notebooks/analysis/mapping_analysis.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "bc372cc5",
   "metadata": {},
   "source": [
    "# MaveDB Mapping Method Analysis\n",
    "This notebook assesses the performance of the mapping method, comparing the match between reference sequences for pre-mapped and post-mapped VRS objects. This notebook also computes the number of unique VRS alleles generated across the examine score sets."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ac0dcb44",
   "metadata": {},
   "source": [
    "## Load Relevant Libraries\n",
    "Run the cell below to load the libraries used in the analysis"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "0e09bf08",
   "metadata": {},
   "outputs": [],
   "source": [
    "import json\n",
    "from Bio.Seq import Seq\n",
    "import pandas as pd"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7ec3d972",
   "metadata": {},
   "source": [
    "## Load List of Examined Score Sets"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "bb03ef26",
   "metadata": {},
   "outputs": [],
   "source": [
    "mave_metadata_dat = pd.read_csv(\"analysis_files/mave_dat.csv\", index_col=[0])\n",
    "score_sets = mave_metadata_dat[\"urn\"].to_list()[:-5]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4d037805",
   "metadata": {},
   "source": [
    "## Determine the MaveDB IDs of Variants with Reference Mismatches\n",
    "Run the cell below to generate a dictionary listing the MaveDB IDs of variants with reference mismatches. The dictionary is keyed by score set URN."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "3f5e44d2",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'The number of examined variant pairs is: 2499044'"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "diff_vars_dict = {}\n",
    "var_count = 0\n",
    "var_count_dict_new = {}\n",
    "\n",
    "for key in score_sets:\n",
    "    if key != \"urn:mavedb:00000072-a-1\" and key != \"urn:mavedb:00000105-a-1\": # No mapping for these score sets\n",
    "        f = open(f\"analysis_files/mappings/{key[11::]}.json\")\n",
    "        dat = json.load(f)\n",
    "        dat = dat[\"mapped_scores\"]\n",
    "\n",
    "        diff_vars = []\n",
    "\n",
    "        for j in range(len(dat)):\n",
    "            if \"members\" not in dat[j][\"pre_mapped\"]:\n",
    "                var_count += 1\n",
    "                seq_pre = dat[j][\"pre_mapped\"][\"vrs_ref_allele_seq\"]\n",
    "                seq_post = dat[j][\"post_mapped\"][\"vrs_ref_allele_seq\"]\n",
    "                seq_pre_rv = str(Seq(seq_pre).reverse_complement())\n",
    "\n",
    "                if seq_pre != seq_post and seq_post != seq_pre_rv:\n",
    "                    diff_vars.append(j)\n",
    "\n",
    "            else:\n",
    "                for k in range(len(dat[j][\"pre_mapped\"][\"members\"])):\n",
    "                    var_count += 1\n",
    "                    seq_pre = dat[j][\"pre_mapped\"][\"members\"][k][\"vrs_ref_allele_seq\"]\n",
    "                    seq_post = dat[j][\"post_mapped\"][\"members\"][k][\"vrs_ref_allele_seq\"]\n",
    "                    seq_pre_rv = str(Seq(seq_pre).reverse_complement())\n",
    "\n",
    "                    if seq_pre != seq_post and seq_post != seq_pre_rv:\n",
    "                        diff_vars.append(j)\n",
    "\n",
    "        diff_vars_dict[key] = diff_vars\n",
    "        var_count_dict_new[key] = var_count\n",
    "f\"The number of examined variant pairs is: {var_count}\""
   ]
  },
  {
   "cell_type": "markdown",
   "id": "57c4d6d7",
   "metadata": {},
   "source": [
    "### The cell below can be uncommented to examine the output of the above cell"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "240e02be",
   "metadata": {},
   "outputs": [],
   "source": [
    "#diff_vars_dict"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2e27dc3f",
   "metadata": {},
   "source": [
    "### Examine Example Mismatch\n",
    "Run the cell below to view an example of reference mismatch discordance due to a MAVE variant that spans an alignment block"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "1c5c6167",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{'pre_mapped': {'id': 'ga4gh:VA.mmZEpS0H-V-dGn4n8V1uPWy18MY18PrR',\n",
       "  'type': 'VariationDescriptor',\n",
       "  'variation': {'id': 'ga4gh:VA.mmZEpS0H-V-dGn4n8V1uPWy18MY18PrR',\n",
       "   'type': 'Allele',\n",
       "   'location': {'id': None,\n",
       "    'type': 'SequenceLocation',\n",
       "    'sequence_id': 'ga4gh:SQ.sVMC1jmTXRvuzBCDJ8aoBmZ_Uu35YFj7',\n",
       "    'interval': {'type': 'SequenceInterval',\n",
       "     'start': {'type': 'Number', 'value': 316},\n",
       "     'end': {'type': 'Number', 'value': 318}}},\n",
       "   'state': {'type': 'LiteralSequenceExpression', 'sequence': 'CT'}},\n",
       "  'vrs_ref_allele_seq': 'TG'},\n",
       " 'post_mapped': {'id': 'ga4gh:VA.th7X__sHHeffw9AGpRyCUfg4jN-hM2mE',\n",
       "  'type': 'VariationDescriptor',\n",
       "  'variation': {'id': 'ga4gh:VA.th7X__sHHeffw9AGpRyCUfg4jN-hM2mE',\n",
       "   'type': 'Allele',\n",
       "   'location': {'id': None,\n",
       "    'type': 'SequenceLocation',\n",
       "    'sequence_id': 'ga4gh:SQ.5ZUqxCmDDgN4xTRbaSjN8LwgZironmB8',\n",
       "    'interval': {'type': 'SequenceInterval',\n",
       "     'start': {'type': 'Number', 'value': 43068507},\n",
       "     'end': {'type': 'Number', 'value': 43068508}}},\n",
       "   'state': {'type': 'LiteralSequenceExpression', 'sequence': 'G'}},\n",
       "  'expressions': [{'type': 'Expression',\n",
       "    'syntax': 'hgvs.g',\n",
       "    'value': 'NC_000021.9:g.43068508C>G',\n",
       "    'syntax_version': None}],\n",
       "  'vrs_ref_allele_seq': 'C'},\n",
       " 'mavedb_id': 'urn:mavedb:00000005-a-6#5601',\n",
       " 'relation': 'SO:is_homologous_to',\n",
       " 'score': 0.0836691871353702}"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "f = open(f'analysis_files/mappings/00000005-a-6.json')\n",
    "dat = json.load(f)\n",
    "dat[\"mapped_scores\"][5601]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bbe07ed2",
   "metadata": {},
   "source": [
    "## Compute the Proportion of Reference Mismatches\n",
    "\n",
    "The cell below computes the proportion of reference mismatches among all MAVE variants that have been mapped."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "de7a5de0",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'There are 24878 instances of reference mismatch. This corresponds to a percentage of 0.996 (24878/2499044)'"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "mm_count = 0\n",
    "for key in diff_vars_dict:\n",
    "    mm_count = mm_count + len(diff_vars_dict[key])\n",
    "f\"There are {mm_count} instances of reference mismatch. This corresponds to a percentage of {round(100*mm_count/var_count,3)} ({mm_count}/{var_count})\""
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0ac2b919",
   "metadata": {},
   "source": [
    "## Compute the Number of Unique Pre-Mapped and Post-Mapped MAVE Variants\n",
    "\n",
    "The cell below computes the number of unique pre-mapped and post-mapped MAVE variants that have been processed using VRS."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "6a22205e",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'The number of unique pre-mapped MAVE variants is 363294. The number of unique post-mapped MAVE variants is 349972.'"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Determine total number of paired VRS alleles in data set\n",
    "allele_count = 0\n",
    "var_count = 0\n",
    "allele_ac_list_pre = []\n",
    "allele_ac_list_post = []\n",
    "for key in score_sets:\n",
    "    if key != \"urn:mavedb:00000072-a-1\" and key != \"urn:mavedb:00000105-a-1\":\n",
    "        f = open(f\"analysis_files/mappings/{key[11::]}.json\")\n",
    "        dat = json.load(f)\n",
    "        dat = dat[\"mapped_scores\"]\n",
    "\n",
    "        for j in range(len(dat)):\n",
    "            var_count += 1\n",
    "            if \"members\" not in dat[j][\"post_mapped\"]:\n",
    "                allele_count += 1\n",
    "                allele_ac_list_pre.append(dat[j][\"pre_mapped\"][\"id\"])\n",
    "                allele_ac_list_post.append(dat[j][\"post_mapped\"][\"id\"])\n",
    "            else:\n",
    "                for k in range(len(dat[j][\"post_mapped\"][\"members\"])):\n",
    "                    allele_count += 1\n",
    "                    allele_ac_list_pre.append(dat[j][\"pre_mapped\"][\"members\"][k][\"id\"])\n",
    "                    allele_ac_list_post.append(dat[j][\"post_mapped\"][\"members\"][k][\"id\"])\n",
    "\n",
    "\n",
    "f\"The number of unique pre-mapped MAVE variants is {len(set(allele_ac_list_pre))}. The number of unique post-mapped MAVE variants is {len(set(allele_ac_list_post))}.\""
   ]
  },
  {
   "cell_type": "markdown",
   "id": "92bc5e87",
   "metadata": {},
   "source": [
    "## Generate VRS Allele ID Dictionary\n",
    "Run the cell below to generate a VRS allele ID dictionary for the MAVE variants. The dictionary is keyed by the post-mapped VRS allele IDs, and each key's contents contain a list of corresponding pre-mapped VRS allele IDs."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "3eb301ed",
   "metadata": {},
   "outputs": [],
   "source": [
    "allele_list_dict = {}\n",
    "for key in score_sets:\n",
    "    if key != \"urn:mavedb:00000072-a-1\" and key != \"urn:mavedb:00000105-a-1\":\n",
    "        f = open(f\"analysis_files/mappings/{key[11::]}.json\")\n",
    "        dat = json.load(f)\n",
    "        dat = dat[\"mapped_scores\"]\n",
    "\n",
    "        for j in range(len(dat)):\n",
    "            if \"members\" not in dat[j][\"post_mapped\"]:\n",
    "                va = dat[j][\"post_mapped\"][\"id\"]\n",
    "                if va not in allele_list_dict:\n",
    "                    allele_list_dict[va] = [dat[j][\"pre_mapped\"][\"id\"]]\n",
    "                else:\n",
    "                    if dat[j][\"pre_mapped\"][\"id\"] in allele_list_dict[va]:\n",
    "                        continue\n",
    "                    else:\n",
    "                        tmp = allele_list_dict[va]\n",
    "                        tmp.append(dat[j][\"pre_mapped\"][\"id\"])\n",
    "                        allele_list_dict[va] = tmp\n",
    "            else:\n",
    "                for k in range(len(dat[j][\"post_mapped\"][\"members\"])):\n",
    "                    va = dat[j][\"post_mapped\"][\"members\"][k][\"id\"]\n",
    "                    if va not in allele_list_dict:\n",
    "                        allele_list_dict[va] = [dat[j][\"pre_mapped\"][\"members\"][k][\"id\"]]\n",
    "                    else:\n",
    "                        if dat[j][\"pre_mapped\"][\"members\"][k][\"id\"] in allele_list_dict[va]:\n",
    "                            continue\n",
    "                        tmp = allele_list_dict[va]\n",
    "                        tmp.append(dat[j][\"pre_mapped\"][\"members\"][k][\"id\"])\n",
    "                        allele_list_dict[va] = tmp"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "85b4e915",
   "metadata": {},
   "source": [
    "### Summary Statistics for VRS Allele ID Dictionary"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "fc87cbbe",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'There are 349972 pre-mapped MAVE variants in the dictionary. 9553 post-mapped MAVE variants have 2 or more corresponding pre-mapped MAVE variants.'"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "count = 0\n",
    "for key in allele_list_dict:\n",
    "    if len(allele_list_dict[key]) > 1:\n",
    "        count += 1\n",
    "f\"There are {len(allele_list_dict)} pre-mapped MAVE variants in the dictionary. {count} post-mapped MAVE variants have 2 or more corresponding pre-mapped MAVE variants.\""
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "venv",
   "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.11.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}