当前位置：首页 > news >正文

nature genetics | scATAC-seq预测scRNA-seq，识别影响基因表达的新染色质区域

news 来源：原创 2025/6/9 22:42:16

–https://doi.org/10.1038/s41588-024-01689-8

Single-cell multi-ome regression models identify functional and disease-associated enhancers and enable chromatin potential analysis

研究团队和单位

Christina S. Leslie–Memorial Sloan Kettering Cancer Center

Kushal K. Dey–Memorial Sloan Kettering Cancer Center

研究简介

多组学单细胞测序技术

多组学单细胞测序技术，特别是同时检测染色质可及性（scATAC-seq）和基因表达（scRNA-seq）的技术，为计算方法的发展提供了新的机遇。这些方法将增强子与基因联系起来，推断基因调控网络，并基于染色质势能的概念解析发育轨迹。染色质势能理论认为，在细胞分化过程中，基因位点的可及性先于基因表达。染色质可及性的研究把基因组和转录组学连接起来。

现有方法的概述

现有的方法主要利用单细胞中ATAC和RNA的联合测量，识别单个可及区域（如峰或开放染色质域，DORCs）与基因表达水平之间的成对相关性，从而实现增强子-基因的关联。

例如，一些方法使用泊松回归来测试峰可及性与基因表达之间的相关性，并建模批次或细胞特异性协变量。

而标准的scATAC-seq分析方法通过聚合基因启动子附近或基因位点的染色质可及性，将数据转换为类似scRNA-seq的读数，从而推断基因表达值。

SCARlink模型的提出

文章构建单细胞ATAC与RNA联合分析（SCARlink），这是一种在单细胞水平上，针对单细胞/单核多组学数据的基因水平预测模型。

SCARlink能够从单细胞中基因位点的染色质可及性预测基因表达，捕捉基因位点内和周围区域的远端元件共同调控基因表达的事实。与独立评估单个峰-基因关联的方法不同，SCARlink使用正则化泊松回归对tile水平数据（500bp滑窗口划分基因组的可及性）进行模型训练，避免了将数据汇总为峰图谱的步骤。

SCARlink模型的优势与应用

SCARlink模型能够预测基因表达，还可以通过Shapley值识别细胞类型特异性的增强子，即对预测特定聚类或注释中细胞基因表达起重要作用的基因组tile。
文章展示了SCARlink在从染色质可及性预测单细胞基因表达方面优于现有方法，并正确识别了细胞类型特异性的增强子，这一点通过启动子捕获Hi-C（PCHi-C）得到了验证。
使用Shapley值从模型中确定的调控区域在非编码全基因组关联研究（GWAS）和表达数量性状位点（eQTL）变异中显著富集。
文章证明使用基因水平模型对一组发育调控基因进行分析，能够稳健地实现染色质势能轨迹推断方法。

研究结果

1.SCARlink 准确预测基因表达并识别假定的增强子

1.1 模型设计与输入数据

SCARlink 使用正则化泊松回归模型从染色质可及性预测基因表达。染色质可及性以非重叠的 500 bp tile 形式作为输入，默认覆盖基因体上游 250 kb 至下游 250 kb 的区域。

这一范围足以捕捉大多数基因的远端基因间区及内含子增强子，且可根据需要调整。由于基因长度不同，每个基因的输入 tile 数量也不同。

例如，基因 CCR7（长度 11,703 bp）和 ZEB2（长度 140,502 bp）的输入 tile 数量分别为 1,024 和 1,282。模型通过约束回归系数为非负值，专注于识别正向调控元件（如增强子），从而提高结果的可解释性。

1.2 数据集与模型性能评估

SCARlink 被应用于多个稀疏程度不同的多组学数据集，包括：

低稀疏数据集：10X Genomics 的外周血单个核细胞（PBMCs）、骨髓单个核细胞（BMMCs）和发育中的人类皮层。
高稀疏数据集：小鼠皮肤、胰腺和垂体。

对每个数据集的前 5,000 个高变异基因（HVG）运行模型，并根据基因表达向量的稀疏性进行过滤（90% 的值为 0的基因过滤掉）。过滤后，各数据集的基因数量分别为：PBMC（1,250）、BMMC（1,655）、小鼠皮肤（393）、人类皮层（1,201）、胰腺（784）和垂体（1,221）。

数据的UMI计数和TSS读数如下图所示：

评估模型性能：通过 Spearman 相关性将预测的基因表达与保留细胞中的真实表达进行比较。

1.3 与其他方法的比较

