notebooks/r3.2.Set-membership-with-Bloom-filters.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 3.2. Set membership with Bloom filters"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Set membership"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import nbimporter\n",
    "import Background as utils\n",
    "\n",
    "# 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",
    "# 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"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* create a Bloom filter:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from bloom_filter import BloomFilter\n",
    "sketch = BloomFilter(max_elements=1000, error_rate=0.1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* add one element and test the sketch:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# get a single k-mer to add\n",
    "testKmer=kmersA[0]\n",
    "\n",
    "# check it is not already in the sketch\n",
    "if testKmer in sketch:\n",
    "    print(\"error: sketch should be empty\")\n",
    "\n",
    "# add the k-mer\n",
    "sketch.add(testKmer)\n",
    "\n",
    "# check the sketch to see if the k-mer was added\n",
    "if testKmer in sketch:\n",
    "    print(\"k-mer {} is in the Bloom filter\" .format(testKmer))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* now, add all the k-mers from geneA to the sketch:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for kmer in kmersA:\n",
    "    sketch.add(kmer)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* use the Bloom filter to see approximately how many k-mers from geneB are present in geneA:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "count = 0\n",
    "for kmer in kmersB:\n",
    "    if kmer in sketch:\n",
    "        count += 1\n",
    "\n",
    "print(\"approximately {} k-mers are shared between gene A and gene B\" .format(count))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* to get the actual number of shared k-mers, we can get the intersection of kmersA and kmersB by looping through one set and checking it against the other."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "count = 0\n",
    "for kmer in kmersA:\n",
    "    if kmer in kmersB:\n",
    "        count += 1\n",
    "\n",
    "print( \"actual number of k-mers shared between gene A and gene B is: {}\" .format(count))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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
}