–https://doi.org/10.1038/s41588-024-01682-1

Tissue-specific enhancer–gene maps from multimodal single-cell data identify causal disease alleles

研究团队和单位

Alkes L. Price–Broad Institute of MIT and Harvard

Soumya Raychaudhuri–Harvard Medical School

研究简介

SCENT（单细胞增强子靶基因图谱）是一种用于分析10x Genomics单细胞多组学ATAC + RNA 多模态数据的模型，揭示增强子染色质可及性与基因表达之间的关联。该模型原理是泊松回归和非参数自举法，绘制出增强子-基因对，并预测与表达相关的增强子可能具有重要功能。

增强子-基因图谱的构建：

SCENT应用于9个多模态数据集，包含超过120,000个单细胞或细胞核，构建了23个细胞类型特异性的增强子-基因图谱。

这些图谱在表达定量基因座（eQTL）和1,143种疾病及性状的全基因组关联研究（GWAS）中高度富集了致病变异。

功能验证：

1.SCENT模型认定的增强子在统计上精细定位的因果变异中表现出富集。

2.SCENT模型结果通过了CRISPR-Flow FISH和H3K27ac HiChIP实验验证。

3.SCENT预测的增强子-基因连接具有较高的功能性和进化保守性。

4.模型成功识别了常见疾病和罕见疾病的潜在致病基因，并将体细胞突变热点与目标基因联系起来。在类风湿性关节炎、1型糖尿病和炎症性肠病等自身免疫性疾病中，SCENT揭示了多个潜在的因果变异和基因。

5.与其他方法的比较：SCENT在识别因果变异方面优于现有的单细胞多模态方法（如ArchR和Signac）和基于bulk组织的方法（如EpiMap和ABC模型）。

模型方法

1.数据处理

单细胞转录组数据：标准质控流程后，基因表达 count 矩阵使用NormalizeData函数进行归一化，使用ScaleData函数进行缩放，并使用Harmony进行批次校正。

单细胞ATAC数据：标准质控流程后，ATAC峰矩阵被二值化，如果计数>0则为1，否则为0。

什么peak对应什么基因呢？

研究是在利用基因组位置，把和基因本体位于同一染色体上，且距离基因本体500kb以内的peak划分为该基因的peak。

2.SCENT方法

泊松回归

SCENT方法原理是泊松回归，核心目标是找到染色质可及性（ATAC-seq数据）和基因表达（RNA-seq数据）之间的关系。具体来说，想探究某个染色质区域（peak峰）的可及性是否会影响某个基因的表达。公式的主要内容如下：

${\beta }_0$ 是偏置项；${\beta }_{peak}$ 、${\beta }_{mito}$ 、${\beta }_{nUMI}$ 、${\beta }_{batch}$ 都是对应的系数；

$E_i$ 是原本测序的基因表达值；

${\lambda }_i$ 是模型预测出来的基因表达值；

$X_{peak}$ 是该基因对应的多个peak值（0/1）；

$X_{mito}$ 是细胞种线粒体基因表达的比例；

$X_{nUMI}$ 是细胞种唯一分子识别符的数量；

$X_{batch}$ 是实验批次效应。

模型训练的目标是找到最佳的参数（${\beta }_0$、${\beta }_{peak}$ 、${\beta }_{mito}$ 、${\beta }_{nUMI}$ 、${\beta }_{batch}$），使得模型的预测值${\lambda }_i$尽可能接近实际观察值$E_i$。

模型拟合：

使用最大似然估计来拟合泊松回归模型。

评估显著性：

为了评估 βpeak 的显著性（即染色质峰对基因表达的影响是否真实），SCENT使用了自举法（Bootstrapping）。自举法的核心思想是通过重采样来模拟数据的分布。

具体来说：

从原始数据中随机抽取细胞（有放回），生成一个新的样本；

在新样本上重新拟合模型，得到一个新的 ${\beta }_{peak’}$

重复这个过程很多次（比如100到50,000次），得到 ${\beta }_{peak’}$ 的分布

通过比较 ${\beta }_{peak’}$ 的分布与零假设（${\beta }_{peak’}$=0）来判断 ${\beta }_{peak}$ 是否显著。

多重检验校正：

由于SCENT同时测试了大量的峰-基因对，可能会产生假阳性结果。因此使用了FDR校正来控制假阳性的比例。

研究结果

1.SCENT模型概览

通过统计方法（如泊松回归）来测试每个峰的可及性是否与目标基因的表达水平显著相关。如果某个峰的可及性与基因表达显著相关，那么这个峰可能是一个功能性的调控元件，能够影响该基因的表达。

作者在文章中认定这些peak是增强子，所以文章中以SCENT增强子去称呼相关性系数高的peak。

