科学网-GLIMES：探索和缓解单细胞数据差异表达分析局限性-张俊鹏的博文

GLIMES：探索和缓解单细胞数据差异表达分析局限性

2025-8-17 18:29

阅读：711

GLIMES：探索和缓解单细胞数据差异表达分析局限性

单细胞转录组学的差异表达（DE）分析为理解细胞对内部和外部刺激的细胞类型特异性反应提供了重要见解。虽然有许多方法可用于从单细胞转录组学中识别差异表达基因，但最近的研究对最先进方法的性能提出了重要关切，包括针对单细胞数据定制的方法和在bulk样本中表现良好的技术。随着群体水平的单细胞研究迅速变得更加可行，强大的和准确的分析方法对于获得有意义的结果至关重要。在这种情况下，目前困扰单细胞数据 DE 分析的四个“诅咒”：过多零值、归一化、供体效应和累积偏差，突出了当前工作流程中的各种局限性和概念缺陷。

零值“诅咒”

Bulk RNA 测序提供了每个基因在异质性细胞类型群体中的平均转录输出。由于样本特性，即使中等测序深度也能提供关于数千个不同基因的信息。相比之下，单细胞 RNA 测序数据要稀疏得多，每个样本表达的基因数量较少，且大量基因的 UMI 计数为零。一个基因的零 UMI 计数可能由三种情况之一引起：真实零，表明该基因未表达；抽样零，表明该基因表达水平较低；或技术零，表明该基因表达水平较高，但未被检测方法捕获。尽管越来越多的证据表明细胞类型异质性是 UMI 数据中观察到的零的主要驱动因素，但单细胞领域普遍认为零主要是由于“丢失”基因（即技术零）造成的不具有信息量的技术伪影。

因此，许多单细胞差异表达研究包括旨在去除所谓零膨胀的预处理步骤。几种流行的预处理方法包括：(1) 通过激进地根据基因的零检测率进行特征选择，例如要求至少在10%的总细胞中存在非零值，并将差异表达分析限制在一个较小的基因集中；(2) 插补零并在插补值上执行差异表达分析；(3) 将零显式建模为额外的一个分量，本质上仅在非零值上执行差异表达分析。

然而，如果零是由于无表达或非常低表达而真实的生物学零，则在单细胞 RNA 测序中忽视或校正零会在任何分析之前丢弃数据集中的大量信息。由于未能考虑细胞类型的异质性，零膨胀的预处理步骤（如归一化和插补）可能会将不需要的噪声引入下游分析，包括差异表达分析。讽刺的是，单细胞差异表达分析中最期望的标记——例如，仅在稀有细胞类型中特异表达的基因——可能会被当前用于处理零的预处理步骤所掩盖。

归一化“诅咒”

"归一化"这一术语在基因组学中已被用来指代多种不同的方法。例如，它可以指在测序文库制备过程中校正 PCR 扩增偏差的过程（即文库大小归一化），也可以指跨不同实验批次协调数据的过程（即批次归一化），或指将数据转换为符合正态分布的过程（即数据分布归一化）。这三种归一化方法均已应用于 bulk RNAs-eq 和单细胞 RNA-seq 数据，旨在最小化不必要的技术变异。为 scRNA-seq 数据的 DE 分析选择适当的归一化技术显然对于保持数据的完整性至关重要，但该领域尚未建立明确的金标准，以明确不同归一化方法应执行的特定条件。

文库大小标准化在bulk RNA 测序分析中至关重要，因为在典型的bulk RNA 测序协议中，由于在文库构建过程中 PCR 引入的扩增水平未知，因此无法追踪 RNA 分子的绝对丰度。在这种情况下，标准化专注于估计并随后校正样本特定的尺寸因子。这个过程使bulk RNA 测序能够估计相对 RNA 丰度。标准化后，样本相对于一个公共参考进行校准，导致大多数基因在样本间显示相似的表达水平。当使用bulk RNA 测序数据执行差异表达分析时，基因被分类为上调或下调，基于大多数基因在组间保持不变的假设。虽然基于尺寸因子的标准化技术适用于bulk RNA 测序，但它并不能有效地转化为单细胞 RNA 测序。单细胞 RNA 测序中的协议，如 10X，采用独特的分子标识符（UMI），以区分真实的 RNA 分子和通过 PCR 产生的那些分子。这使 RNA 水平的绝对定量成为可能。不幸的是，基于大小因子的标准化方法（例如每百万读取映射数，或 CPM），将数据转换为相对丰度，抹去了 UMI 提供的有用信息。此外，由于 CPM 标准化数据中发现的均匀分子数并不能准确代表真实表达水平，CPM 标准化数据无法考虑基因对细胞资源的竞争，最终导致差异表达分析结果不理想。

