人工智能 | 一种细胞类型无关的调控预测多模态 transformer | Cell Genomics

Basic Information

• 英文标题：A multi-modal transformer for cell type-agnostic regulatory predictions
• 中文标题：一种细胞类型无关的调控预测多模态 transformer
• 发表日期：12 February 2025
• 文章类型：Short article
• 所属期刊：Cell Genomics
• 文章作者：Nauman Javed | Bradley E. Bernstein
• 文章链接：https://www.sciencedirect.com/science/article/pii/S2666979X25000187

Highlights

Para_01

1. 多模态变压器学习DNA序列和可及性的表示
2. 一种新的预训练任务对掩码可及性进行训练
3. EpiBERT 提名 caQTLs、调控基序和基因表达
4. EpiBERT 可推广到未观察到的细胞状态

Summary

Para_01

1. 基于序列的深度学习模型已成为解读人类基因组顺式调控语法的强大工具，但无法推广到未观察到的细胞环境中。
2. 在这里，我们介绍了EpiBERT，一个多模态变换器，它通过基于掩码可及性的预训练目标，学习基因组序列和特定细胞类型的染色质可及性的一般化表示。
3. 在预训练之后，EpiBERT可以进行微调以预测基因表达，其准确性可与仅基于序列的Enformer模型相媲美，同时还能推广到未观察到的细胞状态。
4. 所学的表示是可解释的，并且对于预测染色质可及性数量性状位点（caQTLs）、调控基序以及增强子-基因联系是有用的。
5. 我们的工作代表了朝着提高基于序列的深度神经网络在调控基因组学中的泛化能力迈出的一步。

Graphical abstract

Keywords

• deep learning; transformer; gene regulation; genomics; chromatin accessibility; sequence code

Introduction

Para_01

1. 据估计，人类基因组中有10%到20%编码了调控信息，这些信息控制着不同细胞类型和条件下的基因表达。
2. 这些调控序列包含在顺式调节元件（CRE）中，包括靠近基因启动子和可能位于目标基因相当远距离之外的增强子。
3. 通过全基因组分析已经鉴定出数百万个潜在的CRE，这些分析映射了松散包装在可及染色质中的区域、与特征性组蛋白修饰相关联的区域或特定细胞类型中转录因子（TFs）结合的区域。
4. 然而，这些元件仍然没有得到充分的表征。
5. 对于绝大多数CRE，我们还不知道支撑其功能的具体DNA序列，也不知道它们所调控的基因。
6. 尽管测量序列变化对CRE功能影响的实验方法正在发展，但它们缺乏足够的通量来进行全面的调控基因组表征。
7. 因此，能够提名功能序列并预测突变影响的计算模型具有相当大的兴趣。

Para_02

1. 深度神经网络代表了一种有前景的方法来建模和理解CRE功能。
2. 这些模型被训练以从原始基因组序列预测目标信号（例如，基因表达、染色质可及性或转录因子结合）针对特定细胞类型。
3. 一旦训练完成，这些模型可以被解释以预测遗传变异的影响，解析基序语法，或模拟增强子-启动子相互作用。
4. 一个突出的例子是Enformer，这是一个基于大型变压器的模型，它预测基因表达、染色质可及性和其他调控信号，从基因组序列到ENCODE和FANTOM目录中的细胞类型。
5. Enformer采用多任务回归框架，在该框架中，单一基因组序列映射到许多不同的输出（即，对应于各种实验和细胞类型的不同的基因组轨迹）。
6. 虽然这种方法非常强大，但仅基于序列的模型如Enformer的一个缺点是它们必须独立地对每种新细胞类型进行训练，或者至少进行微调，并且无法推广到未观察到的细胞环境。

Para_03

1. 多模态模型通过整合基因组序列与辅助细胞类型特异性信息（例如，染色质可及性）来解决这一局限性。这种框架不依赖于细胞类型，并允许在新的细胞环境中进行预测和调控推断，前提是提供了额外的输入信息。
2. 然而，现有的多模态模型通常在预测跨细胞状态的染色质可及性或基因表达时表现不如仅基于序列的网络，它们的预测分辨率较低（例如，二进制峰值预测或千碱基大小的区间），并且可能需要大量的辅助输入（例如，可及性和三维相互作用数据）。此外，那些以固定辅助信号作为输入的模型无法模拟这些输入本身如何受到基因组序列的影响，从而限制了可解释性。

Para_04

1. 为了应对这些局限性，我们开发了一种多模态神经网络，以一种可泛化的模式预测染色质可及性和细胞类型中的基因表达驱动因素。
2. 我们的模型首先通过掩码预训练学习基因组序列和染色质可及性信号的表示。
3. 我们证明这些表示是可解释的，并且可以用来预测影响染色质可及性的序列变异（染色质可及性数量性状位点[caQTLs]）以及跨细胞类型的调控基序，包括那些在训练过程中未观察到的细胞类型。
4. 在预训练之后，我们将模型微调以预测基因表达，并在保留的基因和细胞类型上实现了与仅基于序列的Enformer模型相当的高准确性。
5. 预训练方法提高了跨细胞类型的表达、功能基序和变异的预测，表明所学表示更准确地预测了上下文特定的CREs及其基因靶标。

Results

EpiBERT architecture and training overview

EpiBERT架构和训练概述

Para_01

1. 我们试图创建一个深度学习模型，能够预测染色质可及性和基因表达的序列决定因素，用于训练过程中未遇到的保留细胞类型。
2. 为此，我们开发了一种多模态神经网络（EpiBERT），该网络通过使用受BERT语言模型启发的预训练目标，整合了基因组序列和局部表观遗传状态（通过染色质可及性测量）（图1A）。
3. 预训练之后，EpiBERT可以微调以用于下游任务，例如预测保留细胞类型的基因表达（图1B）。

图片说明

◉ 图1。EpiBERT训练概述（A）描绘EpiBERT掩码可及性预训练的示意图。在预训练过程中，模型观察特定细胞类型每个输入窗口的掩码ATAC-seq、相应的序列和全局基序富集（STAR方法）。然后，模型预测掩码ATAC-seq信号。预训练后，可以解释该模型以预测调控基序和变异效应。（B）预训练的EpiBERT模型可以微调以预测观察到和未观察到的细胞类型的基因表达。微调后，可以解释该模型以预测增强子-基因（E-G）连接。（C）EpiBERT预训练的高级网络图。在每次训练迭代期间，模型输入包括随机选择的位点的一热编码基因组序列（524 kb），随机掩码（约15%）归一化的ATAC-seq（对应于随机选择的细胞类型的相应位点，4 bp分辨率），以及相对基序富集（STAR方法）。输入通过卷积层和前馈网络（FFN）层传递，合并，并由线性（可扩展神经网络内核）注意力层处理。模型预测掩码区域内的归一化ATAC计数（128 bp分辨率）。训练过程进行2500万次独立迭代/示例（每个对应一个特定的位点和细胞类型）。

Para_01

