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
|
/* This is an example HY-PHY Batch File.
It reads a '#' nucleotide dataset data/actin2.flt, parititions
the data into two blocks - first 100 sites and the rest, then
applies HKY85 to both partitions. The only difference between the two
is the fact that we use different equilibrium frequencies for each partition.
Output is printed out as a Newick Style tree with branch lengths
representing the number of expected substitutions per branch (which
is the default setting for nucleotide models w/o rate variation).
Also, the likelihood ratio statistic is evaluated and the P-value
for the test is reported.
Sergei L. Kosakovsky Pond and Spencer V. Muse
December 1999.
*/
/* 1. Read in the data and store the result in DataSet variables.
The variable nuclotideSequence.sites holds the number of sites
in the data set. */
DataSet nucleotideSequence = ReadDataFile ("data/actin2.flt");
/* 2. Filter the data.
"0-99" specifies first one hundred sites and
(siteIndex>=100)&&(siteIndex<nuclotideSequence.sites)
takes the rest of the data.*/
DataSetFilter filteredData1 = CreateFilter (nucleotideSequence,1,"0-99");
DataSetFilter filteredData2 = CreateFilter (nucleotideSequence,1,
(siteIndex>=100)&&(siteIndex<nucleotideSequence.sites));
/* 3. Collect observed nucleotide frequencies from the filtered data. observedFreqs will
store the vector of frequencies. */
HarvestFrequencies (observedFreqs1, filteredData1, 1, 1, 1);
HarvestFrequencies (observedFreqs2, filteredData2, 1, 1, 1);
/* 4. Define the HKY85 substitution matrix. '*' is defined to be -(sum of off-diag row elements) */
HKY85RateMatrix =
{{*,trvs,trst,trvs}
{trvs,*,trvs,trst}
{trst,trvs,*,trvs}
{trvs,trst,trvs,*}};
/*5. Define the HKY85 model. USE_FUNCTION_FREQUENCIES allows the model
to pull its equilibrium frequencies from the vector passed to the lk function.*/
Model HKY85_1 = (HKY85RateMatrix, observedFreqs1);
Tree actinTree1 = (ZMU60513,ZMU60514,((((((ZMU60511,ZMU60510),OSRAC2),(OSRAC1,SVSOAC1)),((OSRAC3,ZMU60507),ZMU60509)),(MZEACT1G,OSRAC7)),ZMU60508));
Model HKY85_2 = (HKY85RateMatrix, observedFreqs2);
Tree actinTree2 = (ZMU60513,ZMU60514,((((((ZMU60511,ZMU60510),OSRAC2),(OSRAC1,SVSOAC1)),((OSRAC3,ZMU60507),ZMU60509)),(MZEACT1G,OSRAC7)),ZMU60508));
/*7. Set up the likelihood function, maximize, print results */
LikelihoodFunction theLnLik = (filteredData1, actinTree1,filteredData2, actinTree2);
Optimize (paramValues, theLnLik);
fprintf (stdout, "\n ----- RESULTS ----- \n", theLnLik);
|
