Model training - Apr 10 2022 - RQ.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0e462e84",
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "executionInfo": {
     "elapsed": 364,
     "status": "ok",
     "timestamp": 1630409305031,
     "user": {
      "displayName": "Daniel Bojar",
      "photoUrl": "https://lh3.googleusercontent.com/a-/AOh14GgqK9Tu9YrkihvM5n0N7oStKrVaKvnc25sL21EXvg=s64",
      "userId": "10339697633531698497"
     },
     "user_tz": -120
    },
    "id": "0e462e84",
    "outputId": "9ae23fd6-0474-4e16-f0df-102f4012fca8",
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "import torch\n",
    "from torch import nn\n",
    "from torch.utils.data import Dataset\n",
    "from torch.utils.data import DataLoader\n",
    "\n",
    "from sklearn.model_selection import train_test_split"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "da9c46b7",
   "metadata": {},
   "source": [
    "### 1.   Prepare input data\n",
    "-----"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f7ac683f",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Function: determine PHA-L read cut-offs for binary classification \n",
    "def categorize_lectin(data_all, quantile_high, quantile_low, ref_col_loc):\n",
    "    cutoff = np.quantile(data_all.iloc[:,ref_col_loc], [quantile_high, quantile_low], interpolation=\"nearest\").tolist()\n",
    "    print(f\"Cut-off for PHA-L high: {cutoff[0]}; Cut-off for PHA-L low: {cutoff[1]}\")\n",
    "    \n",
    "    high_indices = np.array(data_all.iloc[:,ref_col_loc]>=cutoff[0])\n",
    "    low_indices = np.array(data_all.iloc[:,ref_col_loc]<cutoff[1])\n",
    "    high_low_indices = np.logical_or(high_indices, low_indices)\n",
    "\n",
    "    high_count = high_indices.sum()\n",
    "    low_count = low_indices.sum()\n",
    "    \n",
    "    return cutoff, [high_indices, low_indices, high_low_indices], [high_count, low_count]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a7128392",
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "executionInfo": {
     "elapsed": 384,
     "status": "ok",
     "timestamp": 1630409305399,
     "user": {
      "displayName": "Daniel Bojar",
      "photoUrl": "https://lh3.googleusercontent.com/a-/AOh14GgqK9Tu9YrkihvM5n0N7oStKrVaKvnc25sL21EXvg=s64",
      "userId": "10339697633531698497"
     },
     "user_tz": -120
    },
    "id": "a7128392",
    "outputId": "528d881a-d533-46fd-c6cc-b195b38a5ee2"
   },
   "outputs": [],
   "source": [
    "# Load input file\n",
    "input_df = pd.read_csv('TIL_transformed_data.csv')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6be21275",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Process data: binary classification\n",
    "quantile_high, quantile_low = 0.75, 0.25\n",
    "cutoff, indices, count = categorize_lectin(input_df, quantile_high, quantile_low, -1)\n",
    "\n",
    "input_df.loc[indices[0], \"PHA-L\"] = 1\n",
    "input_df.loc[indices[1], \"PHA-L\"] = 0\n",
    "\n",
    "input_df = input_df.loc[indices[2], :]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f2b2f94e",
   "metadata": {},
   "outputs": [],
   "source": [
    "#y: class array\n",
    "y = input_df['PHA-L'].values \n",
    "#X: transcript data array\n",
    "X = input_df.iloc[:, 1:-1].values"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f8406fe2",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Split training, validation and test set\n",
    "X_train_val, X_test, y_train_val, y_test = train_test_split(\n",
    "    X, y, test_size=0.1, random_state=342, stratify=y)\n",
    "\n",
    "X_train, X_val, y_train, y_val = train_test_split(\n",
    "    X_train_val, y_train_val, test_size=0.2, random_state=2, stratify=y_train_val)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e6be0050",
   "metadata": {},
   "outputs": [],
   "source": [
    "# PHA-L high and PHA-L low counts in each set\n",
    "count_train = [y_train.sum(), len(y_train)-y_train.sum()]\n",
    "count_val = [y_val.sum(), len(y_val)-y_val.sum()]\n",
    "count_test = [y_test.sum(), len(y_test)-y_test.sum()]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7b1ac6f1",
   "metadata": {
    "id": "7b1ac6f1"
   },
   "outputs": [],
   "source": [
    "# Define class of SingleCellDataset\n",
    "class SingleCellDataset(Dataset):\n",
    "    # Initialize\n",
    "    def __init__(self, rna, lectin):\n",
    "        self.transcript = torch.tensor(rna, dtype=torch.float)\n",
    "        self.lectin = torch.tensor(lectin, dtype=torch.float)\n",
    "    \n",
    "    # Total number of cells\n",
    "    def __len__(self):\n",
    "        return self.transcript.shape[0]\n",
    "    \n",
    "    # Index cells\n",
    "    def __getitem__(self, idx):\n",
    "        transcript_value = self.transcript[idx, :]\n",
    "        lectin_value = self.lectin[idx]\n",
    "        return transcript_value, lectin_value"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9be48609",
   "metadata": {
    "id": "9be48609"
   },
   "outputs": [],
   "source": [
    "# Create datasets\n",
    "TrainDataSet = SingleCellDataset(X_train, y_train)\n",
    "ValDataSet = SingleCellDataset(X_val, y_val)\n",
    "TestDataSet = SingleCellDataset(X_test, y_test)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "76ddc59e",
   "metadata": {
    "id": "76ddc59e"
   },
   "outputs": [],
   "source": [
    "# Create dataloaders\n",
    "batch_size = 128\n",
    "\n",
    "train_data_loader = DataLoader(TrainDataSet, batch_size=batch_size, shuffle=True)\n",
    "val_data_loader = DataLoader(ValDataSet, batch_size=batch_size, shuffle=True)\n",
    "test_data_loader = DataLoader(TestDataSet, batch_size=X_test.shape[0], shuffle=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1260f531",
   "metadata": {},
   "source": [
    "### 2.   Model training\n",
    "-----"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4fc87be3",
   "metadata": {
    "id": "4fc87be3"
   },
   "outputs": [],
   "source": [
    "# Define the model class\n",
    "class NeuralNetwork(nn.Module):\n",
    "    def __init__(self, input_size):\n",
    "        super(NeuralNetwork, self).__init__()\n",
    "        self.linear_relu_stack = nn.Sequential(\n",
    "            nn.Linear(input_size, 128),\n",
    "            nn.LeakyReLU(),\n",
    "            nn.Dropout(0.4),\n",
    "            nn.BatchNorm1d(128),\n",
    "            nn.Linear(128, 64),\n",
    "            nn.LeakyReLU(),\n",
    "            nn.Dropout(0.4),\n",
    "            nn.BatchNorm1d(64),\n",
    "            nn.Linear(64, 16),\n",
    "            nn.LeakyReLU(),\n",
    "            nn.Dropout(0.2),\n",
    "            nn.BatchNorm1d(16),\n",
    "            nn.Linear(16, 8),\n",
    "            nn.LeakyReLU(),\n",
    "            nn.BatchNorm1d(8),\n",
    "            nn.Linear(8, 1),\n",
    "            nn.Sigmoid()\n",
    "        )\n",
    "\n",
    "    def forward(self, x):\n",
    "        y = self.linear_relu_stack(x)\n",
    "        return y"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b876fa74",
   "metadata": {},
   "outputs": [],
   "source": [
    "# function: call label based on preset probability cutoff\n",
    "def call_label(pred, prob_cutoff):\n",
    "    pred_label = []\n",
    "    for i in pred.squeeze():\n",
    "        if i >= prob_cutoff:\n",
    "            pred_label.append(1)\n",
    "        else:\n",
    "            pred_label.append(0)\n",
    "    return torch.tensor(pred_label).reshape(len(pred_label),)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9254b550",
   "metadata": {},
   "outputs": [],
   "source": [
    "device = \"cuda\" if torch.cuda.is_available() else \"cpu\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e772614d",
   "metadata": {
    "id": "e772614d"
   },
   "outputs": [],
   "source": [
    "# Training loop\n",
    "def train(dataloader, model, loss_fn, optimizer, scheduler, prob_cutoff, count):\n",
    "    size = len(dataloader.dataset)\n",
    "    num_batches = len(dataloader)\n",
    "    model.train()\n",
    "    train_loss, correct_high, correct_low = 0, 0, 0\n",
    "    high_count, low_count = count[0], count[1]\n",
    "    \n",
    "    for batch, (X, y) in enumerate(dataloader):\n",
    "        X, y = X.to(device), y.to(device)\n",
    "        y=torch.squeeze(y)\n",
    "\n",
    "        # Compute prediction error\n",
    "        pred = model(X)\n",
    "        pred_label = call_label(pred, prob_cutoff)\n",
    "        loss = loss_fn(pred, y.unsqueeze(1))\n",
    "        train_loss += loss.item()\n",
    "\n",
    "        correct_high += (torch.logical_and(pred_label == 1, y == 1)).type(torch.float).sum().item()\n",
    "        correct_low += (torch.logical_and(pred_label == 0, y == 0)).type(torch.float).sum().item()\n",
    "\n",
    "        # Backpropagation\n",
    "        optimizer.zero_grad()\n",
    "        loss.backward()\n",
    "        optimizer.step()\n",
    "\n",
    "        if batch % 10 == 0:\n",
    "            loss, current = loss.item(), batch * len(X)\n",
    "            print(f\"loss: {loss:>7f}  [{current:>5d}/{size:>5d}]\")\n",
    "    \n",
    "    train_loss /= num_batches\n",
    "    print(f\"Avg loss of training set: {train_loss:.4f} \\n\")\n",
    "    print(f\"Accuracy for 'PHA-L high' class of training set: {correct_high}/{high_count} ({100*correct_high/high_count:.4f}%)\")\n",
    "    print(f\"Accuracy for 'PHA-L low' class of training set: {correct_low}/{low_count} ({100*correct_low/low_count:.4f}%)\")\n",
    "    print(f\"Overall accuracy: {correct_high+correct_low}/{high_count+low_count} ({100*(correct_high+correct_low)/(high_count+low_count):.4f}%)\\n\")\n",
    "    \n",
    "    scheduler.step()\n",
    "    \n",
    "    return train_loss, 100*correct_high/high_count, 100*correct_low/low_count"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9b60d917",
   "metadata": {
    "id": "9b60d917"
   },
   "outputs": [],
   "source": [
    "# Validation loop\n",
    "def val(dataloader, model, loss_fn, prob_cutoff, count):\n",
    "    size = len(dataloader.dataset)\n",
    "    num_batches = len(dataloader)\n",
    "    model.eval()\n",
    "    val_loss, correct_high, correct_low = 0, 0, 0\n",
    "    high_count, low_count = count[0], count[1]\n",
    "    \n",
    "    with torch.no_grad():\n",
    "        for X, y in dataloader:\n",
    "            X, y = X.to(device), y.to(device)\n",
    "            y=torch.squeeze(y)\n",
    "            \n",
    "            pred = model(X)\n",
    "            pred_label = call_label(pred, prob_cutoff)\n",
    "            val_loss += loss_fn(pred, y.unsqueeze(1)).item()\n",
    "            correct_high += (torch.logical_and(pred_label == 1, y == 1)).type(torch.float).sum().item()\n",
    "            correct_low += (torch.logical_and(pred_label == 0, y == 0)).type(torch.float).sum().item()\n",
    "    \n",
    "    val_loss /= num_batches\n",
    "    print(f\"Avg loss of test set: {val_loss:.4f} \\n\")\n",
    "    print(f\"Accuracy for 'PHA-L high' class of validation set: {correct_high}/{high_count} ({100*correct_high/high_count:.4f}%)\")\n",
    "    print(f\"Accuracy for 'PHA-L low' class of validation set: {correct_low}/{low_count} ({100*correct_low/low_count:.4f}%)\")\n",
    "    print(f\"Overall accuracy: {correct_high+correct_low}/{high_count+low_count} ({100*(correct_high+correct_low)/(high_count+low_count):.4f}%)\\n\")\n",
    "    \n",
    "    return val_loss, 100*correct_high/high_count, 100*correct_low/low_count"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fd0d3cb0",
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "executionInfo": {
     "elapsed": 4986,
     "status": "ok",
     "timestamp": 1630409968509,
     "user": {
      "displayName": "Daniel Bojar",
      "photoUrl": "https://lh3.googleusercontent.com/a-/AOh14GgqK9Tu9YrkihvM5n0N7oStKrVaKvnc25sL21EXvg=s64",
      "userId": "10339697633531698497"
     },
     "user_tz": -120
    },
    "id": "fd0d3cb0",
    "outputId": "e18e2395-1bb9-4247-d1a9-8abf941eb0fb"
   },
   "outputs": [],
   "source": [
    "# Start training\n",
    "gene_number = X.shape[1]\n",
    "\n",
    "model = NeuralNetwork(input_size=gene_number).to(device)\n",
    "print(model)\n",
    "\n",
    "loss_fn = nn.BCELoss()\n",
    "prob_cutoff = 0.5\n",
    "optimizer = torch.optim.Adam(model.parameters(), lr=0.0001)\n",
    "epochs = 30\n",
    "scheduler = torch.optim.lr_scheduler.CosineAnnealingLR(optimizer, T_max=epochs)\n",
    "\n",
    "train_loss = []\n",
    "train_accuracy_high = []\n",
    "train_accuracy_low = []\n",
    "\n",
    "val_loss = []\n",
    "val_accuracy_high = []\n",
    "val_accuracy_low = []\n",
    "\n",
    "for t in range(epochs):\n",
    "    print(f\"Epoch {t+1}\\n-------------------------------\")\n",
    "    loss_train_epoch, acc_train_epoch_high, acc_train_epoch_low = train(train_data_loader, model, loss_fn, optimizer, scheduler, prob_cutoff, count_train)\n",
    "    train_loss.append(loss_train_epoch)\n",
    "    train_accuracy_high.append(acc_train_epoch_high)\n",
    "    train_accuracy_low.append(acc_train_epoch_low)\n",
    "    \n",
    "    loss_val_epoch, acc_val_epoch_high, acc_val_epoch_low = val(val_data_loader, model, loss_fn, prob_cutoff, count_val)\n",
    "    val_loss.append(loss_val_epoch)\n",
    "    val_accuracy_high.append(acc_val_epoch_high)\n",
    "    val_accuracy_low.append(acc_val_epoch_low)\n",
    "    \n",
    "print(\"Training finished.\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1c907b9c",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot training set and validation set loss\n",
    "plt.plot(np.arange(1, epochs+1), train_loss, label=\"train_loss\")\n",
    "plt.plot(np.arange(1, epochs+1), val_loss, label=\"test_loss\")\n",
    "plt.xlabel('Epoch')\n",
    "plt.ylabel(str(loss_fn))\n",
    "plt.legend(loc=1)\n",
    "plt.xticks(np.arange(0, epochs+2, step=2))\n",
    "plt.xlim(1, epochs)\n",
    "plt.grid()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a85c39b4",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot training set and validation set accuracy\n",
    "plt.plot(np.arange(1, epochs+1), train_accuracy_high, label=\"train_accuracy_high\")\n",
    "plt.plot(np.arange(1, epochs+1), train_accuracy_low, label=\"train_accuracy_low\")\n",
    "plt.plot(np.arange(1, epochs+1), val_accuracy_high, label=\"val_accuracy_high\")\n",
    "plt.plot(np.arange(1, epochs+1), val_accuracy_low, label=\"val_accuracy_low\")\n",
    "plt.xlabel('Epoch')\n",
    "plt.ylabel('Accuracy(%)')\n",
    "plt.legend(loc=4)\n",
    "plt.xticks(np.arange(0, epochs+2, step=2))\n",
    "plt.yticks(np.arange(0, 110, step=10))\n",
    "plt.xlim(1, epochs)\n",
    "plt.ylim(0,105)\n",
    "plt.grid()\n",
    "plt.show()"
   ]
  }
 ],
 "metadata": {
  "colab": {
   "collapsed_sections": [],
   "name": "RNA_lectin_ML_implementation_RQ.ipynb",
   "provenance": []
  },
  "interpreter": {
   "hash": "0c27f3dfbe5c91552ea375c193e935c23c9fbef877fb71378394b3d18f317895"
  },
  "kernelspec": {
   "display_name": "Python 3.8.11 64-bit ('sc_rna_lectin': conda)",
   "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.8.11"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}