Bayesian Phylogenetic Analysis - Revision history

Christoph: /* Examine marginal and joint distributions */

2008-04-14T03:28:38Z

‎Examine marginal and joint distributions

Christoph at 04:57, 13 September 2006

2006-09-13T04:57:21Z

Christoph: /* Examine marginal and joint distributions */ format

2006-06-09T15:55:10Z

‎Examine marginal and joint distributions: format

Christoph: /* Appendix */ format

2006-06-09T15:47:32Z

‎Appendix: format

Christoph: Added more int.links +Cats

2005-12-30T01:33:19Z

Added more int.links +Cats

Christoph: /* Markov chain Monte Carlo sampling */ more int.links

2005-12-29T23:50:29Z

‎Markov chain Monte Carlo sampling: more int.links

Christoph: /* Markov chain Monte Carlo sampling */ more int.links

2005-12-29T23:33:37Z

‎Markov chain Monte Carlo sampling: more int.links

Christoph: Added many internal links

2005-12-29T22:52:00Z

Added many internal links

Show changes

Christoph: format

2005-12-27T19:46:44Z

format

Christoph: format

2005-12-27T18:27:15Z

format

Show changes

@@ Line 281: / Line 281: @@
 Again find a shell window where MrBayes is not running. In the command below I have arbitrarily set the burnin to 50 in order to save you the hassle of entering the value.
-  '''cat''' hcvsmall.nexus.p | '''grep''' -v '^Gen' | '''gawk''' '{if (NR > 50) print $6}' > rCG.data
+  '''grep''' -v '^Gen' hcvsmall.nexus.p | '''gawk''' '{if (NR > 50) print $6}' > rCG.data
 This command takes the parameter file, removes the header line, and for all lines that were sampled after the burnin period it prints the <tt>r(C<->G)</tt> parameter to the file named "<tt>rCG.data</tt>". Now extract a few extra interesting columns in a similar manner.
-  '''cat''' hcvsmall.nexus.p | '''grep''' -v '^Gen' | '''gawk''' '{if (NR > 50) print $9}'  > rAC.data
+  '''grep''' -v '^Gen' hcvsmall.nexus.p | '''gawk''' '{if (NR > 50) print $9}'  > rAC.data
-  '''cat''' hcvsmall.nexus.p | '''grep''' -v '^Gen' | '''gawk''' '{if (NR > 50) print $14}' > rm1.data
+  '''grep''' -v '^Gen' hcvsmall.nexus.p | '''gawk''' '{if (NR > 50) print $14}' > rm1.data
-  '''cat''' hcvsmall.nexus.p | '''grep''' -v '^Gen' | '''gawk''' '{if (NR > 50) print $15}' > rm2.data
+  '''grep''' -v '^Gen' hcvsmall.nexus.p | '''gawk''' '{if (NR > 50) print $15}' > rm2.data
-  '''cat''' hcvsmall.nexus.p | '''grep''' -v '^Gen' | '''gawk''' '{if (NR > 50) print $16}' > rm3.data
+  '''grep''' -v '^Gen' hcvsmall.nexus.p | '''gawk''' '{if (NR > 50) print $16}' > rm3.data
 We can now plot some interesting distributions:

@@ Line 360: / Line 360: @@
   printf "Plotting column ";
   $command = sprintf("head -2 %s | tail -1 | gawk  \'{print \$%d}\'",$parameterfile, $column);
-  '''system''' ($command);
+  system ($command);
   printf "Press ENTER to exit\n";
   <span style="color:blue;"># Plot</span>
-  '''system''' ("cat $parameterfile | grep -v 'ID:'| sed \'s/^Gen/#Gen/\' > tempmb.data");
+  system ("cat $parameterfile | grep -v 'ID:'| sed \'s/^Gen/#Gen/\' > tempmb.data");
   $command = sprintf("plot %s\'tempmb.data\' u %s:%s w li\n",$range,1,$column);
-  '''system''' ("echo \"$command\" > tempmb.gnu");
+  system ("echo \"$command\" > tempmb.gnu");
-  '''system''' ("echo \"pause -1\" >> tempmb.gnu");
+  system ("echo \"pause -1\" >> tempmb.gnu");
-  '''system''' ("gnuplot tempmb.gnu");
+  system ("gnuplot tempmb.gnu");
-  '''system''' ("rm tempmb.gnu tempmb.data");
+  system ("rm tempmb.gnu tempmb.data");
-[[Category:Academic Research]]
+{{Phylogenetics}}
 [[Category:Phylogenetics]]

