Tweak var.genes/PC/resolution selection

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" --pc=0
+Rscript ../r.scripts/sc_PC.Score.Stress.R --p="D17PrF" --r=0.5 --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

r.scripts/sc_LineageSubClust.R

@@ -53,8 +53,8 @@ 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.075,x.high.cutoff=3.5,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=3.5,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.15,x.high.cutoff=3.5,y.cutoff=0.75,do.plot=FALSE)
+sc10x.Group.St <- FindVariableGenes(object=sc10x.Group.St,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.NoStress.Epi <- sc10x.Group.Epi@var.genes
 genes.var.NoStress.St <- sc10x.Group.St@var.genes
 

r.scripts/sc_PC.Score.NE.R

@@ -89,16 +89,16 @@ dev.off()
 
 #find first turning point on a KDE of NE1
 #hist(GetCellEmbeddings(object=sc10x.Group.Epi,reduction.type="NE",dims.use=1),prob=TRUE)
-d1 <- density(GetCellEmbeddings(object=sc10x.Group.Epi,reduction.type="NE",dims.use=1),n=1000)
+d1 <- density(GetCellEmbeddings(object=sc10x.Group.Epi,reduction.type="NE",dims.use=1),n=50)
 #d2 <- density(GetCellEmbeddings(object=sc10x.Group.Epi,reduction.type="NE",dims.use=2))
-#lines(d1)
+#lines(d1,col="blue")
 ts_1 <- ts(d1$y)
 #ts_2 <- ts(d2$y)
 tp1 <- turnpoints(ts_1)
 #tp2 <- turnpoints(ts_2)
 pit1 <- min(d1$x[tp1$pits][d1$x[tp1$pits]>d1$x[tp1$peaks][d1$y[tp1$peaks]==max(d1$y[tp1$peaks])]])
 #pit2 <- d2$x[tp2$pits][length(d2$x[tp2$pits])]
-#abline(v=pit1)
+#abline(v=pit1,col="red")
 
 ALL <- GetCellEmbeddings(object=sc10x.Group.Epi,reduction.type="NE",dims.use=1:2)
 NE <- rownames(ALL[ALL[,1]>pit1,])

r.scripts/sc_PC.Score.Stress.R

@@ -76,7 +76,7 @@ dev.off()
 
 #find first turning point on a KDE of Stress1
 histo <- hist(GetCellEmbeddings(object=sc10x.Group,reduction.type="Stress",dims.use=1),breaks=10,prob=TRUE,plot=TRUE)
-d1 <- density(GetCellEmbeddings(object=sc10x.Group,reduction.type="Stress",dims.use=1),n=1000)
+d1 <- density(GetCellEmbeddings(object=sc10x.Group,reduction.type="Stress",dims.use=1),n=50)
 #d2 <- density(GetCellEmbeddings(object=sc10x.Group,reduction.type="Stress",dims.use=2),n=1000)
 ts_1 <- ts(d1$y)
 #ts_2 <- ts(d2$y)
@@ -99,8 +99,6 @@ plot <- plot+geom_vline(xintercept=pit1,color="red",lwd=2.5)
 plot(plot)
 dev.off()
 
-
-
 #subsample all cells (+Stress) to better visualize their clustering
 if (!is.na(opt$s)){
   rnd <- sample(1:ncol(sc10x.Group@data),opt$sv)
@@ -163,16 +161,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.075,x.high.cutoff=4,y.cutoff=0.075,do.plot=FALSE)
+sc10x.Group.subset <- FindVariableGenes(object=sc10x.Group.subset,mean.function=ExpMean,dispersion.function=LogVMR,x.low.cutoff=0.15,x.high.cutoff=3.5,y.cutoff=0.75,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.75 <- min(which(PC.var.NoStress>=0.75))
+PC.var.NoStress.95 <- min(which(PC.var.NoStress>=0.95))
 if (opt$pc>0){
   PC.Max <- opt$pc
 } else {
-  PC.Max <- PC.var.NoStress.75
+  PC.Max <- PC.var.NoStress.95
 }
 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"))

r.scripts/sc_QC.R

@@ -73,7 +73,7 @@ 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.075,x.high.cutoff=3.5,y.cutoff=0.5,do.plot=FALSE)
+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
 if (opt$cc==TRUE){
   gc()