Commit 180fbd64 authored by Gervaise Henry's avatar Gervaise Henry 🤠
Browse files

Create temp low end cutoff

parent 4682f616
......@@ -29,7 +29,7 @@ scFolders()
sc10x <- scLoad(p=project.name,cellranger=3,ref="GRCh38")
lu=1000
lu=0
results <- scQC(sc10x,lu=lu,sub=FALSE,sp="hu")
sc10x <- results[[1]]
counts.cell.raw <- results[[2]]
......
......@@ -159,27 +159,49 @@ scQC <- function(sc10x,lu=0,sub=FALSE,sp="hu"){
threshPlot <- list()
cell.remove <- NULL
Idents(object=sc10x) <- "samples"
sc10x.sub <- sc10x
for (i in c("percent.mito","nFeature_RNA")){
for (j in unique(sc10x$samples)){
if (i=="nFeature_RNA"){
thresh <- scThresh(scSubset(sc10x,i="samples",g=j),feature=i,direction="b",min=lu,name=j)
if (i %in% c("nCount_RNA","nFeature_RNA")){
thresh <- scThresh(sc10x.sub,feature=i,direction="b",min=lu,name=i)
} else {
thresh <- scThresh(scSubset(sc10x,i="samples",g=j),feature=i,direction="h",name=j)
thresh <- scThresh(sc10x.sub,feature=i,direction="h",name=i)
}
threshPlot[[j]] <- thresh[[3]]
data <- FetchData(subset(sc10x,idents=j),vars=i)
if (i=="nFeature_RNA"){
data <- FetchData(sc10x.sub,vars=i)
if (i %in% c("nCount_RNA","nFeature_RNA")){
cell.remove <- c(cell.remove,rownames(data)[data[,1] < thresh[[2]] | data[,1] > thresh[[1]]])
} else {
cell.remove <- c(cell.remove,rownames(data)[data[,1] > thresh[[1]]])
}
sc10x
}
postscript(paste0(folder,"Plot_qc.hist.",i,".eps"))
do.call("grid.arrange",c(threshPlot,ncol=1,top=i))
plot(thresh[[3]])
dev.off()
}
sc10x.sub <- subset(sc10x,cells=WhichCells(sc10x,cells=cell.remove,invert=TRUE))
sc10x.sub <- subset(sc10x.sub,cells=WhichCells(sc10x,cells=cell.remove,invert=TRUE))
}
# threshPlot <- list()
# cell.remove <- NULL
# Idents(object=sc10x) <- "samples"
# for (i in c("percent.mito","nFeature_RNA")){
# for (j in unique(sc10x$samples)){
# if (i=="nFeature_RNA"){
# thresh <- scThresh(scSubset(sc10x,i="samples",g=j),feature=i,direction="b",min=lu,name=j)
# } else {
# thresh <- scThresh(scSubset(sc10x,i="samples",g=j),feature=i,direction="h",name=j)
# }
# threshPlot[[j]] <- thresh[[3]]
# data <- FetchData(subset(sc10x,idents=j),vars=i)
# if (i=="nFeature_RNA"){
# cell.remove <- c(cell.remove,rownames(data)[data[,1] < thresh[[2]] | data[,1] > thresh[[1]]])
# } else {
# cell.remove <- c(cell.remove,rownames(data)[data[,1] > thresh[[1]]])
# }
# sc10x
# }
# postscript(paste0(folder,"Plot_qc.hist.",i,".eps"))
# do.call("grid.arrange",c(threshPlot,ncol=1,top=i))
# dev.off()
# }
# sc10x.sub <- subset(sc10x,cells=WhichCells(sc10x,cells=cell.remove,invert=TRUE))
postscript(paste0(folder,"Plot_qc.raw.eps"))
density1 <- density(sc10x$nCount_RNA,sc10x$nFeature_RNA,n=1000)
......@@ -231,18 +253,14 @@ scQC <- function(sc10x,lu=0,sub=FALSE,sp="hu"){
scThresh <- function(sc10x,feature,direction,min=0,name){
scale <- data.frame(Score=sc10x[[feature]])
colnames(scale) <- "Score"
h <- hist(scale$Score,breaks=100)
h <- hist(scale$Score,breaks=50)
scale.s.h <- scale
scale.s.h <- data.frame(Score=scale.s.h[scale.s.h$Score>min,])
scale.s.l <- data.frame(Score=scale[scale$Score<h$mids[which.max(h$counts)],])
scale.max.h <- max(scale.s.h$Score)
scale.max.l <- max(scale.s.l$Score)
scale.s.h$Score <- as.integer(scale(scale.s.h$Score,FALSE,max(scale.s.h$Score))*100)
scale.s.l$Score <- as.integer(scale(scale.s.l$Score,FALSE,max(scale.s.l$Score))*100)
thresh.h <- auto_thresh(scale.s.h$Score,method="Triangle")
thresh.l <- auto_thresh(scale.s.l$Score,method="Triangle")
thresh.h <- (thresh.h/100)*scale.max.h
thresh.l <- (thresh.l/100)*scale.max.l
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
}
......
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