Tutorial 4
In this tutorial we delve into some of the more advanced aspects of inheritance tree generation including explicit inheritance relationship specification and time-reversed tree building. We do this by showing how MASTER can be used to simulate trees under the well-known structured coalescent model.
An extension of Kingman's coalescent, the structured coalescent gives the probability density of phylogenetic trees conditional on the number, location and time of the sampled lineage tips. Trees can be drawn from this probability distribution via a backward-time continuous time Markov process.
We set this up by considering a structured population of surviving
lineages, with members of an individual lineage within deme i labeled
. The reverse-time Markov process involves a pair of
reactions. The first describes a coalescent process between two
lineages in the same deme:
The coalescent rate 
is inversely proportional to the
total number of individuals (not lineages) contained in that deme and
is assumed here to be constant. The second reaction describes a
backward-time migration event where a lineage in deme i moves to deme
j.
The backward-time migration rate 
is directly related
to the forward-time migration rate between these two demes.
The MASTER specification of this model is almost identical from a
structural standpoint to that in the previous tutorial. We use exactly the same <run> element, for instance. The major
differences lie in the specification of the model itself and the way
in which the output is generated.
The Markovian process which simulates the tree operates in reverse time. This means that the tree is built starting at the leaves and finishing at the root. A slight complication arises because MASTER's default assignment of parent/child relationships between reactants and products only generates tree-like structures in the forward-time direction - from root to tips. Thus, this default assignment will not work here.
To get around this, we must explicitly specify which reactant is a parent of which child. This is done by appending strings of the form ":ID" where ID is an integer to each of the reactants and products in the reaction descriptor. Reactants and products having the same ID are assumed to have a parent/child relationship. In the coalescent process (the first reaction above), we achieve the desired relationship easily by appending ":1" to all population specifiers in each of the reactions. This has the effect of ensuring MASTER knows that the product lineage is the "child" of both reactant lineages.
<reactionGroup spec='ReactionGroup' reactionGroupName='Coalescence'>
<reaction spec='Reaction' rate="1.0">
2L[0]:1 -> L[0]:1
</reaction>
<reaction spec='Reaction' rate="1.0">
2L[1]:1 -> L[1]:1
</reaction>
</reactionGroup>No explicit inheritance relationship needs to be specified for the migration reactions, as the default assignment (individuals with like types share a lineage) produces the desired relationship in that case.
The initial state for this calculation is entirely made up of lineage seeds. Unlike in the previous tutorial, however, we will here make use of MASTER's ability to seed lineages at different times to allow us to simulate trees conditional on non-contemporaneously sampled tips. Additionally, we use <lineageSeedMultiple> elements to specify the addition of many lineages of the same population type and location simultaneously.
<initialState spec='InitState'>
<lineageSeedMultiple spec='MultipleIndividuals' copies="20" time="0.0">
<population spec='Population' type='@L' location="1"/>
</lineageSeedMultiple>
<lineageSeedMultiple spec='MultipleIndividuals' copies="20" time="0.2">
<population spec='Population' type='@L' location="0"/>
</lineageSeedMultiple>
</initialState>Be aware that the time attributes belong to the <lineageSeedMultiple> (or <lineageSeed>) elements and not the <population> elements.
In the previous tutorial, we used a <lineageEndCondition> element to cause the simulation to stop once the number of extant lineages reached zero. In the case of coalescent tree simulation, things stop becoming interesting when the MRCA (most recent common ancestor) is reached. We thus include an end condition which fires when the number of extant lineages reaches 1.
<lineageEndCondition spec='LineageEndCondition' nLineages="1"/>In this instance we use only the NEXUS output writer. An important addition is the use of the reverseTime="true" attribute, which specifies that the inheritance graph should be read in the opposite direction to that in which it was generated. If this attribute is missing or set to "false" the NEXUS generator will read the graph in a direction in which it is not tree-like. A NEXUS file will still be produced, but the Newick strings it contains will be of the "extended" variety used to describe phylogenetic networks, as described by Cardona, Rosselló and Valiente, BMC Bioinf. (2008).
<output spec='NexusOutput'
fileName='StructuredCoalescentTree_out.nexus'
reverseTime="true"/>The following is the completed XML input file, including the migration reactions omitted from the above description.
<beast version='2.0' namespace='master:master.model:master.steppers:master.conditions:master.postprocessors:master.outputs'>
<run spec='InheritanceTrajectory'
verbosity='2'>
<model spec='Model'>
<populationType spec='PopulationType' typeName='L' id='L' dim="2"/>
<reactionGroup spec='ReactionGroup' reactionGroupName='Coalescence'>
<reaction spec='Reaction' rate="1.0">
2L[0]:1 -> L[0]:1
</reaction>
<reaction spec='Reaction' rate="1.0">
2L[1]:1 -> L[1]:1
</reaction>
</reactionGroup>
<reactionGroup spec='ReactionGroup' reactionGroupName='Migration'>
<reaction spec='Reaction' rate="1.0">
L[0] -> L[1]
</reaction>
<reaction spec='Reaction' rate="1.0">
L[1] -> L[0]
</reaction>
</reactionGroup>
</model>
<initialState spec='InitState'>
<lineageSeedMultiple spec='MultipleIndividuals' copies="20" time="0.0">
<population spec='Population' type='@L' location="1"/>
</lineageSeedMultiple>
<lineageSeedMultiple spec='MultipleIndividuals' copies="20" time="0.2">
<population spec='Population' type='@L' location="0"/>
</lineageSeedMultiple>
</initialState>
<lineageEndCondition spec='LineageEndCondition' nLineages="1"/>
<output spec='NexusOutput' fileName='StructuredCoalescentTree_out.nexus' reverseTime="true"/>
</run>
</beast>The NEXUS file generated by MASTER can again be visualized using FigTree. The figure below, for instance, was generated by using the "appearance" tab to specify that lineages should be coloured by the "location" trait. Note the two sets of lineage introductions at t=0 and t=0.2.
