10  差异可及性分析：为何不能直接套用 scRNA-seq 的 FindMarkers？

10.1 【地基篇】洞悉第一性原理

当我们拥有了 scATAC-seq 的统一 Peak Set 和每个细胞的聚类身份后，一个核心的生物学问题便是：是什么样的染色质开放性差异，定义了 T 细胞与 B 细胞的不同？换言之，我们需要进行差异可及性分析 (Differential Accessibility Analysis)，找出在不同细胞类型间开放程度显著不同的 Peaks。

scRNA-seq 的老手们会立刻想到一个熟悉的工具：FindMarkers (在 Seurat 或 Scanpy 中)。这个函数，通常基于 Wilcoxon 秩和检验或 t 检验，在识别差异表达基因方面表现出色。然而，将 FindMarkers 直接、不加修改地应用于 scATAC-seq 数据，是一个微妙但后果严重的错误。其根源，再次回到了 scATAC-seq 数据的两大本质特征：极度的稀疏性和近乎二值的信号。

让我们回顾一下。scRNA-seq 的数据是计数值 (counts)，一个基因的表达量可以从 0 到数千，其分布近似于负二项分布。Wilcoxon 秩和检验这类非参数检验，通过比较两组数值的秩次 (ranks) 来判断分布是否有差异，非常适合处理这类数据。

而 scATAC-seq 的数据，在构建了 Peak-by-Cell 矩阵后，其值绝大多数是 0，少数是 1 或 2。这是一个二值化 (binary) 的世界。在这个世界里，“秩次”的概念变得非常模糊和不稳定。例如，比较两组充满了 0 和 1 的数据，Wilcoxon 检验的效能 (power) 会急剧下降。更重要的是，scATAC-seq 的差异，其本质不是开放“程度”（从 10 增加到 100），而是开放“比例”（从 1% 的细胞开放增加到 50% 的细胞开放）。

因此，我们需要一种专门为二值化、稀疏数据设计的统计检验方法。直接套用为连续计数值数据设计的 FindMarkers 默认方法，就像用测量体重的秤去测量身高，工具的适用性从根本上就存在问题。

10.2 【构筑篇】从代码到科学决策

10.2.1 第一步：选择正确的统计检验

幸运的是，主流的 scATAC-seq 分析包已经意识到了这个问题，并提供了更合适的统计检验方法。在 Seurat/Signac 的 FindMarkers 函数中，我们必须将 test.use 参数设置为专门为 scATAC-seq 优化的选项。

最常用的两种检验方法是：

  1. 逻辑回归 (Logistic Regression, test.use = 'LR'): 将染色质的开放状态（开放=1，关闭=0）作为因变量，细胞类型作为自变量，构建一个逻辑回归模型。该模型可以直接评估细胞类型是否与一个 Peak 的开放概率显著相关。它还会将细胞的总测序深度（即文库大小）作为一个协变量（latent variable）纳入模型，从而在检验的同时，有效地校正了测序深度不同带来的技术偏差。

  2. Fisher 精确检验 (Fisher’s Exact Test, test.use = 'fisher'): 这是一个经典的用于比较两组分类数据比例差异的检验。对于每一个 Peak，我们可以构建一个 2x2 的列联表：

细胞类型 A 细胞类型 B
Peak 开放 a b
Peak 关闭 c d

Fisher 检验会精确计算观测到当前这个表格或更极端情况的概率，从而判断 Peak 的开放比例在两个细胞类型间是否存在显著差异。

10.2.2 第二步：执行差异可及性分析

以 Signac 为例，正确的差异可及性分析代码示例如下。请注意，这里的核心决策就是选择了 LR 作为检验方法。

# Signac R package example
# seurat_object 中已包含 scATAC-seq 数据
# 假设我们要比较 "Cluster1" 和 "Cluster2"
da_peaks <- FindMarkers(
  object = seurat_object,
  ident.1 = "Cluster1",
  ident.2 = "Cluster2",
  test.use = 'LR',  # 关键决策！
  latent.vars = 'peak_region_fragments'
)

这里的 latent.vars = 'peak_region_fragments' 就是将细胞的总测序深度（在此处用总的 aken Fragments 数量作为代理）作为协变量纳入逻辑回归模型，这是进行稳健差异分析的关键一步。

10.3 【避坑篇】新手常见的思维陷阱

  1. 陷阱一：默认参数的“陷进”。最致命的陷阱就是不假思索地使用 FindMarkers 的默认参数（通常是 test.use = 'wilcox'）来分析 scATAC-seq 数据。这样做，你得到的结果在统计学上是不严谨的，很可能会引入大量的假阳性（将不显著的 Peak 误判为显著）或假阴性（错过真实的差异 Peak）。

  2. 陷阱二：忽略测序深度的影响。即使选择了一些比例检验方法，如果不校正细胞间的测序深度差异，结果同样会产生偏差。一个测序深度更高的细胞类型，天然就会有更多的 Peak 被观测到开放。逻辑回归通过 latent.vars 参数优雅地解决了这个问题，而其他方法可能需要用户手动进行更复杂的归一化。

  3. 陷阱三：对“logFC”的误读。在差异分析的结果中，通常会有一个 avg_log2FC (log fold change) 的值。在 scRNA-seq 中，它代表了平均表达水平的倍数变化。但在 scATAC-seq 中，尤其是当使用逻辑回归时，这个值的解释需要更加小心。它更多地反映了开放细胞比例的变化，而非信号“强度”的变化。例如，一个大的 logFC 可能意味着一个 Peak 从在 5% 的 T 细胞中开放，增加到了在 80% 的 B 细胞中开放。

10.4 【蓝图篇】构建你的分析框架

在你的 scATAC-seq 分析框架中，差异可及性分析这一步需要建立一个清晰、审慎的决策流程。

  1. 数据类型自省: 首先，明确你面对的是二值化的稀疏数据，而不是连续的计数值数据。这个认知是所有正确决策的前提。

  2. 方法论匹配: 基于数据类型，主动放弃不适用的默认检验方法（如 Wilcoxon），选择专为比例数据设计的统计模型。逻辑回归因其能够同时处理二值数据并校正协变量，通常是首选。

  3. 参数的精义: 将 test.use = 'LR' 和 latent.vars = 'peak_region_fragments' 视为进行 scATAC-seq 差异分析的“黄金搭档”。深刻理解这两个参数的组合，是如何从统计模型的层面，保证了你的分析结果的稳健性和可靠性。

  4. 结果的审慎解读: 在解读结果时，始终将差异的本质理解为“开放细胞比例”的变化。利用火山图等可视化手段，结合 p-value 和 logFC，筛选出那些在生物学和统计学上都最显著的差异可及性区域。

最终，通过严谨的统计推断找到的这些差异 Peaks，才是真正定义了不同细胞类型身份的“表观遗传指纹”。它们是下游进行 Motif 富集分析、寻找关键调控转录因子，以及与 scRNA-seq 数据整合的坚实基础。

探索生命科学前沿，提升实战技能！🔥 欢迎加入「生信实战圈」，获取最新技术干货、实战案例与行业动态。📊 点击关注，与同行一起成长！ #生物信息学 #组学数据分析 #生信案例代码分享 #R语言编程