Phylogenetics: BEAST2 Lab

From EEBedia
Jump to: navigation, search
Adiantum.png EEB 5349: Phylogenetics
The goal of this lab exercise is to use BEAST2 to estimate divergence times under a relaxed clock model and to assess whether the strict clock, random local clocks or the uncorrelated lognormal model fits the data best using estimates of the log marginal likelihood.

Getting started

You have the option of performing this lab fully on the cluster or partially on the cluster. If you run it on the cluster, you will need to be able to display graphics locally on your machine. If your machine isn't equipped to do so you will have to run part of this lab (the BEAUTi part) on your own machine.

(q)login to the cluster

Login to the Bioinformatics Facility Dell cluster (the -Y allows graphics to display on your local machine if you have an XWindow client installed):

ssh -Y username@bbcsrv3.biotech.uconn.edu

Then use the qlogin command to find a free node:

qlogin

Once you are transferred to a free node, type

module load java/1.8.0
module load beast/2.5.0

This makes a more recent version of BEAST2 available to you and loads Java 8 (actually version 1.8), which is required by BEAST 2.5.0.

Create a directory

Use the unix mkdir command to create a directory to play in today:

mkdir beastlab
cd beastlab

Download and save the data file

The tutorial uses a data file named primate-mtDNA.nex that is located on the cluster. Here's how to copy it to your beastlab directory:

cp $BEAST/examples/nexus/primate-mtDNA.nex .

