Add code for DS.D17

r.scripts/sc-TissueMapper_RUN.DS_D17.R

+gc()
+library(methods)
+library(optparse)
+library(Seurat)
+library(readr)
+library(fBasics)
+library(pastecs)
+library(qusage)
+library(RColorBrewer)
+library(monocle)
+library(dplyr)
+library(viridis)
+library(reshape2)
+library(NMI)
+
+source("../r.scripts/sc-TissueMapper.R")
+
+#Create folder structure
+setwd("../")
+if (!dir.exists("./analysis")){
+  dir.create("./analysis")
+}
+if (!dir.exists("./analysis/qc")){
+  dir.create("./analysis/qc")
+}
+if (!dir.exists("./analysis/qc/cc")){
+  dir.create("./analysis/qc/cc")
+}
+if (!dir.exists("./analysis/tSNE")){
+  dir.create("./analysis/tSNE")
+}
+if (!dir.exists("./analysis/tSNE/pre.stress")){
+  dir.create("./analysis/tSNE/pre.stress")
+}
+if (!dir.exists("./analysis/pca")){
+  dir.create("./analysis/pca")
+}
+if (!dir.exists("./analysis/pca/stress")){
+  dir.create("./analysis/pca/stress")
+}
+if (!dir.exists("./analysis/violin")){
+  dir.create("./analysis/violin")
+}
+if (!dir.exists("./analysis/violin/stress")){
+  dir.create("./analysis/violin/stress")
+}
+if (!dir.exists("./analysis/table")){
+  dir.create("./analysis/table")
+}
+if (!dir.exists("./analysis/tSNE/post.stress")){
+  dir.create("./analysis/tSNE/post.stress")
+}
+if (!dir.exists("./analysis/cor")){
+  dir.create("./analysis/cor")
+}
+if (!dir.exists("./analysis/tSNE/lin")){
+  dir.create("./analysis/tSNE/lin")
+}
+if (!dir.exists("./analysis/tSNE/epi")){
+  dir.create("./analysis/tSNE/epi")
+}
+if (!dir.exists("./analysis/tSNE/st")){
+  dir.create("./analysis/tSNE/st")
+}
+if (!dir.exists("./analysis/tSNE/merge")){
+  dir.create("./analysis/tSNE/merge")
+}
+if (!dir.exists("./analysis/pca/ne")){
+  dir.create("./analysis/pca/ne")
+}
+if (!dir.exists("./analysis/tSNE/ne")){
+  dir.create("./analysis/tSNE/ne")
+}
+if (!dir.exists("./analysis/violin/ne")){
+  dir.create("./analysis/violin/ne")
+}
+if (!dir.exists("./analysis/tSNE/FINAL")){
+  dir.create("./analysis/tSNE/FINAL")
+}
+if (!dir.exists("./analysis/deg")){
+  dir.create("./analysis/deg")
+}
+if (!dir.exists("./analysis/cca")){
+  dir.create("./analysis/cca")
+}
+if (!dir.exists("./analysis/diy")){
+  dir.create("./analysis/diy")
+}
+if (!dir.exists("./analysis/pseudotime")){
+  dir.create("./analysis/pseudotime")
+}
+
+#Retrieve command-line options
+option_list=list(
+  make_option("--p",action="store",default="DPrF",type='character',help="Project Name"),
+  make_option("--g",action="store",default="ALL",type='character',help="Group To analyze"),
+  make_option("--lg",action="store",default=500,type='integer',help="Threshold for cells with minimum genes"),
+  make_option("--hg",action="store",default=3000,type='integer',help="Threshold for cells with maximum genes"),
+  make_option("--lm",action="store",default=0,type='numeric',help="Threshold for cells with minimum %mito genes"),
+  make_option("--hm",action="store",default=0.1,type='numeric',help="Threshold for cells with maximum %mito genes"),
+  make_option("--lx",action="store",default=0.2,type='numeric',help="x low threshold for hvg selection"),
+  make_option("--hx",action="store",default=5,type='numeric',help="x high threshold for hvg selection"),
+  make_option("--ly",action="store",default=1,type='numeric',help="y low threshold for hvg selection"),
+  make_option("--cc",action="store",default=TRUE,type='logical',help="Scale cell cycle?"),
+  make_option("--cca",action="store",default=50,type='integer',help="Number of CCAs to cacluate"),
+  make_option("--acca",action="store",default=30,type='integer',help="Number of CCAs to align"),
+  make_option("--pc",action="store",default=50,type='integer',help="Number of PCs to cacluate"),
+  make_option("--res.prestress",action="store",default=1,type='numeric',help="Resolution to cluster, pre-stress"),
+  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.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))
+rm(option_list)
+if (opt$lm==0){opt$lm=-Inf}
+
+sc10x.data <- Read10X(data.dir="./analysis/DATA/10x/filtered_gene_bc_matrices/GRCh38/")
+sc10x <- new("seurat",raw.data=sc10x.data)
+cell.codes <- as.data.frame(sc10x@raw.data@Dimnames[[2]])
+colnames(cell.codes) <- "barcodes"
+rownames(cell.codes) <- cell.codes$barcodes
+cell.codes$samples <- "All"
+sc10x <- CreateSeuratObject(raw.data=sc10x.data,meta.data=cell.codes["samples"],min.cells=3,min.genes=200,project="DS.D17")
+#sc10x <- AddMetaData(sc10x,cell.codes["samples"],"samples")
+rm(cell.codes)
+rm(sc10x.data)
+
+if (opt$cc==TRUE){
+  results <- scCellCycle(sc10x)
+  sc10x <- results[[1]]
+  genes.s <- results[[2]]
+  genes.g2m <- results[[3]]
+  rm(results)
+} else {
+  genes.s=""
+  genes.g2m=""
+}
+
+results <- scQC(sc10x,lg=opt$lg,hg=opt$hg,lm=opt$lm,hm=opt$hm)
+sc10x <- results[[1]]
+counts.cell.raw <- results[[2]]
+counts.gene.raw <- results[[3]]
+counts.cell.filtered <- results[[4]]
+counts.gene.filtered <- results[[5]]
+rm(results)
+
+gc()
+if (opt$cc==TRUE){
+  sc10x <- ScaleData(object=sc10x,vars.to.regress=c("nUMI","percent.mito","S.Score","G2M.Score"),display.progress=FALSE,do.par=TRUE,num.cores=45)
+} else {
+  sc10x <- ScaleData(object=sc10x,vars.to.regress=c("nUMI","percent.mito"),display.progress=FALSE,do.par=TRUE,num.cores=45)
+}
+gc()
+
+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]]
+pc.use.prestress <- results[[3]]
+rm(results)
+
+sc10x <- scCluster(sc10x,pc.use=pc.use.prestress,res.use=opt$res.prestress,folder="pre.stress",red="pca")
+
+if (opt$st==TRUE){
+  results <- scStress(sc10x,stg=opt$stg,res.use=opt$res.prestress,cut=opt$cut.stress)
+  sc10x <- results[[1]]
+  counts.cell.filtered.stress <- results[[2]]
+  sc10x.Stress <- results[[3]]
+  rm(results)
+  
+  results <- scPC(sc10x,lx=opt$lx,hx=opt$hx,ly=opt$ly,cc=opt$cc,pc=50,hpc=0.85,file="post.stress",cca=FALSE)
+  sc10x <- results[[1]]
+  genes.hvg.poststress <- results[[2]]
+  pc.use.poststress <- results[[3]]
+  rm(results)
+  
+  sc10x <- scCluster(sc10x,pc.use=pc.use.poststress,res.use=0.5,folder="post.stress",red="pca")
+}
+
+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="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.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)
+
+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 <- 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")
+}
+rm(i)
+rm(g)
+
+
+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")