COVID_19_PANGO_Multi_Input_Variant_Conversion_Tool.ipynb

{
  "nbformat": 4,
  "nbformat_minor": 0,
  "metadata": {
    "colab": {
      "name": "COVID-19 PANGO: Multi-Input Variant Conversion Tool.ipynb",
      "provenance": [],
      "collapsed_sections": [],
      "include_colab_link": true
    },
    "kernelspec": {
      "name": "python3",
      "display_name": "Python 3"
    },
    "language_info": {
      "name": "python"
    }
  },
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "view-in-github",
        "colab_type": "text"
      },
      "source": [
        "<a href=\"https://colab.research.google.com/github/Karine-Moussa/PANGO-Genomic-Conversions/blob/main/COVID_19_PANGO_Multi_Input_Variant_Conversion_Tool.ipynb\" target=\"_parent\"><img src=\"https://colab.research.google.com/assets/colab-badge.svg\" alt=\"Open In Colab\"/></a>"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "7fbCNCGCftFP"
      },
      "source": [
        "# PANGO Mutation Nomenclature: Convert Genetic Location to Genomic Location\n",
        "\n",
        "\n",
        "---\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "i2gZcKJvbSJn"
      },
      "source": [
        "# Instructions\n",
        "\n",
        "1. In colab menu: *Runtime > Runall*\n",
        "\n",
        "2. In the **INPUT** section enter one or more mutations in PANGO format, then hit `ENTER`\n",
        "\n",
        " - `[GENE]:[REF_AMINO_ACID][AMINO_ACID_LOC][ALT_AMINO_ACID]`\n",
        "    - Example: `S:D1118H`\n",
        " -  **Note:** do not include mutation type\n",
        "     - ie: ~~`aa:`~~`S:D1118H`\n",
        " - For multiple mutations, enter mutations in **comma seperated** format\n",
        "\n",
        "---\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "wcRZXUaXfaQU"
      },
      "source": [
        "## INPUT"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "ajAHx0L9MbCf",
        "outputId": "017e0411-85b4-4d6b-c5ce-5917db7cedf1"
      },
      "source": [
        "# Input mutations of interest\n",
        "import re\n",
        "while True:\n",
        "  mut_input = input()\n",
        "  break\n",
        "# Example: S:D1118H, N:P80R, orf1ab:S1188L "
      ],
      "execution_count": null,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "E:I82T\n"
          ],
          "name": "stdout"
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "XWyY-_pQiigw"
      },
      "source": [
        "# MAIN"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "Wc75L6-szZj5"
      },
      "source": [
        "## Check Input Format"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "leqG5FZ4bdN2"
      },
      "source": [
        "# Check inputted mutations. Only proceed with those in the correct format.\n",
        "import re\n",
        "mutations = []\n",
        "mut_input = mut_input.replace(\" \", \"\").split(\",\")\n",
        "for mut in mut_input:\n",
        " # Verify nomenclature\n",
        "  if re.match(\"([a-z 0-9]+):([a-z]+)([0-9]+)([a-z\\*]+)\", mut, flags=re.IGNORECASE):\n",
        "    mutations.append(mut)\n",
        "  else:\n",
        "    print(f\"'%s' invalid format\" % mut)"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "HDeGC55rJ45R"
      },
      "source": [
        "## Set up"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "VmMB3289EnI1",
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "outputId": "f4b27988-62ce-4441-b9d6-aa31d0afd670"
      },
      "source": [
        "# Install python packages\n",
        "!pip -q install pyfaidx # -q suppresses pip output message"
      ],
      "execution_count": null,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "  Building wheel for pyfaidx (setup.py) ... \u001b[?25l\u001b[?25hdone\n"
          ],
          "name": "stdout"
        }
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "VMxI6FzPECyt"
      },
      "source": [
        "# Import python libraries\n",
        "import os.path\n",
        "from os import path\n",
        "import pandas as pd\n",
        "from pyfaidx import Fasta"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "Dr6GQ0pfEwXf",
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "outputId": "fcb3d559-58db-4e4d-8b14-e3162c42124a"
      },
      "source": [
        "# Get Genome\n",
        "if path.exists(\"GCF_009858895.2_ASM985889v3_genomic.fna\") == False:\n",
        "  !wget https://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/009/858/895/GCF_009858895.2_ASM985889v3/GCF_009858895.2_ASM985889v3_genomic.fna.gz\n",
        "  !gunzip GCF_009858895.2_ASM985889v3_genomic.fna.gz"
      ],
      "execution_count": null,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "--2021-05-06 20:38:20--  https://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/009/858/895/GCF_009858895.2_ASM985889v3/GCF_009858895.2_ASM985889v3_genomic.fna.gz\n",
            "Resolving ftp.ncbi.nlm.nih.gov (ftp.ncbi.nlm.nih.gov)... 130.14.250.12, 130.14.250.10, 2607:f220:41e:250::11, ...\n",
            "Connecting to ftp.ncbi.nlm.nih.gov (ftp.ncbi.nlm.nih.gov)|130.14.250.12|:443... connected.\n",
            "HTTP request sent, awaiting response... 200 OK\n",
            "Length: 9591 (9.4K) [application/x-gzip]\n",
            "Saving to: ‘GCF_009858895.2_ASM985889v3_genomic.fna.gz’\n",
            "\n",
            "\r          GCF_00985   0%[                    ]       0  --.-KB/s               \rGCF_009858895.2_ASM 100%[===================>]   9.37K  --.-KB/s    in 0s      \n",
            "\n",
            "2021-05-06 20:38:20 (157 MB/s) - ‘GCF_009858895.2_ASM985889v3_genomic.fna.gz’ saved [9591/9591]\n",
            "\n"
          ],
          "name": "stdout"
        }
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "4HwEblcvFiYx",
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "outputId": "56d56a59-a03d-4d28-a1de-a3572da1d164"
      },
      "source": [
        "# Get genome gff \n",
        "if path.exists(\"GCF_009858895.2_ASM985889v3_genomic.gff\") == False:\n",
        "  !wget https://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/009/858/895/GCF_009858895.2_ASM985889v3/GCF_009858895.2_ASM985889v3_genomic.gff.gz\n",
        "  !gunzip GCF_009858895.2_ASM985889v3_genomic.gff"
      ],
      "execution_count": null,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "--2021-05-06 20:38:21--  https://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/009/858/895/GCF_009858895.2_ASM985889v3/GCF_009858895.2_ASM985889v3_genomic.gff.gz\n",
            "Resolving ftp.ncbi.nlm.nih.gov (ftp.ncbi.nlm.nih.gov)... 130.14.250.11, 130.14.250.13, 2607:f220:41e:250::11, ...\n",
            "Connecting to ftp.ncbi.nlm.nih.gov (ftp.ncbi.nlm.nih.gov)|130.14.250.11|:443... connected.\n",
            "HTTP request sent, awaiting response... 200 OK\n",
            "Length: 2914 (2.8K) [application/x-gzip]\n",
            "Saving to: ‘GCF_009858895.2_ASM985889v3_genomic.gff.gz’\n",
            "\n",
            "\r          GCF_00985   0%[                    ]       0  --.-KB/s               \rGCF_009858895.2_ASM 100%[===================>]   2.85K  --.-KB/s    in 0s      \n",
            "\n",
            "2021-05-06 20:38:21 (121 MB/s) - ‘GCF_009858895.2_ASM985889v3_genomic.gff.gz’ saved [2914/2914]\n",
            "\n"
          ],
          "name": "stdout"
        }
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "TKjhJKcVa6Co"
      },
      "source": [
        "# Load genome, grab DNA like genome[\"NC_045512.2\"][start:end]\n",
        "genome = Fasta('GCF_009858895.2_ASM985889v3_genomic.fna')"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "PADCqsATfWuM"
      },
      "source": [
        "# Load gff as pandas data frame\n",
        "df = pd.read_csv(\"GCF_009858895.2_ASM985889v3_genomic.gff\", comment=\"#\", sep=\"\\t\", header=None)\n",
        "\n",
        "# Collect GFF information on all genes\n",
        "df = df[df[2] == \"gene\"]\n",
        "df[8] = df[8].apply(lambda line: [n.lstrip(\"gene=\").lower() for n in line.split(\";\") if n.startswith(\"gene=\")][0])\n",
        "# Manually add a row for orf1a (266..13483 taken from https://www.ncbi.nlm.nih.gov/nuccore/NC_045512.2) \n",
        "orf1a = pd.DataFrame(pd.DataFrame([['NC_045512.2', 'RefSeq','gene',266,13483,'.','+','.','orf1a']],)) \n",
        "df = df.append(orf1a, ignore_index=True)\n",
        "genes_list = list(df[8])"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "mHmxFkXjJw6J"
      },
      "source": [
        "## Conversion"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "GmWDai-MEJRZ"
      },
      "source": [
        "codon_table = {\"TTT\": \"F\", \"TTC\": \"F\", \"TTA\": \"L\", \"TTG\": \"L\",\n",
        "               \"CTT\": \"L\", \"CTC\": \"L\", \"CTA\": \"L\", \"CTG\": \"L\",\n",
        "               \"ATT\": \"I\", \"ATC\": \"I\", \"ATA\": \"I\", \"ATG\": \"M\",\n",
        "               \"GTT\": \"V\", \"GTC\": \"V\", \"GTA\": \"V\", \"GTG\": \"V\",\n",
        "               \"TCT\": \"S\", \"TCC\": \"S\", \"TCA\": \"S\", \"TCG\": \"S\",\n",
        "               \"CCT\": \"P\", \"CCC\": \"P\", \"CCA\": \"P\", \"CCG\": \"P\",\n",
        "               \"ACT\": \"T\", \"ACC\": \"T\", \"ACA\": \"T\", \"ACG\": \"T\",\n",
        "               \"GCT\": \"A\", \"GCC\": \"A\", \"GCA\": \"A\", \"GCG\": \"A\",\n",
        "               \"TAT\": \"Y\", \"TAC\": \"Y\", \"TAA\": \"*\", \"TAG\": \"*\",\n",
        "               \"CAT\": \"H\", \"CAC\": \"H\", \"CAA\": \"Q\", \"CAG\": \"Q\",\n",
        "               \"AAT\": \"N\", \"AAC\": \"N\", \"AAA\": \"K\", \"AAG\": \"K\",\n",
        "               \"GAT\": \"D\", \"GAC\": \"D\", \"GAA\": \"E\", \"GAG\": \"E\",\n",
        "               \"TGT\": \"C\", \"TGC\": \"C\", \"TGA\": \"*\", \"TGG\": \"W\",\n",
        "               \"CGT\": \"R\", \"CGC\": \"R\", \"CGA\": \"R\", \"CGG\": \"R\",\n",
        "               \"AGT\": \"S\", \"AGC\": \"S\", \"AGA\": \"R\", \"AGG\": \"R\",\n",
        "               \"GGT\": \"G\", \"GGC\": \"G\", \"GGA\": \"G\", \"GGG\": \"G\"}"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "Onn9KshSEMQ-"
      },
      "source": [
        "# Create dna -> aa translate function\n",
        "def translate(seq):\n",
        "    p = \"\"\n",
        "    for loc in range(0, len(seq), 3):\n",
        "        p += codon_table[seq[loc:loc+3]]\n",
        "    return p"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "T6a4NDS4XPkx"
      },
      "source": [
        "# Create mutant parsing function\n",
        "def mutant_parse(mutation):\n",
        "  gene_name, before, loc, after = re.match(\"([a-z 0-9]+):([a-z]+)([0-9]+)([a-z\\*]+)\", mutation, flags=re.IGNORECASE).groups()\n",
        "  loc = int(loc)\n",
        "  return gene_name, before, loc, after\n",
        "gene_name, before, loc, after = mutant_parse(\"S:D80L\")"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "dyd8zOnOoKKG"
      },
      "source": [
        "# Get gene start and end\n",
        "def get_gene_start(x):\n",
        "  if x.lower() in genes_list:\n",
        "    start = int(df[df[8]==x.lower()][3].values[0])\n",
        "  else:\n",
        "    start = int()\n",
        "  return start\n",
        "\n",
        "def get_gene_end(x):\n",
        "  if x.lower() in genes_list:\n",
        "    end = int(df[df[8]==x.lower()][4].values[0])\n",
        "  else:\n",
        "    end = int()\n",
        "  return end"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "DNYKr9oLK8ij"
      },
      "source": [
        "# Locate the amino acid in the genome\n",
        "def locate_aa(gene_data):\n",
        "  start, end = (gene_data[3]-1 + 3*(loc-1), gene_data[3]-1 + 3*(loc-1) + 3)\n",
        "  return(start, end)"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "MQo0s0gzm65B"
      },
      "source": [
        "# Create df with information on each mutation\n",
        "snpaa = pd.DataFrame()\n",
        "snpaa['mutation'] = mutations\n",
        "snpaa[['gene_name','before','loc','after']] = snpaa.apply(\n",
        "    lambda x: pd.Series(mutant_parse(x['mutation'])), axis=1)\n",
        "\n",
        "# add gene start and end to dataframe\n",
        "snpaa['gene_start'] = snpaa['gene_name'].apply(get_gene_start)\n",
        "snpaa['gene_end'] = snpaa['gene_name'].apply(get_gene_end)\n",
        "\n",
        "# add genomic start and end to dataframe\n",
        "snpaa['mut_start'] = (snpaa['gene_start']-1 + 3*(snpaa['loc']-1))\n",
        "snpaa['mut_end'] = (snpaa['gene_start']-1 + 3*(snpaa['loc']-1) + 3)\n",
        "\n",
        "# add translated sequence\n",
        "snpaa['trans_aa'] = snpaa.apply(\n",
        "    lambda x: translate(str(genome[\"NC_045512.2\"][x['mut_start']:x['mut_end']])), axis =1)\n",
        "\n",
        "\n",
        "# add verification column\n",
        "snpaa['pass'] = (snpaa['before'] == snpaa['trans_aa'])"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "j9YkYTNgvIso"
      },
      "source": [
        "# Save df to csv\n",
        "snpaa.to_csv('snpaa_subset.csv')"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "fTLPNjfZcyEe"
      },
      "source": [
        "# OUTPUT"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "kbl26t8-c2Xn",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 80
        },
        "outputId": "58483595-6648-43c2-9569-ac700b6efb3b"
      },
      "source": [
        "snpaa"
      ],
      "execution_count": null,
      "outputs": [
        {
          "output_type": "execute_result",
          "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>mutation</th>\n",
              "      <th>gene_name</th>\n",
              "      <th>before</th>\n",
              "      <th>loc</th>\n",
              "      <th>after</th>\n",
              "      <th>gene_start</th>\n",
              "      <th>gene_end</th>\n",
              "      <th>mut_start</th>\n",
              "      <th>mut_end</th>\n",
              "      <th>trans_aa</th>\n",
              "      <th>pass</th>\n",
              "    </tr>\n",
              "  </thead>\n",
              "  <tbody>\n",
              "    <tr>\n",
              "      <th>0</th>\n",
              "      <td>E:I82T</td>\n",
              "      <td>E</td>\n",
              "      <td>I</td>\n",
              "      <td>82</td>\n",
              "      <td>T</td>\n",
              "      <td>26245</td>\n",
              "      <td>26472</td>\n",
              "      <td>26487</td>\n",
              "      <td>26490</td>\n",
              "      <td>Y</td>\n",
              "      <td>False</td>\n",
              "    </tr>\n",
              "  </tbody>\n",
              "</table>\n",
              "</div>"
            ],
            "text/plain": [
              "  mutation gene_name before  loc  ... mut_start  mut_end  trans_aa   pass\n",
              "0   E:I82T         E      I   82  ...     26487    26490         Y  False\n",
              "\n",
              "[1 rows x 11 columns]"
            ]
          },
          "metadata": {
            "tags": []
          },
          "execution_count": 16
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "Mlo8jw2GvSKI"
      },
      "source": [
        "**Note: The mut_start is the 0 based coordinate of the first nucleotide in the codon**\n",
        "\n",
        "<br></br>\n",
        "\n",
        "To download output, go to the left tab and select the *files* icon > right-click on `snpaa_subset.csv` > select *download*\n",
        "\n",
        "<img width=\"200\" alt=\"image\" src=\"https://i.imgur.com/xAc9Eym.png\">"
      ]
    }
  ]
}