notebooks/05_clustering-analysis.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Clustering analysis"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-04-16T15:57:28.080865Z",
     "start_time": "2020-04-16T15:57:28.076022Z"
    }
   },
   "outputs": [],
   "source": [
    "import os, sys\n",
    "import warnings\n",
    "\n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "\n",
    "from scipy.stats import zscore\n",
    "import scipy.cluster.hierarchy as sch\n",
    "import matplotlib.pyplot as plt\n",
    "from sklearn.cluster import AgglomerativeClustering \n",
    "\n",
    "from tslearn.clustering import TimeSeriesKMeans\n",
    "from tslearn.preprocessing import TimeSeriesScalerMeanVariance, TimeSeriesResampler"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pd.set_option(\"display.max_rows\", 100)\n",
    "pd.set_option(\"display.max_columns\", 100)\n",
    "warnings.filterwarnings('ignore')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%load_ext autoreload\n",
    "%autoreload 2\n",
    "\n",
    "%load_ext watermark\n",
    "%watermark -a 'Fernando Pozo' -u -n -t -z -g -p pandas,numpy\n",
    "\n",
    "!pwd"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def process_logratios(df):\n",
    "    df = df[df.columns[4:]].set_index('peptide')\n",
    "    df['std'] = df.std(axis=1)\n",
    "#     df = df.sort_values('std',ascending = False).drop('std', axis=1).head(int(df.shape[0]/4))\n",
    "    df = df.sort_values('std',ascending = False).drop('std', axis=1)\n",
    "    df_scaled = pd.DataFrame(zscore(df.values, axis=0), index=df.index, columns=df.columns)\n",
    "    return df_scaled\n",
    "\n",
    "\n",
    "def time_series_kmeans_plot(df:list, n_clusters:int, title:str, seed:int=123):\n",
    "    '''\n",
    "    https://tslearn.readthedocs.io/en/stable/auto_examples/clustering/plot_kmeans.html#sphx-glr-auto-examples-clustering-plot-kmeans-py\n",
    "    '''\n",
    "    data = df.values\n",
    "    labels = list(df.columns)\n",
    "    clusters = n_clusters\n",
    "    \n",
    "    X_train = TimeSeriesScalerMeanVariance().fit_transform(data)\n",
    "#     X_train = TimeSeriesResampler(sz=3).fit_transform(X_train)\n",
    "    sz = X_train.shape[1]\n",
    "    km = TimeSeriesKMeans(n_clusters=clusters, metric=\"dtw\", max_iter=10, verbose=True, random_state=seed)\n",
    "    y_pred = km.fit_predict(X_train)\n",
    "\n",
    "    fig, axes = plt.subplots(1, 1, sharex=True, sharey=True, figsize=(20,8), \n",
    "    #                              dpi=1200\n",
    "                                )\n",
    "    fig.suptitle(title, fontsize=16, y=0.98)\n",
    "    for yi in range(clusters):\n",
    "        plt.subplot(3, clusters, yi + 1)\n",
    "        for xx in X_train[y_pred == yi]:\n",
    "            plt.plot(xx.ravel(), \"k-\", alpha=.02)\n",
    "            plt.xticks(np.linspace(0, sz, sz+1))\n",
    "        plt.plot(km.cluster_centers_[yi].ravel(), \"r-\")\n",
    "    #     plt.xlim(0, sz)\n",
    "        plt.ylim(-3, 3)\n",
    "        plt.text(0.55, 0.85,'Cluster %d' % (yi + 1),\n",
    "                 transform=plt.gca().transAxes)\n",
    "    fig.text(0.5, 0.6, 'Timepoints', ha='center', fontsize=14)\n",
    "    fig.text(0.09, 0.77, 'Normalized log ratio', va='center', rotation='vertical', fontsize=14)\n",
    "    return fig, y_pred"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 1. Data Loading"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_wt = pd.read_excel('../data/external/files/September_2018_Julio_WT_TC_1-48_rep_replicate3_20180921s_with_LP.xlsx')\n",
    "df_ko = pd.read_excel('../data/external/files/Julio_KO_PO4_1-30_20170921_Final_PEPTIDE.xlsx')\n",
    "\n",
    "df_wt_log = df_wt[['ptm', 'proteinName', 'proteinDescription', 'FTestQValue', 'peptide',\n",
    "       'log2ratio_0-0', 'log2ratio_1c', 'log2ratio_2c', 'log2ratio_5c', \n",
    "       'log2ratio_10c', 'log2ratio_20c', 'log2ratio_40c', 'log2ratio_80c']]\n",
    "df_wt_log = df_wt_log[df_wt_log['FTestQValue']>0.05]\n",
    "\n",
    "df_ko_log = df_ko[['ptm', 'proteinName', 'proteinDescription', 'FTestQValue', 'peptide',\n",
    "       'log2ratio_KO-0', 'log2ratio_KOc-1', 'log2ratio_KOc-2',\n",
    "       'log2ratio_KOc-5', 'log2ratio_KOc-10', 'log2ratio_KO-20']]\n",
    "df_ko_log = df_ko_log[df_ko_log['FTestQValue']>0.05]\n",
    "\n",
    "df_wt_log_data = process_logratios(df_wt_log)\n",
    "df_wt_log_data.columns = [0, 1, 2, 5, 10, 20, 40, 80]\n",
    "\n",
    "df_ko_log_data = process_logratios(df_ko_log)\n",
    "df_ko_log_data.columns = [0, 1, 2, 5, 10, 20]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 2. Data Exploration"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Exploring the best set of clusters with dendrogram with ward (minimizing variances between them)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data = df_wt_log_data\n",
    "plt.figure(figsize=(20,6))\n",
    "dendrogram = sch.dendrogram(sch.linkage(data, method  = \"ward\"))\n",
    "plt.title('Wild-type peptides dendrogram')\n",
    "plt.xlabel(f'Set of {data.shape[0]} peptides')\n",
    "plt.ylabel('Euclidean distances')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data = df_ko_log_data\n",
    "plt.figure(figsize=(20,6))\n",
    "sch.dendrogram(sch.linkage(data, method  = \"ward\"))\n",
    "plt.title('Knock-out peptides dendrogram')\n",
    "plt.xlabel(f'Set of {data.shape[0]} peptides', fontsize=16)\n",
    "plt.ylabel('Euclidean distances', fontsize=16)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig_wt, y_pred_wt = time_series_kmeans_plot(df_wt_log_data, 3, 'A')\n",
    "df_wt_log_data['cluster_kmeans'] = y_pred_wt + 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig_ko, y_pred_ko = time_series_kmeans_plot(df_ko_log_data, 3, 'B')\n",
    "df_ko_log_data['cluster_kmeans'] = y_pred_ko + 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig_wt.savefig('../reports/figures/figure4A.png', bbox_inches='tight')\n",
    "fig_ko.savefig('../reports/figures/figure4B.png', bbox_inches='tight')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_wt_log_data['cluster'] = AgglomerativeClustering(n_clusters = 6, affinity = 'euclidean', linkage ='ward').fit_predict(df_wt_log_data) + 1\n",
    "df_ko_log_data['cluster'] = AgglomerativeClustering(n_clusters = 6, affinity = 'euclidean', linkage ='ward').fit_predict(df_ko_log_data) + 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_clusters(df, title):\n",
    "    rows = 3\n",
    "    cols = 2\n",
    "    timepoints = df.drop('cluster', axis=1).columns\n",
    "    fig, axes = plt.subplots(rows, cols, sharex=True, sharey=True, figsize=(12,8), \n",
    "                             dpi=1200\n",
    "                            )\n",
    "    fig.suptitle(title, fontsize=20, y=1.05)\n",
    "\n",
    "    for n, ax in enumerate(axes.flat):\n",
    "        n=n+1\n",
    "        data = df[df['cluster'] == n].drop('cluster', axis=1)\n",
    "        for i in range(0, data.shape[0]):\n",
    "            ax.plot(timepoints, data.iloc[i].values, marker='', linewidth=0.5, alpha=0.4, color='orange')\n",
    "            ax.plot(timepoints, data.describe().iloc[1].values, marker='.', markersize=12, \n",
    "                    linestyle='solid', linewidth=0.4, alpha=0.4, color='gray')\n",
    "            ax.set_ylim([-5, 5])\n",
    "            ax.set_title(f'Cluster {n}')\n",
    "    plt.tight_layout()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_clusters(df_wt_log_data, 'Wild-type peptides cluster analysis')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_clusters(df_ko_log_data, 'Knock-out peptides cluster analysis')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "cell_style": "center"
   },
   "source": [
    "# Project contribution <a class=\"anchor\" id=\"project-contribution\"></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Author information**:\n",
    "Fernando Pozo ([@fpozoc](https://gitlab.com/fpozoc))\n",
    "\n",
    "You can find the data driven project jupyter notebook template [here](https://gitlab.com/fpozoc/data-driven-project-template/-/blob/master/notebooks/1.0-nb_template.ipynb)."
   ]
  }
 ],
 "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.7.10"
  },
  "varInspector": {
   "cols": {
    "lenName": 16,
    "lenType": 16,
    "lenVar": 40
   },
   "kernels_config": {
    "python": {
     "delete_cmd_postfix": "",
     "delete_cmd_prefix": "del ",
     "library": "var_list.py",
     "varRefreshCmd": "print(var_dic_list())"
    },
    "r": {
     "delete_cmd_postfix": ") ",
     "delete_cmd_prefix": "rm(",
     "library": "var_list.r",
     "varRefreshCmd": "cat(var_dic_list()) "
    }
   },
   "types_to_exclude": [
    "module",
    "function",
    "builtin_function_or_method",
    "instance",
    "_Feature"
   ],
   "window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}