@@ Line 281: / Line 281: @@
 Again find a shell window where MrBayes is not running. In the command below I have arbitrarily set the burnin to 50 in order to save you the hassle of entering the value.
-<pre>cat hcvsmall.nexus.p | grep -v '^Gen' | gawk '{if (NR > 50) print $6}' > rCG.data</pre>
+ '''cat''' hcvsmall.nexus.p | '''grep''' -v '^Gen' | '''gawk''' '{if (NR > 50) print $6}' > rCG.data
 This command takes the parameter file, removes the header line, and for all lines that were sampled after the burnin period it prints the <tt>r(C<->G)</tt> parameter to the file named "<tt>rCG.data</tt>". Now extract a few extra interesting columns in a similar manner.
-<pre>
+ '''cat''' hcvsmall.nexus.p | '''grep''' -v '^Gen' | '''gawk''' '{if (NR > 50) print $9}'  > rAC.data
-cat hcvsmall.nexus.p | grep -v '^Gen' | gawk '{if (NR > 50) print $9}' > rAC.data
+ '''cat''' hcvsmall.nexus.p | '''grep''' -v '^Gen' | '''gawk''' '{if (NR > 50) print $14}' > rm1.data
-cat hcvsmall.nexus.p | grep -v '^Gen' | gawk '{if (NR > 50) print $14}' > rm1.data
+ '''cat''' hcvsmall.nexus.p | '''grep''' -v '^Gen' | '''gawk''' '{if (NR > 50) print $15}' > rm2.data
-cat hcvsmall.nexus.p | grep -v '^Gen' | gawk '{if (NR > 50) print $15}' > rm2.data
+ '''cat''' hcvsmall.nexus.p | '''grep''' -v '^Gen' | '''gawk''' '{if (NR > 50) print $16}' > rm3.data
 cat hcvsmall.nexus.p | grep -v '^Gen' | gawk '{if (NR > 50) print $16}' > rm3.data
-</pre>
 We can now plot some interesting distributions:
-<pre>histo -w 2 -L"CG":"AC" -m ll rCG.data rAC.data</pre>
+ '''histo''' -w 2 -L"CG":"AC" -m ll rCG.data rAC.data
 We here use the <tt>histo</tt> programme to plot the marginal distributions of rCG and rAC (issue the command <TT>man histo</TT> if you want to learn how to use the <tt>histo</tt> programme). Note that while rAC tends to be larger than rCG, there is still a substantial overlap between the distributions. We therefore need to examine the joint distribution to find the real probability that rAC > rCG.
-<pre>wc -l rCG.data</pre>
+ '''wc''' -l rCG.data
 This gives us a count of how many points remain in the sample files after having discarded the burnin. Write down the number.
-<pre>paste rAC.data rCG.data | gawk '{if ($1>$2) print $0}' | wc -l</pre>
+ '''paste''' rAC.data rCG.data | '''gawk''' '{if ($1>$2) print $0}' | '''wc''' -l
 This counts the number of sample points where rAC > rCG. You can now find the probability that rAC > rCG by dividing the last number by the first. Note how examining the joint distribution provides you with information that you could not get from simply comparing the marginal distributions. This very simple procedure can also be performed using spread-sheet programs, and can obviously be used to answer many different questions.
@@ Line 308: / Line 306: @@
 We will finish off this tutorial by making one final plot:
-<pre>histo -w 0.01 -L"1st":"2nd":"3rd" -m lll rm1.data rm2.data rm3.data</pre>
+ '''histo''' -w 0.01 -L"1st":"2nd":"3rd" -m lll rm1.data rm2.data rm3.data
 This command plots the relative mutation rates at the first, second, and third codon positions. How does this fit with your knowledge of the structure of the genetic code?

@@ Line 319: / Line 319: @@
 * The '''<tt>mb-plot.pl</tt>''' Perl script:
+ <span style="color:purple;">#! /usr/bin/perl</span>
 # Ugly commandline parsing
