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
|
/* 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. The first application of replicate constraint will
set the transition rate at each branch to be equal to the
transversion rate at the same branch. The constraint
can be specified as "branch.trst:=branch.trsv". This is
essentially the same as reducing the HKY85 model to the
F81 model.
The first argument of "ReplicateConstraint" is
the constraint inself. "this1" etc. will be replaced
with the first argument that follows (hivTree), and
this2 - with the 2nd argument (also hivTree). This defines
the tree on which "ReplicateConstraint" operates. It will
recurse down the tree and will attempt to match "?.trst"
and "?.trsv" with the names of actual parameters that
are present at the branch. '?' is the wildcard character
(it matches every name including the empty string).
The output following the optimization will reflect
these new constraints.
*/
ReplicateConstraint ("this1.?.trst:=this2.?.trvs", hivTree, hivTree);
Optimize (res, hivLik);
fprintf (stdout, "\n\n2).HKY85 constrained to F81\n\n", hivLik);
/* 4. If we modify the constraint as follows, it will only
be applied to the branches whose name containts the word "Node".
By default every internal branch is named "Node" followed by a
number, so the constraint will only apply to internal nodes.
Note the use of "ClearConstraints" to remove the constraints
generated in the previous step.
*/
ClearConstraints (hivTree);
ReplicateConstraint ("this1.?Node?.trst:=this2.?Node?.trvs", hivTree, hivTree);
Optimize (res, hivLik);
fprintf (stdout, "\n\n3).HKY85 constrained to F81 on internal branches\n\n", hivLik);
/* 5. Next, we implement a slightly more involved constraint.
We wish to impose the constraint that the ratio of trst and
tsvr for each branch is the same accross the entire tree.
This would be the "classic" HKY85 model, with one shared ratio.
We declare the global variable "R" to be this ratio.
*/
global R = 1;
ClearConstraints (hivTree);
ReplicateConstraint ("this1.?.trvs:=R*this2.?.trst", hivTree, hivTree);
Optimize (res, hivLik);
fprintf (stdout, "\n\n4).HKY85 with shared ratio.\n\n", hivLik);
|
