diff --git a/r.scripts/sc-TissueMapper_RUN.D17_FACS.R b/r.scripts/sc-TissueMapper_RUN.D17_FACS.R index f884af27dad9681b0f722d3ae5b8f7ee4e7e75f6..c991ebae0218f118363fc41d5031d2f2906a2836 100644 --- a/r.scripts/sc-TissueMapper_RUN.D17_FACS.R +++ b/r.scripts/sc-TissueMapper_RUN.D17_FACS.R @@ -107,7 +107,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=0.2,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)) @@ -143,19 +143,6 @@ if (opt$cc==TRUE){ } gc() -#sc10x.D17 <- scSubset(sc10x,"D17","D17") -#sc10x.D27 <- scSubset(sc10x,"D27","D27") -#sc10x.D35 <- scSubset(sc10x,"D35","D35") - -#results <- sc3CCA(sc10x.D17,sc10x.D27,sc10x.D35,"D17","D27","D35",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.D17) -#rm(sc10x.D27) -#rm(sc10x.D35) - results <- scPC(sc10x,lx=opt$lx,hx=opt$hx,ly=opt$ly,cc=opt$cc,pc=50,hpc=0.85,file="pre.stress",cca=FALSE) sc10x <- results[[1]] genes.hvg.prestress <- results[[2]] @@ -192,8 +179,8 @@ names(gene.set1) <- "St" gene.set <- c(gene.set,gene.set1) rm(gene.set1) gc() -min.all <- min(table(sc10x@meta.data[,paste0("res",0.5)])) -results <- scQuSAGE(sc10x,gs=gene.set,res.use=0.5,ds=min.all,nm="Lin",folder="lin") +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="Lin",folder="lin") sc10x <- results[[1]] results.cor.Lin <- results[[2]] results.clust.Lin.id <- results[[3]] @@ -206,12 +193,12 @@ if (any(levels(sc10x@ident)=="Unknown")){ } else { sc10x.St <- scSubset(sc10x,i="Lin",g="St") } -sc10x.Epi <- SetAllIdent(object=sc10x.Epi,id=paste0("res",0.5)) +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",0.5)) -sc10x.St <- SetAllIdent(object=sc10x.St,id=paste0("res",0.5)) +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",0.5)) +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")) @@ -220,7 +207,7 @@ plot <- plot+theme(axis.text.x=element_text(size=20),axis.text.y=element_text(si plot <- plot+guides(colour=guide_legend(override.aes=list(size=10))) plot(plot) dev.off() -postscript(paste0("./analysis/tSNE/epi/tSNE_res",0.5,".eps")) +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))) @@ -248,78 +235,16 @@ gene.set1 <- gene.set1[1] gene.set1 <- as.list(gene.set1) names(gene.set1) <- "OE_KRT13" gene.set <- c(gene.set,gene.set1) -#gene.set1 <- read_delim("./genesets/genes.deg.NE.csv",",",escape_double=FALSE,trim_ws=TRUE,col_names=TRUE) -#gene.set1 <- gene.set1[1] -#gene.set1 <- as.list(gene.set1) -#names(gene.set1) <- "NE" -#gene.set <- c(gene.set,gene.set1) rm(gene.set1) gc() -min.epi <- min(table(sc10x.Epi@meta.data[,paste0("res",0.5)])) -results <- scQuSAGE(sc10x.Epi,gs=gene.set,res.use=0.5,ds=min.epi,nm="Epi.dws.sc",folder="epi") +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.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="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",0.5,".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.Endo.csv",",",escape_double=FALSE,trim_ws=TRUE,col_names=TRUE) gene.set1 <- gene.set1[1] gene.set1 <- as.list(gene.set1) @@ -342,8 +267,8 @@ 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",0.5)])) -results <- scQuSAGE(sc10x.St,gs=gene.set,res.use=0.5,ds=min.st,nm="St.dws.sc",folder="st") +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=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]] @@ -352,72 +277,12 @@ rm(gene.set) 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") - -#sc10x.St <- scMergeSubClust(sc10x.St,i="St.go",g=c("Endo","SM","Fib","Leu"),nm="Merge") - 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") -#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) - -#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 <- 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") @@ -431,91 +296,14 @@ 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") -# 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="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="res0.5",file="ALL",files="D17") -sctSNEbwCol(sc10x.Epi,i="res0.5",file="Epi",files="D17") -sctSNEbwCol(sc10x.St,i="res0.5",file="St",files="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") @@ -535,24 +323,15 @@ for (g in c("Fib","SM","Endo","Leu")){ for (g in c("D17PrPzF_BE","D17PrPzF_LE","D17PrPzF_OE","D17PrPzF_FMSt")){ sctSNEHighlight(sc10x,i="samples",g=g,file="D17") } -#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")) - 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() -#scPseudotime(sc10x.Epi,i="Epi.dws.sub",ds=0) - save(list=ls(pattern="sc10x.Stress"),file="./analysis/sc10x.Stress.Rda") rm(list=ls(pattern="sc10x.Stress")) -save(list=ls(pattern="^genes.deg"),file="./analysis/DEG.Rda") -rm(list=ls(pattern="^genes.deg")) 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")