通过结合活性对 CCCTC 结合因子（CTCF）结合位点变异进行功能注释和优先级排序

本文构建框架识别并优先排序 CTCF 结合位点（CBS）变异，为研究非编码变异提供思路。

　　

一、引言

CCCTC 结合因子（CTCF）通过在数千个基因组位点结合其 DNA 序列基序发挥作用，这些 CTCF 结合位点（CBSs）能够通过调节基因组的三维（3D）组织来促进或抑制基因表达。CBSs 的遗传变异是表型变异的重要驱动因素，在癌症以及孟德尔疾病相关的致病基因表达中频繁出现。然而，目前仅有少数 CBS 变异得到了充分的表征，对于 CBS 破坏在人类群体中的功能影响也缺乏共识。

此前的功能研究在单基因和复杂疾病中对 CBS 变异进行了分析。针对特定基因座的研究表明，CBS 破坏与多种疾病风险相关，如罕见肢体畸形、哮喘和精神分裂症等。小规模测序研究也发现 CBS 变异在包括癌症和阿尔茨海默病等复杂性状中存在负担。但这些研究存在局限性，如特定表型和小样本量的限制。

大规模调查对于 CBS 变异在疾病病理中的功能作用存在不同的结论。一些研究发现 CTCF ChIP 足迹中的变异没有显著富集，而另一些研究则发现与 CBS 破坏相关的全基因组关联研究（GWAS）单核苷酸多态性（SNPs）与某些复杂性状存在显著关系。

由于检测活性 CBS 的困难以及缺乏指导变异注释的遗传密码，大规模优先排序 CBS 变异并非易事。染色质免疫沉淀测序（ChIP-seq）虽能直接检测 CTCF 结合，但存在检测跨度大、需大量生物样本等问题。基于精确权重矩阵（PWMs）检测 CBSs 则存在较高的假阳性率。目前也没有明确的框架将不同细胞类型中检测到的结合模式转化为有用的变异注释，并且在体内提取对转录因子（TF）结合可能有功能影响的变异也具有挑战性。

为解决这些挑战，本文提出了一种用于在全基因组测序（WGS）中大规模识别和优先排序 CBS 变异的分析框架。该框架利用 ENCODE 候选顺式调控元件（cCREs）的数据定义 CBSs，将 CTCF 的跨细胞类型结合模式量化为结合活性这一功能注释，并利用结合活性为 gnomAD v.3 中超过 100 万个单核苷酸变异（SNVs）的等位基因结合测量提供背景信息。

二、材料和方法

