Nature | 对 6.8万个顺式调节元件进行功能测试，高精度预测 cCRE 功能和变异效应

Basic Information

• 英文标题：Massively parallel characterization of transcriptional regulatory elements
• 中文标题：大规模并行表征转录调控元件
• 发表日期：15 January 2025
• 文章类型：Article
• 所属期刊：Nature
• 文章作者：Vikram Agarwal | Nadav Ahituv
• 文章链接：https://www.nature.com/articles/s41586-024-08430-9

Abstract

Para_01

1. 人类基因组包含数百万个候选顺式调节元件（cCREs），这些元件具有特定细胞类型的活性，影响健康和许多疾病状态。
2. 然而，我们缺乏对控制这些cCREs活性及其特定细胞类型特征的序列特征的功能性理解。
3. 在这里，我们使用基于慢病毒的大规模平行报告测定法（lentiMPRAs）测试了超过680,000个序列的调控活性，这些序列代表了三种细胞类型（HepG2、K562和WTC11）中的广泛注释cCRE集，并发现其中41.7%的序列是活跃的。
4. 通过测试两种取向的序列，我们发现启动子具有方向偏倚，并且它们的200个核苷酸核心作为非特定细胞类型‘开开关’，为它们相关的基因提供相似的表达水平。
5. 相比之下，增强子的方向偏倚较弱，但具有更强的组织特异性特征。
6. 利用我们的lentiMPRA数据，我们开发了基于序列的模型来高精度地预测cCRE功能和变异效应，描绘调控基序并建模其组合效应。
7. 在所有三种细胞类型中测试涵盖60,000个cCRE的lentiMPRA文库进一步确定了决定细胞类型特异性的因素。
8. 总体而言，我们的工作提供了三种常用细胞系中功能性CRE的广泛目录，并展示了如何使用大规模功能测量来解析调控语法。

Main

Para_01

1. 序列变异在顺式调节元件（CREs）中是人类疾病的主要原因。
2. 例如，大多数全基因组关联研究（GWAS）指出，携带远端CREs如增强子的非编码单倍型，在常见疾病中起作用。
3. 然而，预测CREs中核苷酸变异的功能效应仍然具有挑战性。
4. 一个主要限制是缺乏对人类基因组中可能数百万个CREs的全面功能界定，其中许多具有组织或细胞类型特异性活性。
5. 这一障碍也限制了开发能够高精度预测组织特异性CRE活性的机器学习工具的能力。

Para_02

1. 基因组规模的生物化学技术的出现，可以全局性地记录开放染色质区域、转录因子结合位点、组蛋白修饰和mRNA表达水平，为研究数百种细胞类型的基因调控和转录景观提供了框架。
2. 这些努力已经导致在人类基因组中发现了数百万个cCREs。
3. 然而，这些方法主要是描述性的，无法确认任何一个给定的cCRE是功能性的。

Para_03

1. 大规模平行报告测定法（MPRAs）通过多重测试数千个序列或变异体的调控活性来克服这些限制。
2. 先前的研究利用了MPRAs或衍生测定法，自我转录活性调控区域测序（STARR-seq），来测试大量cCREs在人类细胞中的调控活性。
3. 然而，这些测定法依赖于瞬时转染，提供了一个附加体（‘基因组外’）的读数，并且主要局限于可以被稳定转染并在大量培养中生长的已建立细胞类型。
4. 为了解决这个问题，我们开发了一种lentiMPRA，它可以实现可重复性和多重性，扩展到难以转染的细胞系，如神经元或类器官，并提供一个‘基因组内’的读数。
5. 由于慢病毒整合是随机的，lentiMPRA测量的是cCREs的功能效应在不同基因组位置上的平均值。
6. lentiMPRA与ENCODE注释和基于序列的模型的相关性更强，并且比附加体MPRA提供了更高的细胞类型特异性预测。
7. 此外，系统地比较lentiMPRA与其他八种MPRA设计发现，它与附加体MPRA和STARR-seq有很强的相关性，但也存在差异。
8. 然而，lentiMPRA的一个局限性是单次实验中可以测试的序列或变异体的数量。

Para_04

1. 我们应用了一种优化的lentiMPRA方法，并确认了该技术在单次实验中测试超过200,000个序列的重现性和可靠性，这涵盖了任何给定人类细胞类型的大部分cCRE。
2. 我们将此方法应用于显著扩展三个ENCODE细胞类型的MPRA数据，包括人肝细胞（HepG2）、淋巴母细胞（K562）和诱导多能干细胞（iPS细胞；WTC11），以检查启动子和增强子的相对取向依赖性。
3. 此外，我们在所有三条细胞系中测试了60,000个序列。
4. 利用这些数据，我们描述了核心启动子区域的活性效应，并训练了可以预测调控和核苷酸变异活性的模型。
5. 我们确定了与细胞类型特异性活性相关的生化和基于序列的特征，并提供了一个包含数千个功能性cCRE的目录，从而加深了我们对基因调控序列中基因型与表型关联的理解。

Optimization of lentiMPRA

Para_01

