From 9d0c65573f782e59c813fd0bbb08faabfc912ae4 Mon Sep 17 00:00:00 2001
From: "Gervaise H. Henry" <gervaise.henry@utsouthwestern.edu>
Date: Wed, 17 Oct 2018 18:07:40 -0500
Subject: [PATCH] Fix CustCol to handle Unknown better and fix PdPbPc.Pb code

---
 r.scripts/sc-TissueMapper_RUN.PdPbPc.Pb.R | 292 ++++------------------
 1 file changed, 51 insertions(+), 241 deletions(-)

diff --git a/r.scripts/sc-TissueMapper_RUN.PdPbPc.Pb.R b/r.scripts/sc-TissueMapper_RUN.PdPbPc.Pb.R
index 8302154..abb3b1e 100644
--- a/r.scripts/sc-TissueMapper_RUN.PdPbPc.Pb.R
+++ b/r.scripts/sc-TissueMapper_RUN.PdPbPc.Pb.R
@@ -116,7 +116,7 @@ if (opt$lm==0){opt$lm=-Inf}
 
 sc10x <- scLoad("Pb")
 
-sc10x <- scSubset(sc10x,i="Glandular",g="Glandular")
+sc10x <- scSubset(sc10x,i="Tumor",g="Tumor")
 
 if (opt$cc==TRUE){
   results <- scCellCycle(sc10x)
@@ -145,18 +145,20 @@ if (opt$cc==TRUE){
 }
 gc()
 
-sc10x.B327 <- scSubset(sc10x,"samples","BPH327PrGF_Via")
-sc10x.B340 <- scSubset(sc10x,"samples","BPH340PrGF_Via")
-sc10x.B342 <- scSubset(sc10x,"samples","BPH342PrF_Via")
+sc10x.C11 <- scSubset(sc10x,"samples","VAMC011PrBlF_Via")
+sc10x.C12 <- scSubset(sc10x,"samples","VAMC012PrRdF_Via")
+sc10x.C13 <- scSubset(sc10x,"samples","VAMC013PrRdF_Via")
+sc10x.C14 <- scSubset(sc10x,"samples","VAMC014PrRdF_Via")
 
-results <- sc3CCA(sc10x.D17,sc10x.D27,sc10x.D35,"B327","B340","B342",cc=opt$cc,cca=opt$cca,acca=opt$acca,lx=opt$lx,hx=opt$hx,ly=opt$ly)
+results <- sc4CCA(sc10x.C11,sc10x.C12,sc10x.C13,sc10x.C14,"C11","C12","C13","C14",cc=opt$cc,cca=opt$cca,acca=opt$acca,lx=opt$lx,hx=opt$hx,ly=opt$ly)
 sc10x <- results[[1]]
 genes.hvg.cca <- results[[2]]
 rm(results)
 
-rm(sc10x.B327)
-rm(sc10x.B340)
-rm(sc10x.B342)
+rm(sc10x.C11)
+rm(sc10x.C12)
+rm(sc10x.C13)
+rm(sc10x.C14)
 
 sc10x <- scCluster(sc10x,pc.use=opt$acca,res.use=opt$res.prestress,folder="pre.stress",red="cca.aligned")
 
@@ -170,11 +172,13 @@ if (opt$st==TRUE){
   sc10x <- scCluster(sc10x,pc.use=opt$acca,res.use=opt$res.poststress,folder="post.stress",red="cca.aligned")
 }
 
-gene.set1 <- read_delim("./genesets/DEG_Epi_5FC.txt","\t",escape_double=FALSE,trim_ws=TRUE,col_names=FALSE)
+gene.set1 <- read_delim("./genesets/genes.deg.Epi.csv",",",escape_double=FALSE,trim_ws=TRUE,col_names=TRUE)
+gene.set1 <- gene.set1[1]
 gene.set1 <- as.list(gene.set1)
 names(gene.set1) <- "Epi"
 gene.set <- c(gene.set1)
-gene.set1 <- read_delim("./genesets/DEG_FMSt_5FC.txt","\t",escape_double=FALSE,trim_ws=TRUE,col_names=FALSE)
+gene.set1 <- read_delim("./genesets/genes.deg.St.csv",",",escape_double=FALSE,trim_ws=TRUE,col_names=TRUE)
+gene.set1 <- gene.set1[1]
 gene.set1 <- as.list(gene.set1)
 names(gene.set1) <- "St"
 gene.set <- c(gene.set,gene.set1)
@@ -217,70 +221,36 @@ plot(plot)
 dev.off()
 rm(plot)
 
-gene.set1 <- read_delim("./genesets/DEG_BE_5FC.txt","\t",escape_double=FALSE,trim_ws=TRUE,col_names=FALSE)
+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.set1 <- read_delim("./genesets/DEG_LE_5FC.txt","\t",escape_double=FALSE,trim_ws=TRUE,col_names=FALSE)
+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)
 names(gene.set1) <- "LE"
 gene.set <- c(gene.set,gene.set1)
