20191008_DDX56_Expression_Renalysis.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "# Re-evaluation of DDX56 expression in GNS lines\n",
    "\n",
    "L.Richards  \n",
    "October 2019  \n",
    "  \n",
    "---"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "library(ggplot2)\n",
    "library(\"Seurat\", \n",
    "        lib =\"/mnt/work1/users/pughlab/projects/BTSCs_scRNAseq/Manuscript/lib/\" ) #v2.3.4\n",
    "library(reshape2)\n",
    "library(vioplot)\n",
    "\n",
    "setwd(\"/mnt/work1/users/pughlab/projects/DDX56_scRNA_PearsonLab\")"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "### load in processed data for G523_L\n",
    "load(\"~/pughlab/projects/BTSCs_scRNAseq/Manuscript_G607removed/CellCyle_Regression/scClustViz/opt_solution/G523_L_res.0.2.RData\")\n",
    "\n",
    "### set ID as DDX56+ or DDX56 -\n",
    "\n",
    "\n",
    "BTSC@meta.data$DDX56_class <- ifelse(BTSC@data[grep(\"DDX56\", rownames(BTSC@data)), ] > 0,\n",
    "                                     \"DDX56+\",\n",
    "                                     \"DDX56-\"\n",
    "                                    \n",
    "                                    )\n",
    "table(BTSC@meta.data$DDX56_class)\n",
    "BTSC@ident <- factor(BTSC@meta.data$DDX56_class)\n",
    "BTSC@ident[1:5]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "### 1.0 Plot additional genes for co-expression with DDX56\n",
    "---\n",
    "\n",
    "https://docs.google.com/spreadsheets/d/1UxN5XvyCmyVXOSmSkb6WU859ClF1q0d5vYOYD4d9jtY/edit#gid=0\n",
    "\n",
    "These genes **are not** in the data:  \n",
    " [1] \"TNFSF10\" \"P14ARF\"  \"IFN\"     \"IRF4\"    \"CHAT\"    \"H1FOO\"   \"IIAE7\"  \n",
    " [8] \"ISG12A\"  \"ISG12\"   \"TRAIL\"  \n",
    " \n",
    " These genes **are** in the data:  \n",
    "  [1] \"IRF3\"   \"IFI27\"  \"IFIT1\"  \"RPL18A\" \"TP53\"   \"RPL18\"  \"AK1\"    \"RPL10A\"  \n",
    " [9] \"RPL31\"  \"RPL36A\" \"RPL37\"  \"RPL39\"  \"RPL5\"   \"SRSF10\" \"RPL37A\" \"SRSF10\"  \n",
    "[17] \"CBFB\"   \"TRMT6\"  \"PES1\"   \"WDHD1\"  \"LYAR\"   \"CBX3\"  \n",
    "\n",
    "I would be weary of the RPL genes, their expression tends to modulate due to technical reasons."
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "### load in gene plot\n",
    "\n",
    "load(\"GeneList_Sept2019.Rdata\")\n",
    "genes <- genes[genes %in% rownames(BTSC@data)]\n",
    "genes\n",
    "\n",
    "\n",
    "\n",
    "### co-expression plots\n",
    "\n",
    "pdf(\"G523_DDX56_coexpression_Sept2019.pdf\")\n",
    "\n",
    "for(i in 1:length(genes)){\n",
    "    \n",
    "#coexpression plor\n",
    "\n",
    "a <- plot(BTSC@data[\"DDX56\", ],\n",
    "         BTSC@data[genes[i], ],\n",
    "         xlab = \"Normalized DDX56 Expression\",\n",
    "         ylab = paste0(\"Normalized \", genes[i], \" Expression\") ,\n",
    "         main = \"G523_L\"\n",
    "        )\n",
    "print(a)\n",
    "    \n",
    "#compare expression between DDX56+/- cells\n",
    "\n",
    "x1 <- BTSC@data[genes[i], BTSC@ident == \"DDX56+\"]\n",
    "x2 <- BTSC@data[genes[i], BTSC@ident == \"DDX56-\"]\n",
    "\n",
    "vioplot(x1, \n",
    "        x2,\n",
    "        names=c(\"DDX56+\", \"DDX56-\"), \n",
    "        col=c(\"white\", \"grey\"),\n",
    "        ylab = \"Normalized Expression\"\n",
    "       )\n",
    "title(paste0(\"G523_L: \", genes[i]))\n",
    "    \n",
    "}\n",
    "\n",
    "#plot all genes as heatmap\n",
    "\n",
    "DoHeatmap(BTSC,\n",
    "          genes.use = genes,\n",
    "          group.by = \"DDX56_class\",\n",
    "          remove.key = FALSE,\n",
    "          cex.col = 0\n",
    "         )\n",
    "\n",
    "\n",
    "dev.off()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "### 2.0 Subtype cells in G523 for TCGA subtypes\n",
    "---\n",
    "\n"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "###load in function to score cells and generate random signatures\n",
    "### scoreAndClassify_AddModuleScore()\n",
    "source(\"randomScore.R\")\n",
    "\n",
    "\n",
    "###load in gene signatures from TCGA verhaak\n",
    "load(\"TCGA_GeneSignatures.Rdata\")\n",
    "str(subtypes)\n",
    "\n",
    "### score cells for each of the four gene signatures\n",
    "\n",
    "      scored <- AddModuleScore(BTSC,\n",
    "                       genes.list = subtypes,\n",
    "                       n.bin = 25,\n",
    "                       seed.use = 1,\n",
    "                       ctrl.size = 100,\n",
    "                       random.seed = 123,\n",
    "                       enrich.name = names(subtypes)     \n",
    "                      )\n",
    "\n",
    "## generate random lists of genes for each of the signatures\n",
    "\n",
    "size <- as.numeric(lapply(subtypes, length)) #extract length of gene signatures\n",
    "\n",
    "random.sigs <- list()\n",
    "\n",
    "for (i in 1:length(size)){\n",
    "    \n",
    "    print(size[i])\n",
    "    print(subtypes[i])\n",
    "\n",
    "    sigs <- replicate(100, lapply(size[i], function(x) sample(rownames(scored@data), size = x)))\n",
    "    names(sigs) <- paste0(\"Random\", 1:100, \"_\", names(subtypes)[i]) \n",
    "        \n",
    "    random.sigs <- do.call(c, list(random.sigs, sigs))\n",
    "                            \n",
    "}\n",
    "\n",
    "save(random.sigs, file = \"RandomGeneSignatures.Rdata\")\n",
    "                                  \n",
    "## score 100 random gene signatures\n",
    "                                  \n",
    "scored <- AddModuleScore(scored,\n",
    "                         genes.list = random.sigs,\n",
    "                       n.bin = 25,\n",
    "                       seed.use = 1,\n",
    "                       ctrl.size = 100,\n",
    "                       random.seed = 123,\n",
    "                       enrich.name = names(random.sigs)     \n",
    "                      )\n",
    "                                  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "##isolate random scores\n",
    "                                  \n",
    "random.scores <- scored@meta.data[grep(\"Random\", colnames(scored@meta.data))]\n",
    "scored@meta.data <- scored@meta.data[,-grep(\"Random\", colnames(scored@meta.data))]\n",
    "                                  \n",
    "\n",
    "### Define top 5% cutoff for classification and plot\n",
    "\n",
    "        \n",
    "        cutoffs <- c()\n",
    "            \n",
    "        for (sig.name in names(subtypes)){\n",
    "         \n",
    "            print(sig.name)\n",
    "            range <- grep(sig.name, colnames(random.scores))\n",
    "            #print(range)\n",
    "            \n",
    "            bb <- as.vector(as.matrix(random.scores[,range]))\n",
    "            print(length(bb))\n",
    "            \n",
    "            d <- density(bb)\n",
    "            cutoff <- qnorm(0.95, mean = mean(bb), sd = sd(bb))\n",
    "            print(cutoff)\n",
    "            \n",
    "            filename <- paste0(sig.name, \"AddModuleScore_RandomListCutoff.pdf\")\n",
    "            pdf(filename)\n",
    "\n",
    "            plot_raw <- plot(d,\n",
    "                 main = paste0(sig.name,\" - Random Signatures\"),\n",
    "                 xlab = \"AddModuleScore()\"\n",
    "                )\n",
    "            legend(\"topright\",\n",
    "            legend = c(paste(\"Mean:\", round(mean(bb),3)),\n",
    "                          paste(\"Sd:\", round(sd(bb), 3)),\n",
    "                          paste(\"Cutoff:\", round(cutoff,3))\n",
    "                     ),\n",
    "               bty = 'n',\n",
    "               cex = 1\n",
    "          )\n",
    "        \n",
    "        abline(v=cutoff, lty = 2, col = \"red\")\n",
    "        dev.off()\n",
    "        cutoffs <- c(cutoffs, cutoff)\n",
    "        \n",
    "        }\n",
    "        \n",
    "        names(cutoffs) <- names(subtypes)\n",
    "        cutoffs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "#### classify the cells\n",
    "\n",
    "meta <- scored@meta.data\n",
    "            \n",
    "for (i in 1:length(subtypes)){\n",
    "    \n",
    "    subset <- meta[ ,grep(names(subtypes[i]), colnames(meta))]\n",
    "    \n",
    "    col.name <- paste0(\"Classify_\", names(subtypes[i]), \"_\", \n",
    "                           round(as.numeric(cutoffs[grep(names(subtypes[i]), names(cutoffs))]), 4)\n",
    "                           )\n",
    "    \n",
    "    meta[col.name] <- subset > as.numeric(cutoffs[grep(names(subtypes[i]), names(cutoffs))])\n",
    "        \n",
    "    }\n",
    "\n",
    "\n",
    "### make a single classification co\n",
    "\n",
    "meta$Class <- do.call(paste, c(meta[17:20], sep=\"-\"))\n",
    "            \n",
    "meta$Class <- gsub(\"FALSE-FALSE-FALSE-FALSE\", \"Unclassified\", meta$Class)\n",
    "meta$Class <- gsub(\"FALSE-FALSE-FALSE-TRUE\", \"CLA\", meta$Class)\n",
    "meta$Class <- gsub(\"FALSE-FALSE-TRUE-FALSE\", \"MES\", meta$Class)\n",
    "meta$Class <- gsub(\"FALSE-FALSE-TRUE-TRUE\", \"MES+CLA\", meta$Class)\n",
    "meta$Class <- gsub(\"FALSE-TRUE-FALSE-FALSE\", \"PRO\", meta$Class)\n",
    "meta$Class <- gsub(\"FALSE-TRUE-FALSE-TRUE\", \"PRO+CLA\", meta$Class)\n",
    "\n",
    "meta$Class <- gsub(\"FALSE-TRUE-TRUE-FALSE\", \"PRO+MES\", meta$Class)\n",
    "meta$Class <- gsub(\"FALSE-TRUE-TRUE-TRUE\", \"PRO+MES+CLA\", meta$Class)\n",
    "meta$Class <- gsub(\"TRUE-FALSE-FALSE-FALSE\", \"NEU\", meta$Class)\n",
    "meta$Class <- gsub(\"TRUE-FALSE-FALSE-TRUE\", \"NEU+CLA\", meta$Class)\n",
    "meta$Class <- gsub(\"TRUE-FALSE-TRUE-FALSE\", \"NEU+MES\", meta$Class)\n",
    "meta$Class <- gsub(\"TRUE-FALSE-TRUE-TRUE\", \"NEU+MES+CLA\", meta$Class)\n",
    "\n",
    "meta$Class <- gsub(\"TRUE-TRUE-FALSE-FALSE\", \"NEU+PRO\", meta$Class)\n",
    "meta$Class <- gsub(\"TRUE-TRUE-FALSE-TRUE\", \"NEU+PRO+CLA\", meta$Class)\n",
    "meta$Class <- gsub(\"TRUE-TRUE-TRUE-FALSE\", \"NEU+PRO+MES\", meta$Class)\n",
    "meta$Class <- gsub(\"TRUE-TRUE-TRUE-TRUE\", \"NEU+PRO+MES+CLA\", meta$Class)\n",
    "\n",
    "meta <- meta[ ,c(\"DDX56_class\", \n",
    "                 \"Neural1\", \n",
    "                 \"Proneural2\",\n",
    "                 \"Mesenchymal3\",\n",
    "                 \"Classical4\",\n",
    "                 \"Classify_Neural_0.05\",\n",
    "                 \"Classify_Proneural_0.047\",\n",
    "                 \"Classify_Mesenchymal_0.0448\",\n",
    "                 \"Classify_Classical_0.0446\",\n",
    "                 \"Class\")]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "### save the data\n",
    "\n",
    "write.csv(meta, file = paste0(\"DDX56_ScoredTCGA\", \".csv\"), row.names = T)\n",
    "            \n",
    "        save(random.sigs,\n",
    "             subtypes,\n",
    "             scored,\n",
    "             random.scores,\n",
    "             meta,\n",
    "             file = \"DDX56_ScoredTCGA.Rdata\"\n",
    "            )"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "R",
   "language": "R",
   "name": "ir"
  },
  "language_info": {
   "codemirror_mode": "r",
   "file_extension": ".r",
   "mimetype": "text/x-r-source",
   "name": "R",
   "pygments_lexer": "r",
   "version": "3.6.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}