interpret-sra-live.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "487ad35c",
   "metadata": {},
   "source": [
    "# Search ~500,000 metagenomes for matches to a genome sequence\n",
    "\n",
    "See [the GitHub repo](https://github.com/sourmash-bio/2022-search-sra-with-mastiff) for more info!\n",
    "\n",
    "The mastiff server is unpublished work by Luiz Irber - please [file an issue](https://github.com/dib-lab/sourmash/issues) if you are interested in how to cite it."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "342d1324",
   "metadata": {},
   "source": [
    "## Python miscellany\n",
    "\n",
    "Import all the Python libraries we need to use"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "0b6af9cd",
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import os\n",
    "import sourmash, screed\n",
    "import io\n",
    "import urllib3"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bd2baaa3",
   "metadata": {},
   "source": [
    "## Download SRA metadata information.\n",
    "\n",
    "This is a precomputed CSV that we download once per instance; it's used to summarize the `ScientificName` (i.e. metagenome type) of matches at the bottom."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "fe12875c",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "File ‘sra.runinfo.csv.gz’ already there; not retrieving.\r\n"
     ]
    }
   ],
   "source": [
    "!wget --no-clobber https://osf.io/download/762mk/ -O sra.runinfo.csv.gz"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "db3e871c",
   "metadata": {},
   "source": [
    "## Load the SRA metadata into a pandas dataframe"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "7a72d081",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Loaded 702013 records from sra.runinfo.csv.gz\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>Unnamed: 0</th>\n",
       "      <th>Run</th>\n",
       "      <th>ScientificName</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>0</td>\n",
       "      <td>SRR18036904</td>\n",
       "      <td>bovine metagenome</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>1</td>\n",
       "      <td>SRR18036905</td>\n",
       "      <td>bovine metagenome</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>2</td>\n",
       "      <td>SRR18036906</td>\n",
       "      <td>bovine metagenome</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>3</td>\n",
       "      <td>SRR18036907</td>\n",
       "      <td>bovine metagenome</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>4</td>\n",
       "      <td>SRR18036908</td>\n",
       "      <td>bovine metagenome</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   Unnamed: 0          Run     ScientificName\n",
       "0           0  SRR18036904  bovine metagenome\n",
       "1           1  SRR18036905  bovine metagenome\n",
       "2           2  SRR18036906  bovine metagenome\n",
       "3           3  SRR18036907  bovine metagenome\n",
       "4           4  SRR18036908  bovine metagenome"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "run_info = pd.read_csv('sra.runinfo.csv.gz')\n",
    "print(f\"Loaded {len(run_info)} records from sra.runinfo.csv.gz\")\n",
    "run_info.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "25a43468",
   "metadata": {},
   "source": [
    "## Generate a sourmash sketch for a query\n",
    "\n",
    "In order to query mastiff with a genome, we need to sketch the genome using [sourmash](https://sourmash.readthedocs.io/). Here use the Python interface for sourmash.\n",
    "\n",
    "The below code loads all of the sequence data using [the screed library](https://screed.readthedocs.io/), and feeds it into a sourmash [FracMinHash sketch](https://dib-lab.github.io/2020-paper-sourmash-gather/).\n",
    "\n",
    "You can replace the `QUERY_SEQUENCE_FILE` with your own query filename, if you like. It should be over 10kb in size (ideally you want ~5-10 hashes to result in order for the query to be robust)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "bf92950e",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "generated 5307 hashes by sketching 5.31291e+06 bp from 'sequences/shewanella.fa.gz'\n"
     ]
    }
   ],
   "source": [
    "# generate sourmash sketch\n",
    "QUERY_SEQUENCE_FILE='sequences/shewanella.fa.gz'\n",
    "\n",
    "total_bp = 0\n",
    "\n",
    "sketch = sourmash.MinHash(0, 21, scaled=1000)\n",
    "with screed.open(QUERY_SEQUENCE_FILE) as records:\n",
    "    for record in records:\n",
    "        sketch.add_sequence(record.sequence, force=True)\n",
    "        total_bp += len(record.sequence)\n",
    "        \n",
    "print(f\"generated {len(sketch)} hashes by sketching {total_bp:g} bp from '{QUERY_SEQUENCE_FILE}'\")\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "61f55d76",
   "metadata": {},
   "source": [
    "## Serialize the sourmash sketch\n",
    "\n",
    "Now, we convert the sourmash sketch into gzipped JSON bytes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "24eb711e",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "serialized sourmash signature into 44044 bytes.\n"
     ]
    }
   ],
   "source": [
    "# serialize to SourmashSignature / gzipped JSON\n",
    "ss = sourmash.SourmashSignature(sketch)\n",
    "\n",
    "buf = io.BytesIO()\n",
    "sourmash.save_signatures([ss], buf, compression=True)\n",
    "print(f\"serialized sourmash signature into {len(buf.getvalue())} bytes.\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1cd2c23a",
   "metadata": {},
   "source": [
    "## Query the mastiff server\n",
    "\n",
    "Next, use `urllib3` to contact the mastiff server, send in the sourmash signature information, and get the results."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "47da93e9",
   "metadata": {},
   "outputs": [],
   "source": [
    "# query mastiff\n",
    "http = urllib3.PoolManager()\n",
    "r = http.request('POST',\n",
    "                 'https://mastiff.sourmash.bio/search',\n",
    "                 body = buf.getvalue(),\n",
    "                 headers = { 'Content-Type': 'application/json'})\n",
    "\n",
    "query_results_text = r.data.decode('utf-8')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "255d2ebd",
   "metadata": {},
   "source": [
    "## Parse the mastiff results as CSV\n",
    "\n",
    "mastiff returns the results in CSV format, with two columns - the SRA accession, and the containment index of the query.\n",
    "\n",
    "The containment index is calculated as the number of hashes in the query that are present in a matching metagenome, divided by the total number of hashes in the query. A containment of 1.0 means every hash in the query sketch is present, and estimates that the original query sequence is entirely present; a containment of 0.5 means that 50% of the query sketch is present, suggesting that only about half of the query sequence is actually in the metagenome."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "e4912aae",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Loaded 5702 mastiff results into a dataframe!\n"
     ]
    }
   ],
   "source": [
    "# load results into a pandas data frame\n",
    "results_wrap_fp = io.StringIO(query_results_text)\n",
    "mastiff0_df = pd.read_csv(results_wrap_fp)\n",
    "\n",
    "print(f\"Loaded {len(mastiff0_df)} mastiff results into a dataframe!\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "414453dc",
   "metadata": {},
   "source": [
    "## Explore the results\n",
    "\n",
    "Look at the best and worst matches."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "f30a3b7a",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>SRA accession</th>\n",
       "      <th>containment</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>ERR1254310</td>\n",
       "      <td>1.000000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>ERR1254284</td>\n",
       "      <td>1.000000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>ERR1254285</td>\n",
       "      <td>0.999812</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>ERR1254311</td>\n",
       "      <td>0.999623</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>SRR606249</td>\n",
       "      <td>0.999246</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "  SRA accession  containment\n",
       "0    ERR1254310     1.000000\n",
       "1    ERR1254284     1.000000\n",
       "2    ERR1254285     0.999812\n",
       "3    ERR1254311     0.999623\n",
       "4     SRR606249     0.999246"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "mastiff0_df.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "ea3d2cbc",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>SRA accession</th>\n",
       "      <th>containment</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>5697</th>\n",
       "      <td>SRR5918872</td>\n",
       "      <td>0.009422</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5698</th>\n",
       "      <td>SRR5207236</td>\n",
       "      <td>0.009422</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5699</th>\n",
       "      <td>SRR5499129</td>\n",
       "      <td>0.009422</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5700</th>\n",
       "      <td>SRR8500488</td>\n",
       "      <td>0.009422</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5701</th>\n",
       "      <td>ERR599067</td>\n",
       "      <td>0.009422</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "     SRA accession  containment\n",
       "5697    SRR5918872     0.009422\n",
       "5698    SRR5207236     0.009422\n",
       "5699    SRR5499129     0.009422\n",
       "5700    SRR8500488     0.009422\n",
       "5701     ERR599067     0.009422"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "mastiff0_df.tail()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "be60b056",
   "metadata": {},
   "source": [
    "## Discard matches below a containment of 0.2\n",
    "\n",
    "Per [our work on containment ANI](https://github.com/sourmash-bio/sourmash/issues/1859), this corresponds to a minimum of 92% average nucleotide identity."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "36ba51c3",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Filtering results at a containment of >= 0.20; 312 of 5702 left.\n"
     ]
    }
   ],
   "source": [
    "THRESHOLD=0.2\n",
    "\n",
    "mastiff_df = mastiff0_df[mastiff0_df['containment'] >= THRESHOLD]\n",
    "print(f\"Filtering results at a containment of >= {THRESHOLD:.2f}; {len(mastiff_df)} of {len(mastiff0_df)} left.\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b480d8b5",
   "metadata": {},
   "source": [
    "## Explore metagenome \"types\" through `Scientific Name`\n",
    "\n",
    "Now, correlate results with the `ScientificName` metadata, which records which \"kind\" of metagenome this is."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "724bcf39",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>containment</th>\n",
       "      <th>ScientificName</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>SRA accession</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>ERR1254310</th>\n",
       "      <td>1.000000</td>\n",
       "      <td>metagenome</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>ERR1254284</th>\n",
       "      <td>1.000000</td>\n",
       "      <td>metagenome</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>ERR1254285</th>\n",
       "      <td>0.999812</td>\n",
       "      <td>metagenome</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>ERR1254311</th>\n",
       "      <td>0.999623</td>\n",
       "      <td>metagenome</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>SRR606249</th>\n",
       "      <td>0.999246</td>\n",
       "      <td>synthetic metagenome</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "               containment        ScientificName\n",
       "SRA accession                                   \n",
       "ERR1254310        1.000000            metagenome\n",
       "ERR1254284        1.000000            metagenome\n",
       "ERR1254285        0.999812            metagenome\n",
       "ERR1254311        0.999623            metagenome\n",
       "SRR606249         0.999246  synthetic metagenome"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "mastiff2_df = mastiff_df.set_index('SRA accession').join(run_info.set_index('Run')['ScientificName'])\n",
    "mastiff2_df.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "f1e11b86",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "9\n"
     ]
    }
   ],
   "source": [
    "# how many have 'null' scientific name?\n",
    "null_df = mastiff2_df[mastiff2_df['ScientificName'].isnull()]\n",
    "print(len(null_df))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "53a0deeb",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Of 312 MAGsearch results, 303 have non-null ScientificName\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>containment</th>\n",
       "      <th>ScientificName</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>SRA accession</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>ERR1254310</th>\n",
       "      <td>1.000000</td>\n",
       "      <td>metagenome</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>ERR1254284</th>\n",
       "      <td>1.000000</td>\n",
       "      <td>metagenome</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>ERR1254285</th>\n",
       "      <td>0.999812</td>\n",
       "      <td>metagenome</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>ERR1254311</th>\n",
       "      <td>0.999623</td>\n",
       "      <td>metagenome</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>SRR606249</th>\n",
       "      <td>0.999246</td>\n",
       "      <td>synthetic metagenome</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "               containment        ScientificName\n",
       "SRA accession                                   \n",
       "ERR1254310        1.000000            metagenome\n",
       "ERR1254284        1.000000            metagenome\n",
       "ERR1254285        0.999812            metagenome\n",
       "ERR1254311        0.999623            metagenome\n",
       "SRR606249         0.999246  synthetic metagenome"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "mastiff3_df = mastiff2_df[~mastiff2_df['ScientificName'].isnull()]\n",
    "print(f\"Of {len(mastiff2_df)} MAGsearch results, {len(mastiff3_df)} have non-null ScientificName\")\n",
    "mastiff3_df.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9803ccc9",
   "metadata": {},
   "source": [
    "## Group and aggregate by metagenome type\n",
    "\n",
    "Explore which types of metagenomes tend to contain this query."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "461ae616",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "metagenome                  111\n",
       "freshwater metagenome        32\n",
       "wastewater metagenome        30\n",
       "oyster metagenome            22\n",
       "aquifer metagenome           17\n",
       "marine metagenome            11\n",
       "food metagenome              11\n",
       "aquatic metagenome            9\n",
       "biofilm metagenome            9\n",
       "seawater metagenome           9\n",
       "groundwater metagenome        6\n",
       "synthetic metagenome          5\n",
       "bioreactor metagenome         3\n",
       "human skin metagenome         3\n",
       "soil metagenome               2\n",
       "money metagenome              2\n",
       "fish metagenome               2\n",
       "fish gut metagenome           2\n",
       "oral metagenome               2\n",
       "aquatic viral metagenome      2\n",
       "Name: ScientificName, dtype: int64"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# what are the top ScientificNames of the matches?\n",
    "mastiff3_df[\"ScientificName\"].value_counts()[:20]"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "py38arm",
   "language": "python",
   "name": "py38arm"
  },
  "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.8.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}