1. EpiBERT的核心架构主要基于Enformer模型，这是一个处理基因组序列的混合卷积神经网络和变换器（图1C和S1A）。
2. 然而，EpiBERT用线性缩放注意力层替换了普通的注意力层，这使得在最小的性能成本下可以处理非常大的序列窗口，适用于基因组任务（图S1B-S1D；STAR方法）。
3. 在预训练过程中，EpiBERT通过类似于BERT中的掩码语言建模的掩码回归预训练目标来学习序列和可及性的嵌入。
4. 它迭代地对约34,000个位点（大小=524kb）的DNA序列以及训练语料库中每种细胞类型的相应可及性谱进行采样（约740次转座酶可及染色质测序[ATAC测序]谱，从ENCODE、CATLAS和GEO收集；STAR方法）。
5. 为了为每种细胞类型提供跨调节背景，我们还提供了归一化的相对基序富集（图S1E；STAR方法）。
6. 在这个过程中，每个524-kb窗口中有15%的部分被掩码，即将随机选择的1.5-kb区域的可及性信号设置为零。
7. 然后，该模型的任务是在这些区域内重建被掩码的信号。
8. 为了模型评估和测试，我们保留了约4,000个524-kb位点（约占整个基因组的10%），这些位点的选择旨在最大程度地减少与预训练期间观察到的区域之间的序列相似性（STAR方法）。
9. 我们还预留了34种具有不同元数据注释的细胞类型，这些细胞类型不同于训练集中的细胞类型，用于验证和测试（STAR方法）。

EpiBERT accurately predicts masked ATAC signal for hold-out genomic intervals

EpiBERT 准确预测留出基因组区间的掩蔽 ATAC 信号

Para_01

1. 我们训练了模型2500万步，这大约相当于在整个训练集上进行了一次迭代。随后，我们在几个阶段评估了它的性能。
2. 我们评估了它恢复隐藏的可及性信号的准确性，在从训练集中随机选择的17种细胞类型的保留基因组区域中（STAR方法）。
3. 我们特别向模型展示了归一化的基序富集度、每个保留位置的序列以及相应的可及性谱，后者在一个特定的1.5-kb窗口内被屏蔽。
4. 我们发现EpiBERT在屏蔽区域内定性地预测了可及性，并且对不同细胞状态之间的差异敏感。
5. 以PPPR4位点为例，该模型在胎儿肾上腺神经元和心室心肌细胞的共享和细胞类型特异性增强子和启动子处恢复了屏蔽信号（图2A）。
6. 为了更全面的评估，我们在保留区域内的约40,000个窗口中进行了屏蔽。这些窗口被选中以涵盖在17种细胞类型之间具有差异可及性的区域，以及随机选择的对照区域（STAR方法）。
7. 然后，我们将任务交给EpiBERT来预测这些屏蔽区域内的可及性。
8. 在所有屏蔽窗口和细胞类型中，EpiBERT定量地重现了真实的可及性信号（图2B）。
9. 我们还评估了EpiBERT捕捉细胞类型特异性可及性模式的能力。
10. 在这里，我们比较了每个屏蔽区域的预测信号和真实信号的分布。
11. 该模型有效地预测了这种跨细胞类型的变异，在所有测试区域中达到了0.79的中位相关性（图2C）。
12. 当对ATAC谱的基因组序列进行局部洗牌时，预测性能显著下降，这表明该模型学习了精细的序列和可及性特征。
13. 与经过微调的Enformer版本以及使用局部信号平均值的简单基线相比，EpiBERT也表现出更好的性能（图2B、2C、S2A和S2B；STAR方法）。

图片说明

◉ 图2.EpiBERT准确地填补了跨保留区域和细胞类型的掩蔽ATAC信号。(A)对于胎儿肾上腺神经元（紫色）和心室心肌细胞（蓝色）的训练细胞类型，EpiBERT信号在PPP4R4位点的保留区域上的填补。ATAC-seq轨迹是每百万片段（FPM）归一化的。灰色条形表示应用了1.5-kb的掩码。每个条形下的标注显示了每个掩码区域的观测、预测（pred.）和真实信号。◉ 基因轨迹对应于折叠的RefSeq22注释。(B和C)直方图显示了在40,000个保留区域上，针对17个随机选择的训练细胞类型之间的真实和预测信号之间的相关性分布（r）(B)，以及针对每个保留区域的17个训练细胞类型之间的相关性分布(C)。◉ (D)与(A)类似，EpiBERT信号在TMEM260位点的保留区域上的填补，涉及保留细胞类型，胰岛β细胞（紫色）和胎儿光感受器（蓝色）。(E和F)直方图显示了在40,000个保留区域上，针对17个保留细胞类型之间的真实和预测信号之间的相关性分布（r）(E)，以及针对每个保留区域的17个保留细胞类型之间的相关性分布(F)。◉ 图(B)，(C)，(E)和(F)中的相关性基于每个掩码中心512 bp内的总信号。虚线和M(r)对应于中位数相关性。线条颜色对应于使用未扰动输入的EpiBERT（紫色），局部洗牌ATAC-seq（棕色），局部洗牌序列（灰色），被破坏的基序富集（绿色），以及经过微调的Enformer模型（蓝色，EnformerFT）（STAR方法）。

EpiBERT accurately imputes masked ATAC signals for hold-out cell types

EpiBERT准确地填补了留出细胞类型的掩码ATAC信号

Para_01

1. 接下来，我们评估了EpiBERT在17种从未见过的细胞类型中的预测准确性。这些细胞类型是在训练过程中没有见过的。
2. 我们使用预训练网络预测了相同掩码区域的可及性。
3. 我们向网络提供了每个保留的细胞类型的相对基序富集情况以及每个保留的524-kb窗口内的细胞类型特异性可及性，不包括掩码区域。
4. EpiBERT再次定性地预测了特定位置的掩码信号，并且对不同细胞类型之间的差异敏感。
5. 定量评估确认，该模型在所有保留区域和保留细胞类型上预测掩码信号的准确性与训练细胞类型观察到的相似，再次优于更简单的基准（图2D-2F和S2B）。
6. EpiBERT的预测准确性对于不同的序列特征具有鲁棒性，并且在不同类别的基因组元件之间没有显著变化（图S2C和S2D）。

Para_02

1. 我们还评估了EpiBERT的样本外预测性能如何受到保留数据集相似性的影响。我们通过计算每个保留细胞类型的全局基序富集相对于所有训练数据集的最大皮尔逊相关系数来计算其独特性（STAR方法）。我们发现模型性能与每个细胞类型与其训练集最近邻相似性的相关性（r = 0.40，p < 0.05）。尽管如此，EpiBERT仍然实现了准确的可及性预测（r > 0.80），即使对于相对不相似的细胞类型（最近邻相似性 < 0.70；图S2E）。

EpiBERT attends to cis sequence and accessibility across input loci

EpiBERT关注顺式序列和输入位置的可及性

Para_01

1. 为了更好地理解EpiBERT的表现，我们消融了各种输入特征和架构组件，并检查了模型嵌入。
2. 首先，我们训练了模型的较小变体，并将每个输入中的三个输入之一替换为随机同分布噪声。
3. 我们发现，消融未屏蔽的ATAC-seq信号或序列输入会导致验证集上的预训练性能显著下降（图S2F）。
4. 相比之下，消融基因组范围的基序富集仅略微降低了性能，这表明该模型主要依赖于局部序列和可及性上下文，而不是转调控信息。
5. 移除自注意力层中的位置编码也降低了模型性能，这与先前的研究一致，证明了位置嵌入对于基因组建模的重要性。
6. 3我们还训练了模型的变体，旨在处理1Mb而非524kb的输入长度，但这并没有显著提高模型性能（图S2G）。

Para_02

