wguo-research / scCancer

A package for automated processing of single cell RNA-seq data in cancer
100 stars 42 forks source link

Malignant cell threshold? (malign.thres <- getBimodalThres(scores = obserScore.smooth)?) #35

Open khl0798 opened 3 years ago

khl0798 commented 3 years ago

why getBimodalThres is used to get the score distribution of malignant cells? getBimodalThres function is based on bimodal distribution. infact the maligantscore‘s actual distribution is not a bimodal distribution?

1635485463(1)

1635485538(1)

Does the red line indicate a bimodal distribution?

getBimodalThres <- function(scores){ x.density <- density(scores) d.x.density <- diff(x.density$y) d.sign <- (d.x.density > 0) + 0

ext.pos <- which(d.sign[2:length(d.sign)] - d.sign[1:(length(d.sign)-1)] != 0)
ext.density <- x.density$y[ext.pos]
y.max <- max(ext.density)
if(length(ext.pos) >= 3){
    del.ix <- c()
    for(ei in 2:length(ext.density)){
        if(abs(ext.density[ei] - ext.density[ei - 1]) < y.max * 0.001){
            del.ix <- c(del.ix, ei - 1, ei)
        }
    }
    sel.ix <- !(1:length(ext.density) %in% unique(del.ix))
    ext.density <- ext.density[sel.ix]
    ext.pos <- ext.pos[sel.ix]
}

if(length(ext.pos) >= 3){
    t.ext.density <- c(0, ext.density, 0)
    ext.height <- sapply(2:(length(ext.pos) + 1), FUN = function(x){
        return(min(abs(t.ext.density[x] - t.ext.density[x-1]), abs(t.ext.density[x] - t.ext.density[(x+1)])))
    })
    ext <- data.frame(x = ext.pos, y = ext.density, height = ext.height)
    max.ix <- order(ext.density, decreasing = T)
    if(ext.height[max.ix[2]] / ext.height[max.ix[1]] > 0.01){
        cut.df <- ext[c(max.ix[2]:max.ix[1]), ]
        threshold <- x.density$x[cut.df[which.min(cut.df$y), ]$x]
    }else{
        threshold <- NULL
    }
}else{
    threshold <- NULL
}

return(threshold)

}