sars-cov-2/1-MapTo3DStructures.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Maps SARS-CoV-2 Mutations to 3D Protein Structures\n",
    "[Work in progress]\n",
    "\n",
    "This notebook maps mutation frequency of SARS-CoV-2 strains onto 3D protein structures with bound antibody fragments in the [Protein Data Bank](https://www.wwpdb.org/). SARS-CoV2-2 strains and their mutations have been aggregated in the [COVID-19-Net Knowledge Graph](https://github.com/covid-19-net/covid-19-community)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import math\n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib.colors\n",
    "import matplotlib.cm as cm\n",
    "import ipywidgets as widgets\n",
    "from ipywidgets import interact, IntSlider, FloatSlider, SelectMultiple\n",
    "from py2neo import Graph\n",
    "import py3Dmol"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "pd.options.display.max_rows = None  # display all rows\n",
    "pd.options.display.max_columns = None  # display all columsns"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "output_file_name = 'mutations3d.csv' # mutations mapped to 3D protein structures"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Connect to COVID-19-Community Knowledge Graph\n",
    "[COVID-19-Net Knowledge Graph](https://github.com/covid-19-net/covid-19-community)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "graph = Graph(\"bolt://132.249.238.185:7687\", user=\"reader\", password=\"demo\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Taxonomy ids for pathogen (SARS-CoV-2) and host (human)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "pathogen_taxonomy_id = 'taxonomy:2697049'\n",
    "host_taxonomy_id = 'taxonomy:9606'\n",
    "with_antibodies = True"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Get list of SARS-CoV-2 proteins"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "query = \"\"\"\n",
    "MATCH (p:Protein{taxonomyId: $pathogen_taxonomy_id})-[t:HAS_TERTIARY_STRUCTURE]->(:Chain)-[:IS_PART_OF_STRUCTURE]->(s:Structure)\n",
    "WHERE t.coverage > 0.2 // eliminate polyprotein\n",
    "RETURN p.name AS protein, p.accession as accession, p.proId as proId, p.sequence as sequence, s.description AS description\n",
    "ORDER BY protein\n",
    "\"\"\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "structures = graph.run(query, pathogen_taxonomy_id=pathogen_taxonomy_id).to_data_frame()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "proteins = structures[['protein', 'accession', 'proId', 'sequence']].drop_duplicates()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "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>protein</th>\n",
       "      <th>accession</th>\n",
       "      <th>proId</th>\n",
       "      <th>sequence</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>2'-O-methyltransferase</td>\n",
       "      <td>uniprot:P0DTD1</td>\n",
       "      <td>uniprot.chain:PRO_0000449633</td>\n",
       "      <td>SSQAWQPGVAMPNLYKMQRMLLEKCDLQNYGDSATLPKGIMMNVAK...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>28</th>\n",
       "      <td>3C-like proteinase</td>\n",
       "      <td>uniprot:P0DTC1</td>\n",
       "      <td>uniprot.chain:PRO_0000449639</td>\n",
       "      <td>SGFRKMAFPSGKVEGCMVQVTCGTTTLNGLWLDDVVYCPRHVICTS...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>434</th>\n",
       "      <td>Envelope small membrane protein</td>\n",
       "      <td>uniprot:P0DTC4</td>\n",
       "      <td>uniprot.chain:PRO_0000449651</td>\n",
       "      <td>MYSFVSEETGTLIVNSVLLFLAFVVFLLVTLAILTALRLCAYCCNI...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>439</th>\n",
       "      <td>Helicase</td>\n",
       "      <td>uniprot:P0DTD1</td>\n",
       "      <td>uniprot.chain:PRO_0000449630</td>\n",
       "      <td>AVGACVLCNSQTSLRCGACIRRPFLCCKCCYDHVISTSHKLVLSVN...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>564</th>\n",
       "      <td>Host translation inhibitor nsp1</td>\n",
       "      <td>uniprot:P0DTC1</td>\n",
       "      <td>uniprot.chain:PRO_0000449635</td>\n",
       "      <td>MESLVPGFNEKTHVQLSLPVLQVRDVLVRGFGDSVEEVLSEARQHL...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>566</th>\n",
       "      <td>Non-structural protein 10</td>\n",
       "      <td>uniprot:P0DTC1</td>\n",
       "      <td>uniprot.chain:PRO_0000449644</td>\n",
       "      <td>AGNATEVPANSTVLSFCAFAVDAAKAYKDYLASGGQPITNCVKMLC...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>601</th>\n",
       "      <td>Non-structural protein 2</td>\n",
       "      <td>uniprot:P0DTC1</td>\n",
       "      <td>uniprot.chain:PRO_0000449636</td>\n",
       "      <td>AYTRYVDNNFCGPDGYPLECIKDLLARAGKASCTLSEQLDFIDTKR...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>603</th>\n",
       "      <td>Non-structural protein 7</td>\n",
       "      <td>uniprot:P0DTC1</td>\n",
       "      <td>uniprot.chain:PRO_0000449641</td>\n",
       "      <td>SKMSDVKCTSVVLLSVLQQLRVESSSKLWAQCVQLHNDILLAKDTT...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>642</th>\n",
       "      <td>Non-structural protein 8</td>\n",
       "      <td>uniprot:P0DTC1</td>\n",
       "      <td>uniprot.chain:PRO_0000449642</td>\n",
       "      <td>AIASEFSSLPSYAAFATAQEAYEQAVANGDSEVVLKKLKKSLNVAK...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>696</th>\n",
       "      <td>Non-structural protein 9</td>\n",
       "      <td>uniprot:P0DTC1</td>\n",
       "      <td>uniprot.chain:PRO_0000449643</td>\n",
       "      <td>NNELSPVALRQMSCAAGTTQTACTDDNALAYYNTTKGGRFVLALLS...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>704</th>\n",
       "      <td>Nucleoprotein</td>\n",
       "      <td>uniprot:P0DTC9</td>\n",
       "      <td>uniprot.chain:PRO_0000449656</td>\n",
       "      <td>MSDNGPQNQRNAPRITFGGPSDSTGSNQNGERSGARSKQRRPQGLP...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>753</th>\n",
       "      <td>ORF3a protein</td>\n",
       "      <td>uniprot:P0DTC3</td>\n",
       "      <td>uniprot.chain:PRO_0000449650</td>\n",
       "      <td>MDLFMRIFTIGTVTLKQGEIKDATPSDFVRATATIPIQASLPFGWL...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>757</th>\n",
       "      <td>ORF7a protein</td>\n",
       "      <td>uniprot:P0DTC7</td>\n",
       "      <td>None</td>\n",
       "      <td>MKIILFLALITLATCELYHYQECVRGTTVLLKEPCSSGTYEGNSPF...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>758</th>\n",
       "      <td>ORF7a protein</td>\n",
       "      <td>uniprot:P0DTC7</td>\n",
       "      <td>uniprot.chain:PRO_0000449654</td>\n",
       "      <td>ELYHYQECVRGTTVLLKEPCSSGTYEGNSPFHPLADNKFALTCFST...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>759</th>\n",
       "      <td>ORF8 protein</td>\n",
       "      <td>uniprot:P0DTC8</td>\n",
       "      <td>None</td>\n",
       "      <td>MKFLVFLGIITTVAAFHQECSLQSCTQHQPYVVDDPCPIHFYSKWY...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>763</th>\n",
       "      <td>ORF8 protein</td>\n",
       "      <td>uniprot:P0DTC8</td>\n",
       "      <td>uniprot.chain:PRO_0000449655</td>\n",
       "      <td>FHQECSLQSCTQHQPYVVDDPCPIHFYSKWYIRVGARKSAPLIELC...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>767</th>\n",
       "      <td>ORF9b protein</td>\n",
       "      <td>uniprot:P0DTD2</td>\n",
       "      <td>uniprot.chain:PRO_0000449657</td>\n",
       "      <td>MDPKISEMHPALRLVDPQIQLAVTRMENAVGRDQNNVGPKVYPIIL...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>771</th>\n",
       "      <td>Proofreading exoribonuclease</td>\n",
       "      <td>uniprot:P0DTD1</td>\n",
       "      <td>uniprot.chain:PRO_0000449631</td>\n",
       "      <td>AENVTGLFKDCSKVITGLHPTQAPTHLSVDTKFKTEGLCVDIPGIP...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>774</th>\n",
       "      <td>RNA-directed RNA polymerase</td>\n",
       "      <td>uniprot:P0DTD1</td>\n",
       "      <td>uniprot.chain:PRO_0000449629</td>\n",
       "      <td>SADAQSFLNRVCGVSAARLTPCGTGTSTDVVYRAFDIYNDKVAGFA...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>797</th>\n",
       "      <td>Spike glycoprotein</td>\n",
       "      <td>uniprot:P0DTC2</td>\n",
       "      <td>None</td>\n",
       "      <td>MFVFLVLLPLVSSQCVNLTTRTQLPPAYTNSFTRGVYYPDKVFRSS...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1573</th>\n",
       "      <td>Spike glycoprotein</td>\n",
       "      <td>uniprot:P0DTC2</td>\n",
       "      <td>uniprot.chain:PRO_0000449646</td>\n",
       "      <td>SQCVNLTTRTQLPPAYTNSFTRGVYYPDKVFRSSVLHSTQDLFLPF...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2349</th>\n",
       "      <td>Spike protein S1</td>\n",
       "      <td>uniprot:P0DTC2</td>\n",
       "      <td>uniprot.chain:PRO_0000449647</td>\n",
       "      <td>SQCVNLTTRTQLPPAYTNSFTRGVYYPDKVFRSSVLHSTQDLFLPF...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2542</th>\n",
       "      <td>Spike protein S2</td>\n",
       "      <td>uniprot:P0DTC2</td>\n",
       "      <td>uniprot.chain:PRO_0000449648</td>\n",
       "      <td>SVASQSIIAYTMSLGAENSVAYSNNSIAIPTNFTISVTTEILPVSM...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2546</th>\n",
       "      <td>Spike protein S2'</td>\n",
       "      <td>uniprot:P0DTC2</td>\n",
       "      <td>uniprot.chain:PRO_0000449649</td>\n",
       "      <td>SFIEDLLFNKVTLADAGFIKQYGDCLGDIAARDLICAQKFNGLTVL...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2552</th>\n",
       "      <td>Uridylate-specific endoribonuclease</td>\n",
       "      <td>uniprot:P0DTD1</td>\n",
       "      <td>uniprot.chain:PRO_0000449632</td>\n",
       "      <td>SLENVAFNVVNKGHFDGQQGEVPVSIINNTVYTKVDGVDVELFENK...</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                                  protein       accession  \\\n",
       "0                  2'-O-methyltransferase  uniprot:P0DTD1   \n",
       "28                     3C-like proteinase  uniprot:P0DTC1   \n",
       "434       Envelope small membrane protein  uniprot:P0DTC4   \n",
       "439                              Helicase  uniprot:P0DTD1   \n",
       "564       Host translation inhibitor nsp1  uniprot:P0DTC1   \n",
       "566             Non-structural protein 10  uniprot:P0DTC1   \n",
       "601              Non-structural protein 2  uniprot:P0DTC1   \n",
       "603              Non-structural protein 7  uniprot:P0DTC1   \n",
       "642              Non-structural protein 8  uniprot:P0DTC1   \n",
       "696              Non-structural protein 9  uniprot:P0DTC1   \n",
       "704                         Nucleoprotein  uniprot:P0DTC9   \n",
       "753                         ORF3a protein  uniprot:P0DTC3   \n",
       "757                         ORF7a protein  uniprot:P0DTC7   \n",
       "758                         ORF7a protein  uniprot:P0DTC7   \n",
       "759                          ORF8 protein  uniprot:P0DTC8   \n",
       "763                          ORF8 protein  uniprot:P0DTC8   \n",
       "767                         ORF9b protein  uniprot:P0DTD2   \n",
       "771          Proofreading exoribonuclease  uniprot:P0DTD1   \n",
       "774           RNA-directed RNA polymerase  uniprot:P0DTD1   \n",
       "797                    Spike glycoprotein  uniprot:P0DTC2   \n",
       "1573                   Spike glycoprotein  uniprot:P0DTC2   \n",
       "2349                     Spike protein S1  uniprot:P0DTC2   \n",
       "2542                     Spike protein S2  uniprot:P0DTC2   \n",
       "2546                    Spike protein S2'  uniprot:P0DTC2   \n",
       "2552  Uridylate-specific endoribonuclease  uniprot:P0DTD1   \n",
       "\n",
       "                             proId  \\\n",
       "0     uniprot.chain:PRO_0000449633   \n",
       "28    uniprot.chain:PRO_0000449639   \n",
       "434   uniprot.chain:PRO_0000449651   \n",
       "439   uniprot.chain:PRO_0000449630   \n",
       "564   uniprot.chain:PRO_0000449635   \n",
       "566   uniprot.chain:PRO_0000449644   \n",
       "601   uniprot.chain:PRO_0000449636   \n",
       "603   uniprot.chain:PRO_0000449641   \n",
       "642   uniprot.chain:PRO_0000449642   \n",
       "696   uniprot.chain:PRO_0000449643   \n",
       "704   uniprot.chain:PRO_0000449656   \n",
       "753   uniprot.chain:PRO_0000449650   \n",
       "757                           None   \n",
       "758   uniprot.chain:PRO_0000449654   \n",
       "759                           None   \n",
       "763   uniprot.chain:PRO_0000449655   \n",
       "767   uniprot.chain:PRO_0000449657   \n",
       "771   uniprot.chain:PRO_0000449631   \n",
       "774   uniprot.chain:PRO_0000449629   \n",
       "797                           None   \n",
       "1573  uniprot.chain:PRO_0000449646   \n",
       "2349  uniprot.chain:PRO_0000449647   \n",
       "2542  uniprot.chain:PRO_0000449648   \n",
       "2546  uniprot.chain:PRO_0000449649   \n",
       "2552  uniprot.chain:PRO_0000449632   \n",
       "\n",
       "                                               sequence  \n",
       "0     SSQAWQPGVAMPNLYKMQRMLLEKCDLQNYGDSATLPKGIMMNVAK...  \n",
       "28    SGFRKMAFPSGKVEGCMVQVTCGTTTLNGLWLDDVVYCPRHVICTS...  \n",
       "434   MYSFVSEETGTLIVNSVLLFLAFVVFLLVTLAILTALRLCAYCCNI...  \n",
       "439   AVGACVLCNSQTSLRCGACIRRPFLCCKCCYDHVISTSHKLVLSVN...  \n",
       "564   MESLVPGFNEKTHVQLSLPVLQVRDVLVRGFGDSVEEVLSEARQHL...  \n",
       "566   AGNATEVPANSTVLSFCAFAVDAAKAYKDYLASGGQPITNCVKMLC...  \n",
       "601   AYTRYVDNNFCGPDGYPLECIKDLLARAGKASCTLSEQLDFIDTKR...  \n",
       "603   SKMSDVKCTSVVLLSVLQQLRVESSSKLWAQCVQLHNDILLAKDTT...  \n",
       "642   AIASEFSSLPSYAAFATAQEAYEQAVANGDSEVVLKKLKKSLNVAK...  \n",
       "696   NNELSPVALRQMSCAAGTTQTACTDDNALAYYNTTKGGRFVLALLS...  \n",
       "704   MSDNGPQNQRNAPRITFGGPSDSTGSNQNGERSGARSKQRRPQGLP...  \n",
       "753   MDLFMRIFTIGTVTLKQGEIKDATPSDFVRATATIPIQASLPFGWL...  \n",
       "757   MKIILFLALITLATCELYHYQECVRGTTVLLKEPCSSGTYEGNSPF...  \n",
       "758   ELYHYQECVRGTTVLLKEPCSSGTYEGNSPFHPLADNKFALTCFST...  \n",
       "759   MKFLVFLGIITTVAAFHQECSLQSCTQHQPYVVDDPCPIHFYSKWY...  \n",
       "763   FHQECSLQSCTQHQPYVVDDPCPIHFYSKWYIRVGARKSAPLIELC...  \n",
       "767   MDPKISEMHPALRLVDPQIQLAVTRMENAVGRDQNNVGPKVYPIIL...  \n",
       "771   AENVTGLFKDCSKVITGLHPTQAPTHLSVDTKFKTEGLCVDIPGIP...  \n",
       "774   SADAQSFLNRVCGVSAARLTPCGTGTSTDVVYRAFDIYNDKVAGFA...  \n",
       "797   MFVFLVLLPLVSSQCVNLTTRTQLPPAYTNSFTRGVYYPDKVFRSS...  \n",
       "1573  SQCVNLTTRTQLPPAYTNSFTRGVYYPDKVFRSSVLHSTQDLFLPF...  \n",
       "2349  SQCVNLTTRTQLPPAYTNSFTRGVYYPDKVFRSSVLHSTQDLFLPF...  \n",
       "2542  SVASQSIIAYTMSLGAENSVAYSNNSIAIPTNFTISVTTEILPVSM...  \n",
       "2546  SFIEDLLFNKVTLADAGFIKQYGDCLGDIAARDLICAQKFNGLTVL...  \n",
       "2552  SLENVAFNVVNKGHFDGQQGEVPVSIINNTVYTKVDGVDVELFENK...  "
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "proteins.head(50)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "protein_list = structures['protein'].unique()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "protein_widget = widgets.Dropdown(options=protein_list, description='Select protein:', value='Spike glycoprotein',\n",
    "                                  style={'description_width': 'initial'}, )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Select SARS-CoV-2 Protein"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "e56b513d43004000bb6fbd74ba242338",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Dropdown(description='Select protein:', index=17, options=(\"2'-O-methyltransferase\", '3C-like proteinase', 'En…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "display(protein_widget)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Protein name: Spike glycoprotein\n"
     ]
    }
   ],
   "source": [
    "protein_name = protein_widget.value\n",
    "print('Protein name:', protein_name)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "subset = structures.query(f'protein == \"{protein_name}\"').copy()\n",
    "subset.fillna('', inplace=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "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>accession</th>\n",
       "      <th>description</th>\n",
       "      <th>proId</th>\n",
       "      <th>protein</th>\n",
       "      <th>sequence</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>797</th>\n",
       "      <td>uniprot:P0DTC2</td>\n",
       "      <td>2-acetamido-2-deoxy-beta-D-glucopyranose, 2-ac...</td>\n",
       "      <td></td>\n",
       "      <td>Spike glycoprotein</td>\n",
       "      <td>MFVFLVLLPLVSSQCVNLTTRTQLPPAYTNSFTRGVYYPDKVFRSS...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>798</th>\n",
       "      <td>uniprot:P0DTC2</td>\n",
       "      <td>2-acetamido-2-deoxy-beta-D-glucopyranose, 2-ac...</td>\n",
       "      <td></td>\n",
       "      <td>Spike glycoprotein</td>\n",
       "      <td>MFVFLVLLPLVSSQCVNLTTRTQLPPAYTNSFTRGVYYPDKVFRSS...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>799</th>\n",
       "      <td>uniprot:P0DTC2</td>\n",
       "      <td>Fab Heavy chain of enhancing antibody, Fab lig...</td>\n",
       "      <td></td>\n",
       "      <td>Spike glycoprotein</td>\n",
       "      <td>MFVFLVLLPLVSSQCVNLTTRTQLPPAYTNSFTRGVYYPDKVFRSS...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>800</th>\n",
       "      <td>uniprot:P0DTC2</td>\n",
       "      <td>Fab Heavy chain of enhancing antibody, Fab lig...</td>\n",
       "      <td></td>\n",
       "      <td>Spike glycoprotein</td>\n",
       "      <td>MFVFLVLLPLVSSQCVNLTTRTQLPPAYTNSFTRGVYYPDKVFRSS...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>801</th>\n",
       "      <td>uniprot:P0DTC2</td>\n",
       "      <td>2-acetamido-2-deoxy-beta-D-glucopyranose, 2-ac...</td>\n",
       "      <td></td>\n",
       "      <td>Spike glycoprotein</td>\n",
       "      <td>MFVFLVLLPLVSSQCVNLTTRTQLPPAYTNSFTRGVYYPDKVFRSS...</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "          accession                                        description proId  \\\n",
       "797  uniprot:P0DTC2  2-acetamido-2-deoxy-beta-D-glucopyranose, 2-ac...         \n",
       "798  uniprot:P0DTC2  2-acetamido-2-deoxy-beta-D-glucopyranose, 2-ac...         \n",
       "799  uniprot:P0DTC2  Fab Heavy chain of enhancing antibody, Fab lig...         \n",
       "800  uniprot:P0DTC2  Fab Heavy chain of enhancing antibody, Fab lig...         \n",
       "801  uniprot:P0DTC2  2-acetamido-2-deoxy-beta-D-glucopyranose, 2-ac...         \n",
       "\n",
       "                protein                                           sequence  \n",
       "797  Spike glycoprotein  MFVFLVLLPLVSSQCVNLTTRTQLPPAYTNSFTRGVYYPDKVFRSS...  \n",
       "798  Spike glycoprotein  MFVFLVLLPLVSSQCVNLTTRTQLPPAYTNSFTRGVYYPDKVFRSS...  \n",
       "799  Spike glycoprotein  MFVFLVLLPLVSSQCVNLTTRTQLPPAYTNSFTRGVYYPDKVFRSS...  \n",
       "800  Spike glycoprotein  MFVFLVLLPLVSSQCVNLTTRTQLPPAYTNSFTRGVYYPDKVFRSS...  \n",
       "801  Spike glycoprotein  MFVFLVLLPLVSSQCVNLTTRTQLPPAYTNSFTRGVYYPDKVFRSS...  "
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "subset.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Get total number of strains"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "query = \"\"\"\n",
    "MATCH (s:Strain)\n",
    "WHERE s.taxonomyId = $pathogen_taxonomy_id AND s.hostTaxonomyId = $host_taxonomy_id\n",
    "RETURN count(s)\n",
    "\"\"\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "strains = graph.evaluate(query, pathogen_taxonomy_id=pathogen_taxonomy_id, \n",
    "                         host_taxonomy_id=host_taxonomy_id)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Total number of human strains: 1857838\n"
     ]
    }
   ],
   "source": [
    "print('Total number of human strains:', strains)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Get variants for selected protein"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "query = \"\"\"\n",
    "MATCH (p:Protein{reviewed: True})-[:HAS_VARIANT]->(v:Variant{variantConsequence:'missense_variant'})<-[:HAS_VARIANT]-(s:Strain)\n",
    "WHERE p.name = $protein_name AND p.taxonomyId = $pathogen_taxonomy_id AND s.hostTaxonomyId = $host_taxonomy_id\n",
    "WITH v.proteinPosition AS residue, count(v.proteinVariant) AS count, \n",
    "     v.proteinVariant + '(' + count(v.proteinVariant) + ')' AS variationId ,\n",
    "     split(v.proteinVariant, ':')[1] + '(' + count(v.proteinVariant) + ')' AS annotation\n",
    "     ORDER by count DESC\n",
    "WITH residue, count, variationId, annotation\n",
    "RETURN residue, collect(variationId) AS variationId, collect(annotation) AS annotation, sum(count) AS count ORDER BY residue\n",
    "\"\"\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Add mutation annotation to each residue"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "variants = graph.run(query, protein_name=protein_name, pathogen_taxonomy_id=pathogen_taxonomy_id, \n",
    "                     host_taxonomy_id=host_taxonomy_id).to_data_frame()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(1242, 4)"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "variants.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "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>annotation</th>\n",
       "      <th>count</th>\n",
       "      <th>residue</th>\n",
       "      <th>variationId</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>[p.2F&gt;L(157), p.2F&gt;S(97), p.2F&gt;V(9), p.2F&gt;Y(2)...</td>\n",
       "      <td>266</td>\n",
       "      <td>2</td>\n",
       "      <td>[S:p.2F&gt;L(157), S:p.2F&gt;S(97), S:p.2F&gt;V(9), QHD...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>[p.3V&gt;G(607), p.3V&gt;I(36), p.3V&gt;F(8), p.3V&gt;A(4)...</td>\n",
       "      <td>657</td>\n",
       "      <td>3</td>\n",
       "      <td>[S:p.3V&gt;G(607), S:p.3V&gt;I(36), S:p.3V&gt;F(8), S:p...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>[p.4F&gt;I(48), p.4F&gt;S(9), p.4F&gt;L(6), p.4F&gt;V(1), ...</td>\n",
       "      <td>65</td>\n",
       "      <td>4</td>\n",
       "      <td>[S:p.4F&gt;I(48), S:p.4F&gt;S(9), S:p.4F&gt;L(6), S:p.4...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>[p.5L&gt;F(43645), p.5L&gt;V(4), p.5L&gt;I(4)]</td>\n",
       "      <td>43653</td>\n",
       "      <td>5</td>\n",
       "      <td>[S:p.5L&gt;F(43645), S:p.5L&gt;V(4), S:p.5L&gt;I(4)]</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>[p.6V&gt;A(765), p.6V&gt;F(498), p.6V&gt;I(96), p.6V&gt;L(...</td>\n",
       "      <td>1369</td>\n",
       "      <td>6</td>\n",
       "      <td>[S:p.6V&gt;A(765), S:p.6V&gt;F(498), S:p.6V&gt;I(96), S...</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                                          annotation  count  residue  \\\n",
       "0  [p.2F>L(157), p.2F>S(97), p.2F>V(9), p.2F>Y(2)...    266        2   \n",
       "1  [p.3V>G(607), p.3V>I(36), p.3V>F(8), p.3V>A(4)...    657        3   \n",
       "2  [p.4F>I(48), p.4F>S(9), p.4F>L(6), p.4F>V(1), ...     65        4   \n",
       "3              [p.5L>F(43645), p.5L>V(4), p.5L>I(4)]  43653        5   \n",
       "4  [p.6V>A(765), p.6V>F(498), p.6V>I(96), p.6V>L(...   1369        6   \n",
       "\n",
       "                                         variationId  \n",
       "0  [S:p.2F>L(157), S:p.2F>S(97), S:p.2F>V(9), QHD...  \n",
       "1  [S:p.3V>G(607), S:p.3V>I(36), S:p.3V>F(8), S:p...  \n",
       "2  [S:p.4F>I(48), S:p.4F>S(9), S:p.4F>L(6), S:p.4...  \n",
       "3        [S:p.5L>F(43645), S:p.5L>V(4), S:p.5L>I(4)]  \n",
       "4  [S:p.6V>A(765), S:p.6V>F(498), S:p.6V>I(96), S...  "
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "variants.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [],
   "source": [
    "variants.sort_values('count', ascending=False, inplace=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [],
   "source": [
    "variants['variationId'] = variants['variationId'].apply(lambda x: ', '.join(x))\n",
    "variants['annotation'] = variants['annotation'].apply(lambda x: ', '.join(x))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [],
   "source": [
    "variants['annotation'] = variants['annotation'].str.replace('p.', '')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Create a color scale based on the log mutation frequency"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [],
   "source": [
    "total = variants['count'].sum()\n",
    "variants['scale'] = variants['count'].apply(np.log10) / math.log10(total)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [],
   "source": [
    "n_colors = 100\n",
    "colors = cm.Reds(np.linspace(0.0, 1.0, n_colors))\n",
    "col = np.empty(n_colors, dtype=object)\n",
    "\n",
    "for i, color in enumerate(colors):\n",
    "    col[i] = matplotlib.colors.rgb2hex(color)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [],
   "source": [
    "variants['color'] = variants['scale'].apply(lambda x: col[round(x*n_colors)])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "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>annotation</th>\n",
       "      <th>count</th>\n",
       "      <th>residue</th>\n",
       "      <th>variationId</th>\n",
       "      <th>scale</th>\n",
       "      <th>color</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>599</th>\n",
       "      <td>614D&gt;G(2141444), 614D&gt;N(274), 614D&gt;A(24), 614D...</td>\n",
       "      <td>2141748</td>\n",
       "      <td>614</td>\n",
       "      <td>S:p.614D&gt;G(2141444), S:p.614D&gt;N(274), S:p.614D...</td>\n",
       "      <td>0.900988</td>\n",
       "      <td>#940b13</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>663</th>\n",
       "      <td>681P&gt;H(1069788), 681P&gt;R(20076), 681P&gt;L(778), 6...</td>\n",
       "      <td>1090816</td>\n",
       "      <td>681</td>\n",
       "      <td>S:p.681P&gt;H(1069788), S:p.681P&gt;R(20076), S:p.68...</td>\n",
       "      <td>0.859286</td>\n",
       "      <td>#a60f15</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>490</th>\n",
       "      <td>501N&gt;Y(1066528), 501N&gt;T(3858), 501N&gt;S(42), 501...</td>\n",
       "      <td>1070490</td>\n",
       "      <td>501</td>\n",
       "      <td>S:p.501N&gt;Y(1066528), S:p.501N&gt;T(3858), S:p.501...</td>\n",
       "      <td>0.858123</td>\n",
       "      <td>#a60f15</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>697</th>\n",
       "      <td>716T&gt;I(1031950), 716T&gt;A(338), 716T&gt;S(6), 716T&gt;...</td>\n",
       "      <td>1032296</td>\n",
       "      <td>716</td>\n",
       "      <td>S:p.716T&gt;I(1031950), S:p.716T&gt;A(338), S:p.716T...</td>\n",
       "      <td>0.855878</td>\n",
       "      <td>#a60f15</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1088</th>\n",
       "      <td>1118D&gt;H(1026532), 1118D&gt;Y(1474), 1118D&gt;V(30), ...</td>\n",
       "      <td>1028076</td>\n",
       "      <td>1118</td>\n",
       "      <td>S:p.1118D&gt;H(1026532), S:p.1118D&gt;Y(1474), QHD43...</td>\n",
       "      <td>0.855624</td>\n",
       "      <td>#a60f15</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                                             annotation    count  residue  \\\n",
       "599   614D>G(2141444), 614D>N(274), 614D>A(24), 614D...  2141748      614   \n",
       "663   681P>H(1069788), 681P>R(20076), 681P>L(778), 6...  1090816      681   \n",
       "490   501N>Y(1066528), 501N>T(3858), 501N>S(42), 501...  1070490      501   \n",
       "697   716T>I(1031950), 716T>A(338), 716T>S(6), 716T>...  1032296      716   \n",
       "1088  1118D>H(1026532), 1118D>Y(1474), 1118D>V(30), ...  1028076     1118   \n",
       "\n",
       "                                            variationId     scale    color  \n",
       "599   S:p.614D>G(2141444), S:p.614D>N(274), S:p.614D...  0.900988  #940b13  \n",
       "663   S:p.681P>H(1069788), S:p.681P>R(20076), S:p.68...  0.859286  #a60f15  \n",
       "490   S:p.501N>Y(1066528), S:p.501N>T(3858), S:p.501...  0.858123  #a60f15  \n",
       "697   S:p.716T>I(1031950), S:p.716T>A(338), S:p.716T...  0.855878  #a60f15  \n",
       "1088  S:p.1118D>H(1026532), S:p.1118D>Y(1474), QHD43...  0.855624  #a60f15  "
      ]
     },
     "execution_count": 29,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "variants.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Get 3D structure for selected protein"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [],
   "source": [
    "# some Electron microscopy structures are of low resolution. Keep only high and medium resolution structures (0 - 4 Å).\n",
    "resolution_threshold = 3.9"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [],
   "source": [
    "query = \"\"\"\n",
    "MATCH (g:Gene)-[:ENCODES]->(p:Protein{name: $protein_name, taxonomyId: $pathogen_taxonomy_id})-[h:HAS_TERTIARY_STRUCTURE]->(c:Chain)-[:IS_PART_OF_STRUCTURE]->(s:Structure)\n",
    "WHERE s.resolution <= $resolution_threshold\n",
    "RETURN p.name AS name, p.start, p.end, c.name AS structureChainId, c.uniprotStart, c.uniprotEnd, c.pdbStart, c.pdbEnd, s.resolution AS resolution, s.description AS description, h.coverage AS coverage\n",
    "ORDER BY resolution, coverage DESC\n",
    "\"\"\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [],
   "source": [
    "chains = graph.run(query, protein_name=protein_name, pathogen_taxonomy_id=pathogen_taxonomy_id, \n",
    "                   resolution_threshold=resolution_threshold).to_data_frame()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [],
   "source": [
    "chains['structureChainId'] = chains['structureChainId'].str[4:]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [],
   "source": [
    "chains['structureId'] = chains['structureChainId'].str[:4]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "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>c.pdbEnd</th>\n",
       "      <th>c.pdbStart</th>\n",
       "      <th>c.uniprotEnd</th>\n",
       "      <th>c.uniprotStart</th>\n",
       "      <th>coverage</th>\n",
       "      <th>description</th>\n",
       "      <th>name</th>\n",
       "      <th>p.end</th>\n",
       "      <th>p.start</th>\n",
       "      <th>resolution</th>\n",
       "      <th>structureChainId</th>\n",
       "      <th>structureId</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>843</th>\n",
       "      <td>[526]</td>\n",
       "      <td>[333]</td>\n",
       "      <td>[526]</td>\n",
       "      <td>[333]</td>\n",
       "      <td>0.152396</td>\n",
       "      <td>REGN10933 antibody Fab fragment heavy chain, R...</td>\n",
       "      <td>Spike glycoprotein</td>\n",
       "      <td>1273</td>\n",
       "      <td>1</td>\n",
       "      <td>3.9</td>\n",
       "      <td>6XDG.E</td>\n",
       "      <td>6XDG</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "    c.pdbEnd c.pdbStart c.uniprotEnd c.uniprotStart  coverage  \\\n",
       "843    [526]      [333]        [526]          [333]  0.152396   \n",
       "\n",
       "                                           description                name  \\\n",
       "843  REGN10933 antibody Fab fragment heavy chain, R...  Spike glycoprotein   \n",
       "\n",
       "     p.end  p.start  resolution structureChainId structureId  \n",
       "843   1273        1         3.9           6XDG.E        6XDG  "
      ]
     },
     "execution_count": 35,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "chains.query('structureId == \"6XDG\"')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {},
   "outputs": [],
   "source": [
    "if with_antibodies:\n",
    "    chains = chains[chains['description'].str.contains('heavy') | \n",
    "                    chains['description'].str.contains('antibody') |\n",
    "                    chains['description'].str.contains('nanobody') |\n",
    "                    chains['description'].str.contains('Fab')]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "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>c.pdbEnd</th>\n",
       "      <th>c.pdbStart</th>\n",
       "      <th>c.uniprotEnd</th>\n",
       "      <th>c.uniprotStart</th>\n",
       "      <th>coverage</th>\n",
       "      <th>description</th>\n",
       "      <th>name</th>\n",
       "      <th>p.end</th>\n",
       "      <th>p.start</th>\n",
       "      <th>resolution</th>\n",
       "      <th>structureChainId</th>\n",
       "      <th>structureId</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>843</th>\n",
       "      <td>[526]</td>\n",
       "      <td>[333]</td>\n",
       "      <td>[526]</td>\n",
       "      <td>[333]</td>\n",
       "      <td>0.152396</td>\n",
       "      <td>REGN10933 antibody Fab fragment heavy chain, R...</td>\n",
       "      <td>Spike glycoprotein</td>\n",
       "      <td>1273</td>\n",
       "      <td>1</td>\n",
       "      <td>3.9</td>\n",
       "      <td>6XDG.E</td>\n",
       "      <td>6XDG</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "    c.pdbEnd c.pdbStart c.uniprotEnd c.uniprotStart  coverage  \\\n",
       "843    [526]      [333]        [526]          [333]  0.152396   \n",
       "\n",
       "                                           description                name  \\\n",
       "843  REGN10933 antibody Fab fragment heavy chain, R...  Spike glycoprotein   \n",
       "\n",
       "     p.end  p.start  resolution structureChainId structureId  \n",
       "843   1273        1         3.9           6XDG.E        6XDG  "
      ]
     },
     "execution_count": 37,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "chains.query('structureId == \"6XDG\"')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "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>c.pdbEnd</th>\n",
       "      <th>c.pdbStart</th>\n",
       "      <th>c.uniprotEnd</th>\n",
       "      <th>c.uniprotStart</th>\n",
       "      <th>coverage</th>\n",
       "      <th>description</th>\n",
       "      <th>name</th>\n",
       "      <th>p.end</th>\n",
       "      <th>p.start</th>\n",
       "      <th>resolution</th>\n",
       "      <th>structureChainId</th>\n",
       "      <th>structureId</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>[527]</td>\n",
       "      <td>[324]</td>\n",
       "      <td>[527]</td>\n",
       "      <td>[324]</td>\n",
       "      <td>0.160251</td>\n",
       "      <td>2-acetamido-2-deoxy-beta-D-glucopyranose, 2-ac...</td>\n",
       "      <td>Spike glycoprotein</td>\n",
       "      <td>1273</td>\n",
       "      <td>1</td>\n",
       "      <td>1.4</td>\n",
       "      <td>7EAM.B</td>\n",
       "      <td>7EAM</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>[527]</td>\n",
       "      <td>[324]</td>\n",
       "      <td>[527]</td>\n",
       "      <td>[324]</td>\n",
       "      <td>0.160251</td>\n",
       "      <td>2-acetamido-2-deoxy-beta-D-glucopyranose, 2-ac...</td>\n",
       "      <td>Spike glycoprotein</td>\n",
       "      <td>1273</td>\n",
       "      <td>1</td>\n",
       "      <td>1.4</td>\n",
       "      <td>7EAM.A</td>\n",
       "      <td>7EAM</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>[1156]</td>\n",
       "      <td>[1147]</td>\n",
       "      <td>[1156]</td>\n",
       "      <td>[1147]</td>\n",
       "      <td>0.007855</td>\n",
       "      <td>B6 antigen-binding (Fab) fragment heavy chain,...</td>\n",
       "      <td>Spike glycoprotein</td>\n",
       "      <td>1273</td>\n",
       "      <td>1</td>\n",
       "      <td>1.4</td>\n",
       "      <td>7M53.A</td>\n",
       "      <td>7M53</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>[528]</td>\n",
       "      <td>[332]</td>\n",
       "      <td>[528]</td>\n",
       "      <td>[332]</td>\n",
       "      <td>0.154753</td>\n",
       "      <td>1,2-ETHANEDIOL, GLYCEROL, Spike protein S1, Sy...</td>\n",
       "      <td>Spike glycoprotein</td>\n",
       "      <td>1273</td>\n",
       "      <td>1</td>\n",
       "      <td>1.7</td>\n",
       "      <td>7MFU.A</td>\n",
       "      <td>7MFU</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td>[528]</td>\n",
       "      <td>[332]</td>\n",
       "      <td>[528]</td>\n",
       "      <td>[332]</td>\n",
       "      <td>0.154753</td>\n",
       "      <td>1,2-ETHANEDIOL, GLYCEROL, Spike protein S1, Sy...</td>\n",
       "      <td>Spike glycoprotein</td>\n",
       "      <td>1273</td>\n",
       "      <td>1</td>\n",
       "      <td>1.7</td>\n",
       "      <td>7MFU.D</td>\n",
       "      <td>7MFU</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "  c.pdbEnd c.pdbStart c.uniprotEnd c.uniprotStart  coverage  \\\n",
       "0    [527]      [324]        [527]          [324]  0.160251   \n",
       "1    [527]      [324]        [527]          [324]  0.160251   \n",
       "2   [1156]     [1147]       [1156]         [1147]  0.007855   \n",
       "4    [528]      [332]        [528]          [332]  0.154753   \n",
       "5    [528]      [332]        [528]          [332]  0.154753   \n",
       "\n",
       "                                         description                name  \\\n",
       "0  2-acetamido-2-deoxy-beta-D-glucopyranose, 2-ac...  Spike glycoprotein   \n",
       "1  2-acetamido-2-deoxy-beta-D-glucopyranose, 2-ac...  Spike glycoprotein   \n",
       "2  B6 antigen-binding (Fab) fragment heavy chain,...  Spike glycoprotein   \n",
       "4  1,2-ETHANEDIOL, GLYCEROL, Spike protein S1, Sy...  Spike glycoprotein   \n",
       "5  1,2-ETHANEDIOL, GLYCEROL, Spike protein S1, Sy...  Spike glycoprotein   \n",
       "\n",
       "   p.end  p.start  resolution structureChainId structureId  \n",
       "0   1273        1         1.4           7EAM.B        7EAM  \n",
       "1   1273        1         1.4           7EAM.A        7EAM  \n",
       "2   1273        1         1.4           7M53.A        7M53  \n",
       "4   1273        1         1.7           7MFU.A        7MFU  \n",
       "5   1273        1         1.7           7MFU.D        7MFU  "
      ]
     },
     "execution_count": 38,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "chains.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {},
   "outputs": [],
   "source": [
    "chains.drop_duplicates(subset=['structureChainId'], inplace=True)\n",
    "chains.sort_values(by=['structureChainId'], inplace=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Map uniprot residue numbers to PDB residue numbers"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {
    "jupyter": {
     "source_hidden": true
    }
   },
   "outputs": [],
   "source": [
    "def uniprot_to_pdb_mapping(row):\n",
    "    mapping = dict()\n",
    "    for (us,ue, ps, pe) in zip(row['c.uniprotStart'], row['c.uniprotEnd'], row['c.pdbStart'], row['c.pdbEnd']):\n",
    "        ps = int(ps)\n",
    "        pe = int(pe)\n",
    "        if (ue-us != pe-ps):\n",
    "            print('length mismatch:', row['structureChainId'], ue-us, pe-ps)\n",
    "        else:\n",
    "            offset = ps - us\n",
    "            for v in range(us, ue+1):\n",
    "                mapping[v] = offset + v\n",
    "                \n",
    "    #print(mapping)\n",
    "    return mapping"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {},
   "outputs": [],
   "source": [
    "chains['mapping'] = chains.apply(lambda row: uniprot_to_pdb_mapping(row), axis=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "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>c.pdbEnd</th>\n",
       "      <th>c.pdbStart</th>\n",
       "      <th>c.uniprotEnd</th>\n",
       "      <th>c.uniprotStart</th>\n",
       "      <th>coverage</th>\n",
       "      <th>description</th>\n",
       "      <th>name</th>\n",
       "      <th>p.end</th>\n",
       "      <th>p.start</th>\n",
       "      <th>resolution</th>\n",
       "      <th>structureChainId</th>\n",
       "      <th>structureId</th>\n",
       "      <th>mapping</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>843</th>\n",
       "      <td>[526]</td>\n",
       "      <td>[333]</td>\n",
       "      <td>[526]</td>\n",
       "      <td>[333]</td>\n",
       "      <td>0.152396</td>\n",
       "      <td>REGN10933 antibody Fab fragment heavy chain, R...</td>\n",
       "      <td>Spike glycoprotein</td>\n",
       "      <td>1273</td>\n",
       "      <td>1</td>\n",
       "      <td>3.9</td>\n",
       "      <td>6XDG.E</td>\n",
       "      <td>6XDG</td>\n",
       "      <td>{333: 333, 334: 334, 335: 335, 336: 336, 337: ...</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "    c.pdbEnd c.pdbStart c.uniprotEnd c.uniprotStart  coverage  \\\n",
       "843    [526]      [333]        [526]          [333]  0.152396   \n",
       "\n",
       "                                           description                name  \\\n",
       "843  REGN10933 antibody Fab fragment heavy chain, R...  Spike glycoprotein   \n",
       "\n",
       "     p.end  p.start  resolution structureChainId structureId  \\\n",
       "843   1273        1         3.9           6XDG.E        6XDG   \n",
       "\n",
       "                                               mapping  \n",
       "843  {333: 333, 334: 334, 335: 335, 336: 336, 337: ...  "
      ]
     },
     "execution_count": 42,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "chains.query('structureId == \"6XDG\"')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Visualize mutation sites\n",
    "\n",
    "Mutations are mapped onto protein chains for available 3D protein structures.\n",
    "\n",
    "Display options:\n",
    "\n",
    "|||\n",
    "|:-|:-|\n",
    "| *show_bio_assembly* | Toggle display of the biologically relevant quaternary structure |\n",
    "| *show_surface* | Toggle surface for protein chain |\n",
    "| *show_short_label* | Toggle display of mutation information<br>{UniProtResidue}{aminoAcid1}>{aminoAcid2}(# observations)<br>Example: 501N>Y(350436)|\n",
    "| *show_long_label* | Toggle display of mutation information<br>{PDBId}.{chainId}.{PDBResidue}: {geneName}.p{UniProtResidue}{aminoAcid1}>{aminoAcid2}(# observations)<br>Example: 6XDG.E.501: S:p.501N>Y(350436) |\n",
    "| *size* | Change size of visualization |\n",
    "| *font* | Change font size of annotations |\n",
    "| *logFreq* | Change minimum threshold to display mutations based on normalized log of mutation frequency [0.0 - 1.0]|\n",
    "| *structure* | Move slider to browse through available structures |"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Example: Move the structure slider to PDB ID:6XDG to see how mutations (e.g., 501N>Y) effect the binding of the Regeneron antibodies "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {
    "jupyter": {
     "source_hidden": true
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "# Setup viewer\n",
    "def view_mutations(df, variants, *args):\n",
    "    chainIds = list(df['structureChainId'])\n",
    "\n",
    "    def view3d(show_bio_assembly, show_surface, show_short_label, show_long_label, size, font, logFreq, i): \n",
    "        pdb_id, chain_id = chainIds[i].split('.')\n",
    "        global viewer1\n",
    "        viewer1 = py3Dmol.view(query='pdb:' + pdb_id, options={'doAssembly': show_bio_assembly}, width=size, height=size)\n",
    "\n",
    "        # polymer style\n",
    "        viewer1.setStyle({'cartoon': {'colorscheme': 'chain', 'width': 0.6, 'opacity':0.8}})\n",
    "\n",
    "        # highlight chain of interest in blue\n",
    "        viewer1.setStyle({'chain': chain_id},{'cartoon': {'color': 'blue'}})\n",
    "        \n",
    "        # non-polymer style\n",
    "        viewer1.setStyle({'hetflag': True}, {'stick':{'radius': 0.3, 'singleBond': False, 'colorscheme': 'greenCarbon'}})\n",
    "        \n",
    "        mapping = df['mapping'].iloc[i]\n",
    "\n",
    "        for row in variants.itertuples():\n",
    "            # get PDB residue mapping from a UniProt residue number\n",
    "            res_num = mapping.get(row.residue, 0)\n",
    "            col = row.color\n",
    "            if res_num > 0 and row.scale > logFreq:\n",
    "                mut_res = {'resi': res_num, 'chain': chain_id}\n",
    "                viewer1.addStyle(mut_res, {'sphere':{'color':col, 'opacity': 1.0}}) \n",
    "\n",
    "                if show_short_label:\n",
    "                    label = row.annotation\n",
    "                if show_long_label:\n",
    "                    label = chainIds[i] + \".\" + str(res_num) + \": \" + row.variationId\n",
    "                if show_short_label or show_long_label:\n",
    "                    viewer1.addLabel(label, {'fontSize':font,'fontColor': 'black','backgroundColor':'ivory', 'opacity': 1.0}, {'resi': res_num, 'chain': chain_id})\n",
    "\n",
    "        description = df['description'].iloc[i]\n",
    "        resolution = df['resolution'].iloc[i]\n",
    "        coverage = df['coverage'].iloc[i]\n",
    "        name = df['name'].iloc[i]\n",
    "        \n",
    "        print(name)\n",
    "        print()\n",
    "        print(f'PDB Id: {pdb_id}, chain Id: {chain_id}, resolution: {resolution}, sequence coverage: {coverage:.2f}')\n",
    "        print(f'description: {description}')\n",
    "        \n",
    "        # print any specified additional columns from the dataframe\n",
    "        for a in args:\n",
    "            print(a + \": \" + df.iloc[i][a])\n",
    "\n",
    "        viewer1.zoomTo({'chain': chain_id})\n",
    "        viewer1.center({'chain': chain_id})\n",
    "        \n",
    "        if show_surface:\n",
    "             viewer1.addSurface(py3Dmol.SES,{'opacity':0.8,'color':'lightblue'},{'chain': chain_id})\n",
    "\n",
    "        return viewer1.show()\n",
    "       \n",
    "    f_widget = IntSlider(value=9, min=5, max=20, description='font size', continuous_update=False)\n",
    "    z_widget = IntSlider(value=750, min=500, max=1200, description='size', continuous_update=False)\n",
    "    s_widget = IntSlider(min=0, max=len(chainIds)-1, description='structure', continuous_update=False)\n",
    "    l_widget = FloatSlider(value=0.8, min=0, max=1, step=0.05, description='logFreq:', \n",
    "                           continuous_update=False, orientation='horizontal', readout=True, readout_format='.2f')\n",
    "    \n",
    "    \n",
    "    return interact(view3d, show_bio_assembly=False, show_surface=False, show_short_label=True, show_long_label=False, size=z_widget, font=f_widget, logFreq=l_widget, i=s_widget)\n",
    "\n",
    "def view_image1():\n",
    "    return viewer1.png()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "cfd649c435f84f8b8f8142ee059fcfd9",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "interactive(children=(Checkbox(value=False, description='show_bio_assembly'), Checkbox(value=False, descriptio…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "view_mutations(chains, variants);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "metadata": {
    "jupyter": {
     "source_hidden": true
    }
   },
   "outputs": [],
   "source": [
    "# https://stackoverflow.com/questions/32468402/how-to-explode-a-list-inside-a-dataframe-cell-into-separate-rows\n",
    "import copy\n",
    "\n",
    "def pandas_explode(df, column_to_explode):\n",
    "    \"\"\"\n",
    "    Similar to Hive's EXPLODE function, take a column with iterable elements, and flatten the iterable to one element \n",
    "    per observation in the output table\n",
    "\n",
    "    :param df: A dataframe to explod\n",
    "    :type df: pandas.DataFrame\n",
    "    :param column_to_explode: \n",
    "    :type column_to_explode: str\n",
    "    :return: An exploded data frame\n",
    "    :rtype: pandas.DataFrame\n",
    "    \"\"\"\n",
    "\n",
    "    # Create a list of new observations\n",
    "    new_observations = list()\n",
    "\n",
    "    # Iterate through existing observations\n",
    "    for row in df.to_dict(orient='records'):\n",
    "\n",
    "        # Take out the exploding iterable\n",
    "        explode_values = row[column_to_explode]\n",
    "        del row[column_to_explode]\n",
    "\n",
    "        # Create a new observation for every entry in the exploding iterable & add all of the other columns\n",
    "        for explode_value in explode_values.items():\n",
    "\n",
    "            # Deep copy existing observation\n",
    "            new_observation = copy.deepcopy(row)\n",
    "\n",
    "            # Add one (newly flattened) value from exploding iterable\n",
    "            new_observation[column_to_explode] = explode_value\n",
    "\n",
    "            # Add to the list of new observations\n",
    "            new_observations.append(new_observation)\n",
    "\n",
    "    # Create a DataFrame\n",
    "    return_df = pd.DataFrame(new_observations)\n",
    "\n",
    "    # Return\n",
    "    return return_df"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Expand chains into residues"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "metadata": {},
   "outputs": [],
   "source": [
    "residues = pandas_explode(chains, 'mapping')\n",
    "residues['uniprotPosition'] = residues['mapping'].apply(lambda x: x[0])\n",
    "residues['pdbPosition'] = residues['mapping'].apply(lambda x: x[1])\n",
    "residues.drop(columns='mapping', inplace=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 47,
   "metadata": {},
   "outputs": [],
   "source": [
    "residues = residues.drop(columns=['c.pdbEnd', 'c.pdbStart', 'c.uniprotEnd', 'c.uniprotStart', 'p.end', 'p.start'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 48,
   "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>coverage</th>\n",
       "      <th>description</th>\n",
       "      <th>name</th>\n",
       "      <th>resolution</th>\n",
       "      <th>structureChainId</th>\n",
       "      <th>structureId</th>\n",
       "      <th>uniprotPosition</th>\n",
       "      <th>pdbPosition</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>0.153181</td>\n",
       "      <td>2-acetamido-2-deoxy-beta-D-glucopyranose, CR30...</td>\n",
       "      <td>Spike glycoprotein</td>\n",
       "      <td>3.084</td>\n",
       "      <td>6W41.C</td>\n",
       "      <td>6W41</td>\n",
       "      <td>333</td>\n",
       "      <td>333</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>0.153181</td>\n",
       "      <td>2-acetamido-2-deoxy-beta-D-glucopyranose, CR30...</td>\n",
       "      <td>Spike glycoprotein</td>\n",
       "      <td>3.084</td>\n",
       "      <td>6W41.C</td>\n",
       "      <td>6W41</td>\n",
       "      <td>334</td>\n",
       "      <td>334</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>0.153181</td>\n",
       "      <td>2-acetamido-2-deoxy-beta-D-glucopyranose, CR30...</td>\n",
       "      <td>Spike glycoprotein</td>\n",
       "      <td>3.084</td>\n",
       "      <td>6W41.C</td>\n",
       "      <td>6W41</td>\n",
       "      <td>335</td>\n",
       "      <td>335</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>0.153181</td>\n",
       "      <td>2-acetamido-2-deoxy-beta-D-glucopyranose, CR30...</td>\n",
       "      <td>Spike glycoprotein</td>\n",
       "      <td>3.084</td>\n",
       "      <td>6W41.C</td>\n",
       "      <td>6W41</td>\n",
       "      <td>336</td>\n",
       "      <td>336</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>0.153181</td>\n",
       "      <td>2-acetamido-2-deoxy-beta-D-glucopyranose, CR30...</td>\n",
       "      <td>Spike glycoprotein</td>\n",
       "      <td>3.084</td>\n",
       "      <td>6W41.C</td>\n",
       "      <td>6W41</td>\n",
       "      <td>337</td>\n",
       "      <td>337</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   coverage                                        description  \\\n",
       "0  0.153181  2-acetamido-2-deoxy-beta-D-glucopyranose, CR30...   \n",
       "1  0.153181  2-acetamido-2-deoxy-beta-D-glucopyranose, CR30...   \n",
       "2  0.153181  2-acetamido-2-deoxy-beta-D-glucopyranose, CR30...   \n",
       "3  0.153181  2-acetamido-2-deoxy-beta-D-glucopyranose, CR30...   \n",
       "4  0.153181  2-acetamido-2-deoxy-beta-D-glucopyranose, CR30...   \n",
       "\n",
       "                 name  resolution structureChainId structureId  \\\n",
       "0  Spike glycoprotein       3.084           6W41.C        6W41   \n",
       "1  Spike glycoprotein       3.084           6W41.C        6W41   \n",
       "2  Spike glycoprotein       3.084           6W41.C        6W41   \n",
       "3  Spike glycoprotein       3.084           6W41.C        6W41   \n",
       "4  Spike glycoprotein       3.084           6W41.C        6W41   \n",
       "\n",
       "   uniprotPosition  pdbPosition  \n",
       "0              333          333  \n",
       "1              334          334  \n",
       "2              335          335  \n",
       "3              336          336  \n",
       "4              337          337  "
      ]
     },
     "execution_count": 48,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "residues.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 49,
   "metadata": {},
   "outputs": [],
   "source": [
    "variants = variants[['residue', 'variationId', 'annotation', 'scale', 'color']]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 50,
   "metadata": {},
   "outputs": [],
   "source": [
    "residues_variants = residues.merge(variants, left_on='uniprotPosition', right_on='residue')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 51,
   "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>coverage</th>\n",
       "      <th>description</th>\n",
       "      <th>name</th>\n",
       "      <th>resolution</th>\n",
       "      <th>structureChainId</th>\n",
       "      <th>structureId</th>\n",
       "      <th>uniprotPosition</th>\n",
       "      <th>pdbPosition</th>\n",
       "      <th>residue</th>\n",
       "      <th>variationId</th>\n",
       "      <th>annotation</th>\n",
       "      <th>scale</th>\n",
       "      <th>color</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>0.153181</td>\n",
       "      <td>2-acetamido-2-deoxy-beta-D-glucopyranose, CR30...</td>\n",
       "      <td>Spike glycoprotein</td>\n",
       "      <td>3.084</td>\n",
       "      <td>6W41.C</td>\n",
       "      <td>6W41</td>\n",
       "      <td>333</td>\n",
       "      <td>333</td>\n",
       "      <td>333</td>\n",
       "      <td>S:p.333T&gt;I(8), S:p.333T&gt;K(4), S:p.333T&gt;A(2)</td>\n",
       "      <td>333T&gt;I(8), 333T&gt;K(4), 333T&gt;A(2)</td>\n",
       "      <td>0.163116</td>\n",
       "      <td>#fdd5c4</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>0.750196</td>\n",
       "      <td>2-acetamido-2-deoxy-beta-D-glucopyranose, 2-ac...</td>\n",
       "      <td>Spike glycoprotein</td>\n",
       "      <td>3.100</td>\n",
       "      <td>6WPS.A</td>\n",
       "      <td>6WPS</td>\n",
       "      <td>333</td>\n",
       "      <td>333</td>\n",
       "      <td>333</td>\n",
       "      <td>S:p.333T&gt;I(8), S:p.333T&gt;K(4), S:p.333T&gt;A(2)</td>\n",
       "      <td>333T&gt;I(8), 333T&gt;K(4), 333T&gt;A(2)</td>\n",
       "      <td>0.163116</td>\n",
       "      <td>#fdd5c4</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>0.750196</td>\n",
       "      <td>2-acetamido-2-deoxy-beta-D-glucopyranose, 2-ac...</td>\n",
       "      <td>Spike glycoprotein</td>\n",
       "      <td>3.100</td>\n",
       "      <td>6WPS.B</td>\n",
       "      <td>6WPS</td>\n",
       "      <td>333</td>\n",
       "      <td>333</td>\n",
       "      <td>333</td>\n",
       "      <td>S:p.333T&gt;I(8), S:p.333T&gt;K(4), S:p.333T&gt;A(2)</td>\n",
       "      <td>333T&gt;I(8), 333T&gt;K(4), 333T&gt;A(2)</td>\n",
       "      <td>0.163116</td>\n",
       "      <td>#fdd5c4</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>0.750196</td>\n",
       "      <td>2-acetamido-2-deoxy-beta-D-glucopyranose, 2-ac...</td>\n",
       "      <td>Spike glycoprotein</td>\n",
       "      <td>3.100</td>\n",
       "      <td>6WPS.E</td>\n",
       "      <td>6WPS</td>\n",
       "      <td>333</td>\n",
       "      <td>333</td>\n",
       "      <td>333</td>\n",
       "      <td>S:p.333T&gt;I(8), S:p.333T&gt;K(4), S:p.333T&gt;A(2)</td>\n",
       "      <td>333T&gt;I(8), 333T&gt;K(4), 333T&gt;A(2)</td>\n",
       "      <td>0.163116</td>\n",
       "      <td>#fdd5c4</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>0.742341</td>\n",
       "      <td>2-acetamido-2-deoxy-beta-D-glucopyranose, 2-ac...</td>\n",
       "      <td>Spike glycoprotein</td>\n",
       "      <td>3.700</td>\n",
       "      <td>6WPT.A</td>\n",
       "      <td>6WPT</td>\n",
       "      <td>333</td>\n",
       "      <td>333</td>\n",
       "      <td>333</td>\n",
       "      <td>S:p.333T&gt;I(8), S:p.333T&gt;K(4), S:p.333T&gt;A(2)</td>\n",
       "      <td>333T&gt;I(8), 333T&gt;K(4), 333T&gt;A(2)</td>\n",
       "      <td>0.163116</td>\n",
       "      <td>#fdd5c4</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   coverage                                        description  \\\n",
       "0  0.153181  2-acetamido-2-deoxy-beta-D-glucopyranose, CR30...   \n",
       "1  0.750196  2-acetamido-2-deoxy-beta-D-glucopyranose, 2-ac...   \n",
       "2  0.750196  2-acetamido-2-deoxy-beta-D-glucopyranose, 2-ac...   \n",
       "3  0.750196  2-acetamido-2-deoxy-beta-D-glucopyranose, 2-ac...   \n",
       "4  0.742341  2-acetamido-2-deoxy-beta-D-glucopyranose, 2-ac...   \n",
       "\n",
       "                 name  resolution structureChainId structureId  \\\n",
       "0  Spike glycoprotein       3.084           6W41.C        6W41   \n",
       "1  Spike glycoprotein       3.100           6WPS.A        6WPS   \n",
       "2  Spike glycoprotein       3.100           6WPS.B        6WPS   \n",
       "3  Spike glycoprotein       3.100           6WPS.E        6WPS   \n",
       "4  Spike glycoprotein       3.700           6WPT.A        6WPT   \n",
       "\n",
       "   uniprotPosition  pdbPosition  residue  \\\n",
       "0              333          333      333   \n",
       "1              333          333      333   \n",
       "2              333          333      333   \n",
       "3              333          333      333   \n",
       "4              333          333      333   \n",
       "\n",
       "                                   variationId  \\\n",
       "0  S:p.333T>I(8), S:p.333T>K(4), S:p.333T>A(2)   \n",
       "1  S:p.333T>I(8), S:p.333T>K(4), S:p.333T>A(2)   \n",
       "2  S:p.333T>I(8), S:p.333T>K(4), S:p.333T>A(2)   \n",
       "3  S:p.333T>I(8), S:p.333T>K(4), S:p.333T>A(2)   \n",
       "4  S:p.333T>I(8), S:p.333T>K(4), S:p.333T>A(2)   \n",
       "\n",
       "                        annotation     scale    color  \n",
       "0  333T>I(8), 333T>K(4), 333T>A(2)  0.163116  #fdd5c4  \n",
       "1  333T>I(8), 333T>K(4), 333T>A(2)  0.163116  #fdd5c4  \n",
       "2  333T>I(8), 333T>K(4), 333T>A(2)  0.163116  #fdd5c4  \n",
       "3  333T>I(8), 333T>K(4), 333T>A(2)  0.163116  #fdd5c4  \n",
       "4  333T>I(8), 333T>K(4), 333T>A(2)  0.163116  #fdd5c4  "
      ]
     },
     "execution_count": 51,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "residues_variants.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 52,
   "metadata": {},
   "outputs": [],
   "source": [
    "residues_variants.to_csv(output_file_name, index=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Now run the next step\n",
    "Map mutations occuring at protein-protein interaction sites: [2-MapToPolymerInteractions.ipynb](2-MapToPolymerInteractions.ipynb)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.7.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}