Species-tree & species-delimitation using bpp (BP&P)

The program bpp by Rannala & Yang (2010; 2015) is a powerful tool for inferring species tree parameters and testing species delimitation hypotheses. It is relatively easy to use, and best of all, it's quite fast, although not highly parallelizable. This notebook describes a streamlined approach to easily setup input files for testing different hypthotheses in bpp, and to do so in a clear programmatic way that makes it easy to perform many tests over many different parameter settings. We also show how to distribute many separate jobs to run in parallel on a cluster.

Notebook setup

This is a Jupyter notebook, a reproducible and executable document. The code in this notebook is Python (2.7), and should be executed in a jupyter-notebook like this one. Execute each cell in order to reproduce our entire analysis. The example data set used in this analysis is from the empirical example ipyrad tutorial.

Install required software

All software required for this notebook can be installed using conda.

In [1]:
## conda install -c ipyrad ipyrad
## conda install -c ipyrad bpp
## conda install -c eaton-lab toytree
In [2]:
import ipyrad.analysis as ipa         ## ipyrad analysis tools
import ipyparallel as ipp             ## parallelization
import pandas as pd                   ## DataFrames
import numpy as np                    ## data generation
import toytree                        ## tree plotting
import toyplot                        ## data plotting

Connect to an ipyparallel cluster

We will use the ipyparallel library to submit jobs to run in parallel on a cluster. You will need to have an ipcluster instance running in a separate terminal on your machine (or ideally, it is running on your HPC cluster). The code below simply connects to that cluster and prints how many CPUs are available for use.

In [3]:
## Connect to a running ipcluster instance
ipyclient = ipp.Client()#profile="alt")

## print information about our cluster
print "Connected to {} cores".format(len(ipyclient))
Connected to 10 cores

Analysis setup

You must define a tree with the "species" names in your analysis which will act either as a fixed-tree or as a guide-tree for your analysis. You must also define an IMAP dictionary which maps sample names to "species" names. And you can optionally also define a MINMAP dictionary which is used to filter RAD loci to include only loci that have at least N samples with data from each species in a locus.

In [16]:
## set the location of our input .loci file
locifile = "./analysis-ipyrad/pedicularis_outfiles/pedicularis.alleles.loci"

## a tree hypothesis (guidetree) (here based on tetrad results)
newick = "((((((rex, lip),rck),tha),cup),(cys,(cya, sup))),prz);"

## a dictionary mapping sample names to 'species' names
imap = {
    "prz": ["32082_przewalskii", "33588_przewalskii"],
    "cys": ["41478_cyathophylloides", "41954_cyathophylloides"],
    "cya": ["30686_cyathophylla"],
    "sup": ["29154_superba"],
    "cup": ["33413_thamno"],
    "tha": ["30556_thamno"],
    "rck": ["35236_rex"],
    "rex": ["35855_rex", "40578_rex"],
    "lip": ["39618_rex", "38362_rex"],  
    }

## optional: loci will be filtered if they do not have data for at
## least N samples/individuals in each species.
minmap = {
    "prz": 2,
    "cys": 2,
    "cya": 1,
    "sup": 1,
    "cup": 1,
    "tha": 1, 
    "rck": 1,
    "rex": 2,
    "lip": 2,
    }
In [17]:
## check your (starting) tree hypothesis
toytree.tree(newick).draw(width=250);
przcyscyasupcuptharckrexlip

The bpp Class object

To simplify the creation of input files for bpp analyses we've created a bpp job generator object that can be accessed from ipa.bpp(). Running bpp requires three input files (.ctl, .imap, and .seq) of which the .ctl file is the most important since it contains the parameters for a run and points to the location of the other two files. The ipa.bpp() object can be used to easily modify parameter settings for a run, to generate the input files, and if desired, to submit the bpp jobs to run on a cluster (your ipyclient cluster).

In [18]:
## create a bpp object to run algorithm 00
b = ipa.bpp(
    name="test",
    data=locifile,
    guidetree=newick, 
    imap=imap,
    minmap=minmap,   
    )
In [45]:
## set some optional params, leaving others at their defaults
b.params.burnin = 2000
b.params.nsample = 20000
b.params.sampfreq = 2
b.params.seed = 33333
b.params.tauprior = (2, 200, 1)