SCARlink 与 ArchR 基因评分、ChrAccR 和 DORC 评分等方法进行了比较：

在高覆盖数据集（PBMC、BMMC、人类皮层）中，SCARlink 显著优于 ArchR 基因评分（单侧符号秩检验，P 值均小于 10^-60）：

在稀疏数据集（小鼠皮肤、胰腺、垂体）中，SCARlink 在小鼠皮肤数据集上表现更优（P < 3.7 × 10^-9），而在胰腺和垂体数据集上与 ArchR 基因评分相当：

在小鼠皮肤数据集中，SCARlink 的预测结果与表达的相关性高于 DORC 评分（P < 3.1 × 10^-18），这可能是因为 SCARlink 同时建模了所有区域染色质可及性的影响：

此外，SCARlink 对细胞数量的下采样具有鲁棒性，在大多数数据集中，使用至少 50% 的细胞数量仍能获得可比的预测结果。

1.4 增强子识别与调控分析

SCARlink 不仅能够预测基因表达，还能通过回归系数识别候选功能增强子。例如：

在 PBMC 数据集中，SCARlink 对 ZEB2 基因的建模识别了其基因组位点上的候选增强子：

在小鼠皮肤 SHARE-seq 数据中，SCARlink 识别了 Lef1 基因的远端调控区域，这些区域未被 DORC 注释，但显示出高回归系数，表明其与转录相关：

这些结果表明，SCARlink 能够有效利用所有区域的染色质可及性信息，识别潜在的增强子-基因关联，并在基因表达预测和调控元件识别方面优于现有方法。

2.Shapley 分析识别细胞类型特异性 SCARlink 增强子

2.1 Shapley 值的计算与意义

通过 SCARlink 生成的回归系数可以反映每个 500 bp tile 的染色质可及性在预测基因表达时的总体重要性。

为了进一步量化每个细胞类型中特定 tile 的贡献，研究计算了每个细胞类型的标准化平均 Shapley 值。这一分析使我们能够识别特定细胞类型中调控基因表达的潜在调控区域。

结果显示，预测的调控元件在基因体内或附近（约 25 kb 范围内）最为富集，而在远端区域的出现频率逐渐降低。

2.2 假设验证：SCARlink 预测区域与三维基因组相互作用的关系

基于活性增强子通常与启动子发生物理相互作用以促进转录的理论，研究假设 SCARlink 预测的调控区域会富集与建模基因启动子的三维相互作用。

为了验证这一假设，我们利用造血细胞的 PCHi-C（启动子捕获 Hi-C）数据，该技术通过全基因组启动子诱饵库识别与启动子相互作用的基因组区域。我们使用广义加性模型分析了 PBMC 数据中 T 细胞亚群、单核细胞和 B 细胞的 SCARlink 预测区域，并将其与 PCHi-C 相互作用进行了比较。

2.3 实例分析：HLA-DQB1 基因的 Shapley 值与 PCHi-C 相互作用的关联

以 HLA-DQB1 基因为例，比较了它的 Shapley 值与 PCHi-C 相互作用的关系。

结果显示，远端 tile 中具有 PCHi-C 相互作用的区域显示出比非相互作用 tile 更高的 Shapley 值，尤其是在 B 细胞中，HLA-DQB1 在该细胞类型中高度表达。这一结果表明，SCARlink 预测的调控区域与三维基因组相互作用高度一致。

2.4 统计分析：相互作用 tile 的 Shapley 值显著更高

为了进一步验证这一发现，研究比较了每个细胞类型中高表达基因的 tile 在有和无 PCHi-C 相互作用两种情况下的 Shapley 值。

结果是具有 PCHi-C 相互作用的 tile 的 Shapley 值显著高于非相互作用 tile。

这一结果支持了 SCARlink 预测的调控区域在功能上与增强子相关的假设，并表明 Shapley 值能够有效识别细胞类型特异性的增强子。

3. 预测的增强子富含精细映射的 GWAS 和 eQTL 变体

3.1 基因相关 tile 的筛选与富集分析

首先根据 Shapley 分数的显著性，为每个基因和细胞类型筛选了一组基因相关 tile。

接着针对 GTEx 组织中精细定位的 eQTLs（后验包含概率 PIP > 0.5）和 UK Biobank 中 82 种性状的 17,769 个精细定位 GWAS 变异（PIP > 0.2），对这些 tile 进行了富集分析。

分析在 PBMC、胰腺和垂体数据集中进行。

3.2 GWAS 变异的富集结果

