Dear Celine,
Try the following:
1). Run a local codon model (see Multimedia File Viewing and Clickable Links are available for Registered Members only!! You need to
for some background) on your codon alignment, by running the "AnalyzeCodonData.bf" in "Basic analyses" from the Analysis->Standard Analyses menu. Select MG94xcustom model, then input whatever nucleotide biases you want (010010 for HKY85, 012345 for GTR - see the above link for notation) and a local option when prompted (this will fit a separate synRate and nonSynRate whose ratio is dN/dS).
2). When the analysis is done, look in the Results menu and select Paramertic Bootstrap. Enter the number of replicates and where to save the summary to
3). Wait a while for the simulation to finish, then look at the results file written out - you can import it into a spreadsheet for example. I am pasting a toy example (3 sequence file and 10 replicates below)
Code:Iteration,Ln-likelihood,Syn_sites,NS_sites,Total_sites,givenTree.b.nonSynRate,givenTree.b.synRate,givenTree.og.nonSynRate,givenTree.og.synRate,givenTree.a.nonSynRate,givenTree.a.synRate,TotalLength(b),SynLength(b),NonSynLength(b),TotalLength(og),SynLength(og),NonSynLength(og),TotalLength(a),SynLength(a),NonSynLength(a)
1,-802.55,85.6211,268.55,366,0.0245751,0.860529,0.108228,0.855417,0.0552774,0.150122,0.140447,0.127596,0.0128502,0.18343,0.126838,0.0565919,0.0511638,0.0222596,0.0289042
2,-889.153,90.0714,265.345,366,0.0477494,0.721171,0.113003,1.49053,0.117413,0.056023,0.137,0.111828,0.0251719,0.290699,0.231128,0.0595713,0.0705831,0.00868716,0.061896
3,-830.321,86.12,268.106,366,0.0560825,0.313611,0.105273,1.678,0.0573788,0.31872,0.0754049,0.0464701,0.0289348,0.302955,0.248642,0.0543137,0.0768309,0.0472273,0.0296036
4,-885.733,91.4508,265.847,366,0.0890521,0.473651,0.100368,0.953327,0.0971588,0.417889,0.121161,0.0748995,0.0462615,0.202892,0.150752,0.05214,0.116555,0.0660817,0.0504728
5,-822.126,84.1483,270.297,366,0.0364101,0.791975,0.107599,1.50065,0.0958834,0.0273557,0.132864,0.114116,0.0187484,0.271635,0.21623,0.0554051,0.0533141,0.00394169,0.0493724
6,-876.111,90.3207,264.25,366,0.0448144,0.923687,0.0986096,1.81218,0.074151,0,0.176121,0.153097,0.0230237,0.351022,0.300361,0.0506613,0.0380956,0,0.0380956
7,-818.456,85.3556,269.879,366,0.0261671,0.884305,0.0901809,1.97265,0.0837307,0,0.142059,0.128351,0.0137076,0.333558,0.286317,0.0472409,0.043862,0,0.043862
8,-838.018,84.7777,269.423,366,0.0236206,0.264488,0.0873169,1.7631,0.118703,0.360056,0.0516829,0.0392838,0.0123991,0.307704,0.261869,0.045835,0.115789,0.0534783,0.0623107
9,-822.697,91.6688,264.08,366,0.0163719,1.00311,0.104712,1.57021,0.0678416,0,0.167871,0.159514,0.0083569,0.303143,0.249693,0.0534492,0.0346292,0,0.0346292
10,-766.61,83.4728,270.162,366,0.0223334,0.58132,0.0916118,1.11329,0.0312125,0.135721,0.0961263,0.0842681,0.0118582,0.210024,0.161382,0.0486426,0.0362469,0.0196742,0.0165727
MLE,-850.063,87.0605,266.185,366,0.0334451,0.607245,0.107903,1.49124,0.0810188,0.186748,0.109673,0.092062,0.0176108,0.282899,0.226082,0.0568173,0.0709734,0.0283122,0.0426612
The last line includes the maximum likelihood estimates and every other one - a bunch of numbers for each replicate. In my example, the tree was (a,b,og) - if you want to estimate the omega for say branch 'b', you would divide givenTree.b.nonSynRate/givenTree.b.synRate (again, see the linked PDF - part II - for an expanation of parameter notation); you could then compute the mean of this value and its variance over all replicates.
HTH,
Sergei
Home
Pages: 1