1. 我们接下来检查了EpiBERT用于染色质可及性预测所关注的特征。
2. 对于一系列示例位点和假定增强子，我们引入了一个掩码，并从每个自注意力层中提取了查询矩阵和键矩阵。
3. 然后，我们计算了所有层和头的平均注意力权重（STAR方法）。
4. 对所得权重的可视化显示了一条强信号沿着对角线分布，表明每个输入位置都关注自身以及附近几个千碱基内的相邻bin（图3A）。
5. 此外，我们观察到强烈的非对角垂直线，对应于由大量查询（行）关注的位置（键）。
6. 这些非对角信号在不同细胞类型之间有所不同，并且经常与特定细胞类型的峰重叠。
7. 人为地移动这些高度关注的区域导致模型关注新的位置（图S3A和S3B）。
8. 因此，对于任何给定位点，EpiBERT根据周围染色质可及性的模式关注不同的序列。
9. 我们还测试了准确性是否受输入窗口大小的影响。
10. 我们使用了与之前相同的测试集，但用随机同分布噪声替换了每个524-kb输入窗口的外围部分（STAR方法）。
11. 这一分析显示，当输入上下文长度减少到524 kb以下时，预测准确性降低，这表明EpiBERT关注整个输入窗口中的序列和可及性特征（图3B）。

图片说明

◉ 图3。EpiBERT预测保留细胞类型中的调控基序实例和caQTLs（A）热图显示了胰腺α细胞（左）和β细胞（右）一个524-kb区域所有层/头的平均注意力权重。上方显示了相应细胞类型的基因轨迹和ATAC-seq信号。（B）条形图显示了EpiBERT预测与真实信号之间的相关性（r），如图2E所示，但使用了不同的输入长度（STAR方法）。（C）EpiBERT（红色）和各种方法预测来自27,500个SNP阴性集的571个dsQTLs23的精确-召回曲线。EpiBERT预测使用了GM12878的ATAC，而Enformer预测使用了GM12878的DNAse-seq。（D）实验与EpiBERT预测效果大小之间的散点图以及Spearman相关系数（ρ）在LCLs中。23（E-H）实验与预测的caQTL效果大小之间的散点图，分别为：（E）批量肝组织，（F）神经祖细胞，（G）胰岛，（H）视网膜。25,26,27,28（I和J）PC3细胞中SAMD4A位点增强子处的基础级别输入∗梯度分数（I）和胰腺β细胞中KTN1位点增强子处的基础级别输入∗梯度分数（J）。每个位点下方是TFModisco输出样本，来自每个细胞类型聚合的输入∗梯度分数的12,500个掩蔽增强子。灰色条对应于PC3细胞中CTCF结合峰（黑条，GEO：GSM3498385）上的1.5-kb掩码或β细胞中JUNB结合峰（黑条，GEO：GSM3387453）上的1.5-kb掩码。图（C）-（H）中的CaQTL预测是基于预训练验证集中表现最好的两个模型的平均值。

EpiBERT predicts regulatory motifs and caQTLs in training and held-out cell types

EpiBERT在训练和预留的细胞类型中预测调控基序和caQTLs

Para_01

1. 我们已经证明EpiBERT能够准确预测跨细胞类型的被屏蔽的ATAC-seq信号，接下来我们探讨我们的模型是否能够识别与染色质可及性相关的序列特征。
2. 首先，我们评估了模型预测与实验上绘制的染色质可及性数量性状位点（QTLs）之间的符合程度。
3. 我们从一个由29个研究淋巴母细胞样细胞系（LCLs）的70个供体中映射的DNase高敏感度QTLs（dsQTLs）集合开始，并将其与一个匹配的阴性背景单核苷酸多态性（SNP）集合进行比较。
4. 我们在LCL GM12878的ATAC-seq图谱中，将每个变异位点为中心的1.5-kb窗口进行屏蔽，并计算了参考序列和替代序列之间预测信号差异的效果大小。
5. 我们发现我们的模型能够有效地将dsQTLs与背景集区分开来（平均精度[AP]=0.51，AUROC=0.91；图3C），优于更简单的基线模型（例如gkm-SVM24），并且接近于专门针对相应细胞类型DNase数据训练的显著更大的仅序列Enformer模型的表现（AP=0.61，AUROC=0.95）。
6. 然后，我们探究是否可以通过计算模型预测与实验测量之间的相关性来预测每个dsQTL的效果大小。
7. 对于LCLs，EpiBERT准确地预测了这些效果大小（Spearman秩相关系数ρ=0.79，皮尔逊相关系数r=0.73），与Enformer的表现相当（Spearman秩相关系数ρ=0.80，皮尔逊相关系数r=0.77；图3D）。
8. 在使用来自87个供体肝脏组织中映射的caQTLs进行的类似基准测试中，我们获得了相似的性能（ρ=0.75，r=0.65；图3E），再次接近Enformer的表现（ρ=0.75，r=0.65），以及在神经祖细胞中的表现（图3F）。
9. 因此，在关注已观察到的细胞类型时，我们的方法和Enformer表现出相似的性能。

Para_02

1. 我们接下来询问EpiBERT是否能够在保留的细胞类型中预测caQTL效应大小，这是一个仅使用序列模型如Enformer无法完成的任务。
2. 采用与上述相同的方法，我们计算了在保留的视网膜28和胰腺岛细胞中映射的caQTL的预测效应大小与测量效应大小之间的相关性。
3. 尽管在训练过程中没有观察到这些细胞类型，EpiBERT仍然准确地预测了caQTL效应大小（图3G-3I）。
4. 这些变异周围序列的基序富集分析确认，它们破坏了这些组织中已知转录因子的基序，证实EpiBERT捕获了每个细胞调控语法的相关组成部分（图S3C和3D）。

Para_03

1. 我们还采用了基于梯度的方法来识别EpiBERT预测的与特定位点和全局可及性模式相关的调控基序。
2. 在这里，我们在保留的细胞类型和地区内屏蔽了感兴趣的增强子，并从底层序列和周围可及性环境中预测它们的可及性。
3. 然后，我们计算了输出预测相对于输入序列的梯度，从而根据核苷酸对其局部可及性预测的重要性进行评分。
4. 这种重要性评分突出了与实验验证的TF结合事件一致的基序（图3I和3J）。
5. 例如，EpiBERT将一个CTCF基序与保留的前列腺癌细胞系PC3中的远端假定增强子的可及性联系起来。
6. 它还在保留的胰腺β细胞类型中的假定增强子内突出了JUN/FOS基序。
7. 在这两种情况下，染色质免疫沉淀测序数据证实了相应TFs对掩蔽增强子的结合。
8. 我们还在选定的保留细胞类型中的12,500个掩蔽可及区域上汇总了基础层面的得分。
9. 随后使用基序发现工具TFModisco来识别相应细胞模型中反复出现的高分基序，结果显示在PC3、胰腺β和视网膜细胞中分别发现了26、24和32个显著基序（图S4-S6；错误发现率<0.05）。
10. 发现的顶级基序包括PC3细胞中的JUN/FOS、ETV1和TEAD1；胰腺β细胞中的RFX6、FOXA1和CREB1；以及视网膜细胞中的SP4、MEF2D和RORA，这些对应于这些细胞类型中已知的特征TF。
11. 这些结果表明，EpiBERT学习到了有意义的细胞类型特异性序列特征。

EpiBERT predicts cell type-specific gene expression from sequence and accessibility

EpiBERT 从序列和可及性预测细胞类型特异性基因表达

Para_01

