6.3 差异表达（DE）分析
6.3.1 Bulk RNA-seq
处理Bulk RNA测序数据时最常见的分析之一是识别差异表达基因。通过比较两种条件下发生变化的基因，例如突变型和野生型或受刺激和不受刺激，可以描述变化背后的分子机制。
如DESeq2和edgeR等方法是为Bulk RNA测序而开发的。此外，还有大量使用RT-qPCR进行了验证的RNA-seq数据集。这些数据可用于对DE查找算法进行基准测试，现有证据表明这些算法的性能相当良好。
6.3.2 单细胞RNA-seq
与Bulk RNA测序相比，在scRNA测序中通常没有明确的实验条件。相反，如上一章（6.2）所示，可以使用无监督聚类方法来识别细胞群。一旦确定了分群，就可以通过比较组间方差的差异来找到差异表达的基因（如SC3中的Kruskal-Wallis检验），或者通过成对的方式比较组之间的基因表达。在下一章中，我们将主要考虑为成对比较而开发的工具。
6.3.3 分布差异
与Bulk RNA测序不同，在单细胞实验中，要比较的每一组通常具有大量样本（即细胞）。因此，我们可以利用每个组中表达量的整体分布来识别组之间的差异，而不是仅仅比较平均表达量，这是Bulk RNA-seq的标准。
比较分布主要有两种方法。首先，我们可以使用现有的统计模型/分布，并对每组的表达使用同一模型拟合，然后进行参数检验，或者检验允许某个特定参数根据组变化的情况下，模型是否拟合得更好。例如，在5.6章中，使用edgeR来检验允许不同批次中的平均表达不同是否显著改善了数据负二项模型的拟合度。
或者，我们可以使用非参数检验，它不假设表达量遵循任何特定分布，例如Kolmogorov-Smirnov检验。非参数检验通常将观察到的表达量转换为秩，并检验一个组的秩分布是否与另一个组的秩分布有显著差异。然而，一些非参数方法在存在大量并列值的情况下会失败，例如单细胞RNA-seq数据中的零值的情况。此外，如果参数检验的条件成立，那么它通常比非参数检验更有效。
6.3.4 单细胞RNA测序数据模型
ScRNA-seq数据最常见的模型是负二项分布模型：

> set.seed(1)
> hist(
      rnbinom(
          1000, 
          mu = 10, 
          size = 100), 
      col = "grey50", 
      xlab = "Read Counts", 
      main = "Negative Binomial"
  )

1000个细胞中单个基因的read数呈负二项分布

均值：μ=mu
方差：σ²=mu+mu²/size
它的参数包括平均表达量（mu）和离散度（size），离散度与方差成反比。负二项分布模型非常适合Bulk RNA-seq数据，并且用于为此类数据设计的大多数统计方法。此外，它能够很好地拟合从UMI标记的数据中获得的count分布（Grun等，2014；Islam等，2011)。
然而，由于零值较高，原始负二项分布模型并不适合全长转录本数据。对于这种类型的数据，已经提出了各种零膨胀负二项分布模型（例如MAST、SCDE）。

> d <- 0.5
> counts <- rnbinom(
      1000, 
      mu = 10, 
      size = 100
  )
> counts[runif(1000) < d] <- 0
> hist(
      counts, 
      col = "grey50", 
      xlab = "Read Counts", 
      main = "Zero-inflated NB"
  )

零膨胀负二项分布

均值：μ=mu·(1-d)
方差：σ²=μ·(1-d)·(1+d·μ+μ/size)
这些模型为负二项分布模型引入了一个新参数d，表示零值。基因表达的零值与基因的平均表达密切相关。不同的零膨胀负二项分布模型使用不同的mu和d之间的关系，有些模型可能将μ和d独立地拟合到每个基因的表达。
最后，有几种方法使用基于转录爆发机械模型的Poisson-Beta分布。该模型得到了强有力的实验支持（Kim和Marioni，2013），它与scRNA-seq数据非常契合，但不如负二项分布模型那么容易使用，而且可供构建的方法比负二项分布模型少得多。

> a <- 0.1
> b <- 0.1
> g <- 100
> lambdas <- rbeta(1000, a, b)
> counts <- sapply(g*lambdas, function(l) {rpois(1, lambda = l)})
> hist(
      counts, 
      col = "grey50", 
      xlab = "Read Counts", 
      main = "Poisson-Beta"
  )

均值：μ=g·a/(a+b)
方差：σ²=g²·a·b/((a+b+1)·(a+b)²)
该模型使用三个参数：a为转录激活率；b为转录抑制率；g为转录活跃时转录本产生的速率。差异表达方法可以检测每个参数在不同组之间的差异，或者只测试一个参数（通常是g）。
所有这些模型都可以进一步扩展，以明确考虑基因表达差异的其他来源，例如根据特定的DE算法的批次效应或文库深度。

往期内容：
重生之我在剑桥大学学习单细胞RNA-seq分析——6. 生物学分析（1）
重生之我在剑桥大学学习单细胞RNA-seq分析——6. 生物学分析（2）