QuSAGE ALL populations at one time for VAMC15x

r.scripts/sc-TissueMapper_RUN.VAMC015PrFx.R → r.scripts/sc-TissueMapper_RUN_PrFx_.VAMC015.R

@@ -10,6 +10,7 @@ library(RColorBrewer)
 library(monocle)
 library(dplyr)
 library(viridis)
+library(readxl)
 
 source("../r.scripts/sc-TissueMapper.R")
 
@@ -107,7 +108,7 @@ option_list=list(
   make_option("--st",action="store",default=TRUE,type='logical',help="Remove stressed cells?"),
   make_option("--stg",action="store",default="dws",type='character',help="Geneset to use for stress ID"),
   make_option("--cut.stress",action="store",default=0.9,type='numeric',help="Cutoff for stress score"),
-  make_option("--res.poststress",action="store",default=1,type='numeric',help="Resolution to cluster, post-stress"),
+  make_option("--res.poststress",action="store",default=0.5,type='numeric',help="Resolution to cluster, post-stress"),
   make_option("--cut.ne",action="store",default=0.999,type='numeric',help="Cutoff for NE score")
 )
 opt=parse_args(OptionParser(option_list=option_list))
@@ -177,65 +178,11 @@ gene.set1 <- gene.set1[1]
 gene.set1 <- as.list(gene.set1)
 names(gene.set1) <- "St"
 gene.set <- c(gene.set,gene.set1)
-rm(gene.set1)
-gc()
-min.all <- min(table(sc10x@meta.data[,paste0("res",opt$res.poststress)]))
-results <- scQuSAGE(sc10x,gs=gene.set,res.use=opt$res.poststress,ds=0,nm="Lin",folder="lin")
-sc10x <- results[[1]]
-results.cor.Lin <- results[[2]]
-results.clust.Lin.id <- results[[3]]
-rm(results)
-rm(gene.set)
-
-sc10x <- SetAllIdent(object=sc10x,id="Lin")
-sc10x.Epi <- scSubset(sc10x,i="Lin",g="Epi")
-if (any(levels(sc10x@ident)=="Unknown")){
-  sc10x.St <- scSubset(sc10x,i="Lin",g=c("St","Unknown"))
-} else {
-  sc10x.St <- scSubset(sc10x,i="Lin",g="St")
-}
-sc10x.Epi <- SetAllIdent(object=sc10x.Epi,id=paste0("res",opt$res.poststress))
-sc10x.Epi <- BuildClusterTree(sc10x.Epi,do.reorder=TRUE,reorder.numeric=TRUE,do.plot=FALSE)
-sc10x.Epi <- StashIdent(object=sc10x.Epi,save.name=paste0("res",opt$res.poststress))
-sc10x.St <- SetAllIdent(object=sc10x.St,id=paste0("res",opt$res.poststress))
-sc10x.St <- BuildClusterTree(sc10x.St,do.reorder=TRUE,reorder.numeric=TRUE,do.plot=FALSE)
-sc10x.St <- StashIdent(object=sc10x.St,save.name=paste0("res",opt$res.poststress))
-
-sc10x.Epi <- RunTSNE(object=sc10x.Epi,reduction.use="pca",dims.use=1:pc.use.poststress,do.fast=TRUE)
-postscript(paste0("./analysis/tSNE/epi/tSNE_Sample.eps"))
-plot <- TSNEPlot(object=sc10x.Epi,group.by="samples",pt.size=2.5,do.return=TRUE,vector.friendly=FALSE)
-plot <- plot+theme(axis.text.x=element_text(size=20),axis.text.y=element_text(size=20),axis.title.x=element_text(size=20),axis.title.y=element_text(size=20),legend.text=element_text(size=20))
-plot <- plot+guides(colour=guide_legend(override.aes=list(size=10)))
-plot(plot)
-dev.off()
-postscript(paste0("./analysis/tSNE/epi/tSNE_res",opt$res.poststress,".eps"))
-plot <- TSNEPlot(object=sc10x.Epi,pt.size=5,do.label=TRUE,label.size=10,do.return=TRUE,vector.friendly=FALSE)
-plot <- plot+theme(axis.text.x=element_text(size=20),axis.text.y=element_text(size=20),axis.title.x=element_text(size=20),axis.title.y=element_text(size=20),legend.text=element_text(size=20))
-plot <- plot+guides(colour=guide_legend(override.aes=list(size=10)))
-plot(plot)
-dev.off()
-rm(plot)
-
-sc10x.St <- RunTSNE(object=sc10x.St,reduction.use="pca",dims.use=1:pc.use.poststress,do.fast=TRUE)
-postscript(paste0("./analysis/tSNE/st/tSNE_Sample.eps"))
-plot <- TSNEPlot(object=sc10x.St,group.by="samples",pt.size=2.5,do.return=TRUE,vector.friendly=FALSE)
-plot <- plot+theme(axis.text.x=element_text(size=20),axis.text.y=element_text(size=20),axis.title.x=element_text(size=20),axis.title.y=element_text(size=20),legend.text=element_text(size=20))
-plot <- plot+guides(colour=guide_legend(override.aes=list(size=10)))
-plot(plot)
-dev.off()
-postscript(paste0("./analysis/tSNE/st/tSNE_res",opt$res.poststress,".eps"))
-plot <- TSNEPlot(object=sc10x.St,pt.size=5,do.label=TRUE,label.size=10,do.return=TRUE,vector.friendly=FALSE)
-plot <- plot+theme(axis.text.x=element_text(size=20),axis.text.y=element_text(size=20),axis.title.x=element_text(size=20),axis.title.y=element_text(size=20),legend.text=element_text(size=20))
-plot <- plot+guides(colour=guide_legend(override.aes=list(size=10)))
-plot(plot)
-dev.off()
-rm(plot)
-
 gene.set1 <- read_delim("./genesets/genes.deg.BE.csv",",",escape_double=FALSE,trim_ws=TRUE,col_names=TRUE)
 gene.set1 <- gene.set1[1]
 gene.set1 <- as.list(gene.set1)
 names(gene.set1) <- "BE"