1. 鉴于染色质可及性和转录之间固有的生物学关系，我们认为预训练模型可以微调以预测转录输出。
2. 我们补充了58个ATAC-seq数据集与配对的RNA表达数据，用于相同的细胞类型（50个用于训练，8个用于保留）。
3. 我们特别使用了RAMPAGE31表达数据（RNA注释和基因表达分析中启动子的映射），它明确地映射了每个基因的转录起始位点（TSSs）和定量输出。
4. 我们首先向预训练模型中添加了一组随机初始化的输出层，包括逐点卷积层和密集层。
5. 然后，我们在另外1000万个步骤中微调模型，针对50个训练细胞类型进行，使用与预训练相同的训练和保留位置集合。
6. 在每一步中，模型观察基因组序列、ATAC-seq谱图以及一个随机采样的训练位置和细胞类型的全局基序富集，并被任务预测该位置的归一化RNA计数。
7. 为了防止过拟合并关联序列特征对局部可及性和基因表达的影响，我们间歇性地屏蔽输入的ATAC谱图，并让模型预测被屏蔽的可及性以及RNA计数（图S7A；STAR方法）。

Para_02

1. 我们评估了EpiBERT在预测保留的基因组区域中的RNA表达方面的性能，并将其与在同一RNA数据集上微调的Enformer模型进行了比较。
2. EpiBERT在训练细胞类型中的保留位点处预测了RNA计数分布（平均r = 0.64；图4A），其准确性可与Enformer（平均r = 0.62）相媲美。
3. 接下来，我们重点关注模型在注释的TSS处特异性预测RNA输出的能力，利用两种独立的指标（STAR方法）。
4. 第一个"跨基因"指标衡量模型在一个细胞类型内预测基因表达的准确性。
5. 在这里，EpiBERT实现了高准确性（平均r = 0.88；图4B），几乎与Enformer（平均r = 0.87）相当，并明显优于简单基线的相关性（基因可及性和表达之间的相关性，平均r = 0.42；STAR方法）。
6. 第二个"跨细胞类型"指标衡量模型在不同细胞类型中预测给定基因表达变异的准确性。
7. 与"跨基因"指标不同，后者可能受到组成型表达管家基因的影响，该指标直接评估模型预测细胞类型特异性表达的能力。
8. 在这项任务中，EpiBERT的表现（平均r = 0.72；图4C）超过了简单的基线（平均r = 0.20）和Enformer（平均r = 0.57）。
9. 根据基因的细胞类型特异性分离基因后，我们观察到EpiBERT相对于Enformer在上下文特异性基因上的预测有了最大的改进（图S7B；STAR方法）。
10. 高"跨细胞类型"表现似乎依赖于在掩码预训练期间学习序列-可及性关系，因为省略这一步骤会导致验证集上的表现显著下降（图S7C）。
11. 我们还确认EpiBERT能够准确地建模CAGE数据集（图S7D和S7E）。
12. 综上所述，这些结果表明，通过掩码ATAC预训练和随后针对RNA预测的微调，EpiBERT有效地学到了一个统一的模型，将ATAC-seq数据和基因组序列结合起来，以预测跨细胞类型的转录输出。

Para_02

1. EpiBERT的一个理论上独特的特征是它能够基于序列、可及性和基因表达之间普遍存在的关系，在训练过程中未观察到的持留细胞类型中推断RNA表达。
2. 为此，我们评估了模型在8种未在训练过程中观察到的持留细胞类型上的性能。
3. 当我们提供序列和可及性信息时，我们的模型可以准确地预测这些细胞类型的表达。
4. 具体而言，EpiBERT准确地再现了RNA计数在持留位点的分布（平均r = 0.62），并且在"跨基因"（平均r = 0.85）和"跨细胞类型"（平均r = 0.60）指标上具有高准确性（图4D-4F）。
5. 这种泛化性能部分依赖于全局基序富集，因为去除这种输入会降低"跨细胞类型"表达预测的准确性。
6. 这些结果证实，经过微调的EpiBERT模型能够泛化到新的细胞状态。

Para_03

1. 迄今为止的结果表明，EpiBERT 准确地从序列和可及性数据预测细胞类型特异性基因表达。这促使我们询问我们的模型是否实际上正在学习可访问元件与其靶基因之间的功能联系。
2. 这个问题引发了我们进一步探索的兴趣。

Para_04

1. 我们专注于一组在K562细胞中经过实验验证的增强子-基因连接。
2. 对于每个基因，我们将序列和周围主要转录起始位点（TSS）的可及性谱（524 kb）提供给EpiBERT，并计算了预测输出相对于这些输入谱的梯度。
3. 检查梯度分数，我们发现高模型梯度区域大致对应于经过实验验证的增强子和预测的假定增强子-基因（E-G）连接（图4G；STAR方法）。
4. 具体而言，EpiBERT为每个测试基因确定的相关基因组位置与经过实验验证的增强子-基因连接集很好地对应（图4H）。
5. 虽然梯度分数在增强子优先级划分上存在噪声，但34这项分析使我们能够将EpiBERT与最先进的预测模型进行比较。
6. 在距转录起始位点较近的距离（<12.5 kb）内，我们发现EpiBERT与活动接触33（ABC）和Enformer模型具有竞争力（STAR方法）。
7. 然而，EpiBERT在优先考虑更远距离的增强子（>12.5 kb）方面表现出优于Enformer的性能，尽管两种方法的表现都不如ABC模型，该模型考虑了经过实验测量的E-G接触频率并且不模拟潜在的序列。
8. 这些结果证实，经过微调的EpiBERT模型有效地捕获了经过验证的E-G连接。

Discussion

Para_01

1. 基于序列的神经网络有可能解决基因组学中的一个基本挑战：预测人类基因组中每个可能的序列变异的影响。
2. 基于Transformer的模型，如Enformer，在预测功能基因组分析的输出方面表现出了非凡的性能，并且有可能产生新的调控见解。
3. 然而，这些模型存在局限性，包括无法推广到未观察到的细胞环境以及由于在多任务框架中整合了许多正交数据类型而导致的次优可解释性。
4. 因此，我们开发了EpiBERT，这是一种多模态Transformer，专门针对序列和染色质可及性进行训练，以预测caQTLs、调控基序、基因表达以及增强子-基因联系，适用于已观察和未观察的细胞类型。

Para_02

1. 我们的方法基于生物学假设，即序列定向相互作用驱动染色质变化，而这些变化反过来控制基因转录。
2. 因此，我们的第一步是训练一个模型，该模型能够有效地学习序列和可及性之间的关系，并且这种关系能够在不同细胞类型之间进行泛化。
3. 我们假设特定区域内的序列和可及性之间的局部关系可以从足够宽的基因组窗口周围的顺式模式中以非细胞类型特异性的方式推断出来。
4. 为了实现这一点，我们利用了一个预训练步骤，在这个步骤中，我们的模型学会了使用周围序列和可及性信号来恢复被掩蔽的ATAC-seq信号。
5. 我们利用高效的Transformer块来处理这么大长度的窗口，显著降低了网络训练和推理的内存和计算需求。
6. 对预训练模型的评估表明，它能够泛化到保留的序列和保留的细胞类型。
7. 模型学到了可以轻松解释的表示，这些表示可以用来恢复caQTLs和调控基序，其准确性与明显更大的仅序列模型相当。
8. 此外，我们的模型能够在未观察到的细胞类型中恢复调控序列、基序和caQTLs，这仅靠序列模型是无法实现的。
9. 在预训练之后，我们将模型微调为使用可及性和序列来预测RNA表达。
10. 在这里，我们的模型展示了相对于仅序列模型的准确基因表达预测能力，并且具有向保留细胞类型的泛化能力。
11. 重要的是，与实验数据的基准测试表明，该模型学习了增强子与其目标基因之间的长程相互作用。
12. 我们的结果代表了朝着提高基于序列的神经网络在调控基因组中的泛化能力迈出的一步。

