plot_Fig_CCLE.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-05-12T15:23:52.979211Z",
     "start_time": "2022-05-12T15:23:43.342533Z"
    }
   },
   "outputs": [],
   "source": [
    "%load_ext rpy2.ipython\n",
    "import warnings\n",
    "warnings.filterwarnings('ignore')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Prep: EMT score calculation\n",
    "- Input files:\n",
    "    - `./CCLE/plot_heatmap/EMT_signature_genes_in_gencode.v31lift37.txt`: file contains 207 EMT signature genes\n",
    "    - `./CCLE/exp/exp_FPKM_EMT_signature_genes_in_gencode.v31lift37.txt`: file contains the expression value of 207 signature genes for all samples. Dimention genes x samples\n",
    "- Output matrix:\n",
    "    - output matrix would be a matrix with two columns (STATISTIC, pval)\n",
    "        - where STATISTIC represent the KS Statistics, -1 -- more epithelial; 1 -- more mysenchymal\n",
    "    - the matrix is also stored in `./CCLE/plot_heatmap/EMT_score_CCLE.txt`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-05-12T15:23:53.278690Z",
     "start_time": "2022-05-12T15:23:52.991615Z"
    }
   },
   "outputs": [],
   "source": [
    "%%R\n",
    "library(zeallot)\n",
    "KS_test <- function(x, y) {\n",
    "    n.x <- length(x)\n",
    "    n.y <- length(y)\n",
    "\n",
    "    w <- c(x, y)\n",
    "\n",
    "    z <- cumsum(ifelse(order(w) <= n.x, 1/n.x, -1/n.y))\n",
    "    z <- z[c(which(diff(sort(w)) != 0), n.x + n.y)] #exclude ties\n",
    "\n",
    "    alternative <- \"two.sided\"\n",
    "    STATISTIC <- switch(alternative, two.sided = z[which.max(abs(z))], \n",
    "                        greater = max(z), less = -min(z))\n",
    "\n",
    "    test = ks.test(x, y)\n",
    "    return(c(STATISTIC, test$p.value))\n",
    "}\n",
    "\n",
    "exp2EMT <- function(exp, signature) {\n",
    "    # exp is a matrix of expression value; colname is sample name\n",
    "    # signature is a vector of EMT signature corresponding to the rows of exp matrix\n",
    "    result = data.frame(STATISTIC = NULL, pval = NULL)\n",
    "    for (sample in colnames(exp)) {\n",
    "        exp_E <- exp[which(signature == 'Epi'), sample]\n",
    "        exp_M <- exp[which(signature == 'Mes'), sample]\n",
    "        c(STATISTIC, pval) %<-% KS_test(exp_E, exp_M)\n",
    "        result = rbind(result, data.frame(STATISTIC= STATISTIC, pval = pval))\n",
    "    }\n",
    "    rownames(result) = colnames(exp)\n",
    "    return(result)\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-05-12T15:23:54.994555Z",
     "start_time": "2022-05-12T15:23:53.282936Z"
    }
   },
   "outputs": [],
   "source": [
    "%%R\n",
    "exp <- read.table('./CCLE/exp/exp_FPKM_EMT_signature_genes_in_gencode.v31lift37.txt', header = T, sep = '\\t', quote = '', stringsAsFactors = FALSE, row.names = 1)\n",
    "EMT <- read.table('./CCLE/plot_heatmap/EMT_signature_genes_in_gencode.v31lift37.txt', header = T, sep = '\\t', quote = '', stringsAsFactors = FALSE, row.names = 7)\n",
    "result = exp2EMT(exp[, -c(1)], EMT[rownames(exp), 'EM_Signature'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-04-29T06:59:45.720076Z",
     "start_time": "2022-04-29T06:59:45.709855Z"
    }
   },
   "source": [
    "### Plot Heatmap\n",
    "- Input files\n",
    "    - EMT score matrix generated by running the code described earlier, or directly from file `./CCLE/plot_heatmap/EMT_score_CCLE.txt`\n",
    "    - `./CCLE/plot_heatmap/EMT_EMBO.txt`: EMT score retreived from reference paper (EMBO)\n",
    "    - `./CCLE/plot_heatmap/EMT_NBT.txt`: EMT classification retreived from reference paper (NBT)\n",
    "    - `./CCLE/plot_heatmap/example_subset_SE.MATS.JC_PSI.txt`: PSI value extracted from rMATS-turbo output; using a subset of events for illustration purpose"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### calculate correlation between PSI value and EMT score\n",
    "- Output:\n",
    "    - `./CCLE/plot_heatmap/CCLE_correlation_EMT-filteredPSI.txt`: correlation between EMT score and PSI value\n",
    "    - There is already one copy in the output folder that users can use directly\n",
    "- Note:\n",
    "    - This step will take several minutes to process if the input file contains lots of events (e.g. ~8min for 50,000 events)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-05-12T15:45:29.520514Z",
     "start_time": "2022-05-12T15:45:09.338065Z"
    }
   },
   "outputs": [],
   "source": [
    "%%R\n",
    "EMT <- read.table('./CCLE/plot_heatmap/EMT_score_CCLE.txt', header = T, sep = '\\t', stringsAsFactor = F, quote = \"\", check.name = F)\n",
    "PSI <- read.table('./CCLE/plot_heatmap/example_subset_SE.MATS.JC_PSI.txt', header = T, sep = '\\t', row.names = 1, quote = \"\", stringsAsFactor = F)\n",
    "\n",
    "result = list()\n",
    "for (ID in rownames(PSI)){\n",
    "    mydf = data.frame(EMT = EMT$STATISTIC, PSI = as.numeric(PSI[ID, rownames(EMT)]))\n",
    "    fit <- lm(EMT ~ PSI, data = mydf)\n",
    "    fit_summary = summary(fit)\n",
    "    result[[ID]] = setNames(c(fit_summary$r.squared, fit$coefficient[2], fit$coefficient[1]), \n",
    "                            c('r.squared', 'coefficient', 'intercept'))\n",
    "}\n",
    "myresult = do.call(rbind, result)\n",
    "\n",
    "write.table(myresult, './CCLE/plot_heatmap/CCLE_correlation_EMT-filteredPSI.txt', sep = '\\t', quote = FALSE, row.names = T, col.names = T)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Filter out events whose PSI value is highly correlated with EMT score (R2 > 0.4)\n",
    "- Output\n",
    "    - `r2_0.4_SE.MATS.JC_PSI.txt`: PSI values of events highly correlated with EMT scores"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-05-12T15:49:24.580281Z",
     "start_time": "2022-05-12T15:49:24.018069Z"
    }
   },
   "outputs": [],
   "source": [
    "%%bash\n",
    "cd ./CCLE/plot_heatmap\n",
    "awk -F '\\t' '($2 > 0.4) {print $0}' CCLE_correlation_EMT-filteredPSI.txt > r2_0.4_CCLE_correlation_EMT-filteredPSI.txt\n",
    "sed -i '1s/^/ID\\t/' r2_0.4_CCLE_correlation_EMT-filteredPSI.txt\n",
    "head -n 1 example_subset_SE.MATS.JC_PSI.txt > r2_0.4_SE.MATS.JC_PSI.txt\n",
    "awk -F '\\t' 'FNR==NR{a[$1]=$1;next} ($1 in a) {print $0}' r2_0.4_CCLE_correlation_EMT-filteredPSI.txt example_subset_SE.MATS.JC_PSI.txt >> r2_0.4_SE.MATS.JC_PSI.txt\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Plot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-05-12T15:52:50.328305Z",
     "start_time": "2022-05-12T15:52:43.330559Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "png \n",
       "  2 \n"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "%%R\n",
    "library(circlize)\n",
    "library(ComplexHeatmap)\n",
    "set.seed(1993)\n",
    "\n",
    "# EMT score and cell line annotation metadata\n",
    "EMT_score <- read.table('./CCLE/plot_heatmap/EMT_score_CCLE.txt', header = T, sep = '\\t', stringsAsFactor = F, quote = \"\", check.name = F)\n",
    "metadata <- read.table('./CCLE/metadata/metadata.txt', header = T, sep = '\\t', stringsAsFactor = FALSE, quote = \"\", check.name = F, row.names = 1)\n",
    "EMT <- cbind(EMT_score, metadata[rownames(EMT_score), ])\n",
    "EMT <- EMT[with(EMT, order(Tissue, STATISTIC)),]\n",
    "\n",
    "# from EMBO paper\n",
    "EMT_score_EMBO <- read.table('./CCLE/plot_heatmap/EMT_EMBO.txt', header = T, sep = '\\t', stringsAsFactor = F, quote = \"\", check.name = F, row.names = 1)\n",
    "EMT_score_EMBO = EMT_score_EMBO[rownames(EMT), ]\n",
    "EMT_EMBO = apply(EMT_score_EMBO, 1, \n",
    "      FUN = function(x){\n",
    "          score = as.numeric(x['Generic.EMT.CL.Ksscore'])\n",
    "          pval = as.numeric(x['Generic.EMT.CL.pv'])\n",
    "          if (is.na(score)) {\n",
    "              EM = NA\n",
    "          } else if (score < 0) {\n",
    "              EM = ifelse(pval < 0.05, 'Epi', 'intermediate Epi')\n",
    "          } else {\n",
    "              EM = ifelse(pval < 0.05, 'Mes', 'intermediate Mes')\n",
    "          }\n",
    "          return(EM)\n",
    "      })\n",
    "\n",
    "\n",
    "# from NBT paper\n",
    "EMT_NBT <- read.table('./CCLE/plot_heatmap/EMT_NBT.txt', header = T, sep = '\\t', stringsAsFactor = F, quote = \"\", check.name = F, row.names = 1)\n",
    "EMT_NBT <- EMT_NBT[rownames(EMT), ]\n",
    "EMT_NBT[which(EMT_NBT$EMT_status == 'Epithelial'), 'EMT_status'] = 'Epi'\n",
    "EMT_NBT[which(EMT_NBT$EMT_status == 'Mesenchymal'), 'EMT_status'] = 'Mes'\n",
    "EMT_NBT[which(EMT_NBT$EMT_status == 'OtherEMT'), 'EMT_status'] = 'Other'\n",
    "\n",
    "################ ht_CCLE ################\n",
    "mydf <- read.table('./CCLE/plot_heatmap/r2_0.4_SE.MATS.JC_PSI.txt', header = T, sep = '\\t', row.names = 1, quote = \"\", stringsAsFactor = F)\n",
    "mydf_CCLE <- mydf[, rownames(EMT)]\n",
    "mat <- t(apply(mydf_CCLE, 1, scale))\n",
    "\n",
    "\n",
    "column_ha = HeatmapAnnotation(#'EMT score' = anno1,\n",
    "                              'EMT score' = anno_points(EMT$STATISTIC, \n",
    "                                                        ylim = c(-1, 1), \n",
    "                                                        size = unit(0.8, \"mm\"),\n",
    "                                                        gp = gpar(col = ifelse(EMT$STATISTIC > 0 & EMT$pval < 0.05, \"dodgerblue4\", ifelse(EMT$STATISTIC < -0 & EMT$pval < 0.05, \"orangered4\", \"grey\"))),\n",
    "                                                        axis_param = list(side = \"left\", at = c(-1, 0, 1), labels = c(\"Epi -1\", \"0\", \"Mes  1\"))\n",
    "                                                         ), \n",
    "                              'Tan et al., 2014' = factor(EMT_EMBO, levels = c('Epi', 'intermediate Epi', 'intermediate Mes', 'Mes', 'NA')), \n",
    "                              'Klijn et al., 2015' = factor(EMT_NBT$EMT_status, levels = c('Epi', 'Other', 'Mes', 'NA')),\n",
    "                              'disease stage' = factor(EMT$disease_stage, levels = c('benign_neoplasia', 'primary', 'metastasis')),\n",
    "                              'tissue' = EMT$Tissue,\n",
    "                              col = list('disease stage' = c(\"benign_neoplasia\" = \"darkgrey\", \"primary\" = \"turquoise\", \"metastasis\" = \"deeppink\"),\n",
    "                                         'Tan et al., 2014' = c('Epi' = 'orangered4', 'intermediate Epi' = 'orangered', 'intermediate Mes' = 'dodgerblue', 'Mes' = 'dodgerblue4', 'NA' = 'white'), \n",
    "                                         'Klijn et al., 2015' = c('Epi' = 'orangered4', 'Other' = 'grey', 'Mes' = 'dodgerblue4', 'NA' = 'white')), \n",
    "                              na_col = \"white\")\n",
    "lgd_EMT = Legend(labels = c('Epi', 'intermediate', 'Mes'), legend_gp = gpar(fill = c('orangered4', 'grey', 'dodgerblue4')), title = \"EMT score\")\n",
    "\n",
    "\n",
    "ht_CCLE  = Heatmap(mat, \n",
    "            show_row_names = FALSE,\n",
    "            na_col = \"gray30\",\n",
    "            column_split = EMT$Tissue,\n",
    "            top_annotation = column_ha,\n",
    "            # right_annotation = row_ha,\n",
    "            cluster_columns = F, \n",
    "            heatmap_legend_param = list(\n",
    "                title = \"CCLE PSI Z-Score\",\n",
    "                legend_height = unit(4, \"cm\"),\n",
    "                #direction = \"horizontal\",\n",
    "                title_position = \"leftcenter-rot\"),\n",
    "            column_title = 'Cancer Cell Line Encyclopedia (CCLE)'\n",
    "        )\n",
    "\n",
    "################ PLOT ################\n",
    "pdf(\"./CCLE/plot_heatmap/plot_heatmap_sortedByEMT_usingCorrelation.pdf\", width = 15, height = 6)\n",
    "lgd_list = list('EMT score' = lgd_EMT)\n",
    "\n",
    "ht = draw(\n",
    "    ht_CCLE,\n",
    "    annotation_legend_list = lgd_list,\n",
    "    annotation_legend_side = \"right\",\n",
    "    heatmap_legend_side = \"right\"\n",
    ")\n",
    "\n",
    "dev.off()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.15"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {},
   "toc_section_display": true,
   "toc_window_display": true
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}