1. 为了扩大lentiMPRA的应用范围，我们修订了已建立的协议，在文库扩增步骤中向被检测序列添加随机条形码以及最小启动子（扩展数据图1a）。随后，通过测序重建元件-条形码关联（扩展数据图1b），并使用MPRAflow进行分析。
2. 为了评估这种修订后的lentiMPRA方法的稳健性，我们设计了两个先导文库（补充图1a）。由于位于峰中心的DNase可及性已被证明是预测调控元件和MPRA活性的良好指标，我们将其作为选择cCRE的主要标准。
3. 第一个先导文库包含HepG2细胞中的9,372个元件，包括：（1）9,172个非启动子重叠的DNase高敏感性峰中心的cCRE；（2）50个阳性控制和50个阴性控制的合成序列（即，具有多个已知转录因子结合位点或没有已知结合位点的序列）；
4. 以及（3）50个阳性控制和50个阴性控制的天然增强子（分别观察到具有高和低增强子活性的序列）。
5. 第二个先导文库包含K562细胞中的7,500个元件，包括：（1）6,394个非启动子重叠的DNase高敏感性峰中心的cCRE；
6. （2）通过将CRISPR干扰（CRISPRi）与单细胞RNA测序测量相结合，识别功能增强子-基因对而确定的290个阳性控制和276个阴性控制；
7. （3）从我们的文库中随机选择的cCRE经过二核苷酸洗牌得到的250个阴性控制；
8. （4）50个阳性控制和200个阴性控制的天然增强子（分别观察到具有高和低增强子活性的序列）；
9. 以及（5）在感兴趣的基因座如α-珠蛋白、β-珠蛋白、GATA1和MYC中手动选择的24个阳性控制和16个阴性控制（补充表1）。

Para_02

1. 这些先导的lentiMPRA文库被用于三重转导细胞，并且条形码在DNA和RNA水平上被测序，如文献16所述。
2. 每个元件的活性分数计算方法是将对应于该元件的所有条形码的RNA分子的归一化计数进行对数2转换，然后除以对应于该元件的所有条形码的DNA分子的归一化数量（扩展数据图1c和补充表2）。
3. 我们观察到，在每个重复实验中，每个元件的中位条形码数量范围为50至250（补充图1b），这提供了每个元件大量独立测量的数据点，以及活性分数（即，对数2转换后的RNA/DNA比率）在重复实验之间高度一致（皮尔逊相关系数为0.88至0.96；补充图1c、d）。
4. 平均三个重复实验，元件活性分数的分布在这组大多数正对照和负对照之间表现出强烈的差异（补充图1e）。
5. 然而，对于基于CRISPRi筛选工作的对照组，我们观察到了一个例外，它们显示出比正对照稍弱的信号，表明在我们的报告基因测定和脱离其表观遗传背景的情况下，它们仍然能够激活转录（补充图1e）。

Para_03

1. 我们接下来分析了两个lentiMPRA文库的功能增强子。
2. 我们发现，在8,960个cCREs中有2,740个（30.6%）比阴性合成对照物19（在HepG2细胞中）更活跃，而在6,315个cCREs中有3,703个（58.6%）比打乱的阴性对照物更活跃（在K562细胞中）（5%假发现率（FDR））。
3. 然而，细胞类型之间活性cCREs比例的差异不应直接比较，因为每次实验中的阴性对照物不同——一个代表打乱的对照物，另一个则不是。
4. 鉴于之前对β-珠蛋白基因座调控元件的广泛研究以及这些序列被纳入我们的K562文库中，我们评估了我们的MPRA结果是否重现了先前研究对于五个已表征的cCREs（称为HS1–5）的发现（扩展数据图1d）。
5. 与之前的工作22,23一致，我们观察到HS2相对于HS1和HS3–5强烈激活转录（扩展数据图1d）。
6. 总之，这些初步实验验证了我们修订后的lentiMPRA方案能够高精度且可重复地测量调控活性。

cCRE characterization with lentiMPRA

Para_01

1. 随着我们的试点文库显示出可重复和稳健的结果，我们接下来着手测试我们修订的lentiMPRA方法是否能够在单次实验中测量超过200,000个序列，这些序列构成了任何给定人类细胞类型的cCREs的主要部分。
2. 使用与我们在试点文库中相似的方案16，我们试图测试所有已知的19,104个蛋白编码基因启动子以及潜在增强子（不靠近启动子的DNase峰）的所有方向组合（图1a）。
3. 在HepG2细胞中，我们测试了所有启动子和66,017个潜在增强子；在K562细胞中，我们测试了所有启动子和87,618个潜在增强子；在WTC11细胞中，由于其转导效率降低，我们测试了7,500个启动子和30,121个潜在增强子中的83,201个（图1b和方法）。
4. 为了进一步探究开放染色质是否是转录激活所必需的，我们在K562文库中另外测试了14,918个异染色质区域，这些区域由ENCODE联盟提名，来自GATA1、MYC、HBE1、LMO2、RBM38、HBA2和BCL11A基因位点两侧各1Mb的区域，这些基因是已知的人类疾病相关基因和红系谱系基因。
5. 总的来说，结合双核苷酸洗牌的阴性对照和其他先前研究中确定的阳性对照和阴性对照18,19,24，我们在HepG2细胞中共设计和测试了164,307个元件，在K562细胞中共设计和测试了243,780个元件，在WTC11细胞中共设计和测试了75,542个元件（图1b和补充表3）。

Fig. 1: lentiMPRA in three cell types.

- 图片说明

◉ a, lentiMPRA 策略用于大规模文库。数千个 cCREs 包括潜在增强子（标记为 DNase 峰）和启动子（位于转录起始位点 (TSS) 中心）被插入到报告质粒中，并且带有条形码。 ◉ 这些文库使用慢病毒感染细胞，整合的 DNA 和转录的 RNA 条形码被测序以量化 cCRE 活性。 ◉ b, HepG2、K562 和 WTC11 文库的组成。每个文库包含数千个潜在增强子和启动子、阴性对照和阳性对照。 ◉ 条形码根据测试的方向着色，并且在条形码上方显示了测试元素的数量，并根据元素类型进行着色。 ◉ c, 元素活性的小提琴图，通过 cCREs 和阴性及阳性对照的 RNA/DNA 比率的对数转换来测量。 ◉ 启动子和增强子分布与随机类别进行了比较，采用单侧 Wilcoxon 秩和检验，随后进行 Bonferroni 校正（*P < 10^-8）。

