ipyrad-analysis toolkit: BUCKy

This notebook uses the Pedicularis example data set from the first empirical ipyrad tutorial. Here I show how to run BUCKy on a large set of loci parsed from the output file with the .alleles.loci ending. All code in this notebook is Python. You can simply follow along and execute this same code in a Jupyter notebook of your own.

Software requirements for this notebook

All required software can be installed through conda by running the commented out code below in a terminal.

## conda install -c BioBuilds mrbayes
## conda install -c ipyrad ipyrad
## conda install -c ipyrad bucky

## import Python libraries
import ipyrad.analysis as ipa
import ipyparallel as ipp

Cluster setup

To execute code in parallel we will use the ipyparallel Python library. A quick guide to starting a parallel cluster locally can be found here, and instructions for setting up a remote cluster on a HPC cluster is available here. In either case, this notebook assumes you have started an ipcluster instance that this notebook can find, which the cell below will test.

## look for running ipcluster instance, and create load-balancer
ipyclient = ipp.Client()
print("{} engines found".format(len(ipyclient)))

6 engines found

Create a bucky analysis object

The two required arguments are the name and data arguments. The data argument should be a .loci file or a .alleles.loci file. The name will be used to name output files, which will be written to {workdir}/{name}/{number}.nexus. Bucky doesn’t deal well with missing data, so loci will only be included if they contain data for all samples in the analysis. By default, all samples found in the loci file will be used, unless you enter a list of names (the samples argument) to subsample taxa, which we do here. It is best to select one individual per species or subspecies. You can set a number of additional parameters in the .params dictionary. Here I use the maxloci argument to limit the total number of loci so that the example analysis will finish faster. But in practice, BUCKy runs quite fast and I would typically just use all of your loci in a real analysis.

## make a list of sample names you wish to include in your BUCKy analysis
samples = [
    "29154_superba",
    "30686_cyathophylla",
    "41478_cyathophylloides",
    "33413_thamno",
    "30556_thamno",
    "35236_rex",
    "40578_rex",
    "38362_rex",
    "33588_przewalskii",
]

## initiate a bucky object
c = ipa.bucky(
    name="buckytest",
    data="analysis-ipyrad/pedic_outfiles/pedic.alleles.loci",
    workdir="analysis-bucky",
    samples=samples,
    minsnps=0,
    maxloci=100,
)

## print the params dictionary
c.params

bucky_alpha           [0.1, 1.0, 10.0]
bucky_nchains         4
bucky_niter           1000000
bucky_nreps           4
maxloci               100
mb_mcmc_burnin        100000
mb_mcmc_ngen          1000000
mb_mcmc_sample_freq   1000
minsnps               0
seed                  224443248

Write data to nexus files

As you will see below, one step of this analysis is to convert the data into nexus files with a mrbayes code block. Let’s run that step quickly here just to see what the converted files look like.

## This will write nexus files to {workdir}/{name}/[number].nex
c.write_nexus_files(force=True)

wrote 100 nexus files to ~/Documents/ipyrad/tests/analysis-bucky/buckytest

An example nexus file

## print an example nexus file
! cat analysis-bucky/buckytest/1.nex

#NEXUS
begin data;
dimensions ntax=9 nchar=64;
format datatype=dna interleave=yes gap=- missing=N;
matrix
30556_thamno            TCCTCGGCAGCCATTAAACCAGTGGAGTATGCACCATGTACCGATCCTGGATAATCAAAACTTG
40578_rex               TCCTCGGCAGCCATTAAACCAGTGGAGTATGCACCATGTACCGATCCTGGATAATCAAAACTTG
38362_rex               TCCTCGGCAGCCATTAAACCGGTGGAGTATGCACCATGTACCGATCCTGGATAATCAAAACTTG
29154_superba           TCCTCGGCAGCCATTAGACCGGTGGAGTATGCACCATGTACCGATCCTGGATAATCAAAACTCG
30686_cyathophylla      TCCTCGGCAGCCATTAGACCGGTGGAATATGCACCATGTACCGATCCTGGATAATCAAAACTCG
33413_thamno            TCCTCGGCAGCCATTAAACCGGTGGAGTATGCACCATGTACTGATCCTGGATAATCAAAACTTG
41478_cyathophylloides  TCCTCGGCAGCCATTAGACCAGTGGAGTATGCACCATGTACCGATCCTGGATAATCAAAACTCG
33588_przewalskii       TCCTCGGCAGCCATTAGACCGGTGGAGTGTGCACCATGCACCGATCCCGGATAATCAAAACTCG
35236_rex               TCCTCGGCAGCCATTAAACCAGTGGAGTATGCACCATGTACCGATCCTGGATAATCAAAACTTG

    ;

begin mrbayes;
set autoclose=yes nowarn=yes;
lset nst=6 rates=gamma;
mcmc ngen=1000000 samplefreq=1000 printfreq=1000000;
sump burnin=100000;
sumt burnin=100000;
end;

Complete a BUCKy analysis

There are four parts to a full BUCKy analysis. The first is converting the data into nexus files. The following are .run_mrbayes(), then .run_mbsum(), and finally .run_bucky(). Each uses the files produced by the previous function in order. You can use the force flag to overwrite existing files. An ipyclient should be provided to distribute the jobs in parallel. The parallelization is especially important for the mrbayes analyses, where more cores will lead to approximately linear speed improvements. An ipyrad.bucky analysis object will run all four steps sequentially by simply calling the .run() command. See the end of the notebook for results.