-gene.set <- c(gene.set1)
+gene.set <- c(gene.set,gene.set1)
 gene.set1 <- read_delim("./genesets/genes.deg.LE.csv",",",escape_double=FALSE,trim_ws=TRUE,col_names=TRUE)
 gene.set1 <- gene.set1[1]
 gene.set1 <- as.list(gene.set1)
@@ -251,21 +198,11 @@ gene.set1 <- gene.set1[1]
 gene.set1 <- as.list(gene.set1)
 names(gene.set1) <- "OE_KRT13"
 gene.set <- c(gene.set,gene.set1)
-rm(gene.set1)
-gc()
-min.epi <- min(table(sc10x.Epi@meta.data[,paste0("res",opt$res.poststress)]))
-results <- scQuSAGE(sc10x.Epi,gs=gene.set,res.use=opt$res.poststress,ds=0,nm="Epi.dws.sc",folder="epi")
-sc10x.Epi <- results[[1]]
-results.cor.Epi.dws <- results[[2]]
-results.clust.Epi.dws.id <- results[[3]]
-rm(results)
-rm(gene.set)
-
 gene.set1 <- read_delim("./genesets/genes.deg.Endo.csv",",",escape_double=FALSE,trim_ws=TRUE,col_names=TRUE)
 gene.set1 <- gene.set1[1]
 gene.set1 <- as.list(gene.set1)
 names(gene.set1) <- "Endo"
-gene.set <- c(gene.set1)
+gene.set <- c(gene.set,gene.set1)
 gene.set1 <- read_delim("./genesets/genes.deg.SM.csv",",",escape_double=FALSE,trim_ws=TRUE,col_names=TRUE)
 gene.set1 <- gene.set1[1]
 gene.set1 <- as.list(gene.set1)