Para_02

1. 我们观察到每种细胞类型中的每个重复样本中有20-50个中位条形码（补充图2a）；至少支持10个条形码的元素（我们的最低阈值）在重复样本中的标准偏差明显降低（补充图2b）。
2. 元素活性评分在重复样本对之间也表现出强烈的符合性，皮尔逊相关系数分别为0.94（HepG2）、0.76（K562）和0.76（WTC11）（补充表4和补充图2c-e）。
3. 在三个重复样本中取平均值后，我们还观察到了在我们的试点和大规模文库中共有的cCREs之间的元素活性评分之间的强相关性（HepG2细胞中的皮尔逊r=0.94和K562细胞中的r=0.81；补充图2f）。
4. 同样地，在HBE1位点可视化大规模K562文库（扩展数据图1d）和其他六个与疾病相关的位点（扩展数据图2和补充图3），证实了与K562试点文库的一致性，大规模文库具有更高的密度，并突出了额外的功能调控元件。
5. 由于文库大小和每个元素测序深度之间的权衡，两个大规模文库的重复样本间相关性低于试点文库。
6. 为了进一步研究这种权衡，我们降低了条形码的数量。
7. 将条形码减少到90%，相对于从完整数据集得出的活性评分，元素活性评分的斯皮尔曼相关性接近完美，而较小比例的条形码则降低了这种相关性（补充图2g）。
8. 考虑更大的条形码减少比例（补充图2h）和更多的重复样本（补充图2i）时，元素活性评分的标准差分布变得更紧密。

Para_03

1. 我们完全处理过的元件活性评分分布，在每种细胞类型中的正负对照之间存在显著差异，大多数启动子和潜在增强子的评分范围介于最大正对照和最小负对照之间（图1c）。
2. 与潜在增强子相比，启动子表现出平均更高的活性评分和双峰分布，在所有细胞类型中潜在增强子呈现出右偏分布（图1c）。
3. 这种双峰分布可能是由于在MPRA中不活跃的启动子表现出很少或没有活性造成的。
4. 接下来，我们分析了所有文库，以实证测量每种元件类型和每种细胞类型中功能性cCRE的比例。
5. 使用每种细胞类型的洗牌对照作为背景集，并将测量的cCRE的两种取向作为前景集，我们发现超过50％的所有启动子序列具有调控活性（HepG2：20,816个中的11,367个（54.6％）；K562：29,376个中的15,362个（52.3％）；WTC11：9,964个中的5,038个（50.6％）；5％FDR）。
6. 对于潜在增强子，我们发现多达42％是活跃的（HepG2：118,433个中的50,714个（42.8％）；K562：169,260个中的69,820个（41.3％）；WTC11：45,942个中的11,861个（25.8％）；5％FDR）。
7. 另外的功率分析表明，我们的检测功能调控元件的能力在条形码降采样后趋于平稳，这表明额外的测序深度不太可能导致估算值发生实质性改变（补充图2j）。

Para_04

1. 为了评估我们潜在的增强子是否可以用于验证显著的增强子-启动子相互作用和/或预测CRISPRi结果，我们将这一组与在三项不同的CRISPRi扰动研究中测试过的那些进行了交集分析。
2. 我们检查了我们的分箱活性评分中有多少比例与显著调节一个启动子相关。
3. 我们观察到，活性最高的分箱相对于低活性分箱，验证的增强子比例几乎翻了一番（补充图4a），这表明我们的MPRA数据集可能被用来预测具有更大CRISPRi效应的序列。
4. 考虑到我们的活性评分与活动接触（ABC）评分相结合，我们在这三个研究中的两个中观察到，我们的评分仅略微提高了区分显著与非显著增强子-启动子相互作用的任务的表现（补充图4b）。
5. 这种性能的小幅提升可能部分可以通过观察到H3K27ac信号来解释，该信号已经在ABC模型中考虑，它具有整合局部和远端调控信息的优势，这可能掩盖了单独考虑元素局部活性的因素。

Promoter properties and orientation effect

Para_01

1. 根据先前研究14中建立的程序，我们利用我们显著扩展的MPRA数据来检查启动子和增强子的相对取向依赖性。
2. 我们量化了cCREs表现出可观察到的取向依赖性的程度。
3. 在所有检查的细胞类型中，与相对于报告者以相反取向克隆的cCREs相比，以相同取向克隆的cCREs在复制对之间的相关性大约高出0.2（图2a）。
4. 在HepG2细胞系（具有最高技术重复性的细胞类型）中，我们在复制之间取平均值后，检测到了大量的cCREs，这些cCREs在一个取向上的活性比在另一个取向上更强（图2b）。
5. 这些发现表明，cCREs的活性在很大程度上，但并非完全，独立于取向。

Fig. 2: Properties and orientation dependence of promoters and potential enhancers.

- 图片说明

◉ a, 蜂群图显示了三个重复样本之间每一对比较的皮尔逊相关值。这些相关性是根据位于相同（正向与正向（FvF）和反向与反向（RvR））或相反（正向与反向（FvR）和反向与正向（RvF））方向上的元件的观察到的cCRE活性值计算得出的。 ◉ b, 正向与反向取向每个cCRE的平均活性评分的散点图。区域根据数据密度从浅蓝色（低）到黄色（高）着色。显示了皮尔逊（r）和斯皮尔曼（ρ）相关值。 ◉ c, 显示每个细胞类型中启动子和潜在增强子的链不对称分布的箱线图。中心线是中位残差值，箱子边缘界定第25和第75百分位数，使用单侧Wilcoxon秩和检验进行评估。 ◉ d, 热图表示启动子和内源基因表达水平（以每百万转录本计（TPM））在RNA测序中的正向（s）和反向（as）取向之间的相关性。圆圈的大小与皮尔逊相关性成比例。 ◉ e, HepG2细胞中正向取向启动子的活性评分与内源基因表达水平的散点图。 ◉ f, g, 表示MPRA活性排名前1,000与后1,000的启动子（f）和潜在增强子（g）中HOCOMOCO v.1230转录因子家族富集的火山图。富集（以优势比衡量）和Benjamini–Hochberg校正的q值通过费舍尔精确检验计算得出。显著家族（超过P值接受阈值的虚线水平线）用一个代表性的转录因子家族成员标注，并括号中包含一般的转录因子家族名称。

