notebooks/r4.3.Outbreak-surveillance.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Outbreak surveillance\n",
    "\n",
    "Let's now perform some outbreak surveillance.\n",
    "\n",
    "We will use sketching to create a dendrogram based on the whole-genome sequencing data for the isolates.\n",
    "\n",
    "\n",
    "For this workflow, we'll be sketching with:\n",
    "\n",
    "* MinHash\n",
    "\n",
    "***\n",
    "\n",
    "We'll be using [mashtree](https://github.com/lskatz/mashtree). This uses [MASH](https://github.com/lskatz/mashtree) to estimate genetic distances for each isolate, based on **MinHash** estimates of Jaccard similarity. It then builds a tree....\n",
    "\n",
    "\n",
    "* first, we need to interleave the read pairs:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash\n",
    "for i in ../data/reads/*_1*;\n",
    "do\n",
    "base=${i%%_1-trimmed.fq.gz}\n",
    "interleave-reads.py ${base}* -o ${base}.fastq.gz\n",
    "done"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* run mashtree to get distances and create a newick tree:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!mashtree.pl --mindepth 5 --sketch-size 10000 --kmerlength 31 ../data/reads/*fastq.gz ../data/GCF_000025565.1_ASM2556v1_genomic.fna.gz > mashtree.dnd\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* now we have a tree, let's check it out:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from ete3 import Tree\n",
    "# Load a tree structure from a newick file.\n",
    "tree = Tree(\"mashtree.dnd\")\n",
    "tree.set_outgroup(\"GCF_000025565.1_ASM2556v1_genomic\")\n",
    "print(tree)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* compare this to the tree in the [paper](https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4001082), built using SNPs from whole genome sequence alignments:\n",
    "\n",
    "![](https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4001082/bin/emss-58061-f0001.jpg)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* and here are the isolate names -> ENA experiment again, for your reference:\n",
    "\n",
    "|isolate name|patient ID|ENA sample|ENA experiment|\n",
    "|-|-|-|-|\n",
    "EC1a|EC1|ERS184249|ERX168346|\n",
    "EC2a|EC2|ERS184250|ERX168347|\n",
    "EC2b|EC2|ERS184251|ERX168341|\n",
    "EC3a|EC3|ERS184252|ERX168340|\n",
    "EC4a|EC4|ERS184245|ERX168345|\n",
    "EC5a|EC5|ERS184246|ERX168339|\n",
    "EC6a|EC6|ERS184247|ERX168343|\n",
    "EC7a|EC7|ERS184248|ERX168344|"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can see that the outbreak cluster, which was identified in the paper as being EC1a, EC2a, EC2b and EC3a, can be identified using our sketch-based tree (the bottom clade). Here's our tree one more time:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(tree)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "***\n",
    "\n",
    "Let's move on to the next workflow in our outbreak analysis: [data mining](r4.4.Data-mining.ipynb)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}