前言
近期匿醒，張澤民教授團(tuán)隊(duì)開發(fā)了一種在cell cluster里面計(jì)算純度的算法ROGUE邪意，什么意思呢骨望？經(jīng)過傳統(tǒng)的降維方式（t-SNE粪牲，UMAP）降維聚類后的cell cluster，其每一個(gè)cell cluster里面并不會(huì)是100%的同一類型的細(xì)胞赢织，里面可能含有一些其他類型的細(xì)胞也分到同一個(gè)cell cluster里面亮靴，因此有必要計(jì)算每一個(gè)cell cluster的純度
文章鏈接：An entropy-based metric for assessing the purity of single cell populations
測(cè)試算法
首先，我們模擬創(chuàng)建一個(gè)單細(xì)胞數(shù)據(jù)中某cell cluster的表達(dá)矩陣
exper = data.frame(abs(matrix(rnorm(10000),nrow=100)))
row.names(exper) = paste(rep('gene'),1:100,sep = '_')
colnames(exper) = paste(rep('cell'),1:100,sep = '_')
然后計(jì)算該矩陣的平均表達(dá)量和香濃熵敌厘，文章是這么定義每個(gè)基因的香濃熵的：
那么從代碼中台猴，我們可以看到，其香濃熵即為每行的均值俱两，我們稱之為真實(shí)entropy饱狂，實(shí)現(xiàn)由函數(shù)：Entropy：
library(tibble)
Entropy <- function(expr, r = 1){
tmp <- log(expr+1)
entropy <- Matrix::rowMeans(tmp) #計(jì)算香濃熵
mean.expr <- log(Matrix::rowMeans(expr)+r) #計(jì)算表達(dá)均值的log值
ent_res <- tibble(
Gene = rownames(expr),
mean.expr = mean.expr,
entropy = entropy
)
return(ent_res)
}
結(jié)果為：
計(jì)算好每個(gè)基因表達(dá)量均值的log值和每個(gè)基因的香濃熵以后，作者利用loess來擬合基因表達(dá)量均值的log值和每個(gè)基因的香濃熵之間的關(guān)系宪彩，實(shí)現(xiàn)由函數(shù)：entropy_fit：
entropy_fit <- function(.x, span = 0.5, mt.method = "fdr"){
.x <- .x %>% dplyr::filter(is.finite(mean.expr)) %>% dplyr::filter(entropy > 0)
fit <- loess(entropy~mean.expr, data = .x, span=span)
prd <- predict(fit, .x$mean.expr)
.x %>%
dplyr::mutate(fit = prd) %>%
dplyr::mutate(ds = fit - entropy) %>%
dplyr::mutate(pv = 1-pnorm(.$ds, mean = mean(.$ds), sd = sd(.$ds))) %>%
dplyr::filter(pv > 0.1) -> tmp
fit <- loess(entropy~mean.expr, data = tmp, span=span)
prd <- predict(fit, .x$mean.expr)
.x %>%
dplyr::mutate(fit = prd) %>%
dplyr::mutate(ds = fit - entropy) %>%
dplyr::filter(is.finite(ds)) %>%
dplyr::mutate(pv = 1-pnorm(.$ds, mean = mean(.$ds), sd = sd(.$ds))) %>%
dplyr::filter(pv > 0.1) -> tmp
fit <- loess(entropy~mean.expr, data = tmp, span=span)
prd <- predict(fit, .x$mean.expr)
.x %>%
dplyr::mutate(fit = prd) %>%
dplyr::mutate(ds = fit - entropy) %>%
dplyr::filter(is.finite(ds)) -> .x
.x <- .x %>% dplyr::mutate(p.value = 1-pnorm(.x$ds, mean = mean(.x$ds), sd = sd(.x$ds)))
p.adj <- p.adjust(.x$p.value, method = mt.method)
.x <- .x %>% dplyr::mutate(p.adj = p.adj) %>% dplyr::arrange(desc(ds))
}
結(jié)果為：
這里用loess（局部加權(quán)回歸）來預(yù)測(cè)基因平均表達(dá)量的log值與熵的關(guān)系休讳，預(yù)測(cè)的模型為fit，而后者需要利用基因表達(dá)量均值的log值作為決策變量尿孔，來根據(jù)模型fit預(yù)測(cè)響應(yīng)變量的值俊柔，我們稱之為預(yù)測(cè)entropy
并且作者定義ds如下：
由代碼可知
dplyr::mutate(fit = prd) %>%
dplyr::mutate(ds = fit - entropy)
ds為預(yù)測(cè)entropy減去真實(shí)entropy，ds越大說明該cluster內(nèi)某基因的熵變很大活合，即真實(shí)表達(dá)水平與預(yù)測(cè)的（與cluster內(nèi)部大部分基因不一樣）相差很大雏婶，說明很有可能是其他類型的細(xì)胞（雜質(zhì)）
最后要計(jì)算的是Rogue值
該值定義如下：
其中而ROGUE介于0-1之間；當(dāng)平臺(tái)是UMI的白指，那么K=45留晚，若平臺(tái)為full-length，則K=500
SE_fun <- function(expr, span = 0.5, r = 1, mt.method = "fdr", if.adj = T){
ent_res <- ROGUE::Entropy(expr, r = r)
ent_res <- ROGUE::entropy_fit(ent_res, span = span, mt.method = mt.method)
if(!isTRUE(if.adj)){
ent_res <- ent_res %>% dplyr::mutate(p.adj = p.value)
}
return(ent_res)
}
## .x即為函數(shù)SE_fun()的返回值
CalculateRogue <- function(.x, platform = NULL, cutoff = 0.05, k = NULL, features = NULL){
if(is.null(k)){
if(is.null(platform)){
warning("Please provide a \"platform\" argument or specify a k value")
}else if(platform == "UMI"){
k = 45
}else if(platform == "full-length"){
k = 500
}else if(!is.null(platform) & !(platform %in% c("UMI","full-length"))){
warning("Please provide valid \"platform\" argument")
}
}else if(!is.null(k)){
k <- k
}
if(!is.null(features)){
.x <- .x %>% dplyr::filter(Gene %in% features)
sig_value <- sum(abs(.x$ds))
Rogue <- 1-sig_value/(sig_value+k)
return(Rogue)
}else{
sig_value <- abs(.x$ds[.x$p.adj < cutoff & .x$p.value < cutoff])
sig_value <- sum(sig_value)
Rogue <- 1-sig_value/(sig_value+k)
return(Rogue)
}
}
其中.x即為函數(shù)SE_fun()的返回值告嘲，在計(jì)算sig_value時(shí)错维，用的是該cell cluster中所有基因的ds值的總和；ROGUE值越高橄唬，說明該cell cluster細(xì)胞純度越高赋焕，反之越低；
Github：https://github.com/PaulingLiu/ROGUE/blob/master/R/ROGUE.R
