sars-spike-protein-analysis.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "d6163d1e",
   "metadata": {},
   "outputs": [],
   "source": [
    "from Bio import SeqIO\n",
    "#make a fasta file with reference genome and other variant genomes, align them with mafft"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "7c368bb0",
   "metadata": {},
   "outputs": [],
   "source": [
    "fasta = list(SeqIO.parse('corona-seq.fasta', 'fasta'))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "4e04cd7f",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Wuhan-Hu-1\n",
      "B.1.1.7|Alpha\n",
      "B.1.351|Beta\n",
      "P.1,\n",
      "B.1.617.2|Delta\n",
      "BA.1|Omicron\n"
     ]
    }
   ],
   "source": [
    "for entry in fasta:\n",
    "    print(entry.id)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "360f19f6",
   "metadata": {},
   "outputs": [],
   "source": [
    "seqs = {}\n",
    "for entry in fasta:\n",
    "    seqs[entry.id] = entry"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "20facd00",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "29903\n",
      "29694\n",
      "29885\n",
      "29859\n",
      "29821\n",
      "29684\n"
     ]
    }
   ],
   "source": [
    "for entry in fasta:\n",
    "    print(len(entry.seq))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "id": "e31a4021",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "##### atgcfreq = 0.999994\n",
      "nseq =  6\n",
      "distance =  ktuples\n",
      "iterate =  0\n",
      "cycle =  2\n",
      "sparsepickup = 0\n",
      "nguidetree = 2\n",
      "nthread = 0\n",
      "sueff_global = 0.100000\n",
      "generating a scoring matrix for nucleotide (dist=200) ... done\n",
      "done\n",
      "done\n",
      "scoremtx = -1\n",
      "Gap Penalty = -1.53, +0.00, +0.00\n",
      "\n",
      "tuplesize = 6, dorp = d\n",
      "\n",
      "\n",
      "Making a distance matrix ..\n",
      "\n",
      "There are 30 ambiguous characters.\n",
      "    1 / 6\n",
      "done.\n",
      "\n",
      "Constructing a UPGMA tree ... \n",
      "    0 / 6\n",
      "done.\n",
      "\n",
      "Progressive alignment 1/2... \n",
      "STEP     5 / 5 f\n",
      "done.\n",
      "\n",
      "Making a distance matrix from msa.. \n",
      "    0 / 6\n",
      "done.\n",
      "\n",
      "Constructing a UPGMA tree ... \n",
      "    0 / 6\n",
      "done.\n",
      "\n",
      "Progressive alignment 2/2... \n",
      "STEP     5 / 5 f\n",
      "done.\n",
      "\n",
      "disttbfast (nuc) Version 7.245 alg=A, model=DNA200 (2), 1.53 (4.59), -0.00 (-0.00), noshift, amax=0.0\n",
      "0 thread(s)\n",
      "\n",
      "\n",
      "Strategy:\n",
      " FFT-NS-2 (Fast but rough)\n",
      " Progressive method (guide trees were built 2 times.)\n",
      "\n",
      "If unsure which option to use, try 'mafft --auto input > output'.\n",
      "For more information, see 'mafft --help', 'mafft --man' and the mafft page.\n",
      "\n",
      "The default gap scoring scheme has been changed in version 7.110 (2013 Oct).\n",
      "It tends to insert more gaps into gap-rich regions than previous versions.\n",
      "To disable this change, add the --leavegappyregion option.\n",
      "\n"
     ]
    }
   ],
   "source": [
    "!mafft --auto corona-seq.fasta > corona-seq-ali.fasta"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "d6db1ed0",
   "metadata": {},
   "outputs": [],
   "source": [
    "fasta = list(SeqIO.parse('corona-seq-ali.fasta', 'fasta'))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "1773554e",
   "metadata": {},
   "outputs": [],
   "source": [
    "seqs = {}\n",
    "for entry in fasta:\n",
    "    seqs[entry.id] = entry"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "c199b68d",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "29916\n",
      "29916\n",
      "29916\n",
      "29916\n",
      "29916\n",
      "29916\n"
     ]
    }
   ],
   "source": [
    "for entry in fasta:\n",
    "    print(len(entry.seq))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "c076b058",
   "metadata": {},
   "outputs": [],
   "source": [
    "def gapped_pos(seq, pos):\n",
    "    non_gap = 0\n",
    "    gaps = 0\n",
    "    \n",
    "    for nt in seq:\n",
    "        if nt != '-':\n",
    "            non_gap += 1\n",
    "        else:\n",
    "            gaps += 1\n",
    "        if non_gap == pos:\n",
    "            return pos + gaps"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "id": "686168c5",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "8"
      ]
     },
     "execution_count": 26,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "gapped_pos(\"aaa-aa-a\", 6)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fd725b8f",
   "metadata": {},
   "outputs": [],
   "source": [
    "#21563..25384 we have the region of the spike protein from the original , we need to right a function to get the \n",
    "#spike protein region from the gapped or the aligned sequences!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "87324142",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "21563"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "gapped_pos (seqs['Wuhan-Hu-1'].seq, 21563)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "0ba8e201",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "25393"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "gapped_pos (seqs['Wuhan-Hu-1'].seq, 25384)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "c3f96b73",
   "metadata": {},
   "outputs": [],
   "source": [
    "spikes = {}\n",
    "for seq in fasta:\n",
    "    spikes[seq.id] = seq.seq[21563-1:25393]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "6ec82e2d",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "atgtttgttt\n",
      "atgtttgttt\n",
      "atgtttgttt\n",
      "atgtttgttt\n",
      "atgtttgttt\n",
      "atgtttgttt\n"
     ]
    }
   ],
   "source": [
    "for spike in spikes:\n",
    "    print(spikes[spike][0:10])\n",
    "    \n",
    "#start codon is there"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "0eb96436",
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_mutations(initial, variant):\n",
    "    seqs = list(zip(initial, variant))\n",
    "    for pos, nt in enumerate(seqs):\n",
    "        if nt[0] != nt[1]:\n",
    "            print(nt[0].upper() + str(pos) + nt[1].upper())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "id": "4f8baa56",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "T202-\n",
      "A203-\n",
      "C204-\n",
      "A205-\n",
      "T206-\n",
      "G207-\n",
      "T428-\n",
      "T429-\n",
      "A430-\n",
      "A1509T\n",
      "C1717A\n",
      "A1849G\n",
      "C2050A\n",
      "C2155T\n",
      "T2952G\n",
      "G3360C\n"
     ]
    }
   ],
   "source": [
    "get_mutations(spikes['Wuhan-Hu-1'],spikes['B.1.1.7|Alpha']) #comparing reference wuhan gneome spike protein \n",
    "#to alpha genome spike protein"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "id": "657c95e2",
   "metadata": {},
   "outputs": [],
   "source": [
    "with open ('spikes.fasta', 'w') as f:\n",
    "    for spike in spikes:\n",
    "        out = spikes[spike].replace('-', '').translate()\n",
    "        f.write('>' + spike + '\\n')\n",
    "        f.write(str(out).upper()+ '\\n')\n",
    "        \n",
    "#this write a fasta file of all spike proteins of different variants\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "id": "e081cc11",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "corona-seq-ali.fasta\t\t\tmpxv-ha-sequences.fasta\r\n",
      "corona-seq.fasta\t\t\tMusic\r\n",
      "Desktop\t\t\t\t\tPictures\r\n",
      "Documents\t\t\t\tPublic\r\n",
      "Downloads\t\t\t\tsars-spike-protein-analysis.ipynb\r\n",
      "google-chrome-stable_current_amd64.deb\tsnap\r\n",
      "miniconda3\t\t\t\tspikes.fasta\r\n",
      "Miniconda3-latest-Linux-x86_64.sh\tTemplates\r\n",
      "mpx-seq-analystis-notebook.ipynb\tVideos\r\n",
      "mpxv-ha-sequences-ali.fasta\r\n"
     ]
    }
   ],
   "source": [
    "!ls"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "id": "0056d6a6",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "##### atgcfreq = 0.302478\n",
      "nseq =  6\n",
      "distance =  local\n",
      "iterate =  16\n",
      "cycle =  1\n",
      "nthread = 0\n",
      "lastonce = 0\n",
      "done.\n",
      "scoremtx = 1\n",
      "    4 / 6\n",
      "\n",
      "##### writing hat3\n",
      "pairlocalalign (aa) Version 7.245 alg=L, model=BLOSUM62, 2.00, -0.10, +0.10, noshift, amax=0.0\n",
      "0 thread(s)\n",
      "minimumweight = 0.000500\n",
      "nthread = 0\n",
      "blosum 62 / kimura 200\n",
      "sueff_global = 0.100000\n",
      "Loading 'hat3' ... \n",
      "done.\n",
      "done.\n",
      "scoremtx = 1\n",
      "Gap Penalty = -1.53, +0.00, +0.00\n",
      "Loading 'hat2' ... done.\n",
      "Constructing a UPGMA tree ... \n",
      "    0 / 6\n",
      "done.\n",
      "\n",
      "Progressive alignment ... \n",
      "STEP     5 /5 c\n",
      "done.\n",
      "tbfast (aa) Version 7.245 alg=A, model=BLOSUM62, 1.53, -0.00, -0.00, noshift, amax=0.0\n",
      "0 thread(s)\n",
      "minimumweight = 0.000500\n",
      "autosubalignment = 0.000000\n",
      "nthread = 0\n",
      "randomseed = 0\n",
      "blosum 62 / kimura 200\n",
      "poffset = 0\n",
      "niter = 16\n",
      "sueff_global = 0.100000\n",
      "Loading 'hat3' ... done.\n",
      "done.\n",
      "scoremtx = 1\n",
      "\n",
      "    0 / 6\n",
      "Segment   1/  1    1-1277\n",
      "STEP 002-004-1  identical.    identical.   \n",
      "Converged.\n",
      "\n",
      "done\n",
      "dvtditr (aa) Version 7.245 alg=A, model=BLOSUM62, 1.53, -0.00, -0.00, noshift, amax=0.0\n",
      "0 thread(s)\n",
      "\n",
      "\n",
      "Strategy:\n",
      " L-INS-i (Probably most accurate, very slow)\n",
      " Iterative refinement method (<16) with LOCAL pairwise alignment information\n",
      "\n",
      "If unsure which option to use, try 'mafft --auto input > output'.\n",
      "For more information, see 'mafft --help', 'mafft --man' and the mafft page.\n",
      "\n",
      "The default gap scoring scheme has been changed in version 7.110 (2013 Oct).\n",
      "It tends to insert more gaps into gap-rich regions than previous versions.\n",
      "To disable this change, add the --leavegappyregion option.\n",
      "\n"
     ]
    }
   ],
   "source": [
    "#now align amino acid sequnces of spike proteins\n",
    "!mafft --auto spikes.fasta > spikes-ali.fasta"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "id": "3e08b829",
   "metadata": {},
   "outputs": [],
   "source": [
    "spikes_aa = list(SeqIO.parse('spikes-ali.fasta', format = 'fasta'))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "id": "d86d9ba8",
   "metadata": {},
   "outputs": [],
   "source": [
    "seqs = {}\n",
    "for entry in spikes_aa:\n",
    "    seqs[entry.id]= entry.seq"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "id": "d9e15165",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[SeqRecord(seq=Seq('MFVFLVLLPLVSSQCVNLTTRTQLPPAYTNSFTRGVYYPDKVFRSSVLHSTQDL...HYT'), id='Wuhan-Hu-1', name='Wuhan-Hu-1', description='Wuhan-Hu-1', dbxrefs=[]),\n",
       " SeqRecord(seq=Seq('MFVFLVLLPLVSSQCVNLTTRTQLPPAYTNSFTRGVYYPDKVFRSSVLHSTQDL...HYT'), id='B.1.1.7|Alpha', name='B.1.1.7|Alpha', description='B.1.1.7|Alpha', dbxrefs=[]),\n",
       " SeqRecord(seq=Seq('MFVFLVLLPLVSSQCVNLTTRTQLPPAYTNSFTRGVYYPDKVFRSSVLHSTQDL...HYT'), id='B.1.351|Beta', name='B.1.351|Beta', description='B.1.351|Beta', dbxrefs=[]),\n",
       " SeqRecord(seq=Seq('MFVFLVLLPLVSSQCVNFTNRTQLPSAYTNSFTRGVYYPDKVFRSSVLHSTQDL...HYT'), id='P.1,', name='P.1,', description='P.1,', dbxrefs=[]),\n",
       " SeqRecord(seq=Seq('MFVFLVLLPLVSSQCVNLRTRTQLPPAYTNSFTRGVYYPDKVFRSSVLHSTQDL...HYT'), id='B.1.617.2|Delta', name='B.1.617.2|Delta', description='B.1.617.2|Delta', dbxrefs=[]),\n",
       " SeqRecord(seq=Seq('MFVFLVLLPLVSSQCVNLTTRTQLPPAYTNSFTRGVYYPDKVFRSSVLHSTQDL...HYT'), id='BA.1|Omicron', name='BA.1|Omicron', description='BA.1|Omicron', dbxrefs=[])]"
      ]
     },
     "execution_count": 26,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "spikes_aa"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "id": "fcae5198",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Seq('MFVFLVLLPLVSSQCVNLTTRTQLPPAYTNSFTRGVYYPDKVFRSSVLHSTQDL...HYT')"
      ]
     },
     "execution_count": 29,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "seqs['Wuhan-Hu-1']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "id": "2aef3d1c",
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_aa_mutations(initial, variant):\n",
    "    seqs = list(zip(initial, variant))\n",
    "    for pos, aa in enumerate(seqs):\n",
    "        if aa[0] != aa[1]:\n",
    "            print(aa[0].upper() + str(pos) + aa[1].upper())\n",
    "\n",
    "#find mutations in spike protein amino acid sequence"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "id": "bc9a1948",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "H68-\n",
      "V69-\n",
      "Y143-\n",
      "N503Y\n",
      "A572D\n",
      "D616G\n",
      "P683H\n",
      "T718I\n",
      "S984A\n",
      "D1120H\n"
     ]
    }
   ],
   "source": [
    "get_aa_mutations (seqs['Wuhan-Hu-1'], seqs ['B.1.1.7|Alpha'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "id": "8efea0e3",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "D79A\n",
      "D217G\n",
      "L243-\n",
      "L244-\n",
      "A245-\n",
      "K419N\n",
      "E486K\n",
      "N503Y\n",
      "D616G\n",
      "A703V\n"
     ]
    }
   ],
   "source": [
    "get_aa_mutations (seqs['Wuhan-Hu-1'], seqs ['B.1.351|Beta'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "id": "8958208f",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "A66V\n",
      "H68-\n",
      "V69-\n",
      "T94I\n",
      "G141D\n",
      "V142-\n",
      "Y143-\n",
      "Y144-\n",
      "N210-\n",
      "L211I\n",
      "-214E\n",
      "-215P\n",
      "-216E\n",
      "G341D\n",
      "S373L\n",
      "S375P\n",
      "S377F\n",
      "K419N\n",
      "N442K\n",
      "G448S\n",
      "S479N\n",
      "T480K\n",
      "E486A\n",
      "Q495R\n",
      "G498S\n",
      "Q500R\n",
      "N503Y\n",
      "Y507H\n",
      "T549K\n",
      "D616G\n",
      "H657Y\n",
      "N681K\n",
      "P683H\n",
      "N766K\n",
      "D798Y\n",
      "N858K\n",
      "Q956H\n",
      "N971K\n",
      "L983F\n"
     ]
    }
   ],
   "source": [
    "get_aa_mutations (seqs['Wuhan-Hu-1'], seqs ['BA.1|Omicron'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "id": "1e208541",
   "metadata": {},
   "outputs": [],
   "source": [
    "#make a list of mutations - we need to modify previous function\n",
    "def get_aa_mutations(initial, variant):\n",
    "    out = []\n",
    "    seqs = list(zip(initial, variant))\n",
    "    for pos, aa in enumerate(seqs):\n",
    "        if aa[0] != aa[1]:\n",
    "            out.append(aa[0].upper() + str(pos) + aa[1].upper())\n",
    "    return out"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "id": "6f78e8af",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "['H68-',\n",
       " 'V69-',\n",
       " 'Y143-',\n",
       " 'N503Y',\n",
       " 'A572D',\n",
       " 'D616G',\n",
       " 'P683H',\n",
       " 'T718I',\n",
       " 'S984A',\n",
       " 'D1120H']"
      ]
     },
     "execution_count": 46,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "get_aa_mutations (seqs['Wuhan-Hu-1'], seqs ['B.1.1.7|Alpha'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 47,
   "id": "9f92247e",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "10"
      ]
     },
     "execution_count": 47,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#How many mutations are there? \n",
    "len(get_aa_mutations (seqs['Wuhan-Hu-1'], seqs ['B.1.1.7|Alpha']))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 48,
   "id": "e63834fb",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "10"
      ]
     },
     "execution_count": 48,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "len(get_aa_mutations (seqs['Wuhan-Hu-1'], seqs ['B.1.351|Beta']))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 49,
   "id": "fd5f1d0f",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "39"
      ]
     },
     "execution_count": 49,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "len(get_aa_mutations (seqs['Wuhan-Hu-1'], seqs ['BA.1|Omicron']))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 54,
   "id": "e4695883",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Wuhan-Hu-1 0\n",
      "B.1.1.7|Alpha 10\n",
      "B.1.351|Beta 10\n",
      "P.1, 12\n",
      "B.1.617.2|Delta 10\n",
      "BA.1|Omicron 39\n"
     ]
    }
   ],
   "source": [
    "for item in seqs:\n",
    "    print(item + ' '+ str(len(get_aa_mutations (seqs['Wuhan-Hu-1'], seqs[item]))))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d4c9413a",
   "metadata": {},
   "outputs": [],
   "source": [
    "#it seems omicron has the highest mutation rate against reference genome!!!"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.10.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}