{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Assembling simulated and empirical data\n",
    "\n",
    "This Jupyter notebook provides a completely reproducible record of all the assembly analyses for the ipyrad manuscript. In this notebook we assemble multiple data sets using five different software programs: ipyrad, pyrad, stacks, aftrRAD, and dDocent. We include executable code to download and install each program, to download each data set, and to run each data set through each program."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Jupyter notebook SSH Tunneling\n",
    "This notebook was executed on the Yale farnam cluster with access to 40 compute cores on a single node. We performed local I/O in the notebook using SSH Tunneling as described in the ipyrad Documentation (http://ipyrad.readthedocs.io/HPC_Tunnel.html).\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Organize the directory structure"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import shutil\n",
    "import glob\n",
    "import sys\n",
    "import os"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "## Set the scratch dir for data analysis\n",
    "WORK_DIR = \"/ysm-gpfs/scratch60/de243/manuscript-analysis\"\n",
    "if not os.path.exists(WORK_DIR):\n",
    "    os.makedirs(WORK_DIR)\n",
    "    \n",
    "## Set a local dir for new (locally) installed software\n",
    "SOFTWARE = os.path.realpath(\"./tmp-software\")\n",
    "if not os.path.exists(SOFTWARE):\n",
    "    os.makedirs(SOFTWARE)\n",
    "    \n",
    "## Make paths to data (we download data into these dirs later)\n",
    "EMPIRICAL_DATA_DIR = os.path.join(WORK_DIR, \"empirical_data\")\n",
    "if not os.path.exists(EMPIRICAL_DATA_DIR):\n",
    "    os.makedirs(EMPIRICAL_DATA_DIR)\n",
    "    \n",
    "SIMULATED_DATA_DIR = os.path.join(WORK_DIR, \"simulated_data\")\n",
    "if not os.path.exists(SIMULATED_DATA_DIR):\n",
    "    os.makedirs(SIMULATED_DATA_DIR)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "## create dirs for results from each software\n",
    "PACKAGES = [\"ipyrad\", \"pyrad\", \"stacks\", \"ddocent\"]\n",
    "rdirs = {}\n",
    "for pack in PACKAGES:\n",
    "    rdir = os.path.join(WORK_DIR, \"assembly-{}\".format(pack))\n",
    "    if not os.path.exists(rdir):\n",
    "        os.makedirs(rdir)\n",
    "    rdirs[pack] = rdir\n",
    "\n",
    "## create subdirectories for results\n",
    "for key, val in rdirs.iteritems():\n",
    "    ## make for sim and real data\n",
    "    real = os.path.join(val, \"REALDATA\")\n",
    "    sim = os.path.join(val, \"SIMDATA\")\n",
    "    for idir in [real, sim]:\n",
    "        if not os.path.exists(idir):\n",
    "            os.makedirs(idir)\n",
    "            \n",
    "    ## extra subdirectories for stacks\n",
    "    if key == \"stacks\":\n",
    "        extras = [os.path.join(real, \"gapped\"), \n",
    "                  os.path.join(real, \"ungapped\"), \n",
    "                  os.path.join(real, \"default\")]\n",
    "        for idir in extras:\n",
    "            if not os.path.exists(idir):\n",
    "                os.makedirs(idir)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "## create simulated data directories\n",
    "SIMDATA = [\"simno\", \"simlo\", \"simhi\", \"simla\"]\n",
    "sdirs = {}\n",
    "for simd in SIMDATA:\n",
    "    sdir = os.path.join(SIMULATED_DATA_DIR, simd)\n",
    "    if not os.path.exists(sdir):\n",
    "        os.makedirs(sdir)\n",
    "    sdirs[simd] = sdir"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Download the empirical RAD *Pedicularis* data set\n",
    "\n",
    "This is a RAD data set for 13 taxa from Eaton and Ree (2013) [open access link](https://www.ncbi.nlm.nih.gov/pmc/articles/PMC3739883/pdf/syt032.pdf). Here we grab the demultiplexed fastq data from a public Dropbox link, but the same data is also hosted on the NCBI SRA as SRP021469.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash -s $EMPIRICAL_DATA_DIR --err stderr --out stdout\n",
    "\n",
    "## the location of the data \n",
    "dataloc=\"https://dl.dropboxusercontent.com/u/2538935/example_empirical_rad.tar.gz\"\n",
    "\n",
    "## grab data from the public link\n",
    "curl -LkO $dataloc \n",
    "tar -xvzf example_empirical_rad.tar.gz\n",
    "rm example_empirical_rad.tar.gz\n",
    "\n",
    "## move data to the EMPIRICAL_DATA_DIR\n",
    "mv ./example_empirical_rad/ $1/Pedicularis"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Download the empirical *Heliconius* reference genome¶\n",
    "Davey, John W., et al. \"Major improvements to the Heliconius melpomene \n",
    "genome assembly used to confirm 10 chromosome fusion events in 6 \n",
    "million years of butterfly evolution.\" G3: Genes| Genomes| Genetics 6.3\n",
    "(2016): 695-708."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash -s $EMPIRICAL_DATA_DIR --err stderr --out stdout\n",
    "\n",
    "## location of the data\n",
    "dataloc=\"http://butterflygenome.org/sites/default/files/Hmel2-0_Release_20151013.tgz\"\n",
    "\n",
    "## grab data from the public link \n",
    "curl -LkO $dataloc\n",
    "tar -zxvf Hmel2-0_Release_20151013.tgz\n",
    "rm Hmel2-0_Release_20151013.tgz\n",
    "\n",
    "## move data to EMPIRICAL_DATA_DIR\n",
    "mv ./Hmel2/ $1/Heliconius"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Install all the software "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### First install a new isolated miniconda directory\n",
    "\n",
    "We could create a separate environment in an existing conda installation, but this is simpler since it works for users whether they have conda installed or not, and we can easily remove all the software at the end when we are done by simply deleting the 'tmp-software' directory. So all of the work we do in this notebook will be done in a completely isolated software environment. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "%%bash -s $SOFTWARE --err stderr --out stdout\n",
    "\n",
    "## Fetch the latest miniconda \n",
    "miniconda=\"https://repo.continuum.io/miniconda/Miniconda-latest-Linux-x86_64.sh\"\n",
    "curl -LkO $miniconda\n",
    "\n",
    "## Install the new miniconda silently into the 'software' directory.\n",
    "## -b = batch mode; -f = force overwrite; -p = install dir\n",
    "bash Miniconda-latest-Linux-x86_64.sh -b -f -p $1\n",
    "\n",
    "## update conda, if this fails your other conda installation may be \n",
    "## getting in the way. Try clearing your ~/.condarc file.\n",
    "$1/bin/conda update -y conda\n",
    "\n",
    "## remove the install file\n",
    "rm Miniconda-latest-Linux-x86_64.sh "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "conda 4.3.4\n"
     ]
    }
   ],
   "source": [
    "%%bash -s $SOFTWARE \n",
    "## check installation\n",
    "$1/bin/conda --version"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Install ipyrad (Eaton & Overcast 2017)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash -s $SOFTWARE --err stderr --out stdout\n",
    "## installs the latest release (silently)\n",
    "$1/bin/conda install -y -c ipyrad ipyrad"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "ipyrad 0.5.15\n"
     ]
    }
   ],
   "source": [
    "%%bash -s $SOFTWARE \n",
    "## check installation\n",
    "$1/bin/ipyrad -v"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Install stacks (Catchen)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash -s $SOFTWARE --err stderr --out stdout\n",
    "## installs the latest release (silently)\n",
    "$1/bin/conda install -y -c bioconda stacks"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "ustacks 1.44\n",
      "\n"
     ]
    }
   ],
   "source": [
    "%%bash -s $SOFTWARE \n",
    "## check installation\n",
    "$1/bin/ustacks -v"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Install ddocent (Puritz)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash -s $SOFTWARE --err stderr --out stdout\n",
    "## install all of the ddocent reqs\n",
    "$1/bin/conda install -y -c conda-forge unzip  \n",
    "$1/bin/conda install -y -c bioconda ddocent=2.2.4 "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "dDocent 2.24 \n"
     ]
    }
   ],
   "source": [
    "%%bash -s $SOFTWARE --err stderr\n",
    "## when running dDocent you have to always set the tmp-software\n",
    "## directory to the top of your PATH\n",
    "export PATH=\"$1/bin:$PATH\"\n",
    "\n",
    "## check installation (no -v argument, print start of stdout)\n",
    "$1/bin/dDocent | head -n 1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Install pyrad (Eaton 2014)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 49,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash -s $SOFTWARE --err stderr --out stdout\n",
    "\n",
    "## Should be unnecessary because numpy and scipy already installed\n",
    "$1/bin/conda install numpy scipy\n",
    "\n",
    "## get muscle binary\n",
    "muscle=\"http://www.drive5.com/muscle/downloads3.8.31/muscle3.8.31_i86linux64.tar.gz\"\n",
    "curl -LkO $muscle \n",
    "tar -xvzf muscle*.tar.gz \n",
    "mv muscle3.8.31_i86linux64 $1/bin/muscle\n",
    "rm muscle3.8.31_i86linux64.tar.gz\n",
    "\n",
    "## Download and install vsearch\n",
    "vsearch=\"https://github.com/torognes/vsearch/releases/download/v2.0.3/vsearch-2.0.3-linux-x86_64.tar.gz\"\n",
    "curl -LkO $vsearch\n",
    "tar xzf vsearch-2.0.3-linux-x86_64.tar.gz\n",
    "mv vsearch-2.0.3-linux-x86_64/bin/vsearch $1/bin/vsearch\n",
    "rm -r vsearch-2.0.3-linux-x86_64/ vsearch-2.0.3-linux-x86_64.tar.gz\n",
    "\n",
    "## Fetch pyrad source from git repository \n",
    "if [ ! -d ./pyrad-git ]; then\n",
    "  git clone https://github.com/dereneaton/pyrad.git pyrad-git\n",
    "fi;\n",
    "\n",
    "## and install to conda using pip\n",
    "cd ./pyrad-git\n",
    "$1/bin/pip install -e . \n",
    "cd .."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 51,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "pyRAD 3.0.66\n"
     ]
    }
   ],
   "source": [
    "%%bash -s $SOFTWARE \n",
    "## check installation\n",
    "$1/bin/pyrad --version"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Simulate data\n",
    "We will use simrrls to generate some simulated RAD-seq data for testing. This is a program that was written by Deren Eaton and is available on github: github.com/dereneaton/simrrls.git. simrrls requires the python egglib module, which is a pain to install in full, but we only need the simulation aspects of it, which are fairly easy to install. See below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 57,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash -s $SOFTWARE --err stderr --out stdout\n",
    "\n",
    "## install gsl\n",
    "$1/bin/conda install gsl -y \n",
    "\n",
    "## fetch egglib cpp and mv into tmp-software/pkgs\n",
    "eggcpp=\"http://mycor.nancy.inra.fr/egglib/releases/2.1.11/egglib-cpp-2.1.11.tar.gz\"\n",
    "curl -LkO $eggcpp\n",
    "tar -xzvf egglib-cpp-2.1.11.tar.gz \n",
    "mv egglib-cpp-2.1.11/ $1/pkgs/\n",
    "\n",
    "## install egglib-cpp\n",
    "cd $1/pkgx/egglib-cpp-2.1.11/\n",
    "sh ./configure --prefix=$1\n",
    "make\n",
    "make install \n",
    "cd -\n",
    "\n",
    "## fetch egglib py module and mv into tmp-software/pkgs\n",
    "eggpy=\"http://mycor.nancy.inra.fr/egglib/releases/2.1.11/egglib-py-2.1.11.tar.gz\"\n",
    "curl -LkO $eggpy \n",
    "tar -xvzf egglib-py-2.1.11.tar.gz\n",
    "mv egglib-py-2.1.11/ $1/pkgs/\n",
    "\n",
    "## install egglib-py\n",
    "cd $1/pkgs/egglib-py-2.1.11/\n",
    "$1/bin/python setup.py build --prefix=$1\n",
    "$1/bin/python setup.py install --prefix=$1\n",
    "cd -\n",
    "\n",
    "## clone simrrls into tmp-software/pkgs/\n",
    "if [ ! -d $1/pkgs/simrrls ]; then\n",
    "  git clone https://github.com/dereneaton/simrrls.git $1/pkgs/simrrls\n",
    "fi\n",
    "\n",
    "## install with tmp-software/bin/pip\n",
    "cd $1/pkgs/simrrls\n",
    "$1/bin/pip install -e .\n",
    "cd -"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 68,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "egglib Version number: 2.1.11\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "simrrls 0.0.13\n"
     ]
    }
   ],
   "source": [
    "%%bash -s $SOFTWARE \n",
    "## check installations\n",
    "echo 'egglib' $($1/bin/egglib | grep Version)\n",
    "$1/bin/simrrls --version"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "### Simulating different RAD-datasets\n",
    "\n",
    "Both pyRAD and stacks have undergone a lot of work since the original pyrad analysis. Because improvements have been made we want to test performance of all the current pipelines and be able to compare current to past performance. We'll follow the original pyRAD manuscript analysis (Eaton 2013) by simulating modest sized datasets with variable amounts of indels. We'll also simulate one much larger dataset. Also, because stacks has since included an option for handling gapped analysis we'll test both gapped and ungapped assembly.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Tuning simrrls indel parameter\n",
    "\n",
    "The -I parameter for simrrls has changed since the initial pyrad manuscript, so the we had to explore new values for this parameter that will approximate the number of indels we are after. I figured out a way to run simrrls and pipe the output to muscle to get a quick idea of the indel variation for different params. This gives a good idea of how many indel bearing seqs are generated.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# simrrls -n 1 -L 1 -I 1 -r1 $RANDOM 2>&1 | \\\n",
    "#                             grep 0 -A 1 | \\\n",
    "#                             tr '@' '>' | \\\n",
    "#                             muscle | grep T | head -n 60"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# for i in {1..50}; \n",
    "#   do simrrls -n 1 -L 1 -I .05 -r1 $RANDOM 2>&1 | \\\n",
    "#                                   grep 0 -A 1 | \\\n",
    "#                                   tr '@' '>' | \\\n",
    "#                                   muscle |  \\\n",
    "#                                   grep T | \\\n",
    "#                                   head -n 40 >> rpt.txt;\n",
    "# done\n",
    "# grep \"-\" rpt.txt | wc -l"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "From experimentation:\n",
    "\n",
    "-I value -- %loci w/ indels\n",
    "\n",
    "    0.02 -- ~10%\n",
    "    0.05 -- ~15%\n",
    "    0.10 -- ~25%\n",
    "\n",
    "The simulated data will live in these directories:\n",
    "\n",
    "    SIM_NO_DIR = WORK_DIR/simulated_data/simno\n",
    "    SIM_LO_DIR = WORK_DIR/simulated_data/simlo\n",
    "    SIM_HI_DIR = WORK_DIR/simulated_data/simhi\n",
    "    SIM_LARGE_DIR = WORK_DIR/simulated_data/simlarge\n",
    "\n",
    "Timing:\n",
    "\n",
    "    10K loci -- ~8MB -- ~ 2 minutes\n",
    "    100K loci -- ~80MB -- ~ 20 minutes\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Call *simrrls* to simulate data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 70,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "%%bash -s $SOFTWARE $SIMULATED_DATA_DIR\n",
    "\n",
    "## call simrrls (No INDELS)\n",
    "$1/bin/simrrls -o $2/simno/simno -ds 2 -L 10000 -I 0 \n",
    "\n",
    "## call simrrls (Low INDELS)\n",
    "$1/bin/simrrls -o $2/simlo/simlo -ds 2 -L 10000 -I 0.02\n",
    "\n",
    "## call simrrls (High INDELS) \n",
    "$1/bin/simrrls -o $2/simhi/simhi -ds 2 -L 10000 -I 0.05\n",
    "\n",
    "## call simrls on Large data set (this takes a few hours, sorry!)\n",
    "## (30x12=360 Individuals at 100K loci) = gzipped 2.2GB\n",
    "$1/bin/simrrls -o $2/simla/Big_i360_L100K -ds 0 -L 100000 -n 30"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "### Demultiplex all four libraries\n",
    "\n",
    "We will start each analysis from the demultiplexed data files, since this is commonly what's available when data are downloaded from NCBI. We use ipyrad to demultiplex the data.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "  Assembly: simno\n",
      "  [####################] 100%  chunking large files  | 0:00:00 | s1 | \n",
      "  [####################] 100%  sorting reads         | 0:00:28 | s1 | \n",
      "  [####################] 100%  writing/compressing   | 0:00:04 | s1 | \n",
      "\n",
      "  Assembly: simlo\n",
      "  [####################] 100%  chunking large files  | 0:00:00 | s1 | \n",
      "  [####################] 100%  sorting reads         | 0:00:28 | s1 | \n",
      "  [####################] 100%  writing/compressing   | 0:00:04 | s1 | \n",
      "\n",
      "  Assembly: simhi\n",
      "  [####################] 100%  chunking large files  | 0:00:00 | s1 | \n",
      "  [####################] 100%  sorting reads         | 0:00:30 | s1 | \n",
      "  [####################] 100%  writing/compressing   | 0:00:04 | s1 | \n",
      "\n",
      "  Assembly: simla\n",
      "  [                    ]   0%  chunking large files  | 0:03:38 | s1 | "
     ]
    }
   ],
   "source": [
    "import ipyrad as ip\n",
    "\n",
    "for name in [\"simno\", \"simlo\", \"simhi\", \"simla\"]:\n",
    "    ## demultiplex the library\n",
    "    data = ip.Assembly(name, quiet=True)\n",
    "    pdir = os.path.join(SIMULATED_DATA_DIR, name)\n",
    "    data.set_params(\"project_dir\", pdir)\n",
    "    data.set_params(\"raw_fastq_path\", os.path.join(pdir, \"*.gz\"))\n",
    "    data.set_params(\"barcodes_path\", os.path.join(pdir, \"*barcodes.txt\"))\n",
    "    data.run(\"1\")\n",
    "    print \"demultiplexed fastq data written to: {}\".format(data.dirs.fastqs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "print \"demultiplexed fastq data written to: {}\".format(data.dirs.fastqs)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Assemble simulated data sets"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### ipyrad : simulated data assembly\n",
    "Results  \n",
    "Data set \tCores \tLoci \tTime  \n",
    "simno \t4 \t10000 \t13:58  \n",
    "simlo \t4 \t10000 \t15:17  \n",
    "simhi \t4 \t10000 \t13:58  \n",
    "simla \t4 \t100000 \t13:38  \n",
    "simla \t80 \t100000 \t...  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "## A function to run ipyrad on a given data set\n",
    "\n",
    "def run_ipyrad(name):\n",
    "    ## create Assembly\n",
    "    data = ip.Assembly(name)\n",
    "\n",
    "    ## set data paths and params\n",
    "    data.set_params(\"project_dir\", \n",
    "        os.path.join(SIMULATED_DATA_DIR, name))\n",
    "    data.set_params(\"sorted_fastq_path\", \n",
    "        os.path.join(SIMULATED_DATA_DIR, name, name+\"_fastqs\", \"*.gz\"))\n",
    "    data.set_params('max_low_qual_bases', 4)\n",
    "    data.set_params('filter_min_trim_len', 69)\n",
    "    data.set_params('max_Ns_consens', (99,99))\n",
    "    data.set_params('max_Hs_consens', (99,99))\n",
    "    data.set_params('max_SNPs_locus', (100, 100))\n",
    "    data.set_params('min_samples_locus', 2)\n",
    "    data.set_params('max_Indels_locus', (99,99))\n",
    "    data.set_params('max_shared_Hs_locus', 99)\n",
    "\n",
    "    ## run ipyrad steps 1-7\n",
    "    data.run(\"1234567\", show_cluster=True)\n",
    "    return data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "%%timeit -n1 -r1 \n",
    "## this records time to run the code in this cell once\n",
    "run_ipyrad(\"simno\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "%%timeit -n1 -r1 \n",
    "## this records time to run the code in this cell once\n",
    "run_ipyrad(\"simlo\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "%%timeit -n1 -r1 \n",
    "## this records time to run the code in this cell once\n",
    "run_ipyrad(\"simhi\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "%%timeit -n1 -r1 \n",
    "## this records time to run the code in this cell once\n",
    "run_ipyrad(\"simla\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### stacks: simulated data assembly"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "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.15"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}