Difference between revisions of "Phylogenetics: Simulating sequence data"

From EEBedia
Jump to: navigation, search
(Submitted for Your Approval: The Felsenstein Zone)
 
(31 intermediate revisions by 2 users not shown)
Line 6: Line 6:
 
|
 
|
 
|}
 
|}
by Kevin Keegan
+
by Kevin Keegan and Paul Lewis
  
 
== Goals ==
 
== Goals ==
Line 24: Line 24:
 
Create an empty text file and add the following lines to it and save it as a .nex file:
 
Create an empty text file and add the following lines to it and save it as a .nex file:
 
  #nexus
 
  #nexus
 
begin paup;
 
    set storebrlens nowarntsave;
 
end;
 
 
begin taxa;
 
    dimensions ntax=4;
 
    taxlabels A B C D;
 
end;
 
 
   
 
   
 
  begin trees;
 
  begin trees;
     tree 1 = [&R] ((A:0.1,B:0.1):0.1,(C:0.1,D:0.1):0.1);
+
     tree 1 = [&R] ((A:0.1,B:0.1):0.05,(C:0.1,D:0.1):0.05);
 
  end;
 
  end;
 
   
 
   
Line 55: Line 46:
 
   
 
   
 
  begin paup;
 
  begin paup;
 +
  set nowarntsave;
 
   quit;
 
   quit;
 
  end;
 
  end;
The initial <tt>paup</tt> block tells PAUP* to store branch lengths of any tree it encounters and to not warn us if there are trees in memory when we try to quit.
 
* ''What is the reason for including <tt>storebrlens</tt>?''
 
The <tt>taxa</tt> block contains the number of taxa and the names of the taxa.
 
* ''Why is this taxa block necessary?''
 
 
The <tt>trees</tt> block contains the description of the true tree that we will use to simulate data. By default, trees are considered unrooted, but the obscure <tt>[&R]</tt> says that this tree is rooted.
 
The <tt>trees</tt> block contains the description of the true tree that we will use to simulate data. By default, trees are considered unrooted, but the obscure <tt>[&R]</tt> says that this tree is rooted.
  
The <tt>beginsim...endsim</tt> loop in the <tt>dnasim</tt> block tells PAUP* to simulate 100 nucleotide data sets using the Jukes-Cantor model with no rate heterogeneity. For each of the 100 data sets simulated, two analyses will be performed: (1) an exhaustive search using parsimony and (2) and exhaustive search using maximum likelihood. The <tt>tally</tt> commands keep track of how many times parsimony and ML infer a tree identical to the true tree used for simulation, and the tallied information is stored in the file <tt>results.txt</tt>, which is best viewed by pasting its contents into an Excel worksheet.
+
The <tt>beginsim...endsim</tt> loop in the <tt>dnasim</tt> block tells PAUP* to simulate 100 nucleotide data sets using the Jukes-Cantor model with no rate heterogeneity. For each of the 100 data sets simulated, two analyses will be performed: (1) an exhaustive search using parsimony and (2) and exhaustive search using maximum likelihood. The <tt>tally</tt> commands keep track of how many times parsimony and ML infer a tree identical to the true tree used for simulation, and the tallied information is stored in the file <tt>results.txt</tt>.  
  
 
For both parsimony and ML, tally calculates the following quantities (where TALLYLABEL is either "parsimony" or "ML-JC"):
 
For both parsimony and ML, tally calculates the following quantities (where TALLYLABEL is either "parsimony" or "ML-JC"):
Line 69: Line 57:
 
* TALLYLABEL_P, the fraction of splits in the true tree that were found in the inferred tree (averaged over number of trees inferred)
 
* TALLYLABEL_P, the fraction of splits in the true tree that were found in the inferred tree (averaged over number of trees inferred)
 
* TALLYLABEL_correct, same as TALLYLABEL_P if no incorrect splits are found in the inferred tree, otherwise 0 (averaged over number of trees inferred)
 