在批次效应标准化中，降维方法识别出在不同批次中具有一致表达模式的基因；这些基因充当锚点，指导数据的对齐和整合。然而，在单细胞 RNA 测序分析中，仅保留高表达或高变异性基因用于估计批次效应和后续整合。因此，与原始 UMI 数据相比，整合后的单细胞 RNA 测序数据集中的基因数量明显减少。

对于数据分布正态化，该领域提供了直接（例如，对数转换）和高级策略（例如，方差稳定转换或 VST）。sctransform是 scRNA-seq 中 VST 的一个显著实现，它采用正则化负二项回归模型，保留 Pearson 残差以供未来的分析步骤，包括差异表达分析。然而，如果底层数据分布与假设模型显著偏离，应用 VST 可能会将偏差引入分析中。

使用总 UMI 计数发现不同细胞类型的文库大小存在显著差异；值得注意的是，巨噬细胞（MP）和分泌性上皮细胞（SE）的 RNA 含量显著高于其他细胞类型（图 1a）。此外，在所有供体中，SE 细胞的平均文库大小均大于肥大细胞（MA）。这些发现与以下认识一致：绝经后输卵管中主要的活跃细胞类型是 MP 和 SE 细胞，而其他细胞类型在绝经后会处于休眠状态。然而，在整合数据中，文库大小分布的差异被缓解，即使在同一细胞类型内部也是如此（图 1a）。虽然整合减少了供体之间的差异，但它掩盖了细胞类型之间的差异。值得一提的是，CPM 标准化使所有细胞类型的文库大小相等；这种标准化可能会掩盖对理解其独特生物学功能至关重要的细胞类型之间的差异。具有细胞类型异质性的样本可能会受到标准化方法的不利影响。在输卵管数据集中，观察到不同细胞类型的 UMI 计数表达模式存在差异；然而，这些差异在插补或某些转换数据集中不太明显（图 1a）。基因表达频率也因细胞类型而异（图 1b）。然而，归一化过程会显著改变非零 UMI（图 1c）和零 UMI 计数的分布。例如，虽然基因频率随原始 UMI 计数增加呈指数下降，但 VST 数据形成了一个更呈钟形的曲线，非零原始 UMI 计数的众数约为 1.5。非零 CPM 归一化数据（乘以 1000）在 0.2 附近达到峰值，并且比 VST 数据更右偏。经过批次整合后，UMI 计数主要低于 5，并且右偏程度减弱。值得注意的是，归一化过程可以使零 UMI 计数获得非零值（CPM 归一化除外）；例如，VST 数据中的零值被调整为负值，并呈左偏分布（图 1d）。相反，整合过程将原始零值转换为紧密围绕零值分布的数值。进一步检查了一个基因的表达分布情况。以 RUNX3 基因为例（图 1e），原始 UMI 计数和 CPM 数据中的分布仍然呈右偏态。相比之下，VST 和整合数据展示了更宽的钟形分布。这些后置数据集（即 VST 和整合数据）中零的处理方式使它们与前者有着本质区别。这种变异性，加上分布偏态的变化，引发了使用标准化值进行差异表达分析的担忧。

图1 归一化对文库大小、零频率和基因计数分布的影响。a 小提琴图展示基于原始 UMI 计数（顶部）和数据整合后（底部）的文库大小，按细胞类型和供体分类。细胞类型缩写：基质（ST）、平滑肌（SM）、纤毛上皮（CE）、分泌上皮（SE）、壁/血管（P/V）、内皮（EN）、淋巴内皮（LE）、T 细胞/自然杀伤细胞（T/NK）、巨噬细胞（MP）、肥大细胞（MA）、B 细胞/浆细胞（B/P）。b 小提琴图说明原始 UMI 数据中基因表达频率（非零计数）的分布。c 直方图表示原始 UMI 数据在各种数据转换下的非零计数分布。d 直方图详细展示原始 UMI 数据中的零计数，比较 VST 与整合数据（其中零值被插补或转换为非零值）。e 直方图显示基因 RUNX3 在不同数据转换下的分布

供体效应“诅咒”

