notebooks/r3.1.Set-similarity-with-MinHash.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 3.1. Set similarity with MinHash"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Recap\n",
    "\n",
    "As we covered in the review, there are several different *flavours* of MinHash (e.g. KMV, Bottom-k). We will use **KMV MinHash** for now, which populates a sketch using the Minimum Values from K random permutations of the input data.\n",
    "\n",
    "It gets a bit confusing when we start talking about k-mers as well as KMV sketches, so from now on we will only use k to refer to k-mers.\n",
    "\n",
    "This notebook will:\n",
    "\n",
    "* demonstrate set similarity\n",
    " - Jaccard similarity\n",
    "\n",
    "* create the MinHash data structure using a Python library\n",
    " - add elements from a set\n",
    " - output a MinHash sketch of the set\n",
    " - estimate Jaccard similarity between sets using MinHash sketches\n",
    " \n",
    "* adjust underlying properties of the MinHash data structure\n",
    " - vary the base hash function and the number of set permutations\n",
    " - compare Jaccard estimates from MinHash to the actual Jaccard set similarity\n",
    "\n",
    "To begin, import the Python modules we will need:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from matplotlib_venn import venn2, venn3\n",
    "from datasketch import MinHash # a MinHash library by Eric Zhu"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We'll also import the functions that we covered in the [Background](Background.ipynb) workbook. This will make things like decomposing sequences to k-mers easier."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import nbimporter\n",
    "import Background as utils"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "***\n",
    "\n",
    "## A scenario\n",
    "\n",
    "To demonstrate set similarity and MinHash for genomic data, let's use a simple scenario. We'll revisit this scenario in the other workbooks, and the [example workflows](r4.1.Workflows-for-genomics.ipynb) are based on this too.\n",
    "\n",
    "### Antimicrobial Resistance\n",
    "\n",
    "**Antimicrobial resistance** (AMR) is a major public health concern and a massive research area, in which genomics plays a big role.\n",
    "\n",
    "The term AMR refers to the ability of a microbe to resist the effects of antimicrobial drugs, which would usually kill or inactivate the microbe. Within AMR, there is antibiotic resistance, which specifically refers to bacteria that are resistant to antibiotic drugs. These resistant bacteria are harder to treat, needing higher doses or different antibiotics. Consequently, treatment regimes can be longer, more expensive or ineffective.\n",
    "\n",
    "We will focus on antibiotic resistance and how we can use genomics and data sketching to research and combat this problem.\n",
    "\n",
    "\n",
    "### Carbapenem resistance\n",
    "\n",
    "Carbapenems are a group of antibiotic drug, belonging to the **β-lactam** class of antibiotics. They are a very effective antibiotic and are usually reserved for treating multidrug-resistant bacterial infections.\n",
    "\n",
    "Resistance to carbapenems is particularly worrying as there are few drugs that can treat infections caused by carbapenem resistant bacteria.\n",
    "\n",
    "Carbapenem resistance is mediated by a particular type of enzyme that is called a **β-lactamase**. This type of enzyme can break down the structure of the antibiotic, inactivating it.\n",
    "\n",
    "Carbapenems were typically stable against β-lactamases but this is no longer the case, with several classes of β-lactamases emerging that can break down carbapenems - these are called **carbapenemases**.\n",
    "\n",
    "\n",
    "### Carbapenemases\n",
    "\n",
    "Hopefully you are still with me - one more definition to go.... Within the carbapenemase group of enzymes, we have several classes (groan). We won't go into this really, but the important thing is that one of these classes contains **IMP-type carbapenemases** (metallo-β-lactamases).\n",
    "\n",
    "These carbapenemase enzymes are encoded by IMP genes and, importantly, are **plasmid-mediated**. So they can be transferred between bacteria, resulting in the spread of AMR.\n",
    "\n",
    "### TL;DR\n",
    "\n",
    "> IMP genes confer resistance to antibiotics that are reserved for treating multi-drug resistant infections. IMP genes are mobile and contibute to the spread of antibiotic resistance.\n",
    "\n",
    "Let's use set similarity and data sketching to look at IMP genes and see if we can differentiate between IMP genes and non-carbapenemase β-lactamases.\n",
    "\n",
    "\n",
    "## Set similarity\n",
    "\n",
    "Okay. Let's start by taking two IMP genes and seeing how similar they are, based on their k-mer composition."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* get the k-mers (k=7) for each gene:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# These AMR genes are downloaded from the ResFinder database (accessed 05/2017):\n",
    "## blaIMP-1_1_DQ522237\n",
    "geneA = \"ATGAGCAAGTTATCTGTATTCTTTATATTTTTGTTTTGCAGCATTGCTACCGCAGCAGAGTCTTTGCCAGATTTAAAAATTGAAAAGCTTGATGAAGGCGTTTATGTTCATACTTCGTTTAAAGAAGTTAACGGGTGGGGCGTTGTTCCTAAACATGGTTTGGTGGTTCTTGTAAATGCTGAGGCTTACCTAATTGACACTCCATTTACGGCTAAAGATACTGAAAAGTTAGTCACTTGGTTTGTGGAGCGTGGCTATAAAATAAAAGGCAGCATTTCCTCTCATTTTCATAGCGACAGCACGGGCGGAATAGAGTGGCTTAATTCTCGATCTATCCCCACGTATGCATCTGAATTAACAAATGAACTGCTTAAAAAAGACGGTAAGGTTCAAGCCACAAATTCATTTAGCGGAGTTAACTATTGGCTAGTTAAAAATAAAATTGAAGTTTTTTATCCAGGCCCGGGACACACTCCAGATAACGTAGTGGTTTGGTTGCCTGAAAGGAAAATATTATTCGGTGGTTGTTTTATTAAACCGTACGGTTTAGGCAATTTGGGTGACGCAAATATAGAAGCTTGGCCAAAGTCCGCCAAATTATTAAAGTCCAAATATGGTAAGGCAAAACTGGTTGTTCCAAGTCACAGTGAAGTTGGAGACGCATCACTCTTGAAACTTACATTAGAGCAGGCGGTTAAAGGGTTAAACGAAAGTAAAAAACCATCAAAACCAAGCAACTAA\"\n",
    "\n",
    "## blaIMP-2_1_AJ243491\n",
    "geneB = \"ATGAAGAAATTATTTGTTTTATGTGTATGCTTCCTTTGTAGCATTACTGCCGCGGGAGCGCGTTTGCCTGATTTAAAAATCGAGAAGCTTGAAGAAGGTGTTTATGTTCATACATCGTTCGAAGAAGTTAACGGTTGGGGTGTTGTTTCTAAACACGGTTTGGTGGTTCTTGTAAACACTGACGCCTATCTGATTGACACTCCATTTACTGCTACAGATACTGAAAAGTTAGTCAATTGGTTTGTGGAGCGCGGCTATAAAATCAAAGGCACTATTTCCTCACATTTCCATAGCGACAGCACAGGGGGAATAGAGTGGCTTAATTCTCAATCTATTCCCACGTATGCATCTGAATTAACAAATGAACTTCTTAAAAAAGACGGTAAGGTGCAAGCTAAAAACTCATTTAGCGGAGTTAGTTATTGGCTAGTTAAAAATAAAATTGAAGTTTTTTATCCCGGCCCGGGGCACACTCAAGATAACGTAGTGGTTTGGTTACCTGAAAAGAAAATTTTATTCGGTGGTTGTTTTGTTAAACCGGACGGTCTTGGTAATTTGGGTGACGCAAATTTAGAAGCTTGGCCAAAGTCCGCCAAAATATTAATGTCTAAATATGTTAAAGCAAAACTGGTTGTTTCAAGTCATAGTGAAATTGGGGACGCATCACTCTTGAAACGTACATGGGAACAGGCTGTTAAAGGGCTAAATGAAAGTAAAAAACCATCACAGCCAAGTAACTAA\"\n",
    "\n",
    "# get the canonical k-mers for each gene\n",
    "kmersA = utils.getKmers(geneA, 7)\n",
    "kmersB = utils.getKmers(geneB, 7)\n",
    "\n",
    "# sanity check and just print a few k-mers out from each gene\n",
    "print(kmersA[0:10])\n",
    "print(kmersB[0:10])\n",
    "\n",
    "# how many k-mers do we have\n",
    "print(\"blaIMP-1 is {} bases long and has {} k-mers (k=7)\" .format(len(geneA), len(kmersA)))\n",
    "print(\"blaIMP-2 is {} bases long and has {} k-mers (k=7)\" .format(len(geneA), len(kmersB)))\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* now we have the two genes decomposed to sets of k-mers, let's see how similar they are\n",
    "\n",
    "* to do this, we will calculate the Jaccard similarity:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# this step converts the Python lists of k-mers into sets, which isn't important for the example but just let's us use some built in functions\n",
    "setA = set(kmersA)\n",
    "setB = set(kmersB)\n",
    "\n",
    "# get the intersection and union of the k-mer sets\n",
    "intersection = len(setA.intersection(setB))\n",
    "union = len(setA.union(setB))\n",
    "\n",
    "# calculate the Jaccard similarity\n",
    "actualJaccardSimilarity = float(intersection)/float(union)\n",
    "\n",
    "# print the values\n",
    "print(\"The intersection of the k-mer sets = {}\\nThe union of the k-mer sets = {}\\n\" .format(intersection, union) )\n",
    "print(\"The Jaccard simularity between genes A and B = {}\\n\" .format(actualJaccardSimilarity))\n",
    "\n",
    "# visualise with a Venn diagram\n",
    "print(\"Venn diagram showing the intersection of the carbapenemase genes:\")\n",
    "venn2([setA, setB], set_labels = ('blaIMP-1', 'blaIMP-2'))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* great, so now we know that are two IMP genes share 299 k-mers and the Jaccard similarity is 0.29\n",
    "\n",
    "* now, how does a non-carbapenemase β-lactamase gene compare?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# CTX-M β-lactamases are also plasmid mediated but they can't break down carbapenemases\n",
    "## blaCTX-M-1_6_DQ915955\n",
    "geneC = \"ATGGTTAAAAAATCACTGCGTCAGTTCACGCTGATGGCGACGGCAACCGTCACGCTGTTGTTAGGAAGTGTGCCGCTGTATGCGCAAACGGCGGACGTACAGCAAAAACTTGCCGAATTAGAGCGGCAGTCGGGAGGAAGACTGGGTGTGGCATTGATTAACACAGCAGATAATTCGCAAATACTTTATCGTGCTGATGAGCGCTTTGCGATGTGCAGCACCAGTAAAGTGATGGCCGTGGCCGCGGTGCTGAAGAAAAGTGAAAGCGAACCGAATCTGTTAAATCAGCGAGTTGAGATCAAAAAATCTGACTTGGTTAACTATAATCCGATTGCGGAAAAGCACGTCGATGGGACGATGTCACTGGCTGAGCTTAGCGCGGCCGCGCTACAGTACAGCGATAACGTGGCGATGAATAAGCTGATTTCTCACGTTGGCGGCCCGGCTAGCGTCACCGCGTTCGCCCGACAGCTGGGAGACGAAACGTTCCGTCTCGACCGTACCGAGCCGACGTTAAACACCGCCATTCCGGGCGATCCGCGTGATACCACTTCACCTCGGGCAATGGCGCAAACTCTGCGTAATCTGACGCTGGGTAAAGCATTGGGTGACAGCCAACGGGCGCAGCTGGTGACATGGATGAAAGGCAATACCACCGGTGCAGCGAGCATTCAGGCTGGACTGCCTGCTTCCTGGGTTGTGGGGGATAAAACCGGCAGCGGTGACTATGGCACCACCAACGATATCGCGGTGATCTGGCCAAAAGATCGTGCGCCGCTGATTCTGGTCACTTACTTCACCCAGCCTCAACCTAAGGCAGAAAGCCGTCGCGATGTATTAGCGTCGGCGGCTAAAATCGTCACCAACGGTTTGTAA\"\n",
    "kmersC = utils.getKmers(geneC, 7)\n",
    "setC = set(kmersC)\n",
    "# visualise with a Venn diagram\n",
    "print(\"Venn diagram showing the intersection of the β-lactamase genes:\")\n",
    "venn3([setA, setB, setC], set_labels = ('blaIMP-1', 'blaIMP-2', 'blaCTX-M-1'))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* we can see that the blaIMP genes share a lot more k-mers than either do with the blaCTX-M gene, which was to be expected!\n",
    "\n",
    "* now let's try approximating these similarities using MinHash"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Estimating set similarity with MinHash\n",
    "\n",
    "We'll use a MinHash library to save us coding up the data structure. This library means we don't have to define our own hash functions (you can do if you want though). All we need to do is specify how many set permutations we want to sample from and then we just provide it the k-mers.\n",
    "\n",
    "* create the MinHash data structure: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# define the number of set permutations to apply to our k-mer sets\n",
    "permutations=100\n",
    "\n",
    "# construct a MinHash data structure for each bla-IMP gene\n",
    "mh1 = MinHash(num_perm=permutations)\n",
    "mh2 = MinHash(num_perm=permutations)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* now we add the k-mers from each gene to their respective MinHash:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# define a helper function to add set elements to a MinHash data structure\n",
    "def add2MinHash(mh, data):\n",
    "    # loop through the list of k-mers\n",
    "    for d in data:\n",
    "        # here, we use the library's update function\n",
    "        # this is where each k-mer is hashed, with a different hash for each set permutation\n",
    "        # if the k-mer hash is the lowest value seen to date, it is added to the MinHash sketch\n",
    "        mh.update(d.encode('utf8'))\n",
    "    return mh\n",
    "\n",
    "# add the k-mers\n",
    "mh1 = add2MinHash(mh1, kmersA)\n",
    "mh2 = add2MinHash(mh2, kmersB)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* now we have MinHash sketches, one for each gene\n",
    "\n",
    "* to compare the sketches and get a Jaccard similarity estimate, we just use the jaccard method from the library!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"Estimated Jaccard similarity for blaIMP-1 and blaIMP-2 is\", mh1.jaccard(mh2))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* not bad, our actual Jaccard similarity between blaIMP-1 and blaIMP-2 was 0.29 and our estimate is 0.24\n",
    "* this estimate is based on only 100 k-mers per sketch, compared to the 735 k-mers per gene for the actual similarity\n",
    "* for the sake of completeness, let's sketch the blaCTX-M gene"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mh3 = MinHash(num_perm=permutations)\n",
    "mh3 = add2MinHash(mh3, kmersC)\n",
    "print(\"Estimated Jaccard similarity for blaIMP-1 and blaCTX-M-1 is\", mh1.jaccard(mh3))\n",
    "print(\"Estimated Jaccard similarity for blaIMP-2 and blaCTX-M-1 is\", mh2.jaccard(mh3))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## MinHash a genomic data stream\n",
    "\n",
    "So MinHash is pretty good at approximating similarity between genes. We haven't really seen the benefits though, as all we did was essential reduce the number of k-mers from 735 to 100, which doesn't save us much in terms of time or space.\n",
    "\n",
    "Let's finish up this notebook by showing how we can apply MinHash to something more useful.\n",
    "\n",
    "A bacteria that we look at in the example workflow is [Enterobacter cloacae](https://en.wikipedia.org/wiki/Enterobacter_cloacae), which is a nosocomial pathogen. Let's use MinHash to compare how similar the reference genome is for this bacteria to another species from the same genus, [Enterobacter mori](https://microbe.directory/species/Enterobacter+mori), which is a plant pathogen.\n",
    "\n",
    "\n",
    "* let's simulate a genomic data stream from an E.cloacae download and MinHash it:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# stream k-mers from E. cloacae\n",
    "dataStream = utils.streamGenome('../data/GCF_000025565.1_ASM2556v1_genomic.fna.gz', 7)\n",
    "\n",
    "# set up the MinHash\n",
    "cloacaeSketch = MinHash(num_perm=128)\n",
    "\n",
    "# sketch! this will take a couple of minutes\n",
    "for kmer in dataStream:\n",
    "    cloacaeSketch.update(kmer.encode('utf8'))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* now do the same for E.mori:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# stream k-mers from E. cloacae\n",
    "dataStream = utils.streamGenome('../data/GCF_000211415.1_ASM21141v1_genomic.fna.gz', 7)\n",
    "\n",
    "# set up the MinHash\n",
    "moriSketch = MinHash(num_perm=128)\n",
    "for kmer in dataStream:\n",
    "    moriSketch.update(kmer.encode('utf8'))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* and compare!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"Estimated Jaccard similarity for E.cloacae and E.mori reference genomes is\", cloacaeSketch.jaccard(moriSketch))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So, these genomes are estimated to be 93% similar, based on k-mer composition. We have used only 128 hashes to approximate each genome, and each hash is derived from a k-mer which is 7 bases long.\n",
    "\n",
    "A good exercise is to try the above again, this time increasing the size of k. Hopefully you should see that the longer the k, the less similar the genomes are. This is because the longer the k-mer, the more species specific it is. So, increase the k-mer length to 31 here and you should be able to differentiate between these two enterobacter species quite nicely.\n",
    "\n",
    "You can also try changing the number of values in the sketches - long sketches give better estimates!\n",
    "\n",
    "We'll revisit MinHash later, this time using some cool bioinformatics software which can do the above more quickly and output more information etc."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "***\n",
    "\n",
    "Return to [this notebook](./notebooks/Background.ipynb) to access the table of contents and select another notebook.\n",
    "\n",
    "***"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "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
}