差异表达edgeR，limma（下）

在我们看来，用于检查基因表达分析的最重要的探索性图表之一便是MDS图或其余类似的图。这种图表使用无监督聚类方法展示出了样品间的相似性和不相似性，能让我们在进行正式的检验之前对于能检测到多少差异表达基因有个大致概念。理想情况下，样本会在不同的实验组内很好的聚类，且可以鉴别出远离所属组的样本，并追踪误差或额外方差的来源。如果存在技术重复，它们应当互相非常接近。

这样的图可以用 limma 中的 plotMDS 函数绘制。第一个维度表示能够最好地分离样品且解释最大比例的方差的引导性的倍数变化（leading-fold-change），而后续的维度的影响更小，并与之前的维度正交。当实验设计涉及到多个因子时，建议在多个维度上检查每个因子。如果在其中一些维度上样本可按照某因子聚类，这说明该因子对于表达差异有影响，在线性模型中应当将其包括进去。反之，没有或者仅有微小影响的因子在下游分析时应当被剔除。

主要说明一下 edgeR 中的 glmQLFTest , exactTest 以及 limma 中的 voom 这几种获取差异基因的不同方式

Limma包基于线性模型建模。 它最初设计用于分析微阵列数据，但最近已扩展到RNA-seq数据。 根据limma用户指南的当前建议是使用edgeR包的TMM标准化和“voom”转换，其本质上将标准化数据取对数并估计它们的均值 - 方差关系以确定在线性建模之前每次观察的权重。 默认情况下，Benjamini-Hochberg程序用于估计FDR 。

首先先建立design矩阵,在此研究中，我们想知道哪些基因在我们研究的两组之间以不同水平表达。在我们的分析中，假设基础数据是正态分布的，为其拟合一个线性模型。在此之前，需要创建一个包含细胞类型design矩阵。

据显示对于RNA-seq计数数据而言，当使用原始计数或当其被转换为log-CPM值时，方差并不独立于均值。使用负二项分布来模拟计数的方法假设均值与方差间具有二次的关系。在limma中，假设log-CPM值符合正态分布，并使用由 voom 函数计算得到的精确权重来调整均值与方差的关系，从而对log-CPM值进行线性建模。

当 *** 作DGEList对象时，voom从x中自动提取文库大小和归一化因子，以此将原始计数转换为log-CPM值。在 voom 中，对于log-CPM值额外的归一化可以通过设定normalize.method参数来进行。

下图左侧展示了这个数据集log-CPM值的均值-方差关系。通常而言，方差是测序实验中的技术差异和不同细胞类型的重复样本之间的生物学差异的结合，而voom图会显示出一个在均值与方差之间递减的趋势。 生物学差异高的实验通常会有更平坦的趋势，其方差值在高表达处稳定。 生物学差异低的实验更倾向于急剧下降的趋势。

不仅如此，voom图也提供了对于上游所进行的过滤水平的可视化检测。如果对于低表达基因的过滤不够充分，在图上表达低的一端，受到非常低的表达计数的影响，可以观察到方差水平的下降。如果观察到了这种情况，应当回到最初的过滤步骤并提高用于该数据集的表达阈值。

edgeR使用经验贝叶斯估计和基于负二项模型的精确检验来确定差异基因。 特别地，经验贝叶斯用于通过在基因之间来调节跨基因的过度离散程度。 使用类似于Fisher精确检验但适应过度分散数据的精确检验用于评估每个基因的差异表达。edgeR 在默认情况下，执行TMM标准化程序以考虑样本之间的不同测序深度，通过Benjamini-Hochberg用于控制FDR 。

精确检验适用于多组实验的精确统计法，使用函数 exactTest 估计差异基因，人们将其称为classic edgeR。

estimateDisp 函数在一组离散网格点上为每个标签(基因)计算一个似然矩阵，然后应用加权似然经验贝叶斯方法获得后验离散度估计。如果没有设计矩阵，它计算每个标签的分位数的条件似然，然后将其最大化。在这种情况下，它类似于函数 estimateCommonDisp 和 estimateTagwiseDisp 。

同样首先利用 calcNormFactors 对因子进行矫正

似然比检验是用广义线性模型（glms）的统计学方法，适用于不同复杂程度的多因素实验，而 QLFTest 则是首选，因为它反映了估计每个基因离散度的不确定性。当重复次数较少时，它可以提供更可靠的错误率控制。

1. 过滤低表达的基因

仅保留在两个样品或更多样本中CPM>1的基因

CPM=Reads/(total reads in the sample/1,000,000)

问题：会受到测序深度的影响

2. 选择一个the most average sample作为reference sample

    1. 对于每一个sample

        1. 除以size做normalization

        2. 选择top25%的基因表达值（75%的基因表达<=这个值）

    2. 对上述取平均值，最接近这个平均值的sample就是reference sample。

    3. 对于非reference sample，分别计算scale factor

        1. 选择gene list

            1. biased genes：log fold change ，过滤掉top/down 30%

            2. highly/lowly transcribed in both samples：geometric mean，过滤掉top/down 5%

            3. 取交集

        2. 对log2计算加权平均

            1. scale factor = 2 ^ (weighted average of log2 ratios)

     4. 对于所有sample的scale factors，center后得到最后的scale factor

1. 去掉不表达或者只在一个sample表达的基因

        DESeq2查看某个基因所有样本均一化read的平均值

2. Log e （defalut：e为底）

3. 为每一个基因

    1. 计算平均值

    2. 每一个值减去平均值

4. 为每一个sample计算median

5. 将median还原成normal number，得到最终的scale factor

一般而言，样本间的变异系数（coefficient of variance，CV）是由两部分组成的，一是技术差异（Technical CV），另一个是生物学差异（Biological coefficient of variance，BCV）。前者是会随着测序通量的提升而消失的，而后者则是样本间真实存在的差异。所以，对于一个基因g而言，它的BCV在样本间足够大的话，就可以认为基因g是一个差异表达基因。

拟合负二项广义对数线性模型，如果某个基因的表达值偏离分布模型，那么该基因为差异表达基因。

为什么使用负二项分布而不是泊松分布？

overdispersion: 真实数据的分布偏离泊松分布（方差=均值），方差明显比均值大

edgeR和后期的DeSeq2使用负二项式广义模型中的NB2模型

早期的DeSeq2：the variance is determined by the smoothed function f of the mean

BH等

参考资料：

1. [Gene expression units explained: RPM, RPKM, FPKM, TPM, DESeq, TMM, SCnorm, GeTMM, and ComBat-Seq](https://www.reneshbedre.com/blog/expression_units.html)

2. statquest: DeSeq2/edgeR共三个视频

3. [广义典型相关分析_广义线性模型（GLM）概述及负二项回归应用举例和R计算_weixin_39629467的博客-CSDN博客](https://blog.csdn.net/weixin_39629467/article/details/111341514)

4. [17. 负二项式模型 — 张振虎的博客 张振虎 文档](https://zhangzhenhu.github.io/blog/glm/source/%E8%B4%9F%E4%BA%8C%E9%A1%B9%E6%A8%A1%E5%9E%8B/content.html)

---