Add NE and QuSAGE to mus DIY

r.scripts/diy_muPrUr_Epi.R

@@ -7,11 +7,23 @@ library(RColorBrewer)
 library(viridis)
 library(scales)
 library(ComplexHeatmap)
-
+library(autothresholdr)
+library(dplyr)
+library(qusage)
 options(bitmapType="cairo")
 
 setwd("../")
 
+source("./r.scripts/sc-TissueMapper_functions.R")
+
+#load("./genesets/sc10x.epi.id.rda")
+#deg.ne <- FindMarkers(sc10x.epi,ident.1="NE",assay="SCT",slot="data",logfc.threshold=log(1.5,2),test.use="MAST",only.pos=TRUE)
+#deg.ne <- as.character(rownames(deg.ne))
+#convert <- read.delim("./genesets/Ensemble.mus-hum.txt")
+#deg.ne <- deg.ne[deg.ne %in% convert$Human.gene.name]
+#deg.ne <- as.character(convert$Gene.name[match(deg.ne,convert$Human.gene.name)])
+#anchor <- c("Chga","Chgb","Calca")
+
 load("./analysis/sc10x.id.rda")
 load("./analysis/singler_objects.RData")
 
@@ -97,7 +109,6 @@ Idents(sc10x.epi) <- "Cell Type"
 Idents(sc10x.epi.pr) <- "Cell Type"
 Idents(sc10x.epi.ur) <- "Cell Type"
 postscript(paste0("./analysis/vis/diy/UMAP_epi.split.pr.eps"))
-DimPlot(sc10x.epi.pr)+scale_color_brewer(palette="Dark2")+theme_cowplot()
 dev.off()
 postscript(paste0("./analysis/vis/diy/UMAP_epi.split.ur.eps"))
 DimPlot(sc10x.epi.ur)+scale_color_brewer(palette="Dark2")+theme_cowplot()