Para_02

1. 为了进一步比较启动子和增强子之间的链不对称特性，我们分析了链不对称分布，定义为一个方向到另一个方向的活性评分的绝对偏差。
2. 与先前的研究27,28一致，我们观察到在所有检查的细胞类型中，启动子相对于潜在的增强子显示出稍微更强的链不对称效应。
3. 这支持了这样的结论：它们可以包含转录因子结合位点（TFBSs），这些位点促进单向转录（或至少比潜在的增强子更单向）。

Para_03

1. 鉴于启动子的方向依赖性略有增强，我们寻求评估由慢病毒MPRA测量的特定方向启动子活性与内源基因的RNA测序之间的关系。
2. 在所有细胞类型对之间，相同方向的MPRA测量显示出比相反方向的测量更高的相关性（图2d和扩展数据图3d）。
3. 此外，当与内源表达水平进行比较时，我们观察到：（1）匹配细胞类型的MPRA测量显示出几乎与不同细胞类型相同的关联；以及（2）对于每种细胞类型，顺向启动子的MPRA测量显示出比反向启动子更高的相关性（图2d和扩展数据图3d）。
4. 为了进一步评估我们的启动子测量是否能解释细胞类型特异性基因表达，我们计算了启动子活性在其跨细胞类型平均活性上的偏差，以及同一细胞类型的内源基因表达水平的相应偏差。
5. 所有细胞类型在这两组测量之间表现出小于或等于0.11的皮尔逊相关性，这支持了我们之前的结论，即增强子、超级增强子和多梳靶向更有可能解释内源转录活性从核心启动子活性的偏差29。
6. 总的来说，这些结果表明核心启动子具有较弱的方向依赖性，并且当单独测试时几乎没有细胞类型特异性。

Para_04

1. 尽管它们对应于启动子的一个非常短的区域，但位于转录起始位点中心的200-bp核心启动子强烈地再现了内源基因的表达水平（皮尔逊相关系数≈0.55；图2e和扩展数据图3a、b）。
2. 由于启动子的开关状（即开/关）状态，表达值落入双峰分布，这略微增加了这些相关性。
3. 去除所有非表达基因导致MPRA测量的启动子活性与内源表达水平之间的相关性降低（皮尔逊相关系数≈0.43；扩展数据图3c-g）。
4. 值得注意的是，在WTC11细胞中，我们发现了一组转录活跃的基因，其启动子在我们的MPRA中是不活跃的（扩展数据图3b、g）。
5. 进一步分析表明，这一观察结果主要可以由WTC11细胞中替代启动子的使用来解释，因为测试的精确启动子的cap分析基因表达（CAGE-seq）信号在三种细胞类型中与MPRA活性一致（皮尔逊相关系数≈0.60；扩展数据图3h-j）。
6. 我们进行了两种互补的分析以进一步了解可能结合到高表达启动子以及促进低表达和潜在增强子的转录因子家族：（1）使用HOCOMOCO v.1230注释的基序进行富集分析（图2f、g，补充方法和补充表5）；以及（2）从头基序发现（扩展数据图4）。
7. 总的来说，这些分析主要确定了富含CpG的基序以及ETV/ETS相关、KLF相关、NFYA/B/C和THAP11转录因子的转录因子结合位点（TFBS）与活跃启动子相关联（图2f和扩展数据图4a）。
8. 这些基序在许多情况下与检测到的高激活潜力与低激活潜力的增强子不同，其中因素如HNF1B、HNF4A、GATA1/2/6和POU5F1-SOX2作为细胞类型特异性激活因子出现，而CTCF作为普遍的抑制因子（图2g和扩展数据图4b）。
9. 尽管我们没有预料到如此短的200-bp启动子片段能反映内源表达水平，但综合来看，这些发现与先前模型一致，显示富含CpG的启动子与增加的基因表达相关；以及位于转录起始位点中心的核心启动子具有弱的细胞类型特异性，信息密度高且强烈预测基因表达水平29。

Sequence-based models predict activity

Para_01

1. 我们接下来着手训练回归模型来预测调控活性。
2. 我们从一个生化模型开始（补充结果），该模型使用了从三种细胞类型中提取的数千个生化特征的汇编（补充表6）。
3. 该模型能够使用我们的数据通过十折交叉验证方法，在所有三种细胞类型中高精度地预测增强子活性（皮尔逊r ≈ 0.72）（补充图5a）。
4. 许多生化特征与元件活性密切相关（补充图5b-e）。
5. 每种细胞类型关联的可变特征数量导致了可能出现偏差。
6. 然而，训练考虑‘通用’特征集的模型，将所有细胞类型的特征合并，仅略微提高了性能（补充图5f）。

Para_02