在三个多组学数据集中，SCARlink 预测的基因相关 tile 在前 15,000 个 tile 中显示出 5.5 倍至 7.5 倍的 GWAS 变异富集，表现优于 ArchR 的标准成对峰-基因关联。随着 PIP 阈值的提高，富集程度也显著增加：

此外，对于 82 种性状中的 79%，SCARlink 的富集程度与 ArchR 相当或更高；在启动子近端和远端区域，SCARlink 的富集表现也优于或与 ArchR 相当：

3.3 疾病遗传力分析

通过分层 LD 评分回归（S-LDSC）方法，研究比较了 SCARlink 和 ArchR 在疾病遗传力分析中的表现。结果显示，SCARlink 在 104 种疾病和性状中的遗传力富集比 ArchR 高 1.2 倍：

在 MAF 和 LD 相关注释条件下，SCARlink 的标准化效应大小（τ*）也显著高于 ArchR（τ* = 0.67 vs. 0.58）。联合遗传力分析进一步表明，SCARlink 的疾病信息量是 ArchR 的 1.9 倍，显示出更强的疾病关联性。

3.4 eQTL 变异的富集结果

对于 GTEx 组织中的精细定位 eQTLs，SCARlink 在 PBMC 中前 20,000 个基因相关 tile 中显示出 12 倍至 20 倍的富集；在 FDR < 0.001 的 tile 中富集程度为 10 倍：

在胰腺数据集中，富集程度也达到 15 倍，显著高于 ArchR 的峰-基因关联。此外，SCARlink 在非匹配 GTEx 组织中的 eQTL 富集程度较低，表明其能够识别组织和细胞类型特异性的调控区域：

3.5 条件独立 eQTL 和 RegulomeDB 变异的富集

SCARlink 在胰腺中前 4,000 个基因相关 tile 中对主要 eQTLs 的富集程度为 10 倍至 21 倍，在 PBMC 中的表现也显著优于 ArchR：

此外，SCARlink 在 RegulomeDB 变异（等级 1a，最严格阈值）中的富集程度也高于 ArchR：

在胰腺和垂体数据集中同样表现出更高的富集：

4. 预测的增强子以细胞类型特异性的方式优先考虑与疾病相关的位点

4.1 哮喘相关变异与 CCR7 基因的关联

研究通过 SCARlink 分析了一个与哮喘相关的精细定位变异 rs112401631（chr17:40608272:T:A）。这个变异与淋巴母细胞系中 CCR7 基因的 eQTL（chr17:40600717:G ）共定位（已有文献支持）。CCR7 基因在 T 细胞归巢至淋巴器官中起关键作用，而 CCR7+ 记忆 CD4+ T 细胞已被证明与哮喘严重程度相关。

在 PBMC 数据中，SCARlink 预测 rs112401631 所在的区域与多种 T 细胞亚型中的 CCR7 基因显著相关，包括：

CD8 效应细胞（P = 3.7 × 10^−08）
CD4 记忆细胞（P = 5.9 × 10^−21）
CD8 初始细胞（P = 7.6 × 10^−18）
CD4 初始细胞（P = 2.2 × 10^−17）

虽然 eQTL 所在的区域未被 SCARlink 预测为与 CCR7 相关，但距离 eQTL 783 个碱基的区域在 CD8 初始 T 细胞中显著相关（P = 2.69 × 10^−07）。该区域包含与 eQTL 处于连锁不平衡的变异 rs1358175（R² = 1）。

此外，GWAS 因果变异周围的 10 kb 窗口包含 T 细胞亚群特有的增强子。这些结果提示，SCARlink 可以帮助确定 GWAS 和 eQTL 共定位背后的潜在因果细胞类型。

4.2 IGF-1 相关变异与 BCL2 基因的关联

另一个例子是精细定位变异 rs12454712（chr18:63178651:T:C），它与循环胰岛素样生长因子 1（IGF-1）浓度和 2 型糖尿病（经 BMI 调整）相关。该变异位于 BCL2 基因的内含子增强子中。IGF-1 通过 BCL2 的活性防止细胞凋亡，而 BCL2 编码一种抗凋亡转录因子。

研究发现该变异位于垂体干细胞（PSCs；P = 1.3 × 10^−10）和生长激素细胞（P = 6.3 × 10^−04）的调控区域，可能暗示其在垂体干细胞分化中的作用。IGF-1 水平异常与胰岛素抵抗和 2 型糖尿病风险相关，而垂体前叶的内分泌细胞（如生长激素细胞）在 IGF-1 的产生中起关键作用。

