notebooks/r4.4.Data-mining.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Data mining\n",
    "\n",
    "Now that we have re-analysed the outbreak data from the [paper](https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4001082), we may wish to augment that original study with new data or ask some new questions.\n",
    "\n",
    "For this workflow, we'll be sketching with:\n",
    "\n",
    "* Bloom filters\n",
    "* HyperLogLog (coming soon)\n",
    "\n",
    "***\n",
    "\n",
    "## Isolates with matching resistome profiles\n",
    "\n",
    "In our [resistome profiling workflow](./r4.3.Resistome-profiling.ipynb), we did a MinHash screen of raw reads and then aligned to AMR gene variation graphs. We ended up with the following resistome profile for isolate EC1a (ERX168346):\n",
    "\n",
    "```\n",
    "(Bla)SHV-12\n",
    "(Sul)SulII\n",
    "(Tmt)DfrA1\n",
    "(Bla)SHV-183\n",
    "(Sul)SulI\n",
    "(AGly)Sat-2A\n",
    "```\n",
    "\n",
    "> If we looked at the alignments against the variation graphs, we would see that the blaSHV-183 is likely a mis-call of blaSHV-12 as they are very similar sequences but blaSHV-12 has more support. We will remove blaSHV-183 from the profile.\n",
    "\n",
    "Now, what if we wanted to find other isolates with a matching resistome profile? Trawling through the [ENA](https://www.ebi.ac.uk/ena) to download the sequence data, re-run resistome profiling experiments etc. will take a long time. Instead, let's use the brilliant [BIGSI](http://www.bigsi.io/) and its index of the ENA.\n",
    "\n",
    "* start with getting the sequences for all the AMR genes that are in our isolate's resistome profile:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# download the AMR gene database we used for resistome profiling\n",
    "!wget https://raw.githubusercontent.com/will-rowe/groot/master/db/full-ARG-databases/arg-annot-db/argannot-args.fna\n",
    "\n",
    "# load the FASTA file into a database\n",
    "import screed\n",
    "screed.make_db(\"./argannot-args.fna\")\n",
    "fadb = screed.ScreedDB(\"./argannot-args.fna\")\n",
    "\n",
    "# store each gene\n",
    "SHV12 = fadb[\"argannot~~~(Bla)SHV-12~~~FJ685654:24-860\"].sequence\n",
    "SulII = fadb[\"argannot~~~(Sul)SulII~~~EU360945:1617-2432\"].sequence\n",
    "DfrA1 = fadb[\"argannot~~~(Tmt)DfrA1~~~JQ794607:474\"].sequence\n",
    "SulI = fadb[\"argannot~~~(Sul)SulI~~~AF071413:6700-7539\"].sequence\n",
    "Sat2A = fadb[\"argannot~~~(AGly)Sat-2A~~~X51546:518-1042\"].sequence\n",
    "\n",
    "# create the resistome profile\n",
    "resistomeProfile = {\"SHV12\":SHV12, \"SulII\":SulII, \"DfrA1\":DfrA1, \"SulI\":SulI, \"Sat2A\":Sat2A}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* now, set up the search function for our calls to the BIGSI ENA index:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import requests\n",
    "def search(seq, threshold):\n",
    "    # set up the api call\n",
    "    url = \"http://api.bigsi.io/search?threshold=%f&seq=%s\" % (float(threshold), seq)\n",
    "    results = requests.get(url).json()\n",
    "    samples = []\n",
    "    # bigsi includes species metadata, derived from Bracken + Kraken analysis of the ENA data\n",
    "    classification = []\n",
    "    for i, j in list(results.values())[0][\"results\"].items():\n",
    "        samples.append(i)\n",
    "        classification.append(j)\n",
    "    return samples, classification\n",
    "\n",
    "# set threshold for proportion of query k-mers contained in results\n",
    "threshold=0.8"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* start with a sanity check and make sure our isolate is returned from the ENA:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# BIGSI uses ERR, so look ours up (see r.4.2.Sample-QC for the table)\n",
    "isolateID = \"ERR193657\"\n",
    "\n",
    "# to save us re-running the BIGSI search (not that it took long), store the results for later\n",
    "results = []\n",
    "classificationDict = {}\n",
    "\n",
    "# search BIGSI for each gene in our resistome profile\n",
    "for gene in resistomeProfile:\n",
    "    print(\"query: {}\" .format(gene))\n",
    "    \n",
    "    # run the search\n",
    "    genomes, classification = search(resistomeProfile[gene], threshold)\n",
    "    \n",
    "    # check the results\n",
    "    found = False\n",
    "    for i, ERR in enumerate(genomes):\n",
    "        classificationDict[ERR] = classification[i][\"species\"]\n",
    "        if (ERR == isolateID):\n",
    "            found = True\n",
    "            results.append(genomes)\n",
    "    \n",
    "    # check if our isolate was in the search results for this gene\n",
    "    if (found == True):\n",
    "        print(\"\\t- our isolate was returned in the results\")\n",
    "    else:\n",
    "        print(\"\\t- our isolate was NOT returned in the results\")\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* great - our isolate was returned for each gene in it's resistome profile. Now let's do something more interesting. Using the search results, find all the isolates in BIGSI that share the resistome profile of our isolate:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# use a set data type\n",
    "matchingIsolates = set(results[0])\n",
    "\n",
    "# get intersection on search results for each gene\n",
    "for isolates in results[1:]:\n",
    "    matchingIsolates.intersection_update(isolates)\n",
    "\n",
    "# another sanity check, make sure our isolate is there\n",
    "check = False\n",
    "for isolate in matchingIsolates:\n",
    "    if isolate == isolateID:\n",
    "        check = True\n",
    "if (check == False):\n",
    "    print(\"fail. our isolate was not returned in the BIGSI search for the resistome profile\")\n",
    "elif (len(matchingIsolates) < 1):\n",
    "    print(\"fail. no isolates matching our resistome profile\")\n",
    "else:\n",
    "    \n",
    "    # print the number of isolates we have that match our resistome profile\n",
    "    print(\"success: {} isolates in ENA matching the resistome profile\" .format(len(matchingIsolates)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* we stored the classifications that BIGSI assigned to each entry in it's index (using kraken/braken), now we can check what we have:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for isolate in matchingIsolates:\n",
    "    print(\"isolate:\")\n",
    "    print(\"\\tERR= {}\" .format(isolate))\n",
    "    print(\"\\tpred.= {}\" .format(classificationDict[isolate]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* there are few *e.cloacae* isolates in there, let's filter them out:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cloacaeIsolates=[]\n",
    "import re\n",
    "for isolate in matchingIsolates:\n",
    "    search = re.search( r'Enterobacter\\scloacae', classificationDict[isolate], re.M|re.I)\n",
    "    if search:\n",
    "        cloacaeIsolates.append(isolate)\n",
    "\n",
    "print(\"success: {} e.cloacae isolates matching the resistome profile\" .format(len(cloacaeIsolates)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* we can try getting some more information on them:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# code adapted from: https://bioinfo.umassmed.edu/bootstrappers/guides/main/python_get_sra_run_ids.html\n",
    "import requests, csv, io\n",
    "def getInfoTableFromSearchTerm(search):\n",
    "        payload = {\"save\": \"efetch\",\"db\": \"sra\",\"rettype\" : \"runinfo\", \"term\" : search };\n",
    "        r = requests.get('http://trace.ncbi.nlm.nih.gov/Traces/sra/sra.cgi', params=payload)\n",
    "        if 200 ==  r.status_code:\n",
    "            if r.text.isspace():\n",
    "                raise Exception(\"Got blank string from \" + str(r.url ))\n",
    "            else:\n",
    "                reader_list = csv.DictReader(io.StringIO(r.text))\n",
    "                infoRows = []\n",
    "                for row in reader_list:\n",
    "                    infoRows.append(row)\n",
    "                if 0 == len(infoRows):\n",
    "                    raise Exception('Found %d entries in SRA for \"%s\" when expecting at least 1' % (len(infoRows), search))\n",
    "                else:        \n",
    "                    return infoRows\n",
    "                return infoRows\n",
    "        else:\n",
    "            raise Exception(\"Error in downloading from \" + str(r.url) + \" got response code \" + str(r.status_code))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for isolate in cloacaeIsolates:\n",
    "    print(isolate)\n",
    "    \n",
    "    # lookup the ENA info\n",
    "    infoTable = getInfoTableFromSearchTerm(isolate)\n",
    "    if (len(infoTable) == 0):\n",
    "        print(\"\\tENA lookup failed - no hits\")\n",
    "\n",
    "    # print some bits of info\n",
    "    for row in infoTable:\n",
    "        print(\"\\tUpload Date:\\t\" + row.get(\"LoadDate\"))\n",
    "        print(\"\\tPlatform:\\t\" + row.get(\"Platform\"))\n",
    "        print(\"\\tBioproject:\\t\" + row.get(\"BioProject\"))\n",
    "        print(\"\\tSubmitting Centre:\\t\" + row.get(\"CenterName\"))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    ">More content coming soon"
   ]
  },
  {
   "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
}