examples/tutorials/biotrainer_intro_compare.ipynb

{
 "cells": [
  {
   "metadata": {},
   "cell_type": "markdown",
   "source": "# Biotrainer Introduction Notebook: Compare results",
   "id": "79db2ec60089f7b9"
  },
  {
   "metadata": {},
   "cell_type": "markdown",
   "source": "In the following tutorial, we will use [biotrainer](https://github.com/sacdallago/biotrainer) to train a baseline model (negative control) for an existing predictor. We will use one_hot_encodings for that. Can we claim that our predictor performs better than the baseline?",
   "id": "5a0d1d36d4c385dc"
  },
  {
   "metadata": {},
   "cell_type": "markdown",
   "source": "## 0. Install Biotrainer",
   "id": "5b00d2a3fe8688b8"
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": [
    "# Create a virtual environment before, then run:\n",
    "!pip install git+https://github.com/sacdallago/biotrainer.git@develop"
   ],
   "id": "5d0613551163a50c",
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "markdown",
   "source": [
    "## 1. The prediction task and predictor\n",
    "\n",
    "In our example, we already trained a regression model on the [FLIP meltome dataset split mixed](https://github.com/J-SNACKKB/FLIP/tree/main/splits/meltome) with [ProtT5](https://github.com/agemagician/ProtTrans) embeddings. So, we created a model that predicts the meltdown temperature for a given protein sequence.\n",
    "\n",
    "A common metric to compare performance of regression models is the `root_mean_squared_error` (rmse). The \"CI\" variable means \"confidence interval\". Both the mean for the rmse and the ci come from bootstrapping the test set with the created model. This is to draw repeatedly from the test set with replacement (e.g. 30 times) and then calculate the performance on that set. You can learn more about bootstrapping here: https://math.mit.edu/~dav/05.dir/class24-prep-a.pdf \n",
    "\n",
    "Now, we already have the result on the test set for ProtT5 given:"
   ],
   "id": "2d42eb88b85165e3"
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": [
    "from biotrainer.protocols import Protocol\n",
    "\n",
    "task_protocol = Protocol.sequence_to_value\n",
    "prott5_rmse_mean = 9.564  # Not real data\n",
    "prott5_rmse_ci = 0.63  # Not real data\n",
    "\n",
    "# Bootstrapping config (sample_size == length of test set)\n",
    "bootstrapping_confidence_level = 0.05\n",
    "bootstrapping_iterations = 30"
   ],
   "id": "909b692e5ed08ffc",
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "markdown",
   "source": [
    "## 2. Downloading the data\n",
    "\n",
    "For this regression task, we can download one simple FASTA file that contains both our sequences and the respective targets (meltdown temperature of the sequence)."
   ],
   "id": "58a7bc3a1ab79d63"
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": [
    "import requests\n",
    "\n",
    "url = \"http://data.bioembeddings.com/public/FLIP/fasta/meltome/mixed_split.fasta\"\n",
    "headers = {\n",
    "    'User-Agent': 'Mozilla/5.0 (Windows NT 10.0; Win64; x64) AppleWebKit/537.36 (KHTML, like Gecko) Chrome/91.0.4472.124 Safari/537.36'\n",
    "}\n",
    "\n",
    "response = requests.get(url, headers=headers)\n",
    "with open(\"meltome_mixed_split.fasta\", \"wb\") as f:\n",
    "    f.write(response.content)"
   ],
   "id": "5b1239a28e13d9cf",
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "markdown",
   "source": [
    "## 3. Defining the biotrainer configuration\n",
    "\n",
    "We already know our sequences and targets are stored in the downloaded fasta file. We also know the protocol and which embedder to use. The rest of the configurations are straight-forward."
   ],
   "id": "27892a8185ce4ec9"
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": [
    "biotrainer_config = {\n",
    "    \"protocol\": task_protocol.name,  # We have a per-sequence regression task, so sequence_to_value\n",
    "    \"sequence_file\": \"meltome_mixed_split.fasta\",  # The path to the downloaded fasta file\n",
    "    \"model_choice\": \"FNN\",  # Fully connected neural netowrk\n",
    "    \"embedder_name\": \"one_hot_encoding\",  # Name of the embedding model\n",
    "    \"num_epochs\": 200,  # Define how long we want to train at max\n",
    "}"
   ],
   "id": "bef44a4ad63d1a67",
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "markdown",
   "source": [
    "## 4. Running the experiment\n",
    "\n",
    "We are now ready to train our baseline prediction model!"
   ],
   "id": "20a2a2f2b9360090"
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": [
    "from biotrainer.utilities.cli import headless_main\n",
    "\n",
    "results = headless_main(biotrainer_config)"
   ],
   "id": "5fb3861c4d855d00",
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "markdown",
   "source": [
    "## 5. Comparing the results\n",
    "\n",
    "After our baseline model finished training, we can now compare it to the ProtT5 result.\n",
    "\n",
    "By the way, did you notice that the model did not train for the given 200 epochs? That is because there is a built-in mechanism called `early stopping` which should be used for any deep learning model. It avoids overfitting by stopping the training after a certain amount of epochs that did not increase the performance on the validation set. This \"certain amount\" can be defined in the training configuration by using `\"patience\": 50 (or some other positive integer value)`.  "
   ],
   "id": "3701d2a0bc2fd607"
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": [
    "# First we get the bootstrapping results\n",
    "bootstrapping_ohe_dict = results[\"test_iterations_results\"][\"bootstrapping\"]\n",
    "ohe_rmse_mean = bootstrapping_ohe_dict[\"rmse\"][\"mean\"]\n",
    "ohe_rmse_ci = bootstrapping_ohe_dict[\"rmse\"][\"error\"]\n",
    "\n",
    "# Now let's compare them!\n",
    "print(f\"ProtT5: RMSE Mean: {prott5_rmse_mean} CI: {prott5_rmse_ci}\")\n",
    "print(f\"OHE: RMSE Mean: {ohe_rmse_mean} CI: {ohe_rmse_ci}\")"
   ],
   "id": "688ac9c9898db58a",
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "markdown",
   "source": [
    "## Final step: z-test\n",
    "\n",
    "Okay, only from looking at this, we can already calculate in our heads, that the confidence intervals do not overlap. That means that there is a significant difference between the two test set performances in our setup.\n",
    "\n",
    "But let's first verify that by using a statistical test, the so called \"Z-Test\". \n",
    "The Z-test essentially asks: 'Given these RMSE scores and their uncertainties from bootstrapping, what's the probability that the difference between ProtT5 and the simpler approach occurred by chance?'\n",
    "\n",
    "It works by:\n",
    "\n",
    "* Measuring how far apart the two RMSE means are\n",
    "* Taking into account how uncertain each measurement is (using the confidence intervals from our bootstrapping)\n",
    "* Comparing this to what we'd expect if there was no real difference between the approaches\n",
    "\n",
    "For example, if ProtT5 gives an RMSE of 3.2±0.2 and the simpler approach gives 4.1±0.2, the Z-test helps us confirm that ProtT5's better performance is 'real' and not just random variation. The test gives us a p-value: if it's less than 0.05, we can say that one approach is significantly better than the other."
   ],
   "id": "e6036c4a6409635c"
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": [
    "import numpy as np\n",
    "from scipy import stats\n",
    "\n",
    "def z_test(mean1, ci1, mean2, ci2):\n",
    "    # Convert CI to standard error for 95% CI\n",
    "    se1 = ci1 / 1.96\n",
    "    se2 = ci2 / 1.96\n",
    "    z_stat = (mean1 - mean2) / np.sqrt(se1**2 + se2**2)\n",
    "    p_value = 2 * (1 - stats.norm.cdf(abs(z_stat)))  # two-tailed\n",
    "    return z_stat, p_value\n",
    "\n",
    "z_stat, p_value = z_test(prott5_rmse_mean, prott5_rmse_ci, ohe_rmse_mean, ohe_rmse_ci)\n",
    "significant = p_value < 0.05\n",
    "\n",
    "print(z_stat, p_value)\n",
    "print(f\"As {p_value} {'<' if significant else '>='} 0.05, we can assume that there is {'a' if significant else 'no'} significant difference between the performance of ProtT5 and one hot encodings!\")"
   ],
   "id": "fbfbf5fe0f31c88d",
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "markdown",
   "source": "",
   "id": "dacc7cf2ead2a883"
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}