Create_Plot_of_CVs.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "\n",
    "import random\n",
    "import sys, os\n",
    "from itertools import combinations\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sns\n",
    "\n",
    "from sklearn.model_selection import train_test_split"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "############################ \n",
    "#Recording coefficient of variation for abundances\n",
    "def recordCV(data_df, coeff_df, random_run_count = 0):\n",
    "    \n",
    "    #We will record the CV scores \n",
    "    \n",
    "    #############################\n",
    "    #Select peptides occuring on each protein\n",
    "    all_proteins = np.unique(data_df['Protein'].values)\n",
    "    print(\"Number of proteins: \", len(all_proteins))\n",
    "    \n",
    "    run_start_index = 11\n",
    "    \n",
    "    all_raw_CVs = []\n",
    "    all_adj_CVs = []\n",
    "    protein_sizes = []\n",
    "    \n",
    "    #For each protein, calculate CV before and after adjustment\n",
    "    for protein in all_proteins:\n",
    "        \n",
    "        sub_df = data_df[data_df['Protein'] == protein]\n",
    "        sub_peptide_indices =  [np.where(data_df.index == i)[0][0] for i in sub_df.index]\n",
    "\n",
    "        protein_sizes.append(len(sub_peptide_indices))\n",
    "        \n",
    "        all_run_raw_CVs = []\n",
    "        all_run_adj_CVs = []\n",
    "        \n",
    "        #Calculate for each run\n",
    "        for random_experiment in range(run_start_index, data_df.shape[1]):\n",
    "            \n",
    "            #For each peptide, record the quantities\n",
    "            quantities = data_df.iloc[sub_peptide_indices, random_experiment]\n",
    "            #print(quantities)\n",
    "            \n",
    "            #Reject sample if any values are nan\n",
    "            if np.all(~np.isnan(quantities.values.ravel())) and \\\n",
    "                   np.all(quantities.values.ravel() != 0.0):\n",
    "\n",
    "                #Calculate the abundances by dividing quantities to peptide coefficients\n",
    "                sibling_coefficients = coeff_df.iloc[sub_peptide_indices, 0]\n",
    "#                 print(\"Sibling intensities: \", quantities)\n",
    "#                 print(\"Sibling coeffs: \", sibling_coefficients)\n",
    "\n",
    "                sibling_abundances = quantities.values.ravel() / sibling_coefficients.values.ravel()\n",
    "#                 print(\"Sibling abundances: \", sibling_abundances)\n",
    "\n",
    "                #Record final CVs\n",
    "                raw_CV = np.std(quantities.values.ravel()) / np.mean(quantities.values.ravel())\n",
    "                adj_CV = np.std(sibling_abundances) / np.mean(sibling_abundances)\n",
    "                \n",
    "                all_run_raw_CVs.append(raw_CV)\n",
    "                all_run_adj_CVs.append(adj_CV)\n",
    "        \n",
    "#                 print(\"Raw run CVs \", np.mean(all_run_raw_CVs))\n",
    "#                 print(\"Adj run CVs \", np.mean(all_run_adj_CVs))\n",
    "    \n",
    "        all_raw_CVs.append(all_run_raw_CVs)\n",
    "        all_adj_CVs.append(all_run_adj_CVs)\n",
    "\n",
    "    all_raw_CVs = np.array(all_raw_CVs)\n",
    "    all_adj_CVs = np.array(all_adj_CVs)\n",
    "\n",
    "#     print(\"Raw CVs \", np.mean(all_raw_CVs))\n",
    "#     print(\"Adj CVs \", np.mean(all_adj_CVs))\n",
    "    \n",
    "    return [all_raw_CVs, all_adj_CVs, protein_sizes]\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "n_runs = 120\n",
    "seq_length = 60\n",
    "\n",
    "#Read dataset\n",
    "data_df = pd.read_csv('preprocess_datasets/preprocessed_datasets/2019_guo_nci60_formatted_peptide_quants.tsv', \n",
    "                      sep = '\\t', index_col = 0)\n",
    "\n",
    "print(\"Peptide df \", data_df.shape)\n",
    "print(\"Peptide df \", data_df.head())\n",
    "\n",
    "#Input 1 is peptide intensity measurements\n",
    "#Also record intensity measurements for pairs\n",
    "q_df = data_df.iloc[:, -n_runs:]\n",
    "\n",
    "#Normalize the intensities such that the sum of elements in each column is equal\n",
    "X = q_df.values\n",
    "X = (X / X.sum(axis=0, keepdims=1)) * X.shape[0]\n",
    "q_df = pd.DataFrame(X, index = q_df.index, columns = q_df.columns)\n",
    "data_df.iloc[:, -n_runs:] = q_df\n",
    "\n",
    "#Input 2 is protein mappings\n",
    "#Convert protein labels to int values\n",
    "protein_labels = data_df['Protein'].values\n",
    "unique_proteins = np.unique(protein_labels)\n",
    "n_proteins = len(unique_proteins)\n",
    "print(\"Number of unique proteins \", n_proteins)\n",
    "int_protein_labels = [np.where(protein_labels[i] == unique_proteins)[0][0] for i in range(protein_labels.shape[0])]\n",
    "int_protein_labels = np.asarray(int_protein_labels)\n",
    "print(\"Protein labels \", int_protein_labels)\n",
    "n_peptides = data_df.shape[0]\n",
    "\n",
    "print(\"No of peptides: \", n_peptides)\n",
    "print(\"No of proteins: \", n_proteins)\n",
    "print(\"No of runs: \", n_runs)\n",
    "\n",
    "#Split the proteins into train/validation/test sets\n",
    "\n",
    "train_proteins, test_proteins = train_test_split((np.arange(len(np.unique(protein_labels)))), \n",
    "                                   test_size=0.2, random_state=12345)\n",
    "\n",
    "#Define train/validation/test peptide pairs\n",
    "train_peptides = np.concatenate([list(np.where(protein_labels == np.unique(protein_labels)[p])[0]) for p in train_proteins])\n",
    "test_peptides = np.concatenate([list(np.where(protein_labels == np.unique(protein_labels)[p])[0]) for p in test_proteins])\n",
    "\n",
    "print(\"No of train/test proteins: %d/%d\" % (len(train_proteins), len(test_proteins)))\n",
    "print(\"No of train/test peptides: %d/%d\" % (len(train_peptides), len(test_peptides)))\n",
    "\n",
    "#Split the runs into train/validation/test sets\n",
    "#Modified code for replicate samples\n",
    "train_runs, test_runs = train_test_split((np.arange(q_df.shape[1] / 2)), \n",
    "                                   test_size=0.2, random_state=12345)\n",
    "\n",
    "train_runs = np.array([[2*i, 2*i+1] for i in train_runs]).astype(int).ravel()\n",
    "test_runs = np.array([[2*i, 2*i+1] for i in test_runs]).astype(int).ravel()\n",
    "\n",
    "print(\"No of train runs \", len(train_runs))\n",
    "print(\"No of test runs \", len(test_runs))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Split the data into training and test set\n",
    "data_df_train = data_df.iloc[train_peptides]\n",
    "data_df_test = data_df.iloc[test_peptides]\n",
    "\n",
    "data_df_train = pd.concat([data_df_train.iloc[:, :-n_runs], data_df_train.iloc[:, data_df.shape[1] - n_runs + train_runs]], axis = 1)\n",
    "data_df_test = pd.concat([data_df_test.iloc[:, :-n_runs], data_df_test.iloc[:, data_df.shape[1] - n_runs + test_runs]], axis = 1)\n",
    "\n",
    "data_df_test"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Define function to create scatter plots\n",
    "#Note that this scatter plot is for peptide ratios\n",
    "def createPlots(scores, n_bins = 100):\n",
    "    \n",
    "    #Create comparison scatter plots\n",
    "    fig, ax = plt.subplots()\n",
    "    fig.set_size_inches(20, 20)\n",
    "\n",
    "    SMALL_SIZE = 40\n",
    "    MEDIUM_SIZE = 50\n",
    "    BIGGER_SIZE = 70\n",
    "\n",
    "    plt.rc('font', size=SMALL_SIZE)          # controls default text sizes\n",
    "    plt.rc('axes', titlesize=MEDIUM_SIZE)     # fontsize of the axes title\n",
    "    plt.rc('axes', labelsize=MEDIUM_SIZE)    # fontsize of the x and y labels\n",
    "    plt.rc('xtick', labelsize=SMALL_SIZE)    # fontsize of the tick labels\n",
    "    plt.rc('ytick', labelsize=SMALL_SIZE)    # fontsize of the tick labels\n",
    "    plt.rc('legend', fontsize=SMALL_SIZE)    # legend fontsize\n",
    "    plt.rc('figure', titlesize=BIGGER_SIZE)  # fontsize of the figure title\n",
    "    \n",
    "    print(\"Creating a plot for all proteins...\")\n",
    "    raw_scores = scores[0].mean(axis = 1)\n",
    "    adj_scores = scores[1].mean(axis = 1) \n",
    "    protein_sizes = scores[2]\n",
    "    \n",
    "    count_raw = 0\n",
    "    count_adj = 0\n",
    "    for i in range(raw_scores.shape[0]):\n",
    "        if raw_scores[i] > adj_scores[i]:\n",
    "            plt.scatter(raw_scores[i], \n",
    "                        adj_scores[i], \n",
    "                        s = 500, alpha = 0.5, color = \"#22a6b3\")\n",
    "            count_raw += 1\n",
    "\n",
    "        elif raw_scores[i] < adj_scores[i]:\n",
    "            plt.scatter(raw_scores[i], \n",
    "                        adj_scores[i], \n",
    "                        s = 500, alpha = 0.5, color = '#22a6b3')\n",
    "            count_adj += 1\n",
    "            \n",
    "    print(\"Count raw \", count_raw)\n",
    "    print(\"Count adj \", count_adj)\n",
    "\n",
    "    plt.plot([0.0, 2.0], [0.0, 2.0], '--', lw = 5, color = 'Black', zorder = -1)\n",
    "    plt.xlim([0.0, 2.0])\n",
    "    plt.ylim([0.0, 2.0])\n",
    "    \n",
    "    plt.xticks([0.0, 0.5, 1.0, 1.5, 2.0])\n",
    "    plt.yticks([0.0, 0.5, 1.0, 1.5, 2.0])\n",
    "    \n",
    "    plt.xlabel('Coefficient of variation for observed abundances')\n",
    "    plt.ylabel('Coefficient of variation for adjusted abundances')\n",
    "    plt.grid()\n",
    "    plt.show()\n",
    "   \n",
    "    print(\"Creating a plot for proteins with 5+ peptides...\")\n",
    "    \n",
    "    #Create comparison scatter plots\n",
    "    fig, ax = plt.subplots()\n",
    "    fig.set_size_inches(20, 20)\n",
    "\n",
    "    SMALL_SIZE = 40\n",
    "    MEDIUM_SIZE = 50\n",
    "    BIGGER_SIZE = 70\n",
    "\n",
    "    plt.rc('font', size=SMALL_SIZE)          # controls default text sizes\n",
    "    plt.rc('axes', titlesize=MEDIUM_SIZE)     # fontsize of the axes title\n",
    "    plt.rc('axes', labelsize=MEDIUM_SIZE)    # fontsize of the x and y labels\n",
    "    plt.rc('xtick', labelsize=SMALL_SIZE)    # fontsize of the tick labels\n",
    "    plt.rc('ytick', labelsize=SMALL_SIZE)    # fontsize of the tick labels\n",
    "    plt.rc('legend', fontsize=SMALL_SIZE)    # legend fontsize\n",
    "    plt.rc('figure', titlesize=BIGGER_SIZE)  # fontsize of the figure title\n",
    "    \n",
    "    raw_scores = scores[0].mean(axis = 1)\n",
    "    adj_scores = scores[1].mean(axis = 1) \n",
    "    protein_sizes = scores[2]\n",
    "    \n",
    "    count_raw = 0\n",
    "    count_adj = 0\n",
    "    for i in range(raw_scores.shape[0]):\n",
    "        if protein_sizes[i] > 5:\n",
    "            if raw_scores[i] > adj_scores[i]:\n",
    "                plt.scatter(raw_scores[i], \n",
    "                            adj_scores[i], \n",
    "                            s = 500, alpha = 0.5, color = \"#22a6b3\")\n",
    "                count_raw += 1\n",
    "\n",
    "            elif raw_scores[i] < adj_scores[i]:\n",
    "                plt.scatter(raw_scores[i], \n",
    "                            adj_scores[i], \n",
    "                            s = 500, alpha = 0.5, color = '#22a6b3')\n",
    "                count_adj += 1\n",
    "            \n",
    "    print(\"Count raw \", count_raw)\n",
    "    print(\"Count adj \", count_adj)\n",
    "\n",
    "    plt.plot([0.0, 2.0], [0.0, 2.0], '--', lw = 5, color = 'Black', zorder = -1)\n",
    "    plt.xlim([0.0, 2.0])\n",
    "    plt.ylim([0.0, 2.0])\n",
    "    \n",
    "    plt.xticks([0.0, 0.5, 1.0, 1.5, 2.0])\n",
    "    plt.yticks([0.0, 0.5, 1.0, 1.5, 2.0])\n",
    "    \n",
    "    plt.xlabel('Coefficient of variation for observed abundances')\n",
    "    plt.ylabel('Coefficient of variation for adjusted abundances')\n",
    "    plt.grid()\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "#Read training coefficients\n",
    "coeff_df = pd.read_csv('trained_models/2019_guo_nci60/2019_guo_nci60_inferred_coefficients.tsv', sep = '\\t', index_col = 0)\n",
    "print(\"Coefficients \", coeff_df.shape)\n",
    "coeff_df = coeff_df.abs()\n",
    "coeff_df.sort_values(by = '0')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "#Create test plots\n",
    "results_test = recordCV(data_df_test, coeff_df.loc[data_df_test.index])\n",
    "createPlots(results_test)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Define function to create histogram of difference in CV\n",
    "def createPlots(scores, n_bins = 25):\n",
    "    \n",
    "    #Create comparison scatter plots\n",
    "    fig, ax = plt.subplots()\n",
    "    fig.set_size_inches(20, 20)\n",
    "\n",
    "    SMALL_SIZE = 40\n",
    "    MEDIUM_SIZE = 50\n",
    "    BIGGER_SIZE = 70\n",
    "\n",
    "    plt.rc('font', size=SMALL_SIZE)          # controls default text sizes\n",
    "    plt.rc('axes', titlesize=MEDIUM_SIZE)     # fontsize of the axes title\n",
    "    plt.rc('axes', labelsize=MEDIUM_SIZE)    # fontsize of the x and y labels\n",
    "    plt.rc('xtick', labelsize=SMALL_SIZE)    # fontsize of the tick labels\n",
    "    plt.rc('ytick', labelsize=SMALL_SIZE)    # fontsize of the tick labels\n",
    "    plt.rc('legend', fontsize=SMALL_SIZE)    # legend fontsize\n",
    "    plt.rc('figure', titlesize=BIGGER_SIZE)  # fontsize of the figure title\n",
    "    \n",
    "    sns.distplot(raw_scores - adj_scores, bins = n_bins, \n",
    "                 kde_kws={\"lw\": 10, \"color\": \"blue\", \"alpha\": 1, },\n",
    "                 hist_kws={\"linewidth\": 1, \"color\": \"blue\"})\n",
    "    plt.vlines(0, ymin = 0, ymax = 6, linewidth = 5, color = 'gray', zorder = -1)\n",
    "    plt.xlabel('Difference in coefficient of variation (observed - adjusted)')\n",
    "    plt.ylim([0, 5.5])\n",
    "    plt.show()\n",
    "   "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "results_test = recordCV(data_df_test, coeff_df.loc[data_df_test.index])\n",
    "createPlots(results_test)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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.6.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}