1. 用 CTCF 结合活性注释代表性 DNase 超敏位点（rDHSs）：从 ENCODE-SCREEN 网络门户下载所有 rDHSs 的 CTCF Z 分数信号矩阵，筛选出在至少一个生物样本中有显著 DNase 信号的 rDHSs（对应完整的人类 cCREs，n = 1,063,878），屏蔽原始 CTCF 结合信号为零（Z 分数等于 - 10）的 rDHS - 生物样本组合，使用 Stouffer 的 Z 分数元分析（公式为 Σzi/N，其中 zi 代表给定生物样本的 Z 分数，N 是该 rDHS 的非屏蔽生物样本总数）生成结合活性注释。
2. 用结合活性注释 CTCF 基序并定义活性分位数：从 JASPAR 数据库下载构建 hg38 中与 CTCF 规范基序（MA0139.2）匹配的所有参考（REF）序列，排除非常染色体和性染色体重叠群、报告的装配间隙、ENCODE 黑名单区域以及蛋白质编码外显子（Gencode v.44）中的序列匹配，最终得到 1,870,772 个序列匹配。将其与 ENCODE 的 rDHS 区域相交，为每个基序分配相交 cCRE 的活性分数，若有多个重叠则取活性分数较高的 cCRE。在某些分析中，将所有注释 rDHSs 的活性分数分布划分为 100 个等大小的箱，以分位数表示结合活性。
3. 比较预测的 CTCF 结合亲和力与结合活性：将每个 CTCF 基序的预测结合亲和力报告为其 PWM 分数，相对于可获得的最小和最大分数。按活性分位数对每个基序进行分组，评估每组的平均相对 PWM 分数，并使用 Pearson 相关性进行统计评估。
4. 比较进化序列保守性与 CTCF 结合活性：使用多种保守性指标评估核苷酸保守性，如 PhyloP100、PhastCons100、GERP++ 和 LINSIGHT 分数。从 UCSC 基因组浏览器下载相关分数的 Bigwig 格式数据，对于在 hg19 构建中生成的分数（GERP++ 和 LINSIGHT），先使用 UCSC LiftOver 工具转换坐标。为具有 rDHS 支持的基序分配保守性分数，计算每个活性十分位数中保守核苷酸的比例，通过随机抽样评估计算的置信度，并使用 Spearman 相关性评估关系和统计显著性。
5. 处理 gnomAD 的遗传数据：下载 gnomAD v.3.1.2 数据库，筛选出 SNVs，去除等位基因计数（AC）和等位基因频率（AF）为零的条目，要求所有 SNVs 质量高（质量控制 [QC]=PASS）且在 gnomAD v.3 中至少一半样本中被覆盖（等位基因数 [AN]>76,000），最终得到 569,860,911 个经过 QC 的 SNVs，并使用 Variant Effect Predictor（VEP）报告每个变异的功能类别。
6. 注释 gnomAD 中所有候选 SNVs 的结合获得和丢失：使用 Bedtools intersect 识别所有与具有 cCRE 支持的 CTCF 结合基序重叠的经过 QC 的 SNVs，使用 R 包 atSNP 计算 ΔPWM 分数和伴随的 p 值，提取变异位点上下游 14bp 的 REF 侧翼序列，报告 REF 和替代（ALT）等位基因的最高对数优势分数，原始 ΔPWM 分数为 REF 和 ALT PWM 分数之差，ΔPWM≤0 的变异为假定的结合获得，ΔPWM>0 的为假定的结合丢失。通过 atSNP 的重要性采样算法计算每个观察到的 ΔPWM 分数的 - log10 p 值排名分数以评估其显著性，对低 p 值（p = ～0）编码虚拟值（10^?7）。
7. 使用缩放的 CADD 分数评估约束：从 gnomAD v.3.1.2 Hail 表中提取每个变异的 PHRED 缩放的综合注释依赖损耗（CADD）分数，按结合获得或丢失以及显著性（p<0.05）对所有变异进行分层，对于每个结合活性十分位数，随机抽样变异 10,000 次，每次迭代计算使用缩放的 CADD 截止值 10 时假定致病变异的比例，并绘制平均值和 95% 置信区间。
8. 实施 MAPS 分数：使用自定义脚本按照 Karzewski 等人的方法实施突变率调整的单例比例（MAPS）分数。首先使用 pySAM 包识别每个 SNV 的三核苷酸上下文，下载 gnomAD v.2 旗舰论文中的三核苷酸突变率数据并与变异的三核苷酸上下文合并。使用所有最严重 VEP 注释为 “同义变异” 的经过 QC 的 SNVs（n = 2,161,831）校准 MAPS 模型，然后预测其他功能变异类别的 MAPS 分数。对于给定类别，计算原始单例比例（单例数（AC = 1）除以变异数），应用校准模型回归出给定平均突变率下该箱中预期的单例数，并计算每个箱的标准误差（SEM）评估置信度。

三、结果

