From bc725ac9749ac2f3a7738aa60ba3f89e52347a04 Mon Sep 17 00:00:00 2001 From: "Gervaise H. Henry" <gervaise.henry@utsouthwestern.edu> Date: Sat, 10 Feb 2018 22:52:54 -0600 Subject: [PATCH] Update folder references --- r.scripts/sc_Cluster.R | 24 ++++++------ r.scripts/sc_D-SampleReorder.R | 10 ++--- r.scripts/sc_DEG.R | 38 +++++++++--------- r.scripts/sc_Demultiplex.R | 16 ++++---- r.scripts/sc_PC.Score.NE.R | 46 +++++++++++----------- r.scripts/sc_PC.Score.Stress.R | 46 +++++++++++----------- r.scripts/sc_QC.R | 56 +++++++++++++-------------- r.scripts/sc_QuSAGE.Epi.R | 40 +++++++++---------- r.scripts/sc_QuSAGE.LGEA.Epi.R | 40 +++++++++---------- r.scripts/sc_QuSAGE.LGEA.St.R | 28 +++++++------- r.scripts/sc_QuSAGE.Lineage.R | 36 ++++++++--------- r.scripts/sc_QuSAGE.St.R | 32 +++++++-------- r.scripts/sc_Seurat.Score.CellCycle.R | 26 ++++++------- r.scripts/sc_Tables.R | 20 +++++----- 14 files changed, 229 insertions(+), 229 deletions(-) diff --git a/r.scripts/sc_Cluster.R b/r.scripts/sc_Cluster.R index 7eb7865..b5d2ff8 100755 --- a/r.scripts/sc_Cluster.R +++ b/r.scripts/sc_Cluster.R @@ -16,23 +16,23 @@ Project.Name <- opt$p PC.Max <- opt$pc #load data -setwd(paste0("./",Project.Name)) -load(paste0("./sc10x.",Project.Name,".Cumulative.RData")) -load(paste0("./sc10x.",Project.Name,".",opt$g,".qc.Rda")) +setwd("../") +load(paste0("./analysis/sc10x.",Project.Name,".Cumulative.RData")) +load(paste0("./analysis/sc10x.",Project.Name,".",opt$g,".qc.Rda")) sc10x.Group <- get(paste0("sc10x.",opt$g,".qc")) rm(list=paste0("sc10x.",opt$g,".qc")) #create folder structure -if (!dir.exists(paste0("./",opt$g))){ - dir.create(paste0("./",opt$g)) +if (!dir.exists(paste0("./analysis/",opt$g))){ + dir.create(paste0("./analysis/",opt$g)) } -if (!dir.exists(paste0("./",opt$g,"/global"))){ - dir.create(paste0("./",opt$g,"/global")) +if (!dir.exists(paste0("./analysis/",opt$g,"/global"))){ + dir.create(paste0("./analysis/",opt$g,"/global")) } #create tSNE of Sample sc10x.Group <- RunTSNE(object=sc10x.Group,dims.use=1:PC.Max,do.fast=TRUE) -png(paste0("./",opt$g,"/global/tSNE_Sample.png"),1000,1000,pointsize=20) +png(paste0("./analysis/",opt$g,"/global/tSNE_Sample.png"),1000,1000,pointsize=20) plot <- TSNEPlot(object=sc10x.Group,group.by="samples",pt.size=2.5,do.return=TRUE) 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))) @@ -44,10 +44,10 @@ for (i in c(5,2,1,0.5,0.2,0.1)){ sc10x.Group <- FindClusters(object=sc10x.Group,reduction.type="pca",dims.use=1:PC.Max,resolution=i,print.output=0,save.SNN=TRUE) sc10x.Group <- BuildClusterTree(sc10x.Group,do.reorder=TRUE,reorder.numeric=TRUE,do.plot=FALSE) sc10x.Group <- StashIdent(object=sc10x.Group,save.name=paste0("res",i)) - png(paste0("./",opt$g,"/global/Tree_Global.res",i,".png"),1000,1000,pointsize=20) + png(paste0("./analysis/",opt$g,"/global/Tree_Global.res",i,".png"),1000,1000,pointsize=20) PlotClusterTree(object=sc10x.Group) dev.off() - png(paste0("./",opt$g,"/global/tSNE_Global.res",i,".png"),1000,1000,pointsize=20) + png(paste0("./analysis/",opt$g,"/global/tSNE_Global.res",i,".png"),1000,1000,pointsize=20) plot <- TSNEPlot(object=sc10x.Group,pt.size=5,do.label=TRUE,label.size=10,do.return=TRUE) 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))) @@ -67,6 +67,6 @@ opt.cluster <- opt rm(opt) #save cluser.Rda file,Cumulative.RData file -save(list=paste0("sc10x.",opt.cluster$g,".cluster"),file=paste0("./sc10x.",Project.Name,".",opt.cluster$g,".cluster.Rda")) +save(list=paste0("sc10x.",opt.cluster$g,".cluster"),file=paste0("./analysis/sc10x.",Project.Name,".",opt.cluster$g,".cluster.Rda")) rm(list=paste0("sc10x.",opt.cluster$g,".cluster")) -save.image(file=paste0("./sc10x.",Project.Name,".Cumulative.RData")) +save.image(file=paste0("./analysis/sc10x.",Project.Name,".Cumulative.RData")) diff --git a/r.scripts/sc_D-SampleReorder.R b/r.scripts/sc_D-SampleReorder.R index c7b0b12..6d4ed1a 100755 --- a/r.scripts/sc_D-SampleReorder.R +++ b/r.scripts/sc_D-SampleReorder.R @@ -13,16 +13,16 @@ rm(option_list) Project.Name <- opt$p #load data -setwd(paste0("./",Project.Name)) -load(paste0("./sc10x.",Project.Name,".",opt$g,".raw.Rda")) +setwd("../") +load(paste0("./analysis/sc10x.",Project.Name,".",opt$g,".raw.Rda")) sc10x.Group <- get(paste0("sc10x.",opt$g,".raw")) rm(list=paste0("sc10x.",opt$g,".raw")) -sc10x.Group <- SetAllIdent(object=sc10x.Goup,id="samples") +sc10x.Group <- SetAllIdent(object=sc10x.Group,id="samples") sc10x.Group@ident <- factor(sc10x.Group@ident,levels=c("D17PrPzF_BE","D17PrPzF_LE","D17PrPzF_OE","D17PrPzF_FMSt","D17PrTzF_Via","D17PrPzF_Via")) assign(paste0("sc10x.",opt$g,".raw"),sc10x.Group) rm(sc10x.Group) #resave Rda file -save(list=paste0("sc10x.",opt$g,".raw"),file=paste0("./sc10x.",Project.Name,".",opt$g,".raw.Rda")) +save(list=paste0("sc10x.",opt$g,".raw"),file=paste0("./analysis/sc10x.",Project.Name,".",opt$g,".raw.Rda")) rm(list=paste0("sc10x.",opt$g,".raw")) -rm(opt) \ No newline at end of file +rm(opt) diff --git a/r.scripts/sc_DEG.R b/r.scripts/sc_DEG.R index 74443e2..9333b15 100755 --- a/r.scripts/sc_DEG.R +++ b/r.scripts/sc_DEG.R @@ -21,24 +21,24 @@ if (opt$st==TRUE){ } #load data -setwd(paste0("./",Project.Name)) -load(paste0("./sc10x.",Project.Name,".Cumulative.RData")) -load(paste0("./sc10x.",Project.Name,".",opt$g,".cluster",sub.file,".IDepi+st+ne+LGEAepi+LGEAst.Rda")) +setwd("../") +load(paste0("./analysis/sc10x.",Project.Name,".Cumulative.RData")) +load(paste0("./analysis/sc10x.",Project.Name,".",opt$g,".cluster",sub.file,".IDepi+st+ne+LGEAepi+LGEAst.Rda")) sc10x.Group.NoStress <- get(paste0("sc10x.",opt$g,".cluster",sub.file,".IDepi+st+ne+LGEAepi+LGEAst")) rm(list=paste0("sc10x.",opt$g,".cluster",sub.file,".IDepi+st+ne+LGEAepi+LGEAst")) -load(paste0("./sc10x.",Project.Name,".",opt$g,".cluster.IDStress.Rda")) +load(paste0("./analysis/sc10x.",Project.Name,".",opt$g,".cluster.IDStress.Rda")) sc10x.Group.Stress <- get(paste0("sc10x.",opt$g,".cluster.IDStress")) rm(list=paste0("sc10x.",opt$g,".cluster.IDStress")) #create folder structure -if (!dir.exists(paste0("./",opt$g))){ - dir.create(paste0("./",opt$g)) +if (!dir.exists(paste0("./analysis/",opt$g))){ + dir.create(paste0("./analysis/",opt$g)) } -if (!dir.exists(paste0("./",opt$g,"/global"))){ - dir.create(paste0("./",opt$g,"/global")) +if (!dir.exists(paste0("./analysis/",opt$g,"/global"))){ + dir.create(paste0("./analysis/",opt$g,"/global")) } -if (!dir.exists(paste0("./",opt$g,"/global/DEG"))){ - dir.create(paste0("./",opt$g,"/global/DEG")) +if (!dir.exists(paste0("./analysis/",opt$g,"/global/DEG"))){ + dir.create(paste0("./analysis/",opt$g,"/global/DEG")) } DEG.Stress <- FindMarkers(object=sc10x.Group.Stress,ident.1="Stress",ident.2="ALL",min.pct=0.5,only.pos=TRUE,thresh.use=0.25,logfc.threshold=log(5)) @@ -80,33 +80,33 @@ DEG.SM2.vs.SM1 <- DEG.SM2.vs.SM1[DEG.SM2.vs.SM1$p_val_adj<0.05,] sc10x.Group.NoStress <- SetAllIdent(object=sc10x.Group.NoStress,id="ALLSubPop") for (i in ls(pattern="^DEG")[!(ls(pattern="^DEG") %in% c("DEG.NE.vs.Epi","DEG.Stress"))]){ - png(paste0("./",opt$g,"/global/DEG/Heatmap.",i,".png"),1000,1000,pointsize=20,type="quartz",antialias="none",family="Arial") + png(paste0("./analysis/",opt$g,"/global/DEG/Heatmap.",i,".png"),1000,1000,pointsize=20,type="quartz",antialias="none",family="Arial") plot(DoHeatmap(sc10x.Group.NoStress,genes.use=row.names(get(i)),slim.col.label=TRUE,group.label.rot=FALSE,group.spacing=0.5,group.order=c("BE","LE","OE","SM1","SM2","Endo"),cex.row=10,group.cex=25)) dev.off() } -png(paste0("./",opt$g,"/global/DEG/Heatmap.DEG.Stress.png"),1000,1000,pointsize=20,type="quartz",antialias="none",family="Arial") +png(paste0("./analysis/",opt$g,"/global/DEG/Heatmap.DEG.Stress.png"),1000,1000,pointsize=20,type="quartz",antialias="none",family="Arial") plot(DoHeatmap(sc10x.Group.Stress,genes.use=row.names(DEG.Stress),slim.col.label=TRUE,group.label.rot=FALSE,group.spacing=0.5,group.order=c("Stress","ALL"),cex.row=10,group.cex=25)) dev.off() sc10x.Group.NoStress <- SetAllIdent(object=sc10x.Group.NoStress,id="ALLSubPop+NE") for (i in ls(pattern="^DEG")[!(ls(pattern="^DEG") %in% c("DEG.Stress"))]){ - png(paste0("./",opt$g,"/global/DEG/Violin.",i,".png"),1000,1000,pointsize=20,type="cairo",antialias="none",family="Arial") + png(paste0("./analysis/",opt$g,"/global/DEG/Violin.",i,".png"),1000,1000,pointsize=20,type="cairo",antialias="none",family="Arial") plot(VlnPlot(sc10x.Group.NoStress,features.plot=row.names(get(i))[1:10],nCol=5)) dev.off() - png(paste0("./",opt$g,"/global/DEG/Ridge.",i,".png"),1000,1000,pointsize=20,type="cairo",antialias="none",family="Arial") + png(paste0("./analysis/",opt$g,"/global/DEG/Ridge.",i,".png"),1000,1000,pointsize=20,type="cairo",antialias="none",family="Arial") plot(RidgePlot(sc10x.Group.NoStress,features.plot=row.names(get(i))[1:10],nCol=5)) dev.off() } -png(paste0("./",opt$g,"/global/DEG/Violin.DEG.Stress.png"),1000,1000,pointsize=20,type="cairo",antialias="none",family="Arial") +png(paste0("./analysis/",opt$g,"/global/DEG/Violin.DEG.Stress.png"),1000,1000,pointsize=20,type="cairo",antialias="none",family="Arial") plot(VlnPlot(sc10x.Group.Stress,features.plot=row.names(DEG.Stress)[1:10],nCol=5)) dev.off() -png(paste0("./",opt$g,"/global/DEG/Ridge.DEG.Stress.png"),1000,1000,pointsize=20,type="cairo",antialias="none",family="Arial") +png(paste0("./analysis/",opt$g,"/global/DEG/Ridge.DEG.Stress.png"),1000,1000,pointsize=20,type="cairo",antialias="none",family="Arial") plot(RidgePlot(sc10x.Group.Stress,features.plot=row.names(DEG.Stress)[1:10],nCol=5)) dev.off() for (i in ls(pattern="^DEG")){ - write.table(get(i),file=paste0("./",opt$g,"/global/DEG/",i,".csv"),row.names=TRUE,col.names=NA,append=FALSE,sep=",") + write.table(get(i),file=paste0("./analysis/",opt$g,"/global/DEG/",i,".csv"),row.names=TRUE,col.names=NA,append=FALSE,sep=",") } -cat(row.names(DEG.Stress),file=paste0("./",opt$g,"/global/DEG/DWS.scStress.txt"),sep="\n") -cat(row.names(DEG.NE.vs.Epi),file=paste0("./",opt$g,"/global/DEG/DWS.scNE.txt"),sep="\n") +cat(row.names(DEG.Stress),file=paste0("./analysis/",opt$g,"/global/DEG/DWS.scStress.txt"),sep="\n") +cat(row.names(DEG.NE.vs.Epi),file=paste0("./analysis/",opt$g,"/global/DEG/DWS.scNE.txt"),sep="\n") diff --git a/r.scripts/sc_Demultiplex.R b/r.scripts/sc_Demultiplex.R index 7e493d7..c9dfcb4 100755 --- a/r.scripts/sc_Demultiplex.R +++ b/r.scripts/sc_Demultiplex.R @@ -16,10 +16,10 @@ rm(option_list) Project.Name <- opt$p #load 10x data and map library_id to samples identity... initial QC min.cells/min.genes loaded based on command line options (default=3/200) -setwd(paste0("./",Project.Name)) -sc10x.data <- Read10X(data.dir="./10x/filtered_gene_bc_matrices_mex/GRCh38/") +setwd("../") +sc10x.data <- Read10X(data.dir="./analysis/DATA/10x/filtered_gene_bc_matrices_mex/GRCh38/") sc10x <- new("seurat",raw.data=sc10x.data) -sc10x.aggr <- read_csv(paste0("./10x/",Project.Name,"-aggr.csv")) +sc10x.aggr <- read_csv(paste0("./analysis/DATA/10x/",Project.Name,"-aggr.csv")) cell.codes <- as.data.frame(sc10x@raw.data@Dimnames[[2]]) colnames(cell.codes) <- "barcodes" rownames(cell.codes) <- cell.codes$barcodes @@ -31,8 +31,8 @@ rm(cell.codes) rm(sc10x.aggr) #map sample categories based on demultiplex.csv and the command line option defined number of categories and save raw.Rda's for each category -sc10x.demultiplex <- read_csv(paste0("./",Project.Name,"-demultiplex.csv")) -Group.List <- c("ALL","samples",colnames(sc10x.demultiplex[,2:(opt$d+1)])) +sc10x.demultiplex <- read_csv(paste0("./analysis/DATA/",Project.Name,"-demultiplex.csv")) +Group.List <- c("samples",colnames(sc10x.demultiplex[,2:(opt$d+1)])) for (i in 2:(opt$d+1)){ sc10x <- SetAllIdent(object=sc10x,id="samples") merge.cluster <- apply(sc10x.demultiplex[,i],1,as.character) @@ -40,7 +40,7 @@ for (i in 2:(opt$d+1)){ sc10x@ident <- plyr::mapvalues(x=sc10x@ident,from=sc10x.demultiplex$Samples,to=merge.cluster) sc10x <- StashIdent(object=sc10x,save.name=colnames(sc10x.demultiplex[,i])) assign(paste0("sc10x.",colnames(sc10x.demultiplex[,i]),".raw"),SubsetData(object=sc10x,ident.use=colnames(sc10x.demultiplex[,i]))) - save(list=paste0("sc10x.",colnames(sc10x.demultiplex[,i]),".raw"),file=paste0("./sc10x.",Project.Name,".",colnames(sc10x.demultiplex[,i]),".raw.Rda")) + save(list=paste0("sc10x.",colnames(sc10x.demultiplex[,i]),".raw"),file=paste0("./analysis/sc10x.",Project.Name,".",colnames(sc10x.demultiplex[,i]),".raw.Rda")) rm(list=paste0("sc10x.",colnames(sc10x.demultiplex[,i]),".raw")) } rm(sc10x.demultiplex) @@ -52,6 +52,6 @@ opt.demultiplex <- opt rm(opt) #save raw.Rda file,Cumulative.RData file -save(sc10x,file=paste0("./sc10x.",Project.Name,".raw.Rda")) +save(sc10x,file=paste0("./analysis/sc10x.",Project.Name,".raw.Rda")) rm(sc10x) -save.image(file=paste0("./sc10x.",Project.Name,".Cumulative.RData")) +save.image(file=paste0("./analysis/sc10x.",Project.Name,".Cumulative.RData")) diff --git a/r.scripts/sc_PC.Score.NE.R b/r.scripts/sc_PC.Score.NE.R index 4b8827d..2176e1d 100755 --- a/r.scripts/sc_PC.Score.NE.R +++ b/r.scripts/sc_PC.Score.NE.R @@ -24,24 +24,24 @@ if (opt$st==TRUE){ } #load data -setwd(paste0("./",Project.Name)) -load(paste0("./sc10x.",Project.Name,".Cumulative.RData")) -load(paste0("./sc10x.",Project.Name,".",opt$g,".cluster",sub.file,".IDepi+st.Rda")) +setwd("../") +load(paste0("./analysis/sc10x.",Project.Name,".Cumulative.RData")) +load(paste0("./analysis/sc10x.",Project.Name,".",opt$g,".cluster",sub.file,".IDepi+st.Rda")) sc10x.Group <- get(paste0("sc10x.",opt$g,".cluster",sub.file,".IDepi+st")) rm(list=paste0("sc10x.",opt$g,".cluster",sub.file,".IDepi+st")) #create folder structure -if (!dir.exists(paste0("./",opt$g))){ - dir.create(paste0("./",opt$g)) +if (!dir.exists(paste0("./analysis/",opt$g))){ + dir.create(paste0("./analysis/",opt$g)) } -if (!dir.exists(paste0("./",opt$g,"/global"))){ - dir.create(paste0("./",opt$g,"/global")) +if (!dir.exists(paste0("./analysis/",opt$g,"/global"))){ + dir.create(paste0("./analysis/",opt$g,"/global")) } -if (!dir.exists(paste0("./",opt$g,"/global",sub.folder))){ - dir.create(paste0("./",opt$g,"/global",sub.folder)) +if (!dir.exists(paste0("./analysis/",opt$g,"/global",sub.folder))){ + dir.create(paste0("./analysis/",opt$g,"/global",sub.folder)) } -if (!dir.exists(paste0("./",opt$g,"/global",sub.folder,"/pc-score"))){ - dir.create(paste0("./",opt$g,"/global",sub.folder,"/pc-score")) +if (!dir.exists(paste0("./analysis/",opt$g,"/global",sub.folder,"/pc-score"))){ + dir.create(paste0("./analysis/",opt$g,"/global",sub.folder,"/pc-score")) } #load resolution @@ -49,9 +49,9 @@ sc10x.Group <- SetAllIdent(object=sc10x.Group,id="Lineage") #subset by NE genes from Table 1 of Vashchenko,Nadezda,and Per-Anders Abrahamsson. "Neuroendocrine differentiation in prostate cancer: implications for new treatment modalities." European urology 47.2 (2005): 147-155 or DWS_scNE if (opt$u==FALSE){ - DEG.NE <- read_delim("../EurUrol.2005.NE.txt","\t",escape_double=FALSE,trim_ws=TRUE,col_names=FALSE) + DEG.NE <- read_delim("./genesets/EurUrol.2005.NE.txt","\t",escape_double=FALSE,trim_ws=TRUE,col_names=FALSE) } else { - DEG.NE <- read_delim("../DWS.scNE.txt","\t",escape_double=FALSE,trim_ws=TRUE,col_names=FALSE) + DEG.NE <- read_delim("./genesets/DWS.scNE.txt","\t",escape_double=FALSE,trim_ws=TRUE,col_names=FALSE) } DEG.NE <- as.list(DEG.NE) sc10x.Group.NE <- as.data.frame(as.matrix(sc10x.Group@data[rownames(sc10x.Group@data) %in% unlist(DEG.NE),])) @@ -64,26 +64,26 @@ sc10x.Group <- SetDimReduction(object=sc10x.Group,reduction.type="NE",slot="cell sc10x.Group <- SetDimReduction(object=sc10x.Group,reduction.type="NE",slot="key",new.data="NE") #generate PCA PC1 (NE1/NE-Score) and PC2 (NE2) and histogram of NE-Scores -png(paste0("./",opt$g,"/global",sub.folder,"/pc-score","/PCS_NE.ALL.png"),1000,1000,pointsize=20) +png(paste0("./analysis/",opt$g,"/global",sub.folder,"/pc-score","/PCS_NE.ALL.png"),1000,1000,pointsize=20) plot <- DimPlot(object=sc10x.Group,reduction.use="NE",pt.size=2.5,do.return=TRUE) 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() -png(paste0("./",opt$g,"/global",sub.folder,"/pc-score","/Hist_NE.ALL.png"),1000,1000,pointsize=20) +png(paste0("./analysis/",opt$g,"/global",sub.folder,"/pc-score","/Hist_NE.ALL.png"),1000,1000,pointsize=20) histo <- hist(GetCellEmbeddings(object=sc10x.Group,reduction.type="NE",dims.use=1),breaks=10,plot=FALSE) plot(histo$mids,histo$density,log="y",type="h",lwd=20,lend=2,xlab="NE Score",ylab="log10 Frequency") dev.off() #subset to Epi only and regenerate figures above sc10x.Group.Epi <- SubsetData(object=sc10x.Group,ident.use="Epi") -png(paste0("./",opt$g,"/global",sub.folder,"/pc-score","/PCS_NE.Epi.png"),1000,1000,pointsize=20) +png(paste0("./analysis/",opt$g,"/global",sub.folder,"/pc-score","/PCS_NE.Epi.png"),1000,1000,pointsize=20) plot <- DimPlot(object=sc10x.Group.Epi,reduction.use="NE",pt.size=5,do.return=TRUE) 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() -png(paste0("./",opt$g,"/global",sub.folder,"/pc-score","/Hist_NE.Epi.png"),1000,1000,pointsize=20) +png(paste0("./analysis/",opt$g,"/global",sub.folder,"/pc-score","/Hist_NE.Epi.png"),1000,1000,pointsize=20) histo <- hist(GetCellEmbeddings(object=sc10x.Group.Epi,reduction.type="NE",dims.use=1),breaks=10,plot=FALSE) plot(histo$mids,histo$density,log="y",type="h",lwd=20,lend=2,xlab="NE Score",ylab="log10 Frequency") dev.off() @@ -99,7 +99,7 @@ if (clust$size[1]>clust$size[2]){ } else { sc10x.Group.Epi <- SetIdent(object=sc10x.Group.Epi,cells.use=attr(clusterCut[clusterCut==as.integer(names(clust$centers[clust$centers==min(clust$centers),]))],"names"),ident.use="NE") } -png(paste0("./",opt$g,"/global",sub.folder,"/pc-score","/PCS_NE.Epi.Group.png"),1000,1000,pointsize=20) +png(paste0("./analysis/",opt$g,"/global",sub.folder,"/pc-score","/PCS_NE.Epi.Group.png"),1000,1000,pointsize=20) plot <- DimPlot(object=sc10x.Group.Epi,reduction.use="NE",pt.size=2.5,do.return=TRUE) 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))) @@ -120,7 +120,7 @@ sc10x.Group.sample <- SetAllIdent(object=sc10x.Group.sample,id="Lineage") sc10x.Group.sample <- SetIdent(object=sc10x.Group.sample,cells.use=clusterCut,ident.use="NE") #generate tSNE plot of subsampled cells + NE -png(paste0("./",opt$g,"/global",sub.folder,"/tSNE_NE.SubSample.png"),1000,1000,pointsize=20) +png(paste0("./analysis/",opt$g,"/global",sub.folder,"/tSNE_NE.SubSample.png"),1000,1000,pointsize=20) plot <- TSNEPlot(object=sc10x.Group.sample,pt.size=5,do.label=FALSE,label.size=15,do.return=TRUE) 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))) @@ -133,11 +133,11 @@ sc10x.Group@ident <- factor(sc10x.Group@ident,levels=c("NE","Epi","St")) #generate violin plot of NE gene expression anchors <- c("CHGA","CHGB","SYP") for (i in anchors){ -png(paste0("./",opt$g,"/global",sub.folder,"/pc-score","/Violin_NE.",i,".png"),1000,1000,pointsize=20) + png(paste0("./analysis/",opt$g,"/global",sub.folder,"/pc-score","/Violin_NE.",i,".png"),1000,1000,pointsize=20) plot <- VlnPlot(object=sc10x.Group.sample,features.plot=i,do.return=TRUE) 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)) plot(plot) -dev.off() + dev.off() } rm(anchors) @@ -174,8 +174,8 @@ opt.ne <- opt rm(opt) #save full .EpiPop+NE.Rda file,Cumulative.RData files -save(list=paste0("sc10x.",opt.ne$g,".cluster",sub.file,".IDepi+st+ne"),file=paste0("./sc10x.",Project.Name,".",opt.ne$g,".cluster",sub.file,".IDepi+st+ne.Rda")) +save(list=paste0("sc10x.",opt.ne$g,".cluster",sub.file,".IDepi+st+ne"),file=paste0("./analysis/sc10x.",Project.Name,".",opt.ne$g,".cluster",sub.file,".IDepi+st+ne.Rda")) rm(list=paste0("sc10x.",opt.ne$g,".cluster",sub.file,".IDepi+st+ne")) rm(sub.folder) rm(sub.file) -save.image(file=paste0("./sc10x.",Project.Name,".Cumulative.RData")) +save.image(file=paste0("./analysis/sc10x.",Project.Name,".Cumulative.RData")) diff --git a/r.scripts/sc_PC.Score.Stress.R b/r.scripts/sc_PC.Score.Stress.R index 1673c8a..40629aa 100755 --- a/r.scripts/sc_PC.Score.Stress.R +++ b/r.scripts/sc_PC.Score.Stress.R @@ -20,21 +20,21 @@ Project.Name <- opt$p PC.Max <- opt$pc #load data -setwd(paste0("./",Project.Name)) -load(paste0("./sc10x.",Project.Name,".Cumulative.RData")) -load(paste0("./sc10x.",Project.Name,".",opt$g,".",opt$id,".Rda")) +setwd("../") +load(paste0("./analysis/sc10x.",Project.Name,".Cumulative.RData")) +load(paste0("./analysis/sc10x.",Project.Name,".",opt$g,".",opt$id,".Rda")) sc10x.Group <- get(paste0("sc10x.",opt$g,".",opt$id)) rm(list=paste0("sc10x.",opt$g,".",opt$id)) #create folder structure -if (!dir.exists(paste0("./",opt$g))){ - dir.create(paste0("./",opt$g)) +if (!dir.exists(paste0("./analysis/",opt$g))){ + dir.create(paste0("./analysis/",opt$g)) } -if (!dir.exists(paste0("./",opt$g,"/pc-score"))){ - dir.create(paste0("./",opt$g,"/pc-score")) +if (!dir.exists(paste0("./analysis/",opt$g,"/pc-score"))){ + dir.create(paste0("./analysis/",opt$g,"/pc-score")) } -if (!dir.exists(paste0("./",opt$g,"/global/NoStress"))){ - dir.create(paste0("./",opt$g,"/global/NoStress")) +if (!dir.exists(paste0("./analysis/",opt$g,"/global/NoStress"))){ + dir.create(paste0("./analysis/",opt$g,"/global/NoStress")) } #load resolution @@ -42,10 +42,10 @@ sc10x.Group <- SetAllIdent(object=sc10x.Group,id="ALL") #subset by Stress genes from CHUANG_OXIDATIVE_STRESS_RESPONSE_UP or DWS_scStress if (opt$u==FALSE){ - DEG.Stress <- read_csv("../DEG_C2.CGP.M10970.txt") + DEG.Stress <- read_csv("./genesets/DEG_C2.CGP.M10970.txt") DEG.Stress <- DEG.Stress[-1,] } else { - DEG.Stress <- read_delim("../DWS.scStress.txt","\t",escape_double=FALSE,trim_ws=TRUE,col_names=FALSE) + DEG.Stress <- read_delim("./genesets/DWS.scStress.txt","\t",escape_double=FALSE,trim_ws=TRUE,col_names=FALSE) } DEG.Stress <- as.list(DEG.Stress) sc10x.Group.Stress <- as.data.frame(as.matrix(sc10x.Group@data[rownames(sc10x.Group@data) %in% unlist(DEG.Stress),])) @@ -58,14 +58,14 @@ sc10x.Group <- SetDimReduction(object=sc10x.Group,reduction.type="Stress",slot=" sc10x.Group <- SetDimReduction(object=sc10x.Group,reduction.type="Stress",slot="key",new.data="Stress") #generate PCA PC1 (Stress1/Stress-Score) and PC2 (Stress2) and histogram of Stress-Scores -png(paste0("./",opt$g,"/pc-score/PCS_Stress.png"),1000,1000,pointsize=20) +png(paste0("./analysis/",opt$g,"/pc-score/PCS_Stress.png"),1000,1000,pointsize=20) plot <- DimPlot(object=sc10x.Group,reduction.use="Stress",pt.size=2.5,do.return=TRUE) plot <- plot+geom_density2d(color="black",bins=25,alpha=0.5) 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() -png(paste0("./",opt$g,"/pc-score/Hist_Stress.png"),1000,1000,pointsize=20) +png(paste0("./analysis/",opt$g,"/pc-score/Hist_Stress.png"),1000,1000,pointsize=20) histo <- hist(GetCellEmbeddings(object=sc10x.Group,reduction.type="Stress",dims.use=1),breaks=10,plot=FALSE) plot(histo$mids,histo$density,log="y",type="h",lwd=20,lend=2,xlab="Stress Score",ylab="log10 Frequency") dev.off() @@ -81,7 +81,7 @@ if (clust$size[1]<clust$size[2]){ } sc10x.Group <- SetAllIdent(object=sc10x.Group,id="ALL") sc10x.Group <- SetIdent(object=sc10x.Group,cells.use=clusterCut,ident.use="Stress") -png(paste0("./",opt$g,"/pc-score/PCS_Stress.Group.png"),1000,1000,pointsize=20) +png(paste0("./analysis/",opt$g,"/pc-score/PCS_Stress.Group.png"),1000,1000,pointsize=20) plot <- DimPlot(object=sc10x.Group,reduction.use="Stress",pt.size=2.5,do.return=TRUE) plot <- plot+geom_density2d(color="black",bins=25,alpha=0.5) 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)) @@ -98,7 +98,7 @@ sc10x.Group.sample <- SetAllIdent(object=sc10x.Group.sample,id="ALL") sc10x.Group.sample <- SetIdent(object=sc10x.Group.sample,cells.use=clusterCut,ident.use="Stress") #generate tSNE plot of subsampled cells + Stress -png(paste0("./",opt$g,"/global/tSNE_Stress.SubSample.png"),1000,1000,pointsize=20) +png(paste0("./analysis/",opt$g,"/global/tSNE_Stress.SubSample.png"),1000,1000,pointsize=20) plot <- TSNEPlot(object=sc10x.Group.sample,pt.size=5,do.label=FALSE,label.size=15,do.return=TRUE) 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))) @@ -109,7 +109,7 @@ sc10x.Group <- SetIdent(object=sc10x.Group,cells.use=clusterCut,ident.use="Stres #generate violin plot of stres gene expression anchors <- c("EGR1","FOS","JUN") for (i in anchors){ - png(paste0("./",opt$g,"/pc-score/Violin_Stress.",i,".png"),1000,1000,pointsize=20) + png(paste0("./analysis/",opt$g,"/pc-score/Violin_Stress.",i,".png"),1000,1000,pointsize=20) plot <- VlnPlot(object=sc10x.Group.sample,features.plot=i,do.return=TRUE) 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)) plot(plot) @@ -135,7 +135,7 @@ sc10x.Group.subset <- SubsetData(object=sc10x.Group,ident.use="ALL") sc10x.Group.subset <- SetAllIdent(object=sc10x.Group.subset,id=paste0("res",opt$r)) sc10x.Group.subset <- BuildClusterTree(sc10x.Group.subset,do.reorder=TRUE,reorder.numeric=TRUE,do.plot=FALSE) sc10x.Group.subset <- StashIdent(object=sc10x.Group.subset,save.name=paste0("res",opt$r,"-Stress")) -png(paste0("./",opt$g,"/global/tSNE_","res",opt$r,"-NoStress.png"),1000,1000,pointsize=20) +png(paste0("./analysis/",opt$g,"/global/tSNE_","res",opt$r,"-NoStress.png"),1000,1000,pointsize=20) plot <- TSNEPlot(object=sc10x.Group.subset,pt.size=5,do.label=FALSE,label.size=15,do.return=TRUE) 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))) @@ -146,7 +146,7 @@ rm(merge.cluster) #regenerate tSNEs for range of resolutions sc10x.Group.subset <- RunTSNE(object=sc10x.Group.subset,dims.use=1:PC.Max,do.fast=TRUE,force.recalc=TRUE) -png(paste0("./",opt$g,"/global/NoStress/tSNE_Sample.png"),1000,1000,pointsize=20) +png(paste0("./analysis/",opt$g,"/global/NoStress/tSNE_Sample.png"),1000,1000,pointsize=20) plot <- TSNEPlot(object=sc10x.Group.subset,group.by="samples",pt.size=2.5,do.return=TRUE) 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))) @@ -156,10 +156,10 @@ for (i in c(5,2,1,0.5,0.2,0.1)){ sc10x.Group.subset <- FindClusters(object=sc10x.Group.subset,reduction.type="pca",dims.use=1:PC.Max,resolution=i,print.output=0,save.SNN=TRUE) sc10x.Group.subset <- BuildClusterTree(sc10x.Group.subset,do.reorder=TRUE,reorder.numeric=TRUE,do.plot=FALSE) sc10x.Group.subset <- StashIdent(object=sc10x.Group.subset,save.name=paste0("res",i)) - png(paste0("./",opt$g,"/global/NoStress/Tree_Global.res",i,".png"),1000,1000,pointsize=20) + png(paste0("./analysis/",opt$g,"/global/NoStress/Tree_Global.res",i,".png"),1000,1000,pointsize=20) PlotClusterTree(object=sc10x.Group.subset) dev.off() - png(paste0("./",opt$g,"/global/NoStress/tSNE_Global.res",i,".png"),1000,1000,pointsize=20) + png(paste0("./analysis/",opt$g,"/global/NoStress/tSNE_Global.res",i,".png"),1000,1000,pointsize=20) plot <- TSNEPlot(object=sc10x.Group.subset,pt.size=5,do.label=TRUE,label.size=10,do.return=TRUE) 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))) @@ -192,8 +192,8 @@ opt.stress <- opt rm(opt) #save full .EpiPop+Stress.Rda file,Cumulative.RData files -save(list=paste0("sc10x.",opt.stress$g,".",opt.stress$id,".IDStress"),file=paste0("./sc10x.",Project.Name,".",opt.stress$g,".",opt.stress$id,".IDStress.Rda")) -save(list=paste0("sc10x.",opt.stress$g,".",opt.stress$id,".NOStress"),file=paste0("./sc10x.",Project.Name,".",opt.stress$g,".",opt.stress$id,".NOStress.Rda")) +save(list=paste0("sc10x.",opt.stress$g,".",opt.stress$id,".IDStress"),file=paste0("./analysis/sc10x.",Project.Name,".",opt.stress$g,".",opt.stress$id,".IDStress.Rda")) +save(list=paste0("sc10x.",opt.stress$g,".",opt.stress$id,".NOStress"),file=paste0("./analysis/sc10x.",Project.Name,".",opt.stress$g,".",opt.stress$id,".NOStress.Rda")) rm(list=paste0("sc10x.",opt.stress$g,".",opt.stress$id,".IDStress")) rm(list=paste0("sc10x.",opt.stress$g,".",opt.stress$id,".NOStress")) -save.image(file=paste0("./sc10x.",Project.Name,".Cumulative.RData")) +save.image(file=paste0("./analysis/sc10x.",Project.Name,".Cumulative.RData")) diff --git a/r.scripts/sc_QC.R b/r.scripts/sc_QC.R index a4294fd..a7acbed 100755 --- a/r.scripts/sc_QC.R +++ b/r.scripts/sc_QC.R @@ -22,24 +22,24 @@ rm(option_list) Project.Name <- opt$p #load data -setwd(paste0("./",Project.Name)) -load(paste0("./sc10x.",Project.Name,".Cumulative.RData")) +setwd("../") +load(paste0("./analysis/sc10x.",Project.Name,".Cumulative.RData")) if (opt$cc==TRUE){ - load(paste0("./sc10x.",Project.Name,".",opt$g,".cc.Rda")) + load(paste0("./analysis/sc10x.",Project.Name,".",opt$g,".cc.Rda")) sc10x.Group <- get(paste0("sc10x.",opt$g,".cc")) rm(list=paste0("sc10x.",opt$g,".cc")) } else { - load(paste0("./sc10x.",Project.Name,".",opt$g,".raw.Rda")) + load(paste0("./analysis/sc10x.",Project.Name,".",opt$g,".raw.Rda")) sc10x.Group <- get(paste0("sc10x.",opt$g,".raw")) rm(list=paste0("sc10x.",opt$g,".raw")) } #create folder structure -if (!dir.exists(paste0("./",opt$g))){ - dir.create(paste0("./",opt$g)) +if (!dir.exists(paste0("./analysis/",opt$g))){ + dir.create(paste0("./analysis/",opt$g)) } -if (!dir.exists(paste0("./",opt$g,"/qc"))){ - dir.create(paste0("./",opt$g,"/qc")) +if (!dir.exists(paste0("./analysis/",opt$g,"/qc"))){ + dir.create(paste0("./analysis/",opt$g,"/qc")) } #filter cells @@ -47,11 +47,11 @@ mito.genes <- grep(pattern="^MT-",x=rownames(x=sc10x.Group@data),value=TRUE) percent.mito <- colSums(sc10x.Group@raw.data[mito.genes,])/colSums(sc10x.Group@raw.data) sc10x.Group <- AddMetaData(object=sc10x.Group,metadata=percent.mito,col.name="percent.mito") sc10x.Group <- SetAllIdent(object=sc10x.Group,id="ALL") -png(paste0("./",opt$g,"/qc/Violin_qc.raw.png"),1000,1000,pointsize=20) +png(paste0("./analysis/",opt$g,"/qc/Violin_qc.raw.png"),1000,1000,pointsize=20) plot <- VlnPlot(object=sc10x.Group,features.plot=c("nGene","nUMI","percent.mito"),nCol=3,do.return=TRUE) plot(plot) dev.off() -png(paste0("./",opt$g,"/qc/Plot_qc.raw.png"),1000,1000,pointsize=20) +png(paste0("./analysis/",opt$g,"/qc/Plot_qc.raw.png"),1000,1000,pointsize=20) par(mfrow=c(1,2)) GenePlot(object=sc10x.Group,gene1="nUMI",gene2="percent.mito") GenePlot(object=sc10x.Group,gene1="nUMI",gene2="nGene") @@ -60,11 +60,11 @@ sc10x.Group.RawCellCount <- length(sc10x.Group@data@Dimnames[[2]]) sc10x.Group.RawGeneCount <- length(sc10x.Group@data@Dimnames[[1]]) sc10x.Group <- FilterCells(object=sc10x.Group,subset.names=c("nGene","percent.mito"),low.thresholds=c(opt$lg,opt$lm),high.thresholds=c(opt$hg,opt$hm)) sc10x.Group <- NormalizeData(object=sc10x.Group,normalization.method="LogNormalize",scale.factor=10000) -png(paste0("./",opt$g,"/qc/Violin_qc.filtered.png"),1000,1000,pointsize=20) +png(paste0("./analysis/",opt$g,"/qc/Violin_qc.filtered.png"),1000,1000,pointsize=20) plot <- VlnPlot(object=sc10x.Group,features.plot=c("nGene","nUMI","percent.mito"),nCol=3,do.return=TRUE) plot(plot) dev.off() -png(paste0("./",opt$g,"/qc/Plot_qc.filtered.png"),1000,1000,pointsize=20) +png(paste0("./analysis/",opt$g,"/qc/Plot_qc.filtered.png"),1000,1000,pointsize=20) par(mfrow=c(1,2)) GenePlot(object=sc10x.Group,gene1="nUMI",gene2="percent.mito") GenePlot(object=sc10x.Group,gene1="nUMI",gene2="nGene") @@ -79,7 +79,7 @@ if (opt$cc==TRUE){ gc() sc10x.Group <- RunPCA(object=sc10x.Group,pc.genes=c(s.genes,g2m.genes),do.print=FALSE,pcs.store=2) sc10x.Group <- SetAllIdent(object=sc10x.Group,id="Phase") - png(paste0("./",opt$g,"/cc/PCA_cc.Norm.png"),1000,1000,pointsize=20) + png(paste0("./analysis/",opt$g,"/cc/PCA_cc.Norm.png"),1000,1000,pointsize=20) plot <- PCAPlot(object=sc10x.Group,do.return=TRUE) 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)) plot <- plot+guides(colour=guide_legend(override.aes=list(size=10))) @@ -97,45 +97,45 @@ rm(mito.genes) rm(percent.mito) #generate qc figures -png(paste0("./",opt$g,"/qc/Plot_PCElbow.png"),1000,1000,pointsize=20) +png(paste0("./analysis/",opt$g,"/qc/Plot_PCElbow.png"),1000,1000,pointsize=20) plot <- PCElbowPlot(object=sc10x.Group,num.pc=50) 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)) plot(plot) dev.off() -png(paste0("./",opt$g,"/qc/Plot_VizPCA.01.09.png"),1000,1000,pointsize=20) +png(paste0("./analysis/",opt$g,"/qc/Plot_VizPCA.01.09.png"),1000,1000,pointsize=20) VizPCA(object=sc10x.Group,pcs.use=1:9) dev.off() -png(paste0("./",opt$g,"/qc/Plot_VizPCA.10.18.png"),1000,1000,pointsize=20) +png(paste0("./analysis/",opt$g,"/qc/Plot_VizPCA.10.18.png"),1000,1000,pointsize=20) VizPCA(object=sc10x.Group,pcs.use=10:18) dev.off() -png(paste0("./",opt$g,"/qc/Plot_VizPCA.19.27.png"),1000,1000,pointsize=20) +png(paste0("./analysis/",opt$g,"/qc/Plot_VizPCA.19.27.png"),1000,1000,pointsize=20) VizPCA(object=sc10x.Group,pcs.use=19:27) dev.off() -png(paste0("./",opt$g,"/qc/Plot_VizPCA.28.36.png"),1000,1000,pointsize=20) +png(paste0("./analysis/",opt$g,"/qc/Plot_VizPCA.28.36.png"),1000,1000,pointsize=20) VizPCA(object=sc10x.Group,pcs.use=28:36) dev.off() -png(paste0("./",opt$g,"/qc/Plot_VizPCA.37.45.png"),1000,1000,pointsize=20) +png(paste0("./analysis/",opt$g,"/qc/Plot_VizPCA.37.45.png"),1000,1000,pointsize=20) VizPCA(object=sc10x.Group,pcs.use=37:45) dev.off() -png(paste0("./",opt$g,"/qc/Plot_VizPCA.46.50.png"),1000,1000,pointsize=20) +png(paste0("./analysis/",opt$g,"/qc/Plot_VizPCA.46.50.png"),1000,1000,pointsize=20) VizPCA(object=sc10x.Group,pcs.use=46:50) dev.off() -png(paste0("./",opt$g,"/qc/Heatmap_PC.01.09.png"),1000,1000,pointsize=20) +png(paste0("./analysis/",opt$g,"/qc/Heatmap_PC.01.09.png"),1000,1000,pointsize=20) PCHeatmap(object=sc10x.Group,pc.use=1:9,cells.use=500,do.balanced=TRUE,label.columns=FALSE,use.full=FALSE) dev.off() -png(paste0("./",opt$g,"/qc/Heatmap_PC.10.18.png"),1000,1000,pointsize=20) +png(paste0("./analysis/",opt$g,"/qc/Heatmap_PC.10.18.png"),1000,1000,pointsize=20) PCHeatmap(object=sc10x.Group,pc.use=10:18,cells.use=500,do.balanced=TRUE,label.columns=FALSE,use.full=FALSE) dev.off() -png(paste0("./",opt$g,"/qc/Heatmap_PC.19.27.png"),1000,1000,pointsize=20) +png(paste0("./analysis/",opt$g,"/qc/Heatmap_PC.19.27.png"),1000,1000,pointsize=20) PCHeatmap(object=sc10x.Group,pc.use=19:27,cells.use=500,do.balanced=TRUE,label.columns=FALSE,use.full=FALSE) dev.off() -png(paste0("./",opt$g,"/qc/Heatmap_PC.28.36.png"),1000,1000,pointsize=20) +png(paste0("./analysis/",opt$g,"/qc/Heatmap_PC.28.36.png"),1000,1000,pointsize=20) PCHeatmap(object=sc10x.Group,pc.use=28:36,cells.use=500,do.balanced=TRUE,label.columns=FALSE,use.full=FALSE) dev.off() -png(paste0("./",opt$g,"/qc/Heatmap_PC.37.45.png"),1000,1000,pointsize=20) +png(paste0("./analysis/",opt$g,"/qc/Heatmap_PC.37.45.png"),1000,1000,pointsize=20) PCHeatmap(object=sc10x.Group,pc.use=37:45,cells.use=500,do.balanced=TRUE,label.columns=FALSE,use.full=FALSE) dev.off() -png(paste0("./",opt$g,"/qc/Heatmap_PC.46.50.png"),1000,1000,pointsize=20) +png(paste0("./analysis/",opt$g,"/qc/Heatmap_PC.46.50.png"),1000,1000,pointsize=20) PCHeatmap(object=sc10x.Group,pc.use=46:50,cells.use=500,do.balanced=TRUE,label.columns=FALSE,use.full=FALSE) dev.off() @@ -157,6 +157,6 @@ opt.qc <- opt rm(opt) #save qc.Rda file,Cumulative.RData file -save(list=paste0("sc10x.",opt.qc$g,".qc"),file=paste0("./sc10x.",Project.Name,".",opt.qc$g,".qc.Rda")) +save(list=paste0("sc10x.",opt.qc$g,".qc"),file=paste0("./analysis/sc10x.",Project.Name,".",opt.qc$g,".qc.Rda")) rm(list=paste0("sc10x.",opt.qc$g,".qc")) -save.image(file=paste0("./sc10x.",Project.Name,".Cumulative.RData")) +save.image(file=paste0("./analysis/sc10x.",Project.Name,".Cumulative.RData")) diff --git a/r.scripts/sc_QuSAGE.Epi.R b/r.scripts/sc_QuSAGE.Epi.R index f62fd25..d93b26b 100755 --- a/r.scripts/sc_QuSAGE.Epi.R +++ b/r.scripts/sc_QuSAGE.Epi.R @@ -25,24 +25,24 @@ if (opt$st==TRUE){ } #load data -setwd(paste0("./",Project.Name)) -load(paste0("./sc10x.",Project.Name,".Cumulative.RData")) -load(paste0("./sc10x.",Project.Name,".",opt$g,".cluster",sub.file,".IDlineage.Rda")) +setwd("../") +load(paste0("./analysis/sc10x.",Project.Name,".Cumulative.RData")) +load(paste0("./analysis/sc10x.",Project.Name,".",opt$g,".cluster",sub.file,".IDlineage.Rda")) sc10x.Group <- get(paste0("sc10x.",opt$g,".cluster",sub.file,".IDlineage")) rm(list=paste0("sc10x.",opt$g,".cluster",sub.file,".IDlineage")) #create folder structure -if (!dir.exists(paste0("./",opt$g))){ - dir.create(paste0("./",opt$g)) +if (!dir.exists(paste0("./analysis/",opt$g))){ + dir.create(paste0("./analysis/",opt$g)) } -if (!dir.exists(paste0("./",opt$g,"/global"))){ - dir.create(paste0("./",opt$g,"/global")) +if (!dir.exists(paste0("./analysis/",opt$g,"/global"))){ + dir.create(paste0("./analysis/",opt$g,"/global")) } -if (!dir.exists(paste0("./",opt$g,"/global",sub.folder))){ - dir.create(paste0("./",opt$g,"/global",sub.folder)) +if (!dir.exists(paste0("./analysis/",opt$g,"/global",sub.folder))){ + dir.create(paste0("./analysis/",opt$g,"/global",sub.folder)) } -if (!dir.exists(paste0("./",opt$g,"/global",sub.folder,"/correlation"))){ - dir.create(paste0("./",opt$g,"/global",sub.folder,"/correlation")) +if (!dir.exists(paste0("./analysis/",opt$g,"/global",sub.folder,"/correlation"))){ + dir.create(paste0("./analysis/",opt$g,"/global",sub.folder,"/correlation")) } #load resolution @@ -78,15 +78,15 @@ for (i in 1:Number.Clusters){ cC[i] <- paste0("Cluster_",i,"-REST") } rm(labels) -gene.set2 <- read_delim("../DEG_BE_2FC.txt","\t",escape_double=FALSE,trim_ws=TRUE,col_names=FALSE) +gene.set2 <- read_delim("./genesets/DEG_BE_2FC.txt","\t",escape_double=FALSE,trim_ws=TRUE,col_names=FALSE) gene.set2 <- as.list(gene.set2) names(gene.set2) <- "BE" gene.set <- c(gene.set2) -gene.set2 <- read_delim("../DEG_LE_2FC.txt","\t",escape_double=FALSE,trim_ws=TRUE,col_names=FALSE) +gene.set2 <- read_delim("./genesets/DEG_LE_2FC.txt","\t",escape_double=FALSE,trim_ws=TRUE,col_names=FALSE) gene.set2 <- as.list(gene.set2) names(gene.set2) <- "LE" gene.set <- c(gene.set,gene.set2) -gene.set2 <- read_delim("../DEG_OE_2FC.txt","\t",escape_double=FALSE,trim_ws=TRUE,col_names=FALSE) +gene.set2 <- read_delim("./genesets/DEG_OE_2FC.txt","\t",escape_double=FALSE,trim_ws=TRUE,col_names=FALSE) gene.set2 <- as.list(gene.set2) names(gene.set2) <- "OE" gene.set <- c(gene.set,gene.set2) @@ -140,7 +140,7 @@ for (j in 1:3){ } else if (j==3) { pop <- "OE" } - png(paste0("./",opt$g,"/global",sub.folder,"/correlation/QuSAGE_",pop,".png"),1000,1000,pointsize=20) + png(paste0("./analysis/",opt$g,"/global",sub.folder,"/correlation/QuSAGE_",pop,".png"),1000,1000,pointsize=20) for (i in Epi){ qs <- get(paste0("qs.results.",i)) if (i==Epi[1]){ @@ -170,7 +170,7 @@ sc10x.Group@ident <- factor(sc10x.Group@ident,levels=c("BE","LE","OE","St")) sc10x.Group <- StashIdent(object=sc10x.Group,save.name="EpiPop") #generate EpiPop tSNE -png(paste0("./",opt$g,"/global",sub.folder,"/tSNE_EpiPop.png"),1000,1000,pointsize=20) +png(paste0("./analysis/",opt$g,"/global",sub.folder,"/tSNE_EpiPop.png"),1000,1000,pointsize=20) plot <- TSNEPlot(object=sc10x.Group,pt.size=5,do.label=TRUE,label.size=10,do.return=TRUE) 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))) @@ -255,7 +255,7 @@ rm(BE) rm(LE) rm(OE) rm(merge.cluster) -png(paste0("./",opt$g,"/global",sub.folder,"/tSNE_EpiSubPop.png"),1000,1000,pointsize=20) +png(paste0("./anaysis/",opt$g,"/global",sub.folder,"/tSNE_EpiSubPop.png"),1000,1000,pointsize=20) plot <- TSNEPlot(object=sc10x.Group,pt.size=5,do.label=TRUE,label.size=10,do.return=TRUE) 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))) @@ -285,10 +285,10 @@ opt.epipop <- opt rm(opt) #save QuSAGE.EpiPop.Rda,IDepi.Rda file,Cumulative.RData files -save(list=ls(pattern="^qs."),file=paste0("./sc10x.",Project.Name,".",opt.lineage$g,sub.file,".QuSAGE.EpiPop.Rda")) +save(list=ls(pattern="^qs."),file=paste0("./analysis/sc10x.",Project.Name,".",opt.lineage$g,sub.file,".QuSAGE.EpiPop.Rda")) rm(list=ls(pattern="^qs.")) -save(list=paste0("sc10x.",opt.epipop$g,".cluster",sub.file,".IDepi"),file=paste0("./sc10x.",Project.Name,".",opt.epipop$g,".cluster",sub.file,".IDepi.Rda")) +save(list=paste0("sc10x.",opt.epipop$g,".cluster",sub.file,".IDepi"),file=paste0("./analysis/sc10x.",Project.Name,".",opt.epipop$g,".cluster",sub.file,".IDepi.Rda")) rm(list=paste0("sc10x.",opt.epipop$g,".cluster",sub.file,".IDepi")) rm(sub.folder) rm(sub.file) -save.image(file=paste0("./sc10x.",Project.Name,".Cumulative.RData")) +save.image(file=paste0("./analysis/sc10x.",Project.Name,".Cumulative.RData")) diff --git a/r.scripts/sc_QuSAGE.LGEA.Epi.R b/r.scripts/sc_QuSAGE.LGEA.Epi.R index 382ff61..dc595a8 100755 --- a/r.scripts/sc_QuSAGE.LGEA.Epi.R +++ b/r.scripts/sc_QuSAGE.LGEA.Epi.R @@ -26,24 +26,24 @@ if (opt$st==TRUE){ } #load data -setwd(paste0("./",Project.Name)) -load(paste0("./sc10x.",Project.Name,".Cumulative.RData")) -load(paste0("./sc10x.",Project.Name,".",opt$g,".cluster",sub.file,".IDepi+st+ne.Rda")) +setwd("../") +load(paste0("./analysis/sc10x.",Project.Name,".Cumulative.RData")) +load(paste0("./analysis/sc10x.",Project.Name,".",opt$g,".cluster",sub.file,".IDepi+st+ne.Rda")) sc10x.Group <- get(paste0("sc10x.",opt$g,".cluster",sub.file,".IDepi+st+ne")) rm(list=paste0("sc10x.",opt$g,".cluster",sub.file,".IDepi+st+ne")) #create folder structure -if (!dir.exists(paste0("./",opt$g))){ - dir.create(paste0("./",opt$g)) +if (!dir.exists(paste0("./analysis/",opt$g))){ + dir.create(paste0("./analysis/",opt$g)) } -if (!dir.exists(paste0("./",opt$g,"/global"))){ - dir.create(paste0("./",opt$g,"/global")) +if (!dir.exists(paste0("./analysis/",opt$g,"/global"))){ + dir.create(paste0("./analysis/",opt$g,"/global")) } -if (!dir.exists(paste0("./",opt$g,"/global",sub.folder))){ - dir.create(paste0("./",opt$g,"/global",sub.folder)) +if (!dir.exists(paste0("./analysis/",opt$g,"/global",sub.folder))){ + dir.create(paste0("./analysis/",opt$g,"/global",sub.folder)) } -if (!dir.exists(paste0("./",opt$g,"/global",sub.folder,"/correlation"))){ - dir.create(paste0("./",opt$g,"/global",sub.folder,"/correlation")) +if (!dir.exists(paste0("./analysis/",opt$g,"/global",sub.folder,"/correlation"))){ + dir.create(paste0("./analysis/",opt$g,"/global",sub.folder,"/correlation")) } #load resolution @@ -79,17 +79,17 @@ for (i in 1:Number.Clusters){ cC[i] <- paste0("Cluster_",i,"-REST") } rm(labels) -gene.set2 <- read_csv("../Basal cells-signature-genes.csv") +gene.set2 <- read_csv("./genesets/Basal cells-signature-genes.csv") gene.set2 <- gene.set2[,2] gene.set2 <- as.list(gene.set2) names(gene.set2) <- "hBC" gene.set <- c(gene.set2) -gene.set2 <- read_csv("../Normal AT2 cells-signature-genes.csv") +gene.set2 <- read_csv("./genesets/Normal AT2 cells-signature-genes.csv") gene.set2 <- gene.set2[,2] gene.set2 <- as.list(gene.set2) names(gene.set2) <- "hAT2" gene.set <- c(gene.set,gene.set2) -gene.set2 <- read_csv("../Club_Goblet cells-signature-genes.csv") +gene.set2 <- read_csv("./genesets/Club_Goblet cells-signature-genes.csv") gene.set2 <- gene.set2[,2] gene.set2 <- as.list(gene.set2) names(gene.set2) <- "hCGC" @@ -144,7 +144,7 @@ for (j in 1:3){ } else if (j==3) { pop <- "hCGC" } - png(paste0("./",opt$g,"/global",sub.folder,"/correlation/QuSAGE_",pop,".png"),1000,1000,pointsize=20) + png(paste0("./genesets/",opt$g,"/global",sub.folder,"/correlation/QuSAGE_",pop,".png"),1000,1000,pointsize=20) for (i in Epi){ qs <- get(paste0("qs.results.",i)) if (i==Epi[1]){ @@ -177,7 +177,7 @@ sc10x.Group@ident <- plyr::mapvalues(x=sc10x.Group@ident,from=Epi,to=merge.clust sc10x.Group <- StashIdent(object=sc10x.Group,save.name="LGEA.Epi") #generate EpiPop tSNE -png(paste0("./",opt$g,"/global",sub.folder,"/tSNE_LGEA.Epi.png"),1000,1000,pointsize=20) +png(paste0("./genesets/",opt$g,"/global",sub.folder,"/tSNE_LGEA.Epi.png"),1000,1000,pointsize=20) plot <- TSNEPlot(object=sc10x.Group,pt.size=5,do.label=TRUE,label.size=10,do.return=TRUE) 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))) @@ -256,7 +256,7 @@ for (i in 1:Number.Clusters){ }}}}} sc10x.Group@ident <- plyr::mapvalues(x=sc10x.Group@ident,from=1:Number.Clusters,to=merge.cluster) sc10x.Group <- StashIdent(object=sc10x.Group,save.name="LGEA.EpiSubPop") -png(paste0("./",opt$g,"/global",sub.folder,"/tSNE_LGEA.EpiSubPop.png"),1000,1000,pointsize=20) +png(paste0("./genesets/",opt$g,"/global",sub.folder,"/tSNE_LGEA.EpiSubPop.png"),1000,1000,pointsize=20) plot <- TSNEPlot(object=sc10x.Group,pt.size=5,do.label=TRUE,label.size=10,do.return=TRUE) 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))) @@ -293,10 +293,10 @@ opt.lgeaepipop <- opt rm(opt) #save QuSAGE.EpiPop.Rda,IDepi.Rda file,Cumulative.RData files -save(list=ls(pattern="^qs."),file=paste0("./sc10x.",Project.Name,".",opt.lineage$g,sub.file,".QuSAGE.LGEA.EpiPop.Rda")) +save(list=ls(pattern="^qs."),file=paste0("./analysis/sc10x.",Project.Name,".",opt.lineage$g,sub.file,".QuSAGE.LGEA.EpiPop.Rda")) rm(list=ls(pattern="^qs.")) -save(list=paste0("sc10x.",opt.lgeaepipop$g,".cluster",sub.file,".IDepi+st+ne+LGEAepi"),file=paste0("./sc10x.",Project.Name,".",opt.lgeaepipop$g,".cluster",sub.file,".IDepi+st+ne+LGEAepi.Rda")) +save(list=paste0("sc10x.",opt.lgeaepipop$g,".cluster",sub.file,".IDepi+st+ne+LGEAepi"),file=paste0("./analysis/sc10x.",Project.Name,".",opt.lgeaepipop$g,".cluster",sub.file,".IDepi+st+ne+LGEAepi.Rda")) rm(list=paste0("sc10x.",opt.lgeaepipop$g,".cluster",sub.file,".IDepi+st+ne+LGEAepi")) rm(sub.folder) rm(sub.file) -save.image(file=paste0("./sc10x.",Project.Name,".Cumulative.RData")) +save.image(file=paste0("./analysis/sc10x.",Project.Name,".Cumulative.RData")) diff --git a/r.scripts/sc_QuSAGE.LGEA.St.R b/r.scripts/sc_QuSAGE.LGEA.St.R index 9bb38f9..f684d6a 100755 --- a/r.scripts/sc_QuSAGE.LGEA.St.R +++ b/r.scripts/sc_QuSAGE.LGEA.St.R @@ -26,18 +26,18 @@ if (opt$st==TRUE){ } #load data -setwd(paste0("./",Project.Name)) -load(paste0("./sc10x.",Project.Name,".Cumulative.RData")) -load(paste0("./sc10x.",Project.Name,".",opt$g,".cluster",sub.file,".IDepi+st+ne+LGEAepi.Rda")) +setwd("../") +load(paste0("./analysis/sc10x.",Project.Name,".Cumulative.RData")) +load(paste0("./analysis/sc10x.",Project.Name,".",opt$g,".cluster",sub.file,".IDepi+st+ne+LGEAepi.Rda")) sc10x.Group <- get(paste0("sc10x.",opt$g,".cluster",sub.file,".IDepi+st+ne+LGEAepi")) rm(list=paste0("sc10x.",opt$g,".cluster",sub.file,".IDepi+st+ne+LGEAepi")) #create folder structure -if (!dir.exists(paste0("./",opt$g))){ - dir.create(paste0("./",opt$g)) +if (!dir.exists(paste0("./anaysis/",opt$g))){ + dir.create(paste0("./analysis/",opt$g)) } -if (!dir.exists(paste0("./",opt$g,"/global",sub.folder,"/correlation"))){ - dir.create(paste0("./",opt$g,"/global",sub.folder,"/correlation")) +if (!dir.exists(paste0("./analysis/",opt$g,"/global",sub.folder,"/correlation"))){ + dir.create(paste0("./analysis/",opt$g,"/global",sub.folder,"/correlation")) } #load resolution @@ -78,7 +78,7 @@ for (i in 1:Number.Clusters){ cC[i] <- paste0("Cluster_",i,"-REST") } rm(labels) -gene.set3 <- read_excel("../journal.pcbi.1004575.s026.XLSX",skip=2) +gene.set3 <- read_excel("./genesets/journal.pcbi.1004575.s026.XLSX",skip=2) gene.set3 <- gene.set3[,c(3,9)] gene.set2 <- gene.set3[gene.set3$`CELL TYPE`=="C1 : Proliferative Fibroblast",] gene.set2 <- as.list(gene.set2[!(is.na(gene.set2$`hu Gene Symbol`)),1]) @@ -161,7 +161,7 @@ for (j in 1:6){ } else if (j==6) { pop <- "MylImm" } - png(paste0("./",opt$g,"/global",sub.folder,"/correlation/QuSAGE_",pop,".png"),1000,1000,pointsize=20) + png(paste0("./analysis/",opt$g,"/global",sub.folder,"/correlation/QuSAGE_",pop,".png"),1000,1000,pointsize=20) for (i in St){ qs <- get(paste0("qs.results.",i)) if (i==St[1]){ @@ -209,7 +209,7 @@ rm(MylImm.cells) rm(UK.cells) #generate LGEA StPops tSNE -png(paste0("./",opt$g,"/global",sub.folder,"/tSNE_LGEA.St.png"),1000,1000,pointsize=20) +png(paste0("./"analysis/,opt$g,"/global",sub.folder,"/tSNE_LGEA.St.png"),1000,1000,pointsize=20) plot <- TSNEPlot(object=sc10x.Group,pt.size=5,do.label=TRUE,label.size=10,do.return=TRUE) 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))) @@ -327,7 +327,7 @@ for (i in 1:Number.Clusters){ }}}}} sc10x.Group@ident <- plyr::mapvalues(x=sc10x.Group@ident,from=1:Number.Clusters,to=merge.cluster) sc10x.Group <- StashIdent(object=sc10x.Group,save.name="LGES.StSubPop") -png(paste0("./",opt$g,"/global",sub.folder,"/tSNE_LGES.StSubPop.png"),1000,1000,pointsize=20) +png(paste0("./analysis/",opt$g,"/global",sub.folder,"/tSNE_LGES.StSubPop.png"),1000,1000,pointsize=20) plot <- TSNEPlot(object=sc10x.Group,pt.size=5,do.label=TRUE,label.size=10,do.return=TRUE) 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))) @@ -364,10 +364,10 @@ opt.lgeastpop <- opt rm(opt) #save QuSAGE.EpiPop.Rda,IDepi.Rda file,Cumulative.RData files -save(list=ls(pattern="^qs."),file=paste0("./sc10x.",Project.Name,".",opt.lineage$g,sub.file,".QuSAGE.LGEA.StPop.Rda")) +save(list=ls(pattern="^qs."),file=paste0("./analysis/sc10x.",Project.Name,".",opt.lineage$g,sub.file,".QuSAGE.LGEA.StPop.Rda")) rm(list=ls(pattern="^qs.")) -save(list=paste0("sc10x.",opt.lgeaepipop$g,".cluster",sub.file,".IDepi+st+ne+LGEAepi+LGEAst"),file=paste0("./sc10x.",Project.Name,".",opt.lgeaepipop$g,".cluster",sub.file,".IDepi+st+ne+LGEAepi+LGEAst.Rda")) +save(list=paste0("sc10x.",opt.lgeaepipop$g,".cluster",sub.file,".IDepi+st+ne+LGEAepi+LGEAst"),file=paste0("./analysis/sc10x.",Project.Name,".",opt.lgeaepipop$g,".cluster",sub.file,".IDepi+st+ne+LGEAepi+LGEAst.Rda")) rm(list=paste0("sc10x.",opt.epipop$g,".cluster",sub.file,".IDepi+st+ne+LGEAepi+LGEAst")) rm(sub.folder) rm(sub.file) -save.image(file=paste0("./sc10x.",Project.Name,".Cumulative.RData")) +save.image(file=paste0("./analysis/sc10x.",Project.Name,".Cumulative.RData")) diff --git a/r.scripts/sc_QuSAGE.Lineage.R b/r.scripts/sc_QuSAGE.Lineage.R index 5a2d306..5bcb1cb 100755 --- a/r.scripts/sc_QuSAGE.Lineage.R +++ b/r.scripts/sc_QuSAGE.Lineage.R @@ -25,24 +25,24 @@ if (opt$st==TRUE){ } #load data -setwd(paste0("./",Project.Name)) -load(paste0("./sc10x.",Project.Name,".Cumulative.RData")) -load(paste0("./sc10x.",Project.Name,".",opt$g,".cluster",sub.file,".Rda")) +setwd("../") +load(paste0("./analysis/sc10x.",Project.Name,".Cumulative.RData")) +load(paste0("./analysis/sc10x.",Project.Name,".",opt$g,".cluster",sub.file,".Rda")) sc10x.Group <- get(paste0("sc10x.",opt$g,".cluster",sub.file)) rm(list=paste0("sc10x.",opt$g,".cluster",sub.file)) #create folder structure -if (!dir.exists(paste0("./",opt$g))){ - dir.create(paste0("./",opt$g)) +if (!dir.exists(paste0("./analysis/",opt$g))){ + dir.create(paste0("./analysis/",opt$g)) } -if (!dir.exists(paste0("./",opt$g,"/global"))){ - dir.create(paste0("./",opt$g,"/global")) +if (!dir.exists(paste0("./analysis/",opt$g,"/global"))){ + dir.create(paste0("./analysis/",opt$g,"/global")) } -if (!dir.exists(paste0("./",opt$g,"/global",sub.folder))){ - dir.create(paste0("./",opt$g,"/global",sub.folder)) +if (!dir.exists(paste0("./analysis/",opt$g,"/global",sub.folder))){ + dir.create(paste0("./analysis/",opt$g,"/global",sub.folder)) } -if (!dir.exists(paste0("./",opt$g,"/global",sub.folder,"/correlation"))){ - dir.create(paste0("./",opt$g,"/global",sub.folder,"/correlation")) +if (!dir.exists(paste0("./analysis/",opt$g,"/global",sub.folder,"/correlation"))){ + dir.create(paste0("./analysis/",opt$g,"/global",sub.folder,"/correlation")) } #load resolution @@ -72,11 +72,11 @@ for (i in 1:Number.Clusters){ cC[i] <- paste0("Cluster_",i,"-REST") } rm(labels) -gene.set2 <- read_delim("../DEG_FMSt_2FC.txt","\t",escape_double=FALSE,trim_ws=TRUE,col_names=FALSE) +gene.set2 <- read_delim("./genesets/DEG_FMSt_2FC.txt","\t",escape_double=FALSE,trim_ws=TRUE,col_names=FALSE) gene.set2 <- as.list(gene.set2) names(gene.set2) <- "St" gene.set <- c(gene.set2) -gene.set2 <- read_delim("../DEG_Epi_2FC.txt","\t",escape_double=FALSE,trim_ws=TRUE,col_names=FALSE) +gene.set2 <- read_delim("./genesets/DEG_Epi_2FC.txt","\t",escape_double=FALSE,trim_ws=TRUE,col_names=FALSE) gene.set2 <- as.list(gene.set2) names(gene.set2) <- "Epi" gene.set <- c(gene.set,gene.set2) @@ -129,7 +129,7 @@ for (j in 1:2){ } else if (j==2) { pop <- "Epi" } - png(paste0("./",opt$g,"/global",sub.folder,"/correlation/QuSAGE_",pop,".png"),1000,1000,pointsize=20) + png(paste0("./analysis/",opt$g,"/global",sub.folder,"/correlation/QuSAGE_",pop,".png"),1000,1000,pointsize=20) for (i in 1:Number.Clusters){ qs <- get(paste0("qs.results.",i)) if (i==1){ @@ -156,7 +156,7 @@ sc10x.Group@ident <- factor(sc10x.Group@ident,levels=c("Epi","St")) sc10x.Group <- StashIdent(object=sc10x.Group,save.name="Lineage") #generate lineage tSNE -png(paste0("./",opt$g,"/global",sub.folder,"/tSNE_Lineage.png"),1000,1000,pointsize=20) +png(paste0("./analysis/",opt$g,"/global",sub.folder,"/tSNE_Lineage.png"),1000,1000,pointsize=20) plot <- TSNEPlot(object=sc10x.Group,pt.size=5,do.label=TRUE,label.size=10,do.return=TRUE) 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))) @@ -183,10 +183,10 @@ opt.lineage <- opt rm(opt) #save QuSAGE.Lineage.Rda,IDlineage.Rda file,Cumulative.RData files -save(list=ls(pattern="^qs."),file=paste0("./sc10x.",Project.Name,".",opt.lineage$g,sub.file,".QuSAGE.Lineage.Rda")) +save(list=ls(pattern="^qs."),file=paste0("./analysis/sc10x.",Project.Name,".",opt.lineage$g,sub.file,".QuSAGE.Lineage.Rda")) rm(list=ls(pattern="^qs.")) -save(list=paste0("sc10x.",opt.lineage$g,".cluster",sub.file,".IDlineage"),file=paste0("./sc10x.",Project.Name,".",opt.lineage$g,".cluster",sub.file,".IDlineage.Rda")) +save(list=paste0("sc10x.",opt.lineage$g,".cluster",sub.file,".IDlineage"),file=paste0("./analysis/sc10x.",Project.Name,".",opt.lineage$g,".cluster",sub.file,".IDlineage.Rda")) rm(list=paste0("sc10x.",opt.lineage$g,".cluster",sub.file,".IDlineage")) rm(sub.folder) rm(sub.file) -save.image(file=paste0("./sc10x.",Project.Name,".Cumulative.RData")) +save.image(file=paste0("./analysis/sc10x.",Project.Name,".Cumulative.RData")) diff --git a/r.scripts/sc_QuSAGE.St.R b/r.scripts/sc_QuSAGE.St.R index 0642644..6932444 100755 --- a/r.scripts/sc_QuSAGE.St.R +++ b/r.scripts/sc_QuSAGE.St.R @@ -24,18 +24,18 @@ if (opt$st==TRUE){ } #load data -setwd(paste0("./",Project.Name)) -load(paste0("./sc10x.",Project.Name,".Cumulative.RData")) -load(paste0("./sc10x.",Project.Name,".",opt$g,".cluster",sub.file,".IDepi.Rda")) +setwd("../") +load(paste0("./analysis/sc10x.",Project.Name,".Cumulative.RData")) +load(paste0("./analysis/sc10x.",Project.Name,".",opt$g,".cluster",sub.file,".IDepi.Rda")) sc10x.Group <- get(paste0("sc10x.",opt$g,".cluster",sub.file,".IDepi")) rm(list=paste0("sc10x.",opt$g,".cluster",sub.file,".IDepi")) #create folder structure -if (!dir.exists(paste0("./",opt$g))){ - dir.create(paste0("./",opt$g)) +if (!dir.exists(paste0("./analysis/",opt$g))){ + dir.create(paste0("./analysis/",opt$g)) } -if (!dir.exists(paste0("./",opt$g,"/global",sub.folder,"/correlation"))){ - dir.create(paste0("./",opt$g,"/global",sub.folder,"/correlation")) +if (!dir.exists(paste0("./analysis/",opt$g,"/global",sub.folder,"/correlation"))){ + dir.create(paste0("./analysis/",opt$g,"/global",sub.folder,"/correlation")) } #load resolution @@ -71,15 +71,15 @@ for (i in 1:Number.Clusters){ cC[i] <- paste0("Cluster_",i,"-REST") } rm(labels) -gene.set2 <- read_delim("../DEG_C5.BP.M11704.txt","\t",escape_double=FALSE,trim_ws=TRUE,col_names=TRUE) +gene.set2 <- read_delim("./genesets/DEG_C5.BP.M11704.txt","\t",escape_double=FALSE,trim_ws=TRUE,col_names=TRUE) gene.set2 <- as.list(gene.set2[-1,]) names(gene.set2) <- "Endo" gene.set <- c(gene.set2) -gene.set2 <- read_delim("../DEG_C5.BP.M10794.txt","\t",escape_double=FALSE,trim_ws=TRUE,col_names=TRUE) +gene.set2 <- read_delim("./genesets/DEG_C5.BP.M10794.txt","\t",escape_double=FALSE,trim_ws=TRUE,col_names=TRUE) gene.set2 <- as.list(gene.set2[-1,]) names(gene.set2) <- "SM" gene.set <- c(gene.set,gene.set2) -gene.set2 <- read_delim("../DEG_C5.BP.M13024.txt","\t",escape_double=FALSE,trim_ws=TRUE,col_names=TRUE) +gene.set2 <- read_delim("./genesets/DEG_C5.BP.M13024.txt","\t",escape_double=FALSE,trim_ws=TRUE,col_names=TRUE) gene.set2 <- as.list(gene.set2[-1,]) names(gene.set2) <- "Fib" gene.set <- c(gene.set,gene.set2) @@ -133,7 +133,7 @@ for (j in 1:3){ } else if (j==3) { pop <- "Fib" } - png(paste0("./",opt$g,"/global",sub.folder,"/correlation/QuSAGE_",pop,".png"),1000,1000,pointsize=20) + png(paste0("./analysis/",opt$g,"/global",sub.folder,"/correlation/QuSAGE_",pop,".png"),1000,1000,pointsize=20) for (i in St){ qs <- get(paste0("qs.results.",i)) if (i==St[1]){ @@ -173,7 +173,7 @@ rm(Fib.cells) rm(UK.cells) #generate ALLPop tSNE -png(paste0("./",opt$g,"/global",sub.folder,"/tSNE_ALLPop.png"),1000,1000,pointsize=20) +png(paste0("./analysis/",opt$g,"/global",sub.folder,"/tSNE_ALLPop.png"),1000,1000,pointsize=20) plot <- TSNEPlot(object=sc10x.Group,pt.size=5,do.label=TRUE,label.size=10,do.return=TRUE) 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))) @@ -291,7 +291,7 @@ for (i in 1:Number.Clusters){ }}}}} sc10x.Group@ident <- plyr::mapvalues(x=sc10x.Group@ident,from=1:Number.Clusters,to=merge.cluster) sc10x.Group <- StashIdent(object=sc10x.Group,save.name="ALLSubPop") -png(paste0("./",opt$g,"/global",sub.folder,"/tSNE_ALLSubPop.png"),1000,1000,pointsize=20) +png(paste0("./analysis/",opt$g,"/global",sub.folder,"/tSNE_ALLSubPop.png"),1000,1000,pointsize=20) plot <- TSNEPlot(object=sc10x.Group,pt.size=5,do.label=TRUE,label.size=10,do.return=TRUE) 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))) @@ -329,10 +329,10 @@ opt.stpop <- opt rm(opt) #save QuSAGE.EpiPop.Rda,IDepi.Rda file,Cumulative.RData files -save(list=ls(pattern="^qs."),file=paste0("./sc10x.",Project.Name,".",opt.lineage$g,sub.file,".QuSAGE.StPop.Rda")) +save(list=ls(pattern="^qs."),file=paste0("./analysis/sc10x.",Project.Name,".",opt.lineage$g,sub.file,".QuSAGE.StPop.Rda")) rm(list=ls(pattern="^qs.")) -save(list=paste0("sc10x.",opt.stpop$g,".cluster",sub.file,".IDepi+st"),file=paste0("./sc10x.",Project.Name,".",opt.stpop$g,".cluster",sub.file,".IDepi+st.Rda")) +save(list=paste0("sc10x.",opt.stpop$g,".cluster",sub.file,".IDepi+st"),file=paste0("./analysis/sc10x.",Project.Name,".",opt.stpop$g,".cluster",sub.file,".IDepi+st.Rda")) rm(list=paste0("sc10x.",opt.stpop$g,".cluster",sub.file,".IDepi+st")) rm(sub.folder) rm(sub.file) -save.image(file=paste0("./sc10x.",Project.Name,".Cumulative.RData")) +save.image(file=paste0("./analysis/sc10x.",Project.Name,".Cumulative.RData")) diff --git a/r.scripts/sc_Seurat.Score.CellCycle.R b/r.scripts/sc_Seurat.Score.CellCycle.R index 4e7ea82..628e30f 100755 --- a/r.scripts/sc_Seurat.Score.CellCycle.R +++ b/r.scripts/sc_Seurat.Score.CellCycle.R @@ -14,37 +14,37 @@ rm(option_list) Project.Name <- opt$p #load data -setwd(paste0("./",Project.Name)) -load(paste0("./sc10x.",Project.Name,".Cumulative.RData")) -load(paste0("./sc10x.",Project.Name,".",opt$g,".raw.Rda")) +setwd("../") +load(paste0("./analysis/sc10x.",Project.Name,".Cumulative.RData")) +load(paste0("./analysis/sc10x.",Project.Name,".",opt$g,".raw.Rda")) sc10x.Group <- get(paste0("sc10x.",opt$g,".raw")) rm(list=paste0("sc10x.",opt$g,".raw")) #create folder structure -if (!dir.exists(paste0("./",opt$g))){ - dir.create(paste0("./",opt$g)) +if (!dir.exists(paste0("./analysis/",opt$g))){ + dir.create(paste0("./analysis/",opt$g)) } -if (!dir.exists(paste0("./",opt$g,"/cc"))){ - dir.create(paste0("./",opt$g,"/cc")) +if (!dir.exists(paste0("./analysis/",opt$g,"/cc"))){ + dir.create(paste0("./analysis/",opt$g,"/cc")) } #score cell cycle and generate figures -cc.genes <- readLines(con="../regev_lab_cell_cycle_genes.txt") +cc.genes <- readLines(con="./genesets/regev_lab_cell_cycle_genes.txt") s.genes <- cc.genes[1:43] g2m.genes <- cc.genes[44:97] sc10x.Group <- NormalizeData(object=sc10x.Group) sc10x.Group <- ScaleData(object=sc10x.Group) sc10x.Group <- CellCycleScoring(object=sc10x.Group,s.genes=s.genes,g2m.genes=g2m.genes,set.ident=TRUE) -png(paste0("./",opt$g,"/cc/Ridge_cc.Raw.png"),1000,1000,pointsize=20) +png(paste0("./analysis/",opt$g,"/cc/Ridge_cc.Raw.png"),1000,1000,pointsize=20) plot <- RidgePlot(object=sc10x.Group,features.plot=c("PCNA","TOP2A","MCM6","MKI67"),y.log=TRUE,nCol=2,do.return=TRUE) plot(plot) dev.off() -png(paste0("./",opt$g,"/cc/Violin_cc.Raw.png"),1000,1000,pointsize=20) +png(paste0("./analysis/",opt$g,"/cc/Violin_cc.Raw.png"),1000,1000,pointsize=20) plot <- VlnPlot(object=sc10x.Group,features.plot=c("PCNA","TOP2A","MCM6","MKI67"),nCol=2,do.return=TRUE) plot(plot) dev.off() sc10x.Group <- RunPCA(object=sc10x.Group,pc.genes=c(s.genes,g2m.genes),do.print=FALSE,pcs.store=2) -png(paste0("./",opt$g,"/cc/PCA_cc.Raw.png"),1000,1000,pointsize=20) +png(paste0("./analysis/",opt$g,"/cc/PCA_cc.Raw.png"),1000,1000,pointsize=20) plot <- PCAPlot(object=sc10x.Group,do.return=TRUE) 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))) @@ -61,6 +61,6 @@ opt.cc <- opt rm(opt) #save qc.Rda file,Cumulative.RData file -save(list=paste0("sc10x.",opt.cc$g,".cc"),file=paste0("./sc10x.",Project.Name,".",opt.cc$g,".cc.Rda")) +save(list=paste0("sc10x.",opt.cc$g,".cc"),file=paste0("./analysis/sc10x.",Project.Name,".",opt.cc$g,".cc.Rda")) rm(list=paste0("sc10x.",opt.cc$g,".cc")) -save.image(file=paste0("./sc10x.",Project.Name,".Cumulative.RData")) +save.image(file=paste0("./analysis/sc10x.",Project.Name,".Cumulative.RData")) diff --git a/r.scripts/sc_Tables.R b/r.scripts/sc_Tables.R index d9d01e7..0267229 100755 --- a/r.scripts/sc_Tables.R +++ b/r.scripts/sc_Tables.R @@ -21,20 +21,20 @@ if (opt$st==TRUE){ } #load data -setwd(paste0("./",Project.Name)) -load(paste0("./sc10x.",Project.Name,".Cumulative.RData")) -load(paste0("./sc10x.",Project.Name,".",opt$g,".cluster",sub.file,".IDepi+st+ne+LGEAepi+LGEAst.Rda")) +setwd("../") +load(paste0("./analysis/sc10x.",Project.Name,".Cumulative.RData")) +load(paste0("./analysis/sc10x.",Project.Name,".",opt$g,".cluster",sub.file,".IDepi+st+ne+LGEAepi+LGEAst.Rda")) sc10x.Group <- get(paste0("sc10x.",opt$g,".cluster",sub.file,".IDepi+st+ne+LGEAepi+LGEAst")) rm(list=paste0("sc10x.",opt$g,".cluster",sub.file,".IDepi+st+ne+LGEAepi+LGEAst")) #create folder structure -if (!dir.exists(paste0("./",opt$g))){ - dir.create(paste0("./",opt$g)) +if (!dir.exists(paste0("./analysis/",opt$g))){ + dir.create(paste0("./analysis/",opt$g)) } -if (!dir.exists(paste0("./",opt$g,"/table"))){ - dir.create(paste0("./",opt$g,"/table")) +if (!dir.exists(paste0("./analysis/",opt$g,"/table"))){ + dir.create(paste0("./analysis/",opt$g,"/table")) } -write.table(table(sc10x.Group@meta.data$samples,sc10x.Group@meta.data$'ALLSubPop+NE')[5:6,],file=paste0("./",opt$g,"/table/Table.ALLSubPop+NE.csv"),row.names=TRUE,col.names=NA,append=FALSE,sep=",") -write.table(round(prop.table(table(sc10x.Group@meta.data$samples,sc10x.Group@meta.data$'ALLSubPop+NE')[5:6,],1)*100,1),file=paste0("./",opt$g,"/table/ProbTable.ALLSubPop+NE.csv"),row.names=TRUE,col.names=NA,append=FALSE,sep=",") -write.table(round(prop.table(table(sc10x.Group@meta.data$samples,sc10x.Group@meta.data$'ALLSubPop+NE')[1:4,],1)*100,1),file=paste0("./",opt$g,"/table/ProbTable.FACS.csv"),row.names=TRUE,col.names=NA,append=FALSE,sep=",") +write.table(table(sc10x.Group@meta.data$samples,sc10x.Group@meta.data$'ALLSubPop+NE')[5:6,],file=paste0("./analysis/",opt$g,"/table/Table.ALLSubPop+NE.csv"),row.names=TRUE,col.names=NA,append=FALSE,sep=",") +write.table(round(prop.table(table(sc10x.Group@meta.data$samples,sc10x.Group@meta.data$'ALLSubPop+NE')[5:6,],1)*100,1),file=paste0("./analysis/",opt$g,"/table/ProbTable.ALLSubPop+NE.csv"),row.names=TRUE,col.names=NA,append=FALSE,sep=",") +write.table(round(prop.table(table(sc10x.Group@meta.data$samples,sc10x.Group@meta.data$'ALLSubPop+NE')[1:4,],1)*100,1),file=paste0("./analysis/",opt$g,"/table/ProbTable.FACS.csv"),row.names=TRUE,col.names=NA,append=FALSE,sep=",") -- GitLab