-gene.set1 <- read_delim("./genesets/DEG_OE_5FC.txt","\t",escape_double=FALSE,trim_ws=TRUE,col_names=FALSE)
+gene.set1 <- read_delim("./genesets/genes.deg.OE1.csv",",",escape_double=FALSE,trim_ws=TRUE,col_names=TRUE)
+gene.set1 <- gene.set1[1]
 gene.set1 <- as.list(gene.set1)
-names(gene.set1) <- "OE"
+names(gene.set1) <- "OE_SCGB"
+gene.set <- c(gene.set,gene.set1)
+gene.set1 <- read_delim("./genesets/genes.deg.OE2.csv",",",escape_double=FALSE,trim_ws=TRUE,col_names=TRUE)
+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=min.epi,nm="Epi.dws",folder="epi")
+results <- scQuSAGE(sc10x.Epi,gs=gene.set,res.use=opt$res.poststress,ds=min.epi,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)
 
-sc10x.Epi <- SetAllIdent(object=sc10x.Epi,id=paste0("res",opt$res.poststress))
-OE <- levels(factor(sc10x.Epi@ident[sc10x.Epi@meta.data$Epi.dws=="OE"]))
-OE.cells <- NULL
-for (i in 1:length(OE)){
-  OE.cells[i] <- list(names(sc10x.Epi@ident[sc10x.Epi@ident==OE[i]]))
-}
-sc10x.Epi <- SetAllIdent(object=sc10x.Epi,id="Epi.dws")
-for (i in 1:length(OE)){
-  sc10x.Epi <- SetIdent(object=sc10x.Epi,cells.use=unlist(OE.cells[i]),ident.use=paste0("OE",i))
-}
-sc10x.Epi <- StashIdent(object=sc10x.Epi,save.name="Epi.dws.sub")
-postscript("./analysis/tSNE/epi/tSNE_Epi.dws.sub.eps")
-plot <- TSNEPlot(object=sc10x.Epi,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()
-rm(plot)
-rm(OE)
-rm(OE.cells)
-rm(i)
-
-gene.set1 <- read_csv("./genesets/Basal cells-signature-genes.csv")
-gene.set1 <- gene.set1[,2]
-gene.set1 <- as.list(gene.set1)
-names(gene.set1) <- "BC"
-gene.set <- c(gene.set1)
-gene.set1 <- read_csv("./genesets/Normal AT2 cells-signature-genes.csv")
-gene.set1 <- gene.set1[,2]
-gene.set1 <- as.list(gene.set1)
-names(gene.set1) <- "AT2"
-gene.set <- c(gene.set,gene.set1)
-gene.set1 <- read_csv("./genesets/Club_Goblet cells-signature-genes.csv")
-gene.set1 <- gene.set1[,2]
-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=min.epi,nm="Epi.dws.sub",folder="lgea")
-rm(gene.set)
-
 sc10x.St <- RunTSNE(object=sc10x.St,reduction.use="cca.aligned",dims.use=1:opt$acca,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)
