Phylogenetics: Simulating sequence data
|EEB 5349: Phylogenetics|
by Kevin Keegan and Paul Lewis
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).
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.
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.
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?
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.