## print params
b.params
Out[45]:
binary          bpp                 
burnin          2000                
cleandata       0                   
delimit_alg     (0, 5)              
finetune        (0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01)
infer_delimit   0                   
infer_sptree    0                   
nsample         20000               
sampfreq        2                   
seed            33333               
tauprior        (2, 200, 1)         
thetaprior      (2, 2000)           
usedata         1                   
In [37]:
## set some optional filters leaving others at their defaults
b.filters.maxloci=100
b.filters.minsnps=0

## print filters
b.filters
Out[37]:
maxloci   100                 
minmap    {'cys': 4, 'rex': 4, 'cup': 2, 'rck': 2, 'cya': 2, 'lip': 4, 'sup': 2, 'tha': 2, 'prz': 4}
minsnps   0                   

Generating files &/or submitting jobs

When you create a bpp object you save it with a variable name (in this example named test) which is simply the name you will use to reference the object. You must also provide a name (prefix for output files) that will be used by either of the two main functions of the bpp object: write_bpp_files() or run(). Both functions make it easy to sample different distributions of loci to include in different replicate bpp analyses. Each replicate will start from a different random seed after the initial seed. If you used a maxloci argument to limit the number of loci that will used in the analysis then you can also use the randomize_order argument to select a different random number of N loci in each rep, otherwise just the first N loci that pass filtering will be sampled.

write_bpp_files()

This writes the .ctl, .seq, and .imap files for the specified run.

In [49]:
## write files (filtered based on test.filters and test.params)
b.write_bpp_files()
input files created for job test (100 loci)

run()

This writes the files for each job and submits the bpp jobs to run on the cluster designated by the ipyclient object. You can efficiently submit many replicate jobs in this way.

In [50]:
## or, write files and submit job to run on ipyclient 
b.run(
    nreps=2, 
    ipyclient=ipyclient, 
    randomize_order=True,
    )
submitted 2 bpp jobs [test_r1] (100 loci)

Accessing job results

When you submit jobs the results files will be stored in the bpp objects .files attribute. In addition, the asychronous result objects (a representation of the running job) from each submitted job will be accessible from the .asyncs attribute of the bpp object. You can view these objects to see if your job has finished or use them to trace errors if an error arises, as we do below.

In [59]:
## files associated with the bpp object (includes replicates)
! ls ~/local/src/ipyrad/tests/analysis-bpp/test*.txt
/home/deren/local/src/ipyrad/tests/analysis-bpp/test.ctl.txt
/home/deren/local/src/ipyrad/tests/analysis-bpp/test.imapfile.txt
/home/deren/local/src/ipyrad/tests/analysis-bpp/test_r0.ctl.txt
/home/deren/local/src/ipyrad/tests/analysis-bpp/test_r0.imapfile.txt
/home/deren/local/src/ipyrad/tests/analysis-bpp/test_r0.mcmc.txt
/home/deren/local/src/ipyrad/tests/analysis-bpp/test_r0.out.txt
/home/deren/local/src/ipyrad/tests/analysis-bpp/test_r0.seqfile.txt
/home/deren/local/src/ipyrad/tests/analysis-bpp/test_r1.ctl.txt
/home/deren/local/src/ipyrad/tests/analysis-bpp/test_r1.imapfile.txt
/home/deren/local/src/ipyrad/tests/analysis-bpp/test_r1.mcmc.txt
/home/deren/local/src/ipyrad/tests/analysis-bpp/test_r1.out.txt
/home/deren/local/src/ipyrad/tests/analysis-bpp/test_r1.seqfile.txt
/home/deren/local/src/ipyrad/tests/analysis-bpp/test.seqfile.txt
In [60]:
## check async objects from the bpp object
for job in b.asyncs:
    if job.ready():
        print 'job finished'
    else:
        print 'job running'
job finished
job finished
In [14]:
## uncomment and run this to block until all running jobs are finished
#ipyclient.wait()

Examples


Algorithm 00 - fixed tree parameter inference

The 00 algorithm means 'infer_sptree=0' and 'infer_delimit=0', thus the tree that you enter will be treated as the fixed species tree and the analysis will infer parameters for the tree under the multispecies coalescent model. This will yield values of $\Theta$ for each branch of the tree, and divergence times ($\tau$) for each split in the tree.

In [38]:
## create two new copies of the bpp object above, but with new names
A00 = b.copy(name="A00")
A00n = b.copy(name="A00-nodata")

