notebooks/r4.2.Resistome-profiling.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Resistome profiling\n",
    "\n",
    "One of the main concerns during the suspected outbreak of *E.cloacae* was the possibility of antimicrobial resistance (AMR), particularly carbapenem resistance.\n",
    "\n",
    "For this workflow, we'll be sketching with:\n",
    "\n",
    "* MinHash\n",
    "* Minimizers\n",
    "\n",
    "***\n",
    "\n",
    "## Taking a quick look with MinHash\n",
    "\n",
    "To quickly tell if there are any known AMR genes in these samples, we can use [GROOT](https://github.com/will-rowe/groot). GROOT is designed for metagenomes but does the job for isolates too. It works by building variation graphs for clusters of genes, then indexing the traverslas in each graph using **MinHash** sketches.\n",
    "\n",
    "We're going to see how you can sketch reads as you download them, this means you could just keep the samples which contain the genes that you are intested in. GROOT also has read QC and trimming built in, so bad reads will be handled.\n",
    "\n",
    "* to begin, download a reference AMR database and index it:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# download the ResFinder database\n",
    "!groot get -d arg-annot\n",
    "\n",
    "# index the database\n",
    "!groot index -i ./arg-annot.90 -l 150 -k 21 -s 24 -o arg-annot-index"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> `-l` specifies the window length to sketch in the graph, which should be similar to the read length\n",
    "\n",
    "> `-k` specifies the k-mer length and `-s` specifies the sketch size\n",
    "\n",
    "> GROOT is quicker if we use multiple cpus to index or sketch (using groot `--processors`)\n",
    "\n",
    "* now we stream the data from the ENA; as reads arrive we align them to the reference graphs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# stream the reads and align them to the graphs\n",
    "!fastq-dump ERX168346 -Z --split-files | groot align -i arg-annot-index -o ERX168346-graphs > ERX168346.bam\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> we use fastq-dump to stream reads from the ENA download into the groot command (via the `-Z` flag and pipe `|`)\n",
    "\n",
    "> the align subcommand produces a BAM file containing all graph traversals for each read. Each BAM file essentially contains the ARG-derived reads.\n",
    "\n",
    "> the gfa variation graphs which had reads align are also kept and can be viewed in Bandage etc.\n",
    "\n",
    "* now, report what AMR genes are present in this sample:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!groot report -i ERX168346.bam --lowCov\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> the `--lowCov` flag is used as we are using GROOT on isolates, not metagenomes. The flag ignores uncovered reads in the first few bases of each gene, which are usually uncovered because there are not enough reads to completely span the gene (partial gene alignments aren't counted by GROOT).\n",
    "\n",
    "This result tells us that our isolates contain AMR genes and this warrants further inspection. Good job we have already QC'd the data and have this waiting for us.\n",
    "\n",
    "The [paper](https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4001082/) that describes this outbreak was particularly interested in beta-lactam resistance. They found *blaSHV-12*, *blaIMP-1* and *blaTEM-1* genes in **isolate EC1a** (ERX168346). This was the only *E.cloacae* isolate in which they found these genes, or did any phenotypic testing for.\n",
    "\n",
    "We managed to find *blaSHV-12* in our MinHash search but not *blaIMP-1* and *blaTEM-1*. If we look at the BAM file, we do find some reads matching these two genes but we didn't get enough coverage to call the genes. We also used the raw data. Let's now try using the cleaned data and doing a full read-alignment.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Mapping reads with MiniMap2\n",
    "\n",
    "The quick look with GROOT told us that these samples contain AMR genes. Let's say that we decided to keep these downloaded samples on disk and then quality checked them using our previous workflow.\n",
    "\n",
    "We will use [MiniMap2](https://github.com/lh3/minimap2) to map our quality checked reads to the reference AMR genes, which we will index with **minimizers**."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* get the fasta sequences for the AMR genes in the ARG-annot database:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!wget https://github.com/raw/will-rowe/groot/master/db/full-ARG-databases/arg-annot-db/argannot-args.fna"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* index the AMR genes and align our quality trimmed reads:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# create the index\n",
    "!minimap2 -x sr -d argannot-args.mmi argannot-args.fna\n",
    "# align the reads and output as a bam file\n",
    "!minimap2 -x sr -t 1 -a argannot-args.mmi ../data/reads/ERX168346_*-trimmed.fq.gz | samtools view -b -o ERX168346-minimap2.bam"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* sort, index and read in the alignment:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pysam\n",
    "pysam.sort(\"-o\", \"ERX168346-minimap2.bam\", \"ERX168346-minimap2.bam\")\n",
    "pysam.index(\"ERX168346-minimap2.bam\")\n",
    "# the rb arguments tell pysam that we are reading in a bam file\n",
    "samfile = pysam.AlignmentFile(\"ERX168346-minimap2.bam\", \"rb\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* let's look at the all the AMR genes that were at least 95% covered:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# check each AMR gene in the database\n",
    "for ref in samfile.header['SQ']:\n",
    "    name=ref['SN']\n",
    "    length=ref['LN']\n",
    "    # create a pilup for the reference AMR gene\n",
    "    pileup=samfile.pileup(name)\n",
    "    coveredBases=0\n",
    "    # see if any reads aligned to this AMR gene\n",
    "    for pos,column in enumerate(pileup,1):\n",
    "        coveredBases+=1\n",
    "\n",
    "    # if >95% of the AMR gene had reads align, print the name of the gene and it's coverage\n",
    "    coverage=(coveredBases/length)*100\n",
    "    if (coverage > 95):\n",
    "        print(\"{} is {}% covered by reads\" .format(name, coverage))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There is a carbapenemase gene (*blaIMP-34*)! But...\n",
    "\n",
    "We are still not finding the *blaTEM* reported in the paper. Also, some of our AMR genes which GROOT found are not showing up here.\n",
    "\n",
    "With regards to no *blaTEM* gene appearing, this may be due to our database being different to the one used in the paper. More likely though is that the paper used a 90% similarity threshold to annotate AMR genes by mapping genes to *de novo* assemblies of the isolates - this approach is probably not the best for resistome profiling but that is a benchmark for another day.\n",
    "\n",
    "With regards to GROOT not finding some of these genes. This is probably because we are using qc filtered reads here, whereas GROOT used the raw data. We can check this theory.\n",
    "\n",
    "* take a closer look at the alignments to blaIMP-34:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# the fetch method creates an iterator around the region we give\n",
    "iter = samfile.fetch(\"argannot~~~(Bla)IMP-34~~~AB715422:1696-2436\")\n",
    "# we can loop over the reads aligned to the specified region\n",
    "for read in iter:\n",
    "    print(read.cigarstring)\n",
    "samfile.close()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can see that the length of the reads aligning to the gene are pretty variable - the trimming has reduced some of these reads from 150bp down to 100bp. So in some cases, our reads had 33% low quality bases. This will have thrown our sketching based mapping, so reads won't have seeded. In GROOT's case, even if the reads seeded, it wouldn't report any alignments that needed 50bp hard clipping.\n",
    "\n",
    "You could go back and run GROOT again, this time using the qc filtered reads. Because of the variable read length, we would need to use `--containment` for our sketch searches.\n",
    "\n",
    "Hopefully this workflow illustrates the pros and cons of resistome profiling using sketching algorithms.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "***\n",
    "\n",
    "Let's move on to the next workflow in our outbreak analysis: [outbreak surveillance](r4.3.Outbreak-surveillance.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
}