最近的综述指出，许多单细胞差异表达分析方法容易产生错误发现。这主要是因为未能考虑生物学重复之间的差异，通常被称为“供体效应”。在单细胞研究中，供体效应与批次效应混杂在一起，因为来自一个生物学样本的细胞通常在同一实验批次中处理。虽然单细胞包含多个样本的研究会在预处理阶段进行批次校正，但在下游分析中进行差异表达检验时，它们通常不会校正供体效应。

然而，尚不清楚单独的批次效应校正是否足以消除与供体相关的效应。为了解决这个问题，作者们研究了批次校正前后不同来源变异的贡献。使用上述相同的输卵管数据集，进一步将 4553 个 T/NK 细胞分为 20 个亚型，使用 HIPPO(图 2a)。借助典型标记物，识别了特定亚型，包括 NK 细胞、CD4+ T 细胞、CD8+ T 细胞和成熟的初始 T 细胞。然后，专注于在所有供体中都观察到的亚型（图 2b–c）。

为了量化不同来源的变异比例，在多个亚型对中为每个基因拟合线性模型，将细胞类型和供体作为协变量。通过所有配对，整合减少了供体变异（图 2d）。然而，在比较同一种细胞类型的两个亚型（12 vs.13）和不同细胞类型的两个亚型（13 vs. 19）时，观察到细胞类型相关变异的比例下降。这表明整合不仅减轻批次效应，还影响目标表型。重要的是，分析表明，即使实施批次校正后，仍有显著比例的基因表现出与供体相关的效应（图 2e）。由于批次效应通常从主要主成分估计，代表基因子集的共识，因此有可能残差供体效应在某些基因上，甚至所有基因上都持续存在。因此，在进行差异表达检验时，必须考虑供体效应，以避免错误发现并获得准确结果，即使在使用批次效应去除方法后也是如此。

图2 输卵管单细胞数据的聚类和变异分析。a UMAP 可视化案例研究 1 中 HIPPO 识别的 20 个簇。b 典型标记界定特定细胞亚型：簇 9、15 和 19 为 NK 细胞；簇 7、10、11、14、16、18 和 20 为 CD4+ T 细胞；簇 4、6、12 和 13 为 CD8+ T 细胞；簇 8 和 17 为成熟初始 T 细胞。c 捐赠者在 20 个识别簇中的分布。d 比较不同配对和数据集中由捐赠者和细胞类型效应引起的变异比例分析。e 散点图比较因捐赠者和细胞类型效应导致的变异比例，涵盖各种配对和数据来源

在单细胞研究中解决供体效应的一种广泛使用的方法是伪bulk分析。这种方法涉及将来自同一供体的细胞聚集起来，并使用 DESeq2 或 edgeR 等工具对合并数据进行分析。虽然这种方法在某些情况下是有效的，但它具有很大的局限性。通过将供体效应视为固定效应，并假设其影响在所有细胞中都是均匀的，它忽略了样本内部的异质性。这种过度简化可能导致分析过于保守，可能会遗漏重要发现。此外，bulk RNA 测序差异表达工具默认进行归一化，这可能会在单细胞研究的背景下产生前面提到的相同缺点。因此，虽然伪bulk分析可能是有用的，但它需要仔细考虑，因为它可能无法完全解决单细胞研究中供体效应带来的挑战。

累积偏差“诅咒”

在单细胞 RNA 测序分析中，通常采用分层、顺序的工作流程进行聚类和差异表达分析。这种方法可能将偏差从一个步骤传递到下一个步骤；这种累积偏差最终会降低检测差异表达基因的能力。

无监督学习，尤其是聚类分析，在单细胞研究中至关重要。它根据基因表达模式对细胞进行分组，从而促进细胞类型的注释。虽然聚类在 CPM 等标准化值上非常有效，但它本质上根据基因特征的相对贡献重新加权基因特征。因此，聚类提供了关于不同细胞类型间基因表达变异的总体视角。对相对表达的依赖也使得聚类对预处理步骤引入的错误和偏差具有相当强的鲁棒性。

另一方面，差异表达分析在基因层面进行，使用聚类过程中的组标签。来自供体或批量处理的偏差对每个基因的影响可能不同。尽管差异表达分析在技术上遵循聚类——鉴于其对组标签的依赖——但所使用的指标不必对两者完全相同。如果使用处理后的表达水平进行差异表达分析，累积的偏差仍可能导致假发现和/或遗漏差异表达基因（DEGs）。

一种替代范式：单细胞表达研究的广义线性混合效应模型（GLIMES）