1. 基于序列的深度学习模型已经展示了相对于生化模型的强大性能，并且已经被用来预测MPRA数据。
2. 我们评估了四个基于序列的模型的性能，这些模型是在我们的MPRA数据上针对三种细胞类型分别训练的：
3. (1) MPRAnn，一个标准的卷积神经网络（CNN）（补充图6）；
4. (2) MPRALegNet，一个基于LegNet架构的CNN，该架构使用类似于EfficientNetV2的卷积块（图3a）；
5. (3) EnformerMPRA，它使用CNN-变压器架构Enformer来生成一组5,313个预测的生化特征，然后使用这些特征对MPRA数据进行套索回归拟合；
6. (4) SeiMPRA，它考虑由Sei预测的21,907个生化特征，拟合了一个类似的套索模型。
7. MPRAnn和MPRALegNet都经过优化过程，以检测可以提高模型性能的超参数和数据增强设置（扩展数据图5a、b，补充方法和补充表1）。
8. 将性能与我们在相同十个保留数据集折上的生化套索回归模型进行比较，我们观察到所有基于序列的模型都优于生化模型，在三种细胞类型的两种中，MPRALegNet表现最佳（图3b和补充表8）。
9. 尽管我们将EnformerMPRA和SeiMPRA包括在内是为了比较目的，但我们警告说，它们的表现可能被夸大了，因为：
10. (1) 它们的特征集比生化模型大八倍以上，并且在额外的细胞类型和生化标记上进行了训练；
11. (2) 因为其几乎在整个基因组上进行了训练，所以有机会观察到测试集中元素相关的生化标记。
12. 此外，由于基于序列的模型可以访问正在测试的精确200-bp序列，而生化信号缺乏这种空间分辨率，因此它们的表现可能优于生化模型。
13. 结合数据的各个折，我们最好的模型MPRALegNet实现了性能（皮尔逊r ≈ 0.83；图3c），这与实验本身的技术噪声相当（即，重复的复制；补充图2c-e）。
14. 然而，下采样分析表明，模型处于仍然可以从额外训练数据中受益的范围内，因为对于每种细胞类型而言，随着训练集大小的增加，性能呈对数线性改善（扩展数据图5c）。
15. ,

Fig. 3: Sequence-based models predict regulatory element activity.

- 图片说明

◉ MPRALegNet 是一种深度卷积神经网络，用于从测试元件的输入序列预测 cCRE 活性。 ◉ 小提琴图显示了基于序列和生物化学模型在十次交叉验证（CV）中的表现，通过单侧配对 t 检验评估相对于另一模型的改进。 ◉ B，生物化学；E，EnformerMPRA；L，MPRALegNet；N，MPRAnn；S，SeiMPRA；NS，不显著。 ◉ 散点图显示了 MPRALegNet 预测与每种细胞类型观察到的元件活性评分之间的关系。 ◉ d，由 TF-MoDISco-lite38 发现的一组富集基序。左侧，跨多种细胞类型检测到的前三种基序。右侧，每种细胞类型检测到的最上层基序。 ◉ e，热图显示了同型转录因子结合位点剂量（n = 1 到 5 个 TFBS）与 K562 细胞中观察到的 MPRALegNet 预测响应之间的关系。 ◉ f，箱线图显示了 STAT1/4/5A/5B 转录因子家族的全剂量依赖分布，以及在乘法或加法模型下预期的效果。每个组中表示的元件数量如图表上方所示。 ◉ g，热图显示了反映超乘法效应（红色）和亚乘法效应（蓝色）的交互项系数，对于在 HepG2 细胞中拥有指定异型转录因子结合位点对的元件。拟合到观测值和预测值的系数分别显示在热图的右上半部分和左下半部分。 ◉ h，HepG2 细胞中，对于仅包含零个或一个对应于 NFYA/C 或 FOXD2 家族的转录因子结合位点，或者同时具有两个家族的一个转录因子结合位点的元件子集，观察到的分布与 MPRALegNet 预测的活性评分之间的关系，调整了由其他转录因子存在引起的潜在混淆效应（方法）。红线是乘法模型下的预期活性评分。箱线中心线和边缘如图 2 所示。须触须显示最极端的数据点，最大范围为四分位距的 1.5 倍。

Para_03

1. 鉴于MPRALegNet表现良好，我们寻求研究它所学到的原则。
2. 我们在全部MPRA数据上进行了体内诱变（ISM），然后使用TF-MoDISco-lite38来识别具有较大预测效应大小的变异中的基序。
3. 这一策略确定了许多预计在所有细胞类型中激活转录的管家因子，包括NRF1、USF1/2、TFEB和TFE3、JUN和与FOS相关的、KLF相关（KLF/SP）、C/EBP相关和ETS相关的转录因子家族；此外，我们还发现了一个REST的基序，这是一种已知的转录抑制子（扩展数据图6）。
4. 值得注意的是，CTCF与转录激活和抑制均有关联，这表明它可能根据序列环境赋予不同的响应。
5. 在所有细胞类型中最常与转录激活相关的前三种转录因子结合位点（TFBS）是KLF相关、ETS相关和CTCF基序；相比之下，最特异性细胞类型的TFBS在HepG2细胞中是HNF4A/G，在K562细胞中是GATA-TAL1二聚体，在WTC11细胞中是POU5F1-SOX2复合元件（图3d）。
6. 总体而言，发现了许多基序，包括一些未知的转录因子结合位点（TFBS），这些基序被经典的基序富集分析忽略（扩展数据图4），支持了TF-MoDISco方法的互补性。
7. 为了验证涉及的转录因子，我们检查了一个MPRA数据集，在该数据集中，通过CRISPR抑制在K562细胞中进行转录因子敲低测试的元件。
8. 我们能够验证GATA1、NRF1、SP1和FOSL1作为调控因子，但STAT1和ATF4的证据不足，可能是由于敲低效率有限或其他转录因子家族成员提供的补偿效应（扩展数据图7a）。

MPRALegNet predicts TFBS combinations

Para_01