* TALLYLABEL_correct, same as TALLYLABEL_P if no incorrect splits are found in the inferred tree, otherwise 0 (averaged over number of trees inferred)
 +
 +
The final paup block sets nowarntsave, which means PAUP will not warn you if you quit without saving stored trees, then quits.
  
 
====Execute the NEXUS File====
 
====Execute the NEXUS File====
  
Log on to the cluster, start a <tt>qlogin</tt> session, and load the current module of PAUP*. If you forget how to load modules, notice the message that was printed right after you logged on to the cluster. Start PAUP*, and execute your NEXUS file.  
+
Log on to the cluster, start a <tt>srun --partition=mcbstudent --qos=mcbstudent --pty bash</tt> session, and load the current module of PAUP*. (Use <tt>module avail</tt> to see available modules, then use <tt>module load xxxx</tt> to start module <tt>xxxx</tt>.) Start PAUP*, and execute your NEXUS file.  
  
Note that PAUP* quits after performing the simulations (because we told it to in that final paup block). Open the <tt>results.txt</tt> file in a spreadsheet viewer to see the results.
+
Note that PAUP* quits after performing the simulations (because we told it to quit in that final paup block). You can also view the <tt>results.txt</tt> file directly in your terminal by typing <tt>column -t results.txt | less</tt> (The <tt>-t</tt> makes the columns align, and the pipe to <tt>less</tt> causes the output to pause after each page of output is shown. Type <tt>q</tt> to get out once you've reached the bottom.).
  
* ''Did both optimality criteria get the tree correct most of the time?'' {{title|Too easy!|answer}}
+
* ''Did both optimality criteria get the tree correct most of the time?''
  
 
==Enter the Felsenstein Zone==
 
==Enter the Felsenstein Zone==
Line 84: Line 74:
 
Execute <tt>paupsimFZ.nex</tt>, then open the new results.txt file in a spreadsheet viewer.
 
Execute <tt>paupsimFZ.nex</tt>, then open the new results.txt file in a spreadsheet viewer.
  
* ''Did both of optimality criteria produce the same tree''? {{title|Nope|answer}}
+
* ''Did both of optimality criteria produce the same tree''?  
  
 
Change the <tt>simdata nchar=10000;</tt> line to <tt>simdata nchar=(10 100 1000 10000);</tt> and change <tt>output=allreps</tt> to <tt>output=meansonly</tt>. Now PAUP* will simulate data sets of four different sequence lengths and summarize the results rather than spitting out a line for every simulation replicate.
 
Change the <tt>simdata nchar=10000;</tt> line to <tt>simdata nchar=(10 100 1000 10000);</tt> and change <tt>output=allreps</tt> to <tt>output=meansonly</tt>. Now PAUP* will simulate data sets of four different sequence lengths and summarize the results rather than spitting out a line for every simulation replicate.
  
* ''Which (parsimony or ML) appears to be statistically consistent''? {{title|Nope|answer}}
+
* ''Which (parsimony or ML) appears to be statistically consistent? Why?''  
 +
 
 +
Add substantial rate heterogeneity (e.g. gamma shape = 0.01) to the simulated data and analyze the data under both parsimony and ML (using a model that assumes rate homogeneity).
 +
 
 +
* ''Is ML statistically consistent when the model is violated in this way? Why?''
 +
 
 +
==Saving Simulated Data==
 +
 
 +
Can you figure out how to change your NEXUS file so that PAUP* simulates one data set and exports it to a file? Start PAUP* and use the "help" command to figure out how to export data to a file. You will also probably want to specify only 1 sequence length (e.g. <tt>nchar=1000</tt>) and 1 simulation replicate. It's a bit easier to comprehend if you tell PAUP* to interleave the sequences when it exports. Do this as a parameter of the <tt>export</tt> function. You can use PAUP*'s <tt>cstatus</tt> command to help answer the questions about proportion of constant sites.
 +
 
 +
* ''Make all branches in the true tree long (e.g. 10). What is the proportion of constant sites? How many substitutions are simulated, on average, per site over the entire tree?''
 +
* ''Make all branches in the true tree short (e.g. 0.001). What is the proportion of constant sites? How many substitutions are simulated, on average, per site over the entire tree?''
 +
* ''Make all branches in the true tree 10 but add significant rate heterogeneity (gamma shape 0.01). What about the proportion of constant sites now? How many substitutions are simulated, on average, per site over the entire tree? To which of the previous 2 simulated data sets is this data set most similar? Can you explain why?''
 +
 
 +
==Strimmer and Rambaut (2002) Study==
 +
 
 +
Download this paper, which I'll refer to as SR from now on:
 +
 
 +
Strimmer K., and Rambaut A. 2002. Inferring confidence sets of possibly misspecified gene trees. Proc. Biol. Sci. 269:137–142. [https://doi.org/10.1098/rspb.2001.1862 https://doi.org/10.1098/rspb.2001.1862]
 +
 
 +
SR simulated data on the tree shown in Figure 1 of their paper and expected the SH (Shimodaira and Hasegawa, 1999) test to reveal that all three possible resolutions of the polytomy were equally supported by the data. Makes sense, doesn't it? What they found instead was that (unless they simulated 4000 sites or more) all 15 (rooted) trees for the four taxa A, B, C, and D were considered equal by the SH test. They concluded that the SH test has a bias making it overly conservative and this bias dissipates as sequence lengths increase. This result motivated Shimodaira to create the AU (Approximately Unbiased) test (Shimodaira, 2002).
 +
 
 +
Can you recreate SR's results for 500 and 5000 sites (see their Table 3)?
 +
 
 +
To do this, you will need to make your true tree equal the tree SR show in their Figure 1, and you need to make the simulation model equal to the one they used (see the bottom right part of SR p. 140). You can delete the <tt>set criterion=...</tt>, <tt>alltrees</tt>, and <tt>tally</tt> commands inside the <tt>beginsim...endsim</tt> loop, replacing these with
 +
set maxtrees=105;
 +
generatetrees all model=equiprobable;
 +
lscores all / shtest autest RELL bootreps=1000;
 +
which generates all 105 possible trees and tests them all using the SH test. To see the output, you'll need to say <tt>monitor=yes</tt> in your <tt>beginsim</tt> command, and you can delete the <tt>resultsfile</tt> statement.
 +
 
 +
(Note: look at the column labeled <tt>SH</tt>, not the column labeled <tt>wtd-SH</tt>.)
  
==Have Some Fun==
+
* ''How many of the 105 trees were not significant using the SH test for 500 sites? 5000 sites?''
 +
* ''Does the AU test produce a different result?''
  
====Time Flies====
+
==Literature Cited==
-lots of time has passed<br />
+
H Shimodaira and M Hasegawa. 1999. Multiple comparisons of log-likelihoods with applications to phylogenetic inference. Molelcular Biology and Evolution 16:1114–1116.
-have them make all of the branch lengths very long and see the nucleotide data
+
  
====Will this Lab Ever End?====
+
H Shimodaira. 2002. An Approximately Unbiased Test of Phylogenetic Tree Selection. Systematic Biology 51:492–508.
-little time has passed<br />
+
-have them make all of the branch lengths teeny and see the nucleotide data
+
  
====You've...Changed====
+
[[Category: Phylogenetics]]
-lots of evolution has occurred for some taxa but not all<br/>
+
-make branch lengths very long for half of the taxa and look at the nucleotide data
+

Latest revision as of 15:44, 30 March 2020

Adiantum.png EEB 5349: Phylogenetics

by Kevin Keegan and Paul Lewis

Goals

The goal of this lab is to gain experience simulating DNA sequence data using PAUP*, which can be useful in testing null hypotheses of interest (parametric bootstrapping) as well as testing the robustness of models to violations of their assumptions and testing the correctness of software and algorithms. PAUP* is a relatively new way to simulate sequence data. The old workhorse for DNA simulations in phylogenetics is Andrew Rambaut's program seq-gen, which is still available (and still as useful as it always was).

Introduction

The development of models and algorithms of any kind requires testing to see how they perform. All models and algorithms make assumptions: they take the infinite complexity of nature and distill them into few components that the maker of the model/algorithm assumes are important. With models of DNA evolution and phylogenetic inference algorithms, one important way of testing the capability of a model/algorithm is by simulating DNA sequence data based on a known phylogeny, and seeing how the model/algorithm performs. If the model/algorithm allows for the recovery of the known or "true" phylogeny then we can rest assured that our model/algorithm is relatively accurate in its distillation of the complexity of the processes it attempts to capture.

Getting Started

We will be using cutting-edge features in PAUP* -- so cutting edge that you will not be able to find any information about these features anywhere online or by using the help command in PAUP*! So don't get confused when you try to look up some of the components of the NEXUS file you will be using. There are some familiar blocks and commands in the NEXUS file though. Feel free to look at past labs or use the help command to refresh your memory.

Simulation Template

Create an empty text file and add the following lines to it and save it as a .nex file:

#nexus

begin trees;
    tree 1 = [&R] ((A:0.1,B:0.1):0.05,(C:0.1,D:0.1):0.05);
end;

begin dnasim;
    simdata nchar=10000;
    lset nst=1 basefreq=equal rates=equal pinvar=0;
    truetree source=memory treenum=1 showtruetree=brlens;
    beginsim nreps=100 seed=12345 monitor=no resultsfile=(name=results.txt replace output=allreps);
        [parsimony]
            set criterion=parsimony;
            alltrees;
            tally parsimony;
        [likelihood under JC]
            set criterion=likelihood;
            alltrees;
            tally 'ML-JC';
    endsim; 
end;

begin paup;
  set nowarntsave;
  quit;
end;

The trees block contains the description of the true tree that we will use to simulate data. By default, trees are considered unrooted, but the obscure [&R] says that this tree is rooted.

The beginsim...endsim loop in the dnasim block tells PAUP* to simulate 100 nucleotide data sets using the Jukes-Cantor model with no rate heterogeneity. For each of the 100 data sets simulated, two analyses will be performed: (1) an exhaustive search using parsimony and (2) and exhaustive search using maximum likelihood. The tally commands keep track of how many times parsimony and ML infer a tree identical to the true tree used for simulation, and the tallied information is stored in the file results.txt.

For both parsimony and ML, tally calculates the following quantities (where TALLYLABEL is either "parsimony" or "ML-JC"):

  • TALLYLABEL_Ntrees, the number of trees tied for being best (ideally 1)
  • TALLYLABEL_P, the fraction of splits in the true tree that were found in the inferred tree (averaged over number of trees inferred)
  • TALLYLABEL_correct, same as TALLYLABEL_P if no incorrect splits are found in the inferred tree, otherwise 0 (averaged over number of trees inferred)

The final paup block sets nowarntsave, which means PAUP will not warn you if you quit without saving stored trees, then quits.

Execute the NEXUS File

Log on to the cluster, start a srun --partition=mcbstudent --qos=mcbstudent --pty bash session, and load the current module of PAUP*. (Use module avail to see available modules, then use module load xxxx to start module xxxx.) Start PAUP*, and execute your NEXUS file.

Note that PAUP* quits after performing the simulations (because we told it to quit in that final paup block). You can also view the results.txt file directly in your terminal by typing column -t results.txt | less (The -t makes the columns align, and the pipe to less causes the output to pause after each page of output is shown. Type q to get out once you've reached the bottom.).

  • Did both optimality criteria get the tree correct most of the time?

Enter the Felsenstein Zone

As you've learned in lecture, parsimony is particularly prone to long branch attraction while maximum likelihood is able to resist joining long edges if the model is correct in the important details. Copy your NEXUS file to create a file named paupsimFZ.nex. Edit the new file and change two edge lengths to 1.0 in order to create a true tree in the Felsenstein zone.

Execute paupsimFZ.nex, then open the new results.txt file in a spreadsheet viewer.

  • Did both of optimality criteria produce the same tree?

Change the simdata nchar=10000; line to simdata nchar=(10 100 1000 10000); and change output=allreps to output=meansonly. Now PAUP* will simulate data sets of four different sequence lengths and summarize the results rather than spitting out a line for every simulation replicate.

  • Which (parsimony or ML) appears to be statistically consistent? Why?

Add substantial rate heterogeneity (e.g. gamma shape = 0.01) to the simulated data and analyze the data under both parsimony and ML (using a model that assumes rate homogeneity).

  • Is ML statistically consistent when the model is violated in this way? Why?

Saving Simulated Data

Can you figure out how to change your NEXUS file so that PAUP* simulates one data set and exports it to a file? Start PAUP* and use the "help" command to figure out how to export data to a file. You will also probably want to specify only 1 sequence length (e.g. nchar=1000) and 1 simulation replicate. It's a bit easier to comprehend if you tell PAUP* to interleave the sequences when it exports. Do this as a parameter of the export function. You can use PAUP*'s cstatus command to help answer the questions about proportion of constant sites.

  • Make all branches in the true tree long (e.g. 10). What is the proportion of constant sites? How many substitutions are simulated, on average, per site over the entire tree?
  • Make all branches in the true tree short (e.g. 0.001). What is the proportion of constant sites? How many substitutions are simulated, on average, per site over the entire tree?
  • Make all branches in the true tree 10 but add significant rate heterogeneity (gamma shape 0.01). What about the proportion of constant sites now? How many substitutions are simulated, on average, per site over the entire tree? To which of the previous 2 simulated data sets is this data set most similar? Can you explain why?

Strimmer and Rambaut (2002) Study

Download this paper, which I'll refer to as SR from now on:

Strimmer K., and Rambaut A. 2002. Inferring confidence sets of possibly misspecified gene trees. Proc. Biol. Sci. 269:137–142. https://doi.org/10.1098/rspb.2001.1862

SR simulated data on the tree shown in Figure 1 of their paper and expected the SH (Shimodaira and Hasegawa, 1999) test to reveal that all three possible resolutions of the polytomy were equally supported by the data. Makes sense, doesn't it? What they found instead was that (unless they simulated 4000 sites or more) all 15 (rooted) trees for the four taxa A, B, C, and D were considered equal by the SH test. They concluded that the SH test has a bias making it overly conservative and this bias dissipates as sequence lengths increase. This result motivated Shimodaira to create the AU (Approximately Unbiased) test (Shimodaira, 2002).

Can you recreate SR's results for 500 and 5000 sites (see their Table 3)?

To do this, you will need to make your true tree equal the tree SR show in their Figure 1, and you need to make the simulation model equal to the one they used (see the bottom right part of SR p. 140). You can delete the set criterion=..., alltrees, and tally commands inside the beginsim...endsim loop, replacing these with

set maxtrees=105;
generatetrees all model=equiprobable;
lscores all / shtest autest RELL bootreps=1000;

which generates all 105 possible trees and tests them all using the SH test. To see the output, you'll need to say monitor=yes in your beginsim command, and you can delete the resultsfile statement.

(Note: look at the column labeled SH, not the column labeled wtd-SH.)

  • How many of the 105 trees were not significant using the SH test for 500 sites? 5000 sites?
  • Does the AU test produce a different result?

Literature Cited

H Shimodaira and M Hasegawa. 1999. Multiple comparisons of log-likelihoods with applications to phylogenetic inference. Molelcular Biology and Evolution 16:1114–1116.

H Shimodaira. 2002. An Approximately Unbiased Test of Phylogenetic Tree Selection. Systematic Biology 51:492–508.