{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"# ipyrad-analysis toolkit: bpp (alg 00)\n",
"\n",
"The ipyrad-analysis `bpp` tool w/ algorithm 00 can be used infer parameters on a fixed species tree. Please see the [full BPP documentation](https://github.com/bpp/bpp/blob/master/bppDOC.pdf) for details about BPP analyses. The ipyrad-analysis wrapper is intended to make it easy to auomate many BPP analyses on large RAD-seq datasets. \n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Load packages"
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [],
"source": [
"# conda install ipyrad toytree ipcoal -c conda-forge -c bioconda\n",
"# conda install bpp -c eaton-lab"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [],
"source": [
"import ipyrad.analysis as ipa\n",
"import pandas as pd\n",
"import numpy as np\n",
"import toytree\n",
"import toyplot\n",
"import ipcoal"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Simulate data under a specified model\n",
"We are interested in estimating the Ne and divergence time values on this tree."
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"
"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"# generate a tree with divergence times in generations\n",
"TREE = toytree.rtree.unittree(ntips=5, treeheight=3e6, seed=123)\n",
"\n",
"# set larger Ne on one lineage than the other\n",
"TREE = TREE.set_node_values(\"Ne\", {i: 5e5 for i in (0, 1, 2, 5, 6)}, default=2e5)\n",
"\n",
"# draw in 'p' style to show Ne\n",
"TREE.draw(ts='p', node_hover=True);"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [],
"source": [
"# create a simulation model for this tree/network with 3 diploids per spp.\n",
"model = ipcoal.Model(tree=TREE, nsamples=6, mut=1e-8)"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"wrote 1000 loci to /tmp/test-bpp.seqs.hdf5\n"
]
}
],
"source": [
"# simulate N loci\n",
"model.sim_loci(nloci=1000, nsites=150)\n",
"\n",
"# write result to a database file\n",
"model.write_loci_to_hdf5(name=\"test-bpp\", outdir=\"/tmp\", diploid=True)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### *True* measurements that will inform our priors\n",
"In BPP we can set a prior on the percent divergence $\\tau$ (pronounced Tau) between the root and tips of the tree. When combined with an estimate of the per-site-per-generation mutation rate ($\\mu$) and generation time (g) this can be converted into units of years. In our `ipcoal` simulation we actually have access to the root ancestral sequence as well as the tip sequences, so we can measure the true percent sequence divergence."
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"0.03360666666666667"
]
},
"execution_count": 6,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# concatenate seqs into single alignment\n",
"concat = np.concatenate(model.seqs, axis=1)\n",
"\n",
"# get ancestral root node sequence\n",
"rootseq = np.concatenate(model.ancestral_seq)\n",
"\n",
"# what is sequence divergence relative to root?\n",
"np.sum(rootseq != concat[0]) / rootseq.size"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"In addition to setting a prior on Tau we must also set a prior on $\\theta$ (pronounced theta), the population mutation rate ($\\theta$=4Ne$\\mu$). Here we know the Ne and $\\mu$ used in the simulation scenario, so again we can calculate it exactly. For real data we will not know these prior values precisely, and so we will want to put a wide prior to incorporate our uncertainty. "
]
},
{
"cell_type": "code",
"execution_count": 13,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"0.008"
]
},
"execution_count": 13,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# theta in simulation\n",
"4 * TREE.treenode.Ne * 1e-8"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### BPP analysis setup"
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {},
"outputs": [],
"source": [
"SEQS = \"/tmp/test-bpp.seqs.hdf5\""
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {},
"outputs": [],
"source": [
"IMAP = {\n",
" \"r0\": [\"r0-0\", \"r0-1\", \"r0-2\"],\n",
" \"r1\": [\"r1-0\", \"r1-1\", \"r1-2\"],\n",
" \"r2\": [\"r2-0\", \"r2-1\", \"r2-2\"],\n",
" \"r3\": [\"r3-0\", \"r3-1\", \"r3-2\"],\n",
" \"r4\": [\"r4-0\", \"r4-1\", \"r4-2\"],\n",
"}"
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {},
"outputs": [],
"source": [
"# create a bpp object to run algorithm 00\n",
"tool = ipa.bpp(\n",
" name=\"test-bpp\",\n",
" workdir=\"/tmp\",\n",
" data=SEQS,\n",
" guidetree=TREE,\n",
" imap=IMAP,\n",
" infer_sptree=0,\n",
" infer_delimit=0,\n",
" maxloci=1000,\n",
" burnin=2e3,\n",
" nsample=2e4,\n",
" sampfreq=5,\n",
")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### The parameters of the object"
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"{'binary': '/tmp/bpp-4.1.4-linux-x86_64/bin/bpp',\n",
" 'infer_sptree': 0,\n",
" 'infer_delimit': 0,\n",
" 'infer_delimit_args': (0, 2),\n",
" 'speciesmodelprior': 1,\n",
" 'seed': 12345,\n",
" 'burnin': 2000,\n",
" 'nsample': 20000,\n",
" 'sampfreq': 5,\n",
" 'thetaprior': (3, 0.002),\n",
" 'tauprior': (3, 0.002),\n",
" 'phiprior': (1, 1),\n",
" 'usedata': 1,\n",
" 'cleandata': 0,\n",
" 'finetune': (0.01, 0.02, 0.03, 0.04, 0.05, 0.01, 0.01),\n",
" 'copied': False}"
]
},
"execution_count": 11,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"tool.kwargs"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Set priors and plot priors of transformed values\n",
"\n",
"Here we will enter our expectations on the generation time and mutation rate of our simulated organisms so that we can toggle the theta and tau priors until the resulting prior distributions on Ne and T seem reasonable. \n",
"\n",
"In our simulation the root divergence time was 3M generations. Assuming a generation time of 1 and prior on root tau of ~3% (as we measured above) yields a prior on the root divergence time with 95% CI from about 0.5-6.5 Ma. So this is a pretty good but not highly informative prior. \n",
"\n",
"Similarly, our simulation used Ne=200K or 500K on different lineages. If we assume a wide prior on the per-site per-generation mutation rate ($\\mu$) with highest density of 1e-8 (as used in our simulations) and a prior on $\\theta$ with highest density at about 0.01 (we calculated 0.008 from the simulated data above) yields a 95% CI on Ne that is again accurate but not highly informative. \n",
"\n",
"These priors incorporate *estimated* knowledge about our organisms, by describing gentime and mutrate as distributions, and puts informative but wide priors on the parameters we hope to infer (theta and tau). \n"
]
},
{
"cell_type": "code",
"execution_count": 14,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"
"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"# toggle these prior settings \n",
"tool.kwargs[\"thetaprior\"] = (2.5, 0.004)\n",
"tool.kwargs[\"tauprior\"] = (3, 0.01)\n",
"\n",
"# draw the distributions when using our assumptions for transforming\n",
"tool.draw_priors(\n",
" gentime_min=0.99, gentime_max=1.01, \n",
" mutrate_min=5e-9, mutrate_max=2e-8,\n",
");"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Run inference\n",
"Distribute jobs in parallel and run replicate jobs that start MCMC from different random seeds. "
]
},
{
"cell_type": "code",
"execution_count": 15,
"metadata": {
"scrolled": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Parallel connection | latituba: 8 cores\n",
"[locus filter] full data: 1000\n",
"[locus filter] post filter: 1000\n",
"[ipa bpp] bpp v4.1.4\n",
"[ipa.bpp] distributed 2 bpp jobs (name=test-bpp, nloci=1000)\n",
"[####################] 100% 3:36:37 | progress on all reps \n",
"Parallel connection closed.\n"
]
}
],
"source": [
"tool.run(auto=True, nreps=2, force=True)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Re-load results\n",
"You can load your results anytime later after your runs have finished by providing the name and workdir for the job to ipa.bpp. This allows you to mess around with plotting and comparing distributions later. You also have the option here to combine the posteriors from multiple replicate runs (that used the same data and priors), or to analyze each separately. "
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [],
"source": [
"import ipyrad.analysis as ipa"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"[ipa.bpp] found 2 existing result files\n",
"[ipa.bpp] summarizing algorithm '00' results\n",
"[ipa.bpp] combining mcmc files\n"
]
}
],
"source": [
"# reload results by creating object with a given name and workdir\n",
"tool = ipa.bpp(name=\"test-bpp\", workdir=\"/tmp\")\n",
"\n",
"# call summarize on the tool\n",
"table, mcmc = tool.summarize_results(\"00\", individual_results=False)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### The \"00\" tables result\n",
"\n",
"The main result of the \"00\" algorithm is a dataframe with the inferred parameters of the multispecies coalescent model. "
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {
"scrolled": false
},
"outputs": [
{
"data": {
"text/html": [
"
\n",
"\n",
"
\n",
" \n",
"
\n",
"
\n",
"
theta_1r0
\n",
"
theta_2r1
\n",
"
theta_3r2
\n",
"
theta_4r3
\n",
"
theta_5r4
\n",
"
theta_6r4r3r2r1r0
\n",
"
theta_7r4r3
\n",
"
theta_8r2r1r0
\n",
"
theta_9r1r0
\n",
"
tau_6r4r3r2r1r0
\n",
"
tau_7r4r3
\n",
"
tau_8r2r1r0
\n",
"
tau_9r1r0
\n",
"
lnL
\n",
"
\n",
" \n",
" \n",
"
\n",
"
mean
\n",
"
0.027995
\n",
"
0.028212
\n",
"
0.025590
\n",
"
0.009736
\n",
"
0.009295
\n",
"
0.007839
\n",
"
0.005385
\n",
"
0.022274
\n",
"
0.017841
\n",
"
0.034456
\n",
"
0.024037
\n",
"
0.026410
\n",
"
0.016234
\n",
"
-418047.448228
\n",
"
\n",
"
\n",
"
median
\n",
"
0.027980
\n",
"
0.028194
\n",
"
0.025576
\n",
"
0.009733
\n",
"
0.009287
\n",
"
0.007836
\n",
"
0.005425
\n",
"
0.022171
\n",
"
0.017826
\n",
"
0.034457
\n",
"
0.024008
\n",
"
0.026409
\n",
"
0.016230
\n",
"
-418046.894000
\n",
"
\n",
"
\n",
"
S.D
\n",
"
0.000885
\n",
"
0.000904
\n",
"
0.000717
\n",
"
0.000304
\n",
"
0.000296
\n",
"
0.000704
\n",
"
0.001142
\n",
"
0.002391
\n",
"
0.001121
\n",
"
0.000394
\n",
"
0.000616
\n",
"
0.000444
\n",
"
0.000313
\n",
"
79.616227
\n",
"
\n",
"
\n",
"
min
\n",
"
0.024424
\n",
"
0.024989
\n",
"
0.022996
\n",
"
0.008656
\n",
"
0.008194
\n",
"
0.005509
\n",
"
0.002602
\n",
"
0.014451
\n",
"
0.013855
\n",
"
0.032841
\n",
"
0.021608
\n",
"
0.024757
\n",
"
0.015049
\n",
"
-418388.687000
\n",
"
\n",
"
\n",
"
max
\n",
"
0.032145
\n",
"
0.032099
\n",
"
0.028488
\n",
"
0.011026
\n",
"
0.010516
\n",
"
0.010827
\n",
"
0.009983
\n",
"
0.032692
\n",
"
0.023367
\n",
"
0.036106
\n",
"
0.026211
\n",
"
0.028206
\n",
"
0.017480
\n",
"
-417729.014000
\n",
"
\n",
"
\n",
"
2.5%
\n",
"
0.026300
\n",
"
0.026481
\n",
"
0.024243
\n",
"
0.009154
\n",
"
0.008728
\n",
"
0.006511
\n",
"
0.003146
\n",
"
0.017910
\n",
"
0.015693
\n",
"
0.033689
\n",
"
0.022892
\n",
"
0.025541
\n",
"
0.015626
\n",
"
-418202.863000
\n",
"
\n",
"
\n",
"
97.5%
\n",
"
0.029785
\n",
"
0.030023
\n",
"
0.027033
\n",
"
0.010342
\n",
"
0.009899
\n",
"
0.009229
\n",
"
0.007568
\n",
"
0.027284
\n",
"
0.020108
\n",
"
0.035223
\n",
"
0.025280
\n",
"
0.027283
\n",
"
0.016861
\n",
"
-417892.551000
\n",
"
\n",
"
\n",
"
2.5%HPD
\n",
"
0.026283
\n",
"
0.026471
\n",
"
0.024219
\n",
"
0.009152
\n",
"
0.008691
\n",
"
0.006477
\n",
"
0.002991
\n",
"
0.017773
\n",
"
0.015733
\n",
"
0.033673
\n",
"
0.022859
\n",
"
0.025528
\n",
"
0.015623
\n",
"
-418201.078000
\n",
"
\n",
"
\n",
"
97.5%HPD
\n",
"
0.029757
\n",
"
0.030004
\n",
"
0.027007
\n",
"
0.010338
\n",
"
0.009860
\n",
"
0.009190
\n",
"
0.007314
\n",
"
0.027099
\n",
"
0.020134
\n",
"
0.035200
\n",
"
0.025239
\n",
"
0.027268
\n",
"
0.016855
\n",
"
-417891.188000
\n",
"
\n",
"
\n",
"
ESS*
\n",
"
12262.897304
\n",
"
11997.844527
\n",
"
13891.059855
\n",
"
6974.445253
\n",
"
6525.824045
\n",
"
436.125051
\n",
"
123.289724
\n",
"
1897.302231
\n",
"
3587.614817
\n",
"
579.708878
\n",
"
147.436912
\n",
"
2641.537188
\n",
"
4378.739376
\n",
"
2087.556267
\n",
"
\n",
"
\n",
"
Eff*
\n",
"
0.613145
\n",
"
0.599892
\n",
"
0.694553
\n",
"
0.348722
\n",
"
0.326291
\n",
"
0.021806
\n",
"
0.006164
\n",
"
0.094865
\n",
"
0.179381
\n",
"
0.028985
\n",
"
0.007372
\n",
"
0.132077
\n",
"
0.218937
\n",
"
0.104378
\n",
"
\n",
" \n",
"
\n",
"
"
],
"text/plain": [
" theta_1r0 theta_2r1 theta_3r2 theta_4r3 theta_5r4 \\\n",
"mean 0.027995 0.028212 0.025590 0.009736 0.009295 \n",
"median 0.027980 0.028194 0.025576 0.009733 0.009287 \n",
"S.D 0.000885 0.000904 0.000717 0.000304 0.000296 \n",
"min 0.024424 0.024989 0.022996 0.008656 0.008194 \n",
"max 0.032145 0.032099 0.028488 0.011026 0.010516 \n",
"2.5% 0.026300 0.026481 0.024243 0.009154 0.008728 \n",
"97.5% 0.029785 0.030023 0.027033 0.010342 0.009899 \n",
"2.5%HPD 0.026283 0.026471 0.024219 0.009152 0.008691 \n",
"97.5%HPD 0.029757 0.030004 0.027007 0.010338 0.009860 \n",
"ESS* 12262.897304 11997.844527 13891.059855 6974.445253 6525.824045 \n",
"Eff* 0.613145 0.599892 0.694553 0.348722 0.326291 \n",
"\n",
" theta_6r4r3r2r1r0 theta_7r4r3 theta_8r2r1r0 theta_9r1r0 \\\n",
"mean 0.007839 0.005385 0.022274 0.017841 \n",
"median 0.007836 0.005425 0.022171 0.017826 \n",
"S.D 0.000704 0.001142 0.002391 0.001121 \n",
"min 0.005509 0.002602 0.014451 0.013855 \n",
"max 0.010827 0.009983 0.032692 0.023367 \n",
"2.5% 0.006511 0.003146 0.017910 0.015693 \n",
"97.5% 0.009229 0.007568 0.027284 0.020108 \n",
"2.5%HPD 0.006477 0.002991 0.017773 0.015733 \n",
"97.5%HPD 0.009190 0.007314 0.027099 0.020134 \n",
"ESS* 436.125051 123.289724 1897.302231 3587.614817 \n",
"Eff* 0.021806 0.006164 0.094865 0.179381 \n",
"\n",
" tau_6r4r3r2r1r0 tau_7r4r3 tau_8r2r1r0 tau_9r1r0 lnL \n",
"mean 0.034456 0.024037 0.026410 0.016234 -418047.448228 \n",
"median 0.034457 0.024008 0.026409 0.016230 -418046.894000 \n",
"S.D 0.000394 0.000616 0.000444 0.000313 79.616227 \n",
"min 0.032841 0.021608 0.024757 0.015049 -418388.687000 \n",
"max 0.036106 0.026211 0.028206 0.017480 -417729.014000 \n",
"2.5% 0.033689 0.022892 0.025541 0.015626 -418202.863000 \n",
"97.5% 0.035223 0.025280 0.027283 0.016861 -417892.551000 \n",
"2.5%HPD 0.033673 0.022859 0.025528 0.015623 -418201.078000 \n",
"97.5%HPD 0.035200 0.025239 0.027268 0.016855 -417891.188000 \n",
"ESS* 579.708878 147.436912 2641.537188 4378.739376 2087.556267 \n",
"Eff* 0.028985 0.007372 0.132077 0.218937 0.104378 "
]
},
"execution_count": 3,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# tables is a summary of the posteriors\n",
"table"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Converting units from the table results\n",
"The results are often difficult to interpret without converting the units into a more easily interpretable type. For example, the theta parameters represent the population effective mutation rate (`4*Ne*u`), but we are typically more interested in just knowing `N`. Similarly, the tau parameters are in units of percent sequence divergence (u/gen) but we instead like to know the divergence time estimates in units of years, or millions of years.\n",
"\n",
"To convert these units we need an estimate of the per-site per-generation mutation-rate (u) and the generation time of our organism (g). Rather than offering a point estimate it is more appropriate in the Bayesian context to provide these values a statistical distribution. This is similar to how we described our priors using a gamma or inverse-gamma distribution earlier, and how we used a distribution of u and g to check the reasonable conversion of the units.\n",
"\n",
"In the example below the converted units can be compared with the true simulation scenario at the beginning of this notebook where we set the divergence times and Ne values on the tree. Among the tip-level taxa we can see that mean Ne estimates are around 2e5 or ~6e5, which is pretty accurate, and the root divergence time is at 3e6, which is correct. "
]
},
{
"cell_type": "code",
"execution_count": 17,
"metadata": {
"scrolled": false
},
"outputs": [],
"source": [
"df = tool.transform(mcmc, 0.99, 1.01, 5e-9, 2e-8).T"
]
},
{
"cell_type": "code",
"execution_count": 18,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"
\n",
"\n",
"
\n",
" \n",
"
\n",
"
\n",
"
mean
\n",
"
median
\n",
"
S.D
\n",
"
min
\n",
"
max
\n",
"
\n",
" \n",
" \n",
"
\n",
"
Ne_1r0
\n",
"
6.16e+05
\n",
"
5.78e+05
\n",
"
2.05e+05
\n",
"
1.89e+05
\n",
"
2.55e+06
\n",
"
\n",
"
\n",
"
Ne_2r1
\n",
"
6.2e+05
\n",
"
5.83e+05
\n",
"
2.06e+05
\n",
"
1.9e+05
\n",
"
2.57e+06
\n",
"
\n",
"
\n",
"
Ne_3r2
\n",
"
5.63e+05
\n",
"
5.29e+05
\n",
"
1.87e+05
\n",
"
1.72e+05
\n",
"
2.28e+06
\n",
"
\n",
"
\n",
"
Ne_4r3
\n",
"
2.14e+05
\n",
"
2.01e+05
\n",
"
7.12e+04
\n",
"
6.57e+04
\n",
"
8.87e+05
\n",
"
\n",
"
\n",
"
Ne_5r4
\n",
"
2.04e+05
\n",
"
1.92e+05
\n",
"
6.8e+04
\n",
"
6.26e+04
\n",
"
8.47e+05
\n",
"
\n",
"
\n",
"
Ne_6r4r3r2r1r0
\n",
"
1.73e+05
\n",
"
1.62e+05
\n",
"
5.93e+04
\n",
"
5.47e+04
\n",
"
7.23e+05
\n",
"
\n",
"
\n",
"
Ne_7r4r3
\n",
"
1.2e+05
\n",
"
1.11e+05
\n",
"
4.54e+04
\n",
"
2.8e+04
\n",
"
5.7e+05
\n",
"
\n",
"
\n",
"
Ne_8r2r1r0
\n",
"
4.91e+05
\n",
"
4.59e+05
\n",
"
1.71e+05
\n",
"
1.45e+05
\n",
"
2.17e+06
\n",
"
\n",
"
\n",
"
Ne_9r1r0
\n",
"
3.92e+05
\n",
"
3.68e+05
\n",
"
1.32e+05
\n",
"
1.3e+05
\n",
"
1.71e+06
\n",
"
\n",
"
\n",
"
div_6r4r3r2r1r0
\n",
"
3.03e+06
\n",
"
2.85e+06
\n",
"
1e+06
\n",
"
9.51e+05
\n",
"
1.24e+07
\n",
"
\n",
"
\n",
"
div_7r4r3
\n",
"
2.11e+06
\n",
"
1.98e+06
\n",
"
7.01e+05
\n",
"
6.42e+05
\n",
"
8.56e+06
\n",
"
\n",
"
\n",
"
div_8r2r1r0
\n",
"
2.32e+06
\n",
"
2.18e+06
\n",
"
7.7e+05
\n",
"
7.2e+05
\n",
"
9.46e+06
\n",
"
\n",
"
\n",
"
div_9r1r0
\n",
"
1.43e+06
\n",
"
1.34e+06
\n",
"
4.74e+05
\n",
"
4.4e+05
\n",
"
5.81e+06
\n",
"
\n",
"
\n",
"
lnL
\n",
"
NaN
\n",
"
NaN
\n",
"
NaN
\n",
"
NaN
\n",
"
NaN
\n",
"
\n",
" \n",
"
\n",
"
"
],
"text/plain": [
" mean median S.D min max\n",
"Ne_1r0 6.16e+05 5.78e+05 2.05e+05 1.89e+05 2.55e+06\n",
"Ne_2r1 6.2e+05 5.83e+05 2.06e+05 1.9e+05 2.57e+06\n",
"Ne_3r2 5.63e+05 5.29e+05 1.87e+05 1.72e+05 2.28e+06\n",
"Ne_4r3 2.14e+05 2.01e+05 7.12e+04 6.57e+04 8.87e+05\n",
"Ne_5r4 2.04e+05 1.92e+05 6.8e+04 6.26e+04 8.47e+05\n",
"Ne_6r4r3r2r1r0 1.73e+05 1.62e+05 5.93e+04 5.47e+04 7.23e+05\n",
"Ne_7r4r3 1.2e+05 1.11e+05 4.54e+04 2.8e+04 5.7e+05\n",
"Ne_8r2r1r0 4.91e+05 4.59e+05 1.71e+05 1.45e+05 2.17e+06\n",
"Ne_9r1r0 3.92e+05 3.68e+05 1.32e+05 1.3e+05 1.71e+06\n",
"div_6r4r3r2r1r0 3.03e+06 2.85e+06 1e+06 9.51e+05 1.24e+07\n",
"div_7r4r3 2.11e+06 1.98e+06 7.01e+05 6.42e+05 8.56e+06\n",
"div_8r2r1r0 2.32e+06 2.18e+06 7.7e+05 7.2e+05 9.46e+06\n",
"div_9r1r0 1.43e+06 1.34e+06 4.74e+05 4.4e+05 5.81e+06\n",
"lnL NaN NaN NaN NaN NaN"
]
},
"execution_count": 18,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"df"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Plot posteriors\n",
"You can plot several posterior distributions on a shared axis using the `.draw_posteriors()` function. This takes a list of tuples as input, where each tuple is the (mean, variance) of a parameter estimate. You can draw the posteriors from a single analysis using this function or combine results from several different analyses if you wish. "
]
},
{
"cell_type": "code",
"execution_count": 43,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"
"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"# get only the divergence time result\n",
"subd = df.loc[[i for i in df.index if \"div_\" in i], :]\n",
"\n",
"# get results as tuple pairs with (mean, var)\n",
"tuples = [\n",
" (i, j) for (i, j) in \n",
" zip(subd[\"mean\"], subd.loc[:, \"S.D\"] ** 2)\n",
"]\n",
"\n",
"# draw the results\n",
"c, a, m = tool.draw_posteriors(\n",
" gamma_tuples=tuples,\n",
" labels=subd.index,\n",
");\n",
"\n",
"# style the axes\n",
"a.y.ticks.labels.angle = -90"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Plot results on a tree\n",
"\n",
"... coming"
]
}
],
"metadata": {
"anaconda-cloud": {},
"kernelspec": {
"display_name": "Python 3",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.6.7"
}
},
"nbformat": 4,
"nbformat_minor": 1
}