Skip to content
Snippets Groups Projects
Commit 81f06892 authored by Gervaise Henry's avatar Gervaise Henry :cowboy:
Browse files

Add custom var.gene and %mito cutoffs and PC choices

parent 6a0e0cb8
Branches
Tags
No related merge requests found
Hide whitespace changes
Inline Side-by-side
bash.scripts/sc_FullAnalysis-D17PrF.sh
+ 1
− 1
View file @ 81f06892
...@@ -14,7 +14,7 @@ Rscript ../r.scripts/sc_Demultiplex.R --p="D17PrF" --d=3 ...@@ -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_Seurat.Score.CellCycle.R --p="D17PrF"
Rscript ../r.scripts/sc_QC.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_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_QuSAGE.Lineage.R --p="D17PrF" --r=0.1 --s=NA
Rscript ../r.scripts/sc_LineageSubClust.R --p="D17PrF" --pc=0 Rscript ../r.scripts/sc_LineageSubClust.R --p="D17PrF" --pc=0
Rscript ../r.scripts/sc_QuSAGE_EpiSubClust.R --p="D17PrF" --r=0.1 --s=NA Rscript ../r.scripts/sc_QuSAGE_EpiSubClust.R --p="D17PrF" --r=0.1 --s=NA
......
r.scripts/sc_Cluster.R
+ 1
− 1
View file @ 81f06892
...@@ -25,7 +25,7 @@ rm(list=paste0("sc10x.",opt$g,".qc")) ...@@ -25,7 +25,7 @@ rm(list=paste0("sc10x.",opt$g,".qc"))
if (opt$pc>0){ if (opt$pc>0){
PC.Max <- opt$pc PC.Max <- opt$pc
} else { } else {
PC.Max <- PC.var.8 PC.Max <- PC.var.75
} }
#create folder structure #create folder structure
......
r.scripts/sc_LineageSubClust.R
+ 12
− 10
View file @ 81f06892
...@@ -53,39 +53,41 @@ tryCatch({ ...@@ -53,39 +53,41 @@ tryCatch({
rm(sc10x.Group) rm(sc10x.Group)
#re-detect variable genes in Epi + St subsets #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.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.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.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 #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.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=25) sc10x.Group.St <- RunPCA(object=sc10x.Group.St,pc.genes=sc10x.Group.St@var.genes,do.print=FALSE,pcs.compute=50)
#generate QC figures #generate QC figures
postscript(paste0("./analysis/",opt$g,"/qc/Plot_PCElbow.Epi.eps")) 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+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) plot(plot)
dev.off() dev.off()
PC.var.epi <- (sc10x.Group.Epi@dr$pca@sdev)^2 PC.var.epi <- (sc10x.Group.Epi@dr$pca@sdev)^2
PC.var.epi <- PC.var.epi/sum(PC.var.epi) PC.var.epi <- PC.var.epi/sum(PC.var.epi)
PC.var.epi <- cumsum(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")) 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+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) plot(plot)
dev.off() dev.off()
PC.var.st <- (sc10x.Group.St@dr$pca@sdev)^2 PC.var.st <- (sc10x.Group.St@dr$pca@sdev)^2
PC.var.st <- PC.var.st/sum(PC.var.st) PC.var.st <- PC.var.st/sum(PC.var.st)
PC.var.st <- cumsum(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){ if (opt$pc>0){
PC.Max.epi <- opt$pc PC.Max.epi <- opt$pc
PC.Max.st <- opt$pc PC.Max.st <- opt$pc
} else { } else {
PC.Max.epi <- PC.var.epi.8 PC.Max.epi <- PC.var.epi.75
PC.Max.st <- PC.var.st.8 PC.Max.st <- PC.var.st.85
} }
#re-generate tSNE dimensions and generage sample figures #re-generate tSNE dimensions and generage sample figures
......
r.scripts/sc_PC.Score.Stress.R
+ 4
− 4
View file @ 81f06892
...@@ -163,16 +163,16 @@ sc10x.Group.subset.FilteredCellCount <- length(sc10x.Group.subset@data@Dimnames[ ...@@ -163,16 +163,16 @@ sc10x.Group.subset.FilteredCellCount <- length(sc10x.Group.subset@data@Dimnames[
rm(merge.cluster) rm(merge.cluster)
#regenerate tSNEs for range of resolutions #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 <- 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=25) 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 <- (sc10x.Group.subset@dr$pca@sdev)^2
PC.var.NoStress <- PC.var.NoStress/sum(PC.var.NoStress) PC.var.NoStress <- PC.var.NoStress/sum(PC.var.NoStress)
PC.var.NoStress <- cumsum(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){ if (opt$pc>0){
PC.Max <- opt$pc PC.Max <- opt$pc
} else { } 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) 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")) postscript(paste0("./analysis/",opt$g,"/global/NoStress/tSNE_Sample.eps"))
......
r.scripts/sc_QC.R
+ 16
− 15
View file @ 81f06892
...@@ -12,7 +12,7 @@ option_list=list( ...@@ -12,7 +12,7 @@ option_list=list(
make_option("--lg",action="store",default=500,type='integer',help="Threshold Cells With Minimum Genes"), 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("--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("--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?") make_option("--cc",action="store",default=TRUE,type='logical',help="Analyze Data With Cell Cycle ID?")
) )
opt=parse_args(OptionParser(option_list=option_list)) opt=parse_args(OptionParser(option_list=option_list))
...@@ -73,7 +73,8 @@ dev.off() ...@@ -73,7 +73,8 @@ dev.off()
sc10x.Group.FilteredCellCount <- length(sc10x.Group@data@Dimnames[[2]]) sc10x.Group.FilteredCellCount <- length(sc10x.Group@data@Dimnames[[2]])
sc10x.Group.FilteredGeneCount <- length(sc10x.Group@data@Dimnames[[1]]) sc10x.Group.FilteredGeneCount <- length(sc10x.Group@data@Dimnames[[1]])
sc10x.Group@data <- sc10x.Group@data[!rownames(sc10x.Group@data) %in% mito.genes,] 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){ if (opt$cc==TRUE){
gc() 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) 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){ ...@@ -92,35 +93,35 @@ if (opt$cc==TRUE){
gc() gc()
} }
sc10x.Group <- SetAllIdent(object=sc10x.Group,id="ALL") 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) sc10x.Group <- ProjectPCA(object=sc10x.Group,do.print=FALSE,pcs.store=25)
rm(mito.genes) rm(mito.genes)
rm(percent.mito) rm(percent.mito)
#generate qc figures #generate qc figures
postscript(paste0("./analysis/",opt$g,"/qc/Plot_PCElbow.eps")) 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+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) plot(plot)
dev.off() dev.off()
PC.var <- (sc10x.Group@dr$pca@sdev)^2 PC.var <- (sc10x.Group@dr$pca@sdev)^2
PC.var <- PC.var/sum(PC.var) PC.var <- PC.var/sum(PC.var)
PC.var <- cumsum(PC.var) PC.var <- cumsum(PC.var)
PC.var.8 <- min(which(PC.var>=0.8)) PC.var.75 <- min(which(PC.var>=0.75))
postscript(paste0("./analysis/",opt$g,"/qc/Plot_VizPCA.01.05.eps")) postscript(paste0("./analysis/",opt$g,"/qc/Plot_VizPCA.01.10.eps"))
VizPCA(object=sc10x.Group,pcs.use=1:5,num.genes=50,nCol=5) VizPCA(object=sc10x.Group,pcs.use=1:10,num.genes=50,nCol=5)
dev.off() dev.off()
postscript(paste0("./analysis/",opt$g,"/qc/Plot_VizPCA.06.10.eps")) postscript(paste0("./analysis/",opt$g,"/qc/Plot_VizPCA.11.20.eps"))
VizPCA(object=sc10x.Group,pcs.use=6:10,num.genes=50,nCol=5) VizPCA(object=sc10x.Group,pcs.use=11:20,num.genes=50,nCol=5)
dev.off() dev.off()
postscript(paste0("./analysis/",opt$g,"/qc/Plot_VizPCA.11.15.eps")) postscript(paste0("./analysis/",opt$g,"/qc/Plot_VizPCA.21.30.eps"))
VizPCA(object=sc10x.Group,pcs.use=11:15,num.genes=50,nCol=5) VizPCA(object=sc10x.Group,pcs.use=21:30,num.genes=50,nCol=5)
dev.off() dev.off()
postscript(paste0("./analysis/",opt$g,"/qc/Plot_VizPCA.16.20.eps")) postscript(paste0("./analysis/",opt$g,"/qc/Plot_VizPCA.31.40.eps"))
VizPCA(object=sc10x.Group,pcs.use=16:20,num.genes=50,nCol=5) VizPCA(object=sc10x.Group,pcs.use=31:40,num.genes=50,nCol=5)
dev.off() dev.off()
postscript(paste0("./analysis/",opt$g,"/qc/Plot_VizPCA.21.25.eps")) postscript(paste0("./analysis/",opt$g,"/qc/Plot_VizPCA.41.50.eps"))
VizPCA(object=sc10x.Group,pcs.use=21:25,num.genes=50,nCol=5) VizPCA(object=sc10x.Group,pcs.use=41:50,num.genes=50,nCol=5)
dev.off() dev.off()
#make generic variables of run specific variables #make generic variables of run specific variables
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment