Manage weird percent.mito=NA in QC/thresh

r.scripts/sc-TissueMapper_functions.R

@@ -110,7 +110,7 @@ scSubset <- function(sc10x,i="ALL",g="ALL"){
 }
 
 
-scQC <- function(sc10x,lu=0,sub=FALSE,sp="hu"){
+scQC <- function(sc10x,sub=FALSE,sp="hu"){
   #QC and filter Seurat object
   
   #Inputs:
@@ -142,6 +142,7 @@ scQC <- function(sc10x,lu=0,sub=FALSE,sp="hu"){
   } else {
     sc10x[["percent.mito"]] <- PercentageFeatureSet(object=sc10x,pattern="^mt-")
   }
+  sc10x <- subset(sc10x,cell=names(which(is.na(sc10x$percent.mito))),invert=TRUE)
   
   #Plot raw stats
   for(x in c("ALL","samples")){
@@ -160,9 +161,9 @@ scQC <- function(sc10x,lu=0,sub=FALSE,sp="hu"){
   cell.remove <- NULL
   Idents(object=sc10x) <- "samples"
   sc10x.sub <- sc10x
-  for (i in c("percent.mito","nFeature_RNA")){
+  for (i in c("nFeature_RNA","percent.mito")){
     if (i %in% c("nCount_RNA","nFeature_RNA")){
-        thresh <- scThresh(sc10x.sub,feature=i,direction="b",min=lu,name=i)
+        thresh <- scThresh(sc10x.sub,feature=i,direction="b",name=i)
       } else {
         thresh <- scThresh(sc10x.sub,feature=i,direction="h",name=i)
       }
@@ -175,7 +176,6 @@ scQC <- function(sc10x,lu=0,sub=FALSE,sp="hu"){
     postscript(paste0(folder,"Plot_qc.hist.",i,".eps"))
     plot(thresh[[3]])
     dev.off()
-    sc10x.sub <- subset(sc10x.sub,cells=WhichCells(sc10x,cells=cell.remove,invert=TRUE))
   }
   
   # threshPlot <- list()
@@ -250,20 +250,18 @@ scQC <- function(sc10x,lu=0,sub=FALSE,sp="hu"){
   return(results)
 }
 
-scThresh <- function(sc10x,feature,direction,min=0,name){
+scThresh <- function(sc10x,feature,direction,name){
   scale <- data.frame(Score=sc10x[[feature]])
   colnames(scale) <- "Score"
+  scale <- data.frame(Score=scale$Score[!is.na(scale$Score)])
   h <- hist(scale$Score,breaks=50)
   scale.s.h <- scale
-  scale.s.h <- data.frame(Score=scale.s.h[scale.s.h$Score>min,])
+  thresh.l <- mean(c(h$mids[which.max(h$counts)],h$mids[-which.max(h$counts)][which.max(h$counts[-which.max(h$counts)])]))
+  scale.s.h <- data.frame(Score=scale.s.h[scale.s.h$Score>thresh.l,])
   scale.max.h <- max(scale.s.h$Score)
   scale.s.h$Score <- as.integer(scale(scale.s.h$Score,FALSE,max(scale.s.h$Score))*100)
   thresh.h <- auto_thresh(scale.s.h$Score,method="Triangle")
   thresh.h <- (thresh.h/100)*scale.max.h
-  thresh.l <- mean(c(h$mids[which.max(h$counts)],h$mids[-which.max(h$counts)][which.max(h$counts[-which.max(h$counts)])]))
-  if (min>0){
-    thresh.l <- min
-  }
   threshPlot <- ggplot(scale,aes(x=Score))+geom_histogram(bins=100,aes(y=..density..))+geom_density()+ggtitle(name)
   if (direction=="h"){
     threshPlot <- threshPlot+geom_vline(xintercept=thresh.h,size=1,color="red")