理論
參考文章為：genesorteR
簡(jiǎn)單理解下茸苇，每個(gè)cell type的marker基因，它們的表達(dá)量一定具有cell type特異性的
假設(shè)單細(xì)胞表達(dá)矩陣為m×n的單細(xì)胞表達(dá)矩陣央星，m個(gè)基因和n個(gè)cell飘诗，并且n個(gè)細(xì)胞劃分到了k個(gè)cell cluster里面，作者通過(guò)貝葉斯公式：
來(lái)反應(yīng)每個(gè)cell cluster中的基因特異性
其中：
- t ∈ { t1 ， t2 致板，...， tk }咏窿，代表不同的cell cluster
- P( ti | gj ) 代表在檢測(cè)到 gene j （gj）有表達(dá)的條件下斟或，觀測(cè)該cell（單個(gè)cell）屬于 cell cluster ti 的概率；其中 gj 代表 gene j
- P( gj | ti ) 代表在 cell cluster ti 的細(xì)胞中檢測(cè)到基因 gj 有表達(dá)的概率
- P( gj ) 代表全部的 n 個(gè)細(xì)胞中都檢測(cè)到 gj 的概率
- P( ti ) 代表全部的 n 個(gè)細(xì)胞屬于 cell cluster ti 的概率
那么試想哪一種基因?qū)儆趍arker呢集嵌？也就是說(shuō) P( ti | gj ) 值越大萝挤，則 gj 為 cell cluster ti 的marker 基因概率越大，根據(jù)定義根欧，當(dāng) gj 有表達(dá)的時(shí)候怜珍，劃分某一個(gè)細(xì)胞為 cell cluster ti 的概率越大，反而越說(shuō)明 gj 為 cell cluster ti 的marker基因
接下來(lái)作者定義：
為cell cluster的特異性分?jǐn)?shù)凤粗，當(dāng) si j 的值為 1 時(shí)酥泛，代表在 cell cluster ti 中的細(xì)胞都檢測(cè)到了 gj，而其他類(lèi)型的細(xì)胞中沒(méi)有檢測(cè)到 gj
定義marker基因特異性分?jǐn)?shù)：
代碼部分
完整的流程為：
library(genesorteR)
data(kidneyTabulaMuris) #three cell types from kidney (Tabula Muris data)
#get specificity scores for each cell type
sg = sortGenes(kidneyTabulaMuris$exp, kidneyTabulaMuris$cellType)
head(sg$specScore) #specificity scores for each gene in each cluster
#define a small set of markers
mm = getMarkers(sg, quant = 0.99)
#cluster genes and make a heatmap
pp = plotMarkerHeat(sg$inputMat, sg$inputClass, mm$markers, clusterGenes=TRUE,
outs = TRUE)
pp$gene_class_info #gene clusters
1.函數(shù) sortGenes
sortGenes = function(x, classLabels, binarizeMethod = "median", TF_IDF = FALSE, returnInput = TRUE, cores = 1) {
# kidneyTabulaMuris$exp 為表達(dá)的單細(xì)胞稀疏矩陣
x = kidneyTabulaMuris$exp
# kidneyTabulaMuris$cellType 為每個(gè)細(xì)胞對(duì)應(yīng)的細(xì)胞類(lèi)型
classLabels = kidneyTabulaMuris$cellType
binarizeMethod = "median"
TF_IDF = FALSE
returnInput = TRUE
cores = 1
# 因子化標(biāo)記每個(gè)細(xì)胞對(duì)應(yīng)的細(xì)胞類(lèi)型
classLabels = as.factor(classLabels)
ww = which(as.vector(table(classLabels)) == 1)
# 將不同的因子類(lèi)型轉(zhuǎn)換為數(shù)字
classLabelsNum = as.integer(classLabels)
# 將上一步的結(jié)果因子化
classLabelsLab = levels(classLabels)
# 轉(zhuǎn)換為 dgCMatrix
x = as(x, "dgCMatrix")
# 計(jì)算稀疏矩陣有表達(dá)(有數(shù)值的,沒(méi)有數(shù)值的不納入計(jì)算)基因表達(dá)量的中位數(shù)
## 將所有細(xì)胞中各個(gè)有表達(dá)(有數(shù)值的,沒(méi)有數(shù)值的不納入計(jì)算)基因按每一列(每個(gè)細(xì)胞)為單位拼成一列向量嫌拣，并計(jì)算中位數(shù)
xbin = binarize(x, method = binarizeMethod)
# 稀疏矩陣求行總和,并去除0表達(dá)的情況
rem = which((Matrix::rowSums(xbin$mat)) == 0)
if (length(rem) > 0) {
xbin$mat = xbin$mat[-rem,]
}
# 求每一個(gè)細(xì)胞類(lèi)別(1,2,3)的概率,每一類(lèi)別的數(shù)量除以總數(shù)量
classProb = getClassProb(classLabelsNum)
# 命名
names(classProb) = classLabelsLab
# 計(jì)算每個(gè)基因的概率,每個(gè)基因在各個(gè)細(xì)胞表達(dá)量的總和除以總的細(xì)胞數(shù)
geneProb = getGeneProb(xbin$mat)
# 計(jì)算每一類(lèi)細(xì)胞(cell cluster)中每個(gè)基因的平均表達(dá)量柔袁，并定義為概率 P( gj | ti )
condGeneCluster = getGeneConditionalCluster_stats(xbin$mat, classProb, classLabels, cores = cores)
# 利用貝葉斯公式計(jì)算每個(gè)基因的后驗(yàn)概率
clusterPostGene = getClusterPostGene(condGeneCluster, geneProb, classProb)
# 計(jì)算特異性分?jǐn)?shù)
specScore = getSpecScore(clusterPostGene, condGeneCluster)
tf_idf = NULL
if (TF_IDF) {
tf_idf = getClusterTFIDF(condGeneCluster)
}
result = list(binary = xbin$mat, cutoff = xbin$cutoff,
removed = rem, geneProb = geneProb,
condGeneProb = condGeneCluster, postClustProb = clusterPostGene,
specScore = specScore, classProb = classProb,
inputMat = x, inputClass = classLabels, tf_idf = tf_idf)
}
小tip幾個(gè)函數(shù)：
# 1 計(jì)算稀疏矩陣有表達(dá)(有數(shù)值的,沒(méi)有數(shù)值的不納入計(jì)算)基因表達(dá)量的中位數(shù)
binarize = function(x, method = "median") {
# 計(jì)算稀疏矩陣有表達(dá)(有數(shù)值的,沒(méi)有數(shù)值的不納入計(jì)算)基因表達(dá)量的中位數(shù)
if (method == "median") {
pi = median(x@x)
} else if (method == "mean") {
pi = mean(x@x)
} else if (method == "naive") {
pi = min(x@x)
} else if (method == "adaptiveMedian") {
mm = Mclust(Matrix::rowSums(x), 1:20, modelNames=c("V"), verbose = FALSE)
if (mm$G == 1) {
stop("Error: you set binarizeMethod to adaptiveMedian but the optimal number of gene groups based on average expression is 1. Please use a different binarization method or use a different gene expression normalization. Also please consider reporting this error to mmibrahim@pm.me or on https://github.com/mahmoudibrahim/genesorteR/issues (preferred).")
} else {
pi = rep(0, mm$G)
for (i in 1:mm$G) {
pi[i] = median(x[which(mm$classification == i),]@x)
}
}
} else if ((is.numeric(method)) & (method >= 0)) {
pi = method
} else {
stop("Unrecognized binarization method! genesorteR stopped.")
}
if (method == "adaptiveMedian") {
mat = list()
for (i in 1:mm$G) {
tx = x[which(mm$classification == i),]
ww = which(tx@x >= pi[i])
dp = diff(tx@p)
colInd = (rep(seq_along(dp),dp))[ww]
rowInd = (tx@i+1)[ww]
genes = rownames(tx)
mat[[i]] = Matrix::sparseMatrix(rowInd[1], colInd[1], dims=tx@Dim)
mat[[i]][cbind(rowInd,colInd)] = 1
mat[[i]] = as(mat[[i]], "dgCMatrix")
rownames(mat[[i]]) = genes
}
mat = do.call(rbind, mat)
x = mat[match(rownames(x),rownames(mat)),]
} else {
ww = which(x@x >= pi)
dp = diff(x@p)
colInd = (rep(seq_along(dp),dp))[ww]
rowInd = (x@i+1)[ww]
genes = rownames(x)
x = Matrix::sparseMatrix(rowInd[1], colInd[1], dims=x@Dim)
x[cbind(rowInd,colInd)] = 1
x = as(x, "dgCMatrix")
rownames(x) = genes
}
return(list(mat = x, cutoff = pi))
}
# 2 求每一個(gè)細(xì)胞類(lèi)別(1,2,3)的概率,每一類(lèi)別的數(shù)量除以總數(shù)量
getClassProb = function(x) {
probs = table(x) / length(x)
return(probs)
}
# 3 計(jì)算每個(gè)基因的概率,每個(gè)基因在各個(gè)細(xì)胞表達(dá)量的總和除以總的細(xì)胞數(shù)
getGeneProb = function(x) {
ncell = ncol(x)
probs = as.vector( (Matrix::rowSums(x)) / ncell )
names(probs) = rownames(x)
return(probs)
}
# 4 計(jì)算每一類(lèi)細(xì)胞(cell cluster)中每個(gè)基因的平均表達(dá)量，并定義為概率 P( gj | ti )
getGeneConditionalCluster_stats = function(mat, classProb, classLabels, cores = 1) {
lab = names(classProb)
if (cores > 1) {
geneProb = do.call(cbind, mclapply(1:length(lab), function(i) Matrix::rowMeans(mat[,which(classLabels == lab[i])]), mc.cores = cores))
} else {
geneProb = do.call(cbind, lapply(1:length(lab), function(i) Matrix::rowMeans(mat[,which(classLabels == lab[i])])))
}
geneProb = as(geneProb, "dgCMatrix")
colnames(geneProb) = lab
return(geneProb)
}
# 5 利用貝葉斯公式計(jì)算每個(gè)基因的后驗(yàn)概率
getClusterPostGene = function(condMat, geneProb, classProb) {
# 計(jì)算貝葉斯公式的分子部分
postProb = t(apply(log(condMat), 1, function(x) x + log(classProb)))
# 計(jì)算貝葉斯公式的分母部分
postProb = exp(apply(postProb, 2, function(x) x - log(geneProb)))
colnames(postProb) = names(classProb)
postProb = as(postProb, "dgCMatrix")
return(postProb)
}
# 6 計(jì)算特異性分?jǐn)?shù)
getSpecScore = function(postMat, condMat) {
# 計(jì)算特異性分?jǐn)?shù)
specScore = exp(log(postMat) + log(condMat))
colnames(specScore) = colnames(condMat)
specScore = as(specScore, "dgCMatrix")
return(specScore)
}
注意這里的特異性分?jǐn)?shù)的計(jì)算异逐，由于變量condGeneCluster代表 P( gj | ti )捶索，因此
specScore = exp(log(postMat) + log(condMat))代表分?jǐn)?shù)si j
2. 基于香農(nóng)熵選出marker基因
library(mclust)
getMarkers = function(gs, quant = 0.99, mutualInfo = FALSE, classEnt = FALSE) {
gs = sg
quant = 0.99
mutualInfo = FALSE
classEnt = FALSE
# 對(duì) sij 分?jǐn)?shù)進(jìn)行標(biāo)準(zhǔn)化
scored = apply(gs$specScore, 2, score)
# 計(jì)算香濃熵
ent = apply(scored, 1, getEntropy)
# 排序并選出最大的值
maxi = apply(scored, 1, max)
ww = which(maxi > quantile(maxi, probs=quant))
m = Mclust(ent[ww], 2, modelNames="E", verbose = FALSE)
# 選出marker基因
markers = names( which(m$classification == (which.min(m$parameters$mean)) ) )
mutInfo = NULL
if (mutualInfo == TRUE) {
mutInfo = apply(gs$binary, 1, function(x) getMutInfo(x,gs$inputClass))
}
classEntropy = NULL
if (classEnt == TRUE) {
classEntropy = apply(scored, 2, getEntropy)
names(classEntropy) = colnames(gs$specScore)
}
return(list(gene_shannon_index = ent, maxScaledSpecScore = maxi, markers = markers,
mutInfo = mutInfo, classEntropy = classEntropy))
}
幾個(gè)函數(shù)：
score = function(x) {
x = ((x-min(x))/(max(x)-min(x)))
return(x)
}
getEntropy = function(x) {
ent = 0
w = which(x > 0)
if (length(w) > 0) {
w = x[w]
ent = -(sum(w * (log(w))))
}
return(ent)
}
