{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## (📗) Cookbook: `ipyrad.analysis.baba` \n",
"\n",
"This notebook demonstrates how to use the `baba` module in `ipyrad` to test for admixture and introgression. The code is written in Python and is best implemented in a Jupyter-notebook like this one. Analyses can be split up among many computing cores. Finally, there are some simple plotting methods to accompany tests which I show below. "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### import packages"
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
" ipyrad v.0.5.15\n",
" ipyparallel v.5.2.0\n",
" numpy v.1.11.2\n"
]
}
],
"source": [
"## imports\n",
"import numpy as np\n",
"import ipyrad as ip\n",
"import ipyrad.analysis.baba as baba\n",
"import ipyparallel as ipp\n",
"\n",
"## print versions\n",
"print \" ipyrad v.{}\".format(ip.__version__)\n",
"print \" ipyparallel v.{}\".format(ipp.__version__)\n",
"print \" numpy v.{}\".format(np.__version__)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### ipcluster setup to run parallel code\n",
"You can find more details about this in our [ipyparallel tutorial]."
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"## start 4 engines by running the commented line below in a separate bash terminal. \n",
"##\n",
"## ipcluster start --n=4"
]
},
{
"cell_type": "code",
"execution_count": 57,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
" host compute node: [4 cores] on oud\n"
]
}
],
"source": [
"## connect to client and print connection summary\n",
"ipyclient = ipp.Client()\n",
"print ip.cluster_info(client=ipyclient)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Get default (example) 12 taxon tree \n",
"The function `baba.Tree()` takes a `newick` string as an argument, or, if no value is passed, it constructs the default 12 taxon tree shown below. The `Tree` Class object holds a representation of the tree, which can be used for plotting, as well as for describing a demographic model for simulating data, as we'll show below. "
]
},
{
"cell_type": "code",
"execution_count": 58,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/html": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"## create a Tree class object with the default tree\n",
"tre = baba.Tree()\n",
"tre.show(width=400, height=300);"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Or get any arbitrary tree"
]
},
{
"cell_type": "code",
"execution_count": 59,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/html": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"## pass a newick string to the Tree class object\n",
"tre = baba.Tree(newick=\"((a,b),c);\")\n",
"tre.show(width=200, height=200);"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### get a tree with an admixture edge\n",
"Shown by an arrow pointing backwards from source to sink (indicating that geneflow should be entered as occurring backwards in time). This is how the program `msprime` will interpret it. "
]
},
{
"cell_type": "code",
"execution_count": 51,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/html": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"## store newick tree and a migration event list [source, sink, start, stop, rate]\n",
"newick = \"(((a,b),c), (d,e));\"\n",
"events = [['e', 'c', 0, 1, 1e-6], \n",
" ['3', '1', 1, 2, 1e-6]]\n",
"\n",
"## add admixture events to tree\n",
"tre = baba.Tree(newick=newick, admix=events)\n",
"tre.show(width=250, height=250, yaxis=True);"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### simulate data on that tree\n"
]
},
{
"cell_type": "code",
"execution_count": 33,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"8 [ 1. 0.]\n",
"7 [ 2. 0.]\n",
"6 [ 3. 0.]\n",
"5 [ 4. 0.]\n",
"4 [ 0.5 1. ]\n",
"3 [ 3.5 1. ]\n",
"2 [ 1.25 2. ]\n",
"1 [ 2.375 3. ]\n",
"0 [ 0. 0.]\n"
]
}
],
"source": [
"for idx in xrange(tre.verts.shape[0]-1, -1, -1):\n",
" print idx, tre.verts[idx-1]\n",
" row = tre.verts[idx-1]\n",
" \n",
" ## get time of the next event\n",
" if row[1] == 0:\n",
" pass\n",
" else:\n",
" tau = row[1]\n",
" \n"
]
},
{
"cell_type": "code",
"execution_count": 43,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"when [ 3. 2. 1. 1.]\n",
"who [0 1 2 3]\n"
]
}
],
"source": [
"## when do coalescent events occur?\n",
"who = tre.verts[:, 1] > 0\n",
"print 'when', tre.verts[who, 1]\n",
"print 'who', np.where(who)[0]\n",
"\n",
"## when do migration events occur\n",
"## ...\n",
"\n",
"for time in when:\n",
" ## does migration change during the next epoch\n"
]
},
{
"cell_type": "code",
"execution_count": 13,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"[[ 2.375 3. ]\n",
" [ 1.25 2. ]\n",
" [ 3.5 1. ]\n",
" [ 0.5 1. ]\n",
" [ 4. 0. ]\n",
" [ 3. 0. ]\n",
" [ 2. 0. ]\n",
" [ 1. 0. ]\n",
" [ 0. 0. ]]\n",
"[[1 8]\n",
" [0 1]\n",
" [2 6]\n",
" [3 4]\n",
" [3 5]\n",
" [2 3]\n",
" [0 2]\n",
" [1 7]]\n"
]
}
],
"source": [
"for node in tre.tree.traverse(\"postorder\"):\n",
" ## check for migration\n",
" if node.name in [i[0] for i in events]:\n",
" node\n",
" #ms.MigrationRateChange(time=0, rate=m_C_B, matrix_index=(1, 2)),\n",
" #ms.MigrationRateChange(time=Taus[1], rate=0),"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Simulate data"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": []
}
],
"metadata": {
"anaconda-cloud": {},
"kernelspec": {
"display_name": "Python [conda root]",
"language": "python",
"name": "conda-root-py"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 2
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython2",
"version": "2.7.12"
}
},
"nbformat": 4,
"nbformat_minor": 1
}