@@ -276,116 +213,23 @@ gene.set1 <- gene.set1[1]
 gene.set1 <- as.list(gene.set1)
 names(gene.set1) <- "Fib"
 gene.set <- c(gene.set,gene.set1)
-gene.set1 <- read_delim("./genesets/genes.deg.Leu.csv",",",escape_double=FALSE,trim_ws=TRUE,col_names=TRUE)
-gene.set1 <- gene.set1[1]
-gene.set1 <- as.list(gene.set1)
-names(gene.set1) <- "Leu"
+genes.leu <- read_excel("genesets/40425_2017_215_MOESM1_ESM.xlsx",sheet="S3. Candidate markers")
+leu <- as.list(unique(genes.leu[,2]))$Cell
+leu.l <- leu[-c(1,3,4,7:9,14,15,17:18,20:21,23:30)]
+genes.leu <- genes.leu[unlist(genes.leu[,2]) %in% unlist(leu),]
+genes.leu.l <- genes.leu[unlist(genes.leu[,2]) %in% unlist(leu.l),]
+gene.set1 <- split(genes.leu.l[,1], genes.leu.l[,2])
+gene.set1 <- lapply(gene.set1,unname)
+gene.set1 <- lapply(gene.set1,unlist)
 gene.set <- c(gene.set,gene.set1)
 rm(gene.set1)
 gc()
-min.st <- min(table(sc10x.St@meta.data[,paste0("res",opt$res.poststress)]))
-results <- scQuSAGE(sc10x.St,gs=gene.set,res.use=opt$res.poststress,ds=0,nm="St.dws.sc",folder="st")
-sc10x.St <- results[[1]]
-results.cor.St.go <- results[[2]]
-results.clust.St.go.id <- results[[3]]
+min.all <- min(table(sc10x@meta.data[,paste0("res",opt$res.poststress)]))
+results <- scQuSAGE(sc10x,gs=gene.set,res.use=opt$res.poststress,ds=min.all,nm="Pop",folder="lin")
+sc10x <- results[[1]]
+results.cor.Lin <- results[[2]]
+results.clust.Lin.id <- results[[3]]
 rm(results)
 rm(gene.set)
 
