理论

参考文章为：genesorteR
简单理解下，每个cell type的marker基因，它们的表达量一定具有cell type特异性的

假设单细胞表达矩阵为m×n的单细胞表达矩阵，m个基因和n个cell，并且n个细胞划分到了k个cell cluster里面，作者通过贝叶斯公式：

来反应每个cell cluster中的基因特异性

其中：

1. t ∈ { t₁ ， t₂ ，...， t_k }，代表不同的cell cluster
2. P( t_i | g_j ) 代表在检测到 gene j （g_j）有表达的条件下，观测该cell（单个cell）属于 cell cluster t_i 的概率；其中 g_j 代表 gene j
3. P( g_j | t_i ) 代表在 cell cluster t_i 的细胞中检测到基因 g_j 有表达的概率
4. P( g_j ) 代表全部的 n 个细胞中都检测到 g_j 的概率
5. P( t_i ) 代表全部的 n 个细胞属于 cell cluster t_i 的概率

那么试想哪一种基因属于marker呢？也就是说 P( t_i | g_j ) 值越大，则 g_j 为 cell cluster t_i 的marker 基因概率越大，根据定义，当 g_j 有表达的时候，划分某一个细胞为 cell cluster t_i 的概率越大，反而越说明 g_j 为 cell cluster t_i 的marker基因

接下来作者定义：

定义marker基因特异性分数：

代码部分

完整的流程为：

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.函数 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)    

  # 将不同的因子类型转换为数字
  classLabelsNum = as.integer(classLabels)
  # 将上一步的结果因子化
  classLabelsLab = levels(classLabels)
  
  # 转换为 dgCMatrix
  x = as(x, "dgCMatrix")
  
  # 计算稀疏矩阵有表达(有数值的,没有数值的不纳入计算)基因表达量的中位数
  ## 将所有细胞中各个有表达(有数值的,没有数值的不纳入计算)基因按每一列(每个细胞)为单位拼成一列向量，并计算中位数
  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)的概率,每一类别的数量除以总数量
  classProb = getClassProb(classLabelsNum)
  # 命名
  names(classProb) = classLabelsLab
  
  # 计算每个基因的概率,每个基因在各个细胞表达量的总和除以总的细胞数
  geneProb = getGeneProb(xbin$mat)
  
  # 计算每一类细胞(cell cluster)中每个基因的平均表达量，并定义为概率 P( gj | ti )
  condGeneCluster = getGeneConditionalCluster_stats(xbin$mat, classProb, classLabels, cores = cores)
  # 利用贝叶斯公式计算每个基因的后验概率
  clusterPostGene = getClusterPostGene(condGeneCluster, geneProb, classProb)
  # 计算特异性分数
  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几个函数：

# 1 计算稀疏矩阵有表达(有数值的,没有数值的不纳入计算)基因表达量的中位数
binarize = function(x, method = "median") {
  # 计算稀疏矩阵有表达(有数值的,没有数值的不纳入计算)基因表达量的中位数
  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)的概率,每一类别的数量除以总数量
getClassProb = function(x) {
  probs = table(x) / length(x)
  
  return(probs)
}

# 3 计算每个基因的概率,每个基因在各个细胞表达量的总和除以总的细胞数
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 计算特异性分数
getSpecScore = function(postMat, condMat) {
  # 计算特异性分数
  specScore = exp(log(postMat) + log(condMat))
  
  colnames(specScore) = colnames(condMat)
  specScore = as(specScore, "dgCMatrix")
  return(specScore)
}

注意这里的特异性分数的计算，由于变量condGeneCluster代表 P( g_j | t_i )，因此specScore = exp(log(postMat) + log(condMat))代表分数s_{i j}

2. 基于香农熵选出marker基因

library(mclust)
getMarkers = function(gs, quant = 0.99, mutualInfo = FALSE, classEnt = FALSE) {
  
  gs = sg
  quant = 0.99
  mutualInfo = FALSE
  classEnt = FALSE
  
  # 对 sij 分数进行标准化
  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))
}

几个函数：

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)
}