Limitations of the study

研究的局限性

Para_03

1. 我们也承认了模型的局限性，并指出了进一步发展的方向。

Para_04

1. 首先，EpiBERT的样本外泛化性能取决于保留单元类型中的调控基序使用情况与训练集中表示的相似性。
2. 因此，在训练过程中从未遇到过的调控语法（例如，罕见且高度不同的发育状态）的细胞状态可能会削弱EpiBERT的表现。
3. 虽然我们向EpiBERT提供了全局基序富集作为细胞身份的代理，但我们的分析表明模型的预测主要依赖于每个输入窗口内的局部序列和染色质状态（即未屏蔽的ATAC-seq信号）。
4. 通过改进细胞状态的表示（例如，结合表达或TF结合数据，如其他工作所述），可以提高EpiBERT的泛化能力。
5. 这种方法还可以通过允许预测特定反式调节因子如何与局部染色质状态特征（如增强子可及性和活性）相互作用来增强模型的解释性。

Para_05

1. 其次，EpiBERT 仅依赖染色质可及性数据来推断每个输入位点的调控状态。
2. 然而，我们认识到可及性数据无法捕捉诸如抑制性表观遗传状态和高级染色质组织等关键调控特征。
3. 因此，整合其他方法类似的数据模态，如组蛋白标记、转录因子结合、DNA 甲基化或染色质拓扑结构，可能会显著提高 EpiBERT 的性能。
4. 此外，尽管我们通过使用高效的变压器块将上下文长度增加到超过 1 Mb，但没有观察到模型性能的显著提升。
5. 这突显了需要其他类型的数据（例如，3D DNA 接触频率）和/或替代训练技术，以便让模型能够学习长距离的调控相互作用。

Para_06

1. 第三，EpiBERT预测染色质可及性和基因表达的分辨率相对较低（128 bp），并且没有纠正输入ATAC-seq数据中存在的实验偏差。
2. 最近的工作9,36已经证明了碱基分辨率预测和明确的Tn5偏差校正对于建模ATAC-seq数据的好处，这两者可能提高EpiBERT的表现，特别是在caQTL预测方面。
3. 最后，这里引入的掩码可及性预训练方法可能是对多模态单细胞基础模型的一个有用补充，目前这些模型并未明确包含基因组序列。
4. 强制这些模型学习基因表达与局部基因组序列和可及性之间的依赖关系，可能会产生更丰富的细胞表示，从而有可能提供越来越精确的调控见解。

Resource availability

Lead contact

主要联系人

Para_01

1. 进一步的信息、资源和试剂的要求应直接联系首席联系人Bradley E. Bernstein（bradley_bernstein@dfci.harvard.edu），并将得到满足。

Material availability

材料可用性

Para_01

1. 本研究没有产生新的独特试剂。

Data and code accessibility

数据和代码的可访问性

Para_01

1. 所有用于数据处理、模型架构和训练的代码可在https://github.com/naumanjaved/EpiBERT 和https://doi.org/10.5281/zenodo.14498038 获取。
2. 模型权重可在gs://epibert_models 和 https://doi.org/10.5281/zenodo.14504049 获取。
3. 训练、测试和验证数据可在gs://epibert 获取。

Acknowledgments

Para_01

1. 我们要感谢Anthony Philippakis、Neno Shoresh、Neva Durand、Jack Lanchantin和Salva Casani提供的有益讨论。
2. 我们要感谢Google TPU研究云提供了TPU访问和支持。
3. 这项工作得到了Gene Regulation Observatory基金、Broad研究所Variant-to-Function倡议以及National Genome Research Institute Impact for Genomic Variation on Function Consortium（UM1 HG011986）的支持。
4. N.J.得到了Sharf Green癌症研究基金的支持。
5. B.E.B.是Dana-Farber癌症研究所Richard和Nancy Lubin家族的终身教授，并且是美国癌症协会的研究教授。

Author contributions

Para_01

1. N.J., T.W., 和 B.E.B. 构想了这项研究。
2. N.J. 收集并处理了输入数据。
3. N.J., T.W., 和 A.S. 设计并实施了模型，在 A.R., A.D., 和 K.C.（共同提出了用于所提出变压器模型的注意力模块）的帮助下。
4. N.J. 进行了模型性能和下游预测的分析。
5. N.J., T.W., 和 B.E.B. 撰写了手稿，并得到了所有作者的反馈。

Declaration of interests

Para_01

1. B.E.B.声明其对外界利益包括Fulcrum Therapeutics, HiFiBio, Arsenal Biosciences, Cell Signaling Technologies和Chroma Medicine.

STAR★Methods

Key resources table

关键资源表

Method details

方法细节

Model architecture and implementation

模型架构和实现

Para_01

1. EpiBERT 是在 Tensorflow 2.13.0.41 上实现的。
2. 详细的模型架构及所有层参数见图 S1A。
3. 简而言之，模型输入包括基因组序列、被掩蔽的 ATAC 图谱和每种细胞类型的基序富集，这些输入通过独立的编码器处理。
4. 序列编码器接受 524288 个碱基对的一热编码基因组序列作为输入，并提取局部基序模式。
5. 与 Enformer 模型类似，我们使用了 7 个卷积和池化块（一个初始卷积茎，后跟由 6 个块组成的卷积塔），每个块包含批标准化层、GELU 激活、一维卷积层和大小为 2、步长为 2 的池化层。
6. 这个卷积塔的输出是对输入基因组序列在 27 = 128bp 分辨率下的表示，结果是 4096 个 bin。
7. 与 Enformer 不同的是，我们去除了卷积塔中的残差逐点卷积块，并使用了一个更简单的 softmax 池化操作。
8. 对于池化窗口内的每个位置 i，在长度 L 和通道 C 下，softmax 池化操作生成权重 Wi = X · w，其中 w 是通过学习到的权重矩阵得到的。
9. 最后，每个池化窗口的输出 Y 通过公式 Y = ∑(j=1)^L X · Si 计算得出，其中 Si 对应于在池化窗口内计算出的每个 Wi 的 softmax。
10. 我们发现这些改变大大降低了模型复杂度和运行时间，对验证指标的影响最小。
11. 我们还观察到，使用更简单的 softmax 池化层相比最大池化层可以得到更平滑的归因图。
12. 我们使用了类似的架构来处理 ATAC 信号编码器，但卷积茎和卷积塔较浅，只有 2 个块。
13. ATAC 卷积塔的输出，类似于序列编码器，是对输入 ATAC-seq 图谱在 128bp 分辨率下的表示，结果也是 4096 个 bin。
14. 对于基序处理，我们使用了一个更简单的前馈网络（FFN），并将其铺平到与序列和 ATAC-seq 编码器输出相同的长度维度（4096 个 bin）。
15. 然后，这三个嵌入层的输出在通道维度上进行拼接，类似于其他网络，随后通过一个由 8 个基于可扩展神经网络内核（SNNK）工作的新型线性注意力层组成的转换器堆栈。
16. 简而言之，受 SNNK 启发的注意力层通过非负随机特征计算"广义"注意力函数，即给定键矩阵 K、查询矩阵 Q 和值矩阵 V，Q' = RELU(QW) 和 K' = RELU(KW)，A^t(t)(Q,K,V) = Q'(K'V)，其中 W 是一个正交随机矩阵。
17. 通过改变矩阵乘法的顺序，我们可以将时间复杂度从 O(n2) 降低到 O(n)，其中 n 是序列长度。
18. 据我们所知，我们是第一个使用这些新型 SNNK 层来替代转换器中的注意力层的作者。

