Add hard cutoff for nCount and a bunch of small stuff

r.scripts/sc-TissueMapper_RAW.R

@@ -62,8 +62,14 @@ rm(results)
 #counts.gene.filtered.3gene <- results[[5]]
 #rm(results)
 
-counts.cell.filtered <- counts.cell.filtered.2mito
-counts.gene.filtered <- counts.gene.filtered.2mito
+sc10x <- lapply(sc10x,function(x) subset(x,subset= nCount_RNA>2500))
+counts.cell.filtered.3gene <- lapply(sc10x,ncol)
+counts.gene.filtered.3gene <- lapply(sc10x,nrow)
+
+counts.cell.filtered <- counts.cell.filtered.3gene
+counts.gene.filtered <- counts.gene.filtered.3gene
+
+save.image(file="./analysis/sc10x.raw_filtered.RData")
 
 sc10x <- sc10x[as.character(unlist(sc10x.groups[sc10x.groups$Keep==1,1]))]
 
@@ -137,23 +143,20 @@ if (length(sc10x)>1){
 }
 gc()
 
-if (ncol(sc10x) > 1000) {
-  pc.calc <- 1000
-} else if (ncol(sc10x) > 500) {
-  pc.calc <- 500
-} else if (ncol(sc10x) > 100) {
+save.image(file="./analysis/sc10x.raw_aligned.RData")
+
+if (ncol(sc10x) > 100) {
+  pc.calc <- 100
+} else if (ncol(sc10x) <= 100) {
   pc.calc <- ncol(sc10x)-1
 }
-results <- scPC(sc10x,pc=pc.calc,hpc=0.9,file="ALL",print="umap")
+results <- scPC(sc10x,pc=pc.calc,hpc=0.9,file=paste0("pre.stress"),print="umap",assay="integrated")
 sc10x <- results[[1]]
-pc.use.poststress <- results[[2]]
+pc.use.prestress <- results[[2]]
 rm(results)
+sc10x <- scCluster(sc10x,res=0.5,red="pca",dim=pc.use.prestress,print="umap",folder=paste0("preStress"),assay="integrated")
 rm(pc.calc)
 
-res <- c(seq(0.1,0.5,0.1),0.75,seq(1,5,1))
-
-sc10x <- scCluster(sc10x,res=res,red="pca",dim=pc.use.poststress,print="umap",folder="ALL")
-
 if (opt$s == "hu"){
   genes.stress <- read_delim("/work/urology/ghenry/RNA-Seq/SingleCell/ANALYSIS/REF/stress/genes.deg.Stress.csv",",",escape_double=FALSE,trim_ws=TRUE,col_names=TRUE)
   genes.stress <- genes.stress[1]
@@ -166,15 +169,36 @@ if (opt$s == "hu"){
   anchor <- c("Egr1","Fos","Jun")
 }
 
-results <- scScore(list(ALL=sc10x),score="Stress",geneset=as.list(genes.stress),cut.pt="renyi",anchor=anchor)
-sc10x.preStress <- results[[1]]
-sc10x <- results[[2]]
+results <- scScore(list(ALL=sc10x),score="Stress",geneset=as.list(genes.stress),cut.pt="triangle",anchor=anchor)
+#sc10x.preStress <- results[[1]]$ALL
+sc10x <- results[[2]]$ALL
 #sc10x <- results[[1]]$ALL
 rm(results)
 #counts.cell.destress <- lapply(sc10x,ncol)
 counts.cell.destress <- ncol(sc10x)
 rm(anchor)
 
+#sc10x <- SCTransform(sc10x,vars.to.regress=c("nFeature_RNA","percent.mito","Stress1"),verbose=FALSE,return.only.var.genes=FALSE,assay="RNA")
+
+save.image(file="./analysis/sc10x.raw_destress.RData")
+
+if (ncol(sc10x) > 1000) {
+  pc.calc <- 1000
+} else if (ncol(sc10x) > 500) {
+  pc.calc <- 500
+} else if (ncol(sc10x) > 100) {
+  pc.calc <- ncol(sc10x)-1
+}
+results <- scPC(sc10x,pc=pc.calc,hpc=0.9,file="ALL",print="umap",assay="integrated")
+sc10x <- results[[1]]
+pc.use.poststress <- results[[2]]
+rm(results)
+rm(pc.calc)
+
+res <- c(seq(0.1,0.5,0.1),0.75,seq(1,5,1))
+
+sc10x <- scCluster(sc10x,res=res,red="pca",dim=pc.use.poststress,print="umap",folder="ALL",assay="integrated")
+
 DefaultAssay(object=sc10x) <- "SCT"
 
 sc10x@meta.data <- sc10x@meta.data[,c("samples","integrated_snn_res.0.1","integrated_snn_res.0.2","integrated_snn_res.0.3","integrated_snn_res.0.4","integrated_snn_res.0.5","integrated_snn_res.0.75","integrated_snn_res.1","integrated_snn_res.2","integrated_snn_res.3","integrated_snn_res.4","integrated_snn_res.5","nCount_RNA","nFeature_RNA","percent.mito","percent.ribo","Stress1")]

r.scripts/sc-TissueMapper_functions.R

@@ -547,7 +547,7 @@ scAlign <-  function(sc10x.l){
     gc()
     #sc10x.l[[i]] <- ScaleData(sc10x.l[[i]],vars.to.regress=c("nFeature_RNA","percent.mito"),verbose = FALSE)
     #sc10x.l[[i]] <- SCTransform(sc10x.l[[i]],variable.features.n=5000,return.only.var.genes=FALSE,vars.to.regress=c("nFeature_RNA","percent.mito","Stress1"),verbose=FALSE,assay="RNA")
-    sc10x.l[[i]] <- SCTransform(sc10x.l[[i]],variable.features.n=10000,return.only.var.genes=FALSE,vars.to.regress=c("nFeature_RNA","percent.mito"),verbose=FALSE,assay="RNA")
+    sc10x.l[[i]] <- SCTransform(sc10x.l[[i]],variable.features.n=5000,return.only.var.genes=FALSE,vars.to.regress=c("nFeature_RNA","percent.mito"),verbose=FALSE,assay="RNA")
     gc()
     #sc10x.l[[i]] <- FindVariableFeatures(sc10x.l[[i]],selection.method="vst",nfeatures=2000,verbose=FALSE)
   }
@@ -573,7 +573,7 @@ scAlign <-  function(sc10x.l){
   #sc10x <- IntegrateData(anchorset=sc10x,dims=1:50)
   
   gc()
-  sc10x <- ScaleData(object=sc10x,do.par=TRUE,num.cores=45,verbose=FALSE,assay="integrated")
+  sc10x <- ScaleData(object=sc10x,verbose=FALSE,assay="integrated")
   gc()
   
   gc()