1. 为了深入了解MPRALegNet学习到的组合TFBS效应的本质，我们检查了每种细胞类型中检测到的最丰富的十种激活TFBS基序。
2. 首先，我们测试了同型（相同）TFBS拷贝数对报告基因表达的影响。
3. 在每种细胞类型中，MPRALegNet能够准确预测包含一个到五个指定TFBS位点的元件的激活谱（图3e和扩展数据图7b-e）。
4. 在大多数情况下，转录因子相对于TFBS剂量显示出接近于乘性（即对数相加）的模式（扩展数据图7c）。
5. 然而，某些转录因子家族，如STAT（K562）和ETS相关（WTC11），在最高剂量时表现出次乘性模式（图3f和扩展数据图7e），这表明存在饱和表达效应。
6. 对于某些剂量，还观察到了超乘性（即合作）效应，例如STAT转录因子家族从一个位点增加到两个位点时的增加（图3f）。

Para_02

1. 接下来，我们评估了异型（不同）转录因子结合位点（TFBS）对的乘性效应偏差，通过考虑只包含以下两种元素的元件子集来量化：（1）单一结合位点针对两个转录因子中的任何一个；或（2）同时存在的两个TFBS实例。调整了由其他转录因子存在可能引起的混淆效应（方法），我们在HepG2细胞中观察到不同TFBS对的超乘性和次乘性效应，分别由正或负的交互项系数表示（图3g）。
2. 这些术语的幅度在预测和观察之间表现出强烈的关联（r = 0.92），表明MPRALegNet学到了更复杂的TFBS对之间的组合属性（扩展数据图7f）。例如，ATF3/FOS-JUN和FOXD2位点的同时存在导致了最强的超乘性效应（图3h）；相反，HNF4A/G和NFYA/C位点的同时存在导致了次乘性效应（扩展数据图7g）。
3. 在K562和WTC11细胞中也观察到了类似的结果（扩展数据图7h-l）。总体而言，MPRALegNet能够在所有细胞类型中建模TFBS组合之间的非线性相互依赖关系。

Predicting fine-mapping and variant effects

Para_01

1. 我们接下来考察了 MPRALegNet 在遗传精细映射和变异效应预测中的效用。
2. 我们检查了与全基因组关联研究（GWAS）目录中得出的领先单核苷酸多态性（SNP）处于连锁不平衡（LD）的所有单核苷酸多态性（SNP）。
3. 最初，我们将研究中的七个分块疾病位点相交（扩展数据图 1d、2 和补充图 3）。
4. 对于每个处于连锁不平衡的 SNP，我们使用我们的 K562 MPRALegNet 模型预测了参考等位基因与替代等位基因之间的差异。
5. 我们发现了一些实例，在这些实例中，预测的效果量非常大。
6. 例如，该模型预测了 RBM38 周围的获得功能（gain-of-function）和失去功能（LOF）SNP（rs2426715、rs376911010 和 rs737092；扩展数据图 8a），以及 LMO2 内含子中一个活性 lentiMPRA 序列内潜在增强子的 LOF（rs75395676；扩展数据图 8b）。

Para_02

1. 为了进一步评估这种精细映射预测策略，我们将模型预测与两个互补任务进行了基准测试。
2. 首先，我们使用了UDACHA43和ADASTRA数据库44中HepG2和K562细胞的染色质可及性（转座酶可及性染色质测序（ATAC-seq）和DNase测序）以及转录因子结合数据（染色质免疫沉淀测序（ChIP–seq））中的六个等位基因特异性变异（ASV）集合来验证变异效应数据的表现。
3. 这些显著的ASV提供了关于变异效应的信息，包括转录因子结合或染色质可及性的偏向性，这可能是偏向参考或替代等位基因的等位基因不平衡。
4. 对于所有六种测试的ASV来源和细胞类型组合，在排除非显著ASV或模型预测不确定性太高的情况前后，我们观察到了观察到的分数和预测分数之间的显著关联（费舍尔精确检验优势比> 1.5且P < 0.05；图4a和补充表9）。
5. 我们得出结论，MPRALegNet成功识别了匹配细胞类型的等位基因特异性调控SNP效应。

Fig. 4: MPRALegNet variant effect prediction.

- 图片说明

◉ 散点图显示了在K562和HepG2细胞中通过ChIP–seq、ATAC–seq和DNase–seq数据预测的变异效应与观察到的等位基因特异性差异。图中标出了预测和观测一致（C；蓝色）、不一致（D；红色）或未被考虑（灰色）的情况数量，因为等位基因特异性变异的错误发现率大于0.05或模型预测过于不确定（P值大于0.05；补充方法）。还标出了相应的比值比（OR）（使用补充表9中提供的2×2列联表进行Fisher精确检验）。 ◉ 来自PKLR增强子的饱和突变数据45。顶部，将参考序列缩放至所有替代突变中的平均效应大小，并注释了8个显著转录因子结合位点中的6个，这些位点与已知基序匹配54。中部，单个变异体的测量效应大小。底部，MPRALegNet预测及其对应的皮尔逊（r）和斯皮尔曼（ρ）相关值。

Para_03