尽管该变异在 PBMC 数据中的分析结果并没有：

SCARlink 正确地将其在 PBMC 中的区域分配为低显著性。此外，GWAS-eQTL 共定位分析显示，eQTL chr18:63179197:G （PP.H4 = 0.9456）位于与 PSCs 中 BCL2 显著相关的区域。

5. SCARlink 预测的基因表达支持染色质潜能分析

5.1 染色质潜能分析的基本思路

研究探索了 SCARlink 识别的调控区域是否在发育过程中先于基因表达变得可访问，从而支持通过染色质潜能推断发育轨迹。简单来说，这种方法就像在细胞嵌入图中画箭头，箭头指向与当前细胞预测表达最相似的邻近细胞的真实基因表达。这些箭头可以帮助我们可视化细胞的分化轨迹。

具体来说，研究为每个“源”细胞计算了一个平滑的 SCARlink 预测基因表达向量，然后找到一组“目标”细胞，其真实基因表达与预测表达最相关。通过计算从源细胞指向目标细胞的平均位置的向量，可以在 FDL 或 UMAP 嵌入图中可视化染色质潜能向量场。

5.2 数据集的基因选择与过滤

在计算染色质潜能时，默认选择了 SCARlink 预测基因表达与真实表达正相关的前 2,000 个高变基因：

小鼠皮肤：434 个基因中过滤掉 19 个（<5%）
BMMC：785 个基因中过滤掉 36 个（<5%）
垂体腺：612 个基因中过滤掉 2 个（<5%）
发育中的人类皮层：1,201 个基因中过滤掉 73 个（6%）

5.3 染色质潜能向量场的表现

SCARlink 的染色质潜能向量场成功重现了小鼠皮肤、BMMC 和垂体腺中已知的分化轨迹：

然而，在发育中的人类皮层中，染色质潜能未能正确识别放射状胶质细胞群为根细胞类型：

进一步分析发现，预测基因表达与真实表达之间的差异在发育轨迹的中间阶段（nIPC/GluN1）最高，随后逐渐减小：

通过层次聚类，我们识别出两个基因簇，其中一个簇显著富集了与胶质细胞分化相关的 GO 通路：

5.4 基因子集分析恢复正确轨迹

对这 470 个与胶质细胞分化相关的基因子集进行染色质潜能分析后，我们成功恢复了正确的发育轨迹：

此外，这些基因的平均预测表达与真实表达之间的差异在轨迹过程中逐渐增加，表明这些位点的染色质在目标基因表达之前开放：

虽然分析证明了染色质潜能作为在多组学数据集中识别分化轨迹策略的实用性，但可能需要预先选择基因子集以获得与已知生物学一致的结果。

此外，研究发现染色质潜能通常在RNA速率表现不一致或失败的情况下识别出发育正确的轨迹：

总结

SCARlink 能够在无需预先计算峰集的情况下，识别细胞类型特异性的基因增强子。通过 Shapley 值分析，能高效解析组织相关 eQTLs 和 GWAS 性状的细胞类型特异性，并计算染色质潜能向量场以追踪发育或分化轨迹。

SCARlink 的设计是一个简单的基因水平模型：使用对数链接函数并约束回归系数为非负的正则化广义线性模型。这种简洁性使得模型能够快速训练和选择，用于预测基因表达，同时高效计算近似的 Shapley 值，从而以细胞类型特异性的方式识别显著相关的基因组区域（tiles）。SCARlink 估算的基因表达还能够从多组学数据中计算染色质潜能。

此外，通过建模正向加性效应，可以获得了一个高度可解释的模型，其中 Shapley 分析识别的显著区域通过染色体构象捕获数据得到了验证，并且在精细定位的 eQTLs 和 GWAS 变异中显著富集。SCARlink 的细胞类型特异性增强子及其与基因的关联可以整合到功能驱动的全转录组关联研究（TWAS）方法中，用于从基因型预测基因表达。

尽管 SCARlink 的广义线性模型非常有效，但更复杂的基因水平模型可能更合适。例如，可以在回归模型中引入区域之间的交互项，甚至使用非线性神经网络架构来完成相同的单细胞基因表达预测任务。

最后，关于 DNA 序列模型在批量表观基因组和 scATAC-seq 数据中的应用已有大量研究，包括在预测批量基因表达方面的应用。在未来可以将 DNA 序列信息整合到 SCARlink 中，跨基因模型共享与每个细胞相关的序列模型，目标是建模增强子中的调控语法及其对目标基因表达的调控影响。

在这里插入图片描述