Rejigger thresholding

r.scripts/sc-TissueMapper_RUN.R

@@ -36,7 +36,7 @@ project.name=opt$p
 
 scFolders()
 
-results <- scLoad(p=project.name,cellranger=3,aggr=FALSE)
+results <- scLoad(p=project.name,cellranger=3,aggr=FALSE,ncell=0,nfeat=0)
 sc10x <- results[[1]]
 sc10x.groups <- results[[2]]
 rm(results)

r.scripts/sc-TissueMapper_functions.R

@@ -68,7 +68,7 @@ scFolders <- function(){
 }
 
 
-scLoad <- function(p,cellranger=3,aggr=TRUE){
+scLoad <- function(p,cellranger=3,aggr=TRUE,ncell=0,nfeat=0){
   #Load and prefilter filtered_gene_bc_matrices_mex output from cellranger
   
   #Inputs:
@@ -101,7 +101,7 @@ scLoad <- function(p,cellranger=3,aggr=TRUE){
       } else {
         sc10x.data[i] <- Read10X(data.dir=paste0("./analysis/DATA/10x/",i,"/filtered_feature_bc_matrix/"))
       }
-      sc10x[i] <- CreateSeuratObject(counts=sc10x.data[[i]],project=p)
+      sc10x[i] <- CreateSeuratObject(counts=sc10x.data[[i]],project=p,min.cells=ncell,min.features=nfeat)
       sc10x[[i]]$samples <- i
     }
   }