1. 接下来，我们进一步验证了我们对变异效应预测的准确性，通过为启动子（F9、LDLR 和 PKLR）和一个增强子（SORT1）生成了 ISM 得分，这些启动子和增强子之前我们在 HepG2 细胞（F9、LDLR 和 SORT1）或 K562 细胞（PKLR）中进行了 MPRA 饱和突变45。
2. 将 MPRALegNet 对 PKLR 启动子的预测与 MPRA 数据进行比较显示，大多数相关的转录因子结合位点（GATA3、KLF9、SP5 和 NFIB）可以被检测到，尽管对于 KLF4 和 GATA2 的预测效应大小相对较小（图 4b）。
3. 总体而言，我们观察到模型预测与观察数据之间在 SORT1 上的相关性为 0.49，在 PKLR 上为 0.65，在 LDLR 上为 0.66，在 F9 上为 0.51（扩展数据图 9），证实了尽管 MPRALegNet 是基于 cCRE 活性的训练，但它仍能部分模拟单个遗传变异的调控效应。
4. 这些结果与 Enformer36 的结果相当（SORT1 上为 0.63，PKLR 上为 0.83，LDLR 上为 0.62，F9 上为 0.28）。综合来看，我们的结果表明了我们的模型如何能够用于预测调控变异的效应。

Characterization of cell-specific factors

Para_01

1. 尽管我们的大规模MPRAs集中于每个细胞类型内的元件活性，但它们并未直接评估每个元件在特定细胞类型中的活性。
2. 因此，我们设计了一个lentiMPRA文库，以测试所有三种细胞类型中的常见一组元件。
3. 这个文库由来自三个细胞系的大约19,000个潜在增强子组成，从以前的大规模MPRA实验中均匀采样，以涵盖广泛的活性范围；
4. 表现出高表达变异以及来自我们以前的大规模MPRA实验中各种细胞类型之间广泛平均表达水平的一组启动子的子集；
5. 双核苷酸洗牌对照；
6. 以及使用先前在HepG2细胞中测试过的合成元件或具有证据表明其具有K562特异性活性的天然元件的一组阳性对照和阴性对照（图5a、补充表10和方法部分）。
7. 元件主要以单向进行测试（启动子为正向，潜在增强子为前向）。

Fig. 5: Assessment of cCRE cell-type-specific activity.

- 图片说明

◉ 联合图书馆在所有三个细胞中进行测试的组成。 ◉ 检测缺乏DNase信号的K562细胞中的活性元件。 ◉ 我们评估了子集中低K562 DNase可及性（垂直红线左侧，由红色方括号指示）的元素在其他两种细胞类型中的DNase信号，然后量化了K562细胞中MPRA活性的平滑核密度估计。 ◉ 主成分分析双图显示第二（PC2）和第三（PC3）主成分，随机绘制每个元素类别最多1,000个数据点。 ◉ 还展示了载荷向量（对应于每种细胞类型）以及每个元素类别最高密度区域的椭圆。 ◉ 小提琴图显示每个元素类别的元件特异性评分分布，并附有关于哪些分布显著大于零的中间值的信息（单侧Wilcoxon符号秩检验，*P < 0.05）。 ◉ 内箱线图按照图3f，h中的描述绘制。 ◉ 训练后的MPRAnn、MPRALegNet、SeiMPRA和EnformerMPRA模型在十个留出数据的交叉验证折叠上的表现，相对于基于生化特征训练的套索回归模型的表现，改善程度通过单侧配对t检验进行评估。

Para_02

1. 我们在所有细胞类型中的每个重复观察到了每种元件10到70个中位条形码（补充图7a）。
2. 元件活性评分在重复对之间表现出强烈的协调性（皮尔逊相关系数为0.98（肝癌细胞系HepG2），0.97（白血病细胞系K562）和0.96（成纤维细胞WTC11）；补充表11和补充图7b-d）。
3. 平均三个重复后，我们观察到共同存在于我们的联合和大规模文库中的cCREs（可确认的活跃调节元件）之间的元件活性评分表现出强烈的同意度（皮尔逊相关系数r=0.90（HepG2），r=0.88（K562），r=0.83（WTC11）；补充图7e）。
4. 在每种细胞类型中，我们观察到阳性对照和阴性对照之间的元件活性评分分布存在强烈差异且弱细胞类型特异性（补充图7f）。
5. 尽管启动子和潜在增强子在所有细胞类型中均显示出显著活性，但来自匹配细胞类型的潜在增强子的活性分布大于来自不匹配细胞类型的潜在增强子（补充图7f）。
6. 虽然启动子和潜在增强子在所有细胞类型中都显示出显著的活性，但是来自匹配细胞类型的潜在增强子的活性分布高于那些来自非匹配细胞类型的潜在增强子（补充图7f）。

Para_03

1. 为了进一步检查细胞类型特异性，我们评估了每个元素类别在每对细胞类型中的行为。
2. 启动子在细胞类型对之间的相关性最强（平均皮尔逊 r = 0.82）；相比之下，潜在增强子在比较从其衍生细胞类型的活性评分与不同细胞类型的活性评分时表现出较弱的相关性（皮尔逊 r = 0.32–0.51（HepG2），r = 0.51–0.65（K562）和 r = 0.64–0.65（WTC11）；补充图8）。
3. 接下来，我们在所有三个细胞系中评估了DNase可及性相对于MPRA活性的关系。
4. 我们观察到强烈的DNase可及性信号并不是MPRA活性的前提。
5. 例如，我们发现尽管K562细胞中的DNase信号几乎不存在，但在HepG2和WTC11细胞中具有高DNase信号的序列仍然可以在K562细胞中导致高MPRA活性（图5b和补充图9a）。
6. 进一步证实这一观察结果，我们注意到5-25%缺乏DNase信号的元件比洗牌的阴性对照更活跃（5% FDR），尽管增加的DNase信号显然与更高比例的活跃元件相关（补充图9b）。
7. 综上所述，我们的结果显示启动子的细胞类型特异性较低，而潜在增强子显示出更强的细胞类型特异性，这与其假定的细胞类型特异性功能一致。

Para_04

