From 84c4e944ba34806811fcdef31d9378bc7bcb51bc Mon Sep 17 00:00:00 2001
From: "Gervaise H. Henry" <gervaise.henry@utsouthwestern.edu>
Date: Fri, 26 Oct 2018 09:40:21 -0500
Subject: [PATCH] Finalize (?) all code for PdPbPc to run

---
 bash.scripts/sc_TissueMapper-PdPbPc.sh    |  17 ++
 r.scripts/sc-TissueMapper_RUN.PdPbPc.Pb.R | 194 ++--------------------
 r.scripts/sc-TissueMapper_RUN.PdPbPc.Pc.R | 178 +-------------------
 r.scripts/sc-TissueMapper_RUN.PdPbPc.Pd.R | 194 ++--------------------
 r.scripts/sc-TissueMapper_RUN.PdPbPc.R    |   4 +-
 5 files changed, 43 insertions(+), 544 deletions(-)
 create mode 100644 bash.scripts/sc_TissueMapper-PdPbPc.sh

diff --git a/bash.scripts/sc_TissueMapper-PdPbPc.sh b/bash.scripts/sc_TissueMapper-PdPbPc.sh
new file mode 100644
index 0000000..7422a5c
--- /dev/null
+++ b/bash.scripts/sc_TissueMapper-PdPbPc.sh
@@ -0,0 +1,17 @@
+#!/bin/bash
+#SBATCH --job-name R_FullAnalysis
+#SBATCH -p 256GB,256GBv1,384GB
+#SBATCH -N 1
+#SBATCH -t 7-0:0:0
+#SBATCH -o job_%j.out
+#SBATCH -e job_%j.out
+#SBATCH --mail-type ALL
+#SBATCH --mail-user gervaise.henry@utsouthwestern.edu
+
+module load R/3.4.1-gccmkl
+
+Rscript ../r.scripts/sc-TissueMapper_RUN.PdPbPc.Pd.R &
+Rscript ../r.scripts/sc-TissueMapper_RUN.PdPbPc.Pb.R &
+Rscript ../r.scripts/sc-TissueMapper_RUN.PdPbPc.Pc.R &
+wait
+Rscript ../r.scripts/sc-TissueMapper_RUN.PdPbPc.R
diff --git a/r.scripts/sc-TissueMapper_RUN.PdPbPc.Pb.R b/r.scripts/sc-TissueMapper_RUN.PdPbPc.Pb.R
index e42948f..4a79f53 100644
--- a/r.scripts/sc-TissueMapper_RUN.PdPbPc.Pb.R
+++ b/r.scripts/sc-TissueMapper_RUN.PdPbPc.Pb.R
@@ -114,7 +114,7 @@ opt=parse_args(OptionParser(option_list=option_list))
 rm(option_list)
 if (opt$lm==0){opt$lm=-Inf}
 
-sc10x <- scLoad("Pb")
+sc10x <- scLoad("Pb",sub="Pb)
 
 sc10x <- scSubset(sc10x,i="Glandular",g="Glandular")
 
@@ -145,20 +145,18 @@ if (opt$cc==TRUE){
 }
 gc()
 
-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")
+sc10x.B327 <- scSubset(sc10x,"samples","BPH327PrGF_Via")
+sc10x.B340 <- scSubset(sc10x,"samples","BPH340PrSF_Via")
+sc10x.B342 <- scSubset(sc10x,"samples","BPH342PrF_Via")
 
-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)
+results <- sc3CCA(sc10x.B327,sc10x.B340,sc10x.B342,"B327","B340","B342",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.C11)
-rm(sc10x.C12)
-rm(sc10x.C13)
-rm(sc10x.C14)
+rm(sc10x.B327)
+rm(sc10x.B340)
+rm(sc10x.B342)
 
 sc10x <- scCluster(sc10x,pc.use=opt$acca,res.use=opt$res.prestress,folder="pre.stress",red="cca.aligned")
 
@@ -168,179 +166,7 @@ if (opt$st==TRUE){
   counts.cell.filtered.stress <- results[[2]]
   sc10x.Stress <- results[[3]]
   rm(results)
-  
-  sc10x <- scCluster(sc10x,pc.use=opt$acca,res.use=opt$res.poststress,folder="post.stress",red="cca.aligned")
 }
 
-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/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)
-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=min.all,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="cca.aligned",dims.use=1:opt$acca,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)
-
-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/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/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_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.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.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)
-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.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/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/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/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=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="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="Merge",i.2="Merge",nm="Merge_Epi.dws_St.go")
-
-sc10x <- scMerge(sc10x,sc10x,sc10x.Epi.NE,i.1="Merge_Epi.dws_St.go",i.2="NE",nm="Merge_Epi.dws_St.go_NE")
-
-sc10x.Epi <- scMerge(sc10x.Epi,sc10x.Epi,sc10x.Epi.NE,i.1="Epi.dws.sub",i.2="NE",nm="Epi.dws.sub_NE")
-
-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")
-
-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","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))
-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="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"))
-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")
-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.Pb <- sc10x
+save(sc10x.Pb,file="./analysis/sc10x.Pb.Rda")
diff --git a/r.scripts/sc-TissueMapper_RUN.PdPbPc.Pc.R b/r.scripts/sc-TissueMapper_RUN.PdPbPc.Pc.R
index f277894..563db10 100644
--- a/r.scripts/sc-TissueMapper_RUN.PdPbPc.Pc.R
+++ b/r.scripts/sc-TissueMapper_RUN.PdPbPc.Pc.R
@@ -114,7 +114,7 @@ opt=parse_args(OptionParser(option_list=option_list))
 rm(option_list)
 if (opt$lm==0){opt$lm=-Inf}
 