1. CBS 变异功能注释框架：通过 CTCF 的 PWM 和实验结合数据发现基因组中所有可能的 CBSs，将 CTCF 的 PWM 与 ENCODE 的 1,063,878 个 rDHSs 重叠，19%（n = 355,418）的 CTCF 基序落在注释的 rDHS 内，27%（n = 290,793）的 rDHSs 至少包含 1 个 CTCF 基序，将这些重叠数据的交集（n = 355,418）定义为 CBSs。研究发现，改变 CTCF ChIP-seq Z 分数的检测阈值会显著影响被称为 “CTCF 结合” 的 rDHSs 数量，且在 ENCODE 框架定义的最小阈值下，许多结合事件仅在单个生物样本中观察到。
2. 开发跨生物背景的 CTCF 结合活性的综合注释：由于基于 CTCF ChIP-seq Z 分数阈值筛选 CBS 存在可变性和重复性问题，因此将 CTCF 在每个 rDHS 的结合活性通过对其在所有生物样本中的结合活性 Z 分数分布进行元分析（平均 = 2.58，标准差 = 10.14）量化为单个指标，即结合活性。为验证结合活性在区分具有稳健活性的 CBSs 中的作用，研究其与可能表明真实结合序列的特征的关系。结果发现，高结合活性与高质量 CTCF 基序的重叠频率显著相关（Pearson R = 0.70，p<2.8×10^?16），与预测的结合能量也有很强的相关性（Pearson R = 0.63，p = 2.9×10^?12）。此外，结合活性分位数与进化序列保守性显著相关，较高结合活性分位数的基序中保守位置的比例显著更高（Pearson R = 0.35，p<2.2×10^?16），且这种相关性在多种保守性指标中一致。
3. 结合活性与 gnomAD 中变异水平等位基因结合预测的整合：识别出 gnomAD 中破坏 CBS 的所有 SNVs（N = 1,253,229），计算其 ΔPWM 分数并分类为结合获得或丢失，发现结合丢失的 SNVs（n = 1,015,133）比结合获得的 SNVs（n = 238,196）更多。通过 atSNP 包的重要性采样程序计算每个 ΔPWM 分数的检验统计量，确定高置信度的结合获得（n = 63,931）和结合丢失（n = 527,785）SNVs（p<0.05）。将这些预测与结合活性整合，发现高置信度的结合丢失和结合获得的 SNVs 在高结合活性分位数中显著更多（结合丢失 Pearson R = 0.65，p = 2.7×10^?13；结合获得 Pearson R = 0.67，p = 2.1×10^?14）。
4. 表征 CTCF 结合活性丧失的选择特征：有害的遗传变异会被自然选择清除，因此可以用 AF 作为评估注释方案捕获功能变异效用的指标。研究通过量化不同选择约束度量与框架中识别的变异之间的关系来评估注释方法，应用了 PHRED 缩放的 CADD 分数和 MAPS 指标。结果发现，被注释为高置信度结合活性丧失的 SNVs 中，预测致病变异（缩放的 CADD≥10）的比例显著更高（Pearson R = 0.53，p<2.2×10^?16）。在高结合活性（95^th分位数及以上）时，高置信度 ΔPWM 分数的预测致病变异频率名义上高于低置信度分数。

对所有 CBS 变异的原始 AFs 分析发现其在 gnomAD 中极为罕见。计算 MAPS 分数发现，CTCF 结合活性的丧失与 MAPS 显著相关，高置信度 SNVs 的总体效应更强（高置信度 Pearson R = 0.74，p = 2.2×10^?16）。结合获得的 SNVs 在高活性水平下虽也显示出较高的 MAPS 分数，但数量较少，导致个体 MAPS 分数的不确定性更高。与不同功能类别的基因间和基因内 SNVs 的 MAPS 分数相比，低活性 CBS 中结合丢失的 SNVs 的 MAPS 分数趋向于可能的中性变异类别（如同义变异和基因间变异），而高影响的结合丢失变异（活性十分位数≥90）的 MAPS 分数高于错义变异。

四、讨论