@@ -296,33 +266,37 @@ plot(plot)
 dev.off()
 rm(plot)
 
-gene.set1 <- read_delim("./genesets/DEG_C5.BP.M11704.txt","\t",escape_double=FALSE,trim_ws=TRUE,col_names=TRUE)
-gene.set1 <- as.list(gene.set1[-1,])
+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.set1 <- read_delim("./genesets/DEG_C5.BP.M10794.txt","\t",escape_double=FALSE,trim_ws=TRUE,col_names=TRUE)
-gene.set1 <- as.list(gene.set1[-1,])
+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)
 names(gene.set1) <- "SM"
 gene.set <- c(gene.set,gene.set1)
-gene.set1 <- read_delim("./genesets/DEG_C5.BP.M13024.txt","\t",escape_double=FALSE,trim_ws=TRUE,col_names=TRUE)
-gene.set1 <- as.list(gene.set1[-1,])
+gene.set1 <- read_delim("./genesets/genes.deg.Fib.csv",",",escape_double=FALSE,trim_ws=TRUE,col_names=TRUE)
+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/DEG_C5.BP.M10124.txt","\t",escape_double=FALSE,trim_ws=TRUE,col_names=TRUE)
-gene.set1 <- as.list(gene.set1[-1,])
+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"
 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=0.2,ds=min.st,nm="St.go",folder="st")
+results <- scQuSAGE(sc10x.St,gs=gene.set,res.use=opt$res.poststress,ds=min.st,nm="St.dws.sc",folder="st")
 sc10x.St <- results[[1]]
 results.cor.St.go <- results[[2]]
 results.clust.St.go.id <- results[[3]]
 rm(results)
 rm(gene.set)
 
-sc10x.Epi.NE <- scNE(sc10x.Epi,neg="EurUro",cut=opt$cut.ne)
+sc10x.Epi.NE <- scNE(sc10x.Epi,neg="dws",cut=opt$cut.ne)
 
 sc10x.Epi <- scMergeSubClust(sc10x.Epi,i="Epi.dws.sub",g=c("BE","LE"),nm="Merge")
 
@@ -338,64 +312,14 @@ sc10x <- SetAllIdent(object=sc10x,id="Lin")
 sc10x <- SetIdent(object=sc10x,cells.use=names(sc10x@ident[sc10x@ident %in% c("St","Unknown")]),ident.use="St")
 sc10x <- StashIdent(object=sc10x,save.name="mLin")
 
