emp_nb_Quercus.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Notebook 3: \n",
    "\n",
    "This is an IPython notebook. Most of the code is composed of bash scripts, indicated by %%bash at the top of the cell, otherwise it is IPython code. This notebook includes code to download, assemble and analyze a published RADseq data set."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": []
    }
   ],
   "source": [
    "### Notebook 3\n",
    "### Data set 3 (American oaks)\n",
    "### Authors: Eaton et al. (2015)\n",
    "### Data Location: NCBI SRA SRP055977"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Download the sequence data\n",
    "Sequence data for this study are archived on the NCBI sequence read archive (SRA). Below I read in SraRunTable.txt for this project which contains all of the information we need to download the data. \n",
    "\n",
    "+ Project SRA: SRP055977\n",
    "+ BioProject ID: PRJNA277574\n",
    "+ Biosample numbers: SAMN03394519 - SAMN03394561\n",
    "+ Runs: SRR1915524 -- SRR1915566\n",
    "+ SRA link: http://trace.ncbi.nlm.nih.gov/Traces/study/?acc=SRP055977"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "%%bash\n",
    "## make a new directory for this analysis\n",
    "mkdir -p empirical_3/fastq/"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "    BioSample_s Library_Name_s  MBases_l  MBytes_l           Organism_s  \\\n",
      "0  SAMN03394519          AR_re       366       240    Quercus arizonica   \n",
      "1  SAMN03394520      BJSL25_re       484       313   Quercus brandegeei   \n",
      "2  SAMN03394521      BJVL19_re       427       275   Quercus brandegeei   \n",
      "3  SAMN03394522          CH_re       411       268  Quercus chrysolepis   \n",
      "4  SAMN03394523     CRL0001_re       339       220     Quercus oleoides   \n",
      "\n",
      "  ReleaseDate_s       Run_s SRA_Sample_s Sample_Name_s  \\\n",
      "0  Mar 13, 2015  SRR1915524    SRS868426         AR_re   \n",
      "1  Mar 13, 2015  SRR1915525    SRS874263     BJSL25_re   \n",
      "2  Mar 13, 2015  SRR1915526    SRS874262     BJVL19_re   \n",
      "3  Mar 13, 2015  SRR1915527    SRS874261         CH_re   \n",
      "4  Mar 13, 2015  SRR1915528    SRS874260    CRL0001_re   \n",
      "\n",
      "                                  geo_loc_name_s Assay_Type_s  AssemblyName_s  \\\n",
      "0              USA: UMN Greenhouse, St. Paul, MN          WGS  <not provided>   \n",
      "1         Mexico: Sierra la Laguna, Baja del Sur          WGS  <not provided>   \n",
      "2            Mexico: Valle Perdido, Baja del Sur          WGS  <not provided>   \n",
      "3              USA: UMN Greenhouse, St. Paul, MN          WGS  <not provided>   \n",
      "4  Costa Rica: RincÃ³n de la Vieja National Park          WGS  <not provided>   \n",
      "\n",
      "  BioProject_s BioSampleModel_s            Center_Name_s Consent_s  \\\n",
      "0  PRJNA277574            Plant  UNIVERSITY OF MINNESOTA    public   \n",
      "1  PRJNA277574            Plant  UNIVERSITY OF MINNESOTA    public   \n",
      "2  PRJNA277574            Plant  UNIVERSITY OF MINNESOTA    public   \n",
      "3  PRJNA277574            Plant  UNIVERSITY OF MINNESOTA    public   \n",
      "4  PRJNA277574            Plant  UNIVERSITY OF MINNESOTA    public   \n",
      "\n",
      "   InsertSize_l LibraryLayout_s    LoadDate_s Platform_s      \n",
      "0             0          SINGLE  Mar 13, 2015   ILLUMINA ...  \n",
      "1             0          SINGLE  Mar 13, 2015   ILLUMINA ...  \n",
      "2             0          SINGLE  Mar 13, 2015   ILLUMINA ...  \n",
      "3             0          SINGLE  Mar 13, 2015   ILLUMINA ...  \n",
      "4             0          SINGLE  Mar 13, 2015   ILLUMINA ...  \n",
      "\n",
      "[5 rows x 27 columns]\n"
     ]
    }
   ],
   "source": [
    "## IPython code\n",
    "import pandas as pd\n",
    "import urllib2\n",
    "import os\n",
    "\n",
    "## open the SRA run table from github url\n",
    "url = \"https://raw.githubusercontent.com/\"+\\\n",
    "      \"dereneaton/RADmissing/master/empirical_3_SraRunTable.txt\"\n",
    "intable = urllib2.urlopen(url)\n",
    "indata = pd.read_table(intable, sep=\"\\t\")\n",
    "\n",
    "## print first few rows\n",
    "print indata.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 47,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def wget_download(SRR, outdir, outname):\n",
    "    \"\"\" Python function to get sra data from ncbi and write to\n",
    "    outdir with a new name using bash call wget \"\"\"\n",
    "    \n",
    "    ## get output name\n",
    "    output = os.path.join(outdir, outname+\".sra\")\n",
    "    \n",
    "    ## create a call string \n",
    "    call = \"wget -q -r -nH --cut-dirs=9 -O \"+output+\" \"+\\\n",
    "           \"ftp://ftp-trace.ncbi.nlm.nih.gov/\"+\\\n",
    "           \"sra/sra-instant/reads/ByRun/sra/SRR/\"+\\\n",
    "           \"{}/{}/{}.sra;\".format(SRR[:6], SRR, SRR)\n",
    "        \n",
    "    ## call bash script\n",
    "    ! $call "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here we pass the SRR number and the sample name to the `wget_download` function so that the files are saved with their sample names. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 48,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "for ID, SRR in zip(indata.Library_Name_s, indata.Run_s):\n",
    "    wget_download(SRR, \"empirical_3/fastq/\", ID)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 49,
   "metadata": {
    "collapsed": false,
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Read 4046890 spots for empirical_3/fastq/AR_re.sra\n",
      "Written 4046890 spots for empirical_3/fastq/AR_re.sra\n",
      "Read 931926 spots for empirical_3/fastq/BJSB3_v.sra\n",
      "Written 931926 spots for empirical_3/fastq/BJSB3_v.sra\n",
      "Read 5352627 spots for empirical_3/fastq/BJSL25_re.sra\n",
      "Written 5352627 spots for empirical_3/fastq/BJSL25_re.sra\n",
      "Read 969575 spots for empirical_3/fastq/BJSL25_v.sra\n",
      "Written 969575 spots for empirical_3/fastq/BJSL25_v.sra\n",
      "Read 4715624 spots for empirical_3/fastq/BJVL19_re.sra\n",
      "Written 4715624 spots for empirical_3/fastq/BJVL19_re.sra\n",
      "Read 817443 spots for empirical_3/fastq/BJVL19_v.sra\n",
      "Written 817443 spots for empirical_3/fastq/BJVL19_v.sra\n",
      "Read 849191 spots for empirical_3/fastq/BZBB1_v.sra\n",
      "Written 849191 spots for empirical_3/fastq/BZBB1_v.sra\n",
      "Read 4539385 spots for empirical_3/fastq/CH_re.sra\n",
      "Written 4539385 spots for empirical_3/fastq/CH_re.sra\n",
      "Read 3742953 spots for empirical_3/fastq/CRL0001_re.sra\n",
      "Written 3742953 spots for empirical_3/fastq/CRL0001_re.sra\n",
      "Read 1012884 spots for empirical_3/fastq/CRL0001_v.sra\n",
      "Written 1012884 spots for empirical_3/fastq/CRL0001_v.sra\n",
      "Read 3242503 spots for empirical_3/fastq/CRL0030_re.sra\n",
      "Written 3242503 spots for empirical_3/fastq/CRL0030_re.sra\n",
      "Read 1204952 spots for empirical_3/fastq/CRL0030_v.sra\n",
      "Written 1204952 spots for empirical_3/fastq/CRL0030_v.sra\n",
      "Read 570148 spots for empirical_3/fastq/CUCA4_v.sra\n",
      "Written 570148 spots for empirical_3/fastq/CUCA4_v.sra\n",
      "Read 90661 spots for empirical_3/fastq/CUMM5_v.sra\n",
      "Written 90661 spots for empirical_3/fastq/CUMM5_v.sra\n",
      "Read 775357 spots for empirical_3/fastq/CUSV6_v.sra\n",
      "Written 775357 spots for empirical_3/fastq/CUSV6_v.sra\n",
      "Read 1130911 spots for empirical_3/fastq/CUVN10_v.sra\n",
      "Written 1130911 spots for empirical_3/fastq/CUVN10_v.sra\n",
      "Read 1953553 spots for empirical_3/fastq/DO_re.sra\n",
      "Written 1953553 spots for empirical_3/fastq/DO_re.sra\n",
      "Read 3744449 spots for empirical_3/fastq/DU_re.sra\n",
      "Written 3744449 spots for empirical_3/fastq/DU_re.sra\n",
      "Read 743556 spots for empirical_3/fastq/EN_re.sra\n",
      "Written 743556 spots for empirical_3/fastq/EN_re.sra\n",
      "Read 3219505 spots for empirical_3/fastq/FLAB109_re.sra\n",
      "Written 3219505 spots for empirical_3/fastq/FLAB109_re.sra\n",
      "Read 347949 spots for empirical_3/fastq/FLAB109_v.sra\n",
      "Written 347949 spots for empirical_3/fastq/FLAB109_v.sra\n",
      "Read 2826213 spots for empirical_3/fastq/FLBA140_re.sra\n",
      "Written 2826213 spots for empirical_3/fastq/FLBA140_re.sra\n",
      "Read 975393 spots for empirical_3/fastq/FLBA140_v.sra\n",
      "Written 975393 spots for empirical_3/fastq/FLBA140_v.sra\n",
      "Read 941880 spots for empirical_3/fastq/FLCK18_v.sra\n",
      "Written 941880 spots for empirical_3/fastq/FLCK18_v.sra\n",
      "Read 791393 spots for empirical_3/fastq/FLCK216_v.sra\n",
      "Written 791393 spots for empirical_3/fastq/FLCK216_v.sra\n",
      "Read 905151 spots for empirical_3/fastq/FLMO62_v.sra\n",
      "Written 905151 spots for empirical_3/fastq/FLMO62_v.sra\n",
      "Read 3976506 spots for empirical_3/fastq/FLSA185_v.sra\n",
      "Written 3976506 spots for empirical_3/fastq/FLSA185_v.sra\n",
      "Read 857771 spots for empirical_3/fastq/FLSF33_v.sra\n",
      "Written 857771 spots for empirical_3/fastq/FLSF33_v.sra\n",
      "Read 908978 spots for empirical_3/fastq/FLSF47_v.sra\n",
      "Written 908978 spots for empirical_3/fastq/FLSF47_v.sra\n",
      "Read 2543324 spots for empirical_3/fastq/FLSF54_re.sra\n",
      "Written 2543324 spots for empirical_3/fastq/FLSF54_re.sra\n",
      "Read 1115926 spots for empirical_3/fastq/FLSF54_v.sra\n",
      "Written 1115926 spots for empirical_3/fastq/FLSF54_v.sra\n",
      "Read 1099920 spots for empirical_3/fastq/FLWO6_v.sra\n",
      "Written 1099920 spots for empirical_3/fastq/FLWO6_v.sra\n",
      "Read 2777404 spots for empirical_3/fastq/HE_re.sra\n",
      "Written 2777404 spots for empirical_3/fastq/HE_re.sra\n",
      "Read 1112292 spots for empirical_3/fastq/HNDA09_v.sra\n",
      "Written 1112292 spots for empirical_3/fastq/HNDA09_v.sra\n",
      "Read 1436228 spots for empirical_3/fastq/LALC2_v.sra\n",
      "Written 1436228 spots for empirical_3/fastq/LALC2_v.sra\n",
      "Read 1154841 spots for empirical_3/fastq/MXED8_v.sra\n",
      "Written 1154841 spots for empirical_3/fastq/MXED8_v.sra\n",
      "Read 1347187 spots for empirical_3/fastq/MXGT4_v.sra\n",
      "Written 1347187 spots for empirical_3/fastq/MXGT4_v.sra\n",
      "Read 1046574 spots for empirical_3/fastq/MXSA3017_v.sra\n",
      "Written 1046574 spots for empirical_3/fastq/MXSA3017_v.sra\n",
      "Read 4482353 spots for empirical_3/fastq/NI_re.sra\n",
      "Written 4482353 spots for empirical_3/fastq/NI_re.sra\n",
      "Read 505504 spots for empirical_3/fastq/SCCU3_v.sra\n",
      "Written 505504 spots for empirical_3/fastq/SCCU3_v.sra\n",
      "Read 879229 spots for empirical_3/fastq/TXGR3_v.sra\n",
      "Written 879229 spots for empirical_3/fastq/TXGR3_v.sra\n",
      "Read 934386 spots for empirical_3/fastq/TXMD3_v.sra\n",
      "Written 934386 spots for empirical_3/fastq/TXMD3_v.sra\n",
      "Read 163427 spots for empirical_3/fastq/TXWV2_v.sra\n",
      "Written 163427 spots for empirical_3/fastq/TXWV2_v.sra\n",
      "Read 76783922 spots total\n",
      "Written 76783922 spots total\n"
     ]
    }
   ],
   "source": [
    "%%bash\n",
    "## convert sra files to fastq using fastq-dump tool\n",
    "## output as gzipped into the fastq directory\n",
    "fastq-dump --gzip -O empirical_3/fastq/  empirical_3/fastq/*.sra\n",
    "\n",
    "## remove .sra files\n",
    "rm empirical_3/fastq/*.sra"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This study includes several re-sequenced individuals. We combine them before beginning the analysis."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 50,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "BJVL19 merged\n",
      "CRL0001 merged\n",
      "BJSL25 merged\n",
      "FLBA140 merged\n",
      "FLAB109 merged\n",
      "CRL0030 merged\n",
      "FLSF54 merged\n"
     ]
    }
   ],
   "source": [
    "##IPython code\n",
    "import glob\n",
    "\n",
    "taxa = [i.split(\"/\")[-1].split('_')[0] for i in glob.glob(\"empirical_3/fastq/*.gz\")]\n",
    "for taxon in set(taxa):\n",
    "    if taxa.count(taxon) > 1:\n",
    "        print taxon, \"merged\"\n",
    "        \n",
    "        ## merge replicate files\n",
    "        ! cat empirical_3/fastq/$taxon\\_v.fastq.gz \\\n",
    "              empirical_3/fastq/$taxon\\_re.fastq.gz \\\n",
    "            > empirical_3/fastq/$taxon\\_me.fastq.gz\n",
    "            \n",
    "        ## remove ind replicate files\n",
    "        ! rm empirical_3/fastq/$taxon\\_v.fastq.gz\n",
    "        ! rm empirical_3/fastq/$taxon\\_re.fastq.gz"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Make a params file"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 51,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "pyRAD 3.0.63\n"
     ]
    }
   ],
   "source": [
    "%%bash\n",
    "pyrad --version"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 52,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "\tnew params.txt file created\n"
     ]
    }
   ],
   "source": [
    "%%bash\n",
    "## remove old params file if it exists\n",
    "rm params.txt \n",
    "\n",
    "## create a new default params file\n",
    "pyrad -n "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "%%bash\n",
    "## substitute new parameters into file\n",
    "sed -i '/## 1. /c\\empirical_3/           ## 1. working directory ' params.txt\n",
    "sed -i '/## 6. /c\\TGCAG                  ## 6. cutters ' params.txt\n",
    "sed -i '/## 7. /c\\20                     ## 7. N processors      ' params.txt\n",
    "sed -i '/## 9. /c\\6                      ## 9. NQual             ' params.txt\n",
    "sed -i '/## 10./c\\.85                    ## 10. clust threshold  ' params.txt\n",
    "sed -i '/## 12./c\\4                      ## 12. MinCov           ' params.txt\n",
    "sed -i '/## 13./c\\10                     ## 13. maxSH            ' params.txt\n",
    "sed -i '/## 14./c\\empirical_3_m4         ## 14. output name      ' params.txt\n",
    "sed -i '/## 18./c\\empirical_3/fastq/*.gz ## 18. data location    ' params.txt\n",
    "sed -i '/## 29./c\\2,2                    ## 29. trim overhang    ' params.txt\n",
    "sed -i '/## 30./c\\p,n,s                  ## 30. output formats   ' params.txt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "==** parameter inputs for pyRAD version 3.0.63  **======================== affected step ==\r\n",
      "empirical_3/           ## 1. working directory \r\n",
      "./*.fastq.gz              ## 2. Loc. of non-demultiplexed files (if not line 18)  (s1)\r\n",
      "./*.barcodes              ## 3. Loc. of barcode file (if not line 18)             (s1)\r\n",
      "vsearch                   ## 4. command (or path) to call vsearch (or usearch)    (s3,s6)\r\n",
      "muscle                    ## 5. command (or path) to call muscle                  (s3,s7)\r\n",
      "TGCAG                  ## 6. cutters \r\n",
      "20                     ## 7. N processors      \r\n",
      "6                         ## 8. Mindepth: min coverage for a cluster              (s4,s5)\r\n",
      "6                      ## 9. NQual             \r\n",
      ".85                    ## 10. clust threshold  \r\n",
      "rad                       ## 11. Datatype: rad,gbs,pairgbs,pairddrad,(others:see docs)(all)\r\n",
      "4                      ## 12. MinCov           \r\n",
      "10                     ## 13. maxSH            \r\n",
      "empirical_3_m4         ## 14. output name      \r\n",
      "==== optional params below this line ===================================  affected step ==\r\n",
      "                       ## 15.opt.: select subset (prefix* only selector)            (s2-s7)\r\n",
      "                       ## 16.opt.: add-on (outgroup) taxa (list or prefix*)         (s6,s7)\r\n",
      "                       ## 17.opt.: exclude taxa (list or prefix*)                   (s7)\r\n",
      "empirical_3/fastq/*.gz ## 18. data location    \r\n",
      "                       ## 19.opt.: maxM: N mismatches in barcodes (def= 1)          (s1)\r\n",
      "                       ## 20.opt.: phred Qscore offset (def= 33)                    (s2)\r\n",
      "                       ## 21.opt.: filter: def=0=NQual 1=NQual+adapters. 2=strict   (s2)\r\n",
      "                       ## 22.opt.: a priori E,H (def= 0.001,0.01, if not estimated) (s5)\r\n",
      "                       ## 23.opt.: maxN: max Ns in a cons seq (def=5)               (s5)\r\n",
      "                       ## 24.opt.: maxH: max heterozyg. sites in cons seq (def=5)   (s5)\r\n",
      "                       ## 25.opt.: ploidy: max alleles in cons seq (def=2;see docs) (s4,s5)\r\n",
      "                       ## 26.opt.: maxSNPs: (def=100). Paired (def=100,100)         (s7)\r\n",
      "                       ## 27.opt.: maxIndels: within-clust,across-clust (def. 3,99) (s3,s7)\r\n",
      "                       ## 28.opt.: random number seed (def. 112233)              (s3,s6,s7)\r\n",
      "2,2                    ## 29. trim overhang    \r\n",
      "p,n,s                  ## 30. output formats   \r\n",
      "                       ## 31.opt.: maj. base call at depth>x<mindepth (def.x=mindepth) (s5)\r\n",
      "                       ## 32.opt.: keep trimmed reads (def=0). Enter min length.    (s2)\r\n",
      "                       ## 33.opt.: max stack size (int), def= max(500,mean+2*SD)    (s3)\r\n",
      "                       ## 34.opt.: minDerep: exclude dereps with <= N copies, def=1 (s3)\r\n",
      "                       ## 35.opt.: use hierarchical clustering (def.=0, 1=yes)      (s6)\r\n",
      "                       ## 36.opt.: repeat masking (def.=1='dust' method, 0=no)      (s3,s6)\r\n",
      "                       ## 37.opt.: vsearch max threads per job (def.=6; see docs)   (s3,s6)\r\n",
      "==== optional: list group/clade assignments below this line (see docs) ==================\r\n"
     ]
    }
   ],
   "source": [
    "cat params.txt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "## Assemble in _pyrad_\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "%%bash\n",
    "pyrad -p params.txt -s 234567 >> log.txt 2>&1  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "%%bash\n",
    "sed -i '/## 12./c\\2                      ## 12. MinCov           ' params.txt\n",
    "sed -i '/## 14./c\\empirical_3_m2         ## 14. output name      ' params.txt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "%%bash\n",
    "pyrad -p params.txt -s 7 >> log.txt 2>&1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Results\n",
    "We are interested in the relationship between the amount of input (raw) data between any two samples, the average coverage they recover when clustered together, and the phylogenetic distances separating samples.  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Raw data amounts\n",
    "The average number of raw reads per sample is 1.36M."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "count         36.000000\n",
      "mean     1955564.750000\n",
      "std      1594388.931033\n",
      "min        84505.000000\n",
      "25%       806803.750000\n",
      "50%      1050559.000000\n",
      "75%      3425603.250000\n",
      "max      5829065.000000\n",
      "Name: passed.total, dtype: float64\n",
      "\n",
      "most raw data in sample:\n",
      "2    BJSL25_me\n",
      "Name: sample , dtype: object\n"
     ]
    }
   ],
   "source": [
    "## read in the data\n",
    "s2dat = pd.read_table(\"empirical_3/stats/s2.rawedit.txt\", header=0, nrows=36)\n",
    "\n",
    "## print summary stats\n",
    "print s2dat[\"passed.total\"].describe()\n",
    "\n",
    "## find which sample has the most raw data\n",
    "maxraw = s2dat[\"passed.total\"].max()\n",
    "print \"\\nmost raw data in sample:\"\n",
    "print s2dat['sample '][s2dat['passed.total']==maxraw]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Look at distributions of coverage\n",
    "pyrad v.3.0.63 outputs depth information for each sample which I read in here and plot. First let's ask which sample has the highest depth of coverage. The std of coverages is pretty low in this data set compared to several others. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": false,
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "summary of means\n",
      "==================\n",
      "count    36.00000\n",
      "mean     17.80300\n",
      "std      13.02276\n",
      "min       3.31500\n",
      "25%       9.02300\n",
      "50%      11.61800\n",
      "75%      25.41200\n",
      "max      51.11900\n",
      "Name: dpt.me, dtype: float64\n",
      "\n",
      "summary of std\n",
      "==================\n",
      "count     36.000000\n",
      "mean      70.155611\n",
      "std       59.520107\n",
      "min        3.718000\n",
      "25%       28.191000\n",
      "50%       47.720500\n",
      "75%       98.996000\n",
      "max      211.783000\n",
      "Name: dpt.sd, dtype: float64\n",
      "\n",
      "summary of proportion lowdepth\n",
      "==================\n",
      "count    36.000000\n",
      "mean      0.427849\n",
      "std       0.157876\n",
      "min       0.211693\n",
      "25%       0.300563\n",
      "50%       0.440962\n",
      "75%       0.497383\n",
      "max       0.860847\n",
      "dtype: float64\n",
      "\n",
      "highest coverage in sample:\n",
      "0    AR_re\n",
      "Name: taxa, dtype: object\n"
     ]
    }
   ],
   "source": [
    "## read in the s3 results\n",
    "s3dat = pd.read_table(\"empirical_3/stats/s3.clusters.txt\", header=0, nrows=39)\n",
    "\n",
    "## print summary stats\n",
    "print \"summary of means\\n==================\"\n",
    "print s3dat['dpt.me'].describe()\n",
    "\n",
    "## print summary stats\n",
    "print \"\\nsummary of std\\n==================\"\n",
    "print s3dat['dpt.sd'].describe()\n",
    "\n",
    "## print summary stats\n",
    "print \"\\nsummary of proportion lowdepth\\n==================\"\n",
    "print pd.Series(1-s3dat['d>5.tot']/s3dat[\"total\"]).describe()\n",
    "\n",
    "## find which sample has the greatest depth of retained loci\n",
    "max_hiprop = (s3dat[\"d>5.tot\"]/s3dat[\"total\"]).max()\n",
    "print \"\\nhighest coverage in sample:\"\n",
    "print s3dat['taxa'][s3dat['d>5.tot']/s3dat[\"total\"]==max_hiprop]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "37.41198396 122.321191848\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "## print mean and std of coverage for the highest coverage sample\n",
    "with open(\"empirical_3/clust.85/AR_re.depths\", 'rb') as indat:\n",
    "    depths = np.array(indat.read().strip().split(\",\"), dtype=int)\n",
    "    \n",
    "print depths.mean(), depths.std()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Plot the coverage for the sample with highest mean coverage\n",
    "Green shows the loci that were discarded and orange the loci that were retained. The majority of data were discarded for being too low of coverage. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<toyplot.mark.BarMagnitudes at 0x7f5de591dad0>"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "text/html": [
       "<div align=\"center\" class=\"toyplot\" id=\"t1082bb5d7c604b1fa934d98fbaf2df7d\"><svg height=\"300.0px\" id=\"t43fff6ddc57040f3af18a9dd87d95d9a\" preserveAspectRatio=\"xMidyMid meet\" style=\"background-color:transparent;fill:rgb(16.1%,15.3%,14.1%);fill-opacity:1.0;font-family:helvetica;font-size:12px;opacity:1.0;stroke:rgb(16.1%,15.3%,14.1%);stroke-opacity:1.0;stroke-width:1.0\" viewBox=\"0 0 350.0 300.0\" width=\"350.0px\" xmlns=\"http://www.w3.org/2000/svg\" xmlns:toyplot=\"http://www.sandia.gov/toyplot\"><g class=\"toyplot-axes-Cartesian\" id=\"t551d9eff8bca4b2b81d94d8a92ad895e\"><clipPath id=\"tf9316bf2480a448d8bc5634bdb97ad7d\"><rect height=\"200.0\" width=\"250.0\" x=\"50.0\" y=\"50.0\"></rect></clipPath><g class=\"toyplot-coordinate-events\" clip-path=\"url(#tf9316bf2480a448d8bc5634bdb97ad7d)\" style=\"cursor:crosshair\"><rect height=\"200.0\" style=\"pointer-events:all;visibility:hidden\" width=\"250.0\" x=\"50.0\" y=\"50.0\"></rect><g class=\"toyplot-mark-BarMagnitudes\" id=\"te2ee6288007b4596ab930fed0d1155a2\" style=\"stroke:rgb(100%,100%,100%);stroke-opacity:1.0;stroke-width:1.0\"><g class=\"toyplot-Series\"><rect class=\"toyplot-Datum\" height=\"0.0\" style=\"fill:rgb(40%,76.1%,64.7%);fill-opacity:1.0;opacity:1.0;stroke:rgb(100%,100%,100%);stroke-opacity:1.0;stroke-width:1.0\" width=\"7.6666666666666714\" x=\"60.0\" y=\"240.0\"></rect><rect class=\"toyplot-Datum\" height=\"180.0\" style=\"fill:rgb(40%,76.1%,64.7%);fill-opacity:1.0;opacity:1.0;stroke:rgb(100%,100%,100%);stroke-opacity:1.0;stroke-width:1.0\" width=\"7.6666666666666572\" x=\"67.666666666666671\" y=\"60.0\"></rect><rect class=\"toyplot-Datum\" height=\"75.984776403425286\" style=\"fill:rgb(40%,76.1%,64.7%);fill-opacity:1.0;opacity:1.0;stroke:rgb(100%,100%,100%);stroke-opacity:1.0;stroke-width:1.0\" width=\"7.6666666666666714\" x=\"75.333333333333329\" y=\"164.01522359657471\"></rect><rect class=\"toyplot-Datum\" height=\"49.571836346336852\" style=\"fill:rgb(40%,76.1%,64.7%);fill-opacity:1.0;opacity:1.0;stroke:rgb(100%,100%,100%);stroke-opacity:1.0;stroke-width:1.0\" width=\"7.6666666666666572\" x=\"83.0\" y=\"190.42816365366315\"></rect><rect class=\"toyplot-Datum\" height=\"40.171265461465254\" style=\"fill:rgb(40%,76.1%,64.7%);fill-opacity:1.0;opacity:1.0;stroke:rgb(100%,100%,100%);stroke-opacity:1.0;stroke-width:1.0\" width=\"7.6666666666666714\" x=\"90.666666666666657\" y=\"199.82873453853475\"></rect><rect class=\"toyplot-Datum\" height=\"39.029495718363478\" style=\"fill:rgb(40%,76.1%,64.7%);fill-opacity:1.0;opacity:1.0;stroke:rgb(100%,100%,100%);stroke-opacity:1.0;stroke-width:1.0\" width=\"7.6666666666666714\" x=\"98.333333333333329\" y=\"200.97050428163652\"></rect><rect class=\"toyplot-Datum\" height=\"36.213130352045681\" style=\"fill:rgb(40%,76.1%,64.7%);fill-opacity:1.0;opacity:1.0;stroke:rgb(100%,100%,100%);stroke-opacity:1.0;stroke-width:1.0\" width=\"7.6666666666666714\" x=\"106.0\" y=\"203.78686964795432\"></rect><rect class=\"toyplot-Datum\" height=\"34.576593720266402\" style=\"fill:rgb(40%,76.1%,64.7%);fill-opacity:1.0;opacity:1.0;stroke:rgb(100%,100%,100%);stroke-opacity:1.0;stroke-width:1.0\" width=\"7.6666666666666572\" x=\"113.66666666666667\" y=\"205.4234062797336\"></rect><rect class=\"toyplot-Datum\" height=\"33.510941960038053\" style=\"fill:rgb(40%,76.1%,64.7%);fill-opacity:1.0;opacity:1.0;stroke:rgb(100%,100%,100%);stroke-opacity:1.0;stroke-width:1.0\" width=\"7.6666666666666714\" x=\"121.33333333333333\" y=\"206.48905803996195\"></rect><rect class=\"toyplot-Datum\" height=\"31.969552806850629\" style=\"fill:rgb(40%,76.1%,64.7%);fill-opacity:1.0;opacity:1.0;stroke:rgb(100%,100%,100%);stroke-opacity:1.0;stroke-width:1.0\" width=\"7.6666666666666572\" x=\"129.0\" y=\"208.03044719314937\"></rect><rect class=\"toyplot-Datum\" height=\"31.569933396764981\" style=\"fill:rgb(40%,76.1%,64.7%);fill-opacity:1.0;opacity:1.0;stroke:rgb(100%,100%,100%);stroke-opacity:1.0;stroke-width:1.0\" width=\"7.6666666666666572\" x=\"136.66666666666666\" y=\"208.43006660323502\"></rect><rect class=\"toyplot-Datum\" height=\"28.125594671741197\" style=\"fill:rgb(40%,76.1%,64.7%);fill-opacity:1.0;opacity:1.0;stroke:rgb(100%,100%,100%);stroke-opacity:1.0;stroke-width:1.0\" width=\"7.6666666666666856\" x=\"144.33333333333331\" y=\"211.8744053282588\"></rect><rect class=\"toyplot-Datum\" height=\"28.125594671741197\" style=\"fill:rgb(40%,76.1%,64.7%);fill-opacity:1.0;opacity:1.0;stroke:rgb(100%,100%,100%);stroke-opacity:1.0;stroke-width:1.0\" width=\"7.6666666666666856\" x=\"152.0\" y=\"211.8744053282588\"></rect><rect class=\"toyplot-Datum\" height=\"27.706945765937178\" style=\"fill:rgb(40%,76.1%,64.7%);fill-opacity:1.0;opacity:1.0;stroke:rgb(100%,100%,100%);stroke-opacity:1.0;stroke-width:1.0\" width=\"7.6666666666666572\" x=\"159.66666666666669\" y=\"212.29305423406282\"></rect><rect class=\"toyplot-Datum\" height=\"26.565176022835402\" style=\"fill:rgb(40%,76.1%,64.7%);fill-opacity:1.0;opacity:1.0;stroke:rgb(100%,100%,100%);stroke-opacity:1.0;stroke-width:1.0\" width=\"7.6666666666666572\" x=\"167.33333333333334\" y=\"213.4348239771646\"></rect><rect class=\"toyplot-Datum\" height=\"25.233111322549973\" style=\"fill:rgb(40%,76.1%,64.7%);fill-opacity:1.0;opacity:1.0;stroke:rgb(100%,100%,100%);stroke-opacity:1.0;stroke-width:1.0\" width=\"7.6666666666666572\" x=\"175.0\" y=\"214.76688867745003\"></rect><rect class=\"toyplot-Datum\" height=\"24.376784015223592\" style=\"fill:rgb(40%,76.1%,64.7%);fill-opacity:1.0;opacity:1.0;stroke:rgb(100%,100%,100%);stroke-opacity:1.0;stroke-width:1.0\" width=\"7.6666666666666856\" x=\"182.66666666666666\" y=\"215.62321598477641\"></rect><rect class=\"toyplot-Datum\" height=\"25.118934348239776\" style=\"fill:rgb(40%,76.1%,64.7%);fill-opacity:1.0;opacity:1.0;stroke:rgb(100%,100%,100%);stroke-opacity:1.0;stroke-width:1.0\" width=\"7.6666666666666572\" x=\"190.33333333333334\" y=\"214.88106565176022\"></rect><rect class=\"toyplot-Datum\" height=\"24.319695528068507\" style=\"fill:rgb(40%,76.1%,64.7%);fill-opacity:1.0;opacity:1.0;stroke:rgb(100%,100%,100%);stroke-opacity:1.0;stroke-width:1.0\" width=\"7.6666666666666572\" x=\"198.0\" y=\"215.68030447193149\"></rect><rect class=\"toyplot-Datum\" height=\"24.262607040913394\" style=\"fill:rgb(40%,76.1%,64.7%);fill-opacity:1.0;opacity:1.0;stroke:rgb(100%,100%,100%);stroke-opacity:1.0;stroke-width:1.0\" width=\"7.6666666666666572\" x=\"205.66666666666666\" y=\"215.73739295908661\"></rect><rect class=\"toyplot-Datum\" height=\"24.205518553758338\" style=\"fill:rgb(40%,76.1%,64.7%);fill-opacity:1.0;opacity:1.0;stroke:rgb(100%,100%,100%);stroke-opacity:1.0;stroke-width:1.0\" width=\"7.6666666666666856\" x=\"213.33333333333331\" y=\"215.79448144624166\"></rect><rect class=\"toyplot-Datum\" height=\"22.98763082778305\" style=\"fill:rgb(40%,76.1%,64.7%);fill-opacity:1.0;opacity:1.0;stroke:rgb(100%,100%,100%);stroke-opacity:1.0;stroke-width:1.0\" width=\"7.6666666666666572\" x=\"221.0\" y=\"217.01236917221695\"></rect><rect class=\"toyplot-Datum\" height=\"23.748810656517605\" style=\"fill:rgb(40%,76.1%,64.7%);fill-opacity:1.0;opacity:1.0;stroke:rgb(100%,100%,100%);stroke-opacity:1.0;stroke-width:1.0\" width=\"7.6666666666666856\" x=\"228.66666666666666\" y=\"216.25118934348239\"></rect><rect class=\"toyplot-Datum\" height=\"22.85442435775451\" style=\"fill:rgb(40%,76.1%,64.7%);fill-opacity:1.0;opacity:1.0;stroke:rgb(100%,100%,100%);stroke-opacity:1.0;stroke-width:1.0\" width=\"7.6666666666666572\" x=\"236.33333333333334\" y=\"217.14557564224549\"></rect><rect class=\"toyplot-Datum\" height=\"23.463368220742154\" style=\"fill:rgb(40%,76.1%,64.7%);fill-opacity:1.0;opacity:1.0;stroke:rgb(100%,100%,100%);stroke-opacity:1.0;stroke-width:1.0\" width=\"7.6666666666666856\" x=\"244.0\" y=\"216.53663177925785\"></rect><rect class=\"toyplot-Datum\" height=\"22.568981921979088\" style=\"fill:rgb(40%,76.1%,64.7%);fill-opacity:1.0;opacity:1.0;stroke:rgb(100%,100%,100%);stroke-opacity:1.0;stroke-width:1.0\" width=\"7.6666666666666288\" x=\"251.66666666666669\" y=\"217.43101807802091\"></rect><rect class=\"toyplot-Datum\" height=\"22.226450999048524\" style=\"fill:rgb(40%,76.1%,64.7%);fill-opacity:1.0;opacity:1.0;stroke:rgb(100%,100%,100%);stroke-opacity:1.0;stroke-width:1.0\" width=\"7.6666666666666856\" x=\"259.33333333333331\" y=\"217.77354900095148\"></rect><rect class=\"toyplot-Datum\" height=\"21.274976213130373\" style=\"fill:rgb(40%,76.1%,64.7%);fill-opacity:1.0;opacity:1.0;stroke:rgb(100%,100%,100%);stroke-opacity:1.0;stroke-width:1.0\" width=\"7.6666666666666856\" x=\"267.0\" y=\"218.72502378686963\"></rect><rect class=\"toyplot-Datum\" height=\"43.063748810656534\" style=\"fill:rgb(40%,76.1%,64.7%);fill-opacity:1.0;opacity:1.0;stroke:rgb(100%,100%,100%);stroke-opacity:1.0;stroke-width:1.0\" width=\"7.6666666666666288\" x=\"274.66666666666669\" y=\"196.93625118934347\"></rect></g></g><g class=\"toyplot-mark-BarMagnitudes\" id=\"t5fd1c6e1bcf34c77aace5defeb3bff5a\" style=\"stroke:rgb(100%,100%,100%);stroke-opacity:1.0;stroke-width:1.0\"><g class=\"toyplot-Series\"><rect class=\"toyplot-Datum\" height=\"0.0\" style=\"fill:rgb(98.8%,55.3%,38.4%);fill-opacity:1.0;opacity:1.0;stroke:rgb(100%,100%,100%);stroke-opacity:1.0;stroke-width:1.0\" width=\"7.6666666666666714\" x=\"60.0\" y=\"240.0\"></rect><rect class=\"toyplot-Datum\" height=\"0.0\" style=\"fill:rgb(98.8%,55.3%,38.4%);fill-opacity:1.0;opacity:1.0;stroke:rgb(100%,100%,100%);stroke-opacity:1.0;stroke-width:1.0\" width=\"7.6666666666666572\" x=\"67.666666666666671\" y=\"240.0\"></rect><rect class=\"toyplot-Datum\" height=\"0.0\" style=\"fill:rgb(98.8%,55.3%,38.4%);fill-opacity:1.0;opacity:1.0;stroke:rgb(100%,100%,100%);stroke-opacity:1.0;stroke-width:1.0\" width=\"7.6666666666666714\" x=\"75.333333333333329\" y=\"240.0\"></rect><rect class=\"toyplot-Datum\" height=\"0.0\" style=\"fill:rgb(98.8%,55.3%,38.4%);fill-opacity:1.0;opacity:1.0;stroke:rgb(100%,100%,100%);stroke-opacity:1.0;stroke-width:1.0\" width=\"7.6666666666666572\" x=\"83.0\" y=\"240.0\"></rect><rect class=\"toyplot-Datum\" height=\"0.0\" style=\"fill:rgb(98.8%,55.3%,38.4%);fill-opacity:1.0;opacity:1.0;stroke:rgb(100%,100%,100%);stroke-opacity:1.0;stroke-width:1.0\" width=\"7.6666666666666714\" x=\"90.666666666666657\" y=\"240.0\"></rect><rect class=\"toyplot-Datum\" height=\"0.0\" style=\"fill:rgb(98.8%,55.3%,38.4%);fill-opacity:1.0;opacity:1.0;stroke:rgb(100%,100%,100%);stroke-opacity:1.0;stroke-width:1.0\" width=\"7.6666666666666714\" x=\"98.333333333333329\" y=\"240.0\"></rect><rect class=\"toyplot-Datum\" height=\"36.213130352045681\" style=\"fill:rgb(98.8%,55.3%,38.4%);fill-opacity:1.0;opacity:1.0;stroke:rgb(100%,100%,100%);stroke-opacity:1.0;stroke-width:1.0\" width=\"7.6666666666666714\" x=\"106.0\" y=\"203.78686964795432\"></rect><rect class=\"toyplot-Datum\" height=\"34.576593720266402\" style=\"fill:rgb(98.8%,55.3%,38.4%);fill-opacity:1.0;opacity:1.0;stroke:rgb(100%,100%,100%);stroke-opacity:1.0;stroke-width:1.0\" width=\"7.6666666666666572\" x=\"113.66666666666667\" y=\"205.4234062797336\"></rect><rect class=\"toyplot-Datum\" height=\"33.510941960038053\" style=\"fill:rgb(98.8%,55.3%,38.4%);fill-opacity:1.0;opacity:1.0;stroke:rgb(100%,100%,100%);stroke-opacity:1.0;stroke-width:1.0\" width=\"7.6666666666666714\" x=\"121.33333333333333\" y=\"206.48905803996195\"></rect><rect class=\"toyplot-Datum\" height=\"31.969552806850629\" style=\"fill:rgb(98.8%,55.3%,38.4%);fill-opacity:1.0;opacity:1.0;stroke:rgb(100%,100%,100%);stroke-opacity:1.0;stroke-width:1.0\" width=\"7.6666666666666572\" x=\"129.0\" y=\"208.03044719314937\"></rect><rect class=\"toyplot-Datum\" height=\"31.569933396764981\" style=\"fill:rgb(98.8%,55.3%,38.4%);fill-opacity:1.0;opacity:1.0;stroke:rgb(100%,100%,100%);stroke-opacity:1.0;stroke-width:1.0\" width=\"7.6666666666666572\" x=\"136.66666666666666\" y=\"208.43006660323502\"></rect><rect class=\"toyplot-Datum\" height=\"28.125594671741197\" style=\"fill:rgb(98.8%,55.3%,38.4%);fill-opacity:1.0;opacity:1.0;stroke:rgb(100%,100%,100%);stroke-opacity:1.0;stroke-width:1.0\" width=\"7.6666666666666856\" x=\"144.33333333333331\" y=\"211.8744053282588\"></rect><rect class=\"toyplot-Datum\" height=\"28.125594671741197\" style=\"fill:rgb(98.8%,55.3%,38.4%);fill-opacity:1.0;opacity:1.0;stroke:rgb(100%,100%,100%);stroke-opacity:1.0;stroke-width:1.0\" width=\"7.6666666666666856\" x=\"152.0\" y=\"211.8744053282588\"></rect><rect class=\"toyplot-Datum\" height=\"27.706945765937178\" style=\"fill:rgb(98.8%,55.3%,38.4%);fill-opacity:1.0;opacity:1.0;stroke:rgb(100%,100%,100%);stroke-opacity:1.0;stroke-width:1.0\" width=\"7.6666666666666572\" x=\"159.66666666666669\" y=\"212.29305423406282\"></rect><rect class=\"toyplot-Datum\" height=\"26.565176022835402\" style=\"fill:rgb(98.8%,55.3%,38.4%);fill-opacity:1.0;opacity:1.0;stroke:rgb(100%,100%,100%);stroke-opacity:1.0;stroke-width:1.0\" width=\"7.6666666666666572\" x=\"167.33333333333334\" y=\"213.4348239771646\"></rect><rect class=\"toyplot-Datum\" height=\"25.233111322549973\" style=\"fill:rgb(98.8%,55.3%,38.4%);fill-opacity:1.0;opacity:1.0;stroke:rgb(100%,100%,100%);stroke-opacity:1.0;stroke-width:1.0\" width=\"7.6666666666666572\" x=\"175.0\" y=\"214.76688867745003\"></rect><rect class=\"toyplot-Datum\" height=\"24.376784015223592\" style=\"fill:rgb(98.8%,55.3%,38.4%);fill-opacity:1.0;opacity:1.0;stroke:rgb(100%,100%,100%);stroke-opacity:1.0;stroke-width:1.0\" width=\"7.6666666666666856\" x=\"182.66666666666666\" y=\"215.62321598477641\"></rect><rect class=\"toyplot-Datum\" height=\"25.118934348239776\" style=\"fill:rgb(98.8%,55.3%,38.4%);fill-opacity:1.0;opacity:1.0;stroke:rgb(100%,100%,100%);stroke-opacity:1.0;stroke-width:1.0\" width=\"7.6666666666666572\" x=\"190.33333333333334\" y=\"214.88106565176022\"></rect><rect class=\"toyplot-Datum\" height=\"24.319695528068507\" style=\"fill:rgb(98.8%,55.3%,38.4%);fill-opacity:1.0;opacity:1.0;stroke:rgb(100%,100%,100%);stroke-opacity:1.0;stroke-width:1.0\" width=\"7.6666666666666572\" x=\"198.0\" y=\"215.68030447193149\"></rect><rect class=\"toyplot-Datum\" height=\"24.262607040913394\" style=\"fill:rgb(98.8%,55.3%,38.4%);fill-opacity:1.0;opacity:1.0;stroke:rgb(100%,100%,100%);stroke-opacity:1.0;stroke-width:1.0\" width=\"7.6666666666666572\" x=\"205.66666666666666\" y=\"215.73739295908661\"></rect><rect class=\"toyplot-Datum\" height=\"24.205518553758338\" style=\"fill:rgb(98.8%,55.3%,38.4%);fill-opacity:1.0;opacity:1.0;stroke:rgb(100%,100%,100%);stroke-opacity:1.0;stroke-width:1.0\" width=\"7.6666666666666856\" x=\"213.33333333333331\" y=\"215.79448144624166\"></rect><rect class=\"toyplot-Datum\" height=\"22.98763082778305\" style=\"fill:rgb(98.8%,55.3%,38.4%);fill-opacity:1.0;opacity:1.0;stroke:rgb(100%,100%,100%);stroke-opacity:1.0;stroke-width:1.0\" width=\"7.6666666666666572\" x=\"221.0\" y=\"217.01236917221695\"></rect><rect class=\"toyplot-Datum\" height=\"23.748810656517605\" style=\"fill:rgb(98.8%,55.3%,38.4%);fill-opacity:1.0;opacity:1.0;stroke:rgb(100%,100%,100%);stroke-opacity:1.0;stroke-width:1.0\" width=\"7.6666666666666856\" x=\"228.66666666666666\" y=\"216.25118934348239\"></rect><rect class=\"toyplot-Datum\" height=\"22.85442435775451\" style=\"fill:rgb(98.8%,55.3%,38.4%);fill-opacity:1.0;opacity:1.0;stroke:rgb(100%,100%,100%);stroke-opacity:1.0;stroke-width:1.0\" width=\"7.6666666666666572\" x=\"236.33333333333334\" y=\"217.14557564224549\"></rect><rect class=\"toyplot-Datum\" height=\"23.463368220742154\" style=\"fill:rgb(98.8%,55.3%,38.4%);fill-opacity:1.0;opacity:1.0;stroke:rgb(100%,100%,100%);stroke-opacity:1.0;stroke-width:1.0\" width=\"7.6666666666666856\" x=\"244.0\" y=\"216.53663177925785\"></rect><rect class=\"toyplot-Datum\" height=\"22.568981921979088\" style=\"fill:rgb(98.8%,55.3%,38.4%);fill-opacity:1.0;opacity:1.0;stroke:rgb(100%,100%,100%);stroke-opacity:1.0;stroke-width:1.0\" width=\"7.6666666666666288\" x=\"251.66666666666669\" y=\"217.43101807802091\"></rect><rect class=\"toyplot-Datum\" height=\"22.226450999048524\" style=\"fill:rgb(98.8%,55.3%,38.4%);fill-opacity:1.0;opacity:1.0;stroke:rgb(100%,100%,100%);stroke-opacity:1.0;stroke-width:1.0\" width=\"7.6666666666666856\" x=\"259.33333333333331\" y=\"217.77354900095148\"></rect><rect class=\"toyplot-Datum\" height=\"21.274976213130373\" style=\"fill:rgb(98.8%,55.3%,38.4%);fill-opacity:1.0;opacity:1.0;stroke:rgb(100%,100%,100%);stroke-opacity:1.0;stroke-width:1.0\" width=\"7.6666666666666856\" x=\"267.0\" y=\"218.72502378686963\"></rect><rect class=\"toyplot-Datum\" height=\"43.063748810656534\" style=\"fill:rgb(98.8%,55.3%,38.4%);fill-opacity:1.0;opacity:1.0;stroke:rgb(100%,100%,100%);stroke-opacity:1.0;stroke-width:1.0\" width=\"7.6666666666666288\" x=\"274.66666666666669\" y=\"196.93625118934347\"></rect></g></g></g><g class=\"toyplot-coordinates\" style=\"visibility:hidden\"><rect height=\"14.0\" style=\"fill:rgb(100%,100%,100%);fill-opacity:1.0;opacity:0.75;stroke:none\" width=\"90.0\" x=\"200.0\" y=\"60.0\"></rect><text style=\"alignment-baseline:middle;font-size:10px;font-weight:normal;stroke:none;text-anchor:middle\" x=\"245.0\" y=\"67.0\"></text></g><line style=\"\" x1=\"60.0\" x2=\"282.3333333333333\" y1=\"250.0\" y2=\"250.0\"></line><g><text style=\"alignment-baseline:middle;baseline-shift:-80%;font-size:10px;font-weight:normal;stroke:none;text-anchor:middle\" x=\"60.0\" y=\"250.0\">0</text><text style=\"alignment-baseline:middle;baseline-shift:-80%;font-size:10px;font-weight:normal;stroke:none;text-anchor:middle\" x=\"136.66666666666666\" y=\"250.0\">10</text><text style=\"alignment-baseline:middle;baseline-shift:-80%;font-size:10px;font-weight:normal;stroke:none;text-anchor:middle\" x=\"213.33333333333331\" y=\"250.0\">20</text><text style=\"alignment-baseline:middle;baseline-shift:-80%;font-size:10px;font-weight:normal;stroke:none;text-anchor:middle\" x=\"290.0\" y=\"250.0\">30</text></g><text style=\"alignment-baseline:middle;baseline-shift:-200%;font-weight:bold;stroke:none;text-anchor:middle\" x=\"175.0\" y=\"250.0\">Depth of coverage (N reads)</text><line style=\"\" x1=\"50.0\" x2=\"50.0\" y1=\"60.0\" y2=\"240.0\"></line><g><text style=\"alignment-baseline:middle;baseline-shift:80%;font-size:10px;font-weight:normal;stroke:none;text-anchor:middle\" transform=\"rotate(-90, 50.0, 240.0)\" x=\"50.0\" y=\"240.0\">0</text><text style=\"alignment-baseline:middle;baseline-shift:80%;font-size:10px;font-weight:normal;stroke:none;text-anchor:middle\" transform=\"rotate(-90, 50.0, 182.91151284490962)\" x=\"50.0\" y=\"182.91151284490962\">3000</text><text style=\"alignment-baseline:middle;baseline-shift:80%;font-size:10px;font-weight:normal;stroke:none;text-anchor:middle\" transform=\"rotate(-90, 50.0, 125.82302568981922)\" x=\"50.0\" y=\"125.82302568981922\">6000</text><text style=\"alignment-baseline:middle;baseline-shift:80%;font-size:10px;font-weight:normal;stroke:none;text-anchor:middle\" transform=\"rotate(-90, 50.0, 68.73453853472884)\" x=\"50.0\" y=\"68.73453853472884\">9000</text></g><text style=\"alignment-baseline:middle;baseline-shift:200%;font-weight:bold;stroke:none;text-anchor:middle\" transform=\"rotate(-90, 50.0, 150.0)\" x=\"50.0\" y=\"150.0\">N loci</text><text style=\"alignment-baseline:middle;baseline-shift:100%;font-size:14px;font-weight:bold;stroke:none;text-anchor:middle\" x=\"175.0\" y=\"50.0\">dataset3/sample=AR_re</text></g></svg><div class=\"toyplot-controls\"><ul class=\"toyplot-mark-popup\" onmouseleave=\"this.style.visibility='hidden'\" style=\"background:rgba(0%,0%,0%,0.75);border:0;border-radius:6px;color:white;cursor:default;list-style:none;margin:0;padding:5px;position:fixed;visibility:hidden\"><li class=\"toyplot-mark-popup-title\" style=\"color:lightgray;cursor:default;padding:5px;list-style:none;margin:0;\"></li><li class=\"toyplot-mark-popup-save-csv\" onmouseout=\"this.style.color='white';this.style.background='steelblue'\" onmouseover=\"this.style.color='steelblue';this.style.background='white'\" style=\"border-radius:3px;padding:5px;list-style:none;margin:0;\">Save as .csv</li></ul><script>\n",
       "(function()\n",
       "{\n",
       "  if(window.CSS !== undefined && window.CSS.supports !== undefined)\n",
       "  {\n",
       "    if(!window.CSS.supports(\"alignment-baseline\", \"middle\"))\n",
       "    {\n",
       "      var re = /\\s*alignment-baseline\\s*:\\s*([^;\\s]*)\\s*/;\n",
       "      var text = document.querySelectorAll(\"#t1082bb5d7c604b1fa934d98fbaf2df7d text\");\n",
       "      for(var i = 0; i != text.length; ++i)\n",
       "      {\n",
       "        var match = re.exec(text[i].attributes.style.value);\n",
       "        if(match)\n",
       "        {\n",
       "          if(match[1] == \"middle\")\n",
       "          {\n",
       "            var style = getComputedStyle(text[i]);\n",
       "            var font_size = style.fontSize.substr(0, style.fontSize.length - 2);\n",
       "            var dy = text[i].dy.baseVal.length ? text[i].dy.baseVal[0].value : 0;\n",
       "            dy += 0.4 * font_size;\n",
       "            text[i].setAttribute(\"dy\", dy);\n",
       "          }\n",
       "        }\n",
       "      }\n",
       "    }\n",
       "    if(!window.CSS.supports(\"baseline-shift\", \"0\"))\n",
       "    {\n",
       "      var re = /\\s*baseline-shift\\s*:\\s*([^;\\s]*)\\s*/;\n",
       "      var text = document.querySelectorAll(\"#t1082bb5d7c604b1fa934d98fbaf2df7d text\");\n",
       "      for(var i = 0; i != text.length; ++i)\n",
       "      {\n",
       "        var match = re.exec(text[i].attributes.style.value);\n",
       "        if(match)\n",
       "        {\n",
       "          var style = getComputedStyle(text[i]);\n",
       "          var font_size = style.fontSize.substr(0, style.fontSize.length - 2);\n",
       "          var percent = 0.01 * match[1].substr(0, match[1].length-1);\n",
       "          var dy = text[i].dy.baseVal.length ? text[i].dy.baseVal[0].value : 0;\n",
       "          dy -= percent * font_size\n",
       "          text[i].setAttribute(\"dy\", dy);\n",
       "        }\n",
       "      }\n",
       "    }\n",
       "  }\n",
       "})();\n",
       "</script><script>\n",
       "(function()\n",
       "{\n",
       "  var data_tables = [{\"data\": [[0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28], [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29], [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], [0.0, 9459.0, 3993.0, 2605.0, 2111.0, 2051.0, 1903.0, 1817.0, 1761.0, 1680.0, 1659.0, 1478.0, 1478.0, 1456.0, 1396.0, 1326.0, 1281.0, 1320.0, 1278.0, 1275.0, 1272.0, 1208.0, 1248.0, 1201.0, 1233.0, 1186.0, 1168.0, 1118.0, 2263.0]], \"title\": \"Bar Data\", \"names\": [\"left\", \"right\", \"baseline\", \"magnitude0\"], \"id\": \"te2ee6288007b4596ab930fed0d1155a2\", \"filename\": \"toyplot\"}, {\"data\": [[0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28], [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29], [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1903.0, 1817.0, 1761.0, 1680.0, 1659.0, 1478.0, 1478.0, 1456.0, 1396.0, 1326.0, 1281.0, 1320.0, 1278.0, 1275.0, 1272.0, 1208.0, 1248.0, 1201.0, 1233.0, 1186.0, 1168.0, 1118.0, 2263.0]], \"title\": \"Bar Data\", \"names\": [\"left\", \"right\", \"baseline\", \"magnitude0\"], \"id\": \"t5fd1c6e1bcf34c77aace5defeb3bff5a\", \"filename\": \"toyplot\"}];\n",
       "\n",
       "  function save_csv(data_table)\n",
       "  {\n",
       "    var uri = \"data:text/csv;charset=utf-8,\";\n",
       "    uri += data_table.names.join(\",\") + \"\\n\";\n",
       "    for(var i = 0; i != data_table.data[0].length; ++i)\n",
       "    {\n",
       "      for(var j = 0; j != data_table.data.length; ++j)\n",
       "      {\n",
       "        if(j)\n",
       "          uri += \",\";\n",
       "        uri += data_table.data[j][i];\n",
       "      }\n",
       "      uri += \"\\n\";\n",
       "    }\n",
       "    uri = encodeURI(uri);\n",
       "\n",
       "    var link = document.createElement(\"a\");\n",
       "    if(typeof link.download != \"undefined\")\n",
       "    {\n",
       "      link.href = uri;\n",
       "      link.style = \"visibility:hidden\";\n",
       "      link.download = data_table.filename + \".csv\";\n",
       "\n",
       "      document.body.appendChild(link);\n",
       "      link.click();\n",
       "      document.body.removeChild(link);\n",
       "    }\n",
       "    else\n",
       "    {\n",
       "      window.open(uri);\n",
       "    }\n",
       "  }\n",
       "\n",
       "  function open_popup(data_table)\n",
       "  {\n",
       "    return function(e)\n",
       "    {\n",
       "      var popup = document.querySelector(\"#t1082bb5d7c604b1fa934d98fbaf2df7d .toyplot-mark-popup\");\n",
       "      popup.querySelector(\".toyplot-mark-popup-title\").innerHTML = data_table.title;\n",
       "      popup.querySelector(\".toyplot-mark-popup-save-csv\").onclick = function() { popup.style.visibility = \"hidden\"; save_csv(data_table); }\n",
       "      popup.style.left = (e.clientX - 50) + \"px\";\n",
       "      popup.style.top = (e.clientY - 20) + \"px\";\n",
       "      popup.style.visibility = \"visible\";\n",
       "      e.stopPropagation();\n",
       "      e.preventDefault();\n",
       "    }\n",
       "\n",
       "  }\n",
       "\n",
       "  for(var i = 0; i != data_tables.length; ++i)\n",
       "  {\n",
       "    var data_table = data_tables[i];\n",
       "    var event_target = document.querySelector(\"#\" + data_table.id);\n",
       "    event_target.oncontextmenu = open_popup(data_table);\n",
       "  }\n",
       "})();\n",
       "</script><script>\n",
       "(function()\n",
       "{\n",
       "  var axes = {\"t551d9eff8bca4b2b81d94d8a92ad895e\": {\"x\": [{\"domain\": {\"bounds\": {\"max\": Infinity, \"min\": -Infinity}, \"max\": 30.0, \"min\": 0}, \"range\": {\"bounds\": {\"max\": Infinity, \"min\": -Infinity}, \"max\": 290.0, \"min\": 60.0}, \"scale\": \"linear\"}], \"y\": [{\"domain\": {\"bounds\": {\"max\": Infinity, \"min\": -Infinity}, \"max\": 9459.0, \"min\": 0.0}, \"range\": {\"bounds\": {\"max\": -Infinity, \"min\": Infinity}, \"max\": 60.0, \"min\": 240.0}, \"scale\": \"linear\"}]}};\n",
       "\n",
       "  function sign(x)\n",
       "  {\n",
       "    return x < 0 ? -1 : x > 0 ? 1 : 0;\n",
       "  }\n",
       "\n",
       "  function _mix(a, b, amount)\n",
       "  {\n",
       "    return ((1.0 - amount) * a) + (amount * b);\n",
       "  }\n",
       "\n",
       "  function _log(x, base)\n",
       "  {\n",
       "    return Math.log(Math.abs(x)) / Math.log(base);\n",
       "  }\n",
       "\n",
       "  function _in_range(a, x, b)\n",
       "  {\n",
       "    var left = Math.min(a, b);\n",
       "    var right = Math.max(a, b);\n",
       "    return left <= x && x <= right;\n",
       "  }\n",
       "\n",
       "  function to_domain(projection, range)\n",
       "  {\n",
       "    for(var i = 0; i != projection.length; ++i)\n",
       "    {\n",
       "      var segment = projection[i];\n",
       "      if(_in_range(segment.range.bounds.min, range, segment.range.bounds.max))\n",
       "      {\n",
       "        if(segment.scale == \"linear\")\n",
       "        {\n",
       "          var amount = (range - segment.range.min) / (segment.range.max - segment.range.min);\n",
       "          return _mix(segment.domain.min, segment.domain.max, amount)\n",
       "        }\n",
       "        else if(segment.scale[0] == \"log\")\n",
       "        {\n",
       "          var amount = (range - segment.range.min) / (segment.range.max - segment.range.min);\n",
       "          var base = segment.scale[1];\n",
       "          return sign(segment.domain.min) * Math.pow(base, _mix(_log(segment.domain.min, base), _log(segment.domain.max, base), amount));\n",
       "        }\n",
       "      }\n",
       "    }\n",
       "  }\n",
       "\n",
       "  // Compute mouse coordinates relative to a DOM object, with thanks to d3js.org, where this code originated.\n",
       "  function d3_mousePoint(container, e)\n",
       "  {\n",
       "    if (e.changedTouches) e = e.changedTouches[0];\n",
       "    var svg = container.ownerSVGElement || container;\n",
       "    if (svg.createSVGPoint) {\n",
       "      var point = svg.createSVGPoint();\n",
       "      point.x = e.clientX, point.y = e.clientY;\n",
       "      point = point.matrixTransform(container.getScreenCTM().inverse());\n",
       "      return [point.x, point.y];\n",
       "    }\n",
       "    var rect = container.getBoundingClientRect();\n",
       "    return [e.clientX - rect.left - container.clientLeft, e.clientY - rect.top - container.clientTop];\n",
       "  };\n",
       "\n",
       "  function display_coordinates(e)\n",
       "  {\n",
       "    var dom_axes = e.currentTarget.parentElement;\n",
       "    var data = axes[dom_axes.id];\n",
       "\n",
       "    point = d3_mousePoint(e.target, e);\n",
       "    var x = Number(to_domain(data[\"x\"], point[0])).toFixed(2);\n",
       "    var y = Number(to_domain(data[\"y\"], point[1])).toFixed(2);\n",
       "\n",
       "    var coordinates = dom_axes.querySelectorAll(\".toyplot-coordinates\");\n",
       "    for(var i = 0; i != coordinates.length; ++i)\n",
       "    {\n",
       "      coordinates[i].style.visibility = \"visible\";\n",
       "      coordinates[i].querySelector(\"text\").textContent = \"x=\" + x + \" y=\" + y;\n",
       "    }\n",
       "  }\n",
       "\n",
       "  function clear_coordinates(e)\n",
       "  {\n",
       "    var dom_axes = e.currentTarget.parentElement;\n",
       "    var coordinates = dom_axes.querySelectorAll(\".toyplot-coordinates\");\n",
       "    for(var i = 0; i != coordinates.length; ++i)\n",
       "      coordinates[i].style.visibility = \"hidden\";\n",
       "  }\n",
       "\n",
       "  for(var axes_id in axes)\n",
       "  {\n",
       "    var event_target = document.querySelector(\"#\" + axes_id + \" .toyplot-coordinate-events\");\n",
       "    event_target.onmousemove = display_coordinates;\n",
       "    event_target.onmouseout = clear_coordinates;\n",
       "  }\n",
       "})();\n",
       "</script></div></div>"
      ],
      "text/plain": [
       "<toyplot.canvas.Canvas at 0x7f5de591d510>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import toyplot\n",
    "import toyplot.svg\n",
    "import numpy as np\n",
    "\n",
    "## read in the depth information for this sample\n",
    "with open(\"empirical_3/clust.85/AR_re.depths\", 'rb') as indat:\n",
    "    depths = np.array(indat.read().strip().split(\",\"), dtype=int)\n",
    "    \n",
    "## make a barplot in Toyplot\n",
    "canvas = toyplot.Canvas(width=350, height=300)\n",
    "axes = canvas.axes(xlabel=\"Depth of coverage (N reads)\", \n",
    "                   ylabel=\"N loci\", \n",
    "                   label=\"dataset3/sample=AR_re\")\n",
    "\n",
    "## select the loci with depth > 5 (kept)\n",
    "keeps = depths[depths>5]\n",
    "\n",
    "## plot kept and discarded loci\n",
    "edat = np.histogram(depths, range(30)) # density=True)\n",
    "kdat = np.histogram(keeps, range(30)) #, density=True)\n",
    "axes.bars(edat)\n",
    "axes.bars(kdat)\n",
    "\n",
    "#toyplot.svg.render(canvas, \"empirical_3_depthplot.svg\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Print final stats table"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {
    "collapsed": false,
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\r\n",
      "\r\n",
      "90871       ## loci with > minsp containing data\r\n",
      "87969       ## loci with > minsp containing data & paralogs removed\r\n",
      "87969       ## loci with > minsp containing data & paralogs removed & final filtering\r\n",
      "\r\n",
      "## number of loci recovered in final data set for each taxon.\r\n",
      "taxon\tnloci\r\n",
      "AR_re     \t54176\r\n",
      "BJSB3_v   \t38030\r\n",
      "BJSL25_me \t61596\r\n",
      "BJVL19_me \t60159\r\n",
      "BZBB1_v   \t35466\r\n",
      "CH_re     \t51058\r\n",
      "CRL0001_me\t62704\r\n",
      "CRL0030_me\t60181\r\n",
      "CUCA4_v   \t25875\r\n",
      "CUMM5_v   \t2520\r\n",
      "CUSV6_v   \t34898\r\n",
      "CUVN10_v  \t43516\r\n",
      "DO_re     \t46865\r\n",
      "DU_re     \t52925\r\n",
      "EN_re     \t28342\r\n",
      "FLAB109_me\t53770\r\n",
      "FLBA140_me\t59903\r\n",
      "FLCK18_v  \t35834\r\n",
      "FLCK216_v \t33674\r\n",
      "FLMO62_v  \t37431\r\n",
      "FLSA185_v \t29650\r\n",
      "FLSF33_v  \t35589\r\n",
      "FLSF47_v  \t37411\r\n",
      "FLSF54_me \t60404\r\n",
      "FLWO6_v   \t35292\r\n",
      "HE_re     \t42687\r\n",
      "HNDA09_v  \t43262\r\n",
      "LALC2_v   \t48666\r\n",
      "MXED8_v   \t41094\r\n",
      "MXGT4_v   \t47152\r\n",
      "MXSA3017_v\t40169\r\n",
      "NI_re     \t44506\r\n",
      "SCCU3_v   \t18442\r\n",
      "TXGR3_v   \t34111\r\n",
      "TXMD3_v   \t37451\r\n",
      "TXWV2_v   \t5109\r\n",
      "\r\n",
      "\r\n",
      "## nloci = number of loci with data for exactly ntaxa\r\n",
      "## ntotal = number of loci for which at least ntaxa have data\r\n",
      "ntaxa\tnloci\tsaved\tntotal\r\n",
      "1\t-\r\n",
      "2\t-\t\t-\r\n",
      "3\t-\t\t-\r\n",
      "4\t7085\t*\t87969\r\n",
      "5\t5234\t*\t80884\r\n",
      "6\t4497\t*\t75650\r\n",
      "7\t3579\t*\t71153\r\n",
      "8\t3074\t*\t67574\r\n",
      "9\t2755\t*\t64500\r\n",
      "10\t2629\t*\t61745\r\n",
      "11\t2561\t*\t59116\r\n",
      "12\t2409\t*\t56555\r\n",
      "13\t2434\t*\t54146\r\n",
      "14\t2477\t*\t51712\r\n",
      "15\t2302\t*\t49235\r\n",
      "16\t2464\t*\t46933\r\n",
      "17\t2394\t*\t44469\r\n",
      "18\t2444\t*\t42075\r\n",
      "19\t2517\t*\t39631\r\n",
      "20\t2580\t*\t37114\r\n",
      "21\t2652\t*\t34534\r\n",
      "22\t2775\t*\t31882\r\n",
      "23\t2944\t*\t29107\r\n",
      "24\t2830\t*\t26163\r\n",
      "25\t2960\t*\t23333\r\n",
      "26\t3060\t*\t20373\r\n",
      "27\t3228\t*\t17313\r\n",
      "28\t3093\t*\t14085\r\n",
      "29\t3047\t*\t10992\r\n",
      "30\t2766\t*\t7945\r\n",
      "31\t2285\t*\t5179\r\n",
      "32\t1581\t*\t2894\r\n",
      "33\t941\t*\t1313\r\n",
      "34\t310\t*\t372\r\n",
      "35\t59\t*\t62\r\n",
      "36\t3\t*\t3\r\n",
      "\r\n",
      "\r\n",
      "## nvar = number of loci containing n variable sites (pis+autapomorphies).\r\n",
      "## sumvar = sum of variable sites (SNPs).\r\n",
      "## pis = number of loci containing n parsimony informative sites.\r\n",
      "## sumpis = sum of parsimony informative sites.\r\n",
      "\tnvar\tsumvar\tPIS\tsumPIS\r\n",
      "0\t4738\t0\t14177\t0\r\n",
      "1\t4963\t4963\t12086\t12086\r\n",
      "2\t5663\t16289\t11896\t35878\r\n",
      "3\t6227\t34970\t10819\t68335\r\n",
      "4\t6540\t61130\t9289\t105491\r\n",
      "5\t6974\t96000\t8034\t145661\r\n",
      "6\t6661\t135966\t6292\t183413\r\n",
      "7\t6482\t181340\t4703\t216334\r\n",
      "8\t6305\t231780\t3517\t244470\r\n",
      "9\t5622\t282378\t2420\t266250\r\n",
      "10\t5119\t333568\t1656\t282810\r\n",
      "11\t4321\t381099\t1104\t294954\r\n",
      "12\t3815\t426879\t772\t304218\r\n",
      "13\t3219\t468726\t507\t310809\r\n",
      "14\t2572\t504734\t297\t314967\r\n",
      "15\t2151\t536999\t165\t317442\r\n",
      "16\t1660\t563559\t107\t319154\r\n",
      "17\t1261\t584996\t61\t320191\r\n",
      "18\t998\t602960\t33\t320785\r\n",
      "19\t753\t617267\t11\t320994\r\n",
      "20\t559\t628447\t10\t321194\r\n",
      "21\t402\t636889\t4\t321278\r\n",
      "22\t287\t643203\t5\t321388\r\n",
      "23\t205\t647918\t3\t321457\r\n",
      "24\t134\t651134\t0\t321457\r\n",
      "25\t112\t653934\t1\t321482\r\n",
      "26\t73\t655832\t0\t321482\r\n",
      "27\t64\t657560\t0\t321482\r\n",
      "28\t33\t658484\t0\t321482\r\n",
      "29\t14\t658890\t0\t321482\r\n",
      "30\t16\t659370\t0\t321482\r\n",
      "31\t7\t659587\t0\t321482\r\n",
      "32\t3\t659683\t0\t321482\r\n",
      "33\t5\t659848\t0\t321482\r\n",
      "34\t6\t660052\t0\t321482\r\n",
      "35\t2\t660122\t0\t321482\r\n",
      "36\t1\t660158\t0\t321482\r\n",
      "37\t1\t660195\t0\t321482\r\n",
      "38\t0\t660195\t0\t321482\r\n",
      "39\t1\t660234\t0\t321482\r\n",
      "total var= 660234\r\n",
      "total pis= 321482\r\n",
      "sampled unlinked SNPs= 83231\r\n",
      "sampled unlinked bi-allelic SNPs= 51854\r\n"
     ]
    }
   ],
   "source": [
    "cat empirical_3/stats/empirical_3_m4.stats\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "\n",
      "145112      ## loci with > minsp containing data\n",
      "142210      ## loci with > minsp containing data & paralogs removed\n",
      "142210      ## loci with > minsp containing data & paralogs removed & final filtering\n",
      "\n",
      "## number of loci recovered in final data set for each taxon.\n",
      "taxon\tnloci\n",
      "AR_re     \t58912\n",
      "BJSB3_v   \t39297\n"
     ]
    }
   ],
   "source": [
    "%%bash\n",
    "head -n 10 empirical_3/stats/empirical_3_m2.stats\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Infer ML phylogeny in _raxml_ as an unrooted tree"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "%%bash\n",
    "## raxml argumement w/ ...\n",
    "raxmlHPC-PTHREADS-AVX -f a -m GTRGAMMA -N 100 -x 12345 -p 12345 -T 20 \\\n",
    "                      -w /home/deren/Documents/RADmissing/empirical_3/ \\\n",
    "                      -n empirical_3_m4 -s empirical_3/outfiles/empirical_3_m4.phy\n",
    "                      "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "%%bash\n",
    "## raxml argumement w/ ...\n",
    "raxmlHPC-PTHREADS-AVX -f a -m GTRGAMMA -N 100 -x 12345 -p 12345 -T 20 \\\n",
    "                      -w /home/deren/Documents/RADmissing/empirical_3/ \\\n",
    "                      -n empirical_3_m2 -s empirical_3/outfiles/empirical_3_m2.phy\n",
    "                      "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {
    "collapsed": false,
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "\n",
      "This is RAxML version 8.0.16 released by Alexandros Stamatakis on March 21 2014.\n",
      "\n",
      "With greatly appreciated code contributions by:\n",
      "Andre Aberer      (HITS)\n",
      "Simon Berger      (HITS)\n",
      "Alexey Kozlov     (HITS)\n",
      "Nick Pattengale   (Sandia)\n",
      "Wayne Pfeiffer    (SDSC)\n",
      "Akifumi S. Tanabe (NRIFS)\n",
      "David Dao         (KIT)\n",
      "Charlie Taylor    (UF)\n",
      "\n",
      "\n",
      "Alignment has 612971 distinct alignment patterns\n",
      "\n",
      "Proportion of gaps and completely undetermined characters in this alignment: 53.74%\n",
      "\n",
      "RAxML rapid bootstrapping and subsequent ML search\n"
     ]
    }
   ],
   "source": [
    "%%bash\n",
    "head -n 20 empirical_3/RAxML_info.empirical_3_m4"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "head: cannot open ‘empirical_3/RAxML_info.empirical_3_m2’ for reading: No such file or directory\n"
     ]
    }
   ],
   "source": [
    "%%bash\n",
    "head -n 20 empirical_3/RAxML_info.empirical_3_m2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Plot the tree in R using `ape`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "%load_ext rpy2.ipython"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZAAAAMgCAMAAAA6NMSiAAADAFBMVEUAAAABAQECAgIDAwMEBAQF\nBQUGBgYHBwcICAgJCQkKCgoLCwsMDAwNDQ0ODg4PDw8QEBARERESEhITExMUFBQVFRUWFhYXFxcY\nGBgZGRkaGhobGxscHBwdHR0eHh4fHx8gICAhISEiIiIjIyMkJCQlJSUmJiYnJycoKCgpKSkqKior\nKyssLCwtLS0uLi4vLy8wMDAxMTEyMjIzMzM0NDQ1NTU2NjY3Nzc4ODg5OTk6Ojo7Ozs8PDw9PT0+\nPj4/Pz9AQEBBQUFCQkJDQ0NERERFRUVGRkZHR0dISEhJSUlKSkpLS0tMTExNTU1OTk5PT09QUFBR\nUVFSUlJTU1NUVFRVVVVWVlZXV1dYWFhZWVlaWlpbW1tcXFxdXV1eXl5fX19gYGBhYWFiYmJjY2Nk\nZGRlZWVmZmZnZ2doaGhpaWlqampra2tsbGxtbW1ubm5vb29wcHBxcXFycnJzc3N0dHR1dXV2dnZ3\nd3d4eHh5eXl6enp7e3t8fHx9fX1+fn5/f3+AgICBgYGCgoKDg4OEhISFhYWGhoaHh4eIiIiJiYmK\nioqLi4uMjIyNjY2Ojo6Pj4+QkJCRkZGSkpKTk5OUlJSVlZWWlpaXl5eYmJiZmZmampqbm5ucnJyd\nnZ2enp6fn5+goKChoaGioqKjo6OkpKSlpaWmpqanp6eoqKipqamqqqqrq6usrKytra2urq6vr6+w\nsLCxsbGysrKzs7O0tLS1tbW2tra3t7e4uLi5ubm6urq7u7u8vLy9vb2+vr6/v7/AwMDBwcHCwsLD\nw8PExMTFxcXGxsbHx8fIyMjJycnKysrLy8vMzMzNzc3Ozs7Pz8/Q0NDR0dHS0tLT09PU1NTV1dXW\n1tbX19fY2NjZ2dna2trb29vc3Nzd3d3e3t7f39/g4ODh4eHi4uLj4+Pk5OTl5eXm5ubn5+fo6Ojp\n6enq6urr6+vs7Ozt7e3u7u7v7+/w8PDx8fHy8vLz8/P09PT19fX29vb39/f4+Pj5+fn6+vr7+/v8\n/Pz9/f3+/v7////isF19AAAgAElEQVR4nO2dd2AUZfrHh2wSkiCJIQQIAhJaDIgEIiUq3uERiYA0\nBS5RigI5iicIKIdHOcVTigh4ooam3lFsIKhIye8oUkMA5QIKQrBEWoAESIfsPr95Z2dm26SR3dkn\n8z6fP6a8O7v7mo+7S/ab53kFIFAheHsChCMkBBkkBBkkBBkkBBkkBBkkBBkkBBkkBBkkBBkkBBkk\nBBkkBBkkBBkkBBkkBBkkBBkkBBkkBBkkBBkkBBkkBBkkBBkkBBkkBBkkBBkkBBkkBBkkBBkkBBkk\nBBkkBBkkBBkkBBkkBBkkBBkkBBkkBBkkBBkkBBkkBBkkBBkkBBkkBBkkBBkkBBkkBBkkBBkkBBkk\nBBkkBBkkBBkkBBkkBBkkBBkkBBkkBBkkBBkkBBkkBBkkBBkkBBkkBBkkBBkkBBkkBBkkBBkkBBkk\nBBkkBBkkBBkkBBkkBBkkBBkkBBkkBBkkBBkkBBmeFtLRxyRy6lbAJPFkfbyHn80AeFhIke85aX/M\nr14RwPTpnn02I+BhIWmNrPuVT7ReAxD/uWefzQh4WMjS3tb9+LkL/gCWer+w4wmzOr717cP1m3/k\n2aeuobhHSFotQZNZz0i7E9Al9XLtk2fDLOzi2IcP7+x63LwrwOyW5zYY7hGyXxgxV4u09l9Kt98M\nyIHEKZ/1YscFtX+HeKYpzC1PbTTcIyRd+EpzvMDngrT/rgXArrBJL7PjPV0Bgi8BFBa55amNhmeF\n7Gti3a8YAmCJCtzAjhe8AHDfa8WH7trnlqc2Gp4VskRgv4X0grHzxJNFwq9s7IlPAY52qtvhC7c8\ns+HwrBCiypAQZJAQZJAQZJAQZJAQZJAQZJAQZJAQZOggxBJXm32v29nH5BOywC1PZ2R0ELI+zv8c\nQIl/Fli+qeuWpzMynhdyM/pwqz0ARxqLxxkt3fJ0RsZdQqZ9qs1v8N4o+NNHACkDAK4/O0u6nDLD\nsnGPkI3agaHI6LyoizBqNsAYdtbeGhJSZlg27nqFLM/UJv8VJmI4QMet4mWBt9jVlBmWg6c/Qy41\nuwHw4UNQ5HuJPZuUElJmWA6eFjJhtrjZ3hjSpOzQKoQyw3LwsJDToVfF7QmhcGk/6dmK2ZYyw3Kg\n39SRQUKQQUKQQUKQQUKQQUKQQUKQQUKQQUKQgUhIgY/JJ7A3K4FrdkjcvNOv+g9ZA0EkZF8zgOyh\nAwCu+LNvWIbPqf5D1kAQCVkySNycC7JA6v3sNHpr9R+yBqKvkBn39CyLXv837HXximt1LDB/nHhw\nw+cKuwd34aK+QmIDYsvkg+ht4hXfxwAkrhIPdraw3oO3cFFfIYlRZd9mfUnMHyu+WWWwg6FskL9w\nEY+Q3ZHiJjcyE8x+heLRIx+yQf7CRTxCFg6Gm3s7sb+ki5yTnTMjhlnhMFzEIyRRMNXuvk56tG5B\nzZ67LA3yFy7iEUJIkBBkkBBkeFTI2sDmLRy44063PJ2R8aiQxcLQZAdaNXTL0xkZjwpZLfzkOEBv\nWRVCQpBBQpCBS4gto4L9T4Q3faHALbOrUeASomZUsDpqc/6pbhPdMrsaBS4hakb1WzPW+exoI7fM\nrkahs5DWOeVRqmZUY95kl+f7sRSEr4xKXyH3lln6JtFdzaiiT7DLM9qyLV8Zlb5CPn5Ws1emwsdK\nRmU2SR/n80YDdxmVvkIqQM2ooNH34tH5xmeBu4wKlRBbRjXpgTN5X7dcxgY5y6hQCbFlVAUTGzcY\nsEMa5CyjQiWEICHoICHIICHIICHIICHIICHIICHIICHIMKCQPb3rN1ks7iMOiptlj3lhBtXBeEJW\nt9lemB6yG875sG+Mx87UfwbVwnBCspuzpeCmLoRN97LTLpt0n0H1wCjkt5kpt8vq4qmz5EeZNVzc\n3AqQ/mKiBoWOGIX0Lz9XLJct938vP0pv6TxCOq5BoSNGIbOFM2W01KyQC+bAfPERStYVWhqy/GTF\n4+wBa1LoiFHIPwTLbc/EEnxS3C7pYcmSUuDkV9hgTQodjSYEhvW6WPjvxj/Axg7sLGYz29ak0NFw\nQq4mhYb3PwEwQ/oDCdNFNlaTQkfDCanpkBBkkBBkkBBkkBBkkBBkkBBkeEHI2Yq+GJlIQqpLVYSU\n1Kr4K0ISUk2qIqRQeKqM9aoUBtccIWoq2dHHJHKq+o/oDSFzK3i0mvMZoqaSRb7n3PSQJKQ6qKlk\nmtuqIUlIFfhjR8fGne+rqeTS3upF1UwnSUgVCAp3aNvZea6aSj4jHUiFkdVMJ0lIFYj4i+O5LZVs\n/6UyVt10koRUAWchaipZ4HNBGatuOklCqoCzEDWV3NdEHatuOklCqoCzEDWVXCKw30J6sbHqppMk\npAo4C/EEJKQKkBBk8CrkKaEFTkx/dMPPqgIwCvlPp2ScBPV1w8+qAjAKQQuvb1loISHIICHI0BBy\nK2CSuF0fD/D4Urc8BwmpAhpCjvnVKwKYPl28Mc0tz0FCqoCGkJVPtF4DEP85nPN1T60DCakCESOd\nunbmwvi5C/4Alnq/wJcd2RXVr50jIVXA3+XPY97rknq59smzYRaYNYZdUf3aORJSBeZMderauSgz\nIAcSp3zWC6BPCrildo6EVIvvWgDsCpv0MlgaHgG31M6RkGqxYgiAJSpwA/zuXwJuqZ1zj5CdQozW\ncqr3Chmu1xpKyNh54maR8Cts7MxO3VA75x4hP4Xdq7WaaqRwzPVaQwlxP+4RUgbGf8tyPyQEGSQE\nGSQEGSQEGSQEGSQEGSQEGVwIGWEyCT4m0wHffQBZjT8HOD68VXCvHwE6+5h8QhboPZ1y4UIIQLE/\nKzl7KhEKOi0E+LzNZ1dPJbeFEv8ssHxTF+D6qAYRr+g/Ky04EXJYavV31O/8kPEW+K05Kx7If6Lg\nSGNxn9ES4PGR2Zl+ONr/GU3I3AWunTE3A7wvtfqDHtF9bgE8+658ccoA8cXx7Cy49UXBhbcfYCOX\nRC35LX/w4AQrwmhCtIrg77TA6H9It04RssRt6zMAXQUhBsawW9uzV8YRQdjKLrAEn4OXp3twfhVi\nNCG+U13XqywAiJHyms+jOonPbTaxPpnwr2ehoyghPVB8zYAl9y1rzc39B35qk+/B+VWI4YTM0Bot\n9D0vbg9GnFwTUQzQlFX759/zryLfS+KBUARJhQCZkdKViR8/5t3Wy3wIOcg+vc9GbIPiBqsAZtz/\nw830HoH70qQXhSik/srS3xMWSVfO6v64B2dXCfgQ8o74U86NXiIezWxrhqIXGtWPP/Jg3tJ+7Dah\nGD5rHtphpbUmZXXAWQ/OrhLwIaQGQUKQQUKQYRwh0/qzIief8R55cP0wjpA7/UNFhIEeeXD9MI6Q\ndoPZlt6yyoOEVB0SggwSggw+hDhFhiX+faRR4WKBj8knsLe7+iW6Az6EOEWGRwJYr8QTdZrCvmYA\n2UMHIMoMeRHiEBmm9Au2AAzqPxCWDBJHzwVZ8GSGRhPiE6lRFvGBU2SYPL9hDqR3nvBPGPa6OHit\njgVPZmg0Ia06aghZCo6RYacdXb6D+O3dtkH0NvH0+xjAkxkaTUhZ2EeGxbWvDd64s8fNgCs3fK6I\no/PHAp7MkBMhDpFhehuYuiRu//eRsJvFhLmRmYgyQ06EOESG7yfBv9r1gxWDYeFguLm30wJAlBly\nIsQhMhy9EL6qdYzVByYKptrd14mjeDJDToTUHEgIMkgIMkgIMkgIMkgIMkgIMkgIMngSYsulspTi\nNnlxNUQ5FU9CbLmUklQpi6upOdW155uFzPLI9CoLT0JsuZSSVCmLq6k5VZ/XLv/o3ZzKmELMl1zL\ndnIstlxKSaqUxdWUnAru/CLn1Tg24rWcyphC7tdaz3W6LZdSkiplcTU1p5osnv6XPYLXcqqaLuRl\nW1HhvbbR3g1dSz9TvlNzKSWpUhZXU3OqpY+cvLzoHmnMWzlVTRfykjBlmpXwR2yjSW00LrXlUkpS\npSyupuZU0WcAfr5bGvNWTlXThbwhKA1A7d+yNIXYciklqVIWV1NzqgYfmH+Of0ca81ZOxZEQNZdS\nkyplcTU1p9rUql7sOu/mVBwJqRmQEGSQEGQYTsiur1NTUx9p4can0BfDCfGTficJdONT6IvhhIQM\npFdImXhDSL3ngT5DyoKEVB0Sggx+hGjkhUozTLniLcuN87pt+BHimheqzTCVijeAnRGlbpzc7cCP\nENe8UG2GqVS8we4GCW6c223hUSHvCUNdl1x+VpjmxqfQFtLPJQ9ZleeaF6rNMJWKt719k+ewvTfL\n2jwqZJXvnaEu3CmMcONTaAvp4JoXbnbNC9VmmHLF24H++V12sMfwZlmbR4VoosdbVtEF50A9t9gl\nL7Q1w7RWvKUNKigMtIaEXixrM6YQDVzzQrUZplzx1lN6JeWyi71Y1saNENe8UG2GKVe8iSxOtu69\nWNbGjRDXvFBthmmteGMkrbDuvVjWxo2QmgIJQQYJQQYJQQYJQYZxhLSMlb4jqTPGjY/uBYwjJET+\njiTOjY/uBYwj5OpPmYw7/+rGR/cCxhEiQ58hVQWDkF/ZH/VGS19XXfa5AYUvNKvTPR0soRrLv+sN\nn0LWdwbITgljh1vaAox/7vqFN1tAZsAt6db9caGRm8DSTvxAGunGqVYOPoVMZz3kfmnFDl8dAZaQ\nqwB5z8An1n8PnK+3xbyqEex5OC8vr9iNU60cnAh5u2GsHX0eWQYlGYmp7JbH3wFL6yFprAbhJeu/\nBz6ZDHD0Lng6Nrj5N+y8tP7zfSKXDeqaDPDtw/Wbf+TG2WvAiZCkWvZNMfvcyf59/Am7wdLoIMCZ\npDvuFu08ov6sT8etLei24eos6Vv5DL/jMPjRm9n+5v1dj5t3BXi2SJcTITNN9meZtUugdF19S7EZ\nsnyl+xfOaGKxhJyw3nx1Qpz10z0nkP3wlw8H6LYfjnSGeOYxzI2z14BLIdKHRXa4OfQ8vN8FzCbx\nIbIjLKfrWP8E6OvuG8Q3sA2ZAIceZOejlkFxcAm8+xwEXxKnX+T66O6ESyEvTgTIHTbZ7P/tydYp\nAJELr2U9+QZ83F26ccf9l/Py8mHIU/m/dpP+5KHd/+DAAwDDV8N9rxUfumufG2evAZdCetQymRpO\nLoIFwRGvii+GvTFBUW+VwlT2y0k3SGDvSw/Bqa4hsVJn5WshpbBoCkCb03C0U90OX7hx8lpwKQQz\nJAQZJAQZJAQZJAQZJAQZJAQZhhNyR+w0DR7wceNzehTjCdHqXSYItVzvqp1S2Y16BcMJ0UbrLUs7\npVJGvbVsG8dCtFMqZdRby7bxJCQn044siNdMqeRRry3bxpMQH4dPlWWhWimVRR310rJtHAmxCI/b\nl4Ee1EyplOzKa8u2cSXkH/ZD2imVPGrx2rJt/ArRTqnkUe8t28avEO2UShn12rJt/ApBCglBBglB\nBidC/uaTmrpd+Lsbn9hTcCKkt/S7YG83PrGn4ETI5U2HD6fTK0QTrwhh0GeINiSkXEgIMvgV4hgY\nfhzLDkZP8XpZG79CHAPDw/XE/emwK0pZG5TE/Qhwqnd4C53TXH6FOAaGucJ1gKTXlLK2W3+LCTXD\n1cgtpWtD9Y0NuRLSc66N446BIYQdg4zG+UpZW+FXUxMAlo8EuBDEhOhX1saREDDZB4bPOAWGXTbB\nwKV2ZW0vvgKwNCYr8zHpw0W/sjaehBTbNcW87lzWlrQ4vWUJqGVt8NA2gBsJQXGJI9iZfmVtPAlx\nwDEwBJg5qdda8WNdLmuDksBcgCwz3Irawk71K2vjVohTYAgf3R0jvg/JZW0AaawBecjG4ikJ0ruT\nfmVt3ApxCgxhn8BeCnJZG8CiUeLmg/DG466zi3Usa+NWCFZICDJICDJICDJICDJICDJICDJICDL4\nErKnd/0mi8V9s0Pi5p1+oNSxKfGUGyd2u3AlZHWb7YXpIbvhij9r3Td8jlrHpsRTkDukUdeVbpxe\n1eFJSHbzX8Tt1IWQej87jd6q1rEp8ZS5+0rznoFunF7V4UnI1Fnywfxx4uaGzxW1jk2Jp1J7yld4\nb502Ywrp0CrWlSfv/16+OXGVuNnJVpSWuy0q8dTTo9uFsdu8uE6bMYUENe3pyihpCbaSdYUQnSEe\nzB9qnc6MJhYlnoK2idlrgqXv2722TpsxhYRP0Bi0BJ8Ut0t6WMx+heLBIx+qdWxKPGUOuASXw6X3\nMK+t08aREBjW62LhvxuLHwyRc7JzZsQUqnVsSjxlCf2yZNxL0rVeW6eNJyFXk0LD+7PMPL1bULPn\nLoNax6bEU7C2UZMXrc2svbZOG09CagQkBBkkBBkkBBkkBBkkBBkkBBkkBBkkBBkkxCVGLPHvw4ZH\nCBdtR9WfaGUhIS4x4pGAe8WDE3Wa2h3pt2Qb90JcY8SUfsEWgEH9B9od6bdkm2GFnMisFMWuMWLy\n/IY5kN55wj/tjtQl2zweJRpVSIp2P2UXhrnGiJ12dPkO4rd322Y7si3Z5vEo0ShCej+ZbEfAk+uF\nl1MqwYoTLjFice1rgzfu7HEz4IrtiD2Ddck2j0eJRhEihLaww/THL4XDlbq/a4yY3gamLonb/32k\n3ZFtyTaPR4lGEVJrtv1o+ITKCnGNEd9Pgn+16wcrBtsd2ZZs83iUqL+Qgkq+u1eWQulRb1uIS4w4\neiF8VesYjJ1nd6Qu2eb5KFF/ITB82lz3kVDdVwg2vCDErVT7LQsbJAQZJAQZJAQZJAQZJAQZJAQZ\nJAQZPApxqTRUR8ASV9sLC4HZw6EQ10pDZQRgfZz/OeminRGlHptzufAnxDUiVEfgZvThVnvY6O4G\nCR6bcvkYWMjvQ9/Q+PIr0zUiVEfgvVHwJ9Z0cW/f5Dns3AulhgYW8prml8NvukaE6khe1EUYJT7Q\ngf75XaQv271QamhgIZ8IJzTuYHaJCG0jrzBjwyFtUEFhoDUV1L/U0JhC7owbIRzSFuIaEaojl5rd\nAPjwIegpvZRy2dX6lxoaU4i/+AP9SluIRqWhMjKBPcZ26W8ZFidbL9a/1NCYQnIvry3rLUuj0lAe\nOR16VTw5IUWQSSusF+tfamhMIQBlfoZgh4Qgg4Qgg4Qgg4Qgg4Qgg4Qgg4Qgg0shSh7V2cfkE7KA\ndQHy8tJsdvAoRMmjSvyzwPJNXQBlaba0mPofueys67XpBodC1DzqCPvaKqMlKEuz3Wj63bZw5511\nvTb94FCImkelDAC4/qx4Ji/NtmI8FJjMTjvrem2g20JtfAn5a3hsbK9YJY8aw75lb29Wl2Z79EvI\njHTeWddrA90WauNLSK/AwYNHqnlUx60A6YG31KXZIn+GLx513lnXawPdFmrjS0hyY7uEqsj3kngg\n3l9ems3sWwh/WeG0k9drA90WauNPiJpHpTVhQ+L95aXZLPVSj9x1w2knr9cGui3UxqEQJaFayv5E\nDoRidWm2D8I6pIHzzrpem34LtXEoBDdGFbJK6NLe9W/lSIjHKUvIrvA/tBcOOV9NQjxOWUJE6C3L\nG1RCyH99bH+3GKDn3G4LDoRsFpKmycTU03NutwUHQv4r7FaG6C3L45AQZJAQZFRHiHNpm7J43q2A\nSeJ+fbyHplw+HAtxKW1TFs875ldPfNDprC7EmhYq2aEe8CvEtbRNWTxv5ROt1wDEf66khXJ2qA98\nCRnTUO1+meda2qYsnjd+7oI/gKXeL0paKGeH+lS48SUk1vYr4oMai+jJi+d1Sb1c++TZMLY4mJQW\nqtmhHhVufAk58M9PFQ5pLKJnXTzvZkAOJE75rBe7g5QWqtmhHhVufAmxw7W0TVk87zvxxbIrbNLL\nIKeFanaoS4Ubt0I0StvkxfNWDBF1RQVuADktVLNDXSrc+BXiWtomL543dp64WST8CkpaqGSHulS4\n8SsEKSQEGSQEGSQEGSQEGSQEGSQEGQYW8rr6vdU6nSdVHQws5Ez/16WWZWOE7TpPqjoYWIhCeW9Z\ntoXxsnz3AWQ1/hwKX2hWp3s6eKshJudCbAvjwVOJUNBpIcD4565feJN9GW9tiHl9VIMI9h28Q8Xh\n29I7oUdKDzkXYlsYD476nR8y3gKWkKsAec+oDTEfH5md6Wd2qjgszsvLmzfW4u7/GAY/QoqeetWl\nI+aPtoXxAHpE97klvlG1HpIm/aStDTFvfVFw4e0HwKXiELYl3ARPZIj8CNmq0RHzH7Yl8gCmCFns\nujNJd9ydamuICUcEYatrxeEPXXPYzv0ZIj9Cjgiu2ZL9wnifR3WSVysrnNHEojTEFH/muW81cak4\nzLn3lPVat2eIfAuxLYwHByNOrokoBrNJfLzsCLUhZlIhsJeGc8XhxHnyI7g9Q+RbiG1hvLMR26C4\nwSrxpbDwWtaTb6gNMeuvLP09YZFzxeGP99yUH8HtGSLfQtSF8XKjl4inM9uaYW9MUNRbpWpDzM+a\nh3ZYaXGuOOzznvIIbs8Q+RaCEBKCDBKCDBKCDBKCDBKCDBKCDBKCDO6F/Mo6z0RvgiPiziTUPiOf\nKsO6w4WQr3NEdmoLWd8ZIDvF2o7s16ZvKqe2YX17YPIg5D/Kl+2LtW6dPlbc/NKKHV5sPVM9VfZ6\n98DkQcj1J/7Jwqi/Cp9o3Rq/DEoyElPFo9wOf7Wop8peTaT0WrCNAyEy7DNkTnALR2J/C2WvHaYq\n/4ERZrDIpxZ1WEmk9FqwjS8hw2oNduTPabVLoHRdfQsUxw+8BZApnyp7UHtg6rVgG19CZvg6j0o9\nlLPDLbcGxRfbnSp7Ww9MvRZs413IixPFD49hky0jH8i3O1X3th6Yei3YxruQHrVMpoaTi04Itdjv\nIavkU2UY1B6Yui3YxrsQdJAQZJAQZJAQZJAQZJAQZJAQZJAQZPAj5KDQK7mTj6fnU234EXLMLzg0\noJbreMRBcbPsMVBbXxb4mHwCe59jt3mhqo0fIQyNt6xzPgXiduxMUFtf7msGkD10ALvRWtUml7Od\n6h3ewvNfL3IvZBMrMmTN/dTWl0sGiQPngixqVZu1nO1q5JbStZ5PD/kTUrDvsB3nZ7GinFsB52yt\nL4e9Lo5cq2NRqtrkcrblIwEuBHm8FyZ/Qto7FLUFPSbtIuxaX0azROr7GFtVm7WcbWlMVuZjUqtl\nj4aH/Anpd9endvy3IVvda8XjttaXUsNYmD8W1Ko2aznbjYSguMQR0sN4MjzkT8if77EfyTKxz/Tk\nV2ytL3ez2s7cyExQqtrkcrYsM9yK2iLdyZPhIe9CNnZg25jNttaXCwfDzb2dFoBa1SaXs4VsLJ6S\nYP1M92R4yLuQGaPFTYHpoq31ZaJgqt19HYBa1SaXs30Q3njcdeudPBke8i4EHTVfyAX2d6I5JAQJ\nf1b+/Vqpd3US4nGO/iVFotYLlbl6qunTTx+M9PScqkVNF6JQubesvuy1FOjhqVQPvoQUpGdm9m3t\n4alUD76EMOgzRBdICDJICDKqI0RZPq+zj8knZAGAJfSYO6dWNUiIunxeiX8WWL6pC5AZcEu6Qe57\nad05N8P0FCREXT7vCFtTL6OlXCSirpYn7xybYXoO/oQkhE2z5//U5fNSBoivg2fFs5f+Kp3LfS+t\nO6dmmACl9Z/vE7lsUNdkgG8frt/cba8a/oQ85JAY+gxRl88bw87biz/sR6w/XbnvpdL+0qkZZobf\ncRj86M1sf/P+rsfNuwLcFbbzJ8QJs7p8Xkfx550eeAssIdYm/nLfS6X9pVMzzOXDAbrthyOdIZ6J\nDHPLfwWQENvyeUW+l8QDoQhO1yllN8hBobxzboYJo5ZBcXAJvPscBIv3Kyxyz38GCbEtn5fWhJ2J\nQj7uLo3LQaG8c26GCe3+BwfEz5Thq+G+14oP3bXPHf8R0gTc9UBe5vaFKMvnLe3HzoRimMranHRT\n+15ad87NMK+FlMKiKQBtTsPRTnU7fOGe/wogIeggIcggIcggIcggIcggIcggIcggIcjgUEhHH/EX\nv3ql1hTK1uzSO4uyucCfkCJfqXxQTqFszS7l8jW16eXOiFLPx1Gu8CckrZG0k1MotQemXL6mNr3c\n3SABPB9HucKfkKW9pZ2cQqk9MOXyNaXp5d6+yXN0iKNcMbaQ2XVaOLP8GRZfjJNTKLXZpboom7Xp\n5YH++V126BBHuWJsISN8nJpeDk7a3F4qfpJTKLXZpboom9T0Mm1QQSELrjweR7libCEv+7sMFfic\nZzs5hVKaXSrla3LTy55SvJvr+TjKFe6E7LtL2skplNLsUilfszW9XJwMOsRRrnAnZInUokFJoeRm\nl7byNbXpZRJ7ZXg8jnKFOyHYISHIICHIICHIICHIICHIICHIICHIICHI4FWIQ2wIhS80q9M9XRn1\nanDIqRDH2BDGP3f9wpstlFG5ek0pYtMXToU4xoaWkKsAec8oo3L1mrzTGT6EHEw/7MA5x9jQ0npI\nGlv/yzoqV68pRWz6rZgnwYWQi4IT4Q6xIcCZpDvuTgWQR+XqNWWn24p5ElwIuSSM+tSBIw6xIaNw\nRhMLtJc7KVqr15SdbivmSXAiZKnjuGNsaDYVAWRHWORRuXpN3jF0WjFPgk8hjrEhRC68lvXkG8qo\nXL0m7xg6rZgnYUAhXzboFCvTyLoagosQx9gQ9sYERb1VqozK1WtKERvotmKehAGFvCJ07ykTaZJG\nXIQgxoBClgtZymFZb1mIISHIML6Qy39PeVN7bXuUGF/Is+y3vVFa90GJ8YUsEY4cEt7xxpxuC+ML\neVu4Qp8h+kNCkEFCkHF7QpSkUO2DqR40OyRu3unnwRmXAd9C5KRQ7YOpHlzxLxZvHj5H/9iQayFK\nUqj2wVQPUu9nt0dv1T825ELIReHRaRocV5JCtQ+mejB/nLi54XNF/9iQCyE5vneEavCOkhSqfTDV\ng8RV4man+Game2zIhZBy7seSQrUPpq0hZnSGeDB/KOgfG/IsREkK1T6Y6oHZr1A8eORDL8SGPAtR\nkkK1D6btIHJOds6MmEIvxIZcC5GTQrUPpu0gvVtQs+cugxdiQ66FYISEIIOEIMPYQkb49OwZJVz2\nxoRuF2MLmUU3ZBIAAA1TSURBVF0vNrapcNEbE7pdjC2EQW9ZXoGEIIOEIKNKQtRl2KyNL9VTNaby\nIlwKUSrZ5MaXyqkaU+2PC43U7c/dneBSiFzJpjS+lE/VdOp8vS3mVY08N9dy4VKIXMmmNL6UT9V0\n6pPJAEdZaYKutWwyfAhZH9jc1gaz5Sa5kk1pfCmfqumUyOm4taBzLZsMH0LeFQYn20iTK9nkxpdy\nYZstnYKrE+KsH/N61rLJ8CHkE+GE3YhcyaY0vpRP1XQKvu6+wfqNu661bDI8CpEr2ZTGl0phm5JO\n7bj/cl6e9XWhZy2bDI9CrJVsauNLpbBNSacS2DvZQ9KVetayyfAoBDUkBBkkBBk8CDmTs4qE6E3Z\nQsZK3U326j2h28X4Qn58eu7cJCFD7wndLsYXwqDPEN0hIcggIciotpARJpPgYzId8N0HkNX489lR\nANeClwNsbfBhLLt99BR3TrdsSIhKsbRs3lOJUNBpITAhc/3mA/R99XA9cfh02BV9gkQSonI4gm2P\n+p0fMt7ChBQ1HTcdMoOyc4XrAEmv6RQk8iHkY+HTw+Uj/szft36z2yO6zy1gQpYPSxkLU5IBwo5B\nRuN8nYJEPoS84twF04WR4qfEP6Rrp0h3nx1VGpXx6dCCeuLPvcsmGGitGtUhSORDSN7bqRUg3iXm\nK3bp51Gd5gITsr4vpD66/DHxOGlxessS0ClI5ENIZSj0ZS0wD0acXBNRLApp03UPHO58X6o4NnNS\nL/bC0CdIJCEKB9mieWcjtkFxg1WiEN+HADL92jMFH90dYwa9gkQSovCO+EPOjV4iHs1sa4bZgvgG\nliOw6mjYJ2wB3YJEEoIMEoIMEoIMEoIMEoIMEoIMEoIMwwh5UG2DlSD85u3ZVAOjCPG3/6bwlLdn\nUw2MIsSOCt+y0C6Zx+BQSLlL5onsjCiFt6VX2o8enac2HAopd8k8gN0NEgCK8/Ly5o21eHKaZWBs\nIVOGu7a+XFneknkAe/smz5EOtiXcBC+UGRpbSB1/19aXseUumXegf36XHWz4h645bKd7maGxhbQc\npnF7eUvmpQ0qKJRqDXPulf+ppneZIX9Cyl0yr6f0YZ4LMHGefLneZYb8CSl3yTyRxcni5sd7bsqn\nepcZ8iek3CXzRJJWiJs+7ymnepcZ8icEOSQEGSQEGSQEGSQEGSQEGSQEGSQEGdwKcUyplP6X5uAT\nAD+1edfjkywT4wr5PTMz8+4yhTilVEr/y1NBpXCg+WZbH8ySOJZS5Q5p1HWlDlMHAwv5j/Q14b1l\nXeSYUqn9L9c8CBtbHgWlfO3W32JCzeLrpvtK856BeszdwEI2CNNTUsL7l3WRY0ql9r+c/PzSxuxv\nVuTytcKvpiaIZ6k95Zt1iKsMK+Rr4ZD1M2RIbY0F2h50TKnU/pcPt+sb+LP10Fq+9iJbVfLp0e3C\npMIEHeIq4wvpEZzsyt8dUirbMm13TDY/NZENKeVrD20TN20Ts9cES3+P4vm4yvhCnm6lcY1jSqX2\nv/yxrhkO1smxla+VBOaKvgIuweVw6dzzcRWnQpxSKqX/5eo/iCddXreVr6W1BdYb/suScS9JF3o+\nruJUiFNKpfS/nCR+lsPqRsVq+dqiUeyytY2avFgs3cHzcRWnQvBCQpBBQpBBQpBBQpBBQpBBQpBB\nQpBhbCH17x48+O5wb0+oShhbyD31W7SoG6JxjUNeWOLfh+1GsDVa9z8R3vSFAj2n64SxhTC03rIc\n88IjASzHOlGnKcDqqM35p7pNVDNC/RfQ41OIY16Y0i/YAjCo/0D4rdkF8fxoIyUj9MICepwI2Wzf\nzu+gxTEvTJ7fMAfSO0/4J4x5k53n+5nljNALC+jxIWSzY8fLXY55YacdXb6D+O3dtkG0lFdltLVl\nhLovoMeHkP3CG3avkB2lDnlhce1rgzfu7HEz4IrZJH2czxutZoT6L6DHi5Ct9lc45oXpbWDqkrj9\n30cCNGLZ+vnGZ5WM0AsL6HEpxDEvfD8J/tWuH6wYDDDpgTN5X7dcpmSE3lhAj0shjnnh6IXwVa1j\nMHae+NKZ2LjBAFYVbc0IvbGAHpdCMENCkEFCkGF8IQMapKRME7726pSqgPGFxEi/Cy716pSqgPGF\nFF3IydlGb1lexEkIgz5DvAkJQQYJQUaVhGg3xARLXG0vtcPkXIh2Q0yA9XHSGnrWEsNrzzcLmSVK\namddq8qzcC5EuyEm3Iw+3GoPKCWGfV67/KOfGfY8nJeXV+zhyXMiZJfQW6OMKvnF69oNMeG9UfCn\nj0ApMbzzi5xXRWNPxwY3/4bdSsvmVQ0NIWm+wRqFhqHNjms3xMyLugijZrMBqcRwsnjNf6Gg24ar\ns9g6VbRsXhXREFIm2g0xpWUPh7NTVmK49JGTlxfdI92aE+jxSkMDCnlfWJuamvpaZYRoN8S81OwG\nwIcsB5FKDKPPAPx8N2zIBDj0oHQ3WjavSkyU/5Thi4ov1W6IOWG2eLydvTtJJYYNPjD/HP8ODHkq\n/9duUktfWjavahSsZ4vbLq7MK0SzIebpUPHfWnBCKJRLDDe1qhe7zgKnuobEyv9Wo2XzboPKfYbg\ng4Qgg4Qgg4Qgg4Qgg4Qgg4Qgg4Qgg18hSLtgciukoi6YDtmUjnArpIIumI7ZlI5wK6SCLpiO2ZSO\nNW08CPnJ3zWe6l5hF0y7bAp0rGnjQUi6EO+S3i6qqAumczalV00bD0KOCetdbq6oC6ZDNsXQq6aN\nVyEVdMF0zKYYetW08Sqkgi6YjtkUQ6+aNl6FoIWEIIOEIMPwQvJycvaQEAQoQvzZL4BveHkyVcDw\nQgLjU2YKa7w8mSpgeCF1p9BnCApICDJICDIqFuIUGVqr2H5lv7pH69tm0QF+hThFhnIV2/rOANkp\nYaAsmJcWU/8jgFO9w1voJIlfIY6RoVLFNn2sePxLKyUxvNH0u23hcDVyS+naUH2CQz6EzE514ZJj\nZKhUscUvg5KMxFQlMVwxHgpM5uUjAS4EMSGl9Z/vE7lsUNdkgG8frt/8I/fPmwch3woaDHKMDOUq\nNksoG/2EjUiJ4aNfQmYkLI3Jynwslg1m+B2HwY/ezPY37+963LwrwP2vGh6EwMfbXV4gOy45VrPJ\nVWyZtUugdF199o27tGBe5M/wxaNwIyEoLnEEu2z5cIBu++FIZ4hndwhz/7y5EKKFY2SoVLFJnyjZ\nbIk864J5voXwlxWQZYZbUVvYdaOWQXFwCbz7HARfAigscv+8uRXiGBkqVWwvsi7jw1hKZV0wr17q\nkbtuQMjG4ikJ0rtTu//BgQcAhq+G+14rPnTXPvfPm1shDpGhWsXWo5bJ1HAy+x/fumDeB2Ed0sRt\neONx19nptZBSWCQ+YJvTcLRT3Q6VqGKsMtwKwQoJQQYJQYZhhWwUnpvLqP0Xb8+kahhWyArlN8Du\n3p5J1TCsEDifI3HHZG9PpGoYV4gMfYYgg4Qgg4QgQ1uIdu9L66iXul/K8ClEu/elPGprdrkzolSv\naarwKUS796U8qja73N0gQadJ2sGFkPOdRzqUT730tmbvSzlDVJpd7u2bPIedez4ltIcLIR8LES3s\nafuUZu9La4aoNLs80D+/i9RAzvMpoT1cCNkl7HAc1O592V79w5KcQHPaoIJCqfJQh5TQHi6FaPe+\nlEflZpc9pe9dckGPlNAeLoVo976UR23NLhcnSwOeTwnt4VKIZu9LZdTW7DJpBdvqkBLaw6UQzJAQ\nZJAQZBheiH/D2Ng2wgZvT6PSGF5Ih449e8aSEFzQWxYySAgySAgytIU4hFSzo8RfAYOXA2xtUAT7\nnwhv+kKBzpOU4VeIY0jFhMz1mw/Q91VYHbU5/1Q39mfXQxp1XanvVDkW4hhSiUKKmo6bDplB2b81\nuyAOHG0E5u4rzXsG6jhPBr9CHEvaRCHLh6WMhSnJMOZNNpDvZ07tKV+qWwdM4ETIx0Jdly6YMxxL\n2mZHlUZlfDq0oN4PEC2lJBlt4enR7cJWsWPdOmACJ0LOd580zZmtjiHV7Kj1fSH10eWPgdkkfZzP\nGw1tE7PXBHt8WTYnuBCihWNIBbPbdN0DhzvflwrQiLWMPd/4rDngElwOlxr86dUBEzgW4hRSzfZ9\nSPwXl1978ec/6YEzeV+3XAaW0C9Lxr0k3apXB0zgWIhTSDVb+AogR2AfGQUTGzcYwP4VsLZRkxet\na9/q1QETOBaCFRKCDBKCDBKCDBKCDBKCDBKCDBKCDBKCDBKCDBKCDBKCDBKCDBKCDBKCDBKCDBKC\nDBKCDBKCDBKCDBKCDBKCDBKCDBKCDBKCDBKCDBKCDBKCDBKCDBKCDBKCDBKCDBKCDBKCDBKCDBKC\nDBKCDBKCDBKCDBKCDBKCDBKCDBKCDBKCDBKCDBKCDBKCDBKCDBKCDBKCDBKCDBKCDBKCDBKCDBKC\nDBKCDBKCDBKCDBKCDBKCDBKCDBKCDBKCDBKCDBKCDBKCDBKCDBKCDBKCDBKCDBKCDBKCDBKCDBKC\nDBKCDBKCDBKCDBKCDBKCDBKCDBKCDBKCDBKCDBKCDBKCDBKCDBKCDBKCDBKCDBKCDBKCDBKCDBKC\nDBKCjP8HRDQfhHWXTKYAAAAASUVORK5CYII=\n"
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "%%R -w 400 -h 800\n",
    "library(ape)\n",
    "tre <- read.tree(\"empirical_3/RAxML_bipartitions.empirical_3\")\n",
    "ltre <- ladderize(tre)\n",
    "plot(ltre, edge.width=2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Meaure phylo distance (GTRgamma distance)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[1] 0.01793173\n"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "%%R\n",
    "mean(cophenetic.phylo(ltre))"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "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": 0
}