为了最小化上述预处理偏差，Wu等人提出一种在实施批次校正、归一化、插补或特征选择之前，对原始 UMI 计数或零比例进行差异表达分析的方法GLIMES（图3，https://github.com/C-HW/GLIMES）。该方法使用广义线性混合模型，保留了数据中的样本特定结构和生物信号。此外，GLIMES可以通过将它们作为固定效应协变量来调整任何潜在的混杂因素，例如年龄、性别或血统。GLIMES能够与其他效应相比，明确考虑生物学重复之间的变异。作者们提供了两个 GLMM 模型，Poisson-glmm 和 Binomial-glmm，分别对 UMI 计数和零比例进行建模。

图3 单细胞分析中既定工作流程与GLIMES比较。左：在当前单细胞分析流程中，从多个供体收集的原始 UMI 计数被整合以消除批次效应，并进行标准化以进行进一步分析。通常在处理后的数据上执行差异表达分析。右：GLIMES直接对原始 UMI 计数执行广义线性混合模型。随机效应可以解释由于样本引起的批次效应。固定效应包含感兴趣组（例如，细胞类型、对照组/治疗组）和其他调整（例如，年龄、性别）。注释的细胞类型可以从现有流程或 HIPPO 算法中获得，该算法根据 UMI 计数的零比例对细胞进行聚类

与现有的基于 GLMM 的软件包（如 Muscat ）不同，GLIMES 将感兴趣组视为固定效应，同时将供体特异性变异和其他批次效应视为随机效应。这一区别很重要，因为 Muscat 主要设计用于跨不同状态或条件进行差异表达分析，其前提是感兴趣组在不同的实验批次中测序。在 Muscat 的实现中，实验单元——定义为供体和感兴趣组的组合——在这种情况下被视为随机效应。

然而，这种设计可能不适用于其他分析，例如跨细胞类型的比较，其中随机效应与研究的结构不一致。此外，Muscat 通过将其作为偏移量纳入其 GLMMs 中的文库大小因子来对计数进行标准化，这更强调相对丰度而不是原始计数。这种方法类似于伪bulk方法的操作，因为 Muscat 在标准化之前将同一供体的计数分组。因此，其性能通常与伪bulk方法相当。

为了评估GLIMES 的性能，严格评估了八种不同的 DE 分析方法：GLIMES 中的两种新方法：Poissonglmm 和 Binomial-glmm；两种传统的伪bulk方法：DESeq2 和 edgeR；以及四种现有的单细胞特异性方法：MAST、Seurat 中的 Wilcox 和两种 Muscat GLMMs（MMvst 和 MMpoisson）。每种方法都在三种情况下应用，即跨细胞类型、组织区域和细胞状态。评估了各种场景，例如组间库大小差异、组内明显异质性以及混杂的批次效应。

Poisson-glmm 对 UMI 计数拟合 GLMM 模型，而 Binomialglmm 对每个基因的零比例拟合 GLMM 模型，并添加供体作为随机效应。Pseudo-bulk DESeq2 应用 VST 和库大小归一化。Pseudo-bulk edgeR 应用库大小归一化。MAST 采用零膨胀负二项模型，使用 log 转换的 CPM 计数，并将细胞检测率作为协变量。Wilcox 检验是非参数的，使用整合归一化计数（Seurat v4）或归一化计数（Seurat v5）。两种 Muscat 模型，MMvst 使用 VST 计数，MMpoisson 使用原始 UMI 计数，均考虑库大小。两种 Muscat 模型均将供体-组组合视为随机效应。

参考文献

[1] Wu CH, Zhou X, Chen M. Exploring and mitigating shortcomings in single-cell differential expression analysis with a new statistical paradigm. Genome Biol. 2025 Mar 17;26(1):58. doi: 10.1186/s13059-025-03525-6.

以往推荐如下：

1. 分子生物标志物数据库MarkerDB

2. 细胞标志物数据库CellMarker 2.0

3. 细胞发育轨迹数据库CellTracer

4. 人类细胞互作数据库：CITEdb

5. EMT标记物数据库：EMTome

6. EMT基因数据库：dbEMT

7. EMT基因调控数据库：EMTRegulome

8. RNA与疾病关系数据库：RNADisease v4.0

9. RNA修饰关联的读出、擦除、写入蛋白靶标数据库：RM2Target

10. 非编码RNA与免疫关系数据库：RNA2Immune

11. 值得关注的宝藏数据库：CNCB-NGDC