## run the complete analysis
c.run(force=True, ipyclient=ipyclient)

wrote 100 nexus files to ~/Documents/ipyrad/tests/analysis-bucky/buckytest
[####################] 100% [mb] infer gene-tree posteriors | 0:41:56 |
[####################] 100% [mbsum] sum replicate runs      | 0:00:02 |
[####################] 100% [bucky] infer CF posteriors     | 0:25:06 |

Alternatively, you can run each step separately

## (1) This will write nexus files to {workdir}/{name}/[number].nex
c.write_nexus_files(force=True)

wrote 100 nexus files to ~/Documents/ipyrad/tests/analysis-bucky/buckytest

## (2) distributes mrbayes jobs across the parallel client
c.run_mrbayes(force=True, ipyclient=ipyclient)

[####################] 100% [mb] infer gene-tree posteriors | 0:47:06 |

## (3) this step is fast, simply summing the gene-tree posteriors
c.run_mbsum(force=True, ipyclient=ipyclient)

[####################] 100% [mbsum] sum replicate runs      | 0:00:07 |

## (4) infer concordance factors with BUCKy. This will run in parallel
## for however many alpha values are in b.params.bucky_alpha list
b.run_bucky(force=True, ipyclient=ipyclient)

[####################] 100% [bucky] infer CF posteriors     | 1:35:16 |

Convenient access to results

View the results in the file [workdir]/[name]/CF-{alpha-value}.concordance. We haven’t yet developed any further ipyrad tools for parsing these results, but hope to do so in the future. The main results you are typically interested in are the Primary Concordance Tree and the Splits in the Primary Concordance Tree.

## print first 50 lines of a results files
! head -n 50 analysis-bucky/buckytest/CF-a1.0.concordance

translate
 1 30686_cyathophylla,
 2 33413_thamno,
 3 33588_przewalskii,
 4 29154_superba,
 5 41478_cyathophylloides,
 6 40578_rex,
 7 30556_thamno,
 8 38362_rex,
 9 35236_rex;

Population Tree:
(((((((1,4),5),3),2),8),6),7,9);

Primary Concordance Tree Topology:
(((((1,4),5),3),7),((2,8),6),9);

Population Tree, With Branch Lengths In Estimated Coalescent Units:
(((((((1:10.000,4:10.000):0.250,5:10.000):0.763,3:10.000):1.798,2:10.000):0.116,8:10.000):0.044,6:10.000):0.000,7:10.000,9:10.000);

Primary Concordance Tree with Sample Concordance Factors:
(((((1:1.000,4:1.000):0.465,5:1.000):0.578,3:1.000):0.754,7:1.000):0.255,((2:1.000,8:1.000):0.231,6:1.000):0.210,9:1.000);

Four-way partitions in the Population Tree: sample-wide CF, coalescent units and Ties(if present)
{1; 4|2,3,6,7,8,9; 5}   0.481, 0.250,
{1,4; 5|2,6,7,8,9; 3}   0.689, 0.763,
{1,4,5; 3|2; 6,7,8,9}   0.890, 1.798,
{1,2,3,4,5,8; 6|7; 9}   0.327, 0.000,
{1,3,4,5; 2|6,7,9; 8}   0.406, 0.116,
{1,2,3,4,5; 8|6; 7,9}   0.362, 0.044,

Splits in the Primary Concordance Tree: sample-wide and genome-wide mean CF (95% credibility), SD of mean sample-wide CF across runs
{1,3,4,5|2,6,7,8,9} 0.754(0.630,0.830) 0.746(0.596,0.862)       0.041
{1,4,5|2,3,6,7,8,9} 0.578(0.470,0.690) 0.572(0.425,0.718)       0.032
{1,4|2,3,5,6,7,8,9} 0.465(0.310,0.600) 0.461(0.276,0.631)       0.059
{1,3,4,5,7|2,6,8,9} 0.255(0.110,0.400) 0.252(0.092,0.423)       0.056
{1,3,4,5,6,7,9|2,8} 0.231(0.040,0.340) 0.229(0.032,0.381)       0.069
{1,3,4,5,7,9|2,6,8} 0.210(0.070,0.350) 0.208(0.059,0.374)       0.030

Splits NOT in the Primary Concordance Tree but with estimated CF > 0.050:
{1,5|2,3,4,6,7,8,9} 0.306(0.140,0.470) 0.304(0.122,0.503)       0.074
{1,2,3,4,5|6,7,8,9} 0.250(0.150,0.390) 0.248(0.121,0.410)       0.030
{1,2,3,4,5,8|6,7,9} 0.194(0.050,0.360) 0.192(0.027,0.389)       0.095
{1,2,3,4,5,7,8|6,9} 0.174(0.000,0.390) 0.173(0.000,0.409)       0.117
{1,2,3,4,5,6|7,8,9} 0.160(0.000,0.420) 0.158(0.000,0.441)       0.138
{1,2,5,6,7,8,9|3,4} 0.146(0.020,0.280) 0.145(0.009,0.300)       0.032
{1,2,3,4,5,6,9|7,8} 0.140(0.040,0.340) 0.139(0.029,0.341)       0.040
{1,3,4,5,9|2,6,7,8} 0.139(0.060,0.260) 0.138(0.038,0.283)       0.051
{1,2,3,4,5,9|6,7,8} 0.133(0.070,0.230) 0.132(0.045,0.255)       0.018
{1,2,3,4,5,8,9|6,7} 0.130(0.000,0.450) 0.130(0.000,0.467)       0.093