SCENT精确地识别了调控区域染色质可及性与单细胞中单个基因表达之间的显著关联。

然而，对于高表达和分散的基因，泊松回归可能不是最优选择：

通过置换细胞条形码，可以打乱ATAC-seq和RNA-seq数据之间的对应关系，从而破坏染色质可及性与基因表达之间的真实关联，等于制造出空的数据集。

在空数据集中，由于ATAC-seq和RNA-seq数据之间的关联已经被破坏，理论上不应该检测到任何显著的关联。然而泊松回归模型在空数据集中仍然检测到了大量显著关联，说明模型的I型错误率过高，即假阳性率失控。

为了解决I型错误失控的问题，作者在后续的分析中引入了双尾非参数自举法（nonparametric bootstrapping），通过重采样细胞来估计统计量的分布，从而更准确地评估显著性。

2.发现新的细胞类型特异性SCENT增强子-基因关联

数据包含9个单细胞多组学数据集，涵盖了13种细胞类型（免疫相关、造血、神经元和垂体）。其中作者从11名类风湿性关节炎患者和1名骨关节炎患者的滑膜组织中自测了一个炎症组织数据集（30893个细胞），并获取了8个公共数据集，包含129672个细胞。

质量控制（QC）后，分析了16,621个基因和1,193,842个顺式开放染色质峰（共4,753,521个峰-基因对，每个基因中位数为28个峰）。

SCENT模型用于每个细胞类型（ncells>500），最后一共构建了23个细胞类型特异性增强子-基因图谱，共包含87,648个峰-基因关联（错误发现率（FDR）<10%）。基因的关联峰数量各不相同（从0到97个，平均4.13个）。

ATAC-seq片段数量越多，可能意味着更高的染色质开放区域检测灵敏度，从而有助于更准确地识别与基因表达相关的染色质区域；RNA分子数量越多，有助于更准确地识别与染色质开放区域相关的基因：

为了评估SCENT选出的peak是否具有功能，作者开展以下几点分析：

• （1）它们是否与传统顺式调控注释共定位；
• （2）它们对表达的影响是否在更接近的峰-基因对中更强；
• （3）它们是否具有较高的序列保守性；
• （4）峰-基因关联是否有实验支持。

2.1 测试SCENT选出的peak峰与ENCODE cCRE的重叠mapping情况

平均98.0%的SCENT峰与cCRE重叠，而随机大小匹配的顺式区域和非SCENT峰的重叠率分别为23.3%和89.0%

还使用18种状态的chromHMM结果对免疫细胞类型中的SCENT峰进行了注释；97.4%的SCENT峰在41个免疫相关样本中与启动子或增强子注释重叠。

2.2 检查增强子-基因关联的强度

作者基于更强的关联会更接近目标基因的转录起始位点（TSS）的假设。随着SCENT峰接近TSS，回归系数βpeak（峰可及性对基因表达的影响大小）变得更大且更正向。

2.3 评估了SCENT峰是否具有更高的phastCons分数

更高的phastCons分数，反映了跨物种的序列保守性；进化保守的调控区域在功能上活跃且富含复杂性状的遗传性。

结果如下图，外显子（Exonic）区域比所有非编码顺式区域进化保守性更高。

SCENT调控区域相对于非编码顺式区域也具有保守性。相比之下，所有顺式ATAC峰与非编码顺式区域的ΔphastCons分数较低。

2.4 测试SCENT峰的目标基因是否富集于实验验证的增强子-基因关联

作者使用了CRISPR-Flow FISH结果，其中包括278个阳性和5,470个阴性增强子-基因关联。与实验相关的组织中的SCENT峰显著富集于阳性关联。

其次，使用了naive T细胞、Th17 T细胞和调节性T细胞中分辨率高达1 kb的H3K27ac HiChIP数据。在我们的6个T细胞数据集中，SCENT峰-基因关联在H3K27ac HiChIP增强子-基因环中富集1.6倍。

功能丧失不耐受概率（pLI）和功能丧失观察/预期上限分数（LOEUF）分析

作者假设，具有最多SCENT峰的基因可能是功能上最受限且对功能丧失突变最不耐受的基因。

研究基于人类群体中缺乏有害变异的概率评估了突变约束性，包括功能丧失不耐受概率（pLI）和功能丧失观察/预期上限分数（LOEUF）。每个基因的SCENT峰数量与基因的平均约束分数显著相关（pLI的β=0.37，P=4.9×10−90，分数越高表示约束性越强；LOEUF的β=−0.35， P=−0.35×10−106，分数越低表示约束性越强）。

研究结果与之前的报告一致，即具有大量调控区域的基因在批量表观基因组数据中富集于功能丧失不耐受基因。