@@ -196,22 +196,38 @@ scQC <- function(sc10x,sp="hu"){
   
   #Calculate cutoffs
   thresh <- list()
-  h <- list()
-  cells.remove <- list()
-  sc10x.temp <- list()
-  thresh.l <- list()
-  cutoff.l.mito <- list()
+  #h <- list()
+  #cells.remove <- list()
+  #sc10x.temp <- list()
   for (i in c("nFeature_RNA","percent.mito")){
-    thresh[[i]] <- scThresh(sc10x,feature=i)
+    if (i == "nFeature_RNA"){
+    h <- list()
+    cells.remove <- list()
+    sc10x.temp <- list()
+    for (j in names(sc10x)){
+    h[[i]] <- hist(data.frame(sc10x[[j]][[i]])$nFeature_RNA,breaks=10,plot=FALSE)
+    cutoff.temp <- mean(c(h[[i]]$mids[which.max(h[[i]]$counts)],h[[i]]$mids[-which.max(h[[i]]$counts)][which.max(h[[i]]$counts[-which.max(h[[i]]$counts)])]))
+    cells.remove[[j]] <- c(cells.remove[[j]],rownames(sc10x[[j]][["nFeature_RNA"]])[sc10x[[j]][[i]][,1] < cutoff.temp])
+    sc10x.temp[[j]] <- subset(sc10x[[j]],cells=setdiff(colnames(sc10x[[j]]),cells.remove[[j]]))
+    }
+    thresh[[i]] <- scThresh(sc10x.temp,feature=i)
+    }
     if (i == "percent.mito"){
+      h <- list()
+      cells.remove <- list()
+      sc10x.temp <- list()
+      thresh.l <- list()
+      cutoff.l.mito <- list()
       for (j in names(sc10x)){
         cutoff.l.mito[[j]] <- NULL
-        h[[i]] <- hist(data.frame(sc10x[[j]][[i]])$percent.mito,breaks=1000,plot=FALSE)
+        h[[i]] <- hist(data.frame(sc10x[[j]][[i]])$percent.mito,breaks=100,plot=FALSE)
         cutoff.temp <- mean(c(h[[i]]$mids[which.max(h[[i]]$counts)],h[[i]]$mids[-which.max(h[[i]]$counts)][which.max(h[[i]]$counts[-which.max(h[[i]]$counts)])]))
         cells.remove[[j]] <- c(cells.remove[[j]],rownames(sc10x[[j]][["percent.mito"]])[sc10x[[j]][[i]][,1] > cutoff.temp])
         sc10x.temp[[j]] <- subset(sc10x[[j]],cells=setdiff(colnames(sc10x[[j]]),cells.remove[[j]]))
         thresh.l[[i]] <- scThresh(sc10x.temp,feature=i,sub="lower")
-        cutoff.l.mito[[j]] <- thresh.l[[i]][[j]]$threshold[thresh.l[[i]][[j]]$method=="Triangle"]
+        #cutoff.l.mito[[j]] <- thresh.l[[i]][[j]]$threshold[thresh.l[[i]][[j]]$method=="Triangle"]
+        cutoff.l.mito[[j]] <- thresh.l[[i]][[j]]$threshold[thresh.l[[i]][[j]]$method=="RenyiEntropy"]
+        thresh[[i]] <- scThresh(sc10x,feature=i)
       }
     }
   }
@@ -253,11 +269,14 @@ scQC <- function(sc10x,sp="hu"){
       plots.v[[j]] <- VlnPlot(object=sc10x.temp,features=i,pt.size=0.1,)+scale_x_discrete(labels=j)+scale_y_continuous(limits=c(0,max.y))+theme(legend.position="none",axis.text.x=element_text(hjust=0.5,angle=0))
       if (i %in% c("nFeature_RNA","percent.mito")){
         if (i == "nFeature_RNA"){
-          cutoff.l <- thresh[[i]][[j]]$threshold[thresh[[i]][[j]]$method=="MinErrorI"]
+          #cutoff.l <- thresh[[i]][[j]]$threshold[thresh[[i]][[j]]$method=="MinErrorI"]
           cutoff.h <- thresh[[i]][[j]]$threshold[thresh[[i]][[j]]$method=="RenyiEntropy"]
+          cutoff.l <- 200
+          #cutoff.h <- thresh[[i]][[j]]$threshold[thresh[[i]][[j]]$method=="Huang2"]
         } else {
-          cutoff.l <- cutoff.l.mito[[j]]
+          #cutoff.l <- cutoff.l.mito[[j]]
           cutoff.h <- thresh[[i]][[j]]$threshold[thresh[[i]][[j]]$method=="Triangle"]
+          cutoff.l <- 0
         }
         plots.v[[j]] <- plots.v[[j]]+geom_hline(yintercept=cutoff.l,size=0.5,color="red")+geom_hline(yintercept=cutoff.h,size=0.5,color="red")
         densities.s[[j]] <- density(sc10x.temp$nCount_RNA,sc10x.temp[[i]][,1],n=1000)
@@ -378,8 +397,7 @@ scThresh <- function(sc10x,feature,sub=FALSE){
     scale.scaled[[i]] <- as.integer((scale[[i]]$Score-min(scale[[i]]$Score))/(max(scale[[i]]$Score)-min(scale[[i]]$Score))*360)
     #scale.scaled[[i]] <- as.integer(scales::rescale(scale[[i]]$Score,to=c(0,1))*360)
     h[[i]] <- hist(scale[[i]]$Score,breaks=100,plot=FALSE)
-    thresh[[i]] <- purrr::map_chr(thresh_methods, ~auto_thresh(scale.scaled[[i]], .)) %>% 
-      tibble(method = thresh_methods, threshold = .)
+    thresh[[i]] <- purrr::map_chr(thresh_methods,~auto_thresh(scale.scaled[[i]],.)) %>% tibble(method = thresh_methods, threshold = .)
     thresh[[i]]$threshold <- as.numeric(thresh[[i]]$threshold)
     thresh[[i]]$threshold <- ((thresh[[i]]$threshold/360)*(max(scale[[i]]$Score)-min(scale[[i]]$Score)))+min(scale[[i]]$Score)
     #thresh[[i]] <- thresh[[i]] %>% mutate(threshold=(scales::rescale(as.numeric(threshold)/360,to=range(scale[[i]]$Score))))