Create_Histograms_of_Sibling_Abundance_Differences.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\n",
    "from matplotlib.colors import LinearSegmentedColormap\n",
    "from random import sample\n",
    "\n",
    "from scipy.stats import kurtosis"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "############################ \n",
    "#Recording sibling peptide intensities and adjusted abundances\n",
    "\n",
    "#Suppressing print statements\n",
    "class HiddenPrints:\n",
    "    \n",
    "    def __enter__(self):\n",
    "        self._original_stdout = sys.stdout\n",
    "        sys.stdout = open(os.devnull, 'w')\n",
    "\n",
    "    def __exit__(self, exc_type, exc_val, exc_tb):\n",
    "        sys.stdout.close()\n",
    "        sys.stdout = self._original_stdout\n",
    "\n",
    "#Super function to sample sibling peptides\n",
    "def randomSiblingRatios(data_df, coeff_df, random_run_count = 0, supress_print = True):\n",
    "    #supresses prints\n",
    "    if not supress_print:\n",
    "        return randomSiblingRatios_sub(data_df, coeff_df, random_run_count)\n",
    "    \n",
    "    #enables all prints for debugging\n",
    "    else:\n",
    "        with HiddenPrints():\n",
    "            return randomSiblingRatios_sub(data_df, coeff_df, random_run_count)\n",
    "        \n",
    "#Function to sample sibling peptides\n",
    "def randomSiblingRatios_sub(data_df, coeff_df, random_run_count = 0):\n",
    "    \n",
    "    #We will record the scores for all peptide pairs for randomly sampled runs\n",
    "    \n",
    "    #############################\n",
    "    #Define all sibling pairs\n",
    "    all_proteins = np.unique(data_df['Protein'].values)\n",
    "    print(\"Number of proteins: \", len(all_proteins))\n",
    "    \n",
    "    #For all proteins, record indices of all sibling peptides\n",
    "    sibling_pairs = []\n",
    "\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",
    "        for i in range(len(sub_peptide_indices)):\n",
    "            for j in range(i + 1, len(sub_peptide_indices)):\n",
    "                pair = [sub_peptide_indices[i], sub_peptide_indices[j]]\n",
    "                sibling_pairs.append(pair)\n",
    "\n",
    "    #Record all sibling pairs\n",
    "    sibling_pairs = np.asarray(sibling_pairs)\n",
    "    print(\"Number of peptide pairs \", sibling_pairs.shape[0])\n",
    "\n",
    "    #############################\n",
    "\n",
    "    run_start_index = 11\n",
    "    #Define all random runs to sample intensities from\n",
    "    random_runs = list(np.arange(run_start_index, data_df.shape[1]))\n",
    "    print(\"Total number of random runs \", len(random_runs))\n",
    "\n",
    "    #Record observed intensities for the peptides\n",
    "    all_sibling_intensitites = []\n",
    "    \n",
    "    #Record adjusted abundances for the peptides\n",
    "    all_sibling_abundances = []\n",
    "    \n",
    "    if random_run_count > 0:\n",
    "        sample_count = random_run_count\n",
    "    else:\n",
    "        sample_count = len(sibling_pairs)\n",
    "    \n",
    "    for n in range(sample_count):\n",
    "        print(\"----------\")\n",
    "        print(data_df.iloc[sibling_pairs[n]]['Protein'])\n",
    "                   \n",
    "        #for random_experiment in range(run_start_index, data_df.shape[1]):\n",
    "        random_experiment = sample(range(run_start_index, data_df.shape[1]), 1)\n",
    "        #For each peptide, record the quantities\n",
    "        quantities = data_df.iloc[sibling_pairs[n], random_experiment]\n",
    "        print(quantities)\n",
    "        print(sibling_pairs[n])\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[sibling_pairs[n], 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 results\n",
    "            all_sibling_intensitites.append(quantities.values.ravel())\n",
    "            all_sibling_abundances.append(list(sibling_abundances))\n",
    "\n",
    "            #Record the other pair as well\n",
    "            all_sibling_intensitites.append(quantities.values.ravel()[::-1])\n",
    "            all_sibling_abundances.append(list(sibling_abundances)[::-1])\n",
    "\n",
    "    all_sibling_intensitites = np.array(all_sibling_intensitites)\n",
    "    all_sibling_intensitites = all_sibling_intensitites.reshape((len(all_sibling_intensitites), 2))\n",
    "\n",
    "    all_sibling_abundances = np.array(all_sibling_abundances)\n",
    "    all_sibling_abundances = all_sibling_abundances.reshape((len(all_sibling_abundances), 2))\n",
    "\n",
    "    print(\"Mean abundance: \", np.mean(all_sibling_abundances))\n",
    "\n",
    "    return [all_sibling_intensitites, all_sibling_abundances]"
   ]
  },
  {
   "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",
    "print(\"Quants before normalization \", X.sum(axis = 0))\n",
    "X = (X / X.sum(axis=0, keepdims=1)) * X.shape[0]\n",
    "print(\"Quants after normalization \", X.sum(axis = 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))\n"
   ]
  },
  {
   "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 createHistogram_peptide_ratios(sibling_ratios, color, n_bins = 300):\n",
    "\n",
    "    #####################\n",
    "    # Create both plots together\n",
    "    results = sibling_ratios[0]\n",
    "    \n",
    "    N = 50\n",
    "    fig, ax = plt.subplots()\n",
    "    fig.set_size_inches(20, 20)\n",
    "\n",
    "    SMALL_SIZE = 60\n",
    "    MEDIUM_SIZE = 80\n",
    "    BIGGER_SIZE = 90\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=MEDIUM_SIZE)    # fontsize of the tick labels\n",
    "    plt.rc('ytick', labelsize=MEDIUM_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",
    "    ratio_scores1 = pd.DataFrame(np.log10(sibling_ratios[0][:, 0] / sibling_ratios[0][:, 1]), \n",
    "                                 columns = ['Ratio'])\n",
    "    ratio_scores2 = pd.DataFrame(np.log10(sibling_ratios[1][:, 0] / sibling_ratios[1][:, 1]), \n",
    "                                 columns = ['Ratio'])\n",
    "    \n",
    "    ratio_scores1 = pd.DataFrame((sibling_ratios[0][:, 0] - sibling_ratios[0][:, 1]), \n",
    "                                 columns = ['Ratio'])\n",
    "    ratio_scores2 = pd.DataFrame((sibling_ratios[1][:, 0] - sibling_ratios[1][:, 1]), \n",
    "                                 columns = ['Ratio'])\n",
    "    \n",
    "    ratio_scores2 = ratio_scores2[~(ratio_scores2 > np.mean(ratio_scores2) + 3 * np.std(ratio_scores2))]\n",
    "    ratio_scores2 = ratio_scores2[~(ratio_scores2 < np.mean(ratio_scores2) - 3 * np.std(ratio_scores2))]\n",
    "    \n",
    "    ratio_scores1 = ratio_scores1[~(ratio_scores1 > np.mean(ratio_scores1) + 3 * np.std(ratio_scores1))]\n",
    "    ratio_scores1 = ratio_scores1[~(ratio_scores1 < np.mean(ratio_scores1) - 3 * np.std(ratio_scores1))]\n",
    "    \n",
    "    sns.distplot(ratio_scores2, bins = n_bins, \n",
    "                 kde_kws={\"lw\": 10, \"color\": \"#8854d0\", \"alpha\": 1, },\n",
    "                 hist_kws={\"linewidth\": 1, \"color\": \"#8854d0\"}, label = 'Adjusted')\n",
    "     \n",
    "    sns.distplot(ratio_scores1, bins = n_bins, \n",
    "                 kde_kws={\"lw\": 10, \"color\": color, \"alpha\": 1, },\n",
    "                 hist_kws={\"linewidth\": 1, \"color\": color}, label = 'Observed')\n",
    "    \n",
    "    print(\"Kurtosis for raw abundances \", kurtosis(ratio_scores1.values[~np.isnan(ratio_scores1.values)].ravel()))\n",
    "    print(\"Kurtosis for adjusted abundances \", kurtosis(ratio_scores2.values[~np.isnan(ratio_scores2.values)].ravel()))\n",
    "    \n",
    "    plt.xticks([-4, -2, 0, 2, 4])\n",
    "    plt.xlim([-5, 5])\n",
    "    #plt.legend()\n",
    "#     plt.xlabel('Difference between \\n sibling peptide abundances')\n",
    "#     plt.ylabel('Density')\n",
    "#     plt.show()\n",
    "    \n",
    "    \n",
    "    "
   ]
  },
  {
   "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": true
   },
   "outputs": [],
   "source": [
    "#Create training plots\n",
    "#randomSiblingRatios(data_df_train, coeff_df.loc[data_df_train.index], supress_print = False, random_run_count = 50)\n",
    "results_train = randomSiblingRatios(data_df_train, coeff_df.loc[data_df_train.index], supress_print = True)\n",
    "createHistogram_peptide_ratios(results_train, color = '#eb4d4b', n_bins = 100)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "#Create test plots\n",
    "#randomSiblingRatios(data_df_test, coeff_df.loc[data_df_test.index], supress_print = False, random_run_count = 50)\n",
    "results_test = randomSiblingRatios(data_df_test, coeff_df.loc[data_df_test.index], supress_print = True)\n",
    "createHistogram_peptide_ratios(results_test, color = '#eb4d4b', n_bins = 100)"
   ]
  },
  {
   "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
}