nexclade-comparison.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "ba1397c2-f1c9-4ac1-9340-d179a3747d19",
   "metadata": {},
   "source": [
    "# Nextclade comparison \n",
    "1. Select a subset of sequences not used in train/val/test for the 8 clades the model was trained on \n",
    "    - [x] consider the data in datasets.json\n",
    "    - [x] drop those files from the GISAID metadata in  /home/disco/Github/GISAID/metadata_tsv_2021_11_11/metadata.tsv\n",
    "    - [x] then, select new files from the metadata for the 8 clades\n",
    "2. Select a new `global-test` dataset\n",
    "    - [x] Extract all the selected sequences in individual fasta files to evaluate the model\n",
    "    - [x] Also, write all the extracted sequences in a single fasta file to run nextclade\n",
    "3. Compare results\n",
    "    - map all predictions to the same nomenclature (GISAID). \n",
    "    - The output from Nextclade is a column called `clade` corresponding to one of the `Nextrain Clade` that can be found in [this](https://covariants.org/variants) table.\n",
    "    - To compare results, map the `Nextrain Clade` (nextclade output) to the `Pangolineage`, that can be found also in the GISAID metadata for each one of the selected sequences.\n",
    "    - Once all predictions are in the GISAID nomenclature, proceed with a comparison of results. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "70e27ce5-66c3-4bc5-8f28-a01f698a9947",
   "metadata": {},
   "outputs": [],
   "source": [
    "import json \n",
    "import pandas as pd\n",
    "from tqdm import tqdm\n",
    "tqdm.pandas()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "f6650735-03ec-4e52-8c4e-351fe79c550a",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "50718"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "with open(\"datasets.json\",\"r\") as fp: \n",
    "    datasets = json.load(fp)\n",
    "    \n",
    "npy_files = []\n",
    "for list_set in datasets.values():\n",
    "    npy_files.extend(list_set)\n",
    "    \n",
    "len(npy_files) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "c526ab6b-0383-4522-89ea-41cbbcdf02bf",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Clean dataset\n",
    "npy_files = [file for file in  npy_files if \"(1)\" not in file]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "817a7235-0920-442e-ad79-8ad00f44d638",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(44073,\n",
       " 'npy-8-mer/hCoV-19/V/hCoV-19_England_SHEF-CE129_2020|2020-03-22|2020-04-29.npy')"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "len(npy_files), npy_files[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "a7f4c482-b9d7-43b4-b703-1d20dd6e2480",
   "metadata": {},
   "outputs": [],
   "source": [
    "npy2fastaid = lambda npy: npy.split(\"/\")[-1].replace(\".npy\",\".fasta\").replace(\"_\",\"/\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "dae17e38-1179-4ab1-b225-4f1fdad5d130",
   "metadata": {},
   "outputs": [],
   "source": [
    "# list of fastaid used in train/val/test\n",
    "fasta_id_used = [npy2fastaid(file) for file in npy_files]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "95e3504d-ed6c-488c-97df-346e6a91a6a3",
   "metadata": {},
   "source": [
    "# load metadata"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "c6eb3a99-cb83-438f-9510-ca05c52af340",
   "metadata": {},
   "outputs": [],
   "source": [
    "CLADES = ['S','L','G','V','GR','GH','GV','GK']#,'GRY']\n",
    "COLS = [\"Virus name\", \"Accession ID\", \"Collection date\", \"Submission date\",\"Clade\", \"Host\", \"Is complete?\"]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "33cba5d1-0bf4-419d-bd04-9712be5496d4",
   "metadata": {},
   "outputs": [],
   "source": [
    "path_metadata = \"/home/disco/Github/GISAID/metadata_tsv_2021_11_11/metadata.tsv\"\n",
    "metadata = pd.read_csv(path_metadata, sep=\"\\t\", usecols=COLS)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "b0d83279-cdb5-4918-a776-aba30c4eb2c3",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(5045607, 7)"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "metadata.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "33b63d0d-6609-4a50-83dd-73216c536466",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Remove NaN in Clades and not-complete sequences\n",
    "metadata.dropna(axis=\"rows\",\n",
    "            how=\"any\",\n",
    "            subset=[\"Is complete?\", \"Clade\"], \n",
    "            inplace=True,\n",
    "            )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "ceff7f15-a861-460b-9778-887fee0ca1b4",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(4956947, 7)"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "metadata.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "b5b94a97-844f-4a51-aa51-d02b90b42a95",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Filter by Clades and Host\n",
    "CLADES = tuple(clade for clade in CLADES)\n",
    "metadata.query(f\"`Clade` in {CLADES} and `Host`=='Human'\", inplace=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "22ee7e47-297b-476b-b6c2-d8ef0ed7bc22",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(3986669, 7)"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "metadata.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "63bbaad2-8e7d-4839-a9ed-e6f251709e1f",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100%|████████████████████████████████| 3986669/3986669 [01:01<00:00, 65119.95it/s]\n"
     ]
    }
   ],
   "source": [
    "# Generate id of sequences in fasta file: \"Virus name|Accession ID|Collection date\"\n",
    "metadata[\"fasta_id\"] = metadata.progress_apply(lambda row: \"|\".join([row[\"Virus name\"],row[\"Collection date\"],row[\"Submission date\"]]), axis=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "id": "3e74cee9-3524-4b4e-b824-8e5e0c0f34a0",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100%|██████████████████████████████████████████████| 8/8 [49:51<00:00, 373.97s/it]\n"
     ]
    }
   ],
   "source": [
    "import random\n",
    "from collections import namedtuple\n",
    "# subsample \n",
    "SAMPLES_PER_CLADE = 5000\n",
    "SampleClade = namedtuple(\"SampleClade\", [\"fasta_id\",\"clade\"])\n",
    "list_fasta_selected = []\n",
    "for clade in tqdm(CLADES):\n",
    "    samples_clade = metadata.query(f\"`Clade` == '{clade}'\")[\"fasta_id\"].tolist()\n",
    "    samples_clade = [sample for sample in samples_clade if sample not in fasta_id_used]\n",
    "    random.shuffle(samples_clade)\n",
    "    # select 'SAMPLES_PER_CLADE' samples for each clade, or all of them if available samples are less than required\n",
    "    list_fasta_selected.extend([SampleClade(fasta_id, clade) for fasta_id in samples_clade[:SAMPLES_PER_CLADE]])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "id": "6643b98a-304b-4591-8325-37168211fc4b",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "False"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "metadata[\"fasta_id\"][0] in fasta_id_used"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "id": "a93d9314-85ca-4068-978b-b3dc78b4ac80",
   "metadata": {},
   "outputs": [],
   "source": [
    "from pathlib import Path\n",
    "Path(\"nextclade-comparison\").mkdir(exist_ok=True, parents=True)\n",
    "pd.DataFrame(list_fasta_selected).to_csv(\"nextclade-comparison/undersample_by_clade.csv\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "id": "b46c4c4e-a026-4cc7-821c-d37a79d07a4b",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "clade\n",
       "G     5000\n",
       "GH    5000\n",
       "GK    5000\n",
       "GR    5000\n",
       "GV    5000\n",
       "L     5000\n",
       "S     5000\n",
       "V     5000\n",
       "dtype: int64"
      ]
     },
     "execution_count": 26,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "pd.DataFrame(list_fasta_selected).groupby(\"clade\").size()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0b53906e-1d09-40a0-a3c0-50fc75b26ff3",
   "metadata": {},
   "source": [
    "___\n",
    "# Compare results\n",
    "- https://covariants.org/variants\n",
    "- https://www.ncbi.nlm.nih.gov/pmc/articles/PMC8184080/ supplementary data contains a table (for the indian sequences) with all the linages (Nextclade, GISAID, Pangolin)\n",
    "- https://cov-lineages.org/lineage_list.html only pangolin\n",
    "\n",
    "### Relation between lineages are saved in dictionaries at `/master-lineages`: \n",
    "- **nextclade2pango**: obtained from covariants.orig/variants\n",
    "- **pango2gisaid**: obtained from the metadata of GISAID"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "3e077667-eed5-4eeb-b1c7-52b39fc29001",
   "metadata": {},
   "outputs": [],
   "source": [
    "import json\n",
    "import pandas as pd\n",
    "pd.set_option(\"display.max_columns\", None)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "id": "7018991a-c712-43c4-84cc-b5911cf30b1a",
   "metadata": {},
   "outputs": [],
   "source": [
    "# load dictionaries from /master-lineages\n",
    "with open(\"master-lineages/nextclade2pango.json\",\"r\") as fp:\n",
    "    nextclade2pango = json.load(fp) # keys are list of pango clades\n",
    "\n",
    "with open(\"master-lineages/pango2gisaid.json\",\"r\") as fp:\n",
    "    pango2gisaid = json.load(fp) # keys are individual gisaid clades\n",
    "    \n",
    "def nextclade2gisaid(nextclade):\n",
    "    list_pango = nextclade2pango.get(nextclade, [\"not found\"])\n",
    "    if \"not found\" in list_pango:\n",
    "        return \"not found\"\n",
    "    # map pango lineages to gisaid clades\n",
    "    list_gisaid = []\n",
    "    for pango in list_pango: \n",
    "        list_gisaid.append(pango2gisaid.get(pango))\n",
    "        \n",
    "    return \",\".join(set(list_gisaid))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "id": "03e6925b-2626-46d3-9c2c-7580b4561d17",
   "metadata": {},
   "outputs": [],
   "source": [
    "# nextclade predictions\n",
    "path_nextclade_preds = \"/home/disco/Github/nextclade/output-comparison-nextclade/nextclade.tsv\"\n",
    "nextclade_preds = pd.read_csv(path_nextclade_preds, sep=\"\\t\", usecols=[\"seqName\",\"clade\"])\n",
    "nextclade_preds.rename({'clade': 'nextclade'}, axis=1, inplace=True) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "id": "76fdefe6-ecf1-4439-880c-9260ce00bb47",
   "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>seqName</th>\n",
       "      <th>nextclade</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>hCoV-19/Angola/CERI-KRISP-K019085/2021|2021-03...</td>\n",
       "      <td>21B (Kappa)</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>hCoV-19/Argentina/INEI018119/2020|2020-06-11|2...</td>\n",
       "      <td>20A</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>hCoV-19/Australia/ACT0179/2021|2021-08-19|2021...</td>\n",
       "      <td>21J (Delta)</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>hCoV-19/Australia/NSW-ICPMR-11015/2021|2021-10...</td>\n",
       "      <td>21J (Delta)</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>hCoV-19/Australia/NSW-ICPMR-11682/2021|2021-10...</td>\n",
       "      <td>21J (Delta)</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                                             seqName    nextclade\n",
       "0  hCoV-19/Angola/CERI-KRISP-K019085/2021|2021-03...  21B (Kappa)\n",
       "1  hCoV-19/Argentina/INEI018119/2020|2020-06-11|2...          20A\n",
       "2  hCoV-19/Australia/ACT0179/2021|2021-08-19|2021...  21J (Delta)\n",
       "3  hCoV-19/Australia/NSW-ICPMR-11015/2021|2021-10...  21J (Delta)\n",
       "4  hCoV-19/Australia/NSW-ICPMR-11682/2021|2021-10...  21J (Delta)"
      ]
     },
     "execution_count": 30,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "nextclade_preds.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 53,
   "id": "4056cfaf-d96b-42bd-8e49-d71690e7cc4c",
   "metadata": {},
   "outputs": [],
   "source": [
    "path_model_preds = \"results_nextclade_comparison.csv\"\n",
    "model_preds = pd.read_csv(path_model_preds)\n",
    "\n",
    "# get ground truth for the sequence\n",
    "model_preds[\"ground_truth\"] = model_preds[\"path_fasta\"].apply(lambda path: path.split(\"/\")[2])\n",
    "\n",
    "# extract seqName to merge with nextclade predictions\n",
    "model_preds[\"seqName\"] = model_preds[\"path_fasta\"].apply(lambda path: path.split(\"/\")[-1].replace(\"_\",\"/\")[:-6])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 54,
   "id": "fe764226-e2b9-4ac3-bb23-140e9de8e3fe",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "There are 36491 sequences shared between nextclade and the model\n",
      "Preds in nextclade 38749\n",
      "Preds in model 38749\n"
     ]
    }
   ],
   "source": [
    "seqs_in_common = len(set(model_preds.seqName).intersection(set(nextclade_preds.seqName)))\n",
    "print(f\"There are {seqs_in_common} sequences shared between nextclade and the model\")\n",
    "print(f\"Preds in nextclade {nextclade_preds.shape[0]}\")\n",
    "print(f\"Preds in model {model_preds.shape[0]}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 55,
   "id": "6a5bf31d-24b8-4e9d-857a-373274bbcffa",
   "metadata": {},
   "outputs": [],
   "source": [
    "#list(set(model_preds.seqName) - set(nextclade_preds.seqName))[:10]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 56,
   "id": "82f6204b-3774-4ea0-91f2-f14b45848ec3",
   "metadata": {},
   "outputs": [],
   "source": [
    "#list(set(nextclade_preds.seqName)-set(model_preds.seqName))[:10]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 57,
   "id": "3b9fd7a8-133d-406f-bef4-c13628e5ad8b",
   "metadata": {},
   "outputs": [],
   "source": [
    "# merge predictions\n",
    "preds = model_preds.merge(nextclade_preds, how=\"inner\", on=\"seqName\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 58,
   "id": "1907dcbb-ee1e-45bc-838c-5cf53d737415",
   "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>Unnamed: 0</th>\n",
       "      <th>path_fasta</th>\n",
       "      <th>prob</th>\n",
       "      <th>pred_class</th>\n",
       "      <th>ground_truth</th>\n",
       "      <th>seqName</th>\n",
       "      <th>nextclade</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>0</td>\n",
       "      <td>data-nextclade-comparison/hCoV-19/GK/hCoV-19_U...</td>\n",
       "      <td>0.899028</td>\n",
       "      <td>GK</td>\n",
       "      <td>GK</td>\n",
       "      <td>hCoV-19/USA/MD-CDC-LC0260562/2021|2021-09-08|2...</td>\n",
       "      <td>21J (Delta)</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>1</td>\n",
       "      <td>data-nextclade-comparison/hCoV-19/GK/hCoV-19_A...</td>\n",
       "      <td>0.874016</td>\n",
       "      <td>GK</td>\n",
       "      <td>GK</td>\n",
       "      <td>hCoV-19/Aruba/AW-RIVM-52828/2021|2021-08-10|20...</td>\n",
       "      <td>21J (Delta)</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>2</td>\n",
       "      <td>data-nextclade-comparison/hCoV-19/GK/hCoV-19_E...</td>\n",
       "      <td>0.774449</td>\n",
       "      <td>GK</td>\n",
       "      <td>GK</td>\n",
       "      <td>hCoV-19/England/MILK-1C6044F/2021|2021-08-22|2...</td>\n",
       "      <td>21J (Delta)</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>3</td>\n",
       "      <td>data-nextclade-comparison/hCoV-19/GK/hCoV-19_F...</td>\n",
       "      <td>0.827125</td>\n",
       "      <td>GK</td>\n",
       "      <td>GK</td>\n",
       "      <td>hCoV-19/France/NAQ-CERBAHC-0717768/2021|2021|2...</td>\n",
       "      <td>21J (Delta)</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>4</td>\n",
       "      <td>data-nextclade-comparison/hCoV-19/GK/hCoV-19_A...</td>\n",
       "      <td>0.879893</td>\n",
       "      <td>GK</td>\n",
       "      <td>GK</td>\n",
       "      <td>hCoV-19/Aruba/AW-RIVM-51660/2021|2021-08-08|20...</td>\n",
       "      <td>21J (Delta)</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   Unnamed: 0                                         path_fasta      prob  \\\n",
       "0           0  data-nextclade-comparison/hCoV-19/GK/hCoV-19_U...  0.899028   \n",
       "1           1  data-nextclade-comparison/hCoV-19/GK/hCoV-19_A...  0.874016   \n",
       "2           2  data-nextclade-comparison/hCoV-19/GK/hCoV-19_E...  0.774449   \n",
       "3           3  data-nextclade-comparison/hCoV-19/GK/hCoV-19_F...  0.827125   \n",
       "4           4  data-nextclade-comparison/hCoV-19/GK/hCoV-19_A...  0.879893   \n",
       "\n",
       "  pred_class ground_truth                                            seqName  \\\n",
       "0         GK           GK  hCoV-19/USA/MD-CDC-LC0260562/2021|2021-09-08|2...   \n",
       "1         GK           GK  hCoV-19/Aruba/AW-RIVM-52828/2021|2021-08-10|20...   \n",
       "2         GK           GK  hCoV-19/England/MILK-1C6044F/2021|2021-08-22|2...   \n",
       "3         GK           GK  hCoV-19/France/NAQ-CERBAHC-0717768/2021|2021|2...   \n",
       "4         GK           GK  hCoV-19/Aruba/AW-RIVM-51660/2021|2021-08-08|20...   \n",
       "\n",
       "     nextclade  \n",
       "0  21J (Delta)  \n",
       "1  21J (Delta)  \n",
       "2  21J (Delta)  \n",
       "3  21J (Delta)  \n",
       "4  21J (Delta)  "
      ]
     },
     "execution_count": 58,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "preds.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 59,
   "id": "70b1744a-bc8f-4c04-8ece-06a10fe70752",
   "metadata": {},
   "outputs": [],
   "source": [
    "# map nextclade prediction to gisaid\n",
    "preds[\"nextclade2gisaid\"] = preds[\"nextclade\"].apply(lambda pred: nextclade2gisaid(pred))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 60,
   "id": "820859b7-238f-4580-9a5c-0194afa395ab",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "nextclade2gisaid\n",
       "G,GR           630\n",
       "GK            7791\n",
       "GK,O          5604\n",
       "GRY            270\n",
       "GV             647\n",
       "O             1064\n",
       "not found    20485\n",
       "dtype: int64"
      ]
     },
     "execution_count": 60,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "preds.groupby([\"nextclade2gisaid\"]).size()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 62,
   "id": "845fe829-34c3-48af-9f4e-6118a0c733e4",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "ground_truth  nextclade2gisaid\n",
       "G             G,GR                   7\n",
       "              GK                   525\n",
       "              GK,O                2753\n",
       "              GRY                  103\n",
       "              GV                    90\n",
       "              O                     26\n",
       "              not found           1066\n",
       "GH            G,GR                 623\n",
       "              GK                   275\n",
       "              GK,O                 911\n",
       "              GV                   482\n",
       "              O                   1036\n",
       "              not found           1207\n",
       "GK            GK,O                   8\n",
       "              GRY                  167\n",
       "              not found           4488\n",
       "GR            GK                  2271\n",
       "              GK,O                1912\n",
       "              GV                    75\n",
       "              not found            219\n",
       "GV            GK                  4720\n",
       "              GK,O                  15\n",
       "              not found            139\n",
       "L             GK,O                   2\n",
       "              not found           4330\n",
       "S             GK,O                   1\n",
       "              O                      1\n",
       "              not found           4474\n",
       "V             GK,O                   2\n",
       "              O                      1\n",
       "              not found           4562\n",
       "dtype: int64"
      ]
     },
     "execution_count": 62,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "preds.groupby([\"ground_truth\",\"nextclade2gisaid\"]).size()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0eb14cd8-2b61-4395-91ea-8bf55bcd7645",
   "metadata": {},
   "source": [
    "___\n",
    "### Metrics model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 66,
   "id": "c03cdb0c-d800-4a72-8f58-34fd362a9249",
   "metadata": {},
   "outputs": [],
   "source": [
    "from collections import namedtuple\n",
    "from sklearn.metrics import precision_recall_fscore_support"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 65,
   "id": "efc3e444-9f3b-4195-b753-636f5c619c94",
   "metadata": {},
   "outputs": [],
   "source": [
    "CLADES = ['S','L','G','V','GR','GH','GV','GK']\n",
    "y_true = preds.ground_truth\n",
    "y_pred = preds.pred_class\n",
    "precision, recall, fscore, support = precision_recall_fscore_support(y_true, y_pred, average=None, labels=CLADES)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 67,
   "id": "e6124c48-d49f-48a4-aa46-16ddc1558f13",
   "metadata": {},
   "outputs": [],
   "source": [
    "list_metrics = []\n",
    "Metrics = namedtuple(\"Metrics\", [\"clade\",\"precision\", \"recall\"])#, \"fscore\", \"support\"]\n",
    "for j,clade in enumerate(CLADES): \n",
    "    list_metrics.append(\n",
    "        Metrics(clade, precision[j], recall[j])\n",
    "    )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 68,
   "id": "97c22bd7-922c-46ad-b321-3e12a5149c9a",
   "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>clade</th>\n",
       "      <th>precision</th>\n",
       "      <th>recall</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>S</td>\n",
       "      <td>0.999548</td>\n",
       "      <td>0.988606</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>L</td>\n",
       "      <td>0.996724</td>\n",
       "      <td>0.983149</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>G</td>\n",
       "      <td>0.858899</td>\n",
       "      <td>0.963020</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>V</td>\n",
       "      <td>0.992771</td>\n",
       "      <td>0.992771</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>GR</td>\n",
       "      <td>0.997117</td>\n",
       "      <td>0.926960</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td>GH</td>\n",
       "      <td>0.992258</td>\n",
       "      <td>0.989413</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>6</th>\n",
       "      <td>GV</td>\n",
       "      <td>0.997910</td>\n",
       "      <td>0.979688</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>7</th>\n",
       "      <td>GK</td>\n",
       "      <td>0.965473</td>\n",
       "      <td>0.959468</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "  clade  precision    recall\n",
       "0     S   0.999548  0.988606\n",
       "1     L   0.996724  0.983149\n",
       "2     G   0.858899  0.963020\n",
       "3     V   0.992771  0.992771\n",
       "4    GR   0.997117  0.926960\n",
       "5    GH   0.992258  0.989413\n",
       "6    GV   0.997910  0.979688\n",
       "7    GK   0.965473  0.959468"
      ]
     },
     "execution_count": 68,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "pd.DataFrame(list_metrics)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d1afe5d9-05d4-4daa-9315-fed3487cd785",
   "metadata": {},
   "source": [
    "___\n",
    "metrics"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d7dc4d36-0083-4195-8bf4-8d2b928a76e3",
   "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.9.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}