fig4_predictive_performance_comparison.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Predictive performance comparison\n",
    "The idea of this notebook is to take a look at the predictive performance on cell lines for all the drugs. The idea is two-fold:\n",
    "<ul>\n",
    "    <li> Assessing that the source top PVs can yield same predictive performance as a direct ridge on the source data. It would mean that the top PVs contain the relevant information for drug response prediction.\n",
    "    <li> Taking a look at which drug gets predicted using both the PV duos and the consensus representation.\n",
    "</ul>\n",
    "We here use all the cell line data for the domain adaptation. Other settings can be imagined as well."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Parameters (to change)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# None for 'rnaseq', 'fpkm' for FPKM\n",
    "type_data = 'rnaseq'\n",
    "normalization = 'TMM'\n",
    "transformation = 'log'\n",
    "mean_center = True\n",
    "std_unit = False\n",
    "filter_mytochondrial = False\n",
    "protein_coding_only = True\n",
    "d_test = [40]\n",
    "n_factors = 70\n",
    "same_pv_pca = True\n",
    "drug_file = 'input/drug_list_small.txt' # To change to drug_list.txt for full-scale analysis\n",
    "n_jobs=5"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os, sys\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sns\n",
    "import scipy\n",
    "from sklearn.model_selection import GroupKFold, GridSearchCV\n",
    "from sklearn.linear_model import ElasticNet, Ridge\n",
    "from sklearn.pipeline import Pipeline\n",
    "from sklearn.preprocessing import StandardScaler\n",
    "from sklearn.externals.joblib import Parallel, delayed\n",
    "import pickle\n",
    "plt.style.use('ggplot')\n",
    "\n",
    "#Import src implementations\n",
    "os.environ['OMP_NUM_THREADS'] = '1'\n",
    "os.environ['KMP_DUPLICATE_LIB_OK']='True'\n",
    "from data_reader.read_data import read_data\n",
    "from data_reader.read_drug_response import read_drug_response\n",
    "from data_reader.read_cna_tumors import read_cna_tumors\n",
    "from normalization_methods.feature_engineering import feature_engineering\n",
    "import precise\n",
    "from precise import DrugResponsePredictor, ConsensusRepresentation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Read all the drug from the file and load all the data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "with open(drug_file,'r') as drug_file_reader:\n",
    "    drug_file_content = drug_file_reader.read()\n",
    "    drug_file_content = drug_file_content.split('\\n')\n",
    "    drug_file_content = [e.split(',') for e in drug_file_content]\n",
    "    \n",
    "    # drug_IDs and tumor tissues are ordered in the same way\n",
    "    drug_IDs = np.array(list(zip(*drug_file_content))[0]).astype(int)\n",
    "    tumor_tissues = np.array(list(zip(*drug_file_content))[1])\n",
    "\n",
    "unique_tumor_tissues = np.unique(tumor_tissues)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "target_raw_data = dict()\n",
    "source_raw_data = dict()\n",
    "target_barcodes = dict()\n",
    "source_names = dict()\n",
    "target_data = dict()\n",
    "source_data = dict()\n",
    "source_data_filtered = dict()\n",
    "source_response_data = dict()\n",
    "source_names_filtered = dict()\n",
    "drug_names = dict()\n",
    "target_primary_site = dict()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Load cell line data\n",
    "# /!\\ Due to some mismatch in the genes available in TCGA, cell line data has to be loaded all the time\n",
    "for tissue_name in unique_tumor_tissues:\n",
    "    print(tissue_name)\n",
    "    \n",
    "    if tissue_name in target_raw_data:\n",
    "        continue\n",
    "        \n",
    "    X_target, X_source, _, s, target_names = read_data('cell_line',\n",
    "                                                       'tumor',\n",
    "                                                       'count',\n",
    "                                                        None,\n",
    "                                                        tissue_name,\n",
    "                                                        filter_mytochondrial)\n",
    "    \n",
    "    target_raw_data[tissue_name] = X_target\n",
    "    source_raw_data[tissue_name] = X_source\n",
    "    target_barcodes[tissue_name] = target_names\n",
    "    source_names[tissue_name] = s"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Normalize the data\n",
    "for tissue_name in unique_tumor_tissues:\n",
    "    print(tissue_name)\n",
    "    if tissue_name in target_data:\n",
    "        continue\n",
    "    \n",
    "    target_data[tissue_name] = feature_engineering(target_raw_data[tissue_name],\n",
    "                                                  normalization,\n",
    "                                                  transformation,\n",
    "                                                  mean_center,\n",
    "                                                  std_unit)\n",
    "    \n",
    "    # source data is not mean-centered as it will be done during cross-validation procedure.\n",
    "    source_data[tissue_name] = feature_engineering(source_raw_data[tissue_name],\n",
    "                                                  normalization,\n",
    "                                                  transformation,\n",
    "                                                  False,\n",
    "                                                  False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Normalize for variance\n",
    "for tissue_name in unique_tumor_tissues:\n",
    "    print(tissue_name)\n",
    "    if tissue_name in target_data:\n",
    "        continue\n",
    "    \n",
    "    target_total_variance = np.sqrt(np.sum(np.var(target_data[tissue_name], 0)))\n",
    "    target_data[tissue_name] = target_data[tissue_name] / target_total_variance * 10**3\n",
    "    \n",
    "    source_total_variance = np.sqrt(np.sum(np.var(source_data[tissue_name], 0)))\n",
    "    source_data[tissue_name] = source_data[tissue_name] / source_total_variance * 10**3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Read drug response\n",
    "for i, (ID, tissue) in enumerate(zip(drug_IDs, tumor_tissues)):\n",
    "    if (ID, tissue) in source_data_filtered:\n",
    "        continue\n",
    "        \n",
    "    x, y, s, name = read_drug_response(ID,\n",
    "                                       source_data[tissue],\n",
    "                                       source_names[tissue],\n",
    "                                      'count')\n",
    "    \n",
    "    source_data_filtered[(ID, tissue)] = x\n",
    "    source_response_data[(ID, tissue)] = y\n",
    "    drug_names[(ID, tissue)] = name\n",
    "    source_names_filtered[(ID, tissue)] = s"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Principal vector test\n",
    "Here we compute the predictive performance for several different drugs using either the osurce, the target of both principal vector. The latter one is still biases towards the source."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Consensus representation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "l1_ratio  = 0\n",
    "\n",
    "for ID, tissue in zip(drug_IDs, tumor_tissues):\n",
    "    print(ID, tissue)\n",
    "    \n",
    "    X_source = source_data_filtered[ID, tissue]\n",
    "    y_source = source_response_data[ID, tissue]\n",
    "    X_target = target_data[tissue]\n",
    "    \n",
    "    pickle_file = 'consensus_drug_%s_tissue_%s_l1_ratio_%s_n_factors_%s.pkl'%(ID,\n",
    "                                                                      tissue,\n",
    "                                                                      l1_ratio,\n",
    "                                                                      n_factors)\n",
    "    if pickle_file in os.listdir('./output/pred_performance/'):\n",
    "        print('%s, %s ALREADY COMPUTED'%(ID, tissue))\n",
    "        continue\n",
    "        \n",
    "    with open('./output/pred_performance/%s'%(pickle_file), 'wb') as f:\n",
    "        pickle.dump(dict(), f, pickle.HIGHEST_PROTOCOL)\n",
    "    \n",
    "    pred_performance = {}\n",
    "    for d in d_test:\n",
    "        print(d)\n",
    "\n",
    "        predictor = DrugResponsePredictor(source_data=source_data[tissue][~np.isin(source_names[tissue], source_names_filtered[(ID, tissue)])],\\\n",
    "                                            method='consensus',\\\n",
    "                                            n_representations = 100,\\\n",
    "                                            target_data=X_target,\\\n",
    "                                            n_pv=d,\\\n",
    "                                            n_factors=n_factors,\\\n",
    "                                            n_jobs=n_jobs,\\\n",
    "                                            mean_center=mean_center,\\\n",
    "                                            std_unit=std_unit,\\\n",
    "                                            l1_ratio=l1_ratio)\n",
    "        predictor.alpha_values = list(np.logspace(-2,10,17))\n",
    "        predictor.verbose = 5\n",
    "        predictor.fit(X_source, y_source, use_data=True)\n",
    "        pred_performance[d] = predictor.compute_predictive_performance(X_source, y_source)\n",
    "        plt.plot(predictor.alpha_values, predictor.regression_model_.cv_results_['mean_test_score'], '+-')\n",
    "        plt.title(pred_performance[d])\n",
    "        plt.xscale('log')\n",
    "        plt.show()\n",
    "    \n",
    "    with open('./output/pred_performance/%s'%(pickle_file), 'wb') as f:\n",
    "        pickle.dump(pred_performance, f, pickle.HIGHEST_PROTOCOL)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### ElasticNet/Ridge comparison"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.model_selection import GroupKFold\n",
    "\n",
    "l1_ratio  = 0.\n",
    "\n",
    "pickle_file = 'elasticnet_drug_l1_ratio_%s_std.pkl'%(l1_ratio)\n",
    "if pickle_file in os.listdir('./output/pred_performance/'):\n",
    "    with open('./output/pred_performance/%s'%(pickle_file), 'rb') as f:\n",
    "        elasticnet_perf = pickle.load(f)\n",
    "\n",
    "for ID, tissue in zip(drug_IDs, tumor_tissues):\n",
    "    print(ID, tissue)\n",
    "    \n",
    "    pickle_file = 'en_std_drug_%s_tissue_%s_l1_ratio_%s_n_factors_%s.pkl'%(ID,\n",
    "                                                                          tissue,\n",
    "                                                                          l1_ratio,\n",
    "                                                                          n_factors)\n",
    "    \n",
    "    if pickle_file in os.listdir('./output/pred_performance/'):\n",
    "        print('%s, %s ALREADY COMPUTED'%(ID, tissue))\n",
    "        continue\n",
    "    \n",
    "    if (ID, tissue) in elasticnet_perf:\n",
    "        continue\n",
    "        \n",
    "    with open('./output/pred_performance/%s'%(pickle_file), 'wb') as f:\n",
    "        pickle.dump(dict(), f, pickle.HIGHEST_PROTOCOL)\n",
    "    \n",
    "    X_source = source_data_filtered[ID, tissue]\n",
    "    y_source = source_response_data[ID, tissue]\n",
    "    X_target = target_data[tissue]\n",
    "\n",
    "    #Parameters for the grid search\n",
    "    alpha_values = np.logspace(-5,10,16)\n",
    "    param_grid ={\n",
    "        'regression__alpha': alpha_values\n",
    "    }\n",
    "\n",
    "    #Grid search setup\n",
    "\n",
    "    k_fold_split = GroupKFold(10)\n",
    "    y_predicted = np.zeros(X_source.shape[0])\n",
    "    \n",
    "    for train_index, test_index in k_fold_split.split(X_source, y_source, y_source):\n",
    "        grid_en = GridSearchCV(Pipeline([\n",
    "                                ('normalization', StandardScaler(with_mean=mean_center, with_std=True)),\n",
    "                                ('regression', ElasticNet(l1_ratio) if l1_ratio > 0 else Ridge())\n",
    "                            ]),\\\n",
    "                            cv=10, n_jobs=30, param_grid=param_grid, verbose=1, scoring='neg_mean_squared_error')\n",
    "        grid_en.fit(X_source[train_index], y_source[train_index])\n",
    "        y_predicted[test_index] = grid_en.predict(X_source[test_index])\n",
    "    \n",
    "    #Fit grid search\n",
    "    grid_en.fit(X_source, y_source)\n",
    "    elasticnet_perf[ID, tissue] = scipy.stats.pearsonr(y_predicted, y_source)[0]\n",
    "    print(elasticnet_perf[ID, tissue])\n",
    "    \n",
    "    with open('./output/pred_performance/%s'%(pickle_file), 'wb') as f:\n",
    "        pickle.dump(elasticnet_perf[ID, tissue], f, pickle.HIGHEST_PROTOCOL)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Load pickle and look at results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "l1_ratio = 0\n",
    "l1_ratio_en = 0.\n",
    "\n",
    "two_pv_results = dict()\n",
    "consensus_pv_results = dict()\n",
    "source_pv_results = dict()\n",
    "target_pv_results = dict()\n",
    "en_results_std = dict()\n",
    "\n",
    "def sort_dictionary(d):\n",
    "    return {e:d[e] for e in sorted(d)}\n",
    "\n",
    "for ID, tissue in zip(drug_IDs, tumor_tissues):\n",
    "    print(ID, tissue)\n",
    "        \n",
    "    # Read results of consensus PVs\n",
    "    pickle_file = 'consensus_drug_%s_tissue_%s_l1_ratio_%s_n_factors_%s.pkl'%(ID,\n",
    "                                                                      tissue,\n",
    "                                                                      l1_ratio,\n",
    "                                                                      n_factors)\n",
    "    with open('./output/pred_performance/%s'%(pickle_file), 'rb') as f:\n",
    "        consensus_pv_results[ID,tissue] = sort_dictionary(pickle.load(f))\n",
    "        \n",
    "    \n",
    "    # Read results of EN\n",
    "    pickle_file = 'en_std_drug_%s_tissue_%s_l1_ratio_%s_n_factors_%s.pkl'%(ID,\n",
    "                                                                          tissue,\n",
    "                                                                          '0.0',\n",
    "                                                                          n_factors)\n",
    "    \n",
    "    with open('./output/pred_performance/%s'%(pickle_file), 'rb') as f:\n",
    "        en_results_std[ID,tissue] = pickle.load(f)\n",
    "        print(en_results[ID, tissue])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for ID, tissue in zip(drug_IDs, tumor_tissues):\n",
    "    # Plot for a specific number of PV\n",
    "    plt.plot([e[0] for e in consensus_pv_results[ID,tissue].items()],\n",
    "             [e[1] for e in consensus_pv_results[ID,tissue].items()],\n",
    "             label='consensus', linewidth=3, alpha=0.5, marker='+')\n",
    "    plt.plot([e[0] for e in source_pv_results[ID,tissue].items()],\n",
    "             [e[1] for e in source_pv_results[ID,tissue].items()],\n",
    "             label='source', linewidth=3, alpha=0.5, marker='+')\n",
    "    plt.plot([e[0] for e in target_pv_results[ID,tissue].items()],\n",
    "             [e[1] for e in target_pv_results[ID,tissue].items()],\n",
    "             label='target', linewidth=3, alpha=0.5, marker='+')\n",
    "    plt.plot([e[0] for e in two_pv_results[ID,tissue].items()],\n",
    "             [e[1] for e in two_pv_results[ID,tissue].items()],\n",
    "             label='2 pv', linewidth=3, alpha=0.5, marker='+')\n",
    "    plt.hlines(en_results[ID,tissue], xmin=0, xmax=plt.xlim()[1], label='Ridge', linewidth=3, alpha=0.7)\n",
    "\n",
    "    plt.title(drug_names[ID, tissue] + ' '+ tissue)\n",
    "    plt.xlabel('Number of Principal Vectors', fontsize=15)\n",
    "    plt.ylabel('Predictive Performance', fontsize=15)\n",
    "    plt.legend()    \n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "n_pv = 40\n",
    "\n",
    "perf_scatter = []\n",
    "for ID, tissue in zip(drug_IDs, tumor_tissues):\n",
    "    #print(ID, tissue)\n",
    "    if n_pv not in consensus_pv_results[ID,tissue]:\n",
    "        print(ID, tissue)\n",
    "        continue\n",
    "    plt.scatter(en_results_std[ID,tissue],\n",
    "                consensus_pv_results[ID,tissue][n_pv],\n",
    "               color='blue', marker='x', alpha=0.7)\n",
    "    perf_scatter.append([en_results_std[ID,tissue], consensus_pv_results[ID,tissue][n_pv]])\n",
    "    \n",
    "plt.xlabel('ElasticNet', fontsize=20)\n",
    "plt.ylabel('Consensus \\n representation', fontsize=20)\n",
    "plt.xticks(fontsize=15, color='black')\n",
    "plt.yticks(fontsize=15, color='black')\n",
    "plt.tight_layout()\n",
    "plt.xlim(0.1,0.8)\n",
    "plt.ylim(0.1,0.8)\n",
    "plt.plot(plt.xlim(), plt.xlim(), linewidth=3, alpha=0.5)\n",
    "#plt.savefig('./figures/fig4_pred_perf_consensus_%s_en_%s.png'%(l1_ratio, l1_ratio_en), dpi=300)\n",
    "plt.show()\n",
    "\n",
    "perf_scatter = np.array(perf_scatter)\n",
    "p = scipy.stats.pearsonr(perf_scatter[:,0], perf_scatter[:,1])\n",
    "print('Pearson Correlation: %s, %s'%(p[0], p[1]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.scatter(perf_scatter[:,1], (perf_scatter[:,0] - perf_scatter[:,1])/perf_scatter[:,0])\n",
    "\n",
    "np.median((perf_scatter[:,0] - perf_scatter[:,1])/perf_scatter[:,0])\n",
    "\n",
    "#for e in en_results:\n",
    "#    print(e, en_results[e], consensus_pv_results[e])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for ID, tissue in zip(drug_IDs, tumor_tissues):\n",
    "    #print(ID, tissue)\n",
    "    if n_pv not in consensus_pv_results[ID,tissue]:\n",
    "        print(ID, tissue)\n",
    "        continue\n",
    "    plt.scatter(en_results[ID,tissue],\n",
    "                en_results_std[ID,tissue],\n",
    "               color='blue', marker='x', alpha=0.7)\n",
    "    #perf_scatter.append([en_results[ID,tissue], consensus_pv_results[ID,tissue][n_pv]])"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python (precise_figures)",
   "language": "python",
   "name": "precise_figures"
  },
  "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.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}