diff --git a/genesets/40425_2017_215_MOESM1_ESM.xlsx b/genesets/40425_2017_215_MOESM1_ESM.xlsx new file mode 100755 index 0000000000000000000000000000000000000000..5a813a3b5b13a22e578c267e8269b26ab51aca94 Binary files /dev/null and b/genesets/40425_2017_215_MOESM1_ESM.xlsx differ diff --git a/r.scripts/sc-TissueMapper_RUN.VAMC015PrFx.R b/r.scripts/sc-TissueMapper_RUN.VAMC015PrFx.R deleted file mode 100644 index c4d16dc38d19bd00c2a9e25aafdaf014012afc85..0000000000000000000000000000000000000000 --- a/r.scripts/sc-TissueMapper_RUN.VAMC015PrFx.R +++ /dev/null @@ -1,391 +0,0 @@ -gc() -library(methods) -library(optparse) -library(Seurat) -library(readr) -library(fBasics) -library(pastecs) -library(qusage) -library(RColorBrewer) -library(monocle) -library(dplyr) -library(viridis) - -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=1,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 <- scLoad("VAMC015PrFx") - -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=opt$res.poststress,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=0,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=0,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=0,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") -} -for (g in c("D17PrPzF_BE","D17PrPzF_LE","D17PrPzF_OE","D17PrPzF_FMSt")){ - sctSNEHighlight(sc10x,i="samples",g=g,file="D17") -} -rm(i) -rm(g) - -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() - -postscript(paste0("./analysis/diy/Feature.KLK3.eps")) -VlnPlot(sc10x,features.plot="KLK3",group.by="Lin",size.title.use=20,point.size.use=0.05) -dev.off() - -sc10x.LE <- scSubset(sc10x,i="Merge_Epi.dws.sc_St.dws.sc",g="LE") -genes.deg.LE.Enza.vs.Ctrl <- scDEG(sc10x.LE,i="samples",g.1="VAMC015PrFx_Enza_Via",g.2="VAMC015PrFx_Ctrl_Via",pct=0.25,t=0) -genes.deg.LE.ODM.vs.Ctrl <- scDEG(sc10x.LE,i="samples",g.1="VAMC015PrFx_ODM_Via",g.2="VAMC015PrFx_Ctrl_Via",pct=0.25,t=0) -genes.deg.LE.ODM.vs.Enza <- scDEG(sc10x.LE,i="samples",g.1="VAMC015PrFx_ODM_Via",g.2="VAMC015PrFx_Enza_Via",pct=0.25,t=0) - -genes.deg.LE.Ctrl.vs.Enza <- scDEG(sc10x.LE,i="samples",g.1="VAMC015PrFx_Ctrl_Via",g.2="VAMC015PrFx_Enza_Via",pct=0.25,t=0) -genes.deg.LE.Ctrl.vs.ODM <- scDEG(sc10x.LE,i="samples",g.1="VAMC015PrFx_Ctrl_Via",g.2="VAMC015PrFx_ODM_Via",pct=0.25,t=0) -genes.deg.LE.Enza.vs.ODM <- scDEG(sc10x.LE,i="samples",g.1="VAMC015PrFx_Enza_Via",g.2="VAMC015PrFx_ODM_Via",pct=0.25,t=0) - -postscript(paste0("./analysis/diy/Feature.LE.KLK3.eps")) -VlnPlot(sc10x.LE,features.plot="KLK3",group.by="samples",size.title.use=20,point.size.use=0.05) -dev.off() - -sc10x.Epi <- scSubset(sc10x,i="Lin",g="Epi") -genes.deg.Epi.Enza.vs.Ctrl <- scDEG(sc10x.Epi,i="samples",g.1="VAMC015PrFx_Enza_Via",g.2="VAMC015PrFx_Ctrl_Via",pct=0.25,t=0) -genes.deg.Epi.ODM.vs.Ctrl <- scDEG(sc10x.Epi,i="samples",g.1="VAMC015PrFx_ODM_Via",g.2="VAMC015PrFx_Ctrl_Via",pct=0.25,t=0) -genes.deg.Epi.ODM.vs.Enza <- scDEG(sc10x.Epi,i="samples",g.1="VAMC015PrFx_ODM_Via",g.2="VAMC015PrFx_Enza_Via",pct=0.25,t=0) - -genes.deg.Epi.Ctrl.vs.Enza <- scDEG(sc10x.Epi,i="samples",g.1="VAMC015PrFx_Ctrl_Via",g.2="VAMC015PrFx_Enza_Via",pct=0.25,t=0) -genes.deg.Epi.Ctrl.vs.ODM <- scDEG(sc10x.Epi,i="samples",g.1="VAMC015PrFx_Ctrl_Via",g.2="VAMC015PrFx_ODM_Via",pct=0.25,t=0) -genes.deg.Epi.Enza.vs.ODM <- scDEG(sc10x.Epi,i="samples",g.1="VAMC015PrFx_Enza_Via",g.2="VAMC015PrFx_ODM_Via",pct=0.25,t=0) - -postscript(paste0("./analysis/diy/Feature.Epi.KLK3.eps")) -VlnPlot(sc10x.Epi,features.plot="KLK3",group.by="samples",size.title.use=20,point.size.use=0.05) -dev.off() - -for (i in c("genes.deg.LE.Ctrl.vs.Enza","genes.deg.LE.Ctrl.vs.ODM","genes.deg.LE.ODM.vs.Enza","genes.deg.LE.Enza.vs.ODM")){ -write.table(get(i),file=paste0(i,".csv"),row.names=TRUE,col.names=NA,append=FALSE,sep=",") -} - -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") diff --git a/r.scripts/sc-TissueMapper_RUN_PrFx_.VAMC015.R b/r.scripts/sc-TissueMapper_RUN_PrFx_.VAMC015.R new file mode 100644 index 0000000000000000000000000000000000000000..2970d345e5a9c7dc4b3180d4190600dc2b9b35fd --- /dev/null +++ b/r.scripts/sc-TissueMapper_RUN_PrFx_.VAMC015.R @@ -0,0 +1,235 @@ +gc() +library(methods) +library(optparse) +library(Seurat) +library(readr) +library(fBasics) +library(pastecs) +library(qusage) +library(RColorBrewer) +library(monocle) +library(dplyr) +library(viridis) +library(readxl) + +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 <- scLoad("VAMC015PrFx") + +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=opt$res.poststress,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) +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.set,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) +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.set,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) +genes.leu <- read_excel("genesets/40425_2017_215_MOESM1_ESM.xlsx",sheet="S3. Candidate markers") +leu <- as.list(unique(genes.leu[,2]))$Cell +leu.l <- leu[-c(1,3,4,7:9,14,15,17:18,20:21,23:30)] +genes.leu <- genes.leu[unlist(genes.leu[,2]) %in% unlist(leu),] +genes.leu.l <- genes.leu[unlist(genes.leu[,2]) %in% unlist(leu.l),] +gene.set1 <- split(genes.leu.l[,1], genes.leu.l[,2]) +gene.set1 <- lapply(gene.set1,unname) +gene.set1 <- lapply(gene.set1,unlist) +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="Pop",folder="lin") +sc10x <- results[[1]] +results.cor.Lin <- results[[2]] +results.clust.Lin.id <- results[[3]] +rm(results) +rm(gene.set) + +sc10x.Epi.NE <- scNE(sc10x.Epi,neg="dws",cut=opt$cut.ne) \ No newline at end of file