Preliminaries
Download the materials for this practical and expand the archive:
$ cd ~/workshop_materials/15_svdq_astral
$ wget phylosolutions.com/tutorials/wpsg2018/svdquartets_tutorial.zip
$ unzip svdquartets_tutorial.zipMove to the directory you just created:
$ cd svdquartets_tutorialIf you have not already done so, install FigTree on your local machine (http://tree.bio.ed.ac.uk/software/figtree/). FigTree is not a critical component of the tutorial, however, so don't sweat it if you run into problems.
In the instructions which follow, the "$" represents the Unix prompt. For example, if the instruction says:
$ cat filename.txt
you would type "cat filename.txt" and hit return. Likewise, "paup>" represents the prompt used by the PAUP* program.
The programs we will use in this practical are installed on the Amazon EC2 instances, and during the workshop you should ignore the following paragraph.
After the workshop, you can run the tutorial locally on your own machines by downloading and installing PAUP* (http://phylosolutions.com/paup-test/), ASTRAL (https://github.com/smirarab/ASTRAL/blob/master/astral-tutorial.md#installation), and RAxML (https://sco.h-its.org/exelixis/web/software/raxml/). Mac and Windows users should use the command-line (console) version of PAUP* rather than the GUI version. The tutorial assumes that PAUP*, ASTRAL, and RAxML can be started on your system using the commands paup, astral, and raxmlHPC, respectively. It will be easiest to set up symbolic links (symlinks) for PAUP* and RAxML; if you do not do this, you may need to modify the commands needed to start the programs from what is given in the tutorial instructions below. If you do not know how to create a symbolic link, "use the google." For ASTRAL, Linux/macOS users can create a script containing the two lines, substituting appropriately for /path/to/astral and modifying the ASTRAL version number if necessary:
#!/bin/sh
java -jar /path/to/astral/astral.4.11.1.jar $@
Windows users would need to create a batch file that runs the "java -jar" command above and passes all arguments from the command line analogously to the " $@" (again, google it).
Part 1: Analysis of a data set simulated from the "anomaly zone"
We will begin with an example used by Liu and Edwards (2009), calling attention to problems that can arise due to the possibility that gene trees for individual loci to conflict with the species tree due to the phenomenon of "incomplete lineage sorting" (ILS). When internal branches of the species tree are very short (or effective population sizes are extremely large), the most probable gene tree can be inconsistent with the species tree, which can cause "concatenation" methods for inferring the tree to fail. See Kubatko and Degnan (2007) and Liu and Edwards (2009) for details.
The data file anomaly_zone.nex represents one replicate of a simulation described in Fig. 2 in Liu and Edwards (2009). This file contains data for 10,000 loci, with 500 sites per locus. We will use my PAUP* program (http://paup.phylosolutions.com) to analyze this data set (PAUP* will subsequently be referred to as PAUP, without the asterisk). It is often convenient to run PAUP from a Nexus-file script containing the commands used to perform the analysis, but we will usually issue the commands interactively here for pedagogic reasons.
Start PAUP and load the data file.
$ paup anomaly_zone.nex -L anomaly.logThe first command argument above is the name of the data file, and the "-L" option opens a log file that will preserve all of the output of your run (you can use any name you like).
We'll begin with a "concatenated" maximum-likelihood analysis to demonstrate the problem. In order to perform ML in PAUP, you first need to estimate a substitution model, which
specifies the rate of substitution from each nucleotide to each other nucleotide, the equilibrium base frequencies, and the distribution of rate variation across sites. PAUP provides a tool to automatically compare models using the AIC or BIC model selection criteria. A reasonable starting tree is needed (but the details of the tree topology have little effect on the chosen model). We'll just do a quick neighbor-joining analysis using LogDet distances (which are a good general purpose distance when you don't know what else to use):
paup> dset distance=logdet;
paup> nj;
paup> automodel;Using the default options, automodel will evaluate a set of 56 models and choose the one with that minimizes the AIC. On completion, the model parameters will be set to their maximum-likelihood estimates, and we can use this as the starting point for a maximum-likelihood search:
paup> set criterion=likelihood;
paup> hsearch;(Note that you can abbreviate PAUP commands as long as they do not become ambiguous. For instance, instead of the above, we could have instead used set crit=like; hs;)
One way to see the tree found by the ML search is to use the describetrees command, with branch lengths drawn proportionally to the estimated number of substitutions/site:
paup> describe/plot=p;Note the relationships among the taxa implied by the topology of this tree. Now we'll do an svdquartets analysis (Chifman and Kubatko, 2014a,b)for the same data. For any command in PAUP, you can see the available options by typing "command-name ?;"
Before we start our next analysis, look at the options available for the svdquartets command:
paup> svdquartets ?;For example, you can see that it is possible to evaluate quartets exhaustively ("evalquartets=all") or sample quartets randomly ("evalquartets=random nquartets=n"). Typically, you will leave most options set at their default settings, and enter values only for those you wish to override. The current default setting for evalquartets is random (this may change in a future version), so we will explicitly request evaluation of all possible quartets. We'll also do a bootstrap analysis to assess confidence in the result
paup> svdq evalq=all bootstrap;The analysis finishes quickly and outputs the estimated tree, which is the one compatible with all 5 of the inferred quartet relationships. In fact, this is the correct tree that generated the data (see Liu and Edwards, 2009). The concatenated ML method estimated the tree incorrectly despite the large number of sites (500,000). Write down the bootstrap proportions for the two internal nodes for discussion later.
ASTRAL (Mirarab et al., 2014, 2015, 2017) is a "summary" method for estimating species trees from a set of gene trees estimated for each locus. Although it is very fast once the gene trees have been estimated, the preliminary phase of gene tree estimation (outside of ASTRAL) can take a long time. We'll use the RAxML program (https://sco.h-its.org/exelixis/web/software/raxml/index.html) to estimate these trees (PAUP could also be used to estimate the trees, but if you have a data set with many tips, RAxML is much faster, and this example will show you how to do it.)
PAUP now provides a wrapper for RAxML that simplifies the process of obtaining the gene trees. I wrote a small Python script to generate the PAUP commands; take a quick look at it:
paup> !cat astral_prep.pyCommands typed at the command line beginning with "!" are passed to a Unix shell, so this command just shows the contents of the file on your terminal. The Python program just loops through the loci and outputs commands that analyze the data for each locus.
You will remain in the shell until you hit the return key at the "!" prompt, so do this now. Now run the python script from the "paup>" prompt to generate a Nexus command file, and hit return to go back to PAUP.
paup> !python astral_prep.py > run_astral.nexThe created Nexus file simply defines a "character partition" specifying the site numbers for each locus and runs a set of three commands for each of the 10,000 loci, e.g.:
include loci.locus_130/only;
raxml;
savetree file=az_tree.tre format=newick append;
I.e., all of the sites except those belonging to a single locus are excluded, RAxML is executed using PAUP's current model settings, and the resulting tree is appended to a growing treefile that will ultimately contain one tree for each locus.
Now run the Nexus command file generated by the Python script (you will be warned that your are going to reset the active data file; this is fine). Running the full analysis would take several hours, so read ahead rather than waiting for the command to complete.
paup> execute run_astral.nex;Abort the command by typing ctrl-C after a few trees have been generated. Scroll back in the output and notice how different the trees are from one locus to the next. This variation comes from two sources: (1) the true variation in gene trees coming from the multispecies coalescent process, and (2) error in estimating the gene trees using a relatively small amount of sequence information per locus.
We will pretend that the analysis actually completed, and then quit the program.
paup> quit;As in all good cooking shows, I have pre-baked the result so that we can pull it out of the oven immediately. It is stored in the file az_10000_trees.tre, which contains one tree for each of the 10,000 trees, written in the "Newick" (nested parentheses) format. It looks like this:
(A,(B,(D,E)),C);
(A,((B,D),E),C);
(A,((B,D),C),E);
(A,(B,(D,E)),C);
(A,((B,E),D),C);
.
.
.
We'll now run this file in ASTRAL. First, look at the options available:
$ astral --helpFor now, we will simply use the default settings to estimate a species tree, using the treefile we generated above as input.
$ astral -i az_10000_trees.treThe analysis runs quickly and outputs a few lines of output. The estimated species tree is found near the end of the output, and looks something like this:
(A,(B,(C,(E,D)1:0.025494729020277836)1:0.033915343495063886));
Like SVDQuartets, ASTRAL infers the (unrooted) species tree correctly. Note the line near the top: "All output trees will be arbitrarily rooted at A". Both the input to and output from ASTRAL are unrooted trees; the rooting implied by the result is meaningless. In our case, the outgroup is tip E, so we would need to reroot the tree appropriately for a figure.
To view the tree, copy the output tree description the clipboard (e.g., ctrl-C), launch FigTree and paste the clipboard contents (e.g., ctrl-V). Select the terminal branch leading to tip E, and click the Reroot tool.
If you want to gain more familiarity with ASTRAL (after our practical session ends), you can run the tutorial at https://github.com/smirarab/ASTRAL/blob/master/astral-tutorial.md. On the EC2 instances, you can use the command "astral" rather than "java -jar astral"5.5.10.jar" as given in the tutorial.
Part 2: Analysis of a real data set
The data set is for members of the family Canidae (dogs, wolves, coyotes, foxes, jackals, etc.). Lindblad-Toh et al. (2005) sequenced a collection of 16 nuclear loci. The alignments used here are from a *BEAST2 tutorial written by Huw Ogilvie (https://github.com/genomescale/starbeast2/releases/download/v0.14.0/StarBEAST2-tutorial.zip). I reformatted them into the Nexus format and renamed the taxa from scientific to common names.
Start PAUP and load the canid sequence data file:
$ paup canids.nex -L canids.logThis file contains data for two individals from each species, so we use a taxon-partition to assign each individual tip to a species. This definition is already in the file and you do not need to type it here. It looks like this:
taxpartition species =
sidestriped_jackal: sidestriped_jackal_a sidestriped_jackal_b,
African_golden_wolf: African_golden_wolf_a African_golden_wolf_b,
coyote: coyote_a coyote_b,
gray_wolf: gray_wolf_a gray_wolf_b,
blackbacked_jackal: blackbacked_jackal_a blackbacked_jackal_b,
Ethiopian_wolf: Ethiopian_wolf_a Ethiopian_wolf_b,
Dhole: Dhole_a Dhole_b,
African_wild_dog: African_wild_dog_a African_wild_dog_b
;Perform an SVDQuartets analysis, referencing this taxon partition:
paup> svdq evalq=all taxpartition=species;Notice that the (North American) gray wolf (Canis lupus) joins with the African golden wolf (Canis anthus) rather than the (North American) coyote (Canis latrans), unlike the result of the *BEAST tutorial (see StarBEAST2-tutorial.pdf. However, this result is not well supported, as we can see by doing a multilocus nonparametric bootstrap, which resamples (with replacement) both loci and sites within loci:
paup> svdq evalq=all taxpartition=species bootstrap=multilocus loci=combined;(The "combined" above is the name of the character partition defining genes in the input file.) Examining the output bootstrap consensus tree, we see that the bootstrap support for the (gray wolf, African golden wolf) grouping is weak, and many bootstrap runs will actually find the (gray wolf, coyote) grouping that is probably correct. Interestingly, the DensiTree plot in the *BEAST tutorial also demonstrates that support for this grouping is very weak.
When you reach this point, quit PAUP (so that it will forget its current settings and we can start fresh):
>paup quit;Turn over your card so that the green side is up. We'll sync up here before starting the last part.
Part 3: Doing your own simulations
PAUP* now supports simulating data under the multispecies coalescent model. We will use the file sim-6tips.nex, and I will walk you through the components of this file:
#NEXUS
begin taxa;
dimensions ntax=6;
taxlabels A B C D E F;
end;
begin trees;
tree 1 = [&R] (A:4.2,((B:1.0,C:1.0):2.2,(D:3.1,(E:3.0,F:3.0):0.1):0.1):1.0);
end;
begin paup;
showtrees/userbrlens;
end;
begin dnasim;
simdata multilocus=y nloci=100 nsitesperlocus=200;
truetree source=memory treenum=1 units=2Ngen
scalebrlen=1
mscoal=y Ne=100000 mu=1e-8
showtruetree=brlens showgenetrees=n storetruetrees=n
seed=1;
lset model=jc nst=1 basefreq=eq; [= Jukes-Cantor model]
beginsim nreps=100 seed=0;
svdq nthreads=2;
tally 'svdq';
endsim;
end;The copy of this file that you have in your directory will not run initially, because I have replaced the nloci=100, nsitesperlocus=200, and scalebrlen=1 with "@" characters. Each of you will run a simulation with a different value for these settings. You can find your settings here: SVDQuartets results
Using a text editor, replace each "@" symbol with your assigned value. If you are not already comfortable using a Unix text editor, I recommend using "nano" (refer to the Unix activity from the first day).
After you have saved the file, execute it in PAUP:
$ paup sim-6tips.nexRecord the accuracy (proportion of correct trees and proportion of correct nodes) at the appropriate location in the Google Sheet. We'll discuss the results when everyone has posted their results.
References
Liu L., Edwards S.V. 2009. Phylogenetic analysis in the anomaly zone. Syst Biol. 58:452–460.