-gene.orthog <- read.delim("./genesets/Ensemble.mus-hum.txt")
-gene.set1 <- read_csv("./genesets/SupTab3_Consensus_Sigs.csv",skip=6)
-gene.set2 <- as.data.frame(gene.set1$Basal[!is.na(gene.set1$Basal)])
-colnames(gene.set2) <- "genes"
-gene.set2 <- as.data.frame(merge(gene.set2,gene.orthog,by.x="genes",by.y="Gene.name")[,4])
-gene.set2 <- as.list(gene.set2)
-names(gene.set2) <- "Basal"
-gene.set <- c(gene.set2)
-gene.set2 <- as.data.frame(gene.set1$Club[!is.na(gene.set1$Club)])
-colnames(gene.set2) <- "genes"
-gene.set2 <- as.data.frame(merge(gene.set2,gene.orthog,by.x="genes",by.y="Gene.name")[,4])
-gene.set2 <- as.list(gene.set2)
-names(gene.set2) <- "Club"
-gene.set <- c(gene.set,gene.set2)
-gene.set2 <- as.data.frame(gene.set1$Ciliated[!is.na(gene.set1$Ciliated)])
-colnames(gene.set2) <- "genes"
-gene.set2 <- as.data.frame(merge(gene.set2,gene.orthog,by.x="genes",by.y="Gene.name")[,4])
-gene.set2 <- as.list(gene.set2)
-names(gene.set2) <- "Ciliated"
-gene.set <- c(gene.set,gene.set2)
-gene.set2 <- as.data.frame(gene.set1$Tuft[!is.na(gene.set1$Tuft)])
-colnames(gene.set2) <- "genes"
-gene.set2 <- as.data.frame(merge(gene.set2,gene.orthog,by.x="genes",by.y="Gene.name")[,4])
-gene.set2 <- as.list(gene.set2)
-names(gene.set2) <- "Tuft"
-gene.set <- c(gene.set,gene.set2)
-gene.set2 <- as.data.frame(gene.set1$Neuroendocrine[!is.na(gene.set1$Neuroendocrine)])
-colnames(gene.set2) <- "genes"
-gene.set2 <- as.data.frame(merge(gene.set2,gene.orthog,by.x="genes",by.y="Gene.name")[,4])
-gene.set2 <- as.list(gene.set2)
-names(gene.set2) <- "Neuroendocrine"
-gene.set <- c(gene.set,gene.set2)
-gene.set2 <- as.data.frame(gene.set1$Ionocyte[!is.na(gene.set1$Ionocyte)])
-colnames(gene.set2) <- "genes"
-gene.set2 <- as.data.frame(merge(gene.set2,gene.orthog,by.x="genes",by.y="Gene.name")[,4])
-gene.set2 <- as.list(gene.set2)
-names(gene.set2) <- "Ionocyte"
-gene.set <- c(gene.set,gene.set2)
-rm(gene.set2)
-gene.set1 <- read_csv("./genesets/SupTab6_Krt13_Hillock.csv",skip=6)
-gene.set1 <- gene.set1[gene.set1$FDR<=0.05 & gene.set1$'log2 fold-change (MAST)'>=1.5,1]
-colnames(gene.set1) <- "genes"
-gene.set1 <- as.data.frame(merge(gene.set1,gene.orthog,by.x="genes",by.y="Gene.name")[,4])
-gene.set1 <- as.list(gene.set1)
-names(gene.set1) <- "Hillock"
-gene.set <- c(gene.set,gene.set1)
-rm(gene.set1)
-rm(gene.orthog)
-gene.set <- lapply(gene.set,droplevels)
-results.cor.Epi.MusLungHierarchy <- scQuSAGEsm(sc10x.Epi,gs=gene.set,ds=min.epi,nm="Epi.dws.sub_NE",folder="MusLungHierarchy")
-rm(gene.set)
+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")
 
-gene.set.c2.all <- read.gmt("./genesets/c2.all.v6.1.symbols.gmt")
-results.cor.c2.ALL <- scQuSAGElg(sc10x,gs=gene.set.c2.all,ds=min.all,nm="Merge_Epi.dws_St.go",folder="c2.all")
-rm(gene.set.c2.all)
+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_St.go")
-sc10x@ident <- factor(sc10x@ident,levels=c("BE","LE","OE1","OE2","Fib","SM","Endo","Leu"))
+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","Unknown"))
 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))
