Add GO correlation to St subclust QuSAGE

r.scripts/sc_QuSAGE_StSubClust.R

@@ -13,6 +13,7 @@ option_list=list(
   make_option("--r",action="store",default=0.1,type='numeric',help="Resolution"),
   make_option("--st",action="store",default=TRUE,type='logical',help="Analyze Data With Stress Pops Removed?"),
   make_option("--s",action="store",default=5000,type='integer',help="Number Of Cells To Downsample To"),
+  make_option("--go",action="store",default=TRUE,type='logical',help="Correlate to GO St?"),
   make_option("--lgea",action="store",default=TRUE,type='logical',help="Correlate to LGEA St?"),
   make_option("--c2",action="store",default=FALSE,type='logical',help="Correlate to c2?"),
   make_option("--c5",action="store",default=FALSE,type='logical',help="Correlate to c5?")
@@ -81,10 +82,28 @@ rm(labels)
 
 #load QuSAGE genesets
 sets <- NULL
+if (opt$go==TRUE){
+  sets=c(sets,"go")
+}
 if (opt$lgea==TRUE){
   sets=c(sets,"lgea")
 }
 
+gene.set2 <- read_delim("./genesets/DEG_C5.BP.M11704.txt","\t",escape_double=FALSE,trim_ws=TRUE,col_names=TRUE)
+gene.set2 <- as.list(gene.set2[-1,])
+names(gene.set2) <- "Endo"
+gene.set <- c(gene.set2)
+gene.set2 <- read_delim("./genesets/DEG_C5.BP.M10794.txt","\t",escape_double=FALSE,trim_ws=TRUE,col_names=TRUE)
+gene.set2 <- as.list(gene.set2[-1,])
+names(gene.set2) <- "SM"
+gene.set <- c(gene.set,gene.set2)
+gene.set2 <- read_delim("./genesets/DEG_C5.BP.M13024.txt","\t",escape_double=FALSE,trim_ws=TRUE,col_names=TRUE)
+gene.set2 <- as.list(gene.set2[-1,])
+names(gene.set2) <- "Fib"
+gene.set.go <- c(gene.set,gene.set2)
+rm(gene.set)
+rm(gene.set2)
+
 gene.set3 <- read_excel("./genesets/journal.pcbi.1004575.s026.XLSX",skip=2)
 gene.set3 <- gene.set3[,c(3,9)]
 gene.set2 <- gene.set3[gene.set3$`CELL TYPE`=="C1 : Proliferative Fibroblast",]
@@ -120,6 +139,9 @@ gene.set.c5 <- read.gmt("./genesets/c5.all.v6.1.symbols.gmt")
 
 #run QuSAGE
 for (i in 1:Number.Clusters){
+  if (opt$go==TRUE){
+    assign(paste0("qs.results.st.go.",i),qusage(eset,unlist(C[i]),unlist(cC[i]),gene.set.go))
+  }
   if (opt$lgea==TRUE){
     assign(paste0("qs.results.st.lgea.",i),qusage(eset,unlist(C[i]),unlist(cC[i]),gene.set.lgea))
   }
@@ -135,6 +157,20 @@ rm(C)
 rm(cC)
 
 #generate correlation figures
+if (opt$go==TRUE){
+  max.x.rg.go <- 0
+  min.x.rg.go <- 0
+  max.y.rg.go <- 0
+  for (i in 1:Number.Clusters){
+    qs <- get(paste0("qs.results.st.go.",i))
+    if (max(qs$path.mean)>max.x.rg.go){
+      max.x.rg.go <- max(qs$path.mean)
+      if (min(qs$path.mean)<min.x.rg.go){
+        min.x.rg.go <- min(qs$path.mean)
+      }}
+    if (max(qs$path.PDF)>max.y.rg.go){
+      max.y.rg.go <- max(qs$path.PDF)
+}}}
 if (opt$lgea==TRUE){
   max.x.rg.lgea <- 0
   min.x.rg.lgea <- 0
@@ -148,8 +184,15 @@ if (opt$lgea==TRUE){
       }}
     if (max(qs$path.PDF)>max.y.rg.lgea){
       max.y.rg.lgea <- max(qs$path.PDF)
-    }}}
+}}}
 for (i in 1:Number.Clusters){
+  if (opt$go==TRUE){
+    qs <- get(paste0("qs.results.st.go.",i))
+    postscript(paste0("./analysis/",opt$g,"/global",sub.folder,"/correlation/QuSAGE_St_",i,".GO.eps"))
+    plotDensityCurves(qs,path.index=1:3,col=1:numPathways(qs),main=paste0("Cluster ",i),xlim=c(min.x.rg.lgea-0.05,max.x.rg.lgea+0.05),ylim=c(0,50*ceiling(max.y.rg.lgea/50)),xlab="Gene Set Activation")
+    legend("topright",legend=c("Enod","SM","Fib"),lty=1,col=1:numPathways(qs),lwd=1,cex=1,ncol=1)
+    dev.off()
+  }
   if (opt$lgea==TRUE){
     qs <- get(paste0("qs.results.st.lgea.",i))
     postscript(paste0("./analysis/",opt$g,"/global",sub.folder,"/correlation/QuSAGE_St_",i,".LGEA.eps"))
@@ -176,7 +219,10 @@ for (i in 1:Number.Clusters){
     tab <- tab[tab$FDR<0.05,]
     tab <- tab[order(-tab[,2]),]
     write.table(tab,file=paste0("./analysis/",opt$g,"/global",sub.folder,"/correlation/QuSAGE_St_",i,".Table.c5.csv"),row.names=TRUE,col.names=NA,append=FALSE,sep=",")
-  }}
+}}
+rm(max.x.rg.go)
+rm(min.x.rg.go)
+rm(max.y.rg.go)
 rm(max.x.rg.lgea)
 rm(min.x.rg.lgea)
 rm(max.y.rg.lgea)