3.SCENT峰中eQTL推定因果变异的富集

为了检测SCENT峰是否在统计学上精细定位了表达数量性状位点（eQTL）的推定因果变异，研究使用了GTEx中49种组织的eQTL数据，并将后验包含概率（PIP）>0.2的变异定义为推定因果变异。

结果显示顺式区域内由ATAC-seq定义的所有可及区域在精细定位变异中均显示出2.7倍的适度富集。SCENT峰在所有数据集中平均显示出9.6倍的精细定位变异富集。

由于许多SCENT峰靠近TSS区域，作者考虑了这种富集是否由TSS邻近性驱动。结果显示TSS距离是因果eQTL变异的最重要特征之一：

此外，SCENT方法能够以细胞类型特异性的方式检测顺式调控区域，并通过整合多个数据集构建了四个主要细胞类型的增强子-基因图谱。而且，细胞类型特异性的SCENT增强子在GTEx相关样本中的eQTL变异中富集程度最高：

这些结果展示了SCENT以细胞类型特异性的方式优先考虑因果eQTL变异的能力。

4.SCENT增强子中GWAS因果变异的富集

SCENT可用于通过将其应用于疾病相关组织的多模态数据来构建疾病特异性增强子-基因图谱。

研究检测了SCENT峰是否可用于优先考虑疾病因果变异。从两个大规模生物库（包括FinnGen中的1,046种疾病性状和UK Biobank中的35种二元性状及59种数量性状）中获取了PIP >0.2的候选因果变异。

结果可见SCENT增强子在FinnGen和UK Biobank中的GWAS因果变异中显著富集：

研究从9个数据集和23种细胞类型中构建了SCENT，仅使用了28个样本，远少于构建EpiMap的833个样本和组织以及构建ABC模型的131个样本和细胞系。

尽管数据集较小，SCENT峰在给定数量的峰-基因链接中显示出比ABC模型和EpiMap更高的GWAS变异富集，而且更严格的PIP阈值增加了富集：

研究假设SCENT识别的推定因果变异能够调控染色质可及性。如果是这样，SCENT增强子与caQTL的交集可能会进一步富集GWAS因果变异。

为了验证这一点，研究使用了具有基因型数据的单细胞ATAC-seq样本，并通过利用等位基因特异性染色质可及性或结合等位基因特异性和个体间差异进行caQTL定位。

分析结果观察到SCENT与caQTL交集区域中的富集高于仅使用SCENT的区域。随着caQTL峰阈值的严格性增加，富集度最高可达333倍。

例如，15q22.33处的一个哮喘GWAS位点包含一个位于SMAD3基因内含子内的SCENT增强子，该增强子携带一个推定因果变异rs172936320。该SCENT增强子具有显著的caQTL效应，来自等位基因特异性效应和染色质可及性的个体间差异。

替代等位基因T降低了染色质可及性，并被报告破坏了保守的AP-1结合位点。等位基因T还降低了SMAD3的表达。SCENT识别的目标基因SMAD3参与TGF-β信号传导，该信号在哮喘中重塑气道。

5. 通过SCENT定义GWAS位点的机制

5.1 研究利用SCENT方法定义疾病的因果机制

研究分析了来自FinnGen、UK Biobank以及类风湿性关节炎、炎症性肠病和1型糖尿病的GWAS队列的精细定位变异。

SCENT将4,124个推定因果变异（PIP >0.1）与1,143个性状中的潜在目标基因联系起来。大多数目标基因靠近因果变异，其中20%是距离因果变异最近的基因：

但30.6%的SCENT关联基因距离因果变异超过300 kb。

5.2 自身免疫位点的分析

由于SCENT轨迹主要来源于免疫细胞类型，研究首先关注自身免疫位点。

6q15位点：

T细胞特异性SCENT增强子内的单个精细定位变异rs72928038（PIP >0.3）。

相关疾病和目标基因：类风湿性关节炎、1型糖尿病、特应性皮炎和甲状腺功能减退，距离该变异最近的基因是BACH2。

实验验证：在T细胞中将保护性等位基因编辑为风险等位基因的实验证实了该变异对BACH2表达的影响。删除rs72928038的小鼠初始CD8 T细胞更容易分化为效应T细胞。

4p15.2位点：

类风湿性关节炎和1型糖尿病种，候选变异21个，每个变异的PIP值较低（<0.14）。SCENT模型优先考虑的变异是来自炎症滑膜组织的关节炎组织数据集中T细胞和成纤维细胞中的单个变异rs35944082，目标基因是 RBPJ。