@@ -403,125 +327,11 @@ 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_St.go")
-scTables(sc10x,i.1="samples",i.2="Merge_Epi.dws_St.go_NE")
-scTables(sc10x,i.1="Merge_Epi.dws_St.go_NE",i.2="Merge_Epi.dws_St.go")
-
-genes.deg.Stress <- scDEG(sc10x.Stress,i="Stress",g.1="Stress",g.2="ALL",pct=0.5,t=5)
-
-genes.deg.Epi <- scDEG(sc10x,i="Lin",g.1="Epi",g.2="St",t=2)
-genes.deg.St <- scDEG(sc10x,i="Lin",g.1="St",g.2="Epi",t=2)
-
-genes.deg.BE <- scDEG(sc10x.Epi.NE,i="Epi.dws.sub",g.1="BE",g.2=c("LE","OE1","OE2"),pct=0.25,t=2)
-genes.deg.LE <- scDEG(sc10x.Epi.NE,i="Epi.dws.sub",g.1="LE",g.2=c("BE","LE","OE1"),pct=0.25,t=2)
-genes.deg.OE1 <- scDEG(sc10x.Epi.NE,i="Epi.dws.sub",g.1="OE1",g.2=c("BE","LE","OE2"),pct=0.25,t=2)
-genes.deg.OE2 <- scDEG(sc10x.Epi.NE,i="Epi.dws.sub",g.1="OE2",g.2=c("BE","LE","OE1"),pct=0.25,t=2)
-
-genes.deg.NE <- scDEG(sc10x.Epi.NE,i="NE",g.1="NE",g.2="ALL",pct=0.01,t=1)
-
-genes.deg.Fib <- scDEG(sc10x.St,i="Merge_Epi.dws_St.go_NE",g.1="Fib",g.2=c("SM","Endo","Leu"),pct=0.25,t=2)
-genes.deg.SM <- scDEG(sc10x.St,i="Merge_Epi.dws_St.go_NE",g.1="SM",g.2=c("Fib","Endo","Leu"),pct=0.25,t=2)
-genes.deg.Endo <- scDEG(sc10x.St,i="Merge_Epi.dws_St.go_NE",g.1="Endo",g.2=c("Fib","SM","Leu"),pct=0.25,t=2)
-genes.deg.Leu <- scDEG(sc10x.St,i="St.go",g.1="Leu",g.2=c("Fib","SM","Endo"),pct=0.25,t=2)
-
-genes.deg.BE.unique <- setdiff(rownames(genes.deg.BE),Reduce(union,list(rownames(genes.deg.LE),rownames(genes.deg.OE1),rownames(genes.deg.OE2),rownames(genes.deg.NE),rownames(genes.deg.Fib),rownames(genes.deg.SM),rownames(genes.deg.Endo),rownames(genes.deg.Leu))))
-genes.deg.LE.unique <- setdiff(rownames(genes.deg.LE),Reduce(union,list(rownames(genes.deg.BE),rownames(genes.deg.OE1),rownames(genes.deg.OE2),rownames(genes.deg.NE),rownames(genes.deg.Fib),rownames(genes.deg.SM),rownames(genes.deg.Endo),rownames(genes.deg.Leu))))
-genes.deg.OE1.unique <- setdiff(rownames(genes.deg.OE1),Reduce(union,list(rownames(genes.deg.BE),rownames(genes.deg.LE),rownames(genes.deg.OE2),rownames(genes.deg.NE),rownames(genes.deg.Fib),rownames(genes.deg.SM),rownames(genes.deg.Endo),rownames(genes.deg.Leu))))
-genes.deg.OE2.unique <- setdiff(rownames(genes.deg.OE2),Reduce(union,list(rownames(genes.deg.BE),rownames(genes.deg.LE),rownames(genes.deg.OE1),rownames(genes.deg.NE),rownames(genes.deg.Fib),rownames(genes.deg.SM),rownames(genes.deg.Endo),rownames(genes.deg.Leu))))
-genes.deg.NE.unique <- setdiff(rownames(genes.deg.NE),Reduce(union,list(rownames(genes.deg.BE),rownames(genes.deg.LE),rownames(genes.deg.OE1),rownames(genes.deg.OE2),rownames(genes.deg.Fib),rownames(genes.deg.SM),rownames(genes.deg.Endo),rownames(genes.deg.Leu))))
-genes.deg.Fib.unique <- setdiff(rownames(genes.deg.Fib),Reduce(union,list(rownames(genes.deg.BE),rownames(genes.deg.LE),rownames(genes.deg.OE1),rownames(genes.deg.OE2),rownames(genes.deg.NE),rownames(genes.deg.SM),rownames(genes.deg.Endo),rownames(genes.deg.Leu))))
-genes.deg.SM.unique <- setdiff(rownames(genes.deg.SM),Reduce(union,list(rownames(genes.deg.BE),rownames(genes.deg.LE),rownames(genes.deg.OE1),rownames(genes.deg.OE2),rownames(genes.deg.NE),rownames(genes.deg.Fib),rownames(genes.deg.Endo),rownames(genes.deg.Leu))))
-genes.deg.Endo.unique <- setdiff(rownames(genes.deg.Endo),Reduce(union,list(rownames(genes.deg.BE),rownames(genes.deg.LE),rownames(genes.deg.OE1),rownames(genes.deg.OE2),rownames(genes.deg.NE),rownames(genes.deg.Fib),rownames(genes.deg.SM),rownames(genes.deg.Leu))))
-genes.deg.Leu.unique <- setdiff(rownames(genes.deg.Leu),Reduce(union,list(rownames(genes.deg.BE),rownames(genes.deg.LE),rownames(genes.deg.OE1),rownames(genes.deg.OE2),rownames(genes.deg.NE),rownames(genes.deg.Fib),rownames(genes.deg.SM),rownames(genes.deg.Endo))))
-
-genes.deg.5 <- c(genes.deg.BE.unique[1:5],genes.deg.LE.unique[1:5],genes.deg.OE1.unique[1:5],genes.deg.OE2.unique[1:5],genes.deg.NE.unique[1:5],genes.deg.Fib.unique[1:5],genes.deg.SM.unique[1:5],genes.deg.Endo.unique[1:5],genes.deg.Leu.unique[1:5])
-genes.deg.5 <- rev(genes.deg.5)
-genes.deg.10 <- c(genes.deg.BE.unique[1:10],genes.deg.LE.unique[1:10],genes.deg.OE1.unique[1:10],genes.deg.OE2.unique[1:10],genes.deg.NE.unique[1:10],genes.deg.Fib.unique[1:10],genes.deg.SM.unique[1:10],genes.deg.Endo.unique[1:10],genes.deg.Leu.unique[1:10])
-genes.deg.10 <- rev(genes.deg.10)
-
-sc10x <- SetAllIdent(object=sc10x,id="Merge_Epi.dws_St.go_NE")
-sc10x <- SetAllIdent(object=sc10x,id="Merge_Epi.dws_St.go_NE")
-sc10x@ident <- factor(sc10x@ident,levels=c("BE","LE","OE1","OE2","NE","Fib","SM","Endo","Leu"))
-postscript("./analysis/deg/Dot.DEG.eps",paper="special",width=10,height=5,horizontal=FALSE)
-DotPlot(sc10x,genes.deg.5,x.lab.rot=TRUE,plot.legend=TRUE,dot.scale=4)
-dev.off()
-postscript("./analysis/deg/Heatmap.DEG.eps",paper="special",width=10,height=5,horizontal=FALSE)
-plot(DoHeatmap(sc10x,genes.use=genes.deg.10,slim.col.label=TRUE,group.label.rot=TRUE,group.spacing=0.25,cex.row=2.5))
-dev.off()
-
-for (i in ls(pattern="^genes.deg*unique")){
-  postscript(paste0("./analysis/deg/Violin.",i,".eps"),paper="special",width=10,height=5,horizontal=FALSE)
-  plot(VlnPlot(sc10x,features.plot=get(i)[1:10],nCol=5,point.size.use=0.1,size.title.use=15,x.lab.rot=TRUE))
-  dev.off()
-  postscript(paste0("./analysis/deg/Ridge.",i,".eps"),paper="special",width=10,height=5,horizontal=FALSE)
-  plot(RidgePlot(sc10x,features.plot=get(i)[1:10],nCol=5,size.x.use=10,size.title.use=15))
-  dev.off()
-  postscript(paste0("./analysis/deg/Heatmap.",i,".eps"),paper="special",width=10,height=5,horizontal=FALSE)
-  plot(DoHeatmap(sc10x,genes.use=get(i),slim.col.label=TRUE,group.label.rot=TRUE,group.spacing=0.25,cex.row=2.5))
-  dev.off()
-}
-
-sc10x.Stress <- SetAllIdent(object=sc10x.Stress,id="Stress")
-postscript("./analysis/deg/Violin.Stress.eps",paper="special",width=10,height=5,horizontal=FALSE)
-plot(VlnPlot(sc10x.Stress,features.plot=rownames(genes.deg.Stress)[1:10],nCol=5,point.size.use=0.1,size.title.use=15,x.lab.rot=TRUE))
-dev.off()
-postscript("./analysis/deg/Ridge.Stress.eps",paper="special",width=10,height=5,horizontal=FALSE)
-plot(RidgePlot(sc10x.Stress,features.plot=rownames(genes.deg.Stress)[1:10],nCol=5,size.x.use=10,size.title.use=15))
-dev.off()
-postscript("./analysis/deg/Heatmap.Stress.eps",paper="special",width=10,height=5,horizontal=FALSE)
-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=",")
-}
-
-sctSNECustCol(sc10x,i="Lin",bl="Epi",rd="St",file="D")
-sctSNECustCol(sc10x,i="mLin",bl="Epi",rd="St",file="D")
-sctSNECustCol(sc10x,i="Merge_Epi.dws_St.go",bl=c("BE","LE","OE1","OE2"),rd=c("Fib","SM","Endo","Leu"),file="D")
-sctSNECustCol(sc10x.Epi,i="Epi.dws.sub",bl=c("BE","LE","OE1","OE2"),rd="",file="D")
-sctSNECustCol(sc10x.St,i="St.go",bl="",rd=c("Fib","SM","Endo","Leu"),file="D")
-
-sctSNEbwCol(sc10x,i="res0.2",file="ALL",files="D")
-sctSNEbwCol(sc10x.Epi,i="res0.2",file="Epi",files="D")
-sctSNEbwCol(sc10x.St,i="res0.2",file="St",files="D")
-sctSNEbwCol(sc10x,i="Merge_Epi.dws_St.go",file="ALL",files="D")
-sctSNEbwCol(sc10x.Epi,i="Epi.dws.sub",file="Epi",files="D")
-sctSNEbwCol(sc10x.St,i="St.go",file="St",files="D")
-
-for (g in c("Epi","St","Unknown")){
-  sctSNEHighlight(sc10x,i="Lin",g=g,file="D")
-}
-sctSNEHighlight(sc10x,i="mLin",g=c("St"),file="D")
-for (g in c("BE","LE","OE1","OE2")){
-  sctSNEHighlight(sc10x,i="Merge_Epi.dws_St.go",g=g,file="D")
-  sctSNEHighlight(sc10x.Epi,i="Epi.dws.sub",g=g,file="D")
-}
-sctSNEHighlight(sc10x.Epi.NE,i="NE",g="NE",file="D")
-for (g in c("Fib","SM","Endo","Leu")){
-  sctSNEHighlight(sc10x,i="Merge_Epi.dws_St.go",g=g,file="D")
-  sctSNEHighlight(sc10x.St,i="St.go",g=g,file="D")
-}
-for (g in c("D17","D27","D35")){
-  sctSNEHighlight(sc10x,i="patient",g=g,file="D")
-}
-sctSNEHighlight(sc10x,i="Pz",g="Pz",file="D")
-sctSNEHighlight(sc10x,i="Tz",g="Tz",file="D")
-rm(i)
-rm(g)
-
-scCustHeatmap(sc10x.Epi,i="Epi.dws.sub",gs=c(genes.deg.BE.unique,genes.deg.LE.unique,genes.deg.OE1.unique,genes.deg.OE2.unique),g=c("BE","LE","OE1","OE2"))
-scCustHeatmap(sc10x.St,i="St.go",gs=c(genes.deg.Fib.unique,genes.deg.SM.unique,genes.deg.Endo.unique,genes.deg.OE2.unique,genes.deg.Leu.unique),g=c("Fib","SM","Endo","Leu"))
+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")
 
-scPseudotime(sc10x.Epi,i="Epi.dws.sub",ds=0)
+sctSNECustCol(sc10x,i="Merge_Epi.dws.sc_St.dws.sc",bl=c("BE","LE","OE_SCGB","OE_KRT13"),rd=c("Fib","SM","Endo","Leu","Unknown"),file="D")
 
 save(list=ls(pattern="sc10x.Stress"),file="./analysis/sc10x.Stress.Rda")
 rm(list=ls(pattern="sc10x.Stress"))
-- 
GitLab