evaluation_on_outer_data.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "54424c79-98c0-48cd-bbf1-6ff9b8079b49",
   "metadata": {},
   "source": [
    "# Select sequences from Clades not used in training"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "2dc2f026-b12a-4fc0-ae79-7895bdb82030",
   "metadata": {},
   "outputs": [],
   "source": [
    "from tqdm.notebook import tqdm\n",
    "import pandas as pd\n",
    "tqdm.pandas()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "95b16bec-c53d-40d0-8d06-86c1ddcb0edd",
   "metadata": {},
   "outputs": [],
   "source": [
    "CLADES = ['GRY','O']\n",
    "COLS = [\"Virus name\", \"Accession ID\", \"Collection date\", \"Submission date\",\"Clade\", \"Host\", \"Is complete?\"]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "cb45abc5-7bab-4175-8d2f-43cec7925ef0",
   "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": 4,
   "id": "8becff46-4653-4f75-9ac2-3f3f036c5606",
   "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": 5,
   "id": "2882f79e-fcf7-477f-81b3-77b34d2d6670",
   "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": 8,
   "id": "a571f236-adf0-4868-9f25-4b6114ab6aa7",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "b536f6909ead4e719c129d20cfe6b458",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "  0%|          | 0/966768 [00:00<?, ?it/s]"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "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": 9,
   "id": "d6f511e7-bba0-4392-b33b-78215eeddc76",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "7c6e363dfb57415fb5976380de1d597b",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "  0%|          | 0/2 [00:00<?, ?it/s]"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import random\n",
    "from collections import namedtuple\n",
    "# subsample \n",
    "SAMPLES_PER_CLADE = 1000\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",
    "    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": 10,
   "id": "da09b976-a20e-4ef1-90c5-5906d2baf5d4",
   "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/outer_clades.csv\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "d8fee4ec-9cc8-4b7c-9e2a-f8715c297fc2",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "clade\n",
       "GRY    1000\n",
       "O      1000\n",
       "dtype: int64"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "pd.DataFrame(list_fasta_selected).groupby(\"clade\").size()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "41eb7ea1-674b-45d1-9933-1954d9e60766",
   "metadata": {},
   "source": [
    "___\n",
    "## Inference and histogram of probabilities"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "id": "5a7da4ac-06b0-4bcf-9178-89464f608efd",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np \n",
    "\n",
    "from pathlib import Path\n",
    "from Bio import SeqIO\n",
    "from supervised_dna.fcgr import FCGR\n",
    "from supervised_dna.utils import (\n",
    "    # array2img,\n",
    "    # find_matches,\n",
    "    # fcgrpos2kmers,\n",
    "    preprocess_seq\n",
    ")\n",
    "fcgr = FCGR(8)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "id": "1a49abcd-1290-4808-9f16-de551d44a042",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1952"
      ]
     },
     "execution_count": 24,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "list_fasta = list((Path(\"data-outer-comparison\").rglob(\"*.fasta\")))\n",
    "len(list_fasta) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "id": "b4a834a5-c4c9-49c0-819f-9a763d52178d",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      " **load model weights_path** : checkpoints/model-02-0.969.hdf5\n",
      "\n",
      "**Model created**\n"
     ]
    }
   ],
   "source": [
    "from supervised_dna import ModelLoader\n",
    "# -1- Load model with best weights\n",
    "loader = ModelLoader()\n",
    "MODEL  = \"resnet50_8mers\"\n",
    "CLADES = ['S','L','G','V','GR','GH','GV','GK']\n",
    "WEIGHTS_PATH = \"checkpoints/model-02-0.969.hdf5\"\n",
    "model  = loader(\n",
    "            model_name=MODEL,\n",
    "            n_outputs=len(CLADES),\n",
    "            weights_path=WEIGHTS_PATH,\n",
    "            ) # get compiled model from ./supervised_dna/models"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "id": "a87b34c3-459a-4a40-a904-500b7e6137b8",
   "metadata": {},
   "outputs": [],
   "source": [
    "# from collections import namedtuple\n",
    "# Results = namedtuple(\"Results\",[\"path_fasta\",\"prob\",\"pred_class\"])\n",
    "# list_results = []\n",
    "# for path_fasta in tqdm(list_test, desc=\"Computing Saliency Maps\"):\n",
    "#     prob, pred_class = compute_analysis(path_fasta)\n",
    "#     list_results.append(Results(path_fasta, prob, pred_class))\n",
    "\n",
    "# pd.DataFrame(list_results).to_csv(\"results_nextclade_comparison.csv\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 48,
   "id": "271116be-0250-4ecb-ae43-3a9981e41ce2",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "e61cc33f419e4e1990c1b4d0097dab46",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "  0%|          | 0/1952 [00:00<?, ?it/s]"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "from tqdm.notebook import tqdm\n",
    "from collections import namedtuple\n",
    "Results = namedtuple(\"Results\",[\"path_fasta\",\"prob\",\"pred_class\"])\n",
    "\n",
    "list_results = []\n",
    "for path_fasta in tqdm(list_fasta):\n",
    "    fasta = next(SeqIO.parse(path_fasta, \"fasta\"))\n",
    "    array_freq = fcgr(sequence=preprocess_seq(str(fasta.seq)))\n",
    "    # preproceesing (divide by 10) and add channel axis\n",
    "    input_model = np.expand_dims(array_freq/10. , axis=-1)\n",
    "    input_model = np.expand_dims(input_model,axis=0)\n",
    "    \n",
    "    # make prediction and get predicted clade\n",
    "    probs      = model.predict(input_model)[0]\n",
    "    pred_class = CLADES[np.argmax(probs)]\n",
    "    prob_class = probs.max()\n",
    "    \n",
    "    # save result\n",
    "    list_results.append(Results(path_fasta, prob_class, pred_class))\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "id": "62aec332-b83c-4d93-b6c1-57234d150e7d",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(0.9999782, 'GR')"
      ]
     },
     "execution_count": 46,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "prob.max(), CLADES[np.argmax(prob)]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 49,
   "id": "bf9df450-af6d-44d7-9fa8-c45aa68a916d",
   "metadata": {},
   "outputs": [],
   "source": [
    "pd.DataFrame(list_results).to_csv(\"outer-comparison/results_pred.csv\")"
   ]
  }
 ],
 "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
}