@@ -154,3 +165,30 @@ write.table(table(sc10x$`Cell Lineage`,sc10x$Region),file="./analysis/vis/diy/Ta
 
 write.table(table(sc10x.epi$`Cell Type`,sc10x.epi$samples),file="./analysis/vis/diy/Table_Epi.Samples.csv",sep=",",quote=FALSE,row.names=TRUE,col.names=NA)
 write.table(table(sc10x.epi$`Cell Type`,sc10x.epi$Region),file="./analysis/vis/diy/Table_Epi.Region.csv",sep=",",quote=FALSE,row.names=TRUE,col.names=NA)
+
+#deg.ne <- rownames(sc10x.epi)[rownames(sc10x.epi) %in% deg.ne]
+#results <- scScore(sc10x.epi,score="NE",geneset=as.list(deg.ne),cut.pt=0.95,anchor=anchor)
+#table(results[[1]]$NE)
+
+E.MTAB.4991_BE <- read.delim("./genesets/BEvsBE.rest.txt", header=TRUE, comment.char="#", stringsAsFactors=FALSE)
+E.MTAB.4991_LE <- read.delim("./genesets/LEvsLE.rest.txt", header=TRUE, comment.char="#", stringsAsFactors=FALSE)
+E.MTAB.4991_OE <- read.delim("./genesets/OEvsOE.rest.txt", header=TRUE, comment.char="#", stringsAsFactors=FALSE)
+
+E.MTAB.4991_BE <- unique(E.MTAB.4991_BE[E.MTAB.4991_BE$Fold.Change > 0 & E.MTAB.4991_BE$FDR.P.val <= 0.05 & (E.MTAB.4991_BE$Group == "Coding" | E.MTAB.4991_BE$Group == "Multiple_Complex") & E.MTAB.4991_BE$Gene.Symbol != "","Gene.Symbol"])
+E.MTAB.4991_BE <- list(E.MTAB.4991_BE)
+names(E.MTAB.4991_BE) <- "Basal"
+E.MTAB.4991_LE <- unique(E.MTAB.4991_LE[E.MTAB.4991_LE$Fold.Change > 0 & E.MTAB.4991_LE$FDR.P.val <= 0.05 & (E.MTAB.4991_LE$Group == "Coding" | E.MTAB.4991_LE$Group == "Multiple_Complex") & E.MTAB.4991_LE$Gene.Symbol != "","Gene.Symbol"])
+E.MTAB.4991_LE <- list(E.MTAB.4991_LE)
+names(E.MTAB.4991_LE) <- "Luminal"
+E.MTAB.4991_OE <- unique(E.MTAB.4991_OE[E.MTAB.4991_OE$Fold.Change > 0 & E.MTAB.4991_OE$FDR.P.val <= 0.05 & (E.MTAB.4991_OE$Group == "Coding" | E.MTAB.4991_OE$Group == "Multiple_Complex") & E.MTAB.4991_OE$Gene.Symbol != "","Gene.Symbol"])
+E.MTAB.4991_OE <- list(E.MTAB.4991_OE)
+names(E.MTAB.4991_OE) <- "LSCmed"
+
+E.MTAB.4991 <- list()
+E.MTAB.4991 <- c(E.MTAB.4991_BE,E.MTAB.4991_LE,E.MTAB.4991_OE)
+
+sc10x.epi$CellType <- sc10x.epi$`Cell Type`
+Idents(sc10x.epi) <- "CellType"
+sc10x.epi$CellType <- Idents(sc10x.epi)
+sc10x.epi$CellType <- factor(sc10x.epi$CellType,levels=c("BE","Ur","VPrLE","DLPrLE","APrLE","ED"))
+results <- scQuSAGE(sc10x.epi,gs=E.MTAB.4991,type="sm",id="CellType",nm="E.MTAB.4991",print="umap")
\ No newline at end of file

r.scripts/huPr_muPr.R

@@ -37,7 +37,7 @@ lin.se <- subset(lin.se,select=!(as.numeric(gsub("_","",substr(colnames(lin.se),
 leu.se <- lin.se
 rm(igd.se)
 
-convert <- read.delim("/work/urology/ghenry/RNA-Seq/SingleCell/PIPELINE/RUN/muPr/genesets/Ensemble.mus-hum.txt")
+convert <- read.delim("./genesets/Ensemble.mus-hum.txt")
 
 if (opt$r == "epi"){
   load("./genesets/sc10x.epi.id.rda")

r.scripts/sc-TissueMapper_functions.R

@@ -847,9 +847,11 @@ scQuSAGE <- function(sc10x,gs,save=FALSE,type,id,ds=0,nm="pops",print="tsne"){
   }
   data <- as.data.frame(as.matrix(GetAssayData(sc10x[,colnames(sc10x) %in% cell.sample])))
   labels <- labels[colnames(sc10x) %in% cell.sample]
-  groups <- sort(unique(labels))
+  #groups <- sort(unique(labels))
+  groups <- paste0("Cluster_",levels(sc10x@active.ident))
   
-  col <- hcl(h=(seq(15,375-375/length(groups),length=length(groups))),c=100,l=65)
+  #col <- hcl(h=(seq(15,375-375/length(groups),length=length(groups))),c=100,l=65)
+  col <- brewer.pal(n=length(groups),name="Dark2")
   
   #Make labels for QuSAGE
   clust <- list()
@@ -970,11 +972,11 @@ scQuSAGE <- function(sc10x,gs,save=FALSE,type,id,ds=0,nm="pops",print="tsne"){
   plot3 <- DimPlot(sc10x,reduction="umap",label=TRUE,repel=TRUE)+theme(legend.position="none")
   if (print=="tsne"){
     postscript(paste0("./analysis/vis/",nm,"/tSNE_",id,"_",nm,".eps"))
-    print(print2)
+    print(plot2)
     dev.off()
   } else if (print=="umap"){
     postscript(paste0("./analysis/vis/",nm,"/UMAP_",id,"_",nm,".eps"))
-    print(print3)
+    print(plot3)
     dev.off()
   } else if (print=="2"){
     postscript(paste0("./analysis/vis/",nm,"/2Vis_",id,"_",nm,".eps"))