-sc10x.Epi.NE <- scNE(sc10x.Epi,neg="dws",cut=opt$cut.ne)
-
-sc10x <- scMerge(sc10x,sc10x.Epi,sc10x.St,i.1="Epi.dws.sc",i.2="St.dws.sc",nm="Merge_Epi.dws.sc_St.dws.sc")
-
-sc10x <- scMerge(sc10x,sc10x,sc10x.Epi.NE,i.1="Merge_Epi.dws.sc_St.dws.sc",i.2="NE",nm="Merge_Epi.dws.sc_St.dws.sc_NE")
-
-sc10x.Epi <- scMerge(sc10x.Epi,sc10x.Epi,sc10x.Epi.NE,i.1="Epi.dws.sc",i.2="NE",nm="Epi.dws.sc_NE")
-
-sc10x <- SetAllIdent(object=sc10x,id="Merge_Epi.dws.sc_St.dws.sc")
-sc10x@ident <- factor(sc10x@ident,levels=c("BE","LE","OE_SCGB","OE_KRT13","Fib","SM","Endo","Leu"))
-postscript("./analysis/tSNE/FINAL/tSNE_FINAL.eps")
-plot <- TSNEPlot(object=sc10x,pt.size=2.5,do.return=TRUE,vector.friendly=FALSE)
-plot <- plot+theme(axis.text.x=element_text(size=20),axis.text.y=element_text(size=20),axis.title.x=element_text(size=20),axis.title.y=element_text(size=20),legend.text=element_text(size=20))
-plot <- plot+guides(colour=guide_legend(override.aes=list(size=10)))
-plot(plot)
-dev.off()
-
-scTables(sc10x,i.1="samples",i.2="Merge_Epi.dws.sc_St.dws.sc")
-scTables(sc10x,i.1="samples",i.2="Merge_Epi.dws.sc_St.dws.sc_NE")
-scTables(sc10x,i.1="Merge_Epi.dws.sc_St.dws.sc_NE",i.2="Merge_Epi.dws.sc_St.dws.sc")
-
-sctSNECustCol(sc10x,i="Lin",bl="Epi",rd="St",file="D17")
-sctSNECustCol(sc10x,i="Merge_Epi.dws.sc_St.dws.sc",bl=c("BE","LE","OE_SCGB","OE_KRT13"),rd=c("Fib","SM","Endo","Leu"),file="D17")
-sctSNECustCol(sc10x.Epi,i="Epi.dws.sc",bl=c("BE","LE","OE_SCGB","OE_KRT13"),rd="",file="D17")
-sctSNECustCol(sc10x.St,i="St.dws.sc",bl="",rd=c("Fib","SM","Endo","Leu"),file="D17")
-
-sctSNEbwCol(sc10x,i=paste0("res",opt$res.poststress),file="ALL",files="D17")
-sctSNEbwCol(sc10x.Epi,i=paste0("res",opt$res.poststress),file="Epi",files="D17")
-sctSNEbwCol(sc10x.St,i=paste0("res",opt$res.poststress),file="St",files="D17")
-sctSNEbwCol(sc10x,i="Merge_Epi.dws.sc_St.dws.sc",file="ALL",files="D17")
-sctSNEbwCol(sc10x.Epi,i="Epi.dws.sc",file="Epi",files="D17")
-sctSNEbwCol(sc10x.St,i="St.dws.sc",file="St",files="D17")
-
-for (g in c("Epi","St","Unknown")){
-  sctSNEHighlight(sc10x,i="Lin",g=g,file="D17")
-}
-for (g in c("BE","LE","OE_SCGB","OE_KRT13")){
-  sctSNEHighlight(sc10x,i="Merge_Epi.dws.sc_St.dws.sc",g=g,file="D17")
-  sctSNEHighlight(sc10x.Epi,i="Epi.dws.c",g=g,file="D17")
-}
-sctSNEHighlight(sc10x.Epi.NE,i="NE",g="NE",file="D17")
-for (g in c("Fib","SM","Endo","Leu")){
-  sctSNEHighlight(sc10x,i="Merge_Epi.dws.sc_St.dws.sc",g=g,file="D17")
-  sctSNEHighlight(sc10x.St,i="St.dws.sc",g=g,file="D17")
-}
-for (g in c("D17PrPzF_BE","D17PrPzF_LE","D17PrPzF_OE","D17PrPzF_FMSt")){
-  sctSNEHighlight(sc10x,i="samples",g=g,file="D17")
-}
-rm(i)
-rm(g)
-
-postscript(paste0("./analysis/diy/Feature.eps"))
-FeaturePlot(sc10x,features.plot=c("PECAM1","CD200","PDPN","EPCAM","DPP4","PSCA","VIM","KRT5","KLK3"),cols.use=c("grey","darkred"),reduction.use="tsne")
-dev.off()
-
-postscript(paste0("./analysis/diy/Feature.KLK3.eps"))
-VlnPlot(sc10x,features.plot="KLK3",group.by="Lin",size.title.use=20,point.size.use=0.05)
-dev.off()
-
-sc10x.LE <- scSubset(sc10x,i="Merge_Epi.dws.sc_St.dws.sc",g="LE")
-genes.deg.LE.Enza.vs.Ctrl <- scDEG(sc10x.LE,i="samples",g.1="VAMC015PrFx_Enza_Via",g.2="VAMC015PrFx_Ctrl_Via",pct=0.25,t=0)
-genes.deg.LE.ODM.vs.Ctrl <- scDEG(sc10x.LE,i="samples",g.1="VAMC015PrFx_ODM_Via",g.2="VAMC015PrFx_Ctrl_Via",pct=0.25,t=0)
-genes.deg.LE.ODM.vs.Enza <- scDEG(sc10x.LE,i="samples",g.1="VAMC015PrFx_ODM_Via",g.2="VAMC015PrFx_Enza_Via",pct=0.25,t=0)
-
-genes.deg.LE.Ctrl.vs.Enza <- scDEG(sc10x.LE,i="samples",g.1="VAMC015PrFx_Ctrl_Via",g.2="VAMC015PrFx_Enza_Via",pct=0.25,t=0)
-genes.deg.LE.Ctrl.vs.ODM <- scDEG(sc10x.LE,i="samples",g.1="VAMC015PrFx_Ctrl_Via",g.2="VAMC015PrFx_ODM_Via",pct=0.25,t=0)
-genes.deg.LE.Enza.vs.ODM <- scDEG(sc10x.LE,i="samples",g.1="VAMC015PrFx_Enza_Via",g.2="VAMC015PrFx_ODM_Via",pct=0.25,t=0)
-
-postscript(paste0("./analysis/diy/Feature.LE.KLK3.eps"))
-VlnPlot(sc10x.LE,features.plot="KLK3",group.by="samples",size.title.use=20,point.size.use=0.05)
-dev.off()
-
-sc10x.Epi <- scSubset(sc10x,i="Lin",g="Epi")
-genes.deg.Epi.Enza.vs.Ctrl <- scDEG(sc10x.Epi,i="samples",g.1="VAMC015PrFx_Enza_Via",g.2="VAMC015PrFx_Ctrl_Via",pct=0.25,t=0)
-genes.deg.Epi.ODM.vs.Ctrl <- scDEG(sc10x.Epi,i="samples",g.1="VAMC015PrFx_ODM_Via",g.2="VAMC015PrFx_Ctrl_Via",pct=0.25,t=0)
-genes.deg.Epi.ODM.vs.Enza <- scDEG(sc10x.Epi,i="samples",g.1="VAMC015PrFx_ODM_Via",g.2="VAMC015PrFx_Enza_Via",pct=0.25,t=0)
-
-genes.deg.Epi.Ctrl.vs.Enza <- scDEG(sc10x.Epi,i="samples",g.1="VAMC015PrFx_Ctrl_Via",g.2="VAMC015PrFx_Enza_Via",pct=0.25,t=0)
-genes.deg.Epi.Ctrl.vs.ODM <- scDEG(sc10x.Epi,i="samples",g.1="VAMC015PrFx_Ctrl_Via",g.2="VAMC015PrFx_ODM_Via",pct=0.25,t=0)
-genes.deg.Epi.Enza.vs.ODM <- scDEG(sc10x.Epi,i="samples",g.1="VAMC015PrFx_Enza_Via",g.2="VAMC015PrFx_ODM_Via",pct=0.25,t=0)
-
-postscript(paste0("./analysis/diy/Feature.Epi.KLK3.eps"))
-VlnPlot(sc10x.Epi,features.plot="KLK3",group.by="samples",size.title.use=20,point.size.use=0.05)
-dev.off()
-
-for (i in c("genes.deg.LE.Ctrl.vs.Enza","genes.deg.LE.Ctrl.vs.ODM","genes.deg.LE.ODM.vs.Enza","genes.deg.LE.Enza.vs.ODM")){
-write.table(get(i),file=paste0(i,".csv"),row.names=TRUE,col.names=NA,append=FALSE,sep=",")
-}
-
-save(list=ls(pattern="sc10x.Stress"),file="./analysis/sc10x.Stress.Rda")
-rm(list=ls(pattern="sc10x.Stress"))
-save(list=ls(pattern="sc10x.Epi"),file="./analysis/sc10x.Epi.Rda")
-rm(list=ls(pattern="^sc10x.Epi"))
-save(list=ls(pattern="sc10x.St"),file="./analysis/sc10x.St.Rda")
-rm(list=ls(pattern="sc10x.St"))
-save(list=ls(pattern="^sc10x"),file="./analysis/sc10x.Rda")
-rm(list=ls(pattern="^sc10x"))
-save.image(file="./analysis/Data.RData")
+sc10x.Epi.NE <- scNE(sc10x.Epi,neg="dws",cut=opt$cut.ne)
\ No newline at end of file