Add ridge and violins for top 10 NE DEGs and update DPrF's

r.scripts/sc-TissueMapper_RUN.D17.R

@@ -98,7 +98,7 @@ option_list=list(
   make_option("--stg",action="store",default="go",type='character',help="Geneset to use for stress ID"),
   #make_option("--hpc.poststress",action="store",default=0.85,type='numeric',help="Max variance cutoff for PCs to use, post-stress"),
   make_option("--res.poststress",action="store",default=0.2,type='numeric',help="Resolution to cluster, post-stress"),
-  make_option("--ds",action="store",default=0,type='integer',help="Number of cells to downsample"),
+  make_option("--ds",action="store",default=10000,type='integer',help="Number of cells to downsample"),
   #make_option("--hpc.epi",action="store",default=0.85,type='numeric',help="Max variance cutoff for PCs to use, Epi"),
   make_option("--res.epi",action="store",default=0.2,type='numeric',help="Resolution to cluster, Epi"),
   #make_option("--hpc.st",action="store",default=0.85,type='numeric',help="Max variance cutoff for PCs to use, St"),
@@ -164,6 +164,7 @@ gene.set1 <- as.list(gene.set1)
 names(gene.set1) <- "Epi"
 gene.set <- c(gene.set,gene.set1)
 rm(gene.set1)
+gc()
 results <- scQuSAGE(sc10x,gs=gene.set,res.use=opt$res.poststress,ds=opt$ds,nm="Lin",folder="lin")
 sc10x <- results[[1]]
 results.cor.Lin <- results[[2]]
@@ -193,6 +194,7 @@ gene.set1 <- as.list(gene.set1)
 names(gene.set1) <- "OE"
 gene.set <- c(gene.set,gene.set1)
 rm(gene.set1)
+gc()
 results <- scQuSAGE(sc10x.Epi,gs=gene.set,res.use=opt$res.epi,ds=opt$ds,nm="Epi.dws",folder="epi")
 sc10x.Epi <- results[[1]]
 results.cor.Epi.dws <- results[[2]]
@@ -240,6 +242,7 @@ gene.set1 <- as.list(gene.set1[-1,])
 names(gene.set1) <- "Leu"
 gene.set <- c(gene.set,gene.set1)
 rm(gene.set1)
+gc()
 results <- scQuSAGE(sc10x.St,gs=gene.set,res.use=opt$res.st,ds=opt$ds,nm="St.go",folder="st")
 sc10x.St <- results[[1]]
 results.cor.St.go <- results[[2]]
@@ -336,6 +339,14 @@ postscript("./analysis/deg/Heatmap.Stress.eps",paper="special",width=10,height=5
 plot(DoHeatmap(sc10x.Stress,genes.use=rownames(genes.deg.Stress),slim.col.label=TRUE,group.label.rot=TRUE,group.spacing=0.25,cex.row=2.5))
 dev.off()
 
+sc10x.Epi.NE <- SetAllIdent(object=sc10x.Epi.NE,id="NE")
+postscript("./analysis/deg/Violin.NE.eps",paper="special",width=10,height=5,horizontal=FALSE)
+plot(VlnPlot(sc10x.Epi.NE,features.plot=rownames(genes.deg.NE)[1:10],nCol=5,point.size.use=0.1,size.title.use=15,x.lab.rot=TRUE))
+dev.off()
+postscript("./analysis/deg/Ridge.NE.eps",paper="special",width=10,height=5,horizontal=FALSE)
+plot(RidgePlot(sc10x.Epi.NE,features.plot=rownames(genes.deg.NE)[1:10],nCol=5,size.x.use=10,size.title.use=15))
+dev.off()
+
 for (i in ls(pattern="^genes.deg")){
   write.table(get(i),file=paste0("./analysis/deg/",i,".csv"),row.names=TRUE,col.names=NA,append=FALSE,sep=",")
 }

r.scripts/sc-TissueMapper_RUN.D2.R

@@ -101,7 +101,7 @@ option_list=list(
   make_option("--stg",action="store",default="dws",type='character',help="Geneset to use for stress ID"),
   #make_option("--hpc.poststress",action="store",default=0.85,type='numeric',help="Max variance cutoff for PCs to use, post-stress"),
   make_option("--res.poststress",action="store",default=0.2,type='numeric',help="Resolution to cluster, post-stress"),
-  make_option("--ds",action="store",default=25000,type='integer',help="Number of cells to downsample"),
+  make_option("--ds",action="store",default=10000,type='integer',help="Number of cells to downsample"),
   #make_option("--hpc.epi",action="store",default=0.85,type='numeric',help="Max variance cutoff for PCs to use, Epi"),
   make_option("--res.epi",action="store",default=0.2,type='numeric',help="Resolution to cluster, Epi"),
   #make_option("--hpc.st",action="store",default=0.85,type='numeric',help="Max variance cutoff for PCs to use, St"),
@@ -153,6 +153,7 @@ rm(results)
 
 rm(sc10x.D17)
 rm(sc10x.D27)
+save.image(file="./analysis/Aligned.RData")
 
 #results <- scPC(sc10x,lx=opt$lx,hx=opt$hx,ly=opt$ly,cc=opt$cc,pc=opt$pc,hpc=opt$hpc,file="pre.stress",cca=TRUE)
 #sc10x <- results[[1]]
@@ -160,11 +161,11 @@ rm(sc10x.D27)
 #pc.use <- results[[3]]
 #rm(results)
 
-pc.use <- 25
+pc.use <- 15
 sc10x <- scCluster(sc10x,pc.use=pc.use,res.use=opt$res.prestress,folder="pre.stress",red="cca.aligned")
 
 if (opt$st==TRUE){
-  results <- scStress(sc10x,stg=opt$stg,res.use=opt$res.prestress,pc.use=pc.use)
+  results <- scStress(sc10x,stg=opt$stg,res.use=opt$res.prestress,pc.use=pc.use,cut=0.95)
   sc10x <- results[[1]]
   counts.cell.filtered.stress <- results[[2]]
   sc10x.Stress <- results[[3]]
@@ -190,6 +191,7 @@ gene.set1 <- as.list(gene.set1)
 names(gene.set1) <- "Epi"
 gene.set <- c(gene.set,gene.set1)
 rm(gene.set1)
+gc()
 results <- scQuSAGE(sc10x,gs=gene.set,res.use=opt$res.poststress,ds=opt$ds,nm="Lin",folder="lin")
 sc10x <- results[[1]]
 results.cor.Lin <- results[[2]]
@@ -227,6 +229,7 @@ gene.set1 <- as.list(gene.set1)
 names(gene.set1) <- "OE_KRT13"
 gene.set <- c(gene.set,gene.set1)
 rm(gene.set1)
+gc()
 results <- scQuSAGE(sc10x.Epi,gs=gene.set,res.use=opt$res.epi,ds=opt$ds,nm="Epi.dws.sc",folder="epi")
 sc10x.Epi <- results[[1]]
 results.cor.Epi.dws.sc <- results[[2]]
@@ -250,6 +253,7 @@ gene.set1 <- as.list(gene.set1)
 names(gene.set1) <- "CGC"
 gene.set<- c(gene.set,gene.set1)
 rm(gene.set1)
+gc()
 results.cor.Epi.lgea <- scQuSAGEsm(sc10x.Epi,gs=gene.set,ds=opt$ds,nm="Epi.dws.sc",folder="lgea")
 rm(gene.set)
 
@@ -276,6 +280,7 @@ gene.set1 <- as.list(gene.set1)
 names(gene.set1) <- "Leu"
 gene.set <- c(gene.set,gene.set1)
 rm(gene.set1)
+gc()
 results <- scQuSAGE(sc10x.St,gs=gene.set,res.use=opt$res.st,ds=opt$ds,nm="St.dws.sc",folder="st")
 sc10x.St <- results[[1]]
 results.cor.St.dws.sc <- results[[2]]
@@ -283,7 +288,7 @@ results.clust.St.dws.sc.id <- results[[3]]
 rm(results)
 rm(gene.set)
 
-sc10x.Epi.NE <- scNE(sc10x.Epi,neg="dws")
+sc10x.Epi.NE <- scNE(sc10x.Epi,neg="dws",cut=0.95)
 
 sc10x.Epi <- scMergeSubClust(sc10x.Epi,i="Epi.dws.sc",g=c("BE","LE","OE_SCGB","OE_KRT13"),nm="Merge")
 
@@ -372,6 +377,14 @@ postscript("./analysis/deg/Heatmap.Stress.eps",paper="special",width=10,height=5
 plot(DoHeatmap(sc10x.Stress,genes.use=rownames(genes.deg.Stress),slim.col.label=TRUE,group.label.rot=TRUE,group.spacing=0.25,cex.row=2.5))
 dev.off()
 
+sc10x.Epi.NE <- SetAllIdent(object=sc10x.Epi.NE,id="NE")
+postscript("./analysis/deg/Violin.NE.eps",paper="special",width=10,height=5,horizontal=FALSE)
+plot(VlnPlot(sc10x.Epi.NE,features.plot=rownames(genes.deg.NE)[1:10],nCol=5,point.size.use=0.1,size.title.use=15,x.lab.rot=TRUE))
+dev.off()
+postscript("./analysis/deg/Ridge.NE.eps",paper="special",width=10,height=5,horizontal=FALSE)
+plot(RidgePlot(sc10x.Epi.NE,features.plot=rownames(genes.deg.NE)[1:10],nCol=5,size.x.use=10,size.title.use=15))
+dev.off()
+
 for (i in ls(pattern="^genes.deg")){
   write.table(get(i),file=paste0("./analysis/deg/",i,".csv"),row.names=TRUE,col.names=NA,append=FALSE,sep=",")
 }

r.scripts/sc-TissueMapper_RUN.D3.R

@@ -101,7 +101,7 @@ option_list=list(
   make_option("--stg",action="store",default="dws",type='character',help="Geneset to use for stress ID"),
   #make_option("--hpc.poststress",action="store",default=0.85,type='numeric',help="Max variance cutoff for PCs to use, post-stress"),
   make_option("--res.poststress",action="store",default=0.2,type='numeric',help="Resolution to cluster, post-stress"),
-  make_option("--ds",action="store",default=25000,type='integer',help="Number of cells to downsample"),
+  make_option("--ds",action="store",default=10000,type='integer',help="Number of cells to downsample"),
   #make_option("--hpc.epi",action="store",default=0.85,type='numeric',help="Max variance cutoff for PCs to use, Epi"),
   make_option("--res.epi",action="store",default=0.2,type='numeric',help="Resolution to cluster, Epi"),
   #make_option("--hpc.st",action="store",default=0.85,type='numeric',help="Max variance cutoff for PCs to use, St"),
@@ -168,6 +168,7 @@ rm(results)
 rm(sc10x.D17)
 rm(sc10x.D27)
 rm(sc10x.D35)
+save.image(file="./analysis/Aligned.RData")
 
 #results <- scPC(sc10x,lx=opt$lx,hx=opt$hx,ly=opt$ly,cc=opt$cc,pc=opt$pc,hpc=opt$hpc,file="pre.stress",cca=TRUE)
 #sc10x <- results[[1]]
@@ -175,11 +176,11 @@ rm(sc10x.D35)
 #pc.use <- results[[3]]
 #rm(results)
 
-pc.use <- 25
+pc.use <- 15
 sc10x <- scCluster(sc10x,pc.use=pc.use,res.use=opt$res.prestress,folder="pre.stress",red="cca.aligned")
 
 if (opt$st==TRUE){
-  results <- scStress(sc10x,stg=opt$stg,res.use=opt$res.prestress,pc.use=pc.use)
+  results <- scStress(sc10x,stg=opt$stg,res.use=opt$res.prestress,pc.use=pc.use,cut=0.95)
   sc10x <- results[[1]]
   counts.cell.filtered.stress <- results[[2]]
   sc10x.Stress <- results[[3]]
@@ -205,6 +206,7 @@ gene.set1 <- as.list(gene.set1)
 names(gene.set1) <- "Epi"
 gene.set <- c(gene.set,gene.set1)
 rm(gene.set1)
+gc()
 results <- scQuSAGE(sc10x,gs=gene.set,res.use=opt$res.poststress,ds=opt$ds,nm="Lin",folder="lin")
 sc10x <- results[[1]]
 results.cor.Lin <- results[[2]]
@@ -242,6 +244,7 @@ gene.set1 <- as.list(gene.set1)
 names(gene.set1) <- "OE_KRT13"
 gene.set <- c(gene.set,gene.set1)
 rm(gene.set1)
+gc()
 results <- scQuSAGE(sc10x.Epi,gs=gene.set,res.use=opt$res.epi,ds=opt$ds,nm="Epi.dws.sc",folder="epi")
 sc10x.Epi <- results[[1]]
 results.cor.Epi.dws.sc <- results[[2]]
@@ -265,6 +268,7 @@ gene.set1 <- as.list(gene.set1)
 names(gene.set1) <- "CGC"
 gene.set<- c(gene.set,gene.set1)
 rm(gene.set1)
+gc()
 results.cor.Epi.lgea <- scQuSAGEsm(sc10x.Epi,gs=gene.set,ds=opt$ds,nm="Epi.dws.sc",folder="lgea")
 rm(gene.set)
 
@@ -291,6 +295,7 @@ gene.set1 <- as.list(gene.set1)
 names(gene.set1) <- "Leu"
 gene.set <- c(gene.set,gene.set1)
 rm(gene.set1)
+gc()
 results <- scQuSAGE(sc10x.St,gs=gene.set,res.use=opt$res.st,ds=opt$ds,nm="St.dws.sc",folder="st")
 sc10x.St <- results[[1]]
 results.cor.St.dws.sc <- results[[2]]
@@ -298,7 +303,7 @@ results.clust.St.dws.sc.id <- results[[3]]
 rm(results)
 rm(gene.set)
 
-sc10x.Epi.NE <- scNE(sc10x.Epi,neg="dws")
+sc10x.Epi.NE <- scNE(sc10x.Epi,neg="dws"cut=0.95)
 
 sc10x.Epi <- scMergeSubClust(sc10x.Epi,i="Epi.dws.sc",g=c("BE","LE","OE_SCGB","OE_KRT13"),nm="Merge")
 
@@ -387,6 +392,14 @@ postscript("./analysis/deg/Heatmap.Stress.eps",paper="special",width=10,height=5
 plot(DoHeatmap(sc10x.Stress,genes.use=rownames(genes.deg.Stress),slim.col.label=TRUE,group.label.rot=TRUE,group.spacing=0.25,cex.row=2.5))
 dev.off()
 
+sc10x.Epi.NE <- SetAllIdent(object=sc10x.Epi.NE,id="NE")
+postscript("./analysis/deg/Violin.NE.eps",paper="special",width=10,height=5,horizontal=FALSE)
+plot(VlnPlot(sc10x.Epi.NE,features.plot=rownames(genes.deg.NE)[1:10],nCol=5,point.size.use=0.1,size.title.use=15,x.lab.rot=TRUE))
+dev.off()
+postscript("./analysis/deg/Ridge.NE.eps",paper="special",width=10,height=5,horizontal=FALSE)
+plot(RidgePlot(sc10x.Epi.NE,features.plot=rownames(genes.deg.NE)[1:10],nCol=5,size.x.use=10,size.title.use=15))
+dev.off()
+
 for (i in ls(pattern="^genes.deg")){
   write.table(get(i),file=paste0("./analysis/deg/",i,".csv"),row.names=TRUE,col.names=NA,append=FALSE,sep=",")
 }