The $BEAST part represents a bash variable (bash is the scripting language you are using to communicate with the cluster's operating system). The module load beast/2.5.0 command above created the BEAST variable, which points to the directory on the cluster's hard drive where beast2 was installed. The cp command stands for copy, so the command above will copy the file primate-mtDNA.nex to the current directory, which should be your newly-created beastlab directory (the period at the end of the cp command is shorthand for the current directory).

Awaken the BEAST2

BEAST2 comprises two main programs: BEAST and BEAUTi. BEAUTi is a graphical application that helps you create the xml file that BEAST runs. If you have logged into the cluster using "ssh -Y" and if you have an XWindow client that can handle the display, you can start BEAUTi by simply typing "beauti" at the prompt. If that does not work, the easiest option is to simply download the latest version of BEAST2 to your own laptop and run BEAUTi from there. Here's the web site:

http://www.beast2.org

(Note that if you are running BEAUTi locally, you will need to download the primate-mtDNA.nex file you copied into your beastlab folder to your laptop in order to import it into BEAUTi.)

Divergence Dating Tutorial

The tutorial at the web site address shown below walks you through the process of dating a primate mtDNA tree but does not show you how to assess model fit using marginal likelihood estimation. After you finish the tutorial, I will show you how to estimate the marginal likelihood of the clock model you used using BEAST2.



Important: read before starting the tutorial

In section 2.3 (Setting the clock model), we will try using three different clock models. Each person will be assigned one of the following:

  1. strict clock (set priors to what tutorial advises)
  2. uncorrelated lognormal relaxed clock (set priors to what tutorial advises, except use Gamma with shape=0.001, scale=1000 prior for ucldMean.cclock)
  3. random local clocks (set priors to what tutorial advises, and use Poisson with Lambda=0.5 for RRateChanges.cclock prior)

Use the clock model and settings assigned to you when you get to section 2.3 of the tutorial


Here is the tutorial

http://beast2-dev.github.io/beast-docs/beast2/DivergenceDating/DivergenceDatingTutorial.html

Note that you will run BEAUTi on your own laptop to generate the xml file, then move the xml file (e.g. myfile.xml) to the cluster and run beast as follows:

beast myfile.xml

Marginal Likelihood Estimation

Steps for setting up your XML file for steppingstone

I would be nice to be able to directly compare different clock models to each other. BEAST2 can estimate the marginal likelihood of a model but unfortunately this currently (as of April 2018, BEAST 2.5.0) involves modifying the xml file directly (there is no menu option in BEAUTi that will do this for you). Open your xml file using your favorite text editor (e.g. BBEdit/TextWrangler on Mac, Notepad++ on Windows) and make the following modifications:

  1. Search for "<run" and change it to "<mcmc"
  2. Search for "</run>" and change it to "</mcmc>"
  3. Before "<mcmc", add the following lines
<run  spec='beast.inference.PathSampler' nrOfSteps='32' alpha='0.3' deleteOldLogs='true' rootdir='ss'>
cd $(dir)
java -cp $(java.class.path) beast.app.beastapp.BeastMain $(resume/overwrite) -java -seed $(seed) beast.xml 
  1. After "</mcmc>" add the following line:
</run>

The sections below show what the file looks like before and after these modifications.

XML file before modification

<?xml version="1.0" encoding="UTF-8" standalone="no"?>
<beast ...>
.
.
.
<run id="mcmc" spec="MCMC" chainLength="6000000">
.
.
.
</run>
</beast>

XML file after modification

<?xml version="1.0" encoding="UTF-8" standalone="no"?>
<beast ...>
.
.
.
<run  spec='beast.inference.PathSampler' nrOfSteps='32' alpha='0.3' deleteOldLogs='true' rootdir='ss'>
cd $(dir)
java -cp $(java.class.path) beast.app.beastapp.BeastMain $(resume/overwrite) -java -seed $(seed) beast.xml 
<mcmc id="mcmc" spec="MCMC" chainLength="6000000">
.
.
.
</mcmc>
</run>
</beast>

Explanation

The steppingstone method (called path sampling in beast documentation) requires samples from a series of MCMC analysis, each exploring a slightly different distribution (and, strangely, none of these equals the posterior distribution!). Each distribution represents a power posterior: a distribution like the posterior except that the likelihood is raised to a power between 0.0 and 1.0. Each power posterior encloses the next one, like matryoshka dolls, allowing accurate estimate of the ratio of the areas of that pair of distributions. Multiplying successive ratios results in cancellation of the numerator of one with the denominator of the next, so that the end result is an estimate of the ratio of the area of the posterior kernel to the area of the prior. Because the area of the prior is 1.0, this overall ratio equals the marginal likelihood.

The nrOfSteps run attribute specifies the number of these MCMC analyses to perform, and the alpha run attribute determines the spacing of the distributions. The XML file you created with BEAUTi specifies everything needed for each component MCMC analysis, so our modifications to the XML file mainly involve wrapping the original run in a shell that causes it to be executed nrOfSteps times. The PathSampler module specified in the spec run attribute handles putting the results together to obtain an estimate of the marginal likelihood.

The only run attributes I haven't mentioned are deleteOldLogs and rootdir. The deleteOldLogs attribute comes into play if you try to repeat an analysis. The PathSampling module will not start an analysis if old log files are still lying around unless you've told it that it is okay to delete old logs. The rootdir attribute specifies the directory used to store the results of each separate MCMC analysis. These have to be stored somewhere until all are finished. Here I've had you specify rootdir=ss, which causes an ss folder to appear inside your beastlab folder. Inside that ss folder you should see nrOfSteps folders appear, numbered step0, step1, ..., step<nrOfSteps-1>.

Running your XML file

Upload this new xml file to your beastlab directory (you may wish to move any output files from the previous run first so they will not be overwritten) on the cluster and run it using beast as before. When beast is finished, it will spit out the log marginal likelihood that it estimated.

Further Information (no need to read this for today's lab)

Read these section only if you are trying to run this lab entirely on your own computer (not the cluster) or if you are looking for more information about marginal likelihood estimation in BEAST2.

Marginal likelihood estimation on your own laptop

You used the cluster to perform the BEAST run that estimated the marginal likelihood. Were you to try this on your own computer, you would discover that you need to install the BEASTLabs and MODEL_SELECTION packages before performing the analysis. This is done using BEAUTi's File > Manage Packages... menu item. The packages are normally installed in the .beast directory inside your home directory on whichever computer you use to run BEAUTi. Rather than have everyone do this separately (and because some of you cannot use graphical applications such as BEAUTi if they are run on the cluster), I moved these packages to a place where BEAST can find them on the cluster (/usr/local/share/beast/2.5).

Run options

For reference, here are all the options that you can include in the <run ...> xml tag. We only used the first two in this lab and let the defaults apply to all other options. Note that the term particles is synonymous with steppingstones or ratios.

  • alpha: alpha parameter of Beta(alpha,1) distribution used to space out steps, default 0.3. If alpha <= 0, uniform intervals are used.
  • nrOfSteps: the number of steps to use, default 8
  • rootdir: root directory for storing particle states and log files (default /tmp)
  • mcmc: MCMC analysis used to specify model and operations in each of the particles
  • chainLength: number of sample to run a chain for a single step (default 100000L)
  • burnInPercentage: burn-In Percentage used for analysing log files (default 50)
  • preBurnin: number of samples that are discarded for the first step, but not the others (default 100000)
  • value: script for launching a job:
    • $(dir) is replaced by the directory associated with the particle]
    • $(java.class.path) is replaced by a java class path used to launch this application
    • $(java.library.path) is replaced by a java library path used to launch this application
    • $(seed) is replaced by a random number seed that differs with every launch
    • $(host) is replaced by a host from the list of hosts
  • hosts: comma separated list of hosts. If there are k hosts in the list, for particle i the term $(host) in the script will be replaced by the (i modulo k) host in the list. Note that whitespace is removed
  • doNotRun: Set up all files but do not run analysis if true. This can be useful for setting up an analysis on a cluster (default false)
  • deleteOldLogs: delete existing log files from root dir (default false)
  • posterior2prior: whether to do steps from posterior to prior or the other way around. Going from posterior to prior is biased towards over estimates, while from prior to posterior the ML estimate is biased towards under estimates (default true)