eCLIP_cross_platform_benchmark.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Cross-platform benchmark: BMF/GraphProt/iDeepE/DeepCLIP"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os \n",
    "import numpy as np\n",
    "from matplotlib import pyplot as plt\n",
    "import pandas as pd\n",
    "from sklearn.metrics import roc_curve, roc_auc_score, average_precision_score, precision_recall_curve, auc\n",
    "import pickle\n",
    "import seaborn as sns\n",
    "from scipy.stats import wilcoxon"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### helper functions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "def read_predictions(pred_file):\n",
    "    pred = np.loadtxt(pred_file)\n",
    "    return pred"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "def read_predictions_graphprot(pred_file):\n",
    "    pred = pd.read_csv(pred_file, sep='\\t', header=None)\n",
    "    return pred.iloc[:,2].values"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "def read_sequence_ids_clip(file_path):\n",
    "    bed_table = pd.read_csv(file_path, sep='\\t', header=None)\n",
    "    return bed_table.iloc[:,4].values"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "def compile_array_by_id(scores, ids):\n",
    "    '''\n",
    "    gets the scores and ids for positive sequences and combines it to \n",
    "    a single score array where each clip peak gets one score\n",
    "    '''\n",
    "    #make a dataframe and split it based on ids\n",
    "    df = pd.DataFrame({'scores':scores, 'ids':ids})    \n",
    "    final_scores = []\n",
    "    for identifier, df in df.groupby('ids'):\n",
    "        #keep the max value in each split\n",
    "        final_score = df.loc[:,'scores'].mean()\n",
    "        final_scores.append(final_score)\n",
    "    return np.array(final_scores)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_new_scores_all(all_scores, clip_factor):\n",
    "    bg_to_ps_ratio=1\n",
    "    clip_all_file = os.path.join('../../cbscratch/eclip_encode/bed_splitted/',f'{clip_factor}_split_all.bed')\n",
    "    ids = read_sequence_ids_clip(clip_all_file)\n",
    "\n",
    "    \n",
    "    final_scores = compile_array_by_id(all_scores, ids)\n",
    "    \n",
    "    no_seqs_pos = int(len(final_scores)/(bg_to_ps_ratio+1))    \n",
    "    y_true = np.append(np.ones(no_seqs_pos), np.zeros(no_seqs_pos*bg_to_ps_ratio))\n",
    "\n",
    "    return y_true, final_scores"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "def create_id_list(ids, bg_to_ps_ratio=10):\n",
    "    id_list = []\n",
    "    id_list.append(ids)\n",
    "    new_ids = ids\n",
    "    for i in range(bg_to_ps_ratio):\n",
    "        new_ids = ids + new_ids[-1] + 1\n",
    "        id_list.append(new_ids)\n",
    "        \n",
    "    return np.concatenate(id_list)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_roc(fpr, tpr, labels, roc_auc, factor_name):\n",
    "    lw = 2\n",
    "    for i in range(len(fpr)):\n",
    "        plt.plot(fpr[i], tpr[i], color=colors[i], lw=lw, label=f'{labels[i]} ROC curve (area = {roc_auc[i]:.2f})')\n",
    "    plt.plot([0, 1], [0, 1], color='navy', lw=lw, linestyle='--')\n",
    "    plt.xlim([0.0, 1.0])\n",
    "    plt.ylim([0.0, 1.05])\n",
    "    plt.xlabel('False Positive Rate')\n",
    "    plt.ylabel('True Positive Rate')\n",
    "    plt.title(f'{factor_name}\\nReceiver operating characteristic example')\n",
    "    plt.legend(loc=\"lower right\")\n",
    "    plt.savefig(os.path.join(roc_plot_dir,f'{factor_name}.pdf'))\n",
    "    plt.close()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Train on HT-SELEX test on eCLIP"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "#directories where prediction files for each method is stored\n",
    "graphprot_path = 'graphprot0/selex_eclip/'\n",
    "ideepe_path = 'ideepe0/selex_eclip/'\n",
    "deepclip_path = 'deepclip0/selex_eclip/'\n",
    "bipartite_cs3_path = 'scripts/param/selex0/benchmark_selex_eclip_cs3/'\n",
    "bipartite_cs4_path = 'scripts/param/selex0/benchmark_selex_eclip_cs4/'\n",
    "bipartite_cs5_path = 'scripts/param/selex0/benchmark_selex_eclip_cs5/'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "all_combos_table = pd.read_csv('../../cbscratch/mapping_selex_eclip_allcombinations.txt', sep='\\t', header=None)\n",
    "selex_factors = all_combos_table.iloc[:,0].values\n",
    "eclip_factors = all_combos_table.iloc[:,1].values\n",
    "combos = [f'{a}_{b}' for a,b in zip(selex_factors, eclip_factors)]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "datasets = ['BMF cs3','BMF cs4','BMF cs5','GraphProt','iDeepE', 'DeepCLIP']\n",
    "colors = ['#6baed6','#3182bd','#08519c','#fb9a99','#b2df8a','#fdbf6f']\n",
    "no_experiments = len(datasets)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "deepclip predictions for KHDRBS1-rep1_KHDRBS1_K562 not found (runtime error)\n",
      "deepclip predictions for ZC3H8-rep0_ZC3H8_K562 not found (runtime error)\n",
      "deepclip predictions for ZC3H8-rep1_ZC3H8_K562 not found (runtime error)\n"
     ]
    }
   ],
   "source": [
    "auc_list = []\n",
    "ap_list = []\n",
    "roc_plot_dir = 'plots/eclip/roc_plots/'\n",
    "\n",
    "for i, combo in enumerate(combos):\n",
    "    clip_factor = eclip_factors[i]\n",
    "    \n",
    "    fpr = []\n",
    "    tpr = []\n",
    "    labels = []\n",
    "    factor_auc = []\n",
    "\n",
    "    bipartite_cs3_file = os.path.join(bipartite_cs3_path, f'{combo}.predictions')\n",
    "    bipartite_cs4_file = os.path.join(bipartite_cs4_path, f'{combo}.predictions')\n",
    "    bipartite_cs5_file = os.path.join(bipartite_cs5_path, f'{combo}.predictions')\n",
    "    graphprot_file = os.path.join(graphprot_path, f'{combo}_eclip.predictions')\n",
    "    ideepe_file = os.path.join(ideepe_path, f'{combo}_predictions.txt')\n",
    "    deepclip_file = os.path.join(deepclip_path, f'{combo}_predictions.tsv')\n",
    "\n",
    "    bipartite_cs3_exists = os.path.isfile(bipartite_cs3_file) \n",
    "    bipartite_cs4_exists = os.path.isfile(bipartite_cs4_file) \n",
    "    bipartite_cs5_exists = os.path.isfile(bipartite_cs5_file) \n",
    "    graphprot_exists = os.path.isfile(graphprot_file) \n",
    "    ideepe_exists = os.path.isfile(ideepe_file) \n",
    "    deepclip_exists = os.path.isfile(deepclip_file) \n",
    "    \n",
    "\n",
    "    if bipartite_cs3_exists:\n",
    "        y_bipartite3 = read_predictions(bipartite_cs3_file)\n",
    "        y_true, y_bipartite3 = get_new_scores_all(y_bipartite3, clip_factor)\n",
    "        \n",
    "        auc_bipartite3 = roc_auc_score(y_true, y_bipartite3)\n",
    "        \n",
    "        #plotting\n",
    "        fpr0, tpr0, _ = roc_curve(y_true, y_bipartite3)\n",
    "        fpr.append(fpr0)\n",
    "        tpr.append(tpr0)\n",
    "        labels.append('bipartite_cs3')\n",
    "        factor_auc.append(auc_bipartite3)\n",
    "        \n",
    "        #print(y_bipartite.shape, fpr0.shape)\n",
    "        \n",
    "        ap_bipartite3 = average_precision_score(y_true, y_bipartite3)\n",
    "        \n",
    "        precision, recall, thresholds = precision_recall_curve(y_true, y_bipartite3)\n",
    "    else:\n",
    "        print(f'bipartite predictions for {combo} not found (runtime error)')\n",
    "        auc_bipartite3 = np.nan\n",
    "        ap_bipartite3 = np.nan\n",
    "    \n",
    "    if bipartite_cs4_exists:\n",
    "        y_bipartite4 = read_predictions(bipartite_cs4_file)\n",
    "        y_true, y_bipartite4 = get_new_scores_all(y_bipartite4, clip_factor)\n",
    "        \n",
    "        auc_bipartite4 = roc_auc_score(y_true, y_bipartite4)\n",
    "        \n",
    "        #plotting\n",
    "        fpr0, tpr0, _ = roc_curve(y_true, y_bipartite4)        \n",
    "        fpr.append(fpr0)\n",
    "        tpr.append(tpr0)\n",
    "        labels.append('bipartite_cs4')\n",
    "        factor_auc.append(auc_bipartite4)\n",
    "        \n",
    "        ap_bipartite4 = average_precision_score(y_true, y_bipartite4)\n",
    "        precision, recall, thresholds = precision_recall_curve(y_true, y_bipartite4)\n",
    "    else:\n",
    "        print(f'bipartite cs4 predictions for {combo} not found (runtime error)')\n",
    "        auc_bipartite4 = np.nan\n",
    "        ap_bipartite4 = np.nan\n",
    "    \n",
    "    if bipartite_cs5_exists:\n",
    "        y_bipartite5 = read_predictions(bipartite_cs5_file)\n",
    "        y_true, y_bipartite5 = get_new_scores_all(y_bipartite5, clip_factor)\n",
    "        \n",
    "        auc_bipartite5 = roc_auc_score(y_true, y_bipartite5)\n",
    "        \n",
    "        #plotting\n",
    "        fpr0, tpr0, _ = roc_curve(y_true, y_bipartite5)        \n",
    "        fpr.append(fpr0)\n",
    "        tpr.append(tpr0)\n",
    "        labels.append('bipartite_cs5')\n",
    "        factor_auc.append(auc_bipartite5)\n",
    "        \n",
    "        ap_bipartite5 = average_precision_score(y_true, y_bipartite5)\n",
    "        precision, recall, thresholds = precision_recall_curve(y_true, y_bipartite5)\n",
    "    else:\n",
    "        print(f'bipartite cs5 predictions for {combo} not found (runtime error)')\n",
    "        auc_bipartite5 = np.nan\n",
    "        ap_bipartite5 = np.nan\n",
    "                \n",
    "    if graphprot_exists:\n",
    "        y_graphprot = read_predictions_graphprot(graphprot_file)\n",
    "        y_true, y_graphprot = get_new_scores_all(y_graphprot, clip_factor)\n",
    "        \n",
    "        #plotting\n",
    "        auc_graphprot = roc_auc_score(y_true, y_graphprot)\n",
    "        fpr0, tpr0, _ = roc_curve(y_true, y_graphprot)\n",
    "        fpr.append(fpr0)\n",
    "        tpr.append(tpr0)\n",
    "        labels.append('graphprot')\n",
    "        factor_auc.append(auc_graphprot)\n",
    "        \n",
    "        ap_graphprot = average_precision_score(y_true, y_graphprot)\n",
    "        precision, recall, thresholds = precision_recall_curve(y_true, y_graphprot)\n",
    "    else:\n",
    "        print(f'graphprot predictions for {combo} not found (runtime error)')\n",
    "        auc_graphprot = np.nan\n",
    "        ap_graphprot = np.nan\n",
    "\n",
    "\n",
    "    if ideepe_exists:\n",
    "        y_ideepe = read_predictions(ideepe_file)\n",
    "        y_true, y_ideepe = get_new_scores_all(y_ideepe, clip_factor)\n",
    "        \n",
    "        #plotting\n",
    "        auc_ideepe = roc_auc_score(y_true, y_ideepe)\n",
    "        fpr0, tpr0, _ = roc_curve(y_true, y_ideepe)\n",
    "        fpr.append(fpr0)\n",
    "        tpr.append(tpr0)\n",
    "        labels.append('iDeepE')\n",
    "        factor_auc.append(auc_ideepe)\n",
    "        \n",
    "        #print(y_ideepe.shape, fpr0.shape)\n",
    "        ap_ideepe = average_precision_score(y_true, y_ideepe)\n",
    "        precision, recall, thresholds = precision_recall_curve(y_true, y_ideepe)\n",
    "\n",
    "    else:\n",
    "        print(f'ideepe predictions for {combo} not found (runtime error)')\n",
    "        auc_ideepe = np.nan\n",
    "        ap_ideepe = np.nan\n",
    "\n",
    "        \n",
    "    if deepclip_exists:\n",
    "        y_deepclip = read_predictions_graphprot(deepclip_file)\n",
    "        y_true, y_deepclip = get_new_scores_all(y_deepclip, clip_factor)\n",
    "        \n",
    "        #plotting\n",
    "        auc_deepclip = roc_auc_score(y_true, y_deepclip)\n",
    "        fpr0, tpr0, _ = roc_curve(y_true, y_deepclip)\n",
    "        fpr.append(fpr0)\n",
    "        tpr.append(tpr0)\n",
    "        labels.append('DeepCLIP')\n",
    "        factor_auc.append(auc_deepclip)\n",
    "        \n",
    "        ap_deepclip = average_precision_score(y_true, y_deepclip)\n",
    "        precision, recall, thresholds = precision_recall_curve(y_true, y_deepclip)\n",
    "    else:\n",
    "        print(f'deepclip predictions for {combo} not found (runtime error)')\n",
    "        auc_deepclip = np.nan\n",
    "        ap_deepclip = np.nan\n",
    "    \n",
    "    #print(combo)\n",
    "    plot_roc(fpr, tpr, labels, factor_auc, combo)\n",
    "    \n",
    "    auc_list.append([auc_bipartite3, auc_bipartite4, auc_bipartite5, auc_graphprot, auc_ideepe, auc_deepclip])  \n",
    "    ap_list.append([ap_bipartite3, ap_bipartite4, ap_bipartite5, ap_graphprot, ap_ideepe, ap_deepclip])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "plotting_dir = 'plots/eclip'\n",
    "pickle.dump([auc_list, ap_list], open(os.path.join(plotting_dir,'metrics.pkl'), 'wb' ))\n",
    "\n",
    "#auc_list, ap_list = pickle.load(open(os.path.join(plotting_dir,'metrics.pkl'), 'rb' ))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Split by factor"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "def combine_reps_plot(metric_list, metric_name):\n",
    "    metric_df = pd.DataFrame(metric_list, columns=datasets)\n",
    "    metric_df['factors'] = [c.split('-')[0] for c in combos]\n",
    "\n",
    "    metric_by_factor = {}\n",
    "    metric_by_factor_range = []\n",
    "    for label, df in metric_df.groupby('factors'):\n",
    "        values = df.replace(0, np.NaN).mean()\n",
    "        metric_by_factor_range.append((df.replace(0, np.NaN).min(),df.replace(0, np.NaN).max()))\n",
    "        metric_by_factor[label] = values\n",
    "\n",
    "    metric_by_factor = pd.DataFrame(metric_by_factor).T.loc[:,datasets].values\n",
    "\n",
    "    metric_by_factor_max = [np.argmax(arr) for arr in np.array(metric_by_factor)]\n",
    "    _ = plt.hist(metric_by_factor_max, bins=np.arange(-0.5,3.5,1), rwidth=0.8)\n",
    "    _ = plt.xticks(np.arange(0,3,1), datasets, rotation=90)\n",
    "    plt.title(f'best at {metric_name} (by factor)')\n",
    "    plt.savefig(os.path.join(plotting_dir,f'{metric_name}_by_factor_summary.pdf'), bbox_inches='tight')\n",
    "    plt.close()\n",
    "    \n",
    "    fig, ax = plt.subplots(figsize=(4,10))\n",
    "    \n",
    "    for i in range(len(datasets)):\n",
    "        y_pos = np.arange(len(metric_by_factor)) - 0.3 + 0.17*i\n",
    "        auc_i = np.array(metric_by_factor)[:,i]\n",
    "        ax.barh(y_pos, auc_i, align='center', height=0.15, label=datasets[i], color=colors[i], edgecolor='black')\n",
    "\n",
    "    best_colors = [colors[i] for i in metric_by_factor_max]\n",
    "    ax.scatter([0.95]*len(metric_by_factor), np.arange(len(metric_by_factor)), color=best_colors, s=50)\n",
    "    ax.vlines(x=0.5, ymin=-1, ymax=len(metric_by_factor)+1, ls='--')\n",
    "\n",
    "    ax.set_yticks(y_pos-0.3)\n",
    "    ax.set_yticklabels(metric_df.loc[:,'factors'].unique())\n",
    "    ax.invert_yaxis()  # labels read top-to-bottom\n",
    "\n",
    "    ax.set_ylim(len(metric_by_factor)+1, -1)\n",
    "\n",
    "    ax.legend(loc=(1.1,0.9))\n",
    "    ax.set_title(f'{metric_name} (by factor)')\n",
    "    \n",
    "    pickle.dump([metric_df.loc[:,'factors'].unique(), metric_by_factor, metric_by_factor_range], open(os.path.join(plotting_dir,f'{metric_name}_byfactor.pkl'), 'wb' ))\n",
    "\n",
    "    plt.savefig(os.path.join(plotting_dir,f'{metric_name}_byfactor_selex_eclip_all.pdf'), bbox_inches='tight')\n",
    "    plt.close()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "metric_list = ap_list\n",
    "metric_name = 'AP'\n",
    "combine_reps_plot(metric_list, metric_name)\n",
    "\n",
    "metric_list = auc_list\n",
    "metric_name = 'AUROC'\n",
    "combine_reps_plot(metric_list, metric_name)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python [conda env:rbp_motif]",
   "language": "python",
   "name": "conda-env-rbp_motif-py"
  },
  "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.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}