CNN model training (Apr_17_2021).ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "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": "7ce87e07",
   "metadata": {
    "id": "7ce87e07"
   },
   "source": [
    "### 1.   Prepare input data\n",
    "-----\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "625290cf",
   "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": 4,
   "id": "194edf7f",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Load input file\n",
    "input_df = pd.read_csv('TIL_transformed_data.csv')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "b483f677",
   "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": "61a6b8dc",
   "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": "1e017399",
   "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": 13,
   "id": "e7d3b06c",
   "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": 15,
   "id": "93c9cb02",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Function: transform data into 1d images\n",
    "def img_transform(input_data):\n",
    "    gene_num = input_data.shape[1]\n",
    "    cell_num = input_data.shape[0]\n",
    "    final_img = input_data.reshape((cell_num, gene_num, 1))\n",
    "\n",
    "    return final_img"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "5f407901",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Transform data into 1d images\n",
    "X_train_img = img_transform(X_train)\n",
    "X_val_img = img_transform(X_val)\n",
    "X_test_img = img_transform(X_test)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "id": "7b1ac6f1",
   "metadata": {
    "id": "7b1ac6f1"
   },
   "outputs": [],
   "source": [
    "# Define class of SingleCellDataset\n",
    "class SingleCellDataset(Dataset):\n",
    "    # Initialize\n",
    "    def __init__(self, input_img, labels, size):\n",
    "        self.label = torch.tensor(labels).squeeze()\n",
    "        self.img = torch.tensor(input_img).reshape((len(self.label), 1, size))\n",
    "    \n",
    "    # Total number of cells\n",
    "    def __len__(self):\n",
    "        return len(self.label)\n",
    "    \n",
    "    # Index images\n",
    "    def __getitem__(self, idx):\n",
    "        img_value = self.img[idx, :]\n",
    "        label_value = self.label[idx]\n",
    "        return img_value, label_value"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "id": "9be48609",
   "metadata": {
    "id": "9be48609"
   },
   "outputs": [],
   "source": [
    "# Create datasets\n",
    "img_size = X_train_img[0].shape[0]\n",
    "\n",
    "TrainDataSet = SingleCellDataset(X_train_img, y_train, img_size)\n",
    "ValDataSet = SingleCellDataset(X_val_img, y_val, img_size)\n",
    "TestDataSet = SingleCellDataset(X_test_img, y_test, img_size)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 48,
   "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": "3702198c",
   "metadata": {
    "id": "3702198c"
   },
   "source": [
    "### 2.   Model training\n",
    "-----"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "id": "de9bdb90",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define the model class\n",
    "class NeuralNetwork(nn.Module):\n",
    "    def __init__(self):\n",
    "        super(NeuralNetwork, self).__init__()\n",
    "        self.conv_stack = nn.Sequential(\n",
    "            nn.Conv1d(1,6,1, stride=1),\n",
    "            nn.BatchNorm1d(6),\n",
    "\n",
    "            nn.Conv1d(6,6,4, stride=4),\n",
    "            nn.BatchNorm1d(6),\n",
    "            nn.LeakyReLU(),\n",
    "            nn.MaxPool1d(2),\n",
    "\n",
    "            nn.Conv1d(6,6,4, stride=4),\n",
    "            nn.BatchNorm1d(6),\n",
    "            nn.LeakyReLU(),\n",
    "            nn.MaxPool1d(2),\n",
    "\n",
    "            nn.Conv1d(6,6,4, stride=4),\n",
    "            nn.BatchNorm1d(6),\n",
    "            nn.LeakyReLU(),\n",
    "            nn.MaxPool1d(2),\n",
    "\n",
    "            nn.Flatten(start_dim=1),\n",
    "            nn.Linear(10*1*6, 16),\n",
    "            nn.LeakyReLU(),\n",
    "            nn.BatchNorm1d(16),\n",
    "\n",
    "            nn.Linear(16, 1),\n",
    "            nn.Sigmoid()\n",
    "        )\n",
    "\n",
    "    def forward(self, x):\n",
    "        y = self.conv_stack(x)\n",
    "        return y"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c69ed129",
   "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": "816f6639",
   "metadata": {},
   "outputs": [],
   "source": [
    "device = \"cuda\" if torch.cuda.is_available() else \"cpu\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "id": "e772614d",
   "metadata": {
    "id": "e772614d"
   },
   "outputs": [],
   "source": [
    "# Training loop\n",
    "def train(dataloader, model, loss_fn, optimizer, scheduler, 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).type(torch.double)\n",
    "\n",
    "        # Compute prediction error\n",
    "        pred = model(X)\n",
    "        loss = loss_fn(pred, y.unsqueeze(1))\n",
    "        train_loss += loss.item()\n",
    "\n",
    "        correct_high += (torch.logical_and(torch.round(pred.squeeze()) == 1, y == 1)).type(torch.float).sum().item()\n",
    "        correct_low += (torch.logical_and(torch.round(pred.squeeze()) == 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\"Correctly classified as high in train set: {correct_high}/{high_count} ({100*correct_high/high_count:>4f}%)\")\n",
    "    print(f\"Correctly classified as low in train set: {correct_low}/{low_count} ({100*correct_low/low_count:>4f}%)\\n\")\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",
    "    scheduler.step()\n",
    "    \n",
    "    return train_loss, 100*correct_high/high_count, 100*correct_low/low_count"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0e9d7998",
   "metadata": {},
   "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.001)\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
}