引言
在生命科学的研究进程中，单细胞 RNA 测序（scRNA-seq）与单细胞染色质可及性测序（scATAC-seq）成为了探索细胞奥秘的关键技术。scRNA-seq 聚焦于转录本，重点捕捉转录 RNA，主要针对约 30,000 个蛋白质编码基因的 3′非翻译区（UTR）上游 100 - 300 nt 区域，目标大小约为 600 万个碱基对。而 scATAC-seq 则借助转座酶 Tn5，剑指整个开放染色质区域，这片区域在细胞基因组中占比 1% - 2%，大小约 3000 - 6000 万个碱基对。
不过，这两种技术在探索基因组的征程中都面临着诸多阻碍，数据稀疏和采样受限的问题尤为突出，而 scATAC-seq 更是深受其扰。scATAC-seq 中的染色质可及性 “峰”，作为数据解读的关键指标，其定义却充满了主观性，这使得不同研究间的数据难以进行有效的比较和整合。这些 “峰” 本质上是测序读数相较于背景水平出现富集的基因组区间，读数丰度高通常暗示着可能存在活跃的候选调控元件（CREs），但峰的位置和宽度极易受到实验方法的影响。在后续的聚类分析、差异可及性（DA）分析以及转录因子基序富集分析等下游研究中，峰内的读数计数是重要的分析单位，然而由于峰定义的不稳定性，导致研究结果的可靠性大打折扣。
为了深入剖析 scATAC-seq 数据处理的困境，研究人员开展了一系列案例研究。他们利用配对的 scRNA-seq 和 scATAC-seq（10× Multiome）数据，以及来自两个独立收集的 10× 外周血单个核细胞（PBMC）数据集的 scATAC-seq 数据，对 Cell Ranger、MACS2 和 MACS3 这三种热门的峰调用程序进行了评估。结果令人震惊，这些程序识别出的染色质可及性峰存在显著差异，而且大多数峰包含多个 CREs，还有许多峰未能准确捕捉到活跃的 CREs。这一发现凸显了传统基于峰的分析方法在 scATAC-seq 研究中的局限性，也为后续新方法的探索埋下了伏笔。
设计
为了攻克 scATAC-seq 分析方法缺乏标准化的难题，研究人员精心打造了 CREscendo 框架。在 scATAC-seq 的研究领域，数据生成和分析方法的多样性如同杂乱的拼图，使得在比较不同细胞类型或实验的染色质可及性数据时，难以得出一致的结论。
为了充分展示基于峰的分析方法存在的技术和方法学局限，研究人员选取了两个来自人类 PBMC 的独立 ATAC 数据集以及一个小鼠皮层数据集展开深入分析。其中，10× Multiome 平台生成的 PBMC 数据集，能够从同一细胞中获取配对的 scRNA-seq 和 scATAC-seq 数据，为研究提供了宝贵的资源。研究人员重点关注 CD14+单核细胞和 CD8+幼稚 T 细胞这两种含量丰富的细胞类型，对 scRNA-seq 的转录组覆盖比例与 scATAC-seq 的基因组覆盖比例进行了直接对比。他们发现，数据稀疏问题对峰调用产生了深远影响，细胞类型之间的差异 CRE 使用情况常常被宽泛的峰所掩盖。
这一结果在后续对 CD14+单核细胞和 CD4+记忆 T 细胞（基于其在数据集中的普遍性以及与 Multiome 数据集的重叠性进行选择）的 10× Chromium X2 scATAC-seq 数据的分析中得到了验证。此外，对小鼠皮层细胞的分析也表明，这种局限性在不同组织来源的 scATAC-seq 数据中普遍存在，有力地证明了这些问题的严重性和普遍性。
结果
- 标准峰调用方法受数据稀疏影响大，难以精准定位单个 CREscATAC-seq 在追求与 scRNA-seq 可比的覆盖度时，面临着巨大的挑战。从理论上讲，由于 scATAC-seq 的目标区域更广，要达到与 scRNA-seq 相似的覆盖度，其文库大小需要接近 scRNA-seq 的 10 倍。然而，在实际研究中，如 10× Multiome 数据集所示，scATAC-seq 的文库大小仅约为 scRNA-seq 的两倍。以 CD14+和 CD8+细胞为例，scATAC-seq 的中位文库大小为 10,000 个片段，而 scRNA-seq 的中位独特分子标识符（UMIs）数量为 5,000 个。这就导致 scATAC-seq 在更广泛的目标区域上覆盖稀疏，每个细胞的平均片段计数较低。
进一步观察发现，scRNA-seq 的 UMI 计数在基因稀疏度方面呈现出 0 到 10 的动态范围，而 scATAC-seq 的峰稀疏度更高，动态范围更窄，98.7% 的峰区域每个细胞的平均片段数少于 1 个。在分析中，常用的将 ATAC 片段计数量化的方法可能会引入误差，因为零值往往可能是由于采样不足造成的，并非真正的染色质不可及。而且，片段计数与峰宽度之间存在显著的相关性，这意味着峰宽度可能会干扰定量分析的准确性。
研究人员深入分析峰大小与注释的候选 CREs（cCREs）之间的关系，他们借助 ENCODE 提供的 SCREEN 中 cCREs 注册表 v.3，其中包含从 1,518 种细胞类型中分离出的信号所识别的 1,063,878 个 cCREs，并进一步细分为候选启动子、近端增强子和远端增强子，还标注了 CTCF 结合状态。在 10× Multiome PBMC 数据集中，Cell Ranger 识别出的峰中位数大小为 659nt，MACS2 识别出的峰中位数大小为 502nt，而增强子的中位数大小为 288nt，启动子为 328nt，可见这些峰普遍大于单个 cCREs。MACS2 识别出的峰中，超过 55.2% 覆盖至少三个 CREs，在转录起始位点（TSSs）附近，这一比例更是高达 75%。
聚焦于 MACS2 峰，研究人员对包含多个 CREs 的峰周围的 Tn5 切割位点频率进行了研究。令人惊讶的是，尽管包含多个 CREs 的峰跨度总计达 53,588,484nt，但只有 25.1% 与 Tn5 切割位点重合，而在这些切割位点中，78.8% 与注释的 CREs 对齐。在分析小鼠皮层 10× Chromium X2 scATAC-seq 数据集时，也出现了类似的情况，包含多个 CREs 的峰中，只有 16.9% 与 Tn5 切割位点重叠，且其中 71.5% 与注释的 CREs 对齐。即使使用更新的峰调用程序 MACS3 分析 PBMC Chromium X2 scATAC-seq 数据集，结果依然相似，MACS3 峰虽然比 Cell Ranger 峰小，但中位数峰宽度仍远超平均 CRE 长度。
综合这些分析结果可以看出，无论使用哪种峰调用方法，scATAC-seq 数据中的染色质可及性峰大多无法精确捕捉到精细的调控元件，而是常常包含相邻 CREs 的区间，特别是在 TSSs 附近，CREs 密集分布的区域。这一结果清晰地揭示了峰调用在识别关键顺式调控元件方面的局限性，也为后续研究指明了改进的方向，即利用调控注释来优化基于峰的分析方法。2. 整合 CRE 注释优化 scATAC-seq 分析，助力解读细胞类型特异性调控特征研究人员推测，借助功能注释，如 ENCODE 的 CRE 注释和基因活性评分，或许能够为染色质可及性峰的解读提供新的视角。基因活性评分作为解读 scATAC-seq 数据的重要指标，通过映射到基因及其启动子（上游 2kb 加上基因体）的所有片段来估算基因的转录活性，分数越高，表明染色质可及性越高，转录潜力越大。
研究人员利用 ENCODE 的基因活性评分和 CRE 注释，对 PBMC multiome 数据集中由 MACS2 识别出的峰进行解读，并借助配对的 scRNA-seq 数据来验证研究结果。他们通过三个具有代表性的例子，展示了 CRE 注释在捕捉调控事件方面的独特优势，这些调控事件是仅依赖峰水平总结和基因活性评分所无法发现的。
在染色体 1 上的一个峰区域，包含六个注释的 CREs，覆盖约 2kb，与 SFT2D2 基因的启动子和部分编码区域重叠。研究发现，CD14+单核细胞和 CD8+幼稚 T 细胞在该区域呈现出截然不同的细胞类型特异性 CRE 使用模式。其中，CRE5 和 CRE6 作为 CD14+单核细胞特有的 CTCF 结合增强子，在 CD8+幼稚细胞中的活性明显较低，而 CD8+幼稚细胞中 CRE2 的活性则较强。来自相同细胞类型的 scRNA-seq 读数进一步证实了 SFT2D2 基因在这两种细胞中的差异表达，这表明 CRE5 和 CRE6 可能在 CD14+单核细胞中发挥增强转录的作用。值得注意的是，这两种细胞在该峰内的平均片段计数相似，SFT2D2 基因的活性评分也相近，这充分说明了仅依靠峰水平或基因活性评分难以准确捕捉到细胞类型特异性的调控差异。
在染色体 1 的另一个区域，一个跨度为 2,274nt 的峰内包含六个 CREs，与 STK40 基因重叠。在这个区域中，CRE6 在 CD14+单核细胞中的可及性显著增加，表现为 Tn5 切割事件的激增，这强烈暗示 CRE6 是驱动 CD14+单核细胞中 STK40 基因表达的关键因素。虽然 CD14+和 CD8+细胞中 STK40 基因的活性评分存在差异，但仅依靠基因活性无法确定具体的调控元件。
在染色体 19 上围绕 VAV1 基因的区域，有五个 MACS2 峰，该区域至少包含四个转录本和 15 个 CREs，每个峰含有 2 - 5 个 CREs。在第三和第四个峰中，CRE7、 - 9、 - 10 和 - 11 呈现出独特的模式，其中 CRE9 在 CD14+单核细胞中具有高可及性，表明它可能是 CD14+单核细胞特异性的增强子。深入研究发现，该区域存在一些与血小板计数、血小板体积和血小板分布宽度相关的全基因组关联研究（GWAS）单核苷酸多态性（SNPs），且这些关联具有显著的 p 值。血小板与单核细胞相互作用形成单核细胞 - 血小板聚集体（MPAs），会改变单核细胞的功能和表型，而这一过程与炎症和血栓形成密切相关。分析表明，CRE9 在 CD14+单核细胞中高度可及，而在其他 PBMC 细胞类型中则不然，这表明其在单核细胞中具有影响血小板表型的特定调控作用。同样，仅关注峰水平的片段计数、读数计数或基因活性评分会忽略这些细微但重要的调控细节。
综上所述，单纯从峰或基因水平解读可及性数据，就如同用一个大网捞细沙，很容易遗漏关于 CREs 的关键信息。而利用 ENCODE CRE 注释则能够为解读 scATAC-seq 信号提供更深入的理解，研究人员将这一方法正式整合到新开发的 CREscendo 框架中。3. CREscendo 框架解析为了挖掘基于峰的分析方法所遗漏的信号，CREscendo 框架应运而生。该框架巧妙地利用调控注释和 Tn5 切割频率，精准捕捉不同细胞类型之间峰内的差异 CRE 使用情况。
对于包含多个 CREs 的峰，CREscendo 借助 ENCODE 的 CRE 注释，将每个峰精细地划分为不同的区域。假设一个峰中有 k 个重叠的 CREs，那么该峰将被划分为 k + 1 个片段，其中 k 个片段分别对应每个 CRE（CRE1，CRE2，…，CREk），还有一个额外的片段用于表示峰内未被任何 CRE 覆盖的碱基对。通过片段文件中每个片段的起始和终止位点，研究人员能够准确识别 Tn5 切割位点，并统计不同细胞类型中每个片段的切割事件数量。
针对每个峰，研究人员构建了一个列联表，用于记录不同细胞类型之间各片段的切割计数情况。随后，运用卡方检验来评估切割频率的统计学显著差异，以此判断是否存在差异 CRE 使用情况。在获得所有峰的检验统计量后，采用错误发现率（FDR）方法对多重比较进行校正，以确保结果的可靠性。此外，为了更深入地剖析每个 CRE 的作用，研究人员还将总体卡方统计量分解为每个 CRE 片段的部分卡方值，从而量化每个 CRE 的贡献。4. CREscendo 挖掘差异 CRE 使用情况，提升结果解读能力研究人员率先将 CREscendo 应用于分析 CD14+单核细胞和 CD8+幼稚细胞（PBMC 10× Multiome 数据）的 MACS2 峰，并与使用 Signac 进行的 DA 分析结果进行对比。在 53,951 个包含多个 CREs 的 MACS2 峰中，有 13,803 个峰的差异 CRE 使用模式未被 Signac 识别。在这些峰中，研究人员发现了一些具有强烈细胞类型特异性调控模式的区域，由于整个峰的变化倍数较小，这些区域在当前基于峰的分析中常常被忽视。例如，CD55、EEF2K、RAC2 和 PDCD4 等基因，在 CRE 水平上表现出显著的信号。
即便 Signac 识别出了 DA 峰，其解释价值也十分有限，因为整个峰的可及性变化往往比较模糊。以靠近 CD248（也称为 TEM1）的一个峰为例，该峰包含四个 CREs，CRE3 和 CRE4 呈现出明显的细胞类型特异性调控模式。然而，由于基因活性和片段计数的差异较小，Signac 将该峰的排名靠后。但在 CREscendo 的分析中，该峰被列为显著峰（卡方统计量为 3,353.924，p < 1e - 16），这表明两种细胞类型在 CRE 使用上存在显著差异。其中，CRE3 在 CD8+幼稚细胞中高度活跃，而 CRE4 在 CD14+单核细胞中高度活跃。ENCODE 注释显示，CRE3 可能是 CD8+幼稚 T 细胞特异性的增强子，相关细胞类型的批量数据也证实了这一发现，在激活的 CD4+和 CD8+αβ T 细胞以及胸腺来源的幼稚 αβ T 细胞中，CRE3 表现出高活性，而在分选的 CD14+单核细胞中，CRE4 活性较高。这与 CD248 的功能相契合，CD248 编码内唾液酸，是一种在幼稚人类 CD8+ T 细胞上表达并调节增殖的基质细胞抗原。
研究人员进一步将差异 CRE 使用分析扩展到 10× Chromium X2 数据集中的 CD14+单核细胞和 CD4+记忆 T 细胞。基于峰的分析还存在一个问题，即不同样本和研究之间的峰及结果缺乏可移植性和可比性。例如，在 Multiome 数据集中靠近 NFKBIZ 的一个峰，在 Chromium X2 数据集中表现为两个不同的峰，且每个峰都包含多个注释的 cCREs。CREscendo 分析显示，所有峰在 CD14+和 CD8+细胞之间都存在显著的差异 CRE 使用情况，但由于不同数据集之间峰注释的差异，给结果的解释带来了挑战。
与 Signac 的 DA 分析相比，CREscendo 识别出了许多 Signac 因整个峰的变化倍数较小而遗漏的具有强烈细胞类型特异性调控模式的区域。例如，在 DIP2A 基因附近的一个峰，CREscendo 将其识别为显著峰（卡方统计量为 21,729.5，p < 1e - 16），表明 CD14+单核细胞和 CD4+记忆 T 细胞在该峰的 CRE 使用上存在显著差异。其中，CRE2 在 CD4+记忆 T 细胞中活性较高，而在 CD14+单核细胞中几乎无活性。ENCODE 注释显示，CRE3 可能是 CD4+记忆 T 细胞特异性的增强子，相关细胞类型的批量数据也支持这一结论，在激活的 CD4+或 CD8+αβ T 细胞以及胸腺来源的幼稚 αβ T 细胞中，CRE2 表现出高活性，而在分选的 CD14+单核细胞中则无活性，凸显了其细胞类型特异性。
研究人员还利用 10× Chromium X2 数据集对小鼠皮层中的 L2/3 IT 细胞（皮层 2/3 层的兴奋性锥体神经元）和少突胶质细胞进行了差异 CRE 使用分析<研究人员还利用 10× chromium x2 数据集对小鼠皮层中的 l2 3 it 细胞（皮层 2 3 层的兴奋性锥体神经元）和少突胶质细胞进行了差异 cre 使用分析。在此次分析中，使用了 37,048 个包含多个 cre 的峰，并与 signac 的 da 分析结果进行对比。结果发现，crescendo 识别出了许多因整个峰变化倍数小而未被 signac 检测到的细胞类型特异性调控区域，像 adcy5、rgs3、sema6d 等基因在 cre 水平上表现出显著信号。在 cell ranger 识别出的包含多个 cre 的峰中，有 1,171 个区域显示出显著的差异 cre 使用情况，其中 25.8% 未被 signac 检测到。而且，crescendo 独有的峰揭示了差异 cre 使用，共涉及 489 个>
进一步比较 L2/3 IT 细胞与星形胶质细胞时，又发现了如 Sybu、Tox2、Nim1k 等在 CRE 水平有强信号的基因。在这组比较中，1,371 个 Cell Ranger 峰显示出显著差异 CRE 使用，18.3% 未被 Signac 检测到，CREscendo 独有的峰影响了 359 个 CRE 。
在上述两组比较中，血管紧张素原（Agt）基因的启动子区域的 CRE 出现了差异活性，在两种比较中均位列差异使用的前 10 个 CRE 中。CREscendo 将其相关事件判定为显著，在 L2/3 IT 与少突胶质细胞比较中，卡方统计量为 1,698.613（p < 1e - 16）；在 L2/3 IT 与星形胶质细胞比较中，卡方统计量为 1,731.345 。相比之下，在 Signac 分析中，Agt 在 L2/3 IT 与少突胶质细胞比较的 25,294 个显著 DA 峰中排名 2,641，在 L2/3 IT 与星形胶质细胞比较的 33,343 个显著 DA 峰中排名 3,731。这一结果凸显了细胞类型间 CRE 使用的显著差异，也体现了 CREscendo 在识别调控模式上的强大能力。根据 ENCODE 的批量 ATAC-seq 和染色质免疫沉淀（ChIP）-seq 数据，Agt 基因的 CRE1 和 CRE2 似乎作为不同的 CRE 发挥功能，它们在不同脑区的活性差异以及转录因子结合谱的不同，都支持了这一点。L2/3 IT 神经元可能依赖 Agt 的特定增强子来调节神经递质、维持突触可塑性或局部产生血管紧张素 II，以支持神经元活动和连接性；少突胶质细胞则可能利用不同的 Agt 增强子，这与其在髓鞘形成、应对神经元信号时减轻氧化应激或炎症的作用相关；而星形胶质细胞作为中枢神经系统中血管紧张素原的主要产生者，其 Agt 活性的差异或增强对于局部产生调节血脑屏障完整性、脑血流量和神经元兴奋性的血管紧张素肽至关重要。
讨论
本研究的分析结果深刻揭示了基于峰的 scATAC-seq 分析方法存在的关键局限性。这种方法描绘出的基因组景观就像是被迷雾笼罩的地图，大型的调控特征虽然能够被识别出来，但细胞类型特异性的 CREs 却如同隐藏在迷雾中的宝藏，难以被精准定位。这种分辨率的不足，不仅降低了下游分析的灵敏度，还使得不同研究中峰定义的严谨性和质量控制面临挑战。由于峰调用的可变性，研究结果缺乏可移植性，解读起来困难重重，严重影响了研究发现的可重复性。
为了应对这些挑战，研究人员大力推荐采用 CREscendo 框架。该框架借助标准化的注释 CRE 参考，对 scATAC-seq 数据进行量化分析，就像给迷雾中的地图加上了清晰的标注，有效减少了因任意峰定义带来的模糊性，成功揭示了基于峰的分析方法常常遗漏的信号。不过，CREscendo 也并非完美无缺，它依然依赖于峰，所以无法完全摆脱基于峰的方法所带来的限制。
尽管存在这些限制，本研究依然坚定地支持从基于峰的 scATAC-seq 量化方法，向基于标准化注释 CRE 参考的方法转变。这一转变不仅仅是技术层面的调整，更是研究理念的重大变革，它推动着该领域朝着像研究基因一样，对具有明确坐标的 CREs 进行更系统的研究。未来研究的一个重要起点，是将 CRE 水平的汇总作为量化 scATAC-seq 信号的指标，并应用于包括聚类和 DA 分析在内的下游分析。展望未来，对于像人类这样注释完善的基因组，采用无峰分析方法是实现 scATAC-seq 分析精准性和可重复性的范式转变。通过直接聚焦于注释的 CREs，研究人员能够获得更具解释性、可重复性和生物学意义的见解，为系统探索基因组的调控密码开辟新的道路。
然而，向以 CRE 为中心的方法转变并非一帆风顺，该领域必须解决诸如注释不完整、计算方法需要全面更新以及研究人员对新方法的接受度等问题。与此同时，基于峰的分析方法在质量控制、探索性发现以及与现有工具的兼容性方面仍然具有不可替代的价值。因此，一种逐步推进的混合过渡方式，结合两种方法的优势，可能是目前最具可行性的前进道路，既能让研究人员充分利用基于 CRE 的方法带来的益处，又能保持基于峰的工作流程的稳健性和灵活性。
研究局限性
虽然标准化能够提高研究的可重复性，但它可能会在一定程度上限制研究人员对新调控元件的探索。因为标准的 CRE 注释可能无法涵盖所有的调控元件，这就导致研究人员在探索新的调控机制时，可能会错过一些非规范或新兴的 CREs，而这些新发现对于深入理解基因调控机制至关重要。
此外，向基于 CRE 的 scATAC-seq 分析方法的转变，在很大程度上依赖于现有数据库（如 ENCODE）的质量和全面性。如果这些数据库存在局限性或偏差，那么基于它们进行的分析也会受到影响，可能会导致研究结果出现偏差，进而影响对研究数据的正确解读。
资源可用性
- 主要联系人：若需要进一步的信息或资源请求，可联系主要联系人 Mengjie Chen，邮箱为 mengjiechen@uchicago.edu 。
- 材料可用性：本研究未生成独特的试剂。
- 数据和代码可用性：研究人员提供了一个 R 包来实现文中讨论的方法，CREscendo 包可在 GitHub 上获取，网址为https://github.com/ChenMengjie/CREscendo 。此外，用于重现本研究所有数据分析的 R 源代码，与教程一同存放在同一 GitHub 目录下，也可在 Zenodo 上获取，链接为https://doi.org/10.5281/zenodo.14788274 。
