analyze_viral_genome_windows_size.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from Bio import SeqIO\n",
    "import re"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Provides utility functions to find whether words like 'flu' are in viral\n",
    "# genome sequence descriptions, copied from\n",
    "# https://stackoverflow.com/a/5320179. Includes tests.\n",
    "\n",
    "def find_whole_word(w, flags=re.IGNORECASE):\n",
    "    return re.compile(r'\\b({0})\\b'.format(w), flags=flags).search\n",
    "\n",
    "find_flu = find_whole_word('flu')\n",
    "assert find_flu('flu') is not None\n",
    "assert find_flu('aksdjflkaj') is None\n",
    "assert find_flu('flue') is None\n",
    "assert find_flu('flu|1234') is not None\n",
    "assert find_flu('flu-like virus') is not None\n",
    "assert find_flu('Human orthopneumovirus (HRSV) Subgroup B, complete genome') is None\n",
    "\n",
    "find_influenza = find_whole_word('influenza')\n",
    "assert find_influenza('influenza') is not None\n",
    "assert find_influenza('Influenza A virus (A/Puerto Rico/8/1934(H1N1)) segment 1, complete sequence') is not None\n",
    "assert find_influenza('Influenza B virus RNA 1, complete sequence') is not None\n",
    "assert find_influenza('Influenza B virus (B/Lee/1940) segment 2, complete sequence.') is not None\n",
    "assert find_influenza('Influenza C virus (C/Ann Arbor/1/50) PB2 gene for polymerase 2, complete cds.') is not None\n",
    "assert find_influenza('Human parainfluenza virus 1, complete genome') is None\n",
    "assert find_influenza('Human orthopneumovirus (HRSV) Subgroup B, complete genome') is None\n",
    "\n",
    "find_human_immunodeficiency_virus = find_whole_word('human immunodeficiency virus')\n",
    "assert find_human_immunodeficiency_virus('human immunodeficiency virus') is not None\n",
    "assert find_human_immunodeficiency_virus('immunodeficiency virus') is None\n",
    "assert find_human_immunodeficiency_virus('Human immunodeficiency virus 1, complete genome') is not None\n",
    "assert find_human_immunodeficiency_virus('Simian immunodeficiency virus (SIV) - SIVagmMAL genomic RNA, complete genome, strain: ZMB') is None\n",
    "\n",
    "find_hiv = find_whole_word('hiv')\n",
    "assert find_hiv('hiv') is not None\n",
    "assert find_hiv('chives') is None\n",
    "assert find_hiv('hiv-1') is not None\n",
    "assert find_hiv('Human immunodeficiency virus 1 (HIV-1)|1234') is not None\n",
    "assert find_hiv('Simian immunodeficiency virus (SIV) - SIVagmMAL genomic RNA, complete genome, strain: ZMB') is None"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Counts the number of unique {dna_window_size}-mers in all the viral\n",
    "# genomes in the Reference Viral Database, a database attempting to\n",
    "# comprehensively catalog viral genomes:\n",
    "# https://hive.biochemistry.gwu.edu/rvdb.\n",
    "#\n",
    "# Notes:\n",
    "# 1. You can download the dataset \"U-RVDBv15.1.fasta\" from the above\n",
    "# website.\n",
    "# 2. Running this chunk may require a machine with lots of RAM. I used an\n",
    "# AWS EC2 m5a.12xlarge machine with 192 GB RAM.\n",
    "\n",
    "dna_window_size = 42\n",
    "\n",
    "# whether to exclude {dna_window_size}-mers exclusively in flu genomes\n",
    "excl_influenza = True\n",
    "# whether to exclude {dna_window_size}-mers exclusively in HIV genomes\n",
    "excl_hiv = True\n",
    "\n",
    "upper_bound_all_rvdb_windows = set() # stores {dna_window_size}-mers\n",
    "num_seqs_seen = 0\n",
    "seqs = SeqIO.parse(\"data/U-RVDBv15.1.fasta\", \"fasta\")\n",
    "\n",
    "# store counts of excluded viral genomes in each category\n",
    "num_flu = 0\n",
    "num_influenza = 0\n",
    "num_human_immunodeficiency_virus = 0\n",
    "num_hiv = 0\n",
    "\n",
    "for seq in seqs:\n",
    "    # report number of viral genomes screened so far\n",
    "    num_seqs_seen += 1\n",
    "    if num_seqs_seen % 100000 == 0:\n",
    "        print(f'{num_seqs_seen} sequences screened')\n",
    "    # exclude flu and HIV sequences, and count these\n",
    "    excl = False\n",
    "    if excl_influenza:\n",
    "        if find_flu(seq.description) is not None:\n",
    "            num_flu += 1\n",
    "            excl = True\n",
    "        if find_influenza(seq.description) is not None:\n",
    "            num_influenza += 1\n",
    "            excl = True\n",
    "    if excl_hiv:\n",
    "        if find_human_immunodeficiency_virus(seq.description) is not None:\n",
    "            num_human_immunodeficiency_virus += 1\n",
    "            excl = True\n",
    "        if find_hiv(seq.description) is not None:\n",
    "            num_hiv += 1\n",
    "            excl = True\n",
    "    if excl:\n",
    "        continue\n",
    "    # add all {dna_window_size}-mers in seq to set, including overlapping\n",
    "    # ones\n",
    "    sequence = str(seq.seq)\n",
    "    for i in range(len(sequence)-(dna_window_size-1)):\n",
    "        upper_bound_all_rvdb_windows.add(sequence[i:(i+dna_window_size)])\n",
    "\n",
    "if excl_influenza:\n",
    "    print(f'{num_flu} sequences with flu')\n",
    "    print(f'{num_influenza} sequences with influenza')\n",
    "if excl_hiv:\n",
    "    print(f'''{num_human_immunodeficiency_virus} sequences with human\n",
    "    immunodeficiency virus''')\n",
    "    print(f'{num_hiv} sequences with hiv')\n",
    "print(f'''{len(upper_bound_all_rvdb_windows)} unique {dna_window_size}-mers\n",
    "      when excl_influenza = {excl_influenza}, excl_hiv = {excl_hiv}''')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In RVDBv15.1, there are\n",
    "\n",
    "- 2 sequences with flu\n",
    "- 661118 sequences with influenza\n",
    "- 814283 sequences with human immunodeficiency virus\n",
    "- 759673 sequences with hiv\n",
    "\n",
    "There are **471495384 unique 39-mers** and **486339373 unique 42-mers** when `excl_influenza = True, excl_hiv = True`."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}