利用貝葉斯的方法獲得cell cluster的marker基因

理論

參考文章為：genesorteR
簡單理解下，每個cell type的marker基因，它們的表達量一定具有cell type特異性的

假設單細胞表達矩陣為m×n的單細胞表達矩陣，m個基因和n個cell，并且n個細胞劃分到了k個cell cluster里面，作者通過貝葉斯公式：



來反應每個cell cluster中的基因特異性

其中：

  1. t ∈ { t1 ， t2 ，...， tk }，代表不同的cell cluster
  2. P( ti | gj ) 代表在檢測到 gene j （gj）有表達的條件下，觀測該cell（單個cell）屬于 cell cluster ti 的概率；其中 gj 代表 gene j
  3. P( gj | ti ) 代表在 cell cluster ti 的細胞中檢測到基因 gj 有表達的概率
  4. P( gj ) 代表全部的 n 個細胞中都檢測到 gj 的概率
  5. P( ti ) 代表全部的 n 個細胞屬于 cell cluster ti 的概率

那么試想哪一種基因?qū)儆趍arker呢？也就是說 P( ti | gj ) 值越大，則 gj 為 cell cluster ti 的marker 基因概率越大，根據(jù)定義，當 gj 有表達的時候，劃分某一個細胞為 cell cluster ti 的概率越大，反而越說明 gj 為 cell cluster ti 的marker基因

接下來作者定義：


為cell cluster的特異性分數(shù)，當 si j 的值為 1 時，代表在 cell cluster ti 中的細胞都檢測到了 gj，而其他類型的細胞中沒有檢測到 gj

定義marker基因特異性分數(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 為表達的單細胞稀疏矩陣
  x = kidneyTabulaMuris$exp
  # kidneyTabulaMuris$cellType 為每個細胞對應的細胞類型
  classLabels = kidneyTabulaMuris$cellType
  binarizeMethod = "median"
  TF_IDF = FALSE
  returnInput = TRUE
  cores = 1
  
  # 因子化標記每個細胞對應的細胞類型
  classLabels = as.factor(classLabels)
  ww = which(as.vector(table(classLabels)) == 1)    

  # 將不同的因子類型轉(zhuǎn)換為數(shù)字
  classLabelsNum = as.integer(classLabels)
  # 將上一步的結(jié)果因子化
  classLabelsLab = levels(classLabels)
  
  # 轉(zhuǎn)換為 dgCMatrix
  x = as(x, "dgCMatrix")
  
  # 計算稀疏矩陣有表達(有數(shù)值的,沒有數(shù)值的不納入計算)基因表達量的中位數(shù)
  ## 將所有細胞中各個有表達(有數(shù)值的,沒有數(shù)值的不納入計算)基因按每一列(每個細胞)為單位拼成一列向量，并計算中位數(shù)
  xbin = binarize(x, method = binarizeMethod)
  
  # 稀疏矩陣求行總和,并去除0表達的情況
  rem = which((Matrix::rowSums(xbin$mat)) == 0)
  if (length(rem) > 0) {
    xbin$mat = xbin$mat[-rem,]
  }
  
  # 求每一個細胞類別(1,2,3)的概率,每一類別的數(shù)量除以總數(shù)量
  classProb = getClassProb(classLabelsNum)
  # 命名
  names(classProb) = classLabelsLab
  
  # 計算每個基因的概率,每個基因在各個細胞表達量的總和除以總的細胞數(shù)
  geneProb = getGeneProb(xbin$mat)
  
  # 計算每一類細胞(cell cluster)中每個基因的平均表達量，并定義為概率 P( gj | ti )
  condGeneCluster = getGeneConditionalCluster_stats(xbin$mat, classProb, classLabels, cores = cores)
  # 利用貝葉斯公式計算每個基因的后驗概率
  clusterPostGene = getClusterPostGene(condGeneCluster, geneProb, classProb)
  # 計算特異性分數(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幾個函數(shù)：

# 1 計算稀疏矩陣有表達(有數(shù)值的,沒有數(shù)值的不納入計算)基因表達量的中位數(shù)
binarize = function(x, method = "median") {
  # 計算稀疏矩陣有表達(有數(shù)值的,沒有數(shù)值的不納入計算)基因表達量的中位數(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 求每一個細胞類別(1,2,3)的概率,每一類別的數(shù)量除以總數(shù)量
getClassProb = function(x) {
  probs = table(x) / length(x)
  
  return(probs)
}

# 3 計算每個基因的概率,每個基因在各個細胞表達量的總和除以總的細胞數(shù)
getGeneProb = function(x) {
  ncell = ncol(x)
  
  probs = as.vector( (Matrix::rowSums(x)) / ncell )
  names(probs) = rownames(x)
  
  return(probs)
}

# 4 計算每一類細胞(cell cluster)中每個基因的平均表達量，并定義為概率 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 利用貝葉斯公式計算每個基因的后驗概率
getClusterPostGene = function(condMat, geneProb, classProb) {
  
  # 計算貝葉斯公式的分子部分
  postProb = t(apply(log(condMat), 1, function(x) x + log(classProb)))
  # 計算貝葉斯公式的分母部分
  postProb = exp(apply(postProb, 2, function(x) x - log(geneProb)))
  
  colnames(postProb) = names(classProb)
  postProb = as(postProb, "dgCMatrix")
  
  return(postProb)  
}

# 6 計算特異性分數(shù)
getSpecScore = function(postMat, condMat) {
  # 計算特異性分數(shù)
  specScore = exp(log(postMat) + log(condMat))
  
  colnames(specScore) = colnames(condMat)
  specScore = as(specScore, "dgCMatrix")
  return(specScore)
}

注意這里的特異性分數(shù)的計算，由于變量condGeneCluster代表 P( gj | ti )，因此specScore = exp(log(postMat) + log(condMat))代表分數(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
  
  # 對 sij 分數(shù)進行標準化
  scored = apply(gs$specScore, 2, score)
  # 計算香濃熵
  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))
}

幾個函數(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)
}
?著作權(quán)歸作者所有,轉(zhuǎn)載或內(nèi)容合作請聯(lián)系作者
【社區(qū)內(nèi)容提示】社區(qū)部分內(nèi)容疑似由AI輔助生成，瀏覽時請結(jié)合常識與多方信息審慎甄別。
平臺聲明：文章內(nèi)容（如有圖片或視頻亦包括在內(nèi)）由作者上傳并發(fā)布，文章內(nèi)容僅代表作者本人觀點，簡書系信息發(fā)布平臺，僅提供信息存儲服務。

相關(guān)閱讀更多精彩內(nèi)容

友情鏈接更多精彩內(nèi)容