-if($#ARGV>3 || $#ARGV<0) {
+ <span style="color:blue;"># Script that plots ".p" file output by MrBayes (v3.0B4). Default is to plot
-        print "Usage: mb-plot.pl <parameterfile> [from to | column | column from to]\n";
+ # LnL in order to figure out how large the "burnin" is.
-        exit;
+ # If one parameter is given, it is assumed to be the column to plot
+ # If two extra parameters are given , they are assumed to be from-to plotrange.
+ # If THREE parameters are given, they are assumed to be column no.
-$range="";
+ # followed by from-to.
 $column=2;
+ # A.G. Pedersen, 2003</span>
 # Otherwise assume everything is just fine...
-$parameterfile=$ARGV[0];
+ <span style="color:blue;"># Ugly commandline parsing</span>
-if ($#ARGV==1) {
+ if($#ARGV>3 || $#ARGV<0) {
-        $column = $ARGV[1];
+         print "Usage: mb-plot.pl <parameterfile> [from to | column | column from to]\n";
+         exit;
-−
+ }
 if ($#ARGV==2) {
-        $from = $ARGV[1];
+ $range="";
-        $to = $ARGV[2];
+ $column=2;
-        $range = sprintf("[%d:%d]",$from,$to);
+ <span style="color:blue;"># Otherwise assume everything is just fine...</span>
+ $parameterfile=$ARGV[0];
-if ($#ARGV==3) {
+ if ($#ARGV==1) {
-        $column = $ARGV[1];
+         $column = $ARGV[1];
-        $from = $ARGV[2];
-        $to = $ARGV[3];
-        $range = sprintf("[%d:%d]",$from,$to);
+ if ($#ARGV==2) {
+         $from = $ARGV[1];
+         $to = $ARGV[2];
-# Print column heading so user knows that she is doing the right thing...
+         $range = sprintf("[%d:%d]",$from,$to);
 printf "Plotting column ";
 $command = sprintf("head -2 %s | tail -1 | gawk  \'{print \$%d}\'",$parameterfile, $column);
-system($command);
+ if ($#ARGV==3) {
-printf "Press ENTER to exit\n";
+         $column = $ARGV[1];
+         $from = $ARGV[2];
-# Plot
+         $to = $ARGV[3];
-system ("cat $parameterfile | grep -v 'ID:'| sed \'s/^Gen/#Gen/\' > tempmb.data");
+         $range = sprintf("[%d:%d]",$from,$to);
 $command = sprintf("plot %s\'tempmb.data\' u %s:%s w li\n",$range,1,$column);
 system ("echo \"$command\" > tempmb.gnu");
-system ("echo \"pause -1\" >> tempmb.gnu");
+ <span style="color:blue;"># Print column heading so user knows that she is doing the right thing...</span>
-system ("gnuplot tempmb.gnu");
+ printf "Plotting column ";
-system("rm tempmb.gnu tempmb.data");
+ $command = sprintf("head -2 %s | tail -1 | gawk  \'{print \$%d}\'",$parameterfile, $column);
-</pre>
+ '''system''' ($command);
 [[Category:Academic Research]]
 [[Category:Phylogenetics]]
 [[Category:Statistics]]

@@ Line 29: / Line 29: @@
 We will now use Bayesian phylogenetic analysis.
-===Start MrBayes===
+===MrBayes===
 We will use the programme, '''[[MrBayes]]'''. For now we will just test the programme by issuing ''sample'' commands given below without wondering what they mean. They will be explained later.
@@ Line 84: / Line 84: @@
 </pre>
-==Posterior probability of trees==
+==[[Posterior probability|Posterior probability]] of trees==
 We will continue to use the programme "'''MrBayes'''" to perform '''Bayesian phylogenetic analysis'''. MrBayes is a programme that, like PAUP*, can be controlled by giving commands at a command line prompt. In fact, there is a substantial overlap between the commands used to control MrBayes and the PAUP command language. This should be a help when you are trying to understand how to use the program.
@@ Line 102: / Line 102: @@
 * '''Load your sequences:'''
-<pre>execute foo_mitDNA.nexus</pre>
+<pre>execute foo_mtDNA.nexus</pre>
 This file contains [[mitochondrial DNA]] (mtDNA) sequences from 5 different primates. Note that MrBayes accepts input in [[NEXUS file format|NEXUS format]], and that this is the same command that was used to load sequences in PAUP*. In general, you can use many of the PAUP commands in MrBayes also.
@@ Line 152: / Line 152: @@
 Tree topology is also a parameter in our model, and exactly like the other parameters we also get samples from tree-space. One tree is printed per line in the parenthetical format used by most (all?) [[Phylogeny software|phylogeny software]]. Since there are 5 taxa in the present data set, the tree-space consists of only 15 different possible trees. Since we have taken 3000 sample points, there must be many lines containing the same tree topology.
-===Determine "burn-in"===
+===Determine "[[Burnin|burn-in]]"===
 We now have to determine the "'''burnin'''" (i.e., the number of sample points we have to exclude from our further analysis). Here we will do that by inspecting the likelihood values of the trees saved by MrBayes.
@@ Line 198: / Line 198: @@
 <pre>sumt contype=halfcompat showtreeprobs=no burnin=XX</pre>
-(Again remember to put your actual number of samples to discard instead of the <tt>XX</tt> in the command line above.) Now examine the consensus tree that is plotted to screen: On the branches that are resolved, you will notice that numbers have been plotted. These are '''clade-credibility values''', and are in fact the '''posterior probability''' that the ''clade is real'' (based on the present data set). These numbers are different from the ''bootstrap values'' that we briefly investigated previously and that we will return to later: unlike bootstrap support (which have no clear statistical meaning) these are actual probabilities. Furthermore, they have been found using a full probabilistic model, instead of neighbour joining, and have still finished in a reasonable amount of time. These features make Bayesian phylogeny very useful for assessing hypotheses about '''monophyly'''.
+(Again remember to put your actual number of samples to discard instead of the <tt>XX</tt> in the command line above.) Now examine the consensus tree that is plotted to screen: On the branches that are resolved, you will notice that numbers have been plotted. These are '''clade-credibility values''', and are in fact the '''posterior probability''' that the ''clade is real'' (based on the present data set). These numbers are different from the ''[[Bootstrapping|bootstrap values]]'' that we briefly investigated previously and that we will return to later: unlike bootstrap support (which have no clear statistical meaning) these are actual probabilities. Furthermore, they have been found using a full probabilistic model, instead of neighbour joining, and have still finished in a reasonable amount of time. These features make Bayesian phylogeny very useful for assessing hypotheses about '''[[monophyly]]'''.
 Based on your values, how much support does the present data set give to the "out-of-Africa hypothesis"?
@@ Line 376: / Line 376: @@
 [[Category:Academic Research]]

@@ Line 132: / Line 132: @@
 <pre>mcmc ngen=300000 samplefreq=100 nchains=4 savebrlens=yes</pre>
-What you are doing here is to use the method known as '''MCMCMC''' ("[[Metropolis-coupled Markov chain Monte Carlo]]" or '''MC<sup>3</sup>''') to empirically determine the '''posterior probability distribution of trees''', branch lengths and substitution parameters. Recall that in the Bayesian framework this is how we learn about parameter values: instead of finding the best point estimates, we typically want to quantify the probability of the entire range of possible values. The run will take about 5 minutes to finish (estimated time to finish is shown in the last column of output).
+What you are doing here is to use the method known as '''MCMCMC''' ("[[Metropolis-coupled Markov chain Monte Carlo]]" or '''(MC)<sup>3</sup>''') to empirically determine the '''posterior probability distribution of trees''', branch lengths and substitution parameters. Recall that in the Bayesian framework this is how we learn about parameter values: instead of finding the best point estimates, we typically want to quantify the probability of the entire range of possible values. The run will take about 5 minutes to finish (estimated time to finish is shown in the last column of output).
 Let us examine the command in detail. First, <tt>nchain=4</tt> means that the MCMCMC sampling uses 4 parallel chains: one "cold" from which sampling takes place, and three "heated" that move around in the parameter space more quickly to find additional peaks in the probability distribution. The options <tt>ngen=300000 samplefreq=100</tt> lets the search run for 300,000 steps ("'''generations'''") and saves a tree once every 100 rounds (meaning that a total of 3,000 trees will be saved).
@@ Line 138: / Line 138: @@
 Usually, you would want to sample from more generations (at least 1,000,000 or so, but the more the better). As you might have guessed, <tt>savebrlens=yes</tt> means we want the branch length information saved along with the sampled trees.
-When the run starts you will see reports about the progress of the four chains. Each line of output lists the '''generation number''' and the '''log likelihoods''' of the current tree/parameter combination for each of the four chains. The '''cold chain''' is the one enclosed in ''brackets'' [], while the '''[[Heated chains|heated chains]]''' are enclosed in ''parentheses'' (). Occasionally the chains will swap so one of the heated chains now becomes cold (and sampling then takes place from this chain).
+When the run starts you will see reports about the progress of the four chains. Each line of output lists the '''generation number''' and the '''[[Likelihood-ratio test|log likelihoods]]''' of the current tree/parameter combination for each of the four chains. The '''cold chain''' is the one enclosed in ''brackets'' [], while the '''[[Heated chains|heated chains]]''' are enclosed in ''parentheses'' (). Occasionally the chains will swap so one of the heated chains now becomes cold (and sampling then takes place from this chain).
 * '''Have a look at the resulting sample files:'''
@@ Line 146: / Line 146: @@
 <pre>vi foo_mtDNA.nexus.p</pre>
-This file contains one line for each sampled point. Each ''row'' corresponds to a certain sample time (or generation). Each ''column'' contains the sampled values of one specific parameter. (The first line contains headings telling what the different columns are.) Note how the values of most parameters change a lot during the initial "'''burnin'''" period, before they settle near their most probable values. Now, have a look at the file containing sampled trees (<tt>.t</tt> extension):
+This file contains one line for each sampled point. Each ''row'' corresponds to a certain sample time (or generation). Each ''column'' contains the sampled values of one specific parameter. (The first line contains headings telling what the different columns are.) Note how the values of most parameters change a lot during the initial "'''[[Burnin|burnin]]'''" period, before they settle near their most probable values. Now, have a look at the file containing sampled trees (<tt>.t</tt> extension):
 <pre>vi foo_mtDNA.nexus.t</pre>
-Tree topology is also a parameter in our model, and exactly like the other parameters we also get samples from tree-space. One tree is printed per line in the parenthetical format used by most (all?) phylogeny software. Since there are 5 taxa in the present data set, the tree-space consists of only 15 different possible trees. Since we have taken 3000 sample points, there must be many lines containing the same tree topology.
+Tree topology is also a parameter in our model, and exactly like the other parameters we also get samples from tree-space. One tree is printed per line in the parenthetical format used by most (all?) [[Phylogeny software|phylogeny software]]. Since there are 5 taxa in the present data set, the tree-space consists of only 15 different possible trees. Since we have taken 3000 sample points, there must be many lines containing the same tree topology.
 ===Determine "burn-in"===

@@ Line 162: / Line 162: @@
 <pre>mb-plot.pl primatemitDNA.nexus.p 0 50000</pre>
-This will show you the <tt>lnL</tt> values from generation 0 to generation 50000. Feel free to use other limits if you want to zoom further in or out. Now, look at the ''x''-axis to determine the number of generations required before the plateau is reached. Note down this number and close the window.
+This will show you the <tt>LnL</tt> values from generation 0 to generation 50000. Feel free to use other limits if you want to zoom further in or out. Now, look at the ''x''-axis to determine the number of generations required before the plateau is reached. Note down this number and close the window.
 ===Determine number of sample points to discard===
@@ Line 257: / Line 257: @@
 <pre>sump burnin=XX</pre>
-The <tt>sump</tt> command works much like the <tt>sumt</tt> command for the non-tree parameters. First, you get a plot of the <tt>lnL</tt> (not as nice as when using <tt>mb-plot.pl</tt> though). Secondly, the posterior probability distribution of each parameter is summarized by giving the '''mean''', '''variance''', '''median''', and '''95% credible interval'''.
+The <tt>sump</tt> command works much like the <tt>sumt</tt> command for the non-tree parameters. First, you get a plot of the <tt>LnL</tt> (not as nice as when using <tt>mb-plot.pl</tt> though). Secondly, the posterior probability distribution of each parameter is summarized by giving the '''mean''', '''variance''', '''median''', and '''95% credible interval'''.
 Based on the reported posterior means, does it seem that <tt>r(C<->G)</tt> is different from <tt>r(A<->C)</tt>?. Here <tt>r(C<->G)</tt> is the relative substitution rate between C and G. (Note that all substitution parameters are given as relative values, with the '''GT rate''' arbitrarily set at a value of 1.00).
@@ Line 319: / Line 319: @@
+#
 # Script that plots ".p" file output by MrBayes (v3.0B4). Default is to plot
-# lnL in order to figure out how large the "burnin" is.
+# LnL in order to figure out how large the "burnin" is.
 # If one parameter is given, it is assumed to be the column to plot
 # If two extra parameters are given , they are assumed to be from-to plotrange.