diff --git a/bash.scripts/sc_FullAnalysis-D17PrF.sh b/bash.scripts/sc_FullAnalysis-D17PrF.sh
index 8200e96d107dcbb9605251311ed6dc4bbb8056d7..af7d83b9c17b95c2f47636ac3f42f6d81d86adfa 100755
--- a/bash.scripts/sc_FullAnalysis-D17PrF.sh
+++ b/bash.scripts/sc_FullAnalysis-D17PrF.sh
@@ -14,7 +14,7 @@ Rscript ../r.scripts/sc_Demultiplex.R --p="D17PrF" --d=3
Rscript ../r.scripts/sc_Seurat.Score.CellCycle.R --p="D17PrF"
Rscript ../r.scripts/sc_QC.R --p="D17PrF"
Rscript ../r.scripts/sc_Cluster.R --p="D17PrF" --pc=0
-Rscript ../r.scripts/sc_PC.Score.Stress.R --p="D17PrF" --r=0.5 --pc=0
+Rscript ../r.scripts/sc_PC.Score.Stress.R --p="D17PrF" --pc=0
Rscript ../r.scripts/sc_QuSAGE.Lineage.R --p="D17PrF" --r=0.1 --s=NA
Rscript ../r.scripts/sc_LineageSubClust.R --p="D17PrF" --pc=0
Rscript ../r.scripts/sc_QuSAGE_EpiSubClust.R --p="D17PrF" --r=0.1 --s=NA
diff --git a/r.scripts/sc_Cluster.R b/r.scripts/sc_Cluster.R
index 251a8dc3284c09e0c8d80e3bf9231232024b5f0d..5e13d530a886f38572c8527b8c74cdd9751033f8 100755
--- a/r.scripts/sc_Cluster.R
+++ b/r.scripts/sc_Cluster.R
@@ -25,7 +25,7 @@ rm(list=paste0("sc10x.",opt$g,".qc"))
if (opt$pc>0){
PC.Max <- opt$pc
} else {
- PC.Max <- PC.var.8
+ PC.Max <- PC.var.75
}
#create folder structure
diff --git a/r.scripts/sc_LineageSubClust.R b/r.scripts/sc_LineageSubClust.R
index 58fce37acc31a641496bb535b98a34272a870ee2..d217ca9d615561ac667753f9f85112fa0ebf700b 100755
--- a/r.scripts/sc_LineageSubClust.R
+++ b/r.scripts/sc_LineageSubClust.R
@@ -53,39 +53,41 @@ tryCatch({
rm(sc10x.Group)
#re-detect variable genes in Epi + St subsets
-sc10x.Group.Epi <- FindVariableGenes(object=sc10x.Group.Epi,mean.function=ExpMean,dispersion.function=LogVMR,x.low.cutoff=0.0125,x.high.cutoff=3,y.cutoff=0.5,do.plot=FALSE)
-sc10x.Group.St <- FindVariableGenes(object=sc10x.Group.St,mean.function=ExpMean,dispersion.function=LogVMR,x.low.cutoff=0.0125,x.high.cutoff=3,y.cutoff=0.5,do.plot=FALSE)
+sc10x.Group.Epi <- FindVariableGenes(object=sc10x.Group.Epi,mean.function=ExpMean,dispersion.function=LogVMR,x.low.cutoff=0.075,x.high.cutoff=4,y.cutoff=0.075,do.plot=FALSE)
+sc10x.Group.St <- FindVariableGenes(object=sc10x.Group.St,mean.function=ExpMean,dispersion.function=LogVMR,x.low.cutoff=0.075,x.high.cutoff=4,y.cutoff=0.075,do.plot=FALSE)
+genes.var.NoStress.Epi <- sc10x.Group.Epi@var.genes
+genes.var.NoStress.St <- sc10x.Group.St@var.genes
#re-PCA on new variable genes
-sc10x.Group.Epi <- RunPCA(object=sc10x.Group.Epi,pc.genes=sc10x.Group.Epi@var.genes,do.print=FALSE,pcs.compute=25)
-sc10x.Group.St <- RunPCA(object=sc10x.Group.St,pc.genes=sc10x.Group.St@var.genes,do.print=FALSE,pcs.compute=25)
+sc10x.Group.Epi <- RunPCA(object=sc10x.Group.Epi,pc.genes=sc10x.Group.Epi@var.genes,do.print=FALSE,pcs.compute=50)
+sc10x.Group.St <- RunPCA(object=sc10x.Group.St,pc.genes=sc10x.Group.St@var.genes,do.print=FALSE,pcs.compute=50)
#generate QC figures
postscript(paste0("./analysis/",opt$g,"/qc/Plot_PCElbow.Epi.eps"))
-plot <- PCElbowPlot(object=sc10x.Group.Epi,num.pc=25)
+plot <- PCElbowPlot(object=sc10x.Group.Epi,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.epi <- (sc10x.Group.Epi@dr$pca@sdev)^2
PC.var.epi <- PC.var.epi/sum(PC.var.epi)
PC.var.epi <- cumsum(PC.var.epi)
-PC.var.epi.8 <- min(which(PC.var.epi>=0.8))
+PC.var.epi.75 <- min(which(PC.var.epi>=0.75))
postscript(paste0("./analysis/",opt$g,"/qc/Plot_PCElbow.St.eps"))
-plot <- PCElbowPlot(object=sc10x.Group.St,num.pc=25)
+plot <- PCElbowPlot(object=sc10x.Group.St,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.st <- (sc10x.Group.St@dr$pca@sdev)^2
PC.var.st <- PC.var.st/sum(PC.var.st)
PC.var.st <- cumsum(PC.var.st)
-PC.var.st.8 <- min(which(PC.var.st>=0.8))
+PC.var.st.85 <- min(which(PC.var.st>=0.85))
if (opt$pc>0){
PC.Max.epi <- opt$pc
PC.Max.st <- opt$pc
} else {
- PC.Max.epi <- PC.var.epi.8
- PC.Max.st <- PC.var.st.8
+ PC.Max.epi <- PC.var.epi.75
+ PC.Max.st <- PC.var.st.85
}
#re-generate tSNE dimensions and generage sample figures
diff --git a/r.scripts/sc_PC.Score.Stress.R b/r.scripts/sc_PC.Score.Stress.R
index 176d89b583061195862719425c8acc6492f748fc..9ebbd2803331f5807073ae5996ba6fbebf494da4 100755
--- a/r.scripts/sc_PC.Score.Stress.R
+++ b/r.scripts/sc_PC.Score.Stress.R
@@ -163,16 +163,16 @@ sc10x.Group.subset.FilteredCellCount <- length(sc10x.Group.subset@data@Dimnames[
rm(merge.cluster)
#regenerate tSNEs for range of resolutions
-sc10x.Group.subset <- FindVariableGenes(object=sc10x.Group.subset,mean.function=ExpMean,dispersion.function=LogVMR,x.low.cutoff=0.0125,x.high.cutoff=3,y.cutoff=0.5,do.plot=FALSE)
-sc10x.Group.subset <- RunPCA(object=sc10x.Group.subset,pc.genes=sc10x.Group.subset@var.genes,do.print=FALSE,pcs.compute=25)
+sc10x.Group.subset <- FindVariableGenes(object=sc10x.Group.subset,mean.function=ExpMean,dispersion.function=LogVMR,x.low.cutoff=0.075,x.high.cutoff=4,y.cutoff=0.075,do.plot=FALSE)
+sc10x.Group.subset <- RunPCA(object=sc10x.Group.subset,pc.genes=sc10x.Group.subset@var.genes,do.print=FALSE,pcs.compute=50)
PC.var.NoStress <- (sc10x.Group.subset@dr$pca@sdev)^2
PC.var.NoStress <- PC.var.NoStress/sum(PC.var.NoStress)
PC.var.NoStress <- cumsum(PC.var.NoStress)
-PC.var.NoStress.8 <- min(which(PC.var.NoStress>=0.8))
+PC.var.NoStress.75 <- min(which(PC.var.NoStress>=0.75))
if (opt$pc>0){
PC.Max <- opt$pc
} else {
- PC.Max <- PC.var.NoStress.8
+ PC.Max <- PC.var.NoStress.75
}
sc10x.Group.subset <- RunTSNE(object=sc10x.Group.subset,dims.use=1:PC.Max,do.fast=TRUE,force.recalc=TRUE)
postscript(paste0("./analysis/",opt$g,"/global/NoStress/tSNE_Sample.eps"))
diff --git a/r.scripts/sc_QC.R b/r.scripts/sc_QC.R
index ca21d8df97f5b7aafa0c3278f0ecf811b7b5aa1c..a2c2722cc9ecb74bfd98f2af67a5a1f4e19f06cd 100755
--- a/r.scripts/sc_QC.R
+++ b/r.scripts/sc_QC.R
@@ -12,7 +12,7 @@ option_list=list(
make_option("--lg",action="store",default=500,type='integer',help="Threshold Cells With Minimum Genes"),
make_option("--hg",action="store",default=2500,type='integer',help="Threshold Cells With Maximum Genes"),
make_option("--lm",action="store",default=0,type='numeric',help="Threshold Cells With Minimum %Mitochodiral genes"),
- make_option("--hm",action="store",default=0.05,type='numeric',help="Threshold Cells With Maximum %Mitochodiral genes"),
+ make_option("--hm",action="store",default=0.1,type='numeric',help="Threshold Cells With Maximum %Mitochodiral genes"),
make_option("--cc",action="store",default=TRUE,type='logical',help="Analyze Data With Cell Cycle ID?")
)
opt=parse_args(OptionParser(option_list=option_list))
@@ -73,7 +73,8 @@ dev.off()
sc10x.Group.FilteredCellCount <- length(sc10x.Group@data@Dimnames[[2]])
sc10x.Group.FilteredGeneCount <- length(sc10x.Group@data@Dimnames[[1]])
sc10x.Group@data <- sc10x.Group@data[!rownames(sc10x.Group@data) %in% mito.genes,]
-sc10x.Group <- FindVariableGenes(object=sc10x.Group,mean.function=ExpMean,dispersion.function=LogVMR,x.low.cutoff=0.0125,x.high.cutoff=3,y.cutoff=0.5,do.plot=FALSE)
+sc10x.Group <- FindVariableGenes(object=sc10x.Group,mean.function=ExpMean,dispersion.function=LogVMR,x.low.cutoff=0.075,x.high.cutoff=4,y.cutoff=0.075,do.plot=FALSE)
+genes.var <- sc10x.Group@var.genes
if (opt$cc==TRUE){
gc()
sc10x.Group <- ScaleData(object=sc10x.Group,vars.to.regress=c("nUMI","percent.mito","S.Score","G2M.Score"),display.progress=FALSE,do.par=TRUE,num.cores=50)
@@ -92,35 +93,35 @@ if (opt$cc==TRUE){
gc()
}
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=25)
+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)
rm(mito.genes)
rm(percent.mito)
#generate qc figures
postscript(paste0("./analysis/",opt$g,"/qc/Plot_PCElbow.eps"))
-plot <- PCElbowPlot(object=sc10x.Group,num.pc=25)
+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.8 <- min(which(PC.var>=0.8))
-postscript(paste0("./analysis/",opt$g,"/qc/Plot_VizPCA.01.05.eps"))
-VizPCA(object=sc10x.Group,pcs.use=1:5,num.genes=50,nCol=5)
+PC.var.75 <- min(which(PC.var>=0.75))
+postscript(paste0("./analysis/",opt$g,"/qc/Plot_VizPCA.01.10.eps"))
+VizPCA(object=sc10x.Group,pcs.use=1:10,num.genes=50,nCol=5)
dev.off()
-postscript(paste0("./analysis/",opt$g,"/qc/Plot_VizPCA.06.10.eps"))
-VizPCA(object=sc10x.Group,pcs.use=6:10,num.genes=50,nCol=5)
+postscript(paste0("./analysis/",opt$g,"/qc/Plot_VizPCA.11.20.eps"))
+VizPCA(object=sc10x.Group,pcs.use=11:20,num.genes=50,nCol=5)
dev.off()
-postscript(paste0("./analysis/",opt$g,"/qc/Plot_VizPCA.11.15.eps"))
-VizPCA(object=sc10x.Group,pcs.use=11:15,num.genes=50,nCol=5)
+postscript(paste0("./analysis/",opt$g,"/qc/Plot_VizPCA.21.30.eps"))
+VizPCA(object=sc10x.Group,pcs.use=21:30,num.genes=50,nCol=5)
dev.off()
-postscript(paste0("./analysis/",opt$g,"/qc/Plot_VizPCA.16.20.eps"))
-VizPCA(object=sc10x.Group,pcs.use=16:20,num.genes=50,nCol=5)
+postscript(paste0("./analysis/",opt$g,"/qc/Plot_VizPCA.31.40.eps"))
+VizPCA(object=sc10x.Group,pcs.use=31:40,num.genes=50,nCol=5)
dev.off()
-postscript(paste0("./analysis/",opt$g,"/qc/Plot_VizPCA.21.25.eps"))
-VizPCA(object=sc10x.Group,pcs.use=21:25,num.genes=50,nCol=5)
+postscript(paste0("./analysis/",opt$g,"/qc/Plot_VizPCA.41.50.eps"))
+VizPCA(object=sc10x.Group,pcs.use=41:50,num.genes=50,nCol=5)
dev.off()
#make generic variables of run specific variables