-sc10x <- scLoad("Pc")
+sc10x <- scLoad("Pc",sub="Pc")
 
 sc10x <- scSubset(sc10x,i="Tumor",g="Tumor")
 
@@ -168,179 +168,7 @@ if (opt$st==TRUE){
   counts.cell.filtered.stress <- results[[2]]
   sc10x.Stress <- results[[3]]
   rm(results)
-  
-  sc10x <- scCluster(sc10x,pc.use=opt$acca,res.use=opt$res.poststress,folder="post.stress",red="cca.aligned")
 }
 
-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/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)
-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=min.all,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="cca.aligned",dims.use=1:opt$acca,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)
-
-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/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/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_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.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.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)
-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.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/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/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/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=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="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="Merge",i.2="Merge",nm="Merge_Epi.dws_St.go")
-
-sc10x <- scMerge(sc10x,sc10x,sc10x.Epi.NE,i.1="Merge_Epi.dws_St.go",i.2="NE",nm="Merge_Epi.dws_St.go_NE")
-
-sc10x.Epi <- scMerge(sc10x.Epi,sc10x.Epi,sc10x.Epi.NE,i.1="Epi.dws.sub",i.2="NE",nm="Epi.dws.sub_NE")
-
-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")
-
-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","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))
-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="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"))
-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")
-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.Pc <- sc10x
+save(sc10x,file="./analysis/sc10x.Pc.Rda")
diff --git a/r.scripts/sc-TissueMapper_RUN.PdPbPc.Pd.R b/r.scripts/sc-TissueMapper_RUN.PdPbPc.Pd.R
index dddb26a..0aef238 100644
--- a/r.scripts/sc-TissueMapper_RUN.PdPbPc.Pd.R
+++ b/r.scripts/sc-TissueMapper_RUN.PdPbPc.Pd.R
@@ -114,7 +114,7 @@ opt=parse_args(OptionParser(option_list=option_list))
 rm(option_list)
 if (opt$lm==0){opt$lm=-Inf}
 
-sc10x <- scLoad("Pd")
+sc10x <- scLoad("Pd",sub="Pd")
 
 sc10x <- scSubset(sc10x,i="ALL",g="ALL")
 
@@ -145,20 +145,18 @@ if (opt$cc==TRUE){
 }
 gc()
 
-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")
+sc10x.D17 <- scSubset(sc10x,"patient","D17")
+sc10x.D27 <- scSubset(sc10x,"patient","D27")
+sc10x.D35 <- scSubset(sc10x,"patient","D35")
 
-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)
+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.C11)
-rm(sc10x.C12)
-rm(sc10x.C13)
-rm(sc10x.C14)
+rm(sc10x.D17)
+rm(sc10x.D27)
+rm(sc10x.D35)
 
 sc10x <- scCluster(sc10x,pc.use=opt$acca,res.use=opt$res.prestress,folder="pre.stress",red="cca.aligned")
 
@@ -168,179 +166,7 @@ if (opt$st==TRUE){
   counts.cell.filtered.stress <- results[[2]]
   sc10x.Stress <- results[[3]]
   rm(results)
-  
-  sc10x <- scCluster(sc10x,pc.use=opt$acca,res.use=opt$res.poststress,folder="post.stress",red="cca.aligned")
 }
 
-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/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)
-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=min.all,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="cca.aligned",dims.use=1:opt$acca,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)
-
-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/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/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_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.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.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)
-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.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/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/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/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=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="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="Merge",i.2="Merge",nm="Merge_Epi.dws_St.go")
-
-sc10x <- scMerge(sc10x,sc10x,sc10x.Epi.NE,i.1="Merge_Epi.dws_St.go",i.2="NE",nm="Merge_Epi.dws_St.go_NE")
-
-sc10x.Epi <- scMerge(sc10x.Epi,sc10x.Epi,sc10x.Epi.NE,i.1="Epi.dws.sub",i.2="NE",nm="Epi.dws.sub_NE")
-
-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")
-
-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","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))
-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="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"))
-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")
-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.Pd <- sc10x
+save(sc10x,file="./analysis/sc10x.Pd.Rda")
diff --git a/r.scripts/sc-TissueMapper_RUN.PdPbPc.R b/r.scripts/sc-TissueMapper_RUN.PdPbPc.R
index 43fabef..ab650ee 100755
--- a/r.scripts/sc-TissueMapper_RUN.PdPbPc.R
+++ b/r.scripts/sc-TissueMapper_RUN.PdPbPc.R
@@ -118,7 +118,9 @@ opt=parse_args(OptionParser(option_list=option_list))
 rm(option_list)
 if (opt$lm==0){opt$lm=-Inf}
 
-load("./analysis/sc10x.Rda")
+load("./analysis/sc10x.Pd.Rda")
+load("./analysis/sc10x.Pb.Rda")
+load("./analysis/sc10x.Pc.Rda")
 
 for (i in c("Pb","Pc","Pd")){
   if (!dir.exists(paste0("./analysis/tSNE/post.stress/",i))){
-- 
GitLab