1. 我们接下来试图探究每个元件在特定细胞类型中的活性。
2. 我们使用我们在三种细胞类型中的元件活性评分矩阵进行了主成分分析，并去除了主要的主成分，该主成分代表了不同细胞类型之间元件活性的‘普遍’信号。
3. 对主成分2和3（PC2和PC3）的分析显示启动子在WTC11细胞中表达略有偏向，而对照组和潜在增强子在它们来源的细胞类型中表现出更强的偏向性（图5c）。
4. 我们计算了一个元件特异性评分，该评分衡量了每个元件在其所有细胞类型中的平均活性之间的偏差。
5. 这些评分重现了不同元件类别活性富集或耗竭的预期模式，其中HepG2和K562对照在各自的细胞类型中显示出强烈的相对活性；潜在增强子在各自的细胞类型中也显示出强烈的相对活性；而启动子和阴性对照相对于其他细胞类型在WTC11细胞中表现出较弱的更强活性（图5d和补充图10）。
6. 阴性对照在WTC11细胞中表现出较强活性的一个可能解释是干细胞倾向于表现出更广泛的开放染色质状态，这使得它们相对于其他细胞类型更容易受到背景转录水平升高的影响。

Para_05

1. 我们接下来评估了生化模型和基于序列的模型在预测我们的lentiMPRA元件特异性评分方面的表现。
2. 与先前的结果一致，MPRALegNet的多任务版本的表现优于生化模型，而EnformerMPRA的表现优于MPRALegNet和生化模型对于三种细胞类型中的每一种（EnformerMPRA的皮尔逊相关系数约为0.81；图5e和补充表12）。
3. 使用TF-MoDISco-lite38，我们识别了MPRALegNet学习到的细胞类型特异性基序，检测到21个基序与HepG2和K562细胞中的细胞类型特异性活性有关，并在WTC11细胞中检测到12个基序（扩展数据图10）。
4. 与每个细胞类型的前三个排名最高的细胞类型特异性结合基序相关的个别转录因子也表现出预期细胞类型中的强细胞类型特异性表达；此外，CTCF在WTC11中显示出弱富集表达（补充图11）。
5. 需要注意的是，这里测试的转录因子并不代表识别相同基序的所有转录因子家族成员的全面集合，而且未测试的其他家族成员可能进一步解释观察到的基序的细胞类型特异性。

Discussion

Para_01

1. 大规模的MPRA数据集可用于其他细胞系6,7,8,9。
2. 然而，它们主要通过附加型STARR-seq进行测试，需要大量的细胞，提供附加型读出，并倾向于使用强启动子以增加检测活性的能力7。
3. 相比之下，我们改进的lentiMPRA提供了具有‘基因组内’读出的大规模功能数据集。
4. 系统性地、无偏地测试数千个cCREs的能力，使我们能够以高置信度为每种细胞类型识别预测性的生化和基于序列的特征。
5. 然而，尽管测试的序列数量很高，但由于我们的选择标准、未使用的额外cCRE注释分析、标记或工具、所用生化分析的技术问题以及其他因素，许多额外的cCREs可能在我们的注释和随后的分析中被遗漏了。

Para_02

1. 我们测试了所有19,104个已知的蛋白质编码基因启动子在HepG2和K562细胞中的两个方向，并在WTC11细胞中测试了其中的7,500个。
2. 除了观察到与先前研究一致的启动子活性链向偏倚27,28之外，我们还广泛地表征了生成这些开关所需的基于序列的信息。
3. 我们发现位于转录起始位点（TSS）中心的200-bp片段含有足够的序列信息来提供这种开关，并且足以驱动表达，方式类似于它们所关联的基因。
4. 对这些活跃的核心启动子进行测序显示，它们富含已知具有普遍功能的CpG丰富的基序（扩展数据图4a）。
5. 它们还包括已知与转录起始复合物相互作用的KLF相关转录因子家族，以及提供普遍启动子活性的其他转录因子48，以及在普遍表达的启动子中富集的ETS相关家族49。
6. 我们还观察到NF-Y家族的富集，该家族已知与CCAAT框和不含TATA盒的真核启动子相互作用50。
7. 值得注意的是，我们的lentiMPRA设计测试了启动子与一个32bp长的最小启动子一起，这可能会影响启动子活性。
8. 然而，这种方法使我们能够测试数十万增强子和数千启动子，并在同一实验中进行比较。
9. 我们的结果与先前报告27,28相似，显示出启动子的方向偏倚和已知提供普遍启动子表达的基序富集，支持了这种想法：将这个32-bp序列添加到我们检测的启动子上可能没有影响我们的发现。

Para_03

1. 遵循先前的工作10,14,51，我们展示了基于序列的模型能够从MPRA中更优越地预测功能序列。
2. MPRALegNet使我们能够区分出许多对这些预测重要的基序，无论是普遍的还是细胞类型特异性的，并且可以模拟它们的组合效应。
3. 值得注意的是，在所有三种细胞类型中，启动子和增强子中最主要的富集转录因子基序之一是条纹KLF相关转录因子。
4. 条纹转录因子被认为可以在启动子和增强子中提供共可及性并增加其他与转录相关的因子的驻留时间52，并且在最近一项大规模的lentiMPRA研究11中也被发现富集于活性调节元件中，这与其通用功能一致。
5. 尽管MPRALegNet是在三种细胞类型上训练的，但在变异效应预测任务中，非细胞类型特定模型和细胞类型特定模型的类似性能36以及在多种细胞类型中观察到相同变异的相似测量效应量45支持其在其他细胞类型中的应用。
6. 虽然MPRALegNet仅在变异效应预测上与Enformer竞争，但其参数复杂度大约减少了200倍，从Enformer的大约2.49亿减少到大约130万，这提供了一种计算效率高且实用的方法，可以快速地在全基因组范围内预测变异效应。
7. 我们在补充讨论中讨论了这项工作的几个局限性。