## set params on the 'nodata' object to not use seq data
A00n.params.usedata = 0
In [39]:
## submit a few replicate jobs from different random subsets of loci 
A00.run(nreps=3, randomize_order=True, ipyclient=ipyclient, force=True)
A00n.run(nreps=1, randomize_order=True, ipyclient=ipyclient, force=True)
submitted 3 bpp jobs [A00_r2] (100 loci)
submitted 1 bpp jobs [A00-nodata_r0] (100 loci)
In [16]:
## wait for jobs to finish
ipyclient.wait()

Summarize results tables for algorithm 00

Different bpp algorithms produce different types of results files. For algorithm 00 the mcmc results file is simply a table of $\Theta$ and $\tau$ values so we can just parse it as a CSV file to summarize results. The same results will be available in the .out.txt file, but I find that parsing the results this way is a bit easier to view.

The 'no-data' results

We expect that when we use no data the priors will be returned. In this case we see no significant differences in theta values among samples, which is good. The tau values will be distributed according to the prior on the root divergence, with the other nodes distributed according to a dirichlet process. It will be of interest to compare these divergence time estimates to the ones produced by the analysis with sequence data, to assess whether the sequence data causing the divergence times estimates to diverge significantly from their priors.

In [41]:
## set options to print prettier tables
## (this suppresses scientific notation)
pd.set_option('display.float_format', lambda x: '%.4f' % x)
In [42]:
## parse the mcmc table
table1 = pd.read_csv(A00n.files.mcmcfiles[0], sep="\t", index_col=0)
table1.describe().T
Out[42]:
count mean std min 25% 50% 75% max
theta_1cup 20000.0000 0.0011 0.0007 0.0001 0.0006 0.0009 0.0014 0.0064
theta_2cya 20000.0000 0.0010 0.0007 0.0001 0.0005 0.0008 0.0013 0.0075
theta_3cys 20000.0000 0.0013 0.0008 0.0001 0.0008 0.0011 0.0017 0.0067
theta_4lip 20000.0000 0.0011 0.0006 0.0001 0.0006 0.0009 0.0014 0.0050
theta_5prz 20000.0000 0.0003 0.0004 0.0000 0.0001 0.0001 0.0003 0.0042
theta_6rck 20000.0000 0.0011 0.0007 0.0001 0.0005 0.0009 0.0014 0.0067
theta_7rex 20000.0000 0.0012 0.0008 0.0001 0.0006 0.0010 0.0015 0.0081
theta_8sup 20000.0000 0.0011 0.0007 0.0001 0.0005 0.0009 0.0014 0.0055
theta_9tha 20000.0000 0.0010 0.0007 0.0001 0.0005 0.0008 0.0014 0.0054
theta_10rexliprckthacupcyscyasupprz 20000.0000 0.0009 0.0007 0.0000 0.0004 0.0008 0.0012 0.0053
theta_11rexliprckthacupcyscyasup 20000.0000 0.0011 0.0007 0.0000 0.0005 0.0009 0.0014 0.0051
theta_12rexliprckthacup 20000.0000 0.0010 0.0007 0.0001 0.0005 0.0008 0.0013 0.0072
theta_13rexliprcktha 20000.0000 0.0010 0.0008 0.0000 0.0005 0.0009 0.0014 0.0066
theta_14rexliprck 20000.0000 0.0011 0.0008 0.0000 0.0005 0.0009 0.0014 0.0064
theta_15rexlip 20000.0000 0.0010 0.0007 0.0001 0.0005 0.0009 0.0014 0.0058
theta_16cyscyasup 20000.0000 0.0009 0.0007 0.0000 0.0003 0.0007 0.0012 0.0066
theta_17cyasup 20000.0000 0.0009 0.0007 0.0000 0.0004 0.0008 0.0013 0.0049
tau_10rexliprckthacupcyscyasupprz 20000.0000 0.0100 0.0073 0.0010 0.0049 0.0082 0.0128 0.0597
tau_11rexliprckthacupcyscyasup 20000.0000 0.0088 0.0061 0.0006 0.0043 0.0073 0.0113 0.0494
tau_12rexliprckthacup 20000.0000 0.0075 0.0055 0.0005 0.0033 0.0062 0.0099 0.0417
tau_13rexliprcktha 20000.0000 0.0063 0.0047 0.0005 0.0028 0.0052 0.0083 0.0377
tau_14rexliprck 20000.0000 0.0050 0.0040 0.0005 0.0022 0.0039 0.0063 0.0299
tau_15rexlip 20000.0000 0.0032 0.0031 0.0001 0.0011 0.0022 0.0042 0.0211
tau_16cyscyasup 20000.0000 0.0065 0.0054 0.0002 0.0024 0.0049 0.0092 0.0461
tau_17cyasup 20000.0000 0.0044 0.0047 0.0001 0.0009 0.0025 0.0064 0.0329