随着基因组测序的广泛应用和调控功能的大规模检测，为非编码 SNVs 的系统表征提供了机会。本文利用这些数据量化了 76,156 多个 WGS 样本中 CBS 变异的功能影响，发现 CTCF 结合活性的丧失与高约束水平相关，在人类群体中很少发生。通过该方法生成了功能优先排序的 CBS 变异目录，并为 TF 结合位点变异的可解释注释提供了可扩展的蓝图。

ENCODE 的 cCRE 数据集和 PWMs 的广泛可用性使得 CTCF 结合序列的广泛映射成为可能，但如何利用这些多维数据有效注释变异仍是挑战。本文的方法为 CTCF 基序分配了跨 214 个生物背景的结合活性综合度量，从多方面展示了结合活性的实用性。一方面，高结合活性与进化保守性等功能特征相关，表明结合活性能有效捕获可能具有功能后果的序列。另一方面，结合活性为等位基因结合预测的有效解释提供了必要背景，研究发现低结合活性 CBS 中的高置信度 ΔPWM 分数在群体频率上类似于可能的中性变异类别（如同义变异），这表明大量基于 PWM 测量的 ΔPWM 分数可能是假阳性，进一步的注释工作应结合结合活性的功能测定。

研究观察到对导致 CTCF 结合活性丧失的 SNVs 存在强烈的负选择信号，且对结合活性影响最大的变异的 MAPS 分数高于错义变异。CBSs 在 3D 基因组组织中具有关键功能，如参与拓扑相关结构域（TADs）的形成和绝缘作用，变异可能破坏这些功能，导致基因表达异常。对本研究中涉及的 CBSs 进行进一步的功能表征，将有助于深入理解高影响变异的机制。

研究也存在一些局限性。注释方法需要对 TF 的结合模式进行统一处理和整合，目前这种深度的数据仅适用于少数 TF，如 CTCF。随着数据的扩展，有望注释更多的 TF 和 SNVs。结合活性的度量没有纳入组织特异性活性，可能会低估在研究较少的组织中高度活跃的 CBSs。在整合 CTCF 的 PWM 与结合活性时，对高质量基序的阈值筛选存在偏差，这种偏差对结合获得和结合丢失 SNVs 的鉴定以及两者功能特性差异的影响尚不清楚。此外，研究方法不是通过单一分数注释变异，而是结合结合活性和 ΔPWM 类别评估其影响，如何在单个综合分数中权衡这些指标尚不明确，但这是一个有前景的未来研究方向。

本文的方法虽然专门针对 CBS 变异，但有望扩展到其他 TF。该框架由大规模遗传数据（gnomAD）、实验结合注释和 PWMs 三个部分组成。虽然 CTCF 是目前唯一直接整合到 cCRE 数据集中的 TF，但许多其他蛋白质的结合数据也存在，可以类似的方式进行整理。PWMs 广泛适用于数百个 TF，且简单、计算成本低，适合大规模变异注释。因此，该框架在未来非编码变异注释中具有广泛的应用潜力。

目前存在许多非编码变异注释器，如 CADD、DeepSea 和 CATO 等。CADD 基于变异与许多基因组注释的关系生成致病性分数，但对测量的变异缺乏生物学见解；DeepSea 等等位基因结合预测器需要事先设定输入变异集，在没有本文提供的结合活性等额外信息的情况下，难以确定要靶向的 CBSs 和变异。此外，结合活性与机器学习输出的整合可能会带来益处，但由于缺乏 CBS 变异的金标准功能数据，评估存在挑战。

综上所述，本文提出了一种优先排序 CBS 变异的框架，利用现有的大规模群体基因组测序数据集和广泛的基因调控注释，为识别功能性 CBS 变异提供了帮助，并为未来非编码变异注释奠定了基础。

五、数据和代码可用性

构建研究数据集、进行分析和生成图表所使用的所有代码都可在 Ruderfer Lab CBS_FunctionalAnnotation GitHub 存储库中作为 Snakemake 工作流程获取，该存储库还包含两个总结文件，详细说明了研究中评估的所有 CBS 和变异的活性分数和 ΔPWM 分数。