RBPJ转录因子对NOTCH信号传导至关重要，NOTCH信号传导与类风湿性关节炎组织炎症相关。Rbpj敲低导致T细胞分化异常并破坏调节性T细胞表型。

5.3 垂体数据的分析

11p14.1位点，与多种妇科性状（子宫内膜异位症、月经过多、卵巢囊肿和绝经年龄）相关，其目标基因是FSHB，在垂体中特异性表达，促进卵巢卵泡生成。

罕见的FSHB变异会导致常染色体隐性遗传的低促性腺激素性性腺功能减退症。来自其他组织的多模态数据和bulk组织方法未能优先考虑这一变异，因为它们遗漏了最相关的疾病组织——垂体。

6. SCENT增强子中的罕见疾病变异和体细胞突变

6.1 罕见非编码变异的分析

目前，仅在约30-40%的孟德尔疾病患者中能够识别出因果突变，许多病例中的变异被注释为意义未明的变异（VUS），尤其是非编码变异。

研究检查了ClinVar报告的临床非良性非编码变异（总计400,300个变异）与SCENT增强子的重叠情况，ClinVar变异富集中，SCENT增强子中平均包含的ClinVar变异数量是相同基因组长度ATAC区域的2.0倍：

此外，ClinVar变异的密度分别是ENCODE cCREs和所有非编码区域的3.2倍和12倍。SCENT为33,618个非编码ClinVar变异定义了3,724个目标基因。

研究发现了40个与LDLR基因相关的非编码变异，这些变异导致家族性高胆固醇血症1；3个与IL10RA相关的非编码变异，导致常染色体隐性早发性炎症性肠病28（如下图所示）；以及一个与ATM基因相关的内含子变异rs1591491477，导致遗传性癌症易感综合征。

6.2 体细胞突变的分析

在白血病的17个非编码突变热点中，来自血液相关细胞类型的SCENT增强子包含了12个热点。SCENT增强子-基因链接将这些热点与已知的驱动基因（例如，BACH2、BCL6、BCR、CXCR4（下图所示）和IRF8）联系起来。

对比ABC模型：在某些情况下，SCENT为这些突变热点提名的目标基因与原研究中使用的ABC模型不同。SCENT将白血病中chr14:105568663-106851785的一个体细胞突变热点与免疫球蛋白重链相关基因（如IGHA1）联系起来，这可能比ABC模型提名的ADAM6更具生物学相关性。

总结

1.SCENT方法优势

1. I型错误控制：SCENT模型展示了良好的I型错误控制能力，优于常用的统计模型。

2. 因果变异富集：SCENT映射的增强子在eQTL和GWAS中显示出高富集的推定因果变异，并优于之前的方法。

3. 样本效率：尽管使用的样本数量显著较少，SCENT定义的增强子在因果变异富集方面与基于大块组织的方法相当甚至更高。

4. 单细胞水平建模：SCENT在单细胞水平上进行建模，避免了通过将细胞聚合成单个样本来模糊关联。

2.SCENT的局限性

1. 增强子数量较少：与其他资源相比，SCENT增强子-基因图谱中的增强子数量相对较少（图2a）。

2. 扩展潜力：细胞数量与显著SCENT峰-基因链接数量之间的线性关系（补充图8）表明，将SCENT应用于更大的数据集将扩展当前的增强子-基因图谱。

3. bulk组织方法的限制：基于大块组织的增强子-基因图谱可能更难以扩展，因为需要更多的样本。

4. 因果机制的局限性：SCENT专注于基因顺式调控机制以精细定位疾病因果等位基因。可能存在其他因果机制，例如通过反式调控效应、剪接效应或转录后效应起作用的等位基因。为了证明SCENT优先考虑的等位基因的因果关系，需要使用基因编辑技术进行实验验证。

5. 计算强度：由于泊松回归和自举法的使用，SCENT的计算强度高于之前的方法（例如，SCENT需要1.5×10^7 CPU秒，Signac需要2.5×10^5 CPU秒，ArchR需要2.2×10^2 CPU秒）。

6. 优化：在SCENT中实现了多线程和并行化选项，可以线性加快计算速度，但需要额外的计算资源。对于非常大的数据集，算法改进（如下采样或聚合细胞）可能是有用的。

3. SCENT的实用性

1. 疾病组织相关图谱：SCENT的真正实用性在于它能够构建与疾病组织相关的增强子-基因图谱。

2. 数据获取：多模态单细胞数据可以从广泛的人类原代组织中获取，包括那些难以解离的组织，因为这些数据可以从核材料中获取而无需组织解离。

3. 应用示例：FSHB位点：理解妇科性状中的FSHB位点特别需要垂体图谱；RBPJ位点：类风湿性关节炎中的RBPJ位点特别需要滑膜组织图谱。