Results with data

You can see here that the theta and tau values are quite different from the priors, which is good and tells us that our sequence data is informative. If you check the .out file produced by bpp you can also see "effective-sample-size (ess)" which is indicative of whether or not we've run the mcmc long enough.

In [43]:
## parse the mcmc table
table2_r0 = pd.read_csv(A00.files.mcmcfiles[0], sep="\t", index_col=0)
table2_r1 = pd.read_csv(A00.files.mcmcfiles[1], sep="\t", index_col=0)
table2_r1.describe().T
Out[43]:
count mean std min 25% 50% 75% max
theta_1cup 20000.0000 0.0015 0.0005 0.0004 0.0011 0.0014 0.0017 0.0040
theta_2cya 20000.0000 0.0007 0.0003 0.0001 0.0005 0.0007 0.0009 0.0024
theta_3cys 20000.0000 0.0007 0.0002 0.0004 0.0006 0.0007 0.0008 0.0020
theta_4lip 20000.0000 0.0009 0.0002 0.0003 0.0007 0.0008 0.0010 0.0021
theta_5prz 20000.0000 0.0035 0.0005 0.0021 0.0031 0.0035 0.0038 0.0059
theta_6rck 20000.0000 0.0011 0.0003 0.0004 0.0008 0.0010 0.0013 0.0029
theta_7rex 20000.0000 0.0023 0.0006 0.0009 0.0019 0.0022 0.0026 0.0055
theta_8sup 20000.0000 0.0010 0.0004 0.0002 0.0007 0.0009 0.0012 0.0031
theta_9tha 20000.0000 0.0016 0.0005 0.0005 0.0012 0.0015 0.0019 0.0046
theta_10rexliprckthacupcyscyasupprz 20000.0000 0.0061 0.0026 0.0001 0.0046 0.0065 0.0080 0.0140
theta_11rexliprckthacupcyscyasup 20000.0000 0.0023 0.0007 0.0005 0.0017 0.0022 0.0027 0.0059
theta_12rexliprckthacup 20000.0000 0.0039 0.0012 0.0009 0.0030 0.0038 0.0047 0.0097
theta_13rexliprcktha 20000.0000 0.0040 0.0018 0.0002 0.0027 0.0041 0.0053 0.0112
theta_14rexliprck 20000.0000 0.0015 0.0009 0.0001 0.0008 0.0013 0.0020 0.0067
theta_15rexlip 20000.0000 0.0017 0.0010 0.0001 0.0009 0.0015 0.0023 0.0084
theta_16cyscyasup 20000.0000 0.0024 0.0014 0.0003 0.0014 0.0022 0.0031 0.0125
theta_17cyasup 20000.0000 0.0032 0.0010 0.0006 0.0025 0.0031 0.0038 0.0082
tau_10rexliprckthacupcyscyasupprz 20000.0000 0.0082 0.0019 0.0049 0.0067 0.0077 0.0091 0.0158
tau_11rexliprckthacupcyscyasup 20000.0000 0.0040 0.0005 0.0026 0.0037 0.0040 0.0043 0.0060
tau_12rexliprckthacup 20000.0000 0.0017 0.0005 0.0008 0.0014 0.0017 0.0021 0.0039
tau_13rexliprcktha 20000.0000 0.0009 0.0002 0.0004 0.0007 0.0009 0.0010 0.0017
tau_14rexliprck 20000.0000 0.0009 0.0002 0.0004 0.0007 0.0008 0.0010 0.0017
tau_15rexlip 20000.0000 0.0008 0.0002 0.0003 0.0006 0.0008 0.0009 0.0016
tau_16cyscyasup 20000.0000 0.0029 0.0006 0.0008 0.0026 0.0030 0.0033 0.0048
tau_17cyasup 20000.0000 0.0009 0.0004 0.0002 0.0006 0.0008 0.0011 0.0031
lnL 20000.0000 -11881.7791 21.4493 -11968.4540 -11895.5505 -11880.9480 -11866.9993 -11799.6660