Para_02

1. 最后，为了注入位置信息，我们使用了旋转向量编码（RoPE），它统一了相对和绝对位置编码。
2. 具体来说，通过旋转经过仿射变换的关键向量或查询向量来实现相对位置嵌入，如果xm是关键/查询向量，m是在给定序列中的向量索引，则x'm = RdΘ,m xm，其中：
3. Rθ,md=(cosmθ1−sinmθ100⋯00sinmθ1cosmθ100⋯0000cosmθ2−sinmθ2⋯0000sinmθ2cosmθ2⋯00⋮⋮⋮⋮⋱⋮⋮0000⋯cosmθd2−sinmθd20000⋯sinmθd2cosmθd2)

Para_03

1. 因此，如果键和查询向同一方向移动相同的量，绝对嵌入会发生变化，然而相对位置不会改变，因为它仅依赖于向量之间的角度。与标准的位置嵌入不同，旋转嵌入应用于每个变换器层。

Para_04

1. 变压器堆栈之后是一个逐点卷积层和最终的密集层。
2. 为了微调，我们添加了一个额外的输出逐点卷积层和最终的密集层来预测RAMPAGE-seq计数。

Data collection and processing

数据收集与处理

Para_01

1. 对于保留的序列，我们使用了与Enformer相同的序列分割方法。简而言之，Enformer的作者将基因组分割成1Mb区域，并构建了一个链接序列对之间存在超过100kb序列比对的图。
2. 然后，这些分区组件被随机分为训练集、测试集和验证集。
3. 我们使用了这些数据分割，并将每个片段的总长度从196k扩展到524k，并从验证/测试集中移除了任何与训练序列重叠的扩展序列。
4. 这导致了34021个训练窗口，2160个验证窗口和1840个测试窗口。

Para_02

1. 我们汇集了一组来自ENCODE、CATLAS和GEO的配对末端ATAC-seq数据集（表S1）。
2. 对于ENCODE，我们直接下载了与hg38对齐的ATAC-seq bams。
3. 对于GEO数据集，我们下载了原始fastqs，并使用ENCODE ATAC-seq管道和默认参数将其对齐到hg38。
4. 然后使用bedtools处理配对末端bams以生成片段文件。
5. 对于CATLAS数据集，我们直接从CATLAS数据库下载了片段文件，这些文件对应于注释的scATAC-seq集群内所有细胞的伪批量数据。
6. 我们将数据集分为训练集（约95％）、验证集（2.5％）和测试集（2.5％），分别对应738、17和17个数据集。
7. 我们根据正交验证数据的可用性选择了保留的数据集。
8. 我们确定了在视网膜细胞和胰岛细胞中识别的caQTLs集合。
9. 然后我们选择了所有注释的视网膜和胰腺数据类型，并将这些数据作为保留数据集。
10. 接着我们添加了随机选择的细胞类型，直到达到34个保留的细胞类型，这大约占总数据集数量的5％。

Para_03

1. 每个数据集使用MACS240调用了峰，并使用基于读取深度的FDR阈值进行了过滤，该阈值依赖于Zhang和Hocker等人2021年的研究。
2. 为了识别可重复的峰，我们为每个数据集随机抽取了50%的片段文件生成了两个伪复制。
3. 然后，我们只保留了那些在两个伪复制中都被调用且至少有50%重叠的峰。

Para_04

1. 为了生成信号文件，我们从每个数据集的每个片段末端提取了Tn5切割位点（片段5'端后4bp和3'端前5bp）。
2. 然后，信号bedgraph文件通过将每个bedGraph条目除以缩放因子进行片段每百万（FPM）归一化，该缩放因子计算公式为：缩放因子=（总片段数）/1,000,000。
3. 接着，对每个输入窗口内的值进行4bp箱的求和处理。
4. 我们在将ATAC-seq谱图输入模型之前，通过对每个值进行软剪切处理来转换它们，转换公式为：x'=min(150,max(0,x))+sqrt(RELU(x-150))。
5. 输出目标向量是通过在128bp分辨率下从信号bedgraph文件中求和值生成的，并且通过对得到的向量进行软剪切处理来进一步处理，处理公式为：x'=min(2000,max(0,x))+sqrt(RELU(x-2000))。

Para_05

1. 为了生成基序富集，我们通过取所有预训练数据集中ATAC-seq峰的并集，并随机选择在至少50%的所有数据集中存在的50000个峰（用于定义存在于大多数细胞类型中的峰）来生成一种与细胞类型无关的背景峰集。
2. 然后，我们为每个数据集选择了50000个随机峰，并使用MEME套件中的SEA44工具和Vierstra等人定义的693个共识基序模型（2020年发表），计算它们相对于与细胞类型无关的背景的基序富集。
3. 对于每种细胞类型，我们提取了每个共识基序模型的富集q值，得到了一个大小为693×1的基序富集向量。
4. 然后，在将该向量输入模型之前，我们对其进行了最小-最大缩放。
5. 为了计算图S1E所示的距离矩阵，我们计算了每个数据集的全局基序富集向量之间的皮尔逊相关系数。

Para_06

1. 对于图 S2E，我们将最近邻相似性度量定义为每个保留的数据集与训练细胞类型内的数据集之间的最大皮尔逊相关系数。
2. 在计算最近邻相似性和测试集性能之间的相关性时，我们仅限于正常组织（因为这些组织不会有我们的模型无法解释的基因组重排和染色体异常，例如K562，它是三倍体）。
3. 报告的p值是使用scipy.stats.pearsonr函数计算的。

Para_07

1. 为了对RAMPAGE-seq进行微调，我们从ENCODE门户下载了与hg38对齐的预处理RAMPAGE-seq bams。
2. 然后我们将每个bam文件配对到同一组织/供体或细胞系对应的ATAC-seq数据集，并使用Picard-Tools中的CrosscheckFingerprints46工具验证共享基因型（表S2）。
3. 对于少量的RAMPAGE-seq实验，我们在ENCODE门户中找不到精确的供体匹配，因此我们选择了相同组织类型的随机ATAC-seq实验作为相应的可及性谱。

Model pre-training and evaluation

模型预训练和评估

Para_01

1. 在每个模型迭代期间的预训练过程中，EpiBERT观察一个524kb窗口内的序列，从随机选择的训练细胞类型中获取该窗口对应的掩蔽ATAC-seq谱图，以及同一细胞类型的全局基序富集。
2. 输入的每个ATAC-seq谱图中有15%被掩蔽，通过将信号设置为0来实现，覆盖随机的1536bp区域。
3. 这些区域的选择是为了包含至少一个在524kb窗口内存在的已调用ATAC-seq峰，并且不会延伸过窗口的边缘。
4. 受间隔掩蔽方法的启发，我们选择了1536bp作为掩蔽大小，迫使网络依赖远距离相互作用来推断掩蔽信号，而不是依赖局部信号。
5. 网络的任务是在每个掩蔽的bin内预测标准化的ATAC-seq计数。
6. 我们使用了与Enformer相同的泊松负对数似然损失，但仅计算掩蔽区域的损失。
7. 简而言之，对于每个具有预测标准化ATAC-seq计数ypred和真实计数ytrue的掩蔽bin i，损失L定义为L = ∑i (ytrue − ypred) log(ypred)。
8. 预训练进行了380,000步，在TPU v3-64上批量大小为64，总共训练了约2500万个样本，这大约对应于预训练语料库的总大小。
9. 这相当于大约4.5天的训练时间。
10. 我们使用了Tensorflow 2.13.0版本的AdamW优化器实现，初始学习率在前5000步内逐渐上升到1.0e−04，随后通过余弦衰减至1.0e−05，权重衰减为1.0e−05，梯度范数裁剪为1.0，并使用了其他参数的默认值。
11. 训练数据通过随机移动输入序列（最多±2bp），输入序列的反向互补以及以50%的概率反转输入ATAC谱图来增强，同时在输入ATAC谱图上添加均值和标准差为1.0e−04的正高斯噪声。

