1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
|
/* This is an example HY-PHY Batch File.
In this file, we will illustrate how to use the "ReplicateConstraint" command to
apply a constraint to the tree branches.
Sergei L. Kosakovsky Pond and Spencer V. Muse
June 2001.
*/
/* 1. Include a file which reads the data and the tree,
defines a model and the likelihood function */
#include "shared.bf";
/* 2. We begin by obtaining the MLE estimates for the
tree with each branch having 2 rates - transitions
and transversions, and printing out the result.
LIKELIHOOD_FUNCTION_OUTPUT = 3 formats the output
as a list of parameters (incl. constraints) for easier
understanding of the effect of constraints */
Optimize (res, hivLik);
LIKELIHOOD_FUNCTION_OUTPUT = 3;
fprintf (stdout, "1).Unconstrained HKY85\n\n", hivLik);
/* 3. ReplicateConstrain can aslo be used to apply a
constraint to a subtree rooted at some internal node.
To do that, we pass the name of the internal node as
the "this" argument replacement as in the following example.
The tree in the example data set has an internal node named
"iNode1".
*/
ReplicateConstraint ("this1.?.trst:=this2.?.trvs", hivTree.iNode1, hivTree.iNode1);
Optimize (res, hivLik);
fprintf (stdout, "\n\n2).HKY85 constrained to F81 on the \"iNode1\" subtree \n\n", hivLik);
/* 4. If you wish to apply the constraint to an internal node itself only,
and not to its descendants, embed the name of the node in the constraint itself.
*/
ClearConstraints(hivTree);
ReplicateConstraint ("this1.iNode1.trst:=this2.iNode1.trvs", hivTree.iNode1, hivTree.iNode1);
Optimize (res, hivLik);
fprintf (stdout, "\n\n3).HKY85 constrained to F81 on the \"iNode1\" branch \n\n", hivLik);
/* 5. ReplicateConstraint can also set up a constraint between
two different subtrees (of the same structure though).
The tree in the example file, has two subtrees of the same structure (,()),
one rooted an iNode1, and the other at iNode2. We set up all the parameters
in subtree at iNode1 to be equal to corresponding parameters at iNode2.
*/
ClearConstraints(hivTree);
ReplicateConstraint ("this1.?.?:=this2.?.?", hivTree.iNode1, hivTree.iNode2);
Optimize (res, hivLik);
fprintf (stdout, "\n\n4).HKY85 with \"iNode1\" subtree constrained equal to \"iNode2\" subtree\n\n", hivLik);
/* 6. A slight twist on the example above, is to make
the parameters in one subtree proportional, rather than equal, to
the corresponding parameters in the second subtree.
*/
global R = 1;
ClearConstraints(hivTree);
ReplicateConstraint ("this1.?.?:=R*this2.?.?", hivTree.iNode1, hivTree.iNode2);
Optimize (res, hivLik);
fprintf (stdout, "\n\n5).HKY85 with \"iNode1\" subtree constrained proportional to \"iNode2\" subtree\n\n", hivLik);
/* 7. Here's an example of a more complicated constraint:
on the subtree starting at "iNode2", transitions and transversion rates are the same,
and equal to the average of the transitions and transversion rates of the corresponding
branch in the subtree starting "iNode1".
*/
ClearConstraints(hivTree);
ReplicateConstraint ("this1.?.trst:=.5*(this2.?.trst+this3.?.trvs)", hivTree.iNode1, hivTree.iNode2, hivTree.iNode2);
ReplicateConstraint ("this1.?.trvs:=this2.?.trst", hivTree.iNode1, hivTree.iNode1);
Optimize (res, hivLik);
fprintf (stdout, "\n\n5).HKY85 with \"iNode1\" subtree constrained to the average of those in the \"iNode2\" subtree\n\n", hivLik);
|