Draw some histograms

Below we draw the posterior parameter distributions for the analyses with and without sequence data. In a real analysis you will probably want to run the analyses for longer, but you can see here that the results depart significantly from the prior expectation.

In [44]:
## draw barplots
canvas = toyplot.Canvas(width=800, height=800)

for aidx in range(table1.shape[1]):
    ## get data
    param = table1.columns[aidx]
    hist1 = np.histogram(table1[param], density=True)
    hist2 = np.histogram(table2_r0[param], density=True)
    hist3 = np.histogram(table2_r1[param], density=True)

    ## build plot
    nx = int(np.ceil(np.sqrt(table1.shape[1])))
    axes = canvas.cartesian(grid=(nx, nx, aidx), margin=25)
    axes.bars(hist1, opacity=0.75)
    axes.bars(hist2, opacity=0.75)
    axes.bars(hist3, opacity=0.75)

    ## style axes
    axes.label.text = param[:15]
    axes.label.style["font-size"] = "12px"
    axes.x.ticks.labels.show = False
    axes.y.ticks.labels.show = False
    
## save plot to pdf
import toyplot.pdf
toyplot.pdf.render(canvas, "analysis-bpp/bpp-A00.pdf")
canvas
Out[44]:
theta_1cuptheta_2cyatheta_3cystheta_4liptheta_5prztheta_6rcktheta_7rextheta_8suptheta_9thatheta_10rexliprtheta_11rexliprtheta_12rexliprtheta_13rexliprtheta_14rexliprtheta_15rexliptheta_16cyscyastheta_17cyasuptau_10rexliprcktau_11rexliprcktau_12rexliprcktau_13rexliprcktau_14rexliprcktau_15rexliptau_16cyscyasuptau_17cyasup

Algorithm 10 - species tree inference

The algorithm 10 aims to infer the correct species tree from the data by implemented a tree search method, thus the input tree is treated only as a starting tree. Based on the results above I'm also going to increase the priors on the divergence times (tau) by about 10X.

In [46]:
## create new bpp objects
A10 = A00.copy("A10")
A10n = A00.copy("A10-nodata")

## set new params
A10.params.infer_sptree = 1
A10.params.infer_delimit = 0
A10.params.tauprior = (2, 200, 1)
A10n.params.infer_sptree = 1
A10n.params.infer_delimit = 0
A10n.params.tauprior = (2, 200, 1)

## also set no-data on the 'n' copy
A10n.params.usedata = 0
In [47]:
## submit job reps to the cluster
A10.run(nreps=2, randomize_order=True, ipyclient=ipyclient)
A10n.run(nreps=1, randomize_order=True, ipyclient=ipyclient)
submitted 2 bpp jobs [A10_r1] (100 loci)
submitted 1 bpp jobs [A10-nodata_r0] (100 loci)
In [48]:
## block until finished
ipyclient.wait()
Out[48]:
True

Plot the distribution of species trees from algorithm 10

Here we use the multitree() object of toytree to plot the posterior distribution of trees. The code for drawing 'cloudtrees' here is still in development and a littl clunky, but will improve in the future. From the example here (which has not been run long enough, btw), you can see clearly that the two replicates yield quite different trees. In this case the replicates are starting from different random seeds, but probably more importantly they are sampling a different random subset of 100 RAD loci that met our filtering requirements (e.g., minsamples, minSNPs).

In [52]:
## load trees slicing out every Nth: (start:end:by)
trees1 = toytree.multitree(
    newick=A10.files.mcmcfiles[0],
    treeslice=(500, 5000, 30),
    )

## build canvas
canvas = toyplot.Canvas(width=800, height=400)
axes1 = canvas.cartesian(grid=(1, 2, 0), padding=25)
axes2 = canvas.cartesian(grid=(1, 2, 1), padding=25)

## draw cloud tree
trees1.draw_cloudtree(
    axes=axes1,
    use_edge_lengths=True,
    orient='right',
    edge_style={"opacity": 0.025},
);

## draw consensus (doesn't save edge lengths currently)
cons = trees1.get_consensus_tree()
cons.root(wildcard="prz")
cons.draw(
    axes=axes2,
    node_labels=cons.get_node_values("support"),
);

## style axes
axes1.y.show = False
axes2.y.show = False