Commit 24967dcc authored by Gervaise Henry's avatar Gervaise Henry 🤠
Browse files

Make lower cutoff for nCount with filtered cells

parent 5d1064b7
......@@ -221,6 +221,7 @@ scQC <- function(sc10x,sp="hu",feature="nFeature_RNA"){
thresh[[i]] <- scThresh(sc10x.temp,feature=i)
}
if (i == "percent.mito"){
thresh[[i]] <- scThresh(sc10x,feature=i)
h <- list()
cells.remove <- list()
sc10x.temp <- list()
......@@ -228,31 +229,37 @@ scQC <- function(sc10x,sp="hu",feature="nFeature_RNA"){
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=100,plot=FALSE)
temp <- data.frame(sc10x[[j]][[i]])$percent.mito
h[[i]] <- hist(temp,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=="RenyiEntropy"]
thresh[[i]] <- scThresh(sc10x,feature=i)
}
}
if (i == "percent.ribo"){
thresh[[i]] <- scThresh(sc10x,feature=i)
}
if (i == "nCount_RNA"){
thresh[[i]] <- scThresh(sc10x,feature=i)
h <- list()
cells.remove <- list()
sc10x.temp <- list()
thresh.l <- list()
cutoff.l.ribo <- list()
cutoff.l.count <- list()
for (j in names(sc10x)){
cutoff.l.ribo[[j]] <- NULL
h[[i]] <- hist(data.frame(sc10x[[j]][[i]])$percent.ribo,breaks=100,plot=FALSE)
thresh[[i]] <- scThresh(sc10x,feature=i)
cutoff.l.count[[j]] <- NULL
temp <- data.frame(sc10x[[j]][[i]])$nCount_RNA
h[[i]] <- hist(temp,breaks=seq(min(temp),max(temp),length.out=6),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]][["nCount_RNA"]])[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.count[[j]] <- thresh.l[[i]][[j]]$threshold[thresh.l[[i]][[j]]$method=="RenyiEntropy"]
}
}
if (i == "nCount_RNA"){
thresh[[i]] <- scThresh(sc10x,feature=i)
}
}
#Plot raw stats
......@@ -310,7 +317,8 @@ scQC <- function(sc10x,sp="hu",feature="nFeature_RNA"){
cutoff.h <- thresh[[i]][[j]]$threshold[thresh[[i]][[j]]$method=="RenyiEntropy"]
cutoff.l <- 0
} else if (i == "nCount_RNA") {
cutoff.l <- thresh[[i]][[j]]$threshold[thresh[[i]][[j]]$method=="MinErrorI"]
#cutoff.l <- thresh[[i]][[j]]$threshold[thresh[[i]][[j]]$method=="MinErrorI"]
cutoff.l <- cutoff.l.count[[j]]
cutoff.h <- max(sc10x[[j]][[i]])
}
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")
......
Markdown is supported
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