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 paup; set storebrlens nowarntsave; end; begin taxa; dimensions ntax=4; taxlabels A B C D; end; 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; quit; end;
The initial paup 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 storebrlens?
The taxa block contains the number of taxa and the names of the taxa.
- Why is this taxa block necessary?
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, which is best viewed by pasting its contents into an Excel worksheet. You can also view the results directly in your terminal by typing column -t results.txt | less
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)
Execute the NEXUS File
Log on to the cluster, start a qlogin 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.
Note that PAUP* quits after performing the simulations (because we told it to in that final paup block). Open the results.txt file in a spreadsheet viewer to see the results.
- 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. You will also probably want to specify only 1 sequence length and 1 simulation replicate. It's a bit easier to comprehend if you tell PAUP* to interleave the sequence when it displays it. Do this as a parameter of the export function. If you add taxa to your tree and taxa block the effects of the the following changes will be more pronounced.
- Make all branches in the true tree long (e.g. 10) and see if the simulated data is what you expect
- Make all branches in the true tree short (e.g. 0.001) and see if the simulated data is what you expect
- Make all branches in the true tree 0.1 but add significant rate heterogeneity (gamma shape 0.01) and see if the simulated data is what you expect
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 1000 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 modify the taxa block accordingly), and you may 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
generatetrees all model=equiprobable; lscores all / shtest RELL;
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.
- How many of the 105 trees were not significant using the SH test for 1000 sites? 5000 sites?
Change to the AU test and see if that produces a different result.