diff --git a/bash.scripts/sc_FullAnalysis-DPrF.sh b/bash.scripts/sc_FullAnalysis-DPrF.sh new file mode 100644 index 0000000000000000000000000000000000000000..49f95538d3dd8510c0cf4984d16b4bd84cc77470 --- /dev/null +++ b/bash.scripts/sc_FullAnalysis-DPrF.sh @@ -0,0 +1,28 @@ +#!/bin/bash +#SBATCH --job-name R_FullAnalysis +#SBATCH -p 256GB,256GBv1,384GB +#SBATCH -N 1 +#SBATCH -t 7-0:0:0 +#SBATCH -o job_%j.out +#SBATCH -e job_%j.out +#SBATCH --mail-type ALL +#SBATCH --mail-user gervaise.henry@utsouthwestern.edu + +module load R/3.4.1-gccmkl + +#Rscript ../r.scripts/sc_Demultiplex.R --p="DPrF" --d=3 +#Rscript ../r.scripts/sc_Seurat.Score.CellCycle.R --p="DPrF" +#Rscript ../r.scripts/sc_QC.R --p="DPrF" + +Rscript ../r.scripts/sc_DIY.DPrF.aggr.R + +Rscript ../r.scripts/sc_Cluster.R --p="DPrF" --pc=0 +Rscript ../r.scripts/sc_PC.Score.Stress.R --p="DPrF" --r=1 --pc=0 --u=TRUE +Rscript ../r.scripts/sc_QuSAGE.Lineage.R --p="DPrF" --gs="sc" --r=0.5 --s=NA +Rscript ../r.scripts/sc_LineageSubClust.R --p="DPrF" --pc=0 +Rscript ../r.scripts/sc_QuSAGE_EpiSubClust.R --p="DPrF" --r=1 --s=NA --sc=TRUE +Rscript ../r.scripts/sc_QuSAGE_StSubClust.R --p="DPrF" --r=1 --s=NA --sc=TRUE +Rscript ../r.scripts/sc_MergeSubClust.R --p="DPrF" --e="sc" --s="sc" +Rscript ../r.scripts/sc_PC.Score.NE.R --p="DPrF" --u=TRUE +#Rscript ../r.scripts/sc_DEG.R --p="DPrF" +Rscript ../r.scripts/sc_Tables.R --p="DPrF" diff --git a/r.scripts/sc_DIY.DPrF.aggr.R b/r.scripts/sc_DIY.DPrF.aggr.R new file mode 100644 index 0000000000000000000000000000000000000000..c0c37b613be09026ef04a6948393a246a1bf8dcf --- /dev/null +++ b/r.scripts/sc_DIY.DPrF.aggr.R @@ -0,0 +1,101 @@ +gc() +.libPaths("/home2/ghenry/R/x86_64-pc-linux-gnu-library/3.4") +library(methods) +library(Seurat) + +load("../analysis/DATA/sc10x.D17PrF.ALL.qc.Rda") +sc10x.D17 <- sc10x.ALL.qc +load("../analysis/DATA/sc10x.D27PrF.ALL.qc.Rda") +sc10x.D27 <- sc10x.ALL.qc +rm(sc10x.ALL.qc) + +if (!dir.exists("../analysis/ALL")){ + dir.create("../analysis/ALL") +} +if (!dir.exists("../analysis/ALL/qc")){ + dir.create("../analysis/ALL/qc") +} +if (!dir.exists("../analysis/CCA")){ + dir.create("../analysis/CCA") +} + +sc10x.D17 <- FindVariableGenes(sc10x.D17,do.plot=FALSE) +sc10x.D27 <- FindVariableGenes(sc10x.D27,do.plot=FALSE) +var.genes.D17 <- head(rownames(sc10x.D17@hvg.info),1000) +var.genes.D27 <- head(rownames(sc10x.D27@hvg.info),1000) +var.genes.Comb <- unique(c(var.genes.D17,var.genes.D27)) +var.genes.Comb <- intersect(var.genes.Comb,rownames(sc10x.D17@scale.data)) +var.genes.Comb <- intersect(var.genes.Comb,rownames(sc10x.D27@scale.data)) + +sc10x.D17@meta.data$patient <- "D17" +sc10x.D27@meta.data$patient <- "D27" +sc10x.Group <- RunCCA(sc10x.D17,sc10x.D27,genes.use=var.genes.Comb,num.cc=50,add.cell.id1="D17",add.cell.id2="D27") +rm(sc10x.D17) +rm(sc10x.D27) + +postscript(paste0("../analysis/CCA/CCA.eps"),paper="special",width=10,height=5,horizontal=FALSE) +DimPlot(sc10x.Group,reduction.use="cca",group.by="patient",do.return=FALSE,pt.size=0.05) +dev.off() +postscript(paste0("../analysis/CCA/Violin.CCA.Raw.eps"),paper="special",width=10,height=5,horizontal=FALSE) +plot1 <- VlnPlot(sc10x.Group,features.plot="CC1",group.by="patient",do.return=TRUE,point.size.use=0.05) +plot2 <- VlnPlot(sc10x.Group,features.plot="CC2",group.by="patient",do.return=TRUE,point.size.use=0.05) +plot_grid(plot1,plot2) +dev.off() +rm(plot1) +rm(plot2) +postscript(paste0("../analysis/CCA/Bicor.eps"),paper="special",width=10,height=5,horizontal=FALSE) +MetageneBicorPlot(sc10x.Group,dims.eval=1:50,grouping.var="patient") +dev.off() + +gc() +sc10x.Group <- AlignSubspace(sc10x.Group,reduction.type="cca",grouping.var="patient",dims.align=1:30) +gc() + +postscript(paste0("../analysis/CCA/Violin.CCA.Aligned.eps"),paper="special",width=10,height=5,horizontal=FALSE) +plot1 <- VlnPlot(sc10x.Group,features.plot="ACC1",group.by="patient",do.return=TRUE,point.size.use=0.05) +plot2 <- VlnPlot(sc10x.Group,features.plot="ACC2",group.by="patient",do.return=TRUE,point.size.use=0.05) +plot_grid(plot1,plot2) +dev.off() +rm(plot1) +rm(plot2) + +sc10x.Group <- FindVariableGenes(object=sc10x.Group,mean.function=ExpMean,dispersion.function=LogVMR,x.low.cutoff=0.15,x.high.cutoff=3.5,y.cutoff=0.75,do.plot=FALSE) +genes.var <- sc10x.Group@var.genes +sc10x.Group <- SetAllIdent(object=sc10x.Group,id="ALL") +sc10x.Group <- RunPCA(object=sc10x.Group,pc.genes=sc10x.Group@var.genes,do.print=FALSE,pcs.compute=50) +sc10x.Group <- ProjectPCA(object=sc10x.Group,do.print=FALSE,pcs.store=25) + + +#generate qc figures +postscript("../analysis/ALL/qc/Plot_PCElbow.eps") +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() +PC.var <- (sc10x.Group@dr$pca@sdev)^2 +PC.var <- PC.var/sum(PC.var) +PC.var <- cumsum(PC.var) +PC.var.75 <- min(which(PC.var>=0.75)) +postscript("../analysis/ALL/qc/Plot_VizPCA.01.10.eps") +VizPCA(object=sc10x.Group,pcs.use=1:10,num.genes=50,nCol=5) +dev.off() +postscript("../analysis/ALL/qc/Plot_VizPCA.11.20.eps") +VizPCA(object=sc10x.Group,pcs.use=11:20,num.genes=50,nCol=5) +dev.off() +postscript("../analysis/ALL/qc/Plot_VizPCA.21.30.eps") +VizPCA(object=sc10x.Group,pcs.use=21:30,num.genes=50,nCol=5) +dev.off() +postscript("../analysis/ALL/qc/Plot_VizPCA.31.40.eps") +VizPCA(object=sc10x.Group,pcs.use=31:40,num.genes=50,nCol=5) +dev.off() +postscript("../analysis/ALL/qc/Plot_VizPCA.41.50.eps") +VizPCA(object=sc10x.Group,pcs.use=41:50,num.genes=50,nCol=5) +dev.off() + +#save qc.Rda file,Cumulative.RData file +assign("sc10x.ALL.qc",sc10x.Group) +rm(sc10x.Group) +save(sc10x.ALL.qc,file="../analysis/sc10x.DPrF.ALL.qc.Rda") +rm(sc10x.ALL.qc) +rm(plot) +save.image(file=paste0("../analysis/sc10x.DPrF.Cumulative.RData"))