SARS-CoV-2.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Imports & constants & tunable parameters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 53,
   "metadata": {},
   "outputs": [],
   "source": [
    "import math\n",
    "from os import listdir\n",
    "from re import search\n",
    "\n",
    "# affinity for susceptibility (negative/positive) - indices and directories\n",
    "NEG, POS, TEST = 0, 1, 2\n",
    "\n",
    "# sum of significance threshold parameter\n",
    "ALPHA = 0.5\n",
    "\n",
    "# weight function for a mutation given total # mutations in sequence\n",
    "def weight(n):\n",
    "  return (1 / (math.log(n) + 2))\n",
    "\n",
    "# length of human ACE2 sequence\n",
    "SEQ_LEN = 805"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Parse sequences from FASTA files\n",
    "\n",
    "The sequences are divided into:\n",
    "<ol start=\"0\">\n",
    "  <li>Known insusceptible sequences</li>\n",
    "  <li>Known susceptible sequences</li>\n",
    "  <li>Test sequences for predictions</li>\n",
    "</ol>\n",
    "\n",
    "Sequences are stored in tuples as (species name, sequence)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 54,
   "metadata": {},
   "outputs": [],
   "source": [
    "# complete list of data in format (organism, sequence) for each affinity\n",
    "seqs = [[], [], []]\n",
    "\n",
    "for aff in [NEG, POS, TEST]:\n",
    "  for filename in listdir('./data/%d' % aff):\n",
    "    with open(\"./data/%d/%s\" % (aff, filename)) as f:\n",
    "      lines = f.readlines()\n",
    "      name = search(r'\\[(.*?)\\]',lines[0][1:]).group(1)\n",
    "      sequence = ''.join(map(lambda k: k.strip(), lines[1:])).strip()\n",
    "      seqs[aff].append((name, sequence))\n",
    "\n",
    "# seqs"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Consensus Sequence\n",
    "\n",
    "The sequence made up of the most common residues at each site among the susceptible sequences"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 55,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'MSGSFWLLLSLVAVTAAQSTTEELAKTFLEKFNHEAEDLSYQSSLASWNYNTNITDENVQKMNEAGAKWSAFYEEQSKLAKTYPLQEIQNLTVKRQLQALQQSGSSVLSADKSKRLNTILNTMSTIYSTGKVCNPNNPQECLLLEPGLDDIMENSKDYNERLWAWEGWRSEVGKQLRPLYEEYVVLKNEMARANNYEDYGDYWRGDYEVEGVDGYDYSRDQLIEDVERTFTEIKPLYEHLHAYVRAKLMDAYPSYISPTGCLPAHLLGDMWGRFWTNLYSLTVPFGQKPNIDVTDAMVNQSWDARRIFKEAEKFFVSVGLPNMTQGFWENSMLTEPGDGRKVVCHPTAWDLGKGDFRIKMCTKVTMDDFLTAHHEMGHIQYDMAYAAQPFLLRNGANEGFHEAVGEIMSLSAATPNHLKSIGLLSPDFSEDNETEINFLLKQALTIVGTLPFTYMLEKWRWMVFKGEIPKEQWMQKWWEMKREIVGVVEPVPHDETYCDPASLFHVANDYSFIRYYTRTIYQFQFQEALCQIAKHEGPLHKCDISNSTEAGQKLLQMLSLGKSEPWTLALENVVGAKNMDVRPLLNYFEPLFTWLKEQNRNSFVGWSTDWSPYADQSIKVRISLKSALGDKAYEWNDNEMYLFRSSVAYAMREYFLKVKNQTIPFGEEDVRVSDLKPRISFNFFVTAPKNVSDIIPRTEVEEAIRMSRSRINDAFRLDDNSLEFLGIQPTLEPPYQPPVTIWLIVFGVVMGVVVVGIVLLIFTGIRDRRKKNQARSEENPYASVDLSKGENNPGFQNTDDVQTSF'"
      ]
     },
     "execution_count": 55,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "consensus_seq = \"\"\n",
    "\n",
    "for i in range(SEQ_LEN):\n",
    "  acids = []\n",
    "  for (name, seq) in seqs[POS]:\n",
    "    acids.append(seq[i])\n",
    "  consensus_seq += max(acids, key=acids.count)\n",
    "\n",
    "consensus_seq"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Identifying influential mutations\n",
    "\n",
    "Find mutations that exist exclusively in insusceptible sequences, assign a weight to each based on the total count, and sum the weights for each negative sequence. The mutations whose weights exceed a defined threshold (alpha) are recorded."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 56,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[(31, 'D', 0.6863135611816127),\n",
       " (41, 'A', 0.5),\n",
       " (66, 'A', 0.5385357335328163),\n",
       " (83, 'F', 0.5226018801516173),\n",
       " (113, 'N', 0.525115887738091),\n",
       " (353, 'H', 0.8362883189700047),\n",
       " (426, 'S', 0.5226018801516173),\n",
       " (679, 'V', 0.525115887738091)]"
      ]
     },
     "execution_count": 56,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "\n",
    "mutts = {}\n",
    "\n",
    "for (name, seq) in seqs[NEG]:\n",
    "  n_mutations = 0\n",
    "  ix_mutation = []\n",
    "  for (i, ch_i) in enumerate(seq):\n",
    "    if (ch_i) == \"-\": continue\n",
    "    has_match = False\n",
    "    for (_, seq_pos) in seqs[POS]:\n",
    "      if (seq_pos[i] == ch_i): \n",
    "        has_match = True\n",
    "        break\n",
    "    if (not has_match): \n",
    "      n_mutations += 1\n",
    "      mut = \"%c%d\" % (ch_i, i + 1)\n",
    "      ix_mutation.append(mut) \n",
    "  for i in ix_mutation:\n",
    "    # computed weighted value = 1 / (log(n) + 2)\n",
    "    wv = weight(n_mutations)\n",
    "    mutts[i] = wv if i not in mutts else mutts[i] + wv\n",
    "\n",
    "result = []\n",
    "for (m, f) in mutts.items():\n",
    "  if (f >= ALPHA): \n",
    "    m, i = m[0], int(m[1:])\n",
    "    result.append((i, m, f))\n",
    "result.sort(key=lambda k: k[0])\n",
    "\n",
    "result"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Important sites\n",
    "\n",
    "Print residues at influential sites for susceptible (including the consensus sequence) and then insusceptible sequences."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 57,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "  31  41  66  83 113 353 426 679\n",
      "  K   Y   G   Y   S   K   P   I   (+) Felis catus\n",
      "  K   Y   G   Y   S   K   P   I   (+) Chlorocebus sabaeus\n",
      "  K   Y   G   Y   S   K   P   I   (+) Rhinolophus macrotis\n",
      "  K   Y   E   Y   S   K   Y   I   (+) Oryctolagus cuniculus\n",
      "  K   H   G   Y   S   K   P   I   (+) Callithrix jacchus\n",
      "  K   Y   G   Y   S   K   P   I   (+) Macaca fascicularis\n",
      "  K   Y   R   Y   S   K   P   -   (+) Odocoileus virginianus texanus\n",
      "  K   Y   G   Y   S   K   P   I   (+) Panthera tigris\n",
      "  K   Y   G   Y   R   K   P   I   (+) Mustela putorius furo\n",
      "  K   Y   R   H   T   K   P   I   (+) Cynopterus sphinx\n",
      "  K   Y   G   Y   R   K   P   I   (+) Mustela lutreola biedermanni\n",
      "  K   Y   G   Y   S   K   P   I   (+) Homo sapiens\n",
      "  K   Y   G   Y   S   K   P   I   (+) Consensus sequence\n",
      "\n",
      "  -   A   -   -   -   -   -   -   (-) Mutation A41\n",
      "  N   Y   A   Y   R   K   P   I   (-) Procyon lotor\n",
      "  D   Y   -   -   -   -   -   -   (-) Mutation D31\n",
      "  -   -   -   -   -   H   -   -   (-) Mutation H353\n",
      "  N   Y   A   F   N   H   S   V   (-) Mus musculus\n",
      "  K   Y   G   Y   N   K   P   V   (-) Manis pentadactyla\n",
      "  K   Y   A   F   N   H   S   V   (-) Rattus norvegicus\n",
      "  D   H   G   F   S   K   S   I   (-) Rhinolophus ferrumequinum\n"
     ]
    }
   ],
   "source": [
    "indices = list(map(lambda k: k[0], result))\n",
    "\n",
    "def print_influence_pts(seq, name, aff):\n",
    "  p_str = \"  \"\n",
    "  for i in indices: p_str += seq[i-1] + \"   \"\n",
    "  p_str += \"(%s) %s\" % (aff, name)\n",
    "  print(p_str)\n",
    "\n",
    "print(''.join([\"%4d\" % (i) for i in indices]))\n",
    "for (name, seq) in seqs[POS]: print_influence_pts(seq, name, \"+\") \n",
    "print_influence_pts(consensus_seq, \"Consensus sequence\", \"+\")\n",
    "print()\n",
    "for (name, seq) in seqs[NEG]: print_influence_pts(seq, name, \"-\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Predictions based on results\n",
    "\n",
    "Predict the affinity of test sequences based on the presence or absence of influential mutations"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 58,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "  31  41  66  83 113 353 426 679\n",
      "  F   E   A   S   S   K   P   V   (+) Meleagris gallopavo\n",
      "  K   Y   R   Y   S   K   P   I   (-) Capra hircus\n",
      "  K   Y   G   F   S   K   P   I   (+) Mutation NFS82\n",
      "  -   -   -   -   -   K   P   I   (+) Columba livia\n",
      "  K   Y   A   Y   N   K   S   V   (-) Mesocricetus auratus\n",
      "  K   Y   G   Y   N   R   P   I   (-) Nyctereutes procyonoides\n",
      "  F   E   A   S   S   K   P   V   (+) Coturnix japonica\n",
      "  F   E   A   S   S   K   P   V   (+) Gallus gallus\n",
      "\n"
     ]
    }
   ],
   "source": [
    "def predict(seq):\n",
    "  for (i, m, _) in result:\n",
    "    if seq[i] == m: return \"-\"\n",
    "  return \"+\"\n",
    "\n",
    "print(''.join([\"%4d\" % (i) for i in indices]))\n",
    "for (name, seq) in seqs[TEST]: print_influence_pts(seq, name, predict(seq)) \n",
    "print()\n"
   ]
  }
 ],
 "metadata": {
  "interpreter": {
   "hash": "31f2aee4e71d21fbe5cf8b01ff0e069b9275f58929596ceb00d14d90e3e16cd6"
  },
  "kernelspec": {
   "display_name": "Python 3.7.3 64-bit",
   "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.7.3"
  },
  "orig_nbformat": 4
 },
 "nbformat": 4,
 "nbformat_minor": 2
}