Para_02

1. 为了跟踪模型在训练过程中的性能，我们在验证间隔内对所有验证细胞类型中被屏蔽区域的模型预测进行了皮尔逊相关性计算，每8500步（约500k样本）进行一次。
2. 我们在验证集上使用的屏蔽率、屏蔽大小和数据增强与训练期间相同。
3. 我们确保每次模型验证时使用相同的随机种子，以确保掩码和增强不会在不同轮次之间发生变化。
4. 我们选择了在验证集上皮尔逊相关性最高的两个模型作为最终模型，用于预训练模型评估（包括caQTL和motif分析），并将皮尔逊相关性最高的模型用于下游微调。

Para_03

1. 我们也尝试了使用更长的输入序列长度来训练模型。
2. 我们在多达128块TPU v5e芯片上进行了各种长度的预训练，最长达到1百万个碱基对。
3. 在这个长度上，验证指标没有观察到改进，因此研究中大多数情况选择了51.2万个碱基对的长度作为成本和性能之间的实际权衡。
4. 在1百万个碱基对的长度上，所有验证指标都没有观察到损失，因此这种方法应该适用于这个长度及以上的应用，并且内存和计算成本呈线性增加。

Para_04

1. 为了评估预训练模型，我们收集了17个随机选择的训练细胞类型和34个与测试区间重叠的保留细胞类型的ATAC-seq峰的并集。
2. 由于计算限制，我们仅选择了17个训练细胞类型而不是超过700个数据集的全部。
3. 然后，我们随机选择了其中的20,000个区域，并提取了中心1536bp作为测试的正样本。
4. 我们还包含了20,000个未被任何17个训练或34个保留细胞类型中的ATAC-seq峰识别的1536bp随机背景区域作为负样本。
5. 总共，这产生了40,000个测试区间。
6. 然后，我们将每个1536bp窗口对称扩展以创建524kb的基因组窗口。
7. 对于用于评估的每个测试基因组窗口和细胞类型，我们向预训练模型提供了基因组序列、覆盖中心1536bp的ATAC-seq谱图以及该评估细胞类型的全局基序富集。
8. 这个过程类似于模型训练/验证，但不使用任何随机序列移动、反向互补、dropout或高斯噪声添加。
9. 然后，我们提取了基因组窗口中心512bp内的预测信号作为对该区间的模型预测。
10. 为了计算在打乱序列输入上的性能，我们向训练好的模型提供了未损坏的ATAC-seq谱图和基序富集，但在相邻的16bp窗口内打乱了序列输入以保留粗粒度的序列特征（例如GC含量），但破坏了序列基序。
11. 对于打乱的ATAC-seq输入，训练好的模型提供了未损坏的序列和基序富集，但ATAC-seq输入是在相邻的512bp窗口内打乱的。
12. 对于打乱的基序输入，模型提供了未损坏的序列和ATAC-seq信号，并且基序富集在特征上进行了随机打乱。

Ablations

消融研究

Para_01

1. 由于计算限制，我们使用了一个较小版本的EpiBERT进行模型消融研究，称为EpiBERT-small，其变压器层的数量和每个卷积层的通道数都是原来的一半。
2. 对于序列消融，输入基因组序列被随机打乱。
3. 对于ATAC-seq消融，输入信号被替换为来自0到150之间均匀分布的随机噪声，这是允许的最大ATAC信号。
4. 对于基序富集消融，输入被替换为零。
5. 对于旋转位置编码（RoPE）消融，RoPE嵌入被简单移除，没有用于修改每个SNNK-注意层中的查询和关键张量。
6. 每个模型变体训练了1000万步。

Fine-tuning

微调

Para_01

1. 为了模型微调，我们在预训练模型中添加了一个额外的逐点卷积层和输出密集层来预测标准化的RAMPAGE-seq计数。
2. 我们利用了在预训练期间使用的相同训练、验证和测试区间。
3. 在微调过程中的每一步，模型观察一个随机选择的524kb窗口的基因组序列以及随机细胞类型和相应全局基序富集的掩蔽ATAC-seq谱。
4. 我们使用了与预训练期间相同的掩码程序，但将掩码大小减少到1024bp，掩码率降低到7.5%。
5. 我们采用了与预训练期间相同的数据增强方法，包括向输入的ATAC谱添加随机高斯噪声、+/− 2bp的随机移位以及随机反向互补。
6. 微调模型的输出包括对输入窗口内标准化ATAC-seq和RAMPAGE-seq计数的预测。
7. 我们使用了一种组合损失函数 L = λLATAC + LRNA，其中 LATAC = ∑i (yt rue − ypred) log(y pred) 对于掩码的bins i，LRNA = ∑k (yt rue − ypred) log(y pred) 对于整个区间的所有bins k，λ 是一个经验确定的参数设置为0.001，它降低了ATAC损失项的权重。
8. 我们使用了64的批量大小和两个独立的AdamW优化器进行微调。
9. 第一个优化器使用1.0e−04的学习率更新模型的基础权重，该基础权重包括所有模型层（不包括用于标准化ATAC和RAMPAGE-seq预测的逐点卷积和密集层的输出头）。
10. 第二个优化器使用更高的学习率5.0e−04，并应用于输出的逐点卷积和密集层。
11. 对于两个优化器，我们在前5000步使用线性预热，随后使用alpha因子为0.05的余弦衰减。
12. 我们还使用梯度范数裁剪至0.05，并使用了其他参数的默认值，包括 ϵ = 1.0e−08。

Para_02

1. 为了评估模型在训练过程中的性能，我们计算了三个性能指标。
2. 首先，我们在每个基因的转录起始位点（TSS）生成了524kb的基因组窗口，这些基因位于预先计算的验证区间内，并且不与任何训练区间重叠（结果得到1731个基因/窗口）。
3. 窗口以每个基因最长异构体的TSS为中心。
4. ‘profile level’指标是通过预测分布和真实分布之间在所有验证区间和细胞类型上的归一化RAMPAGE计数的皮尔逊相关系数来计算的。
5. 为了更直接地追踪基因表达预测，我们将每个基因的三个相邻bin中重叠的预测和真实RAMPAGE计数相加。
6. 然后对这些汇总计数进行log1p归一化，并分别针对每种细胞类型标准化结果，使得验证基因的真实值和预测值的集合均值为0，标准差为1。
7. 然后我们计算了‘跨基因’指标，它是通过每种细胞类型上所有验证基因的真实值和预测值之间的皮尔逊相关系数来计算的，然后将这些系数平均。
8. ‘跨细胞类型’指标是通过每种基因上所有细胞类型的真实值和预测值之间的皮尔逊相关系数来计算的，然后将这些系数平均。
9. 在微调过程中，我们每隔4250步（约272k个样本）计算这三种指标，覆盖所有的训练细胞类型。
10. 我们继续微调模型，直到在至少5个周期内‘跨基因’指标和‘profile level’指标没有增加。
11. 我们选择了在大约1000万个样本（约6次遍历整个微调语料库）时达到最高相关性的模型。

