Exercise 1:The objective of this activity is to use CODEML to evaluate the likelihood of the GstD1 sequences for a variety of ω values. Plot log-likelihood scores against the values of ω and determine the maximum likelihood estimate of ω. Check your finding by running CODEML’s hill-climbing algorithm.
#run this command to execute:
The MLE of ω is 0.04.
Now we are going to have PAML estimate ω. In the control file, set fix_omega=0.
------------------
2 (Sim) ... 1 (Mel)
lnL = -756.568200
0.12527 2.53119 0.06698
t= 0.1253 S= 45.2 N= 554.8 dN/dS= 0.0670 dN= 0.0204 dS= 0.3041
---------------------
Exercise 2: In this exercise you will investigate the sensitivity of your estimate of ω to the transition/transversion ratio (κ), and to the assumed model for codon frequencies (π’s). After you collect the required data you will determine which assumptions yield the largest and smallest values of S, and what effect it as on the value of ω.
Here is the control file:
Test#1 Codon bias = none; Ts/Tv bias = none
CodonFreq = 0 * 0:equal,
fix_kappa = 1 * 1:kappa fixed,
2 (Sim) ... 1 (Mel)
lnL = -927.178200
0.10683 0.27421
t= 0.1068 S= 152.9 N= 447.1 dN/dS= 0.2742 dN= 0.0213 dS= 0.0776
Test#2 Codon bias = none; Ts/Tv bias = Yes (ML)
*Highlighted is the ML estimate of of sequences kappa (ts/tv)
CodonFreq = 0 * 0:equal,
fix_kappa = 0 * 1:kappa fixed, 0:kappa to be estimated
2 (Sim) ... 1 (Mel)
lnL = -926.280300
0.10533 1.88371 0.32004
t= 0.1053 S= 165.8 N= 434.2 dN/dS= 0.3200 dN= 0.0221 dS= 0.0691
Test#3 Codon bias = yes (F3x4); Ts/Tv bias = none
CodonFreq = 2 * 2 F3X4,
fix_kappa = 1 * 1:kappa fixed,
2 (Sim) ... 1 (Mel)
lnL = -844.507294
0.10684 0.11807
t= 0.1068 S= 70.6 N= 529.4 dN/dS= 0.1181 dN= 0.0189 dS= 0.1605
Test#4 Codon bias = yes (F3x4); Ts/Tv bias = Yes (ML)
CodonFreq = 2 * 2 F3X4,
fix_kappa = 0 * 0:kappa to be estimated
2 (Sim) ... 1 (Mel)
lnL = -842.214727
0.10686 2.71034 0.12661
t= 0.1069 S= 73.4 N= 526.6 dN/dS= 0.1266 dN= 0.0193 dS= 0.1526
Test #5 Codon bias = yes (F61); Ts/Tv bias = none
CodonFreq = 3 * 3 F61
fix_kappa = 1 * 1:kappa fixed,
2 (Sim) ... 1 (Mel)
lnL = -758.551863
0.12089 0.06278
t= 0.1209 S= 40.5 N= 559.5 dN/dS= 0.0628 dN= 0.0201 dS= 0.3198
Test #6 Codon bias = yes (F61); Ts/Tv bias = Yes (ML)
CodonFreq = 3 * 3 F61
fix_kappa = 0 * 0:kappa to be estimated
2 (Sim) ... 1 (Mel)
lnL = -756.568200
0.12527 2.53121 0.06698
t= 0.1253 S= 45.2 N= 554.8 dN/dS= 0.0670 dN= 0.0204 dS= 0.3040
Summarized in a table
The assumption that yields the smallest value of S is F61 + k=1. It changes the value of omega by an order of magnitude.
Exercise 3: The objective of this exercise is to use three LRTs to evaluate the following hypotheses: (1) the mutation rate of Ldh-C has increased relative to Ldh-A, (2) a burst of positive selection for functional divergence occurred following the duplication event that gave rise to Ldh-C, and (3) there was a long term shift in selective constraints following the duplication event that gave rise to Ldh-C.
Here is the control file:
You should change the name of the main result file (via outfile= in the control file) or you will overwrite your previous results.
#run this command to execute:
codeml ex1_codeml.ctl
This is what the control file looks like:
seqfile = seqfile.txt * sequence data filename
outfile = results_0.010.txt * main result file name
noisy = 9 * 0,1,2,3,9: how much rubbish on the screen
verbose = 1 * 1:detailed output
runmode = -2 * -2:pairwise
seqtype = 1 * 1:codons
CodonFreq = 3 * 0:equal, 1:F1X4, 2:F3X4, 3:F61
model = 0 *
NSsites = 0 *
icode = 0 * 0:universal code
fix_kappa = 0 * 1:kappa fixed, 0:kappa to be estimated
kappa = 2 * initial or fixed kappa
fix_omega = 1 * 1:omega fixed, 0:omega to be estimated
* omega = 0.001 * 1st fixed omega value [change this]
* EXCERCISE 1
*alternate fixed omega values
*omega = 0.005 * 2nd fixed value
*omega = 0.01 * 3rd fixed value
*omega = 0.05 * 4th fixed value
*omega = 0.10 * 5th fixed value
*omega = 0.20 * 6th fixed value
*omega = 0.40 * 7th fixed value
*omega = 0.80 * 8th fixed value
*omega = 1.60 * 9th fixed value
omega = 2.00 * 10th fixed value
Now we are going to have PAML estimate ω. In the control file, set fix_omega=0.
------------------
2 (Sim) ... 1 (Mel)
lnL = -756.568200
0.12527 2.53119 0.06698
t= 0.1253 S= 45.2 N= 554.8 dN/dS= 0.0670 dN= 0.0204 dS= 0.3041
---------------------
Exercise 2: In this exercise you will investigate the sensitivity of your estimate of ω to the transition/transversion ratio (κ), and to the assumed model for codon frequencies (π’s). After you collect the required data you will determine which assumptions yield the largest and smallest values of S, and what effect it as on the value of ω.
Here is the control file:
seqfile = seqfile.txt * sequence data filename
outfile = results.txt * main result file name
noisy = 9 * 0,1,2,3,9: how much rubbish on the screen
verbose = 1 * 1:detailed output
runmode = -2 * -2:pairwise
seqtype = 1 * 1:codons
CodonFreq = 0 * 0:equal, 1:F1X4, 2:F3X4, 3:F61 [change this]
model = 0 *
NSsites = 0 *
icode = 0 * 0:universal code
fix_kappa = 1 * 1:kappa fixed, 0:kappa to be estimated [change this]
kappa = 1 * initial or fixed kappa
fix_omega = 0 * 1:omega fixed, 0:omega to be estimated
omega = 0.5 * initial omega value
* EXCERCISE 2
* Codon bias = none; Ts/Tv bias = none
* Codon bias = none; Ts/Tv bias = Yes (ML)
* Codon bias = yes (F3x4); Ts/Tv bias = none
* Codon bias = yes (F3x4); Ts/Tv bias = Yes (ML)
* Codon bias = yes (F61); Ts/Tv bias = none
* Codon bias = yes (F61); Ts/Tv bias = Yes (ML)
Test#1 Codon bias = none; Ts/Tv bias = none
CodonFreq = 0 * 0:equal,
fix_kappa = 1 * 1:kappa fixed,
2 (Sim) ... 1 (Mel)
lnL = -927.178200
0.10683 0.27421
t= 0.1068 S= 152.9 N= 447.1 dN/dS= 0.2742 dN= 0.0213 dS= 0.0776
Test#2 Codon bias = none; Ts/Tv bias = Yes (ML)
*Highlighted is the ML estimate of of sequences kappa (ts/tv)
CodonFreq = 0 * 0:equal,
fix_kappa = 0 * 1:kappa fixed, 0:kappa to be estimated
2 (Sim) ... 1 (Mel)
lnL = -926.280300
0.10533 1.88371 0.32004
t= 0.1053 S= 165.8 N= 434.2 dN/dS= 0.3200 dN= 0.0221 dS= 0.0691
Test#3 Codon bias = yes (F3x4); Ts/Tv bias = none
CodonFreq = 2 * 2 F3X4,
fix_kappa = 1 * 1:kappa fixed,
2 (Sim) ... 1 (Mel)
lnL = -844.507294
0.10684 0.11807
t= 0.1068 S= 70.6 N= 529.4 dN/dS= 0.1181 dN= 0.0189 dS= 0.1605
Test#4 Codon bias = yes (F3x4); Ts/Tv bias = Yes (ML)
CodonFreq = 2 * 2 F3X4,
fix_kappa = 0 * 0:kappa to be estimated
2 (Sim) ... 1 (Mel)
lnL = -842.214727
0.10686 2.71034 0.12661
t= 0.1069 S= 73.4 N= 526.6 dN/dS= 0.1266 dN= 0.0193 dS= 0.1526
Test #5 Codon bias = yes (F61); Ts/Tv bias = none
CodonFreq = 3 * 3 F61
fix_kappa = 1 * 1:kappa fixed,
lnL = -758.551863
0.12089 0.06278
t= 0.1209 S= 40.5 N= 559.5 dN/dS= 0.0628 dN= 0.0201 dS= 0.3198
Test #6 Codon bias = yes (F61); Ts/Tv bias = Yes (ML)
CodonFreq = 3 * 3 F61
fix_kappa = 0 * 0:kappa to be estimated
2 (Sim) ... 1 (Mel)
lnL = -756.568200
0.12527 2.53121 0.06698
t= 0.1253 S= 45.2 N= 554.8 dN/dS= 0.0670 dN= 0.0204 dS= 0.3040
Summarized in a table
| Assumptions | K | S | N | dS | dN | omega | lnL |
| Fequal + k=1 | 1 | 152.9 | 447.1 | 0.021 | 0.078 | 0.27 | -927.1782 |
| Fequal + k=estimated | 1.88 | 165.8 | 434.2 | 0.07 | 0.022 | 0.32 | -926.2803 |
| F3X4 + k=1 | 1 | 70.6 | 529.4 | 0.16 | 0.019 | 0.118 | -844.507294 |
| F3X4 + k=estimated | 2.71 | 73.4 | 526.6 | 0.15 | 0.02 | 0.127 | -842.214727 |
| F61 + k=1 | 1 | 40.5 | 559.5 | 0.32 | 0.02 | 0.063 | -758.551863 |
| F61 + k=estimated | 2.53 | 45.2 | 554.8 | 0.3 | 0.02 | 0.067 | -756.5682 |
The assumption that yields the smallest value of S is F61 + k=1. It changes the value of omega by an order of magnitude.
Exercise 3: The objective of this exercise is to use three LRTs to evaluate the following hypotheses: (1) the mutation rate of Ldh-C has increased relative to Ldh-A, (2) a burst of positive selection for functional divergence occurred following the duplication event that gave rise to Ldh-C, and (3) there was a long term shift in selective constraints following the duplication event that gave rise to Ldh-C.
Here is the control file:
seqfile = seqfile.txt * sequence data filename
treefile = treeH0.txt * tree structure file name
outfile = results.txt * main result file name
noisy = 9 * 0,1,2,3,9: how much rubbish on the screen
verbose = 1 * 1:detailed output
runmode = 0 * 0:user defined tree
seqtype = 1 * 1:codons
CodonFreq = 2 * 0:equal, 1:F1X4, 2:F3X4, 3:F61
model = 0 * 0:one omega ratio for all branches
* 1:separate omega for each branch
* 2:user specified dN/dS ratios for branches
NSsites = 0 *
icode = 0 * 0:universal code
fix_kappa = 0 * 1:kappa fixed, 0:kappa to be estimated
kappa = 2 * initial or fixed kappa
fix_omega = 0 * 1:omega fixed, 0:omega to be estimated
omega = 0.2 * initial omega
* comments:
* H0 in Table 3: model = 0
* H1 in Table 3: model = 2
* H2 in Table 3: model = 2
* H3 in Table 3: model = 2
H0: homogeneous selection pressure over the tree
(X02152Hom,U07178Sus,(M22585rab,((NM017025Rat,U13687Mus),(((AF070995C,(X04752Mus,U07177Rat)),(U95378Sus,U13680Hom)),(X53828OG1,U28410OG2)))));
Here is the output (lnL highlighted in yellow):
TREE # 1: (1, 2, (5, ((3, 4), (((6, (8, 9)), (7, 10)), (11, 12))))); MP score: -1
lnL(ntime: 21 np: 23): -6018.633010 +0.000000
13..1 13..2 13..14 14..5 14..15 15..16 16..3 16..4 15..17 17..18 18..19 19..6 19..20 20..8 20..9 18..21 21..7 21..10 17..22 22..11 22..12
0.165986 0.145030 0.060255 0.252090 0.030239 0.230093 0.133935 0.101768 0.125411 0.692354 0.249468 0.205980 0.294515 0.131299 0.213038 0.174959 0.155890 0.172753 0.321380 0.290395 0.356759 2.401783 0.136656
Note: Branch length is defined as number of nucleotide substitutions per codon (not per neucleotide site).
tree length = 4.50360
(1: 0.165986, 2: 0.145030, (5: 0.252090, ((3: 0.133935, 4: 0.101768): 0.230093, (((6: 0.205980, (8: 0.131299, 9: 0.213038): 0.294515): 0.249468, (7: 0.155890, 10: 0.172753): 0.174959): 0.692354, (11: 0.290395, 12: 0.356759): 0.321380): 0.125411): 0.030239): 0.060255);
(X02152Hom: 0.165986, U07178Sus: 0.145030, (M22585rab: 0.252090, ((NM017025Rat: 0.133935, U13687Mus: 0.101768): 0.230093, (((AF070995C: 0.205980, (X04752Mus: 0.131299, U07177Rat: 0.213038): 0.294515): 0.249468, (U95378Sus: 0.155890, U13680Hom: 0.172753): 0.174959): 0.692354, (X53828OG1: 0.290395, U28410OG2: 0.356759): 0.321380): 0.125411): 0.030239): 0.060255);
Detailed output identifying parameters
kappa (ts/tv) = 2.40178
omega (dN/dS) = 0.13666
dN & dS for each branch
branch t N S dN/dS dN dS N*dN S*dS
13..1 0.166 729.3 269.7 0.1367 0.0205 0.1497 14.9 40.4
13..2 0.145 729.3 269.7 0.1367 0.0179 0.1308 13.0 35.3
13..14 0.060 729.3 269.7 0.1367 0.0074 0.0543 5.4 14.7
14..5 0.252 729.3 269.7 0.1367 0.0311 0.2273 22.7 61.3
14..15 0.030 729.3 269.7 0.1367 0.0037 0.0273 2.7 7.4
15..16 0.230 729.3 269.7 0.1367 0.0283 0.2075 20.7 55.9
16..3 0.134 729.3 269.7 0.1367 0.0165 0.1208 12.0 32.6
16..4 0.102 729.3 269.7 0.1367 0.0125 0.0918 9.1 24.7
15..17 0.125 729.3 269.7 0.1367 0.0155 0.1131 11.3 30.5
17..18 0.692 729.3 269.7 0.1367 0.0853 0.6242 62.2 168.3
18..19 0.249 729.3 269.7 0.1367 0.0307 0.2249 22.4 60.7
19..6 0.206 729.3 269.7 0.1367 0.0254 0.1857 18.5 50.1
19..20 0.295 729.3 269.7 0.1367 0.0363 0.2655 26.5 71.6
20..8 0.131 729.3 269.7 0.1367 0.0162 0.1184 11.8 31.9
20..9 0.213 729.3 269.7 0.1367 0.0262 0.1921 19.1 51.8
18..21 0.175 729.3 269.7 0.1367 0.0216 0.1577 15.7 42.5
21..7 0.156 729.3 269.7 0.1367 0.0192 0.1406 14.0 37.9
21..10 0.173 729.3 269.7 0.1367 0.0213 0.1558 15.5 42.0
17..22 0.321 729.3 269.7 0.1367 0.0396 0.2898 28.9 78.1
22..11 0.290 729.3 269.7 0.1367 0.0358 0.2618 26.1 70.6
22..12 0.357 729.3 269.7 0.1367 0.0440 0.3217 32.1 86.7
tree length for dN: 0.55489
tree length for dS: 4.06048
H1: episodic change in selection pressure in Ldh-C (only in the branch that immediately follows the gene duplication event).
Here is the tree file for this:
(X02152Hom,U07178Sus,(M22585rab,((NM017025Rat,U13687Mus),(((AF070995C,(X04752Mus,U07177Rat)),(U95378Sus,U13680Hom))#1,(X53828OG1,U28410OG2)))));
Output:
TREE # 1: (1, 2, (5, ((3, 4), (((6, (8, 9)), (7, 10)), (11, 12))))); MP score: -1
lnL(ntime: 21 np: 24): -6017.572460 +0.000000
13..1 13..2 13..14 14..5 14..15 15..16 16..3 16..4 15..17 17..18 18..19 19..6 19..20 20..8 20..9 18..21 21..7 21..10 17..22 22..11 22..12
0.166139 0.145654 0.060220 0.253769 0.027868 0.234240 0.134451 0.101620 0.136088 0.615641 0.250544 0.206999 0.295643 0.131224 0.213679 0.176223 0.155609 0.173597 0.319034 0.291976 0.357768 2.397123 0.131877 0.192033
Note: Branch length is defined as number of nucleotide substitutions per codon (not per nucleotide site).
tree length = 4.44798
(1: 0.166139, 2: 0.145654, (5: 0.253769, ((3: 0.134451, 4: 0.101620): 0.234240, (((6: 0.206999, (8: 0.131224, 9: 0.213679): 0.295643): 0.250544, (7: 0.155609, 10: 0.173597): 0.176223): 0.615641, (11: 0.291976, 12: 0.357768): 0.319034): 0.136088): 0.027868): 0.060220);
(X02152Hom: 0.166139, U07178Sus: 0.145654, (M22585rab: 0.253769, ((NM017025Rat: 0.134451, U13687Mus: 0.101620): 0.234240, (((AF070995C: 0.206999, (X04752Mus: 0.131224, U07177Rat: 0.213679): 0.295643): 0.250544, (U95378Sus: 0.155609, U13680Hom: 0.173597): 0.176223): 0.615641, (X53828OG1: 0.291976, U28410OG2: 0.357768): 0.319034): 0.136088): 0.027868): 0.060220);
Detailed output identifying parameters
kappa (ts/tv) = 2.39712
dN & dS for each branch
branch t N S dN/dS dN dS N*dN S*dS
13..1 0.166 729.4 269.6 0.1319 0.0199 0.1512 14.5 40.8
13..2 0.146 729.4 269.6 0.1319 0.0175 0.1326 12.8 35.7
13..14 0.060 729.4 269.6 0.1319 0.0072 0.0548 5.3 14.8
14..5 0.254 729.4 269.6 0.1319 0.0305 0.2310 22.2 62.3
14..15 0.028 729.4 269.6 0.1319 0.0033 0.0254 2.4 6.8
15..16 0.234 729.4 269.6 0.1319 0.0281 0.2132 20.5 57.5
16..3 0.134 729.4 269.6 0.1319 0.0161 0.1224 11.8 33.0
16..4 0.102 729.4 269.6 0.1319 0.0122 0.0925 8.9 24.9
15..17 0.136 729.4 269.6 0.1319 0.0163 0.1239 11.9 33.4
17..18 0.616 729.4 269.6 0.1920 0.0961 0.5004 70.1 134.9
18..19 0.251 729.4 269.6 0.1319 0.0301 0.2281 21.9 61.5
19..6 0.207 729.4 269.6 0.1319 0.0249 0.1884 18.1 50.8
19..20 0.296 729.4 269.6 0.1319 0.0355 0.2691 25.9 72.6
20..8 0.131 729.4 269.6 0.1319 0.0158 0.1195 11.5 32.2
20..9 0.214 729.4 269.6 0.1319 0.0257 0.1945 18.7 52.4
18..21 0.176 729.4 269.6 0.1319 0.0212 0.1604 15.4 43.3
21..7 0.156 729.4 269.6 0.1319 0.0187 0.1417 13.6 38.2
21..10 0.174 729.4 269.6 0.1319 0.0208 0.1580 15.2 42.6
17..22 0.319 729.4 269.6 0.1319 0.0383 0.2904 27.9 78.3
22..11 0.292 729.4 269.6 0.1319 0.0351 0.2658 25.6 71.7
22..12 0.358 729.4 269.6 0.1319 0.0430 0.3257 31.3 87.8
tree length for dN: 0.55619
tree length for dS: 3.98922
dS tree:
(X02152Hom: 0.151246, U07178Sus: 0.132597, (M22585rab: 0.231020, ((NM017025Rat: 0.122398, U13687Mus: 0.092510): 0.213241, (((AF070995C: 0.188443, (X04752Mus: 0.119460, U07177Rat: 0.194523): 0.269141): 0.228084, (U95378Sus: 0.141659, U13680Hom: 0.158035): 0.160426): 0.500425, (X53828OG1: 0.265802, U28410OG2: 0.325696): 0.290434): 0.123889): 0.025370): 0.054822);
dN tree:
(X02152Hom: 0.019946, U07178Sus: 0.017487, (M22585rab: 0.030466, ((NM017025Rat: 0.016142, U13687Mus: 0.012200): 0.028122, (((AF070995C: 0.024851, (X04752Mus: 0.015754, U07177Rat: 0.025653): 0.035494): 0.030079, (U95378Sus: 0.018682, U13680Hom: 0.020841): 0.021157): 0.096098, (X53828OG1: 0.035053, U28410OG2: 0.042952): 0.038302): 0.016338): 0.003346): 0.007230);
w ratios as labels for TreeView:
(X02152Hom #0.1319 : 0.019946, U07178Sus #0.1319 : 0.017487, (M22585rab #0.1319 : 0.030466, ((NM017025Rat #0.1319 : 0.016142, U13687Mus #0.1319 : 0.012200) #0.1319 : 0.028122, (((AF070995C #0.1319 : 0.024851, (X04752Mus #0.1319 : 0.015754, U07177Rat #0.1319 : 0.025653) #0.1319 : 0.035494) #0.1319 : 0.030079, (U95378Sus #0.1319 : 0.018682, U13680Hom #0.1319 : 0.020841) #0.1319 : 0.021157)#1 #0.1920 : 0.096098, (X53828OG1 #0.1319 : 0.035053, U28410OG2 #0.1319 : 0.042952) #0.1319 : 0.038302) #0.1319 : 0.016338) #0.1319 : 0.003346) #0.1319 : 0.007230);
H2: Long term shift in selection pressure in Ldh-C only; Ldh-C has a permanent change in selection pressure (as compared to its ancestors) whereas Ldh-A remains subject to the ancestral level of selection pressure.
(X02152Hom,U07178Sus,(M22585rab,((NM017025Rat,U13687Mus),(((AF070995C #1,(X04752Mus #1,U07177Rat #1)#1)#1,(U95378Sus #1,U13680Hom #1)#1)#1,(X53828OG1,U28410OG2)))));
Output:
TREE # 1: (1, 2, (5, ((3, 4), (((6, (8, 9)), (7, 10)), (11, 12))))); MP score: -1
lnL(ntime: 21 np: 24): -5985.635237 +0.000000
13..1 13..2 13..14 14..5 14..15 15..16 16..3 16..4 15..17 17..18 18..19 19..6 19..20 20..8 20..9 18..21 21..7 21..10 17..22 22..11 22..12
0.171224 0.151409 0.058980 0.269203 0.017888 0.262039 0.138108 0.103579 0.168914 0.560187 0.231662 0.203731 0.272797 0.130747 0.204423 0.162430 0.151155 0.167095 0.366793 0.312139 0.396584 2.403414 0.075516 0.239425
Note: Branch length is defined as number of nucleotide substitutions per codon (not per neucleotide site).
tree length = 4.50108
(1: 0.171224, 2: 0.151409, (5: 0.269203, ((3: 0.138108, 4: 0.103579): 0.262039, (((6: 0.203731, (8: 0.130747, 9: 0.204423): 0.272797): 0.231662, (7: 0.151155, 10: 0.167095): 0.162430): 0.560187, (11: 0.312139, 12: 0.396584): 0.366793): 0.168914): 0.017888): 0.058980);
(X02152Hom: 0.171224, U07178Sus: 0.151409, (M22585rab: 0.269203, ((NM017025Rat: 0.138108, U13687Mus: 0.103579): 0.262039, (((AF070995C: 0.203731, (X04752Mus: 0.130747, U07177Rat: 0.204423): 0.272797): 0.231662, (U95378Sus: 0.151155, U13680Hom: 0.167095): 0.162430): 0.560187, (X53828OG1: 0.312139, U28410OG2: 0.396584): 0.366793): 0.168914): 0.017888): 0.058980);
Detailed output identifying parameters
kappa (ts/tv) = 2.40341
dN & dS for each branch
branch t N S dN/dS dN dS N*dN S*dS
13..1 0.171 729.3 269.7 0.0755 0.0133 0.1756 9.7 47.3
13..2 0.151 729.3 269.7 0.0755 0.0117 0.1552 8.5 41.9
13..14 0.059 729.3 269.7 0.0755 0.0046 0.0605 3.3 16.3
14..5 0.269 729.3 269.7 0.0755 0.0208 0.2760 15.2 74.4
14..15 0.018 729.3 269.7 0.0755 0.0014 0.0183 1.0 4.9
15..16 0.262 729.3 269.7 0.0755 0.0203 0.2687 14.8 72.5
16..3 0.138 729.3 269.7 0.0755 0.0107 0.1416 7.8 38.2
16..4 0.104 729.3 269.7 0.0755 0.0080 0.1062 5.8 28.6
15..17 0.169 729.3 269.7 0.0755 0.0131 0.1732 9.5 46.7
17..18 0.560 729.3 269.7 0.2394 0.1005 0.4198 73.3 113.2
18..19 0.232 729.3 269.7 0.2394 0.0416 0.1736 30.3 46.8
19..6 0.204 729.3 269.7 0.2394 0.0366 0.1527 26.7 41.2
19..20 0.273 729.3 269.7 0.2394 0.0490 0.2045 35.7 55.1
20..8 0.131 729.3 269.7 0.2394 0.0235 0.0980 17.1 26.4
20..9 0.204 729.3 269.7 0.2394 0.0367 0.1532 26.8 41.3
18..21 0.162 729.3 269.7 0.2394 0.0291 0.1217 21.3 32.8
21..7 0.151 729.3 269.7 0.2394 0.0271 0.1133 19.8 30.6
21..10 0.167 729.3 269.7 0.2394 0.0300 0.1252 21.9 33.8
17..22 0.367 729.3 269.7 0.0755 0.0284 0.3761 20.7 101.4
22..11 0.312 729.3 269.7 0.0755 0.0242 0.3200 17.6 86.3
22..12 0.397 729.3 269.7 0.0755 0.0307 0.4066 22.4 109.7
tree length for dN: 0.56113
tree length for dS: 4.04015
dS tree:
(X02152Hom: 0.175561, U07178Sus: 0.155244, (M22585rab: 0.276022, ((NM017025Rat: 0.141606, U13687Mus: 0.106203): 0.268676, (((AF070995C: 0.152691, (X04752Mus: 0.097991, U07177Rat: 0.153209): 0.204454): 0.173624, (U95378Sus: 0.113287, U13680Hom: 0.125233): 0.121737): 0.419845, (X53828OG1: 0.320046, U28410OG2: 0.406630): 0.376085): 0.173193): 0.018341): 0.060474);
dN tree:
(X02152Hom: 0.013258, U07178Sus: 0.011723, (M22585rab: 0.020844, ((NM017025Rat: 0.010693, U13687Mus: 0.008020): 0.020289, (((AF070995C: 0.036558, (X04752Mus: 0.023462, U07177Rat: 0.036682): 0.048951): 0.041570, (U95378Sus: 0.027124, U13680Hom: 0.029984): 0.029147): 0.100521, (X53828OG1: 0.024169, U28410OG2: 0.030707): 0.028400): 0.013079): 0.001385): 0.004567);
w ratios as labels for TreeView:
(X02152Hom #0.0755 : 0.013258, U07178Sus #0.0755 : 0.011723, (M22585rab #0.0755 : 0.020844, ((NM017025Rat #0.0755 : 0.010693, U13687Mus #0.0755 : 0.008020) #0.0755 : 0.020289, (((AF070995C#1 #0.2394 : 0.036558, (X04752Mus#1 #0.2394 : 0.023462, U07177Rat#1 #0.2394 : 0.036682)#1 #0.2394 : 0.048951)#1 #0.2394 : 0.041570, (U95378Sus#1 #0.2394 : 0.027124, U13680Hom#1 #0.2394 : 0.029984)#1 #0.2394 : 0.029147)#1 #0.2394 : 0.100521, (X53828OG1 #0.0755 : 0.024169, U28410OG2 #0.0755 : 0.030707) #0.0755 : 0.028400) #0.0755 : 0.013079) #0.0755 : 0.001385) #0.0755 : 0.004567);
Time used: 0:56
H3: Long term shift in selection on both Ldh-C and Ldh-A; those lineages are subject to selection pressures different from each other and from the ancestor.
Tree file:
(X02152Hom,U07178Sus,(M22585rab,((NM017025Rat,U13687Mus),(((AF070995C #1,(X04752Mus #1,U07177Rat #1)#1)#1,(U95378Sus #1,U13680Hom #1)#1)#1,(X53828OG1 #2,U28410OG2 #2)#2))));
Output:
TREE # 1: (1, 2, (5, ((3, 4), (((6, (8, 9)), (7, 10)), (11, 12))))); MP score: -1
lnL(ntime: 21 np: 25): -5984.111015 +0.000000
13..1 13..2 13..14 14..5 14..15 15..16 16..3 16..4 15..17 17..18 18..19 19..6 19..20 20..8 20..9 18..21 21..7 21..10 17..22 22..11 22..12
0.173367 0.152089 0.058499 0.273176 0.018609 0.265853 0.138699 0.104431 0.178874 0.556050 0.231889 0.203607 0.272766 0.130870 0.204220 0.161922 0.150629 0.167529 0.341576 0.302793 0.379370 2.399088 0.064218 0.240341 0.094898
Note: Branch length is defined as number of nucleotide substitutions per codon (not per neucleotide site).
tree length = 4.46682
(1: 0.173367, 2: 0.152089, (5: 0.273176, ((3: 0.138699, 4: 0.104431): 0.265853, (((6: 0.203607, (8: 0.130870, 9: 0.204220): 0.272766): 0.231889, (7: 0.150629, 10: 0.167529): 0.161922): 0.556050, (11: 0.302793, 12: 0.379370): 0.341576): 0.178874): 0.018609): 0.058499);
(X02152Hom: 0.173367, U07178Sus: 0.152089, (M22585rab: 0.273176, ((NM017025Rat: 0.138699, U13687Mus: 0.104431): 0.265853, (((AF070995C: 0.203607, (X04752Mus: 0.130870, U07177Rat: 0.204220): 0.272766): 0.231889, (U95378Sus: 0.150629, U13680Hom: 0.167529): 0.161922): 0.556050, (X53828OG1: 0.302793, U28410OG2: 0.379370): 0.341576): 0.178874): 0.018609): 0.058499);
Detailed output identifying parameters
kappa (ts/tv) = 2.39909
dN & dS for each branch
branch t N S dN/dS dN dS N*dN S*dS
13..1 0.173 729.4 269.6 0.0642 0.0117 0.1824 8.5 49.2
13..2 0.152 729.4 269.6 0.0642 0.0103 0.1600 7.5 43.1
13..14 0.058 729.4 269.6 0.0642 0.0040 0.0616 2.9 16.6
14..5 0.273 729.4 269.6 0.0642 0.0185 0.2874 13.5 77.5
14..15 0.019 729.4 269.6 0.0642 0.0013 0.0196 0.9 5.3
15..16 0.266 729.4 269.6 0.0642 0.0180 0.2797 13.1 75.4
16..3 0.139 729.4 269.6 0.0642 0.0094 0.1459 6.8 39.4
16..4 0.104 729.4 269.6 0.0642 0.0071 0.1099 5.1 29.6
15..17 0.179 729.4 269.6 0.0642 0.0121 0.1882 8.8 50.7
17..18 0.556 729.4 269.6 0.2403 0.1000 0.4162 73.0 112.2
18..19 0.232 729.4 269.6 0.2403 0.0417 0.1736 30.4 46.8
19..6 0.204 729.4 269.6 0.2403 0.0366 0.1524 26.7 41.1
19..20 0.273 729.4 269.6 0.2403 0.0491 0.2041 35.8 55.0
20..8 0.131 729.4 269.6 0.2403 0.0235 0.0979 17.2 26.4
20..9 0.204 729.4 269.6 0.2403 0.0367 0.1528 26.8 41.2
18..21 0.162 729.4 269.6 0.2403 0.0291 0.1212 21.2 32.7
21..7 0.151 729.4 269.6 0.2403 0.0271 0.1127 19.8 30.4
21..10 0.168 729.4 269.6 0.2403 0.0301 0.1254 22.0 33.8
17..22 0.342 729.4 269.6 0.0949 0.0319 0.3357 23.2 90.5
22..11 0.303 729.4 269.6 0.0949 0.0282 0.2976 20.6 80.2
22..12 0.379 729.4 269.6 0.0949 0.0354 0.3728 25.8 100.5
tree length for dN: 0.56167
tree length for dS: 3.99725
dS tree:
(X02152Hom: 0.182422, U07178Sus: 0.160033, (M22585rab: 0.287445, ((NM017025Rat: 0.145944, U13687Mus: 0.109885): 0.279740, (((AF070995C: 0.152387, (X04752Mus: 0.097947, U07177Rat: 0.152845): 0.204148): 0.173554, (U95378Sus: 0.112736, U13680Hom: 0.125384): 0.121188): 0.416167, (X53828OG1: 0.297568, U28410OG2: 0.372824): 0.335682): 0.188217): 0.019581): 0.061554);
dN tree:
(X02152Hom: 0.011715, U07178Sus: 0.010277, (M22585rab: 0.018459, ((NM017025Rat: 0.009372, U13687Mus: 0.007057): 0.017964, (((AF070995C: 0.036625, (X04752Mus: 0.023541, U07177Rat: 0.036735): 0.049065): 0.041712, (U95378Sus: 0.027095, U13680Hom: 0.030135): 0.029127): 0.100022, (X53828OG1: 0.028239, U28410OG2: 0.035380): 0.031856): 0.012087): 0.001257): 0.003953);
w ratios as labels for TreeView:
(X02152Hom #0.0642 : 0.011715, U07178Sus #0.0642 : 0.010277, (M22585rab #0.0642 : 0.018459, ((NM017025Rat #0.0642 : 0.009372, U13687Mus #0.0642 : 0.007057) #0.0642 : 0.017964, (((AF070995C#1 #0.2403 : 0.036625, (X04752Mus#1 #0.2403 : 0.023541, U07177Rat#1 #0.2403 : 0.036735)#1 #0.2403 : 0.049065)#1 #0.2403 : 0.041712, (U95378Sus#1 #0.2403 : 0.027095, U13680Hom#1 #0.2403 : 0.030135)#1 #0.2403 : 0.029127)#1 #0.2403 : 0.100022, (X53828OG1#2 #0.0949 : 0.028239, U28410OG2#2 #0.0949 : 0.035380)#2 #0.0949 : 0.031856) #0.0642 : 0.012087) #0.0642 : 0.001257) #0.0642 : 0.003953);
Time used: 1:10
In addition, carry out likelihood ratio tests (LRT) of the hypotheses below. See the lecture notes for additional details about LRTs. Use 1 degree of freedom to obtain the P-value for each LRT.
Summary table of hypotheses:
| Model | ω_A0 | ω_A1 | ω_C1 | ω_C0 | lnL | LRT |
| H0 | 0.1367 | 0.1367 | 0.1367 | 0.1367 | -6018.633 | |
| H1 | 0.1319 | 0.1319 | 0.1319 | 0.192 | -6017.572 | |
| H2 | 0.0755 | 0.0755 | 0.2394 | 0.2394 | -5985.635 | |
| H3 | 0.0949 | 0.0642 | 0.2403 | 0.2403 | -5984.111 |
You should change the name of the main result file (via outfile= in the control file) or you will overwrite your previous results.
Exercise 4: The objective of this exercise is to use a series of LRTs to test for sites evolving under positive selection in the nef gene. If you find significant evidence for positive selection, then identify the involved sites by using empirical Bayes methods.
seqfile = seqfile.txt * sequence data filename
* treefile = treefile_M0.txt * SET THIS for tree file with ML branch lengths under M0
* treefile = treefile_M1.txt * SET THIS for tree file with ML branch lengths under M1
* treefile = treefile_M2.txt * SET THIS for tree file with ML branch lengths under M2
* treefile = treefile_M3.txt * SET THIS for tree file with ML branch lengths under M3
* treefile = treefile_M7.txt * SET THIS for tree file with ML branch lengths under M7
* treefile = treefile_M8.txt * SET THIS for tree file with ML branch lengths under M8
outfile = results.txt * main result file name
noisy = 9 * lots of rubbish on the screen
verbose = 1 * detailed output
runmode = 0 * user defined tree
seqtype = 1 * codons
CodonFreq = 2 * F3X4 for codon ferquencies
model = 0 * one omega ratio for all branches
* NSsites = 0 * SET THIS for M0
* NSsites = 1 * SET THIS for M1
* NSsites = 2 * SET THIS for M2
* NSsites = 3 * SET THIS for M3
* NSsites = 7 * SET THIS for M7
* NSsites = 8 * SET THIS for M8
icode = 0 * universal code
fix_kappa = 1 * kappa fixed
* kappa = 4.43491 * SET THIS to fix kappa at MLE under M0
* kappa = 4.39117 * SET THIS to fix kappa at MLE under M1
* kappa = 5.08964 * SET THIS to fix kappa at MLE under M2
* kappa = 4.89033 * SET THIS to fix kappa at MLE under M3
* kappa = 4.22750 * SET THIS to fix kappa at MLE under M7
* kappa = 4.87827 * SET THIS to fix kappa at MLE under M8
fix_omega = 0 * omega to be estimated
omega = 5 * initial omega
* ncatG = 3 * SET THIS for 3 site categories under M3
* ncatG = 10 * SET THIS for 10 of site categories under M7 and M8
fix_blength = 2 * fixed branch lengths from tree file
sdsf
No comments:
Post a Comment