Add DIY script for CCA aggregation of D17 and D27 and increase resolution of downstream analysis

bash.scripts/sc_FullAnalysis-DPrF.sh

+#!/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"

r.scripts/sc_DIY.DPrF.aggr.R

+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"))