oldTutorial SIR2 Treelikelihood

XML Generation: Calculating a structured tree likelihood

The goal of these tutorials is to introduce, step-by-step, PhyDyn's (BEAST2) components, analyses and tools. The tutorial sequence also suggests a work-flow, which starts with the definition and execution (test) of the population model; the final output of this initial step (Tutorial) is an XML PopModelODE specification that, when plugged to a TrajectoryOutput object, generates a trajectory time-series. The next step is to introduce phylogenetic trees and compute their likelihood with respect to a population model trajectory using STreeLikelihoodODE class.

Simulation: Phylogenetic trees for our two-deme SIR model

Phylogenetic trees are key components to BEAST analyses. In the context of PhyDyn, we want to assess the likelihood of a structured tree given a population process trajectory or history. A BEAST analysis can use a fixed tree e.g. a maximum likelihood tree generated by another analysis, or work through a sequence of trees generated via sampling. For this example (and the following tutorial) we will work with a 100-tip Newick tree generated by sampling on a population process using MASTER. The ree is stored in sir2Master.nwk.

Our example Newick tree, sir2Master.nwk, was generated by executing the genSIR2Trees.R script:

  > Rscript genSIR2Trees.R <num-trees>

The following steps are executed until <num-trees> are generated:

• The script calls BEAST with MASTER file sir2Master.xml as input. MASTER is a BEAST package so, in order for this to work, MASTER should be loaded beforehand (standalone MASTER versions can be used by changing the master variable in the script).
• MASTER generates an annotated NEXUS file. The script processes the tree (trees with less than 200 tips are discarded, internal node annotations are removed) and writes it back in Newick format.
• Trees are named sir2Master_i.nwk, with i varying from 1 to <num-trees>.

The MASTER model describes the same two-model population model we are working with but specified in terms of Chemical Master Equations (CME). The specification instructs MASTER to start the simulation with initial values 999,1 and 0 for S, I0 and I1, respectively, and sample 200 demes (tree tips) at time 20. The resulting tree is stored in Nexus format by MASTER, processed by the script (trees with less than 200 tips are discarded, internal node annotations are removed) and written back in Newick format.

MASTER produces an annotated phylogenetic tree. In our context, it produces a structured phylogenetic tree with tips labeled with deme information (population names). For example, the following extract from sir2Master.nwk

((95_I0:0.009020825955,96_I0:0.009020825955):1.633546797,97_I1:1.642567623):1.326547235

shows that tips 95,96 and 97 correspond to demes I0,I0 and I1, respectively.

XML Generation

Copy the PhyDyn file sir2Trajectory0.xml (or your sir2TrajectoryFinal file created in the previous tutorial) and name it sir2Generate.xml.

  > cp sir2Trajectory0.xml sir2Generate.xml.

Open the new file in a text editor. Replace the TrajectoryOut component

<run spec="TrajectoryOut" model='@twodeme' file="sir2Final.csv" />

with

<run spec="XMLGenerator" xmlType="likelihood" model='@twodeme'
     outputFile="sir2Likelihood.xml" treeFile="sir2Master.nwk"
     adjustTipHeights="true" createDateTrait="false" createTypeTrait="false" />

Execute the file by running BEAST2/Phydyn:

  > beast sir2Generate.xml

Upon execution, PhyDyn's XMLGenerator creates the sir2Likelihood.xml file with the following components:

• id="twodeme-tree": A Tree XML element (spec TreeParser) with the contents of the Newick input tree file sirMaster.nwk verbatim.
• id="twodeme-data": A TaxonSet element with taxa ids extracted from the annotated tree tips.
• id="twodeme":The population model specification (PopModelODE) included in the original file.
• id="stlh": A StreeLikelihoodODE object with default parameters, and references to the TreeParser and PopModelODE objects.
• A LikelihoodOut XML element pointing to the likelihood element stlh:

     <run spec='LikelihoodOut' stlikelihood='@stlh' />

The new file is ready to go:

> beast sir2Likelihood.xml 
File: sir2Likelihood.xml seed: 127 threads: 1
model-name = twodeme;
.... model definition
logP=-1128.9024695637336