Para_03

1. 为了模型评估，我们计算了与上述相同的三个指标，但利用了所有落在测试区间内的基因的TSS中心间隔，并且这些间隔没有任何与训练区间重叠（导致2106个基因/窗口）。我们分别计算了在训练过程中观察到的50种细胞类型和8种保留的细胞类型的这些指标。
2. We separately computed these metrics for the 50 cell types observed during training and the 8 hold out cell types.

Enformer fine-tuning

Enformer微调

• 待补充
• 待补充

Attention matrix calculations

注意力矩阵计算

• 待补充

caQTL predictions

caQTL 预测

Para_01

1. 为了获得LCL GM12878的平均精度曲线和效应大小分析，我们从Lee等人的研究（2015年）下载了处理过的数据集，该数据集包含了Degner等人（2012年）所映射的每个dsQTL的替代和参考等位基因、hg18坐标、FDR值、rsID以及GKM-SVM评分，并且包含了一个由27.5k个阴性SNP组成的背景集合。
2. 我们将rsID输入到ENSEMBL VEP工具中，将其从hg18坐标转换为hg38。
3. 对于Enformer预测，我们在测试集中每个SNP处以196kb的窗口为中心，并生成一个参考序列和带有替代等位基因的突变序列。
4. 然后，我们对模型进行了4次前向传递，每次传递分别针对参考/突变序列及其反向互补序列。
5. 随后，我们计算了Enformer效应大小，作为在中心4个bin（输出bin 446、447、448、449）中预测信号的差异，这对应于GM12878中的DNase（第12个输出头），并在两条链上取平均值。
6. 对于EpiBERT预测，我们在测试集中每个SNP处以524kb的窗口为中心，并生成一个参考序列和带有替代等位基因的突变序列。
7. 对于ATAC-seq输入，我们使用了ENCODE数据库中的处理后的bedGraph文件（数据处理见上文）（数据集编号为ENCFF962FMH）。
8. 我们在每个变异体处以1536bp的掩码为中心。
9. 然后，我们使用了在预训练过程中验证集皮尔逊相关系数最高的两个EpiBERT检查点。
10. 对于每个检查点，我们对模型进行了4次前向传递，每次传递分别针对参考/突变序列及其反向互补序列。
11. 然后，我们计算了EpiBERT效应大小，作为在中心4个bin（索引2044-2047）中预测ATAC-seq信号的差异，并在两条链上和检查点之间取平均值。
12. 平均精度值和曲线是使用sci-kit-learn中的precision_recall_curve和average_precision_score函数计算出来的。
13. 对于效应大小分析，我们下载了Degner等人（2012年）提供的已计算的效应大小。
14. 对于部分dsQTL，我们发现处理表中记录的参考等位基因和替代等位基因被交换了，之后根据VEP工具中的注释参考等位基因进行了纠正。

Para_02

1. 类似的方法被用于计算EpiBERT效应量，该方法适用于批量肝脏caQTLs、批量胰腺岛caQTLs和视网膜组织caQTLs。
2. 对于每个数据集，我们下载了RASQUAL计算的效应量，并根据作者的指示进行了调整（例如，减去0.50）。
3. 我们将分析限制在距相关峰中心128bp以内且FDR小于0.01的caQTLs上。
4. 对于视网膜组织，我们进一步筛选了存在于杆状细胞和其他一种细胞类型中的caQTLs，并使用这两种细胞类型的最大效应量进行后续计算。
5. 对于批量肝脏组织，我们使用ENCODE ATAC-seq数据集ENCFF005CQM生成EpiBERT输入的ATAC谱。
6. 为了生成这个数据集的Enformer评分，我们使用了第448个输出头（对应于53岁女性成人肝脏右叶的DNase）。
7. 对于胰腺岛细胞，我们使用了GEO数据库中的ATAC-seq数据集（访问号GSE16725052）。
8. 对于图S3C和S3D，我们提取了每个caQTL位点中心的20bp区域，并使用MEME套件中的SEA44工具计算与洗牌背景相比匹配HOCOMOCO数据库的富集基序。
9. 对于图S3C和S3D，我们提取了每个caQTL位点中心的20bp区域，并使用MEME套件中的SEA44工具计算与洗牌背景相比匹配HOCOMOCO数据库的富集基序。

Gradient attribution score analysis

梯度归因得分分析

Para_01

1. 为了进行梯度归因评分分析，我们收集了一组包含12.5k个随机选择的ATAC-seq峰，分别来自PC3细胞（ENCODE访问编号#ENCFF024FNF）和胰腺β细胞（CATLAS细胞类型β1）。
2. 我们在每个峰中心放置了一个524kb的窗口，并屏蔽了中心的1536bp。
3. 然后我们计算了模型的前向传递，并计算了输出预测在中心4个bin（索引2044-2047）相对于输入序列的梯度。
4. 接着，我们提取了梯度得分的中心512个元素（对应于中心4个bin）以及相应的512bp的一热编码基因组序列。
5. 最后，对于这12.5k个输入峰中的每一个，长度为512的梯度集和相应的单热编码序列被存储为numpy数组，并输入到TF-MoDISco-Lite。
6. 得到的seqlets与HOCOMOCO v11.53中的基序进行了匹配。

Enhancer analysis

增强子分析

Para_01

1. 处理后的数据，包括坐标、重要性、TSS-增强子距离和预计算的ABC评分，已从Fulco等人2019年的工作中下载。33
2. 对于每个基因，我们在原始ABC手稿中标注的TSS处为中心设置了一个196kb（用于Enformer预测）或524kb窗口（用于EpiBERT预测）。33
3. 对于Enformer，我们随后计算了输入相对于中心12个bin（第441到452个bin）在输出头#4828（K562的CAGE-seq）预测输出处的输入*梯度得分。
4. 我们专注于输入*梯度得分，因为在Enformer论文中这是最佳评分指标。
5. 然后，我们将重叠于指示候选增强子坐标的bin中的输入*梯度得分相加。
6. 对于EpiBERT，我们向经过微调的模型提供了K562对应的标准化ATAC-seq输入（ENCODE访问编号#ENCFF445RUW）。
7. 然后，我们预测了相应的RAMPAGE-seq谱，并计算了输入相对于中心12个bin（第441到452个bin）在预测输出处的输入*梯度得分。
8. 然后，我们将重叠于指示候选增强子坐标的bin中的输入*梯度得分相加。
9. 使用sci-kit-learn中的precision_recall_curve和average_precision_score函数，计算了每个距离的平均精确度值和曲线。

Quantification and statistical analysis

量化和统计分析

Para_01

1. 所有模型代码和软件的使用细节详见上述方法部分。
2. 对于图S2E，我们使用scipy.stats.pearsonr函数计算了n等于26的相关性和相应的p值。
3. 对于图3C和3H，我们使用了sci-kit-learn中的precision_recall_curve和average_precision_score函数。
4. 对于图3D至3H，我们使用了scipy.stats.spearmanr和scipy.stats.pearsonr函数来计算相关性，样本量在每个图中标明。
5. 对于图3I和3J以及图S4至S6中的TFModisco结果，我们使用了n等于10000，并且使用了tfmodisco-lite包的其他默认设置。
6. 对于图S3C和S3D，我们提取了每个caQTL位点周围的中心20bp区域，并使用MEME套件中的SEA44工具（默认设置）来计算与洗牌背景相比匹配HOCOMOCO数据库的富集基序。

Supplemental information