Day3/TissueClassification_with_solution.ipynb

{
  "cells": [
    {
      "attachments": {},
      "cell_type": "markdown",
      "metadata": {
        "id": "pyoD2k8uLwpD"
      },
      "source": [
        "# Classifying tissues from gene expression profiles\n",
        "\n",
        "A problem clinicians often encounter is the primary origin of cancer. Often, metastatic tumors are never fully diagnosed even after various clinical tests to determine the primary site and are thus labeled cancers of unknown origin (CUP).\n",
        "\n",
        "RNA-seq data seems like a promising solution to the problem, since tissues differ substantially in their gene expression profiles.\n",
        "\n",
        "## Exercise:\n",
        "\n",
        "In this exercise, we will learn how to classify a gene expression count sample into its respective tissue. If there are more than two tissues, we call this a multi-class classification problem.\n",
        "\n",
        "In the first part of the exercise you will have to choose an appropiate loss function for the problem and you will have to decide upon the neural network architecture. In the first 10 minutes of the exercise, discuss the following with your classmates:\n",
        "\n",
        "1. Which loss function should you use? (Hint: Look at Anders lecture for a recap)\n",
        "2. Which architecture (number of layers and hidden units) should you use?\n",
        "\n",
        "\n",
        "\n",
        "After you made yourself familiar with the code and the results from our first trained neural network, try out a few things. Here are some guiding questions:\n",
        "\n",
        "1.   Does the model overfit?\n",
        "2.   Try to improve the model by changing model width (number of hidden units in a layer) and model depth (number of hidden layers).\n",
        "3.   Does the model perform equally well for all tissues?\n",
        "\n",
        "\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "JJZa15x-8Xw9"
      },
      "source": [
        "## The data"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "HwXHYUp4S_KA"
      },
      "source": [
        "### Loading the data"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "4nvjR93IJ9Fe"
      },
      "source": [
        "We have prepared a small subset from the Genotype-Tissue Expression (GTEx) database. It contains 3382 samples with expression counts from 1000 genes for 7 different tissues. You can see the raw counts below."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "H3RqFaNf8Xw8",
        "outputId": "ce8f655c-bbef-4c79-d526-6a06ead63dc1"
      },
      "outputs": [],
      "source": [
        "!wget https://zenodo.org/record/7822717/files/gtex_1000.csv.gz"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 574
        },
        "id": "PbZASIo98Xw-",
        "outputId": "8b075349-84ff-4908-e867-bd4d5bc4ef13"
      },
      "outputs": [],
      "source": [
        "import pandas as pd\n",
        "\n",
        "gtex_raw = pd.read_table(\"gtex_1000.csv.gz\",sep='\\t')\n",
        "gtex_raw"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "kycfVDJS8XxA"
      },
      "source": [
        "As you can see in the las column of the dataframe, the tissues are enconded as strings, and pytorch does not support this. To address so, we will simply encode the tissues as numbers. Luckily, sklearn has a function for that! See the before and after below."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "jwUgoUEX8XxA",
        "outputId": "e6690070-58d6-4d8e-d487-f60a88a962c1"
      },
      "outputs": [],
      "source": [
        "# tissues before transform\n",
        "print('tissues before \\n', gtex_raw['tissue'])\n",
        "\n",
        "# transforming strings to integer labels\n",
        "gtex = gtex_raw.copy()\n",
        "#Create a set of unique labels\n",
        "label_set = set(gtex['tissue'])  \n",
        "# Create a dictionary mapping labels to numbers\n",
        "label_dict = {label: i for i, label in enumerate(label_set)}  \n",
        "# Encode each label using the dictionary\n",
        "gtex['tissue'] = [label_dict[label] for label in gtex['tissue']]\n",
        "# tissues after transform\n",
        "print('tissues after \\n', gtex['tissue'])"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "GdkJZlU48XxB"
      },
      "source": [
        "### An initial exploration of the data using PCA"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "0XTGTirk8XxB"
      },
      "source": [
        "Before diving \"deep\" into a neural network, we will explore the data a bit using a principal component analysis (PCA). This will serve as a good baseline to see how easy it is to (linearly) separate the tissues from the data. In this way, we can have an initial guess how complex our model should be.\n",
        "\n",
        "Due to the extreme range of the expression counts, we log-transform the data before visualizing it in a PCA to get a \"clearer visual\"."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "BzL11DkG8XxB"
      },
      "outputs": [],
      "source": [
        "import numpy as np\n",
        "from sklearn.decomposition import PCA\n",
        "\n",
        "# log scale the data\n",
        "x = gtex.loc[:, gtex.columns.drop('tissue')].values\n",
        "x = np.log(x+0.01) # add pseudo-counts\n",
        "\n",
        "# fit pca \n",
        "pca = PCA(n_components=2)\n",
        "principal_components = pca.fit_transform(x)\n",
        "df_pca = pd.DataFrame(data = principal_components\n",
        "             , columns = ['principal component 1', 'principal component 2'])"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 632
        },
        "id": "q62AGGKJXRou",
        "outputId": "954866fd-276d-4e81-b5ed-7b18a37d6814"
      },
      "outputs": [],
      "source": [
        "import matplotlib.pyplot as plt\n",
        "import seaborn as sns\n",
        "sns.set(font_scale = 2)\n",
        "\n",
        "# add the original tissue strings to the PCA data\n",
        "df_pca['tissue'] = gtex_raw['tissue'].values\n",
        "\n",
        "fig, ax = plt.subplots(figsize=(10,10))\n",
        "sns.scatterplot(data=df_pca, x='principal component 1', y='principal component 2', hue='tissue', ax=ax)\n",
        "ax.legend(bbox_to_anchor=(1.05, 1.), loc='upper left')\n",
        "plt.show()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "PHqKJx5Y8XxC"
      },
      "source": [
        "From the PCA, it becomes apparent that tissues have different gene expression profiles. Let's try to create a neural network that classifies the samples given the gene expression values."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "EmQAW0Fg8XxD"
      },
      "source": [
        "### Preliminaries: preparing the data for pytorch"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "VtKwAHoj8XxD"
      },
      "source": [
        "Pytorch has certain requirements that datasets need to fulfill. We don't have to bother with the details, though. The Dataset module takes care of things for us, we just need to use it to specify what the input and what the output is, and how to get samples from it.\n",
        "\n",
        "In more detail, since datasets can be of different types, we need to create new class functions for making the dataset (\\__init__\\()), telling us how many samples there are (\\__len__\\()), and accessing samples by index (\\__getitem__\\())."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "oxv_nemu8XxD"
      },
      "outputs": [],
      "source": [
        "import torch\n",
        "from torch.utils.data import Dataset, DataLoader\n",
        "\n",
        "class GeneExpressionDataset(Dataset):\n",
        "    '''\n",
        "    Creates a Dataset class for gene expression dataset\n",
        "    gene_dim is the number of genes (features)\n",
        "    The rows of the dataframe contain samples, and the \n",
        "    columns contain gene expression values \n",
        "    and the class label (tissue) at label_position.\n",
        "    '''\n",
        "    def __init__(self, gtex, label_position=-1):\n",
        "        '''\n",
        "        Args:\n",
        "            gtex: pandas dataframe containing input and output data\n",
        "            label_position: column id of the class labels\n",
        "        '''\n",
        "        self.dataframe = gtex\n",
        "        self.label_position = label_position\n",
        "\n",
        "        self.label = torch.tensor(self.dataframe.iloc[:,label_position].to_numpy())\n",
        "        self.data = torch.tensor(self.dataframe.drop(self.dataframe.columns[[self.label_position]], axis=1).values).float()\n",
        "        \n",
        "    def __len__(self):\n",
        "        return(len(self.dataframe))\n",
        "    \n",
        "    def __getitem__(self, idx):\n",
        "        # get expression and labels\n",
        "        expression = self.data[idx,:]\n",
        "        label = self.label[idx]\n",
        "        return expression, label"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "bf-8F4uoJvPl"
      },
      "source": [
        "## The neural network"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "GsftRoRl8XxE"
      },
      "source": [
        "Once we have defined our dataset, we need to define our neural network. "
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "Z8pyJRroGrBi"
      },
      "outputs": [],
      "source": [
        "import torch.nn as nn\n",
        "import torch.nn.functional as F\n",
        "\n",
        "class Net(nn.Module):\n",
        "    def __init__(self):\n",
        "        super().__init__()\n",
        "        self.fc1 = nn.Linear(in_features=1000,out_features=200)\n",
        "        self.fc2 = nn.Linear(in_features=200,out_features=7)\n",
        "    def forward(self, x):\n",
        "        x = F.relu(self.fc1(x))\n",
        "        x = self.fc2(x)\n",
        "        return x"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "m0z7UmsVJ3xi"
      },
      "source": [
        "## Training"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "HAFhhQo-8XxE"
      },
      "source": [
        "Finally, we can move to the main part. We first create a train and test set from the data. This is important for knowing how the model performs on unseen samples from the same \"distribution\".\n",
        "\n",
        "We then use these two sets to make so-called dataloaders. These essentially take care of splitting the data into randomly grouped subsets (they are called mini-batches).\n",
        "\n",
        "After that we initialize the network and define the criterion, aka the loss function, and the optimizer that will update the model weights for us."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "O_ZimOkq8XxE"
      },
      "outputs": [],
      "source": [
        "from sklearn.model_selection import train_test_split\n",
        "from torch import optim\n",
        "\n",
        "train, test = train_test_split(gtex,test_size=0.2,random_state=42)\n",
        "\n",
        "gtex_train = GeneExpressionDataset(train)\n",
        "gtex_test = GeneExpressionDataset(test)\n",
        "batch_size = 32\n",
        "train_loader = torch.utils.data.DataLoader(gtex_train, batch_size=batch_size, num_workers=0, shuffle=True)\n",
        "test_loader = torch.utils.data.DataLoader(gtex_test, batch_size=batch_size, num_workers=0, shuffle=True)\n",
        "\n",
        "epochs = 20\n",
        "\n",
        "net = Net()\n",
        "criterion = nn.CrossEntropyLoss()\n",
        "optimizer = optim.Adam(net.parameters(), lr=1e-6)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {},
      "outputs": [],
      "source": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "fPBuorkhIcSQ"
      },
      "source": [
        "Now we can finally start training!"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "0ru3IrRM8XxF",
        "outputId": "1db444b6-d353-4850-8487-e431bfc5dff2"
      },
      "outputs": [],
      "source": [
        "# in these lists we store loss and accuracy per epoch to keep track of the progress\n",
        "train_loss = []\n",
        "test_loss = []\n",
        "train_acc = []\n",
        "test_acc  = []\n",
        "\n",
        "for epoch in range(epochs):\n",
        "    # for each epoch we collect the stats as averages over the mini-batches\n",
        "    epoch_train_loss = 0.0\n",
        "    epoch_test_loss = 0.0\n",
        "    epoch_train_acc = 0.0\n",
        "    epoch_test_acc = 0.0\n",
        "    \n",
        "    # train\n",
        "    for inputs,labels in train_loader:\n",
        "        \n",
        "        # zero the parameter gradients\n",
        "        optimizer.zero_grad()\n",
        "        # forward + loss\n",
        "        outputs = net(inputs)\n",
        "        loss = criterion(outputs, labels)\n",
        "        epoch_train_loss += loss.item()\n",
        "        # backward\n",
        "        loss.backward()\n",
        "        # optimize\n",
        "        optimizer.step()\n",
        "\n",
        "        # calculate train accuracy\n",
        "        correct_pred = outputs.argmax(dim=1) == labels\n",
        "        acc = correct_pred.sum() / len(correct_pred)\n",
        "        epoch_train_acc += acc.item()\n",
        "\n",
        "    # test\n",
        "    for inputs,labels in test_loader:\n",
        "\n",
        "        # forward only\n",
        "        outputs = net(inputs)\n",
        "        loss = criterion(outputs, labels)\n",
        "        epoch_test_loss += loss.item()\n",
        "        \n",
        "        correct_pred = outputs.argmax(dim=1) == labels\n",
        "        acc = correct_pred.sum() / len(correct_pred)\n",
        "        epoch_test_acc += acc.item()\n",
        "        \n",
        "    # average the stats and append to lists\n",
        "    epoch_train_loss /= len(train_loader)\n",
        "    epoch_test_loss /= len(test_loader)\n",
        "    epoch_train_acc /= len(train_loader)\n",
        "    epoch_test_acc /= len(test_loader)\n",
        "    train_loss.append(epoch_train_loss)\n",
        "    test_loss.append(epoch_test_loss)\n",
        "    train_acc.append(epoch_train_acc)\n",
        "    test_acc.append(epoch_test_acc)\n",
        "                                \n",
        "    print(\"Epoch:\", epoch)\n",
        "    print(\"train loss:\", round(epoch_train_loss, 4),\"test loss:\", round(epoch_test_loss, 4))\n",
        "    print(\"train acc:\", round(epoch_train_acc, 4),\"test accuracy:\", round(epoch_test_acc, 4))"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "-vI6x1iC8XxF"
      },
      "source": [
        "Lets now plot the loss curves"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 306
        },
        "id": "ccsv71UE8XxF",
        "outputId": "42ea86f8-bdff-4137-ccfc-93e4a7287916"
      },
      "outputs": [],
      "source": [
        "plt.plot(np.arange(epoch+1),train_loss,label=\"Train\")\n",
        "plt.plot(np.arange(epoch+1),test_loss,label=\"Test\")\n",
        "plt.legend()\n",
        "plt.xlabel(\"epoch\")\n",
        "plt.ylabel(\"Loss\")\n",
        "plt.show()"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 306
        },
        "id": "AnRW4gbO8XxG",
        "outputId": "a187eb35-71f3-49a1-9b4a-de02c3cd0610"
      },
      "outputs": [],
      "source": [
        "plt.plot(np.arange(epoch+1),train_acc,label=\"Train\")\n",
        "plt.plot(np.arange(epoch+1),test_acc,label=\"Test\")\n",
        "plt.legend()\n",
        "plt.xlabel(\"epoch\")\n",
        "plt.ylabel(\"Accuracy\")\n",
        "plt.show()"
      ]
    }
  ],
  "metadata": {
    "colab": {
      "collapsed_sections": [
        "JJZa15x-8Xw9",
        "HwXHYUp4S_KA",
        "GdkJZlU48XxB",
        "EmQAW0Fg8XxD",
        "bf-8F4uoJvPl"
      ],
      "name": "tissue_classification.ipynb",
      "provenance": []
    },
    "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.8.5"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 0
}