combo_biome.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# One big biome table"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Since I am hoping to make comparisons of wood and leaf endophyte environmental patterns, I need to combine these datasets early in the biomformatics pipeline, to make them as comparable as possible. We'll try to stick to the [usearch (uparse)](http://drive5.com/usearch/) pipeline for the process, as much as possible."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Table of contents"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[Work environment](#environment)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[Merge paired ends](#mergepairs)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[Quality filtering](#qf)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[Convert form fastq to fasta format](#fastq2fasta)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[Convert form fastq to fasta format](#fastq2fasta)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[Demultiplex leaf reads](#demult)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[Remove primers](#removeprimers)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[Floating primers](#defloat)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[Chimera checking](#chimera)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[Combining fasta files from all 3 studies](#combine)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[OTU clustering](#otuclust)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[Customizing UNITE database](#UNITE)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[Assign taxonomy](#asstax)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[Make and tidy up biom table](#makebiom)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[Add metadata](#addmetadata)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='environment'></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Work environment"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Working directory, on my machine:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": []
    }
   ],
   "source": [
    "cd /home/daniel/Documents/Taiwan_data/combined/combo_biome"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We'll be using the [usearch (uparse)](http://drive5.com/usearch/) pipeline, version usearch v8.1.1861_i86linux32 on my personal machine, and the equivalent 64-bit version on our computing cluster [ACISS](http://aciss-computing.uoregon.edu/). These programs are abbreviated (soft linked) to \"usearch81\" and \"usearch\", respectively. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='mergepairs'></a> "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Merging paired-end reads"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's re-pair all readsets in the same manner, except for Roo's stromatal readset, which he aligned by hand. \n",
    "\n",
    "First, the leaf study reads include a split 6+6 bp barcode scheme for identifying reads, so these need to be clipped from one read and combined on the other. I wrote a python script for this, let's see if it works:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": []
    }
   ],
   "source": [
    "./BCunsplit4.py Roo_R2.fastq Roo_R1.fastq"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This outputs two files, \"rearranged_Roo_R2.fastq\" and \"rearranged_Roo_R2.fastq\"\n",
    "\n",
    "Next we trim a little to make sure we doing our alignments with high quality base calls. The sites for trimming are decided by looking at the raw reads [(see below)](#quality), and finding where quality begins to drop off. \n",
    "To trim, we'll use the [FASTX-toolkit](http://hannonlab.cshl.edu/fastx_toolkit/)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": []
    }
   ],
   "source": [
    "## wood\n",
    "fastx_trimmer -l 255 -i woodR1.fastq -o woodR1_trimmed.fastq\n",
    "fastx_trimmer -l 210 -i woodR2.fastq -o woodR2_trimmed.fastq\n",
    "\n",
    "## leaves. These lengths were decided by Roo\n",
    "fastx_trimmer -l 263 -i rearranged_Roo_R1.fastq -o Roo_R1_trimmed.fastq\n",
    "fastx_trimmer -l 170 -i rearranged_Roo_R2.fastq -o Roo_R2_trimmed.fastq"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Do the actual pairing:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": []
    }
   ],
   "source": [
    "## wood, let's pair both trimmed and untrimmed just to compare:\n",
    "usearch -fastq_mergepairs woodR1.fastq -reverse woodR2.fastq  -fastqout woodtrimmedmerged.fastq -notrunclabels\n",
    "usearch -fastq_mergepairs woodR1.fastq -reverse woodR2.fastq -fastqout woodmerged.fastq -notrunclabels"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The \"-notrunclabels\" tag above asks usearch to keep the entire label of the forward reads, which is necessary because the wood reads, which are more recently sequenced than the leaf reads, contain sample info in their identifier lines. The leaves do not require this, their sample info is still in the sequence itself, to be use to [demultiplex](#demult) them, later."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "** typical report from these:**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "usearch v8.1.1803_i86linux64, 74.2Gb RAM, 12 cores\n",
    "(C) Copyright 2013-15 Robert C. Edgar, all rights reserved.\n",
    "http://drive5.com/usearch\n",
    "\n",
    "Licensed to: ronh@molbio.uoregon.edu\n",
    "\n",
    "<br>03:02 925Mb  100.0% 95.6% merged<br>5567799  Pairs (5.6M)<br>5323274  Merged (5.3M, 95.61%)<br>2067236  Alignments with zero diffs (37.13%)<br>0  Fwd tails Q <= 2 trimmed (0.00%)<br>15  Rev tails Q <= 2 trimmed (0.00%)<br>244525  No alignment found (4.39%)<br>0  Alignment too short (< 16) (0.00%)<br>4512849  Staggered pairs (81.05%) merged & trimmed<br>239.63  Mean alignment length\n",
    "259.16  Mean merged read length<br>3.81  Mean fwd expected errors<br>6.27  Mean rev expected errors<br>1.57  Mean merged expected errors"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "These numbers look pretty good. Many of the erroneous reads will be taken out below, in a [quality filtering](#qf) step"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": []
    }
   ],
   "source": [
    "## leaves\n",
    "usearch -fastq_mergepairs Roo_R2.fastq -reverse Roo_R1.fastq -fastqout leafmerged.fastq\n",
    "usearch -fastq_mergepairs Roo_R2_trimmed.fastq -reverse Roo_R1_trimmed.fastq -fastqout leaftrimmedmerged.fastq"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "What did trimming do to our reads? Let's take a look. To plot these we'll use the [fastx wrapper](http://hannonlab.cshl.edu/fastx_toolkit/commandline.html#fastq_quality_boxplot_usage) for [gnuplot](http://www.gnuplot.info/), which I've tinkered with just a little to change up the crowded axes of the original setup. Fastx requires that we first compile the quality data from the fastq files:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": []
    }
   ],
   "source": [
    "## wood quality stats:\n",
    "fastx_quality_stats -i woodR1.fastq -o woodR1_fastxstats.txt\n",
    "fastx_quality_stats -i woodR2.fastq -o woodR2_fastxstats.txt\n",
    "fastx_quality_stats -i woodR1_trimmed.fastq -o woodR1_fastxstats.txt\n",
    "fastx_quality_stats -i woodR2_trimmed.fastq -o woodR2_fastxstats.txt\n",
    "fastx_quality_stats -i woodmerged.fastq -o woodmerged_fastxstats.txt\n",
    "fastx_quality_stats -i woodtrimmedmerged.faestq -o woodtrimmedmerged_fastxstats.txt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "## leaf quality stats:\n",
    "fastx_quality_stats -i Roo_R1.fastq -o Roo_R1_fastxstats.txt\n",
    "fastx_quality_stats -i Roo_R2.fastq -o Roo_R2_fastxstats.txt\n",
    "fastx_quality_stats -i Roo_R1_trimmed.fastq -o Roo_R1_trimmed_fastxstats.txt\n",
    "fastx_quality_stats -i Roo_R2_trimmed.fastq -o Roo_R2_trimmed_fastxstats.txt\n",
    "fastx_quality_stats -i leafmerged.fastq -o leafmerged_fastxstats.txt\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Make the graphics:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": []
    }
   ],
   "source": [
    "## wood graphics\n",
    "./dan_plot.sh -i woodR1_fastxstats.txt -o woodR1_quality.png\n",
    "./dan_plot.sh -i woodR2_fastxstats.txt -o woodR2_quality.png\n",
    "./dan_plot.sh -i  woodR1_trimmed_fastxstats.txt -o woodR1_trimmed_quality.png\n",
    "./dan_plot.sh -i  woodR2_trimmed_fastxstats.txt -o woodR2_trimmed_quality.png\n",
    "./dan_plot.sh -i woodmerged_fastxstats.txt -o woodmerged_quality.png\n",
    "./dan_plot.sh -i  woodtrimmedmerged_fastxstats.txt -o woodtrimmedmerged_quality.png\n",
    "\n",
    "## leaf graphics\n",
    "./dan_plot.sh -i  Roo_R1_trimmed_fastxstats.txt -o Roo_R1_trimmed_quality.png\n",
    "./dan_plot.sh -i  Roo_R2_trimmed_fastxstats.txt -o Roo_R2_trimmed_quality.png\n",
    "./dan_plot.sh -i leaftrimmedmerged_fastxstats.txt -o leaftrimmedmerged_quality.png"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='quality'></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The **untrimmed wood R1** reads look like this:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img src='woodR1_quality.png'>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Compare to **trimmed wood R1** reads:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img src='woodR1_trimmed_quality.png'>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The **untrimmed wood R2** reads look like this:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img src='woodR2_quality.png'>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Compare to **trimmed wood R2** reads:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img src='woodR2_trimmed_quality.png'>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And the **untrimmed, merged wood** file looks like this:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img src='woodmerged_quality.png'>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Obviously some problems here. So compare to the **merged, trimmed wood** reads:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img src='woodtrimmedmerged_quality.png'>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Looks much better, but still a large dip in quality around 15 bp. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Roo has already decided the trimming sites for his leaf data [(see above)](#mergepairs). Skip to the trimmed leaf R1 reads:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img src='Roo_R1_trimmed_quality.png'>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The trimmed leaf R2 reads:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img src='Roo_R2_trimmed_quality.png'>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Trimmed, merged leaf reads:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img src='leaftrimmedmerged_quality.png'>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In both leaves and wood, I think we've vastly improved the situation by merging. The nice thing about using the usearch merging algorithms is that they use a bayesian approach to calculating the q-scores of paired reads, so that agreements on base calls between R1 and R2 improve Q scores, and disagreements reduce them. The reduction by disagreement is proportional to the confidence of the two base calls at a site, so if the disagreement occurs at the end of a read, where quality is lower, the higher (more reliable) Q has a much greater influence on the final q-score of a base call. Check out the explanation [here](http://drive5.com/usearch/manual/exp_errs.html)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='qf'></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Quality filtering"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Continuing with the usearch pipeline, let's do some quality filtering. We'll use the [expected error approach](http://drive5.com/usearch/manual/exp_errs.html). We can set error cutoff of 1% of all bases in a read, meaning that a read of length 400 bp is thrown out if it likely contains 4 or more erroneous bases. I think this is permissable, given our OTU clustering will ultimately be done at 95% similarity."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": []
    }
   ],
   "source": [
    "## wood sequences\n",
    "usearch -fastq_filter woodtrimmedmerged.fastq -fastq_maxee_rate .01 -fastqout wood_filtered.fastq -notrunclabels"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "usearch v8.1.1803_i86linux64, 74.2Gb RAM, 12 cores\n",
    "(C) Copyright 2013-15 Robert C. Edgar, all rights reserved.\n",
    "http://drive5.com/usearch\n",
    "\n",
    "Licensed to: ronh@molbio.uoregon.edu\n",
    "\n",
    "01:04 857Mb  100.0% Filtering, 85.4% passed\n",
    "   5323274  FASTQ recs (5.3M)              \n",
    "   4548698  Converted (4.5M, 85.4%)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": []
    }
   ],
   "source": [
    "## leaves\n",
    "usearch -fastq_filter leaftrimmedmerged.fastq -fastq_maxee_rate 0.01 -fastqout leaf_filtered.fastq"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "usearch v8.1.1803_i86linux64, 529Gb RAM, 32 cores\n",
    "(C) Copyright 2013-15 Robert C. Edgar, all rights reserved.\n",
    "http://drive5.com/usearch\n",
    "\n",
    "Licensed to: ronh@molbio.uoregon.edu\n",
    "\n",
    "07:01 2.4Gb  100.0% Filtering, 90.7% passed\n",
    "  16701565  FASTQ recs (16.7M)             \n",
    "  15145323  Converted (15.1M, 90.7%)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can inspect these graphically as above, to see if avg read quality has improved. Not placing this here because the graphs look basically the same before and after filtering. But notice that we drop 15% of wood reads, and 9% of the leaf reads. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "<a id='fastq2fasta'></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Convert from fastq to fasta format"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now that we've merged paired ends and done some quality filtering to ensure that we've hopefully mostly eliminated sequencer error (hah! see [index bleed](#bleed) below), let's convert to fasta as required by most downstream steps. Using FASTX toolkit again:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": []
    }
   ],
   "source": [
    "fastq_to_fasta -n -i leaf_filtered.fastq -o  leaf.fasta\n",
    "fastq_to_fasta -n -i wood_filtered.fastq -o  wood.fasta"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The \"-n\" flag tells fastx to retain sequences with \"N\" basecalls. Otherwise, these are removed, by default. Since ~1/2 of our leaf reads contain an \"N\", we need these. A single N basecall is an acceptable loss of information, OTU clustering and taxonomic assignments should be to deal with this. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id=\"demult\"></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "## Demultiplex leaf reads"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Roo's leaf reads were prepped at an earlier date than the wood reads. At the time of their sequencing, the standard method for denoting sample identities was to look for the presence of 12 bp golay barcodes, the ones that we cut and pasted when we [merged paired end reads](#mergepairs). For probably the only time in this pipeline, we will use a [qiime](http://qiime.org/) script, [\"demultiplex_fasta.py\"](http://qiime.org/scripts/demultiplex_fasta.html) that was made to parse samples by these golay barcodes."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This script requires a mapping file that lists the barcodes and their accompanying sample info, plus a \"linkerprimersequence\". I use a map file supplied by Roo. I do not know what a \"linkerprimersequence\" is, I believe this is supplied by the illumina software. The script seems to prefer .tsv format to .csv, and looks like this:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "#SampleID\tBarcodeSequence\tLinkerPrimerSequence\r\n",
      "1Leaf\tACCCATATATCC\tGCTGCGTTCTTCATCGATGC\r\n",
      "2Leaf\tACCCATAAGACG\tGCTGCGTTCTTCATCGATGC\r\n",
      "3Leaf\tTCGCCAGAACCA\tGCTGCGTTCTTCATCGATGC\r\n",
      "4Leaf\tACCCATATCAAA\tGCTGCGTTCTTCATCGATGC\r\n",
      "5Leaf\tACCCATATAGTA\tGCTGCGTTCTTCATCGATGC\r\n",
      "6Leaf\tACCCATCTACAG\tGCTGCGTTCTTCATCGATGC\r\n",
      "7Leaf\tACCCATCATACC\tGCTGCGTTCTTCATCGATGC\r\n",
      "8Leaf\tACCCATCATTAT\tGCTGCGTTCTTCATCGATGC\r\n",
      "9Leaf\tACCCATCTATCT\tGCTGCGTTCTTCATCGATGC\r\n"
     ]
    }
   ],
   "source": [
    "head leaf_sample_map.tsv"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": []
    }
   ],
   "source": [
    "demultiplex_fasta.py -m leaf_sample_map.tsv -f leaf.fasta -o ./leaf_demult"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This produces a folder with a log file and a file called \"demultiplexed_seqs.fna\". I will rename this to \"leaf.fna\" and bring it into our working directory. We can remove primers from the sequences now. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='removeprimers'></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Remove primers"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "PCR Primers in our wood reads were not included in these sequences, because they are used as part of the sequencing primers. The leaf reads, however, still contain our forward and reverse PCR primers. To remove these, we just clip the appropriate number of BPs from each end. We'll use the FASTX toolkit again. We will trim these reads further [below](#trim3) when combining the three studies (stromata, leaf endophyte, and wood endophyte) into a single fasta file. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "How long are our primers?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "22\r\n",
      "20\r\n"
     ]
    }
   ],
   "source": [
    "##ITS1F\n",
    "expr length \"CTTGGTCATTTAGAGGAAGTAA\"\n",
    "## ITS2\n",
    "expr length \"GCTGCGTTCTTCATCGATGC\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": []
    }
   ],
   "source": [
    "fastx_trimmer -f 23 -i leaf.fna | fastx_trimmer -t 20 -o leaf_noprim.fna"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "On my personal machine this took about 20 minutes."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='defloat'></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Floating primers"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In both of our read sets, \"floating\" primer sequences appear. This happens in other studies, as indicated by [Bálint et al. (2014)](http://onlinelibrary.wiley.com/doi/10.1002/ece3.1107/abstract;jsessionid=FBCBBBE428CAA870889926051DBC9927.f04t02). As advised by these authors, I remove the sequences that contain these floating primers with a script."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": []
    }
   ],
   "source": [
    "## ITS1F\n",
    "grep CTTGGTCATTTAGAGGAAGTAA wood.fasta | wc -l"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "129"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": []
    }
   ],
   "source": [
    "## ITS2\n",
    "grep GCTGCGTTCTTCATCGATGC wood.fasta | wc -l"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "18"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Not many. This doesn't include the reverse complements that also occur. To remove these floating primers:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": []
    }
   ],
   "source": [
    "## wood\n",
    "./floatingprimers.py wood.fasta wood_defloat.fasta CTTGGTCATTTAGAGGAAGTAA GCTGCGTTCTTCATCGATGC\n",
    "## leaves\n",
    "./floatingprimers.py leaf_noprim.fna leaf_defloat.fna CTTGGTCATTTAGAGGAAGTAA GCTGCGTTCTTCATCGATGC"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Arguments for this script: floatingprimers.py input_fasta_file output_fasta_file forward_primer reverse_primer. \n",
    "\n",
    "Reverse compliments are checked automatically from the forward and reverse primer sequences that are given as arguments.\n",
    "\n",
    "When doing this, we lose 230 reads from the wood sequences (0.005%), and 312 from the leaves (0.002%). There were many more in the raw reads, but I think they were removed by the quality filtering."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='chimera'></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Chimera checking"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Time to look for chimeras. We'll use the [uchime](http://www.drive5.com/usearch/manual/cmd_uchime_ref.html) algorithm, another step in the uparse/usearch pipeline. This is actually just the first of two checks for chimeras, the other being part of the [otu clustering](#otus).\n",
    "\n",
    "Get the latest [UNITE](https://unite.ut.ee/repository.php) database for usearch:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "--2016-08-21 18:27:34--  https://unite.ut.ee/sh_files/uchime_reference_dataset_01.01.2016.zip\r\n",
      "Resolving unite.ut.ee (unite.ut.ee)... 2001:bb8:2002:500:ec4:7aff:fe0a:37b2, 193.40.5.164\r\n",
      "Connecting to unite.ut.ee (unite.ut.ee)|2001:bb8:2002:500:ec4:7aff:fe0a:37b2|:443... connected.\r\n",
      "HTTP request sent, awaiting response... 200 OK\r\n",
      "Length: 8233914 (7.9M) [application/zip]\r\n",
      "Saving to: ‘uchime_reference_dataset_01.01.2016.zip’\r\n",
      "\r\n",
      "\r",
      " 0% [                                       ] 0           --.-K/s              \r",
      " 1% [                                       ] 89,823       386KB/s             \r",
      " 2% [>                                      ] 220,895      384KB/s             \r",
      " 7% [=>                                     ] 597,727      766KB/s             \r",
      "14% [====>                                  ] 1,203,935   1.19MB/s             \r",
      "27% [=========>                             ] 2,268,895   1.86MB/s             \r",
      "43% [================>                      ] 3,612,383   2.53MB/s             \r",
      "61% [======================>                ] 5,054,175   3.08MB/s             \r",
      "78% [=============================>         ] 6,495,967   3.51MB/s             \r",
      "96% [====================================>  ] 7,954,143   3.86MB/s             \r",
      "100%[======================================>] 8,233,914   3.95MB/s   in 2.0s   \r\n",
      "\r\n",
      "2016-08-21 18:27:37 (3.95 MB/s) - ‘uchime_reference_dataset_01.01.2016.zip’ saved [8233914/8233914]\r\n",
      "\r\n"
     ]
    }
   ],
   "source": [
    "wget https://unite.ut.ee/sh_files/uchime_reference_dataset_01.01.2016.zip"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Unzip this, then use the ITS1 only database to check for chimeras:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": []
    }
   ],
   "source": [
    "## wood\n",
    "usearch -uchime_ref wood_defloat.fasta -db /home6/dthomas/combobiom/uchime_reference_dataset_01.01.2016/ITS1_ITS2_datasets/uchime_sh_refs_dynamic_develop_985_01.01.2016.ITS1.fasta -nonchimeras wood_notchim.fasta -strand plus -uchimeout woodchim_log.txt -notrunclabels"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "usearch v8.1.1861_i86linux64, 74.2Gb RAM, 12 cores\n",
    "(C) Copyright 2013-15 Robert C. Edgar, all rights reserved.\n",
    "http://drive5.com/usearch\n",
    "\n",
    "Licensed to: ronh@molbio.uoregon.edu\n",
    "\n",
    "00:19 1.7Gb  100.0% Reading wood_defloat.fasta\n",
    "00:19 1.7Gb    0.1% Reading /home6/dthomas/combobiom/uchime_reference_dataset_01.01.2016/ITS1_ITS2_datasets/uchime_sh_refs_dynamic_dev00:19 1.7Gb  100.0% Reading /home6/dthomas/combobiom/uchime_reference_dataset_01.01.2016/ITS1_ITS2_datasets/uchime_sh_refs_dynamic_develop_985_01.01.2016.ITS1.fasta\n",
    "00:19 1.7Gb    0.1% Masking                                                                                                           00:19 1.7Gb  100.0% Masking   \n",
    "00:19 1.7Gb  100.0% Word stats\n",
    "00:19 1.7Gb  100.0% Alloc rows\n",
    "00:20 1.7Gb  100.0% Build index\n",
    "15:13 14.7Gb  100.0% Found 17692/4548468 chimeras (0.4%), 3623750 not classified (79.7%)\n",
    "16:03 14.7Gb  100.0% Writing hits                                                       \n",
    "17:07 14.7Gb  100.0% Writing 4530776 non-chimeras and unclassifieds"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": []
    }
   ],
   "source": [
    "## leaf\n",
    "usearch -uchime_ref leaf_defloat.fna -db /home6/dthomas/combobiom/uchime_reference_dataset_01.01.2016/ITS1_ITS2_datasets/uchime_sh_refs_dynamic_develop_985_01.01.2016.ITS1.fasta -nonchimeras leaf_notchim.fna -strand plus -uchimeout leafchim_log.txt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "usearch v8.1.1861_i86linux64, 74.2Gb RAM, 12 cores\n",
    "(C) Copyright 2013-15 Robert C. Edgar, all rights reserved.\n",
    "http://drive5.com/usearch\n",
    "\n",
    "Licensed to: ronh@molbio.uoregon.edu\n",
    "\n",
    "00:36 3.5Gb  100.0% Reading leaf_defloat.fna\n",
    "00:36 3.5Gb    0.1% Reading /home6/dthomas/combobiom/uchime_reference_dataset_01.01.2016/ITS1_ITS2_datasets/uchime_sh_refs_dynamic_de$\n",
    "00:36 3.5Gb  100.0% Reading /home6/dthomas/combobiom/uchime_reference_dataset_01.01.2016/ITS1_ITS2_datasets/uchime_sh_refs_dynamic_de$\n",
    "elop_985_01.01.2016.ITS1.fasta\n",
    "00:36 3.5Gb    0.1% Masking                                                                                                           00:36 3.5Gb  100.0% Masking   \n",
    "00:37 3.5Gb  100.0% Word stats\n",
    "00:37 3.5Gb  100.0% Alloc rows\n",
    "00:37 3.5Gb  100.0% Build index\n",
    "40:45 32.6Gb  100.0% Found 53858/10820734 chimeras (0.5%), 7478748 not classified (69.1%)\n",
    "42:49 32.6Gb  100.0% Writing hits                                                        \n",
    "44:48 32.6Gb  100.0% Writing 10766876 non-chimeras and unclassifieds"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "check leaf chimera outputs, calculate loss below as for wood"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "How many reads were chimeric? In the wood:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": []
    }
   ],
   "source": [
    "grep '>' wood_defloat.fasta | wc -l\n",
    "grep '>' wood_notchim.fasta | wc -l"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "4548468<br>\n",
    "4530776<br>\n",
    "Looks like we lost 17,692 reads (~0.4%) of the wood reads from the previous step.\n",
    "\n",
    "Leaves:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": []
    }
   ],
   "source": [
    "grep '>' leaf_defloat.fna | wc -l\n",
    "grep '>' leaf_notchim.fna | wc -l "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "10820734<br>\n",
    "10766876<br>\n",
    "We lost 53,858 reads (0.5%) of leaf reads from the previous step."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='combine'></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Combining fasta files from all 3 studies"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So, the whole point of this exercise is get our three types of sequence files (Leaf endophytes, wood endophytes, and stromata) into a single fasta file that can be used to create OTU clusters. To help this, we'll simplify our leaf and wood identifier labels, extract ITS1 from our stromata, trim our illumina reads to ITS1 region manually, and concatenate the three fasta files. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Simplifying leaf endophytes identifiers:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note: I omitted a step in this section. The utax command, which we use to generate our [biom table](#makebiom), only allows alpha-numeric characters in names of samples. This is an issue for us, the leaf reads have periods in some of the names. Stromata sequences also require a single sample name to differentiate them from the other types of sequences. Corrected [here](#underscores) -dan."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Leaf read identifiers look like this:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": []
    }
   ],
   "source": [
    "head -n 1 leaf_notchim.fna"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\\>78Leaf_5"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Use gnu [SED](https://www.gnu.org/software/sed/) to strip these down and reformat them:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": []
    }
   ],
   "source": [
    "sed s/_.*//g leaf_notchim.fna > leaf_relab.fna"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Leaf identifiers just contain sample number and study (host) info. They look like this:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": []
    }
   ],
   "source": [
    "head -n 1 leaf_relab.fna"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\\>78Leaf"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Simplifying wood endophytes identifiers:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Wood identifiers look this:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": []
    }
   ],
   "source": [
    "head -n 1 wood_notchim.fasta"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\\>M01498:244:000000000-ANT97:1:1101:17999:1109 1:N:0:160"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here's a series of regexes that seems to work for simplifying them:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": []
    }
   ],
   "source": [
    "sed s/^\\>.*[0-9]\\ //g wood_notchim.fasta | sed s/1:[YN]:0:/\\>/g | sed '/^>/ s/$/wood/g' > wood_relab.fasta"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Looks like this:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": []
    }
   ],
   "source": [
    "head -n 1 wood_relab.fasta"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\\>160wood"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Extract ITS1 from stromata full ITS sequences"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "These stromatal reads are hand-aligned and curated by Roo, but their origins vary, as do their primer sites. To make them comparable during otu clustering, we'll extract the ITS1 region from them. The ITS1 region from stromata reads will be extracted using [Bengtsson-Palme et al.'s (Nilsson's) ITS extractor](http://microbiology.se/software/itsx/). \n",
    "\n",
    "The installation of the ITS extractor was theoretically simple but a little bit of a pain in the ass. It was simple in that you simply need to download a compressed, archived package from the [ITS extractor website](http://microbiology.se/software/itsx/), and it's main dependency, the [HMMER package](http://hmmer.org/). \n",
    "\n",
    "Each of these packages contains binaries that worked with my linux setups, but initial attempts at using the command line programs returned a lot of errors, so had do several things: delete old databases [(details here)](http://microbiology.se/2013/07/08/metaxa-and-hmmer-3-1b/) and use \"--reset T\" flag to create new dbs. I also had use the \"hmmpress\" command from the HMMER package with every one of these .hmm files, as a superuser, to create databases with correct formats. \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": []
    }
   ],
   "source": [
    "./ITSx -i ../Final_Stromata_Ref_Sequences.fasta -o strom --preserve T --allow_single_domain -t F"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Some notes about the settings used here:\n",
    "\n",
    "--preserve T  == keep our identifier lines, don't let ITSx create new ones \n",
    "\n",
    "--allow_single_domain == if ITSx can only find one conserved region (18s, 5.8, or 28s) to anchor into, this is enough (usually require two)\n",
    "\n",
    "--t F == look at fungal reads only\n",
    "\n",
    "--multi_thread 12 == use 12 cores, cuz we got 'em."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "ITS1 from all 51 stromatal sequences were extracted. Problems reported by ITSx:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "A_purpureonitens_consesus\tEnd of SSU sequence not found; Start of LSU sequence not found\r\n",
      "H_perforatum_China\tEnd of SSU sequence not found\r\n",
      "X_enterogena_Ecuador\tStart of LSU sequence not found\r\n",
      "X_flabelliforme_australia_thai\tStart of LSU sequence not found\r\n",
      "X_flabelliforme_china\tStart of LSU sequence not found\r\n"
     ]
    }
   ],
   "source": [
    "cat strom.problematic.txt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='trim3'></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Trimming leaf and wood illumina reads to ITS1 region"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we need to get rid of SSU and 5.8s region sequences from our Illumina reads. Manual inspection of our stromata sequences, and the ITSx mapping of a random subset of our leaf and wood illumina datasets both show a constant distance of 46 base pairs from the end of our forward (ITS1F) primer and the beginning of the ITS region, and 30 bp after the end of the ITS1 region we hit our ITS2 primer site. \n",
    "\n",
    "Each run of ITSx generates a positions file which tells us where SSU, 5.8s, and LSU regions (not present in our illumina reads) are located:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": []
    }
   ],
   "source": [
    "## get 100 randomly selected reads from our illumina fastas, after removing linebreaks\n",
    "./subset_fasta.py leaf_relab_nolb.fna 100 leaf_sub.fna\n",
    "./subset_fasta.py wood_relab_nolb.fasta 100 wood_sub.fasta\n",
    "cat leaf_sub.fna wood_sub.fasta > combo_sub.fasta\n",
    "## have ITSx program take a look at them...\n",
    "./ITSx -i ../combo_sub.fasta -o combo_sub --preserve T --allow_single_domain -t F"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "130Leaf\t330 bp.\tSSU: 1-46\tITS1: 47-300\t5.8S: No end\tITS2: Not found\tLSU: Not found\tBroken or partial sequence, only partial 5.8S! \r\n",
      "38Leaf\t315 bp.\tSSU: 1-46\tITS1: 47-285\t5.8S: No end\tITS2: Not found\tLSU: Not found\tBroken or partial sequence, only partial 5.8S! \r\n",
      "127Leaf\t315 bp.\tSSU: 1-46\tITS1: 47-285\t5.8S: No end\tITS2: Not found\tLSU: Not found\tBroken or partial sequence, only partial 5.8S! \r\n",
      "102Leaf\t262 bp.\tSSU: 1-46\tITS1: 47-226\t5.8S: No end\tITS2: Not found\tLSU: Not found\tBroken or partial sequence, only partial 5.8S! \r\n",
      "95Leaf\t257 bp.\tSSU: 1-46\tITS1: 47-227\t5.8S: No end\tITS2: Not found\tLSU: Not found\tBroken or partial sequence, only partial 5.8S! \r\n",
      "96Leaf\t315 bp.\tSSU: 1-46\tITS1: 47-227\t5.8S: No end\tITS2: Not found\tLSU: Not found\tBroken or partial sequence, only partial 5.8S! \r\n",
      "24Leaf\t256 bp.\tSSU: 1-46\tITS1: 47-226\t5.8S: No end\tITS2: Not found\tLSU: Not found\tBroken or partial sequence, only partial 5.8S! \r\n",
      "117Leaf\t254 bp.\tSSU: 1-46\tITS1: 47-224\t5.8S: No end\tITS2: Not found\tLSU: Not found\tBroken or partial sequence, only partial 5.8S! \r\n",
      "15Leaf\t243 bp.\tSSU: 1-46\tITS1: 47-213\t5.8S: No end\tITS2: Not found\tLSU: Not found\tBroken or partial sequence, only partial 5.8S! \r\n",
      "32Leaf\t250 bp.\tSSU: 1-46\tITS1: 47-220\t5.8S: No end\tITS2: Not found\tLSU: Not found\tBroken or partial sequence, only partial 5.8S! \r\n",
      "245wood\t237 bp.\tSSU: 1-46\tITS1: 47-204\t5.8S: No end\tITS2: Not found\tLSU: Not found\tBroken or partial sequence, only partial 5.8S! \r\n",
      "246wood\t254 bp.\tSSU: 1-46\tITS1: 47-224\t5.8S: No end\tITS2: Not found\tLSU: Not found\tBroken or partial sequence, only partial 5.8S! \r\n",
      "247wood\t244 bp.\tSSU: 1-46\tITS1: 47-219\t5.8S: No end\tITS2: Not found\tLSU: Not found\tBroken or partial sequence, only partial 5.8S! \r\n",
      "248wood\t224 bp.\tSSU: 1-46\tITS1: 47-194\t5.8S: No end\tITS2: Not found\tLSU: Not found\tBroken or partial sequence, only partial 5.8S! \r\n",
      "250wood\t268 bp.\tSSU: 1-46\tITS1: 47-238\t5.8S: No end\tITS2: Not found\tLSU: Not found\tBroken or partial sequence, only partial 5.8S! \r\n",
      "251wood\t226 bp.\tSSU: 1-46\tITS1: 47-196\t5.8S: No end\tITS2: Not found\tLSU: Not found\tBroken or partial sequence, only partial 5.8S! \r\n",
      "253wood\t219 bp.\tSSU: 1-47\tITS1: 48-189\t5.8S: No end\tITS2: Not found\tLSU: Not found\tBroken or partial sequence, only partial 5.8S! \r\n",
      "254wood\t296 bp.\tSSU: 1-98\tITS1: 99-266\t5.8S: No end\tITS2: Not found\tLSU: Not found\tBroken or partial sequence, only partial 5.8S! \r\n",
      "256wood\t285 bp.\tSSU: 1-46\tITS1: 47-255\t5.8S: No end\tITS2: Not found\tLSU: Not found\tBroken or partial sequence, only partial 5.8S! \r\n",
      "257wood\t256 bp.\tSSU: 1-46\tITS1: 47-226\t5.8S: No end\tITS2: Not found\tLSU: Not found\tBroken or partial sequence, only partial 5.8S! \r\n"
     ]
    }
   ],
   "source": [
    "head combo_sub.positions.txt; tail combo_sub.positions.txt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And so on. Almost all reads examined contain 46 bp of SSU DNA, then ITS1 region starts. Also notice that the ITS1 region ends 30 bp before the end of the read, so our ITS2 primers are located 30 bp into the 5.8s region. Knowing this, we can trim these ends to just the ITS1 region:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "## not sure if fastx tools can handle linebreaks. So remove them with a script. \n",
    "\n",
    "## wood\n",
    "./fasta_remove_linebreaks.py wood_relab.fasta wood_relab_nolb.fasta\n",
    "\n",
    "## leaves\n",
    "./fasta_remove_linebreaks.py leaf_relab.fna leaf_relab_nolb.fna"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": []
    }
   ],
   "source": [
    "## The wood:\n",
    "fastx_trimmer -f 47 -i wood_relab_nolb.fasta | fastx_trimmer -t 30 -o wood.ITS1.fasta\n",
    "\n",
    "## the leaves:\n",
    "fastx_trimmer -f 47 -i leaf_relab_nolb.fna | fastx_trimmer -t 30 -o leaf.ITS1.fasta"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Combining wood and leaf fasta files"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='underscores'></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we combine our leaf and wood endophyte fastas into a single fasta. We will do this in two ways: \n",
    "\n",
    "(1) all reads from all three libraries will be combined, this is our data, which will await taxonomic assignments and sample data. For this let's add a uniform \"Strom\" label to the strom' to all the stromatal sequences, this is important for making the biom table [below](#makebiom). We'll also change out the periods in the leaf reads for underscores, this messes with the utax command. This probably should have been done when we [simplified illumina identifiers](#combine). Anyway, sed to rescue again, on the ITS1 region of stromata sequences."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": []
    }
   ],
   "source": [
    "sed '/^>/  s/>/>Strom;/' strom.ITS1.fasta > strom_add_to_combo.fasta"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "sed '/>/ s/\\./_/' leaf.ITS1.fasta -i"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now combine them:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "cat wood.ITS1.fasta leaf.ITS1.fasta strom_add_to_combo.fasta > combo.fasta"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "(2) However we also need to generate an OTU repset as part of our tool box for making taxonomic assignments. As part of this, leaf and wood reads need to be combined, dereplicated, then stromatal sequences added after removing singletons from our illumina sets. For our OTU representative set, we wait to add in our stromata, because we don't want to lose these in the singleton removal [below](#removesings)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": []
    }
   ],
   "source": [
    "cat leaf.ITS1.fasta wood.ITS1.fasta > comboLW.fasta"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='otuclust'></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## OTU clustering"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We're following the [UPARSE pipeline](http://drive5.com/usearch/manual/uparse_pipeline.html) recommendations as much as possible. They recommend, for otu clustering, to remove singletons from our reads. They can be added back in later, during the final steps of [making our biome table](), [according to drive5](http://drive5.com/usearch/manual/mapreadstootus.html). But here, as we generate our OTU rep set, they can form spurious OTUs. So out they go. We're not removing our carefully curated stromatal sequences, though, as the point of excluding singletons here is that they have a high probability of being erroneous. Also, we want our stromata to inform the clustering process.\n",
    "\n",
    "First we [dereplicate](http://drive5.com/usearch/manual/cmd_derep_fulllength.html) our reads, takes a few seconds:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": []
    }
   ],
   "source": [
    "usearch -derep_fulllength comboLW.fasta -fastaout combo_derepLW.fasta -sizeout"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='removesings'></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "usearch v8.1.1861_i86linux64, 74.2Gb RAM, 12 cores\n",
    "(C) Copyright 2013-15 Robert C. Edgar, all rights reserved.\n",
    "http://drive5.com/usearch\n",
    "\n",
    "Licensed to: ronh@molbio.uoregon.edu\n",
    "\n",
    "00:36 3.8Gb  100.0% Reading comboLW.fasta 00:03 353Mb    8.5% Reading comboLW.fasta\n",
    "00:45 7.8Gb 15295207 seqs, 2304847 uniques, 1836384 singletons (79.7%)\n",
    "00:45 7.8Gb Min size 1, median 1, max 416662, avg 6.64\n",
    "01:01 4.8Gb  100.0% Writing combo_derepLW.fasta"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next we [sort](http://drive5.com/usearch/manual/cmd_sortbysize.html) these dereplicated reads by size. In the process, we remove singletons."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": []
    }
   ],
   "source": [
    "usearch -sortbysize combo_derepLW.fasta -fastaout combo_sorted.fasta -minsize 2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "usearch v8.1.1861_i86linux64, 74.2Gb RAM, 12 cores\n",
    "(C) Copyright 2013-15 Robert C. Edgar, all rights reserved.\n",
    "http://drive5.com/usearch\n",
    "\n",
    "Licensed to: ronh@molbio.uoregon.edu\n",
    "\n",
    "00:36 3.8Gb  100.0% Reading comboLW.fasta[[B^[[B^[[B^[[B^[[B^[[B^[[B^[[B00:03 353Mb    8.5% Reading comboLW.fasta\n",
    "00:45 7.8Gb 15295207 seqs, 2304847 uniques, 1836384 singletons (79.7%)\n",
    "00:45 7.8Gb Min size 1, median 1, max 416662, avg 6.64\n",
    "01:01 4.8Gb  100.0% Writing combo_derepLW.fasta\n",
    "[dthomas@cn113 combobiom]$ usearch -sortbysize combo_derepLW.fasta -fastaout combo_sorted.fasta -minsize 2\n",
    "usearch v8.1.1861_i86linux64, 74.2Gb RAM, 12 cores\n",
    "(C) Copyright 2013-15 Robert C. Edgar, all rights reserved.\n",
    "http://drive5.com/usearch\n",
    "\n",
    "Licensed to: ronh@molbio.uoregon.edu\n",
    "\n",
    "00:06 639Mb  100.0% Reading combo_derepLW.fasta\n",
    "00:06 605Mb Getting sizes                      \n",
    "00:09 624Mb Sorting 468463 sequences\n",
    "00:11 626Mb  100.0% Writing output"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "How many of our sequences were singletons? "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": []
    }
   ],
   "source": [
    "grep '>' combo_derepLW.fasta | wc -l\n",
    "grep '>' combo_sortedLW.fasta | wc -l"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "2304847<br>\n",
    "464475"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Wow, a lot of singletons. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Time to add the stromata back into this fasta file. Let's change the stromata labels to match the dereplicated, sorted illumina read labels. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": []
    }
   ],
   "source": [
    "## sed, add \"size 100000 to end of each identifier\n",
    "sed '/^>/ s/$/;size=1000000;/' strom.ITS1.fasta > strom.ITS1.relab.fasta"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      ">A_aff_atroroseum;size=1;\r\n",
      "GCGAGTTAGCAAAACTCCAAAACCCTTTGTGAACCTTACCGTCGTTGCCTCGGCGTGTGCCGCGGCTACCCTGGAGTAGTTACCCTGGACAGGTTACCCTATAGGGGCTACCCTGGAGGGGTTCCTACCCTGGAAGCCGGCACCCGGCCCGCCAAAGGACCCGTACAAAATTCTGTCTTACCAGTGTATCTCTGAATGCTTCAACTGAAATAAGTTA\r\n"
     ]
    }
   ],
   "source": [
    "head strom.ITS1.relab.fasta -n 2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": []
    }
   ],
   "source": [
    "cat strom.ITS1.relab.fasta combo_sortedLW.fasta > combo_sorted.fasta"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And then give this combined file to usearch to chew on. We'll do a .97 and .95 similarity radius. Since we're going outside of the normal 0.97 radius, and because Roo's other analyses in this paper were done using an older version of usearch  (uclust) at the .95 radius, we'll use the \"usearch -cluster_smallmem\" command. I believe this is the closest descendant of the older algorithms used by usearch, and the newer algorithms haven't been tested for other radii. See discussion [here](http://drive5.com/usearch/manual/uparse_otu_radius.html)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "usearch -cluster_smallmem combo_sorted.fasta -id 0.95 -centroids otus_95_combo.fasta -sizein -sizeout -sortedby size\n",
    "usearch -cluster_smallmem combo_sorted.fasta -id 0.97 -centroids otus_97_combo.fasta -sizein -sizeout -sortedby size"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "usearch v8.1.1861_i86linux64, 74.2Gb RAM, 12 cores\n",
    "(C) Copyright 2013-15 Robert C. Edgar, all rights reserved.\n",
    "http://drive5.com/usearch\n",
    "\n",
    "Licensed to: ronh@molbio.uoregon.edu\n",
    "\n",
    "00:30  61Mb  100.0% 10302 clusters, max size 1099864, avg 6247.7\n",
    "00:30  61Mb  100.0% Writing centroids to otus_95_combo.fasta    \n",
    "                                                            \n",
    "      Seqs  464526 (464526 ()\n",
    "  Clusters  10302 (10302 (10302 (103)\n",
    "  Max size  1099864 (1099864 (1099864 (2 (103))\n",
    "  Avg size  6247.7\n",
    "  Min size  2\n",
    "Singletons  0, 0.0% of seqs, 0.0% of clusters\n",
    "   Max mem  61Mb\n",
    "      Time  30.0s\n",
    "Throughput  15.5k seqs/sec."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "That was freakishly fast, a couple seconds each with the computing cluster. But I can't see any problems in the output. I think removing all of the singletons for this process vastly decreases the processing time (we went from ~ 1 million unique reads to ~1/4 million)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='UNITE'></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "# Customize UNITE database"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[UNITE](http://www2.dpes.gu.se/project/unite/UNITE_intro.htm) is an attempt by mycologists to get some quality control on fungal accessions in large public databases like GenBank. We will use it to identify our sequences. \n",
    "\n",
    "It's handy, but last I checked, there were some problems with sequences only identified to very low taxonomic resolution, e. g. accessions identified only to phylum. The cause problems because BLAST may display these matches over other very close matches that are more completely identified, if they match ever-so-slightly-better to the poorly identified sequences.\n",
    "\n",
    "I will use usearch's utax algorithm. So let's get the UNITE database made expressly for this:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "wget https://unite.ut.ee/sh_files/utax_reference_dataset_31.01.2016.zip"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "I had lot of trouble getting my python interpreter to work on the database due issues with special characters like umlauts and other symbols (especially \"ë\" and \"×\" that are used in some names of fungi in UNITE. My versions of sed and iconv and my python interpeter would not work with UNITE names. My $LANG = en_US.UTF-8, which seems like it should be able to handle these characters, but character encoding is a complex subject which I do not understand. Most of the database is in simple unicode8, though, so I was able to use my favorite text editor VIM to substitute these characters out:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "unzip utax_reference_dataset_31.01.2016.zip\n",
    "vim utax_reference_dataset_31.01.2016.fasta"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Inside vim command mode (press ':'):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    ":%s/ë/e/g\n",
    ":%s/×/x/g"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This fixed it for me. Good luck with your unicode. Maybe mac users won't have this issue. Exit vim, back in the shell. Here I've written another python script, \"UNITE_with_class2.py\", for only keeping those UNITE records with class identification or better. This is in the repository. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "./UNITE_with_class2.py utax_reference_dataset_31.01.2016.fasta utax_ref_class.fasta"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We'll need to append our curated stromata to this database. The usearch pipeline has a standard heading needed to create a searchable database from fasta files, described [here](http://drive5.com/usearch/manual/tax_annot.html). Since the UNITE database as we downloaded it is already formatted to these specs, we just need to reformat our stromata sequence fasta. More SED magic to add in most of the taxonomy:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": []
    }
   ],
   "source": [
    "sed '/^>/ s/>/,s:/g' Final_Stromata_Ref_Sequences.fasta | sed '/,s:/ s/^/>RefStrom;tax=d:Fungi,p:Ascomycota,c:Sordariomycetes,o:Xylariales,f:Xylariaceae/' > strom_add_to_UNITE.fasta"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Expand our genus names:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": []
    }
   ],
   "source": [
    "sed 's/s:A_/g:Annulohypoxylon,s:/g' strom_add_to_UNITE.fasta -i\n",
    "sed 's/s:B_/g:Biscogniauxia,s:/g' strom_add_to_UNITE.fasta -i\n",
    "sed 's/s:H_/g:Hypoxylon,s:/g' strom_add_to_UNITE.fasta -i\n",
    "sed 's/s:K_/g:Kretzschmaria,s:/g' strom_add_to_UNITE.fasta -i\n",
    "sed 's/s:N_/g:Nemania,s:/g' strom_add_to_UNITE.fasta -i\n",
    "sed 's/s:W_/g:Whalleya,s:/g' strom_add_to_UNITE.fasta -i\n",
    "sed 's/s:X_/g:Xylaria,s:/g' strom_add_to_UNITE.fasta -i"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Our new headers look like this:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      ">RefStrom;tax=d:Fungi,p:Ascomycota,c:Sordariomycetes,o:Xylariales,f:Xylariaceae,g:Annulohypoxylon,s:aff_atroroseum\r\n",
      ">RefStrom;tax=d:Fungi,p:Ascomycota,c:Sordariomycetes,o:Xylariales,f:Xylariaceae,g:Annulohypoxylon,s:aff_stygium\r\n",
      ">RefStrom;tax=d:Fungi,p:Ascomycota,c:Sordariomycetes,o:Xylariales,f:Xylariaceae,g:Annulohypoxylon,s:atroroseum\r\n",
      ">RefStrom;tax=d:Fungi,p:Ascomycota,c:Sordariomycetes,o:Xylariales,f:Xylariaceae,g:Annulohypoxylon,s:bovei_var_microspora\r\n",
      ">RefStrom;tax=d:Fungi,p:Ascomycota,c:Sordariomycetes,o:Xylariales,f:Xylariaceae,g:Annulohypoxylon,s:moriforme\r\n",
      ">RefStrom;tax=d:Fungi,p:Ascomycota,c:Sordariomycetes,o:Xylariales,f:Xylariaceae,g:Xylaria,s:sp_nov_1_long\r\n",
      ">RefStrom;tax=d:Fungi,p:Ascomycota,c:Sordariomycetes,o:Xylariales,f:Xylariaceae,g:Xylaria,s:sp_nov_1_short\r\n",
      ">RefStrom;tax=d:Fungi,p:Ascomycota,c:Sordariomycetes,o:Xylariales,f:Xylariaceae,g:Xylaria,s:sp_nov_2\r\n",
      ">RefStrom;tax=d:Fungi,p:Ascomycota,c:Sordariomycetes,o:Xylariales,f:Xylariaceae,g:Xylaria,s:telfairii\r\n",
      ">RefStrom;tax=d:Fungi,p:Ascomycota,c:Sordariomycetes,o:Xylariales,f:Xylariaceae,g:Xylaria,s:telfairii_Ecuador\r\n"
     ]
    }
   ],
   "source": [
    "sed '1~2p' strom_add_to_UNITE.fasta -n | head -n 5 \n",
    "sed '1~2p' strom_add_to_UNITE.fasta -n | tail -n 5"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Compare this to the already-formatted UNITE identifiers:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      ">EU821669|SH188517.07FU;tax=d:Fungi,p:Basidiomycota,c:Agaricomycetes,o:Agaricales,f:Cortinariaceae,g:Cortinarius,s:Cortinarius_balaustinus_SH188517.07FU;\r\n",
      ">AF201716|SH204808.07FU;tax=d:Fungi,p:Ascomycota,c:Sordariomycetes,o:Xylariales,f:Xylariaceae,g:Annulohypoxylon,s:Annulohypoxylon_truncatum_SH204808.07FU;\r\n",
      "sed: couldn't write 145 items to stdout: Broken pipe\r\n"
     ]
    }
   ],
   "source": [
    "sed '1~2p' utax_ref_class.fasta -n | head -n 2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Okay, add this to our UNITE database:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": []
    }
   ],
   "source": [
    "cat strom_add_to_UNITE.fasta utax_ref_class.fasta > utax-strom_ref.fasta"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And convert this \"classy\" database to a format that usearch/utax likes:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "usearch -makeudb_utax utax-strom_ref.fasta -output utax-strom_ref.udb -report dbmake_report.txt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='asstax'></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Assign taxonomy to otus"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Time to put some names on the otu reference set. To do this with the usearch pipeline, we need to give all of our otu clusters a unique name/number. Since we used [usearch -cluster_smallmem](http://www.drive5.com/usearch/manual/cmd_cluster_smallmem.html) instead of [usearch -cluster_otus](http://www.drive5.com/usearch/manual/cmd_cluster_otus.html) due to our .95 similarity radius for our otus, [see here for discussion](http://www.drive5.com/usearch/manual/uparse_otu_radius.html), we need to add our otu names ourselves. As everyone knows, I like to do things the hard way. Anyway, wrote a python script for this, called \"addOTUtag.py\", it's in the repository. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": []
    }
   ],
   "source": [
    "./addOTUtag.py otus_95_combo.fasta OTU otus_95_combo_relabel.fasta\n",
    "./addOTUtag.py otus_97_combo.fasta OTU otus_97_combo_relabel.fasta"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now use the database created above to assign taxonomy. We'll use the [usearch -utax](http://drive5.com/usearch/manual/cmd_utax.html) command."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": []
    }
   ],
   "source": [
    "usearch -utax otus_95_combo_relabel.fasta -db utax-strom_ref.udb -strand both -fastaout otus_95_combo_asstax.fasta\n",
    "usearch -utax otus_97_combo_relabel.fasta -db utax-strom_ref.udb -strand both -fastaout otus_97_combo_asstax.fasta"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "What does this look like? "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": []
    }
   ],
   "source": [
    "head otus_95_combo_asstax.fasta"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\\>OTU10:X_atrosphaerica;size=1000916;tax=d:Fungi(1.0000),p:Ascomycota(0.9997),c:Sordariomycetes(0.9995),o:Xylariales(0.9998),f:Xylariaceae(0.9997),g:Xylaria(0.9994),s:atrosphaerica;\n",
    "AAGAGTGTTATAACTCCTAAACCCATGTGAATATACCGTATGTTGCTTCGGTAGGCTGGATCTATCCCGTGGGGTTCCCT\n",
    "TGGGAGATGTGGAAAAGCCTGCCGGCGGCCTATTTAATCCTGTTTATTTATTTTTGAATTCTGAGGTTATAATTAATCAG\n",
    "TTA\n",
    "\\>OTU1:X_hypoxylon_Oregon;size=1099864;tax=d:Fungi(1.0000),p:Ascomycota(1.0000),c:Sordariomycetes(0.9963),o:Xylariales(0.9992),f:Xylariaceae(0.9994),g:Xylaria(0.9990),s:hypoxylon_Oregon;\n",
    "AAGAGTTAATTACAACTCCCAAACCCATGTGAACTTACCTTCTGTTGCCTCGGCAGGTCGTGTTTACCCTGTGAGGTCCT\n",
    "ACCCTGTAGGACCCTACCTGGTAGACACGGGTACGCCTGCCGGTGGCCCATGAAACTCTGTTAATTCTATGTTATTCTGA\n",
    "ATCTATAACTAAATAAGTTA\n",
    "\\>OTU4:N_sp_1;size=1004185;tax=d:Fungi(1.0000),p:Ascomycota(0.9995),c:Sordariomycetes(0.9995),o:Xylariales(0.9998),f:Xylariaceae(0.9996),g:Nemania(0.9604),s:sp_1;\n",
    "GAGAGTTCTCATAAAACTCCCAAACCCATGTGAACATACCTTCTGTTGCCTCGGCAGGTCGGCGCCGGCGCTCACCGGCG\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='makebiom'></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Make and tidy up biom tables"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So now we take these OTU clusters, with taxonomic assignments, and map our reads to them to create a biom table. Using the [usearch -usearch_global](http://drive5.com/usearch/manual/mapreadstootus.html) command. This command outputs the [first generation format](http://biom-format.org/documentation/biom_format.html)  of biom tables, a json table that is human readable. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "usearch -usearch_global combo.fasta -db otus_97_combo_asstax.fasta -strand both -id 0.97 -biomout combo_otu_97.biom\n",
    "usearch -usearch_global combo.fasta -db otus_95_combo_asstax.fasta -strand both -id 0.95 -biomout combo_otu_95.biom"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This makes our two biome tables. Log files for these report that 99.5% of our quality filtered reads mapped to an OTU in the 97% radius OTU rep set and 99.9% mapped to an OTU in our 95% radius OTU rep set. \n",
    "\n",
    "What do these look like? "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 61,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "{\r\n",
      "\t\"id\":\"combo_otu_95.biom\",\r\n",
      "\t\"format\": \"Biological Observation Matrix 1.0\",\r\n",
      "\t\"format_url\": \"http://biom-format.org\",\r\n",
      "\t\"generated_by\": \"usearch8.1.1861\",\r\n",
      "\t\"type\": \"OTU table\",\r\n",
      "\t\"date\": \"Thu Sep  8 21:40:16 2016\",\r\n",
      "\t\"matrix_type\": \"sparse\",\r\n",
      "\t\"matrix_element_type\": \"float\",\r\n",
      "\t\"shape\": [10269,232],\r\n",
      "\t\"rows\":[\r\n",
      "\t\t{\"id\":\"OTU17:114Leaf\", \"metadata\":{\"taxonomy\":\"d:Fungi(1.0000),p:Ascomycota(1.0000),c:Sordariomycetes(0.9972),o:Xylariales(0.9992),f:Xylariaceae(0.9993),g:Xylaria(0.9990),s:hypoxylon_Oregon\"}},\r\n",
      "\t\t{\"id\":\"OTU7811:160wood\", \"metadata\":{\"taxonomy\":\"d:Fungi(1.0000),p:Ascomycota(0.9855),c:Dothideomycetes(0.9551),o:Pleosporales(0.7729)\"}},\r\n"
     ]
    }
   ],
   "source": [
    "head -n 13 combo_otu_95.biom"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 62,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "{\r\n",
      "\t\"id\":\"combo_otu_97.biom\",\r\n",
      "\t\"format\": \"Biological Observation Matrix 1.0\",\r\n",
      "\t\"format_url\": \"http://biom-format.org\",\r\n",
      "\t\"generated_by\": \"usearch8.1.1861\",\r\n",
      "\t\"type\": \"OTU table\",\r\n",
      "\t\"date\": \"Thu Sep  8 21:40:29 2016\",\r\n",
      "\t\"matrix_type\": \"sparse\",\r\n",
      "\t\"matrix_element_type\": \"float\",\r\n",
      "\t\"shape\": [12558,232],\r\n",
      "\t\"rows\":[\r\n",
      "\t\t{\"id\":\"OTU16:114Leaf\", \"metadata\":{\"taxonomy\":\"d:Fungi(1.0000),p:Ascomycota(1.0000),c:Sordariomycetes(0.9972),o:Xylariales(0.9992),f:Xylariaceae(0.9993),g:Xylaria(0.9990),s:hypoxylon_Oregon\"}},\r\n",
      "\t\t{\"id\":\"OTU9350:160wood\", \"metadata\":{\"taxonomy\":\"d:Fungi(1.0000),p:Ascomycota(0.9855),c:Dothideomycetes(0.9551),o:Pleosporales(0.7729)\"}},\r\n"
     ]
    }
   ],
   "source": [
    "head -n 13 combo_otu_97.biom"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### clean up unidentified OTUs"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can check initial, obvious issues with our biome table by using the \"biom validate-table command\", part of the [biom-format python package](http://biom-format.org/index.html):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 63,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Invalid format 'Biological Observation Matrix 1.0', must be '1.0.0'\r\n",
      "'id' in {'id': '', 'metadata': {'taxonomy': 'd:Fungi(1.0000),p:Basidiomycota(0.9808),c:Agaricomycetes(0.9300),o:Agaricales(0.9049),f:Marasmiaceae(0.8875),g:Marasmius(0.7400),s:Marasmius_rotula_SH190961.07FU'}} appears empty\r\n",
      "Bad value at idx 0: [0, 0, 70749]\r\n",
      "Timestamp does not appear to be ISO 8601\r\n",
      "The input file is not a valid BIOM-formatted file.\r\n"
     ]
    }
   ],
   "source": [
    "biom validate-table -i combo_otu_95.biom"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "All of this appears to be pretty trivial, a missing id and a \"Bad value\". The \"Bad value\" happens often when I make biom tables, and the biom authors seem to think it is likely minor, due to having a floating decimal instead of an integer at some place in the table, discussion [here](https://github.com/biocore/biom-format/issues/701). The missing OTU is weird, but this seems to also happen every time I make a biom table. Let's see if we can find it:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 64,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\t\t{\"id\":\"\", \"metadata\":{\"taxonomy\":\"d:Fungi(1.0000),p:Basidiomycota(0.9808),c:Agaricomycetes(0.9300),o:Agaricales(0.9049),f:Marasmiaceae(0.8875),g:Marasmius(0.7400),s:Marasmius_rotula_SH190961.07FU\"}},\r\n"
     ]
    }
   ],
   "source": [
    "grep '\"id\":\"\"' combo_otu_95.biom"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can search our OTU taxonomic classifications to figure out what OTU this was:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 65,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      ">OTU1166:3Leaf;size=661;tax=d:Fungi(1.0000),p:Basidiomycota(0.9808),c:Agaricomycetes(0.9300),o:Agaricales(0.9049),f:Marasmiaceae(0.8875),g:Marasmius(0.7400),s:Marasmius_rotula_SH190961.07FU;\r\n"
     ]
    }
   ],
   "source": [
    "grep 'd:Fungi(1.0000),p:Basidiomycota(0.9808),c:Agaricomycetes(0.9300),o:Agaricales(0.9049),f:Marasmiaceae(0.8875),g:Marasmius(0.7400),s:Marasmius_rotula_SH190961.07FU' otus_95_combo_asstax.fasta"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So one OTU pops up, it's name is \"OTU1166:3Leaf\". We need to plug this into our json file. What would I do without SED?:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 66,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": []
    }
   ],
   "source": [
    "sed '/\"id\":\"\"/ s/\"id\":\"\",/\"id\":\"OTU1166:3Leaf\",/' combo_otu_95.biom -i"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Did that work?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 67,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0\r\n"
     ]
    }
   ],
   "source": [
    "grep '\"id\":\"\"' combo_otu_95.biom | wc -l"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 68,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Invalid format 'Biological Observation Matrix 1.0', must be '1.0.0'\r\n",
      "Bad value at idx 0: [0, 0, 70749]\r\n",
      "Timestamp does not appear to be ISO 8601\r\n",
      "The input file is not a valid BIOM-formatted file.\r\n"
     ]
    }
   ],
   "source": [
    "biom validate-table -i combo_otu_95.biom"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Grep and biom-format can't find any other empty ids. Still not sure why this happens, some minor bug in the usearch code? Something weird in my data? As mentioned above, I think the \"Bad value\" is a minor issue, but let's look at it:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 70,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\t\t{\"id\":\"67Leaf\", \"metadata\":null},\r\n",
      "\t\t{\"id\":\"112_1Leaf\", \"metadata\":null},\r\n",
      "\t\t{\"id\":\"113_2Leaf\", \"metadata\":null},\r\n",
      "\t\t{\"id\":\"74Leaf\", \"metadata\":null},\r\n",
      "\t\t{\"id\":\"113_1Leaf\", \"metadata\":null},\r\n",
      "\t\t{\"id\":\"126_2Leaf\", \"metadata\":null},\r\n",
      "\t\t{\"id\":\"91Leaf\", \"metadata\":null},\r\n",
      "\t\t{\"id\":\"Strom\", \"metadata\":null}\r\n",
      "\t],\r\n",
      "\t\"data\": [\r\n",
      "\t\t[0,0,70749],\r\n",
      "\t\t[0,3,3],\r\n",
      "\t\t[0,4,8],\r\n",
      "\t\t[0,5,58],\r\n",
      "\t\t[0,13,50],\r\n",
      "\t\t[0,14,87],\r\n",
      "\t\t[0,17,6],\r\n",
      "\t\t[0,21,12],\r\n",
      "\t\t[0,26,15],\r\n",
      "\t\t[0,27,67],\r\n",
      "\t\t[0,34,15],\r\n"
     ]
    }
   ],
   "source": [
    "grep \\\\[0,0,70749\\\\] -A 10 -B 10 combo_otu_95.biom"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Looks fine to me. This is a plentiful OTU (\"OTU17\", in sample 114leaf, 99,865 reads) and the sample (160w) has slightly more reads in it than this amount (71,220 total reads), so it seems like this value is possible. Not sure. Our .97 biom table also had an unidentified OTU and a Bad Value at the [0,0] (first data entry). I fixed the unidentified OTU in the same manner as above, don't know why the bad value report. We'll keep going for now, but we should check on this \"Bad value\" OTU and sample downstream, cuz something might be fishy.  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Change biom taxonomy metadata format"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Another clean-up step for adapting the uparse-outputted biom tables for downstream applications is reformatting the taxonomy info. Formatting these JSON tables is a little tricky, but there is a good example that I will follow [here](http://biom-format.org/documentation/adding_metadata.html). I wrote a script for this, \"format_tax.py\", in the repository. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 86,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": []
    }
   ],
   "source": [
    "./format_tax.py combo_otu_95.biom combo_95_relab.biom\n",
    "./format_tax.py combo_otu_97.biom combo_97_relab.biom"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And the taxonomy metadata now looks like this:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 87,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\t\"rows\":[\r\n",
      "\t\t{\"id\":\"OTU17:114Leaf\", \"metadata\":{\"taxonomy\": [\"k__Fungi\", \"p__Ascomycota\", \"c__Sordariomycetes\", \"o__Xylariales\", \"f__Xylariaceae\", \"g__Xylaria\", \"s__hypoxylon_Oregon\"]}},\r\n",
      "\t\t{\"id\":\"OTU7811:160wood\", \"metadata\":{\"taxonomy\": [\"k__Fungi\", \"p__Ascomycota\", \"c__Dothideomycetes\", \"o__Pleosporales\"]}},\r\n",
      "\t\t{\"id\":\"OTU2412:53Leaf\", \"metadata\":{\"taxonomy\": [\"k__Fungi\", \"p__Ascomycota\", \"c__Eurotiomycetes\", \"o__Chaetothyriales\"]}},\r\n",
      "\t\t{\"id\":\"OTU494:114Leaf\", \"metadata\":{\"taxonomy\": [\"k__Fungi\", \"p__Ascomycota\", \"c__Dothideomycetes\", \"o__Pleosporales\"]}},\r\n",
      "\t\t{\"id\":\"OTU66:65Leaf\", \"metadata\":{\"taxonomy\": [\"k__Fungi\", \"p__Ascomycota\", \"c__Dothideomycetes\", \"o__Capnodiales\", \"f__Mycosphaerellaceae\", \"g__Mycosphaerella\", \"s__Mycosphaerella_tassiana_SH216250.07FU\"]}},\r\n"
     ]
    }
   ],
   "source": [
    "grep 'rows' -A 5 combo_95_relab.biom"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='addmetadata'></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Add metadata"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So, theoretically, these biom tables should now parse well in downstream applications, like phyloseq. Let's add some information about our samples first. For the moment, we'll just get real sample numbers onto our columns. Up to this point, the wood samples have been identified using the numbers given by illumina software when they dumultiplexed them. These numbers are meaningless to our study, except to keep track of samples. We can map them by [adding metadata with biom scripts](http://biom-format.org/documentation/adding_metadata.html). The leaf reads were given correct sample names in their identifiers when we [demultiplexed them](#demult), but we'll need to fill this info into their metadata, also.\n",
    "\n",
    "For this we need a \"metadata mapping file.\" For the moment, I will use a very simple map, looks like this:  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "#SampleID\tSampleNumber\tLibrary\tSorC\r\n",
      "160wood\tDc-X\tW\tControl\r\n",
      "161wood\tDc-PosG\tW\tControl\r\n",
      "162wood\tDc-PosI\tW\tControl\r\n",
      "163wood\tDc-Neg\tW\tControl\r\n",
      "164wood\t1\tW\tSample\r\n",
      "165wood\t2\tW\tSample\r\n",
      "166wood\t3\tW\tSample\r\n",
      "167wood\t4\tW\tSample\r\n",
      "168wood\t5\tW\tSample\r\n",
      "127Leaf\t127\tL\tSample\r\n",
      "128Leaf\t128\tL\tSample\r\n",
      "129Leaf\t129\tL\tSample\r\n",
      "130Leaf\t130\tL\tSample\r\n",
      "131Leaf\t131\tL\tSample\r\n",
      "132Leaf\t132\tL\tSample\r\n",
      "133Leaf\t133\tL\tSample\r\n",
      "NC_1Leaf\tNC_1\tL\tControl\r\n",
      "NC_2Leaf\tNC_2\tL\tControl\r\n",
      "Strom\tS\tS\tControl\r\n"
     ]
    }
   ],
   "source": [
    "head combo_meta.tsv; tail combo_meta.tsv"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Note:** Tab-separated values seem to work better for adding metadata using biom scripts than comma-separated values, at least in my environment/versions. If this is an issue, a quick SED command can convert csv format to tsv:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": []
    }
   ],
   "source": [
    "sed 's/,/\\t/g' combo_meta.csv combo_meta.tsv "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So we use this map to add column metadata:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 97,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": []
    }
   ],
   "source": [
    "biom add-metadata -i combo_95_relab.biom -o combo_95_wMeta.biom -m combo_meta.tsv --output-as-json\n",
    "\n",
    "biom add-metadata -i combo_97_relab.biom -o combo_97_wMeta.biom -m combo_meta.tsv --output-as-json"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Did this work? A preliminary check with biom scripts."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\r\n",
      "The input file is a valid BIOM-formatted file.\r\n"
     ]
    }
   ],
   "source": [
    "biom validate-table -i combo_95_wMeta.biom"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\r\n",
      "The input file is a valid BIOM-formatted file.\r\n"
     ]
    }
   ],
   "source": [
    "biom validate-table -i combo_97_wMeta.biom"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "Okay! Looks like we have a biom table to run some stats on. We'll do this in a different notebook, just to keep things organized. "
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Bash",
   "language": "bash",
   "name": "bash"
  },
  "language_info": {
   "codemirror_mode": "shell",
   "file_extension": ".sh",
   "mimetype": "text/x-sh",
   "name": "bash"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}