Basic Information
英文标题：Predicting RNA-seq coverage from DNA sequence as a unifying model of gene regulation 中文标题：从 DNA 序列预测 RNA-seq 覆盖度作为基因调控的统一模型 发表日期：08 January 2025 文章类型：Article 所属期刊：Nature Genetics 文章作者：Johannes Linder | David R. Kelley 文章链接：https://www.nature.com/articles/s41588-024-02053-6
Abstract
Para_01
基于基因组学数据训练的基于序列的机器学习模型通过提供描述其对顺式调控代码影响的功能预测，提高了遗传变异的解释。 然而，由于建模挑战，当前工具无法预测RNA-seq表达谱。 在这里，我们介绍了Borzoi，这是一种从DNA序列中学习预测细胞类型特异性和组织特异性RNA-seq覆盖度的模型。 利用Borzoi预测的覆盖率统计信息，我们可以隔离并准确评分DNA变异对多个调节层的影响，包括转录、剪接和多聚腺苷酸化。 在数量性状位点上进行评估时，Borzoi与针对个体调控功能训练的最先进模型相比具有竞争力，且通常表现更好。 通过对派生统计信息应用归因方法，我们提取了驱动正常组织中RNA表达和转录后调控的顺式调节基序。 在不同物种、条件和检测特定调控方面的RNA-seq数据的广泛可用性，突显了这种方法在解析从DNA序列到调控功能的映射中的潜力。
Main
Para_01
长期以来，遗传学的一个目标是准确预测修改人类基因组中的每个三十亿个核苷酸对基因调控活性的影响，包括染色质可及性、转录激活、剪接和多聚腺苷酸化。 这样的预测将极大地提高研究人员解释致病突变和优先考虑全基因组关联研究（GWAS）中涉及位点的功能变异的能力，甚至可以通过功能导向的发现和精细定位来改进GWAS本身。
Para_02
机器学习模型被训练用于从DNA序列预测功能，已经成功地描述了调控语法并解释了遗传变异的影响。 迄今为止，此类模型主要集中在测量活性与局部测序读数成比例的测定上。 例如，转录因子（TF）染色质免疫沉淀与测序（ChIP–seq）或DNase I高敏感位点测序（DNase–seq）以及使用转座酶可及性染色质测序（ATAC–seq）读取表明了在读取对齐的位置发生了TF结合事件或DNA是可访问的。 这使得可以使用相对较短的周围序列区域进行准确预测，通常为500–2,000 bp4,5,6,7,8,9,10。
Para_03
相比之下，最受欢迎的测序分析方法，RNA 测序（RNA-seq），不具有这一特性；RNA-seq 的读段与转录本对齐将取决于包含基因外显子和相关顺式调节元件的更大区域的序列。 与基因 3' 端对齐的读段可能距离影响测定信号强度的启动子和增强子数万至数十万核苷酸远。 此外，RNA-seq 覆盖模式整合了多层基因调控；即转录、剪接、终止或聚腺苷酸化以及 RNA 稳定性。 这些特性使得从序列预测 RNA-seq 覆盖变得具有挑战性。
Para_04
先前的模型仅尝试在总结单一统计量的基因表达后处理RNA测序。 通过处理一个以转录起始位点（TSS）为中心的大区域，几个模型可以预测标准化的基因计数。 这种方法依赖于准确的TSS注释，并忽略了异构体的复杂性。 其他模型预测了基因表达的帽分析（CAGE），该方法测量了帽状RNA（代表TSS）5'端的表达，但不捕捉单个外显子的覆盖情况。 同样，基于序列的转录后调控模型依赖于基因组注释和从RNA测序中提取的转换测量来隔离每个调控机制（例如，拼接的百分比用于剪接）。 然而，此类指标不可避免地难以描述复杂的剪接结果、未注释的新事件或转录、剪接与（内含子）多聚腺苷酸化之间的复杂且有时竞争的关系。
Para_05
建模RNA-seq覆盖度将有几个好处。 首先，RNA-seq远比以前建模的测定方法丰富。 尽管同时建模多个调控层更具挑战性，但它蕴含着巨大的潜力；层与层之间的相互作用很常见，并且同时考虑它们可能会改善对每个调控过程的模型。 迄今为止，模型（例如，那些基于ChIP或ATAC训练的模型）主要集中在单一调控层上。 其次，有大量的RNA-seq数据可用，描述了多种物种中的各种细胞和组织状态。 基于多物种数据训练的模型已被证明可以提高性能，但与RNA-seq相比，染色质谱分析和CAGE基因表达测定在更少的物种上进行过。
Para_06
鉴于哺乳动物基因经常跨越数十万核苷酸，有效的RNA-seq建模需要处理非常大的序列和能够在大距离上传播信息的算法。 最近关于使用自注意力的Enformer模型的工作已经展示了一条实现这一目标的路径。 因此，我们着手通过基于潜在DNA序列的方法来建模RNA-seq以及额外的表观遗传测定的覆盖范围，而不依赖于基因注释的先验知识。 我们开发了一个名为Borzoi的模型，该模型能够有效地学习多层基因调控。 通过将归因方法应用于训练数据中个体RNA-seq实验的预测覆盖模式，Borzoi推导出主要的细胞类型特异性或状态特异性转录因子基序以及全基因组的核苷酸对基因结构和表达的影响图谱。 我们的模型在识别远端增强子和预测遗传变异对基因表达的影响等下游任务上相比Enformer表现更好，并且引入了新的能力来预测变异对剪接和多聚腺苷化的影响，这些能力与最先进的技术相当或超越。 我们预计这个工具包将加速确定许多未解的人类遗传关联影响特征的机制的进程。
Results
RNA-seq model design
RNA测序模型设计
Para_01
RNA测序是对转录和通常被处理的RNA进行碱基分辨率读出。 因此，在碱基分辨率上建模RNA测序覆盖率是理想的。 然而，哺乳动物基因的长跨度意味着我们必须处理非常长的序列来覆盖所有的外显子和相关的调控元件。 计算限制在这些考虑之间创造了一种权衡。 我们倾向于使用较长的序列，尽管会牺牲一些分辨率，选择预测32碱基对箱中覆盖率的524千碱基序列。 训练样本是从跨越人类和小鼠基因组的镶嵌的524千碱基窗口中提取的，因此每个窗口包含位于不同位置的基因。
Para_02
我们的神经网络模型称为Borzoi，在图1a中有所展示。 我们使用核心的Enformer架构，该架构包括一系列卷积和下采样块，随后是一系列在128 bp分辨率下工作的自注意力块27,28。 自注意力是一种关键操作，允许每一对位置交换信息。 从这一点开始，我们利用U-net架构将分辨率提高回32 bp29,30。 对于每次序列长度扩展（以及分辨率增加），我们将来自自注意力块的位置向量上采样，并将其与由初始卷积塔产生的相应大小的功能图结合（见方法部分）。 为了从表示128 bp的嵌入过渡到表示32 bp的嵌入，我们执行这个过程两次，每次上采样的比例为两倍。
Fig. 1: Borzoi: a neural network for predicting RNA-seq coverage from sequence.
Borzoi神经网络架构由若干卷积和下采样层组成，随后是堆叠的自注意力层，这些层使用相对位置编码，在128 bp分辨率下运行，类似于Enformer架构。然后输出被反复上采样并通过额外的匹配U-net连接的卷积层处理，以在32 bp分辨率下进行预测。带有'+'符号的连接表示通过残差卷积将前一层的输出与新层的输入相结合。 通过四个模型复制的预测结果平均得到保留的测试基因INSR（GTEx‘脂肪组织’）的RNA-seq覆盖度预测。 "压缩"尺度是指应用于训练数据的变换尺度（方法）。 当预测CAGE、RNA-seq、DNase-seq或ChIP-seq时，在保留的测试数据各个bin级别上的皮尔逊相关性（n=覆盖轨道的数量）。预测结果是通过对四个模型复制的预测结果取平均值得到。 当比较预测的与测量的外显子范围内RNA覆盖度总和时，基因级别的皮尔逊相关性（n=测序实验的数量）。 在对RNA覆盖轨道进行分位数标准化并减去各轨道间平均基因表达后，基因级别的皮尔逊相关性（n=测序实验的数量）。
Para_03
我们选择使用来自ENCODE的均匀处理的RNA测序数据，提供了866个人类和279个小鼠数据集，这些数据集跨越了多种生物样本，包括细胞系、成人人体组织和发育中的小鼠。 我们还包括了由recount3项目处理的每个基因型-组织表达（GTEx）组织的两到三个重复样本。 为了帮助模型识别重要的调控元件，我们将这些数据与Enformer模型的数千个训练数据集配对，其中包括CAGE、DNase测序、ATAC测序和ChIP测序轨迹（方法部分）。 为了评估模型性能的变化并实现集成，我们训练了四个随机初始化的重复模型。 我们在一组从人类基因组和同源小鼠区域中随机保留的序列上评估了模型的性能。
Borzoi accurately predicts RNA-seq and other assays
Borzoi准确预测了RNA-seq和其他检测方法
Para_01
尽管仅基于底层DNA序列建模RNA-seq覆盖率存在挑战，Borzoi仍能预测外显子-内含子覆盖率模式，即使对于包含许多外显子的长基因也是如此，如图1b中的190kb基因INSR所示。 测试集预测与RNA-seq覆盖率相匹配，使用一个模型副本时，人类样本的平均皮尔逊相关系数为0.74。 当使用整个集合的预测平均值时，皮尔逊相关系数增加到0.75（图1c）。 由于数据处理方法的不同，很难直接将性能与Enformer进行比较（方法部分）。 然而，在重叠数据集上的测试准确性总体相似（扩展数据图1a-e），但有两个例外：DNase的平均皮尔逊相关系数低于Enformer，而CAGE的则高于Enformer。
Para_02
为了研究基因层面的预测，我们将覆盖度聚合并进行对数归一化处理。 当比较预测和测量的基因层面覆盖度值时，我们在保留的基因上观察到平均皮尔逊相关系数为0.87（每个模型复制为0.86）（图1d和补充图1a–d）。 在对实验中的预测进行分位数标准化并减去每个基因的平均表达量（以便该值表示超出平均值的剩余表达量）之后，我们观察到平均皮尔逊相关系数为0.58（每个复制为0.55）（图1e），这表明该模型解释了在轨道之间观察到的大量变异（如组织特异性和细胞类型特异性差异）。 最后，我们注意到Borzoi准确地预测了转录结构内的变化；通过对测试集中前20%的基因进行评估，这些基因在外显子和内含子范围内覆盖度的方差最高，在所有基因和样本中，预测和测量的RNA覆盖度（在bin层面）之间的平均皮尔逊相关系数为0.88（补充图1e）。
Para_03
在补充信息中，我们展示模型依赖于众所周知的调控特征来进行预测，并且模型的注意力矩阵全面捕捉了基因结构（扩展数据图2和补充图2）。 ,
Inference of tissue-specific expression and isoform usage
组织特异性表达和异构体使用推断
Para_01
基因表达是一个由许多调控步骤控制的多方面过程，包括转录起始、剪接和多聚腺苷酸化，这些步骤可能表现出组织特异性效应。 为了研究Borzoi预测组织特异性表达的能力，我们重点关注了五个GTEx组织：全血、肝脏、大脑、肌肉和食道。 我们首先注意到，Borzoi能够准确地预测在预留测试基因上的组织特异性基因表达覆盖度（例如，图2a中展示的与血液相关的基因ADGRE1；另见补充图3a、b）。 我们将一个组织相对于其他四个组织的平均覆盖度的基因水平覆盖度预测值与实际测量值进行了比较，观察到当使用四个模型复制体的集合时，Spearman相关系数范围从0.52到0.75（图2b）。
Fig. 2: Predicting tissue-specific patterns of RNA-seq coverage in normal tissues.
例如，使用Borzoi在五个GTEx组织中对血液特异性基因ADGRE1进行组织特异性基因表达预测。每种RNA测序实验的预测和测量覆盖率在与外显子重叠的bin中汇总（蓝色阴影区域；'max'和'sum'表示最大和总覆盖率）。每个覆盖率轨迹下方显示了外显子注释（GENCODE v.41）。 b, 对于留出的测试基因（n = 1,940），在给定组织中的聚合覆盖率与另外四个组织的平均覆盖率之间的预测和测量倍数变化进行了比较。蓝色和红色点分别代表复制和集成模型的表现。柱状图高度代表平均相关性。插图显示了血液的预测（颜色条表示高斯核密度估计）。 c, 基因SGK1的替代TSS异构体预测示例。TSS使用率通过重叠每个替代起始位点的bin之间的覆盖率比来估算（比率标注在每个轨迹上方）。 d, 在给定组织和其余四个组织的剩余覆盖率比率（COVR）之间，预测和测量的TSS覆盖率比率倍数变化进行了比较（n = 337个留出的至少有两个TSS的基因）。 e, 基因RWDD1的3'UTR APA异构体预测示例。远端位点使用率通过重叠最远端和最近端多聚腺苷酸化位点的bin之间的覆盖率比来估算。 f, 对于给定组织和其余四个组织的APA覆盖率比率之间的预测和测量倍数变化进行了比较（n = 994个留出的至少有两个位点的基因）。
Para_02
基因通常具有替代的TSSs，在不同组织中使用程度不同。 例如，SGK1包含一个上游TSS，在大脑中高度表达而在血液中不表达（图2c；有关更多例子，请参见扩展数据图3a）。 我们计算了从我们的集成预测中得出的最5'端和最3'端TSS的使用比率（方法），并发现与实验测量结果相关（Spearman相关系数R = 0.85；补充图3c），FANTOM5 TSS使用比例（补充图3d）以及组织特异性TSS使用比率的折叠变化（在保留的基因上，Spearman相关系数R = 0.29−0.50；图2d和补充图3e）。
Para_03
3′ 非翻译区（UTR）包含称为多聚腺苷酸化信号（PAS）的调控区域，这些信号可以通过选择性多聚腺苷酸化（APA）生成具有不同 3′ 末端的多个异构体。 例如，RWDD1 在大脑中表现出对最远端的PAS的偏用（Fig. 2e；有关更多例子，请参见扩展数据图3b）。 预测的组织混合的远端至近端多聚腺苷酸化覆盖率比与来自GTEx的测量值高度相关（斯皮尔曼相关系数= 0.81；补充图3f）和PolyADB v.3（参考文献43,44）（补充图3g）。 预测的组织特异性覆盖率比变化显示出与GTEx组织之间测量的变化之间的中度相关性（斯皮尔曼相关系数= 0.23-0.41；图2f和补充图3h）。
Para_04
在补充信息中，我们展示尽管Borzoi能够从匹配的阴性样本中竞争性地识别剪接点，但该模型并未很好地学会预测不同组织间的可变剪接（扩展数据图4；参见讨论）。
Borzoi identifies regulatory motifs driving RNA expression
Borzoi识别驱动RNA表达的调控基序
Para_01
Borzoi通过将归因方法应用于预测的RNA-seq覆盖统计信息，能够直接表征组织特异性顺式调控TF基序。 专注于前一部分分析的五个GTEx组织，我们选择了每个组织中转录每百万（TPM）相对于其他组织最大折叠变化的1,000个基因，并计算了每个基因的组织特异性聚集外显子覆盖梯度。 这些显著性得分描述了每个核苷酸对预测表达的贡献。 例如，对于基因CFHR2，在最大肝特异性显著性的位置处的梯度突出了CEBPA/B和HNF4A/G的基序命中（图3a）。 我们发现，梯度评分在重复实验之间普遍相似，并且与计算机模拟饱和诱变（ISM）紧密匹配（补充图4）。
Fig. 3: Identifying transcriptional cis-regulatory motifs through tissue-specific attribution.
a, 在肝脏特异性基因CFHR2的最大显著性模式下，对五个GTEx组织的梯度归因（平均集合显著性）。所有序列标志图的y轴标尺相同（最小和最大值显示在右上角）。展示了来自HOCOMOCO（v.11）的可能基序命中及其位置权重矩阵（PWMs）。注释的Tomtom E值代表了基序匹配的显著性。插图，比较肝脏和肌肉覆盖轨迹的核苷酸级显著性。 b, MoDISco从对应于四个GTEx组织的梯度显著性中识别出的一组基序簇。显示了MoDISco PWMs、与HOCOMOCO最佳匹配的PWMs以及属于给定簇的seqlets的组织特异性梯度显著性的分布（n=seqlets的数量）。使用双侧Wilcoxon检验计算P值，该检验在具有最大和第二大第95百分位数值的组织之间的梯度显著性之间进行。P值范围从CEBPA/B/D（不显著）的0.075到SPI1/B的5.7×10^-305。E值代表由Tomtom计算的基序匹配的显著性。左下角，比较假定的CEBPA/B/D在全血和肝脏之间的seqlet显著性。每个点的颜色表示目标基因的测量log(TPM)差异。 c, 梯度显著性seqlets属于基序簇的成对GTEx组织之间的平均差异与相应TF基因的测量log(TPM)差异之间的比较。属于同一TF亚家族（HOCOMOCO）的基因的中位TPM被平均。
Para_02
接下来，对于每组1,000个组织特异性基因，我们选择了相应的梯度，并减去了所有其他组织的平均梯度，从而获得了残差组织特异性评分。 我们对所有五个组织基因集运行了TF-MoDISco，这是一种从头开始的基序聚类工具51，并使用Tomtom MEME套件和HOCOMOCO（v.11）52,53将基序聚类与其最有可能的数据库匹配进行对齐。 一些最高得分的基序及其在基因中的显著性分布如图3b所示（另见补充图5a、b）。 我们检测到每个组织已知的调节因子，例如血液中的SPI1/B和IRF4/8，肝脏中的HNF4A/G和HNF1A，大脑中的SOX9和REST以及肌肉中的MYOD1和MEF2D。 组织间共享的基序通常倾向于调节不同的位置（图3b插图）。 类似地，我们同样重现了食道和K562已知的调节基序（补充图5c-e）。
Para_03
最后，我们汇总了匹配每个转录因子的序列单元对之间梯度显著性差异，获得了一个标量分数，该分数描述了某个特定转录因子在一个组织相对于另一个组织的重要性。 这些分数与对应转录因子的观察到的TPM倍数变化高度相关（图3c和补充图5f，g）。例如，当比较血液和肌肉中的转录因子显著性时，Spearman相关系数达到了0.77。 请注意，抑制元件如REST在与大脑进行比较时应该处于非对角位置，因此我们不期望得到完美的相关性。
Improved context use for gene expression prediction
改进的基因表达预测上下文使用
Para_01
我们接下来评估了Borzoi识别和优先处理远端增强子-基因相互作用的能力，这对于细胞和组织特异性调控至关重要。 对于每个目标基因，我们在K562 RNA-seq样本中计算了外显子覆盖预测的输入梯度，突出显示驱动基因表达预测的调控元件。 从梯度显著性得出的统计量，平均到模型集合中，与高通量CRISPR筛选的测量结果进行了比较。 与Enformer相比，Borzoi可以评分距离基因可达两倍远的位点，即262kb，并且我们使用外显子注释而不是转录起始位点注释，后者通常更能抵抗替代异构体的影响。 图4a和图4b展示了HBE1和MYC基因的梯度归属，其中Borzoi正确地识别了近端（距转录起始位点小于20,000 bp）和远端（距转录起始位点大于200,000 bp）增强子，尽管存在假阳性。 Fig. 4a,b 显示了HBE1和MYC基因的梯度归属，在其中Borzoi正确地识别了近端（距转录起始位点小于20,000 bp）和远端（距转录起始位点大于200,000 bp）增强子，尽管假阳性也存在。
Fig. 4: Predicting the impact of context and distal regulatory elements on gene expression.
a, 输入524 kb范围内HBE1的外显子聚集梯度（四次模型复制曲线）。注释了被测量到调控（绿色）或不调控（红色）HBE1的CRE区域。 底部显示了以最远端增强子为中心的192 bp窗口内的输入门控梯度（最大值和最小值显示在右上角）。 b, MYC的外显子聚集梯度。 c, 当使用Borzoi或Enformer梯度在一个每个CRE位置周围局部窗口内计算的统计量来分类是否它调控目标基因时的平均精度（AUPRC）（来自先前出版物65的数据测量）。每个距离箱下的正样本数量和总样本数量如下所示。总样本数量为：<15K) n = 144, (15K–45K) n = 277, (45K–98K) n = 500 和 (98K–262K) n = 1,220。95%置信区间通过1,000倍自助法估计。 d, 使用Borzoi或Enformer梯度来分类Gasperini等人（2019）数据中的调控和非调控CRE时的AUPRCs。总样本数量为：<15K) n = 1,230, (15K–45K) n = 2,445, (45K–98K) n = 4,058 和 (98K–262K) n = 10,051。95%置信区间通过1,000倍自助法估计。 e, 左侧，基于Borzoi DNase覆盖范围预测的与测量的TRIP报告构建表达水平。颜色对应于DamID LMNB1测量。右侧，基于不同评分预测TRIP表达的平均Spearman相关系数（20倍交叉验证（CV））。
Para_02
当比较Borzoi、Enformer和一个距离到TSS基线在分类以前工作中测量的正负增强子-基因相互作用的能力时60,61,62,63,64,65，我们发现Borzoi在所有距离上具有更优的平均精确率（AUPRC）和接收者操作特征曲线下的面积（AUROC）（图4c和扩展数据图5a）。 类似的结果也出现在Gasperini等人（2019）的数据上58（图4d和扩展数据图5b）。 与最近的工作一致66，我们发现在正负样本中，随着TSS距离的增加，平均预测的表达变化百分比普遍呈下降趋势（补充图6a）。 我们在补充图6b-e中更详细地研究了跨转录本的覆盖模式。 通过消融实验，我们发现除了RNA-seq之外，加入DNase-seq和ATAC-seq训练数据可以提高性能（补充图7a-c）。
Para_03
为了进一步证明模型对更广泛的基因组背景的依赖性以进行预测，我们分析了通过TRIP测定法整合到数千个基因组位置的七个不同启动子的表达数据67,68。 我们从多种类的覆盖轨迹中预测了活性评分，包括DNase、组蛋白修饰、CAGE和RNA-seq（补充图8a,b和方法部分）。 总体而言，来自DNase轨迹的评分与测量的表达水平最为一致（图4e和补充图8c；20倍交叉验证，ARHGEF9启动子的斯皮尔曼相关系数R=0.58）。 这些预测比LMNB1 DamID-seq的相关性更高，后者测量核纤层相互作用，并构成一个强有力的基线。
Borzoi prioritizes genetic variants that influence expression
Borzoi优先考虑影响表达的遗传变异。
Para_01
准确预测遗传变异对基因表达的影响对于理解人类群体中遗传关联的调控机制至关重要。 在这里，我们评估了Borzoi区分精细映射的GTEx表达数量性状位点（eQTL）与一组匹配的阴性样本的能力，同时控制了TSS距离。 例如，图5a显示了GTEx全血中SHTN1基因的RNA-seq覆盖预测，包括参考序列和替换单核苷酸多态性（SNP）rs1905542的替代等位基因的改变序列。 我们还展示了在GTEx个体中携带每个等位基因的实际覆盖情况。 由于CEBP结合基序的形成，Borzoi正确地预测了SHTN1表达的上调69,70,71,72（有关更多实例，请参见补充图9a和扩展数据图6a,b）。
Fig. 5: Borzoi predictions of variant effects align with eQTL results and negative selection.
例如eQTL rs1905542。显示的是预测的RNA-seq全血覆盖轨迹，包括参考（蓝色）和替代（红色）等位基因，以及32名纯合子携带参考等位基因和32名杂合子或纯合子携带替代等位基因的全血中测量到的聚集RNA-seq覆盖。 重叠外显子的bin区域用浅蓝色表示。每个等位基因的外显子聚集覆盖度及其比率被注释。 底部显示了具有相同比例y轴的ISM图谱，以及可能的基序命中和Tomtom基序E值。 b，使用Borzoi或Enformer对距离匹配的负样本进行分类时，每个GTEx组织的AUROC。 每个模型的平均AUROC被注释。 c，随着距转录起始位点距离的变化，特定组织的GTEx eQTL分类性能比较。 每个小提琴图展示了中位AUPRC、四分位间距和1.5倍四分位间距作为须线。 P值通过双侧Wilcoxon检验计算（n=49个组织）。 d，左图，使用Borzoi或Enformer预测与观察到的GTEx eQTL效应大小之间的Spearman相关系数。 所使用的差异对数总覆盖率统计量（'SUM'；方法）。每个模型的平均Spearman相关系数被注释。 右图，Borzoi在全血中的预测与观察到的eQTL效应大小对比。 e，左图，通过对常见变异（AF>0.05）从gnomAD中分类单例变异获得的ROC曲线。 右图，平均AUROC及误差线，表示95%置信区间，通过1000次自助法估计（十折交叉验证）。所有变异均采自ENCODE候选顺式调节元件（cCREs）。 AUROC得分显示在图例中。
Para_02
Borzoi预测一个大序列区域的覆盖度，从中必须提炼出变异效应分数。 对于RNA-seq轨迹，我们计算外显子间差异覆盖度的对数倍变化和L2范数（方法部分）。 使用Borzoi的集成模型和L2得分，在区分eQTL方面优于Enformer及其原始总和聚合（平均AUROC = 0.794 vs 0.747；图5b、c）。 即使只使用单一模型或切换回原始总和统计量，Borzoi仍然优于Enformer（AUROC = 0.788；AUROC = 0.772）。 当将效应大小预测与精细映射的eQTL系数进行比较时，Borzoi比Enformer表现出更高的Spearman相关性（平均R = 0.334 vs R = 0.227；图5d和补充图9b、c）。 即使只使用单一模型，Borzoi也优于Enformer（平均R = 0.292）。 在消融实验中，我们发现除了RNA-seq之外，还训练了DNase-seq和ATAC-seq数据以及小鼠数据，显著提高了预测效果（补充图9d）。 我们进一步评估了该模型在eQTL周围基因中优先识别真实eGenes的能力（补充图9e）。 该模型相对于TSS距离基线，表现最好时也只是稍微好一些。
Para_03
为了进一步测试来自Borzoi的变异评分的实用性，我们调查了该模型在GnomAD数据库中区分常见变异（通常是良性的）与匹配的单例集（在单一个体中观察到的罕见变异）的程度，这些罕见变异相对于致病性有所富集。 为了比较，我们考虑了CADD（v.1.6）评分。 限制在ENCODE候选顺式调节元件内，Borzoi和CADD表现出相同的辨别能力（平均AUROC = 0.55；图5e和补充图9f）。 将它们的评分结合起来产生了最高的准确性（平均AUROC = 0.57）。
Para_04
在补充信息中，我们展示Borzoi在预测启动子和增强子中的非编码调控突变方面的表现与Enformer相当，这种比较是通过大规模平行报告分析（MPRAs）来衡量的（补充图10）。
Functional polyadenylation variant interpretation
功能多聚腺苷酸化变异解释
Para_01
我们首先使用归因方法研究了Borzoi预测的3′非翻译区覆盖度，以了解哪些序列特征影响预测形状（图6a）。 已知多聚腺苷酸化调节因子的基序（例如，CFIm、CPSF、CstF）从预测远端多聚腺苷酸化比率的归因分数中浮现出来（图6b）。 尽管我们在3′非翻译区归因中通常没有发现mRNA半衰期的决定因素，但我们确实观察到基因外显子覆盖率的密码子聚集梯度显著性与先前发表的研究中的MPRA测量结果之间存在相关性（皮尔逊R=0.59）（补充图11a）。 我们还注意到，在3′非翻译区中，由于缓冲效应，窗口洗牌的ISM是一种更可靠的归因方法（补充图11b）。
Fig. 6: Predicting APA and 3′ polyadenylation QTLs.
预测和测量的SRSF11（GTEx组织混合样本）远端PAS的RNA-seq覆盖度。图中说明了基于多聚腺苷酸化中心的覆盖比率（COVR）的计算方法。底部显示了基于梯度显著性、ISM和ISM洗牌的归因评分（右下角显示最小值和最大值）。 从混合的GTEx覆盖比率梯度中获得的著名APA调节因子的MoDISco权重矩阵（PWMs），该数据是针对Gasperini基因集58计算得出的。 变异rs114880747的预测RNA-seq覆盖度（GTEx组织混合样本），以及在两个携带参考等位基因的个体和两个杂合子个体中的测量覆盖度（三种组织）。图表中标注了变异和参考COVR统计数据之间的对数值。 基于归因评分（底部；使用相同的y轴比例绘制）表明获得了CstF基序。 没有和带有变异rs80618986的个体的预测和测量覆盖度（每个个体三种组织）。归因评分（底部）表明获得了HNRNPA1基序。 当根据预测的RNA-seq覆盖比率统计量（组织混合样本）分类精细映射的GTEx paQTLs时，AUPRC随最近的3'UTR PAS距离阈值减小而变化。每个点代表一次置换检验（n=100；虚线表示平均值；方法部分）。 比较Borzoi、APARENT2和APARENT2+PolyADB的变异预测的paQTL分类AUPRC。每个点代表一次置换检验（n=100；虚线表示平均值）。 100次置换测试的平均paQTL分类AUPRC，随最近的PAS距离阈值减小而变化。‘A2+S+Borzoi’代表所有模型的集成。
Para_02
我们接下来研究了Borzoi区分精细映射的3′ QTLs（来自eQTL目录79,80；多聚腺苷酸化QTLs (paQTLs)；n = 1,058）和一组经过表达匹配的阴性对照的能力，同时控制了PAS距离。 我们计算变异效应评分作为组织混合GTEx轨迹中任何3′切割位点预测覆盖率比值的最大绝对变化。 我们关注组织混合的预测，因为有限数量的QTLs阻止了特定组织的分析。 图6c-d和补充图11c-d展示了两个paQTLs的覆盖度预测。 与携带替代等位基因的GTEx个体的RNA-seq轨迹相比，Borzoi正确预测了每个变异引起的位点使用的变化。 扩展数据图7a展示了更多例子。
Para_03
从预测的RNA-seq轨迹得出的变异效应评分能够区分与匹配的阴性样本之间的paQTLs，并且在接近最近的多聚腺苷酸化位点（PAS）的距离上准确性呈单调增加（图6e；AUPRC = 0.64-0.74）。 与APARENT2模型22预测的变异得分相比，Borzoi始终表现出更高的准确性（图6f）。 然而，当根据PolyADB中的参考异构体百分比对APARENT2的预测进行缩放时，性能差距减小，这表明上下文是一个重要的决定因素。 我们进一步将APARENT2和Saluki23（方法部分）组成的3'UTR广泛集合进行了比较。 Borzoi在较长距离上表现更好（dAUPRC > 0.050在2,000 bp处），而在接近PAS的位置，其性能与APARENT2和Saluki的表现更为接近（dAUPRC = 0.025在50 bp处）（图6g）。 在较近距离上，所有模型预测（Borzoi、APARENT2和Saluki）的平均排名超过了任一模型的单独表现。
Functional splicing variant interpretation
功能性剪接变异解释
Para_01
我们基于预测的跨越剪接连接点的外显子到内含子覆盖率比，定义了一个以剪接为中心的归因评分（图7a）。 当对Gasperini集合中的基因从组织混合的外显子到内含子覆盖率梯度运行MoDISco时，我们发现了已知的剪接调控基序（图7b）。 使用ISMs解释重复样式的剪接基序时，缓冲效应对结果的影响较小（补充图12a）。
Fig. 7: Classifying sQTLs and intronic paQTLs from RNA-seq coverage predictions.
a, 预测和测量的RNA-seq覆盖度跨越SRSF11基因（GTEx混合组织）的一个外显子。图中说明了计算外显子与内含子覆盖率统计（COVR）的方法。 b, 通过运行MoDISco处理混合GTEx覆盖率比值梯度获得的假定剪接调节因子的PWMs。 c, 对于变异rs55695858预测的RNA-seq覆盖度（GTEx组织睾丸），以及五个具有参考等位基因的个体和五个杂合个体在睾丸中的测量覆盖度（sQTL在睾丸中显著）。图中标注了变异与参考COVR统计数据之间的对数比率。下面显示了归因分数（y轴按相同比例绘制）。 d, 在不同距离阈值下，Borzoi、Pangolin和这两个模型组合在分类精细映射的剪接QTL方面的变异效应预测比较。每个点代表每个模型在给定GTEx组织中的AUPRC指标（中位AUPRC用虚线表示）。 e, 随着距离阈值减小到最近的剪接位点，Pangolin、Borzoi及其组合的平均AUPRC变化情况。
Para_02
我们从eQTL目录中精选了精细映射的剪接QTL（sQTL），并构建了表达匹配和剪接距离匹配的阴性对照（n = 4,105）。 这一相对较大的变异集合允许进行组织特异性分析。 变异效应评分是通过基因跨度内各bin之间的相对覆盖度的最大绝对差异来计算的。 图7c展示了示例sQTL（rs55695858）的RNA-seq覆盖预测（有关更多示例，请参见补充图12b和扩展数据图8a、b）。 此外，还展示了五个GTEx个体在携带或不携带替代等位基因情况下的测量覆盖度。 该变异削弱了一个替代的3'剪接位点，从而上调了相应外显子的延伸。 在比较Borzoi和Pangolin16对因果sQTL与匹配阴性对照进行分类的任务时，Pangolin具有轻微优势（图7d-e和补充图12c；dAUPRC = 0.01，在距已注释剪接位点≤10,000 bp的所有SNP上评估）。 大多数远端SNP是新产生的剪接增益突变，Pangolin可以基于变异等位基因处的局部预测效果轻松对其进行分类，而Borzoi的剪接增益预测似乎不太精确。 相比之下，Borzoi在接近接合点的距离上表现更好（图7d-e；dAUPRC = 0.02，在距已注释接合点≤200 bp的变异上评估）。 重要的是，两种模型的平均排名预测均优于任一模型单独使用（dAUPRC > 0.02）。
Intronic polyadenylation variant interpretation
内含子多聚腺苷酸化变异解释
Para_01
候选多聚腺苷酸化位点经常出现在内含子中，导致多聚腺苷酸化位点与周围的剪接位点之间的竞争。 在这种情况下，内含子要么被剪切掉，要么被保留并进行多聚腺苷酸化。 我们很好奇Borzoi是否已经了解了这种不同调控功能之间的竞争，因此我们从eQTL目录中筛选出更靠近内含子多聚腺苷酸化位点而不是3'非翻译区位点的SNP，并构建了新的匹配内含子多聚腺苷酸化距离的表达控制阴性样本。 Borzoi对精细映射的因果内含子paQTLs预测效果良好，平均AUPRC为0.725（扩展数据图9a、b和补充图13a）。
Discussion
Para_01
本文提出了一种新的基于序列的机器学习模型Borzoi，该模型从大量的RNA测序实验中学习预测测序覆盖度。 Borzoi通过多层调控（包括转录、剪接和多聚腺苷酸化）实现了变异评分和解释，并展示了与最先进的模型在分类精细映射QTL方面的竞争性能。 当对模型复制的预测进行平均时，Borzoi的表现进一步提升。 通过将序列归因方法应用于由预测覆盖度轨迹得出的统计量，Borzoi提供了驱动RNA表达和转录后调控的组织特异性增强子解释。 通过多次消融研究，我们发现，在RNA测序之外还使用DNase测序和ATAC测序数据进行训练，可以持续提高测试集的准确性，并且与eQTL测量和增强子-基因连接数据有更好的一致性。 这一观察表明，最近的多组学数据集（单细胞中同时测量可及性和表达量的数据集）作为联合训练数据将是很有价值的。 变异预测质量仅略微受变异是否出现在训练期间见过的基因组序列中的影响，这意味着遗传学家在使用该模型时可以忽略这一因素。
Para_02
建模 RNA 测序覆盖仍然存在挑战，Borzoi 在预测这些数据方面远非完美。 例如，尽管在不同组织中准确预测了被保留基因的差异 5'（TSS）和 3'（APA）异构体，但大多数组织特异性剪接事件并未被该模型很好地捕获，反而倾向于预测平均 RNA 测序形状。 此外，我们在 Borzoi 的序列归因中未发现与 mRNA 半衰期相关的序列元件。 在存在测序偏差的情况下，解开这些调控层特别困难。 例如，当读段在转录物的 3' 末端以更高的密度对齐时，以及其他干扰因素（如 GC 偏差）导致我们试图根据预测的覆盖率对可变使用的剪接位点进行分类时出现了假阳性。 我们还强调选择适当的归因方法来解释模型的重要性。 虽然输入梯度和 ISM 对剪接接合点和增强子-启动子区域产生了高质量的归因，但我们发现由于缓冲效应，窗口洗牌 ISM 对于 3' UTR 更有效。
Para_03
对于打算在遗传变异分析中使用Borzoi的研究人员，我们建议使用本文中得出的与特定目标基因相关的变异效应评分来优先考虑变异。 这些评分包括（1）目标基因的预测外显子聚集覆盖度对数值变化（用于检测丰度差异），（2）任何3'切割位点之间预测的最大覆盖度对数比差异（用于检测多聚腺苷酸化差异），以及（3）基因体内部任何覆盖区间内预测的最大标准化差异（用于检测剪接差异）。 如果事先不知道目标基因，我们建议使用一种不依赖于基因的统计量，例如基于总L2范数的统计量，来量化整个输出窗口内覆盖模式可能发生的变化。
Para_04
在未来的工作中，我们设想了几个改进的方向。 我们认为，从基于RNA测序的其他测定方法中增加训练数据将进一步提高模型质量；例如，交联和免疫沉淀测序来测量RNA结合蛋白，核糖体图谱来测量翻译，以及测量mRNA半衰期的时间序列。 同样地，我们预计，在那些调节蛋白被扰动的实验上进行训练将改善模型的整体性能，并通过将特定调节因子与序列基序联系起来实现因果推断。 数据量是机器学习成功的关键因素，我们认为，增加更多哺乳动物的RNA测序数据是一个可行的方法，可以增加训练数据和提高模型质量。 相关地，利用个体人类基因组和来自像GTEx这样的群体测序项目的匹配RNA测序数据进行训练，可能有助于进一步提高变异效应预测。 最后，我们渴望整合新的高效注意力模块，以扩大感受野到兆碱基规模，并在更精细的分辨率下进行预测。
Para_05
总的来说，我们开发了一种神经网络模型，用于从序列预测RNA覆盖，并在多个变异解释任务中展示了其性能。 直接对RNA测序进行建模开启了研究广泛实验分析的大门，提高了我们理解遗传变异对基因调控过程影响的能力。
Methods
Para_01
本研究进行的实验不需要特定伦理委员会的批准。
Training data
训练数据
Para_01
训练数据包括大量的人类和小鼠RNA-seq实验。 为了帮助模型使用重要的序列特征来做出RNA覆盖度预测，我们还在训练数据中包含了Enformer和Basenji模型研究的实验测定方法。 这包括FANTOM5联盟提供的经过精心策划的人类和小鼠CAGE测定集合，我们认为这将有助于模型关联TSS使用和强度与多个（替代）TSS之间的RNA-seq覆盖度。 还包括ENCODE和表观基因组路线图的DNase-seq和ChIP–seq数据，以及CATlas的伪批量单细胞ATAC–seq数据，这使模型关注远端调控元件。 相对于之前的分析，我们对数据处理略有不同。 首先，我们在32 bp分辨率下汇总了比对后的读取计数。 其次，我们将CAGE比对的读取按链分开，要求模型预测正向和反向互补链的覆盖度。
Para_02
我们从 ENCODE 收集了 867 个人类和 278 个小鼠的 RNA-seq 覆盖轨迹。 该集合包括来自多种组织和细胞类型的样本，并且测量涵盖了人类和小鼠发育过程中的整个谱系。 可下载的轨迹代表了通过 STAR 对齐程序标准化的独特映射读取的覆盖度102。 大多数实验使用了一种协议来实现分链分析，从而创建了一个正向和反向覆盖轨迹。 我们将 Borzoi 训练为直接预测这些连续覆盖值，预测基于 32 bp 基因组箱。 由于 RNA-seq 具有相对较大的动态范围，我们将每个覆盖轨迹的箱值都通过 3/4 的指数进行归一化。 如果在指数运算后箱值仍然大于 384，则我们对残差值应用一个额外的平方根变换。 这些操作有效地限制了非常高度表达的基因对模型训练损失的影响。 下面的公式总结了应用于目标张量 y 的第 j 个箱中组织 t 的转换：
Para_03
我们将在正文中称这一组转换为‘压缩比例’。 参数的选择使得大多数基因的箱型值小于1000（这是一个相当大的最大值，标准的tensorflow数据类型可以很好地处理它）。 对于大多数下游任务，例如，当由于突变而从预测值计算对数折叠变化时，我们首先通过将逆变换应用于预测值来取消标准化（因此在‘计数’空间中操作）。 一个例外是在可视化测试序列的参考预测时，在这种情况下，除了在384处的残差指数化外，所有变换都被反转，否则接近阈值的小噪声量会被放大。
Para_04
我们通过来自 GTEx 全组织样本的 89 条轨迹补充了训练数据，这些数据由 recount3 项目（GTEx 版本 8 发布）统一处理。 recount3 将 49 个 GTEx 组织聚类为 30 个元组织，结合了高度相关的生理区域（如大脑区域）。 对于每个元组织，我们通过对所有样本的基因表达谱进行 k-means 聚类（k=3）来选择一个子集作为训练数据（尽管一些元组织在 k=2 时就聚类在一起了）。 对于每个聚类，我们选择了与所有聚类成员平均距离最小的样本。 在处理这些数据时，没有考虑 recount3 中的链路信息，这意味着 GTEx 训练轨迹是非链向的，而大多数其他 RNA-seq 轨迹是链向的。 对于这些轨迹，除了上述的指数转换外，我们还通过将对齐片段计数除以其平均长度的倒数来对其进行缩放，以将每个片段视为单一事件。
Para_05
我们将人类（hg38）和小鼠（mm10）染色体打碎，并将这些片段随机分为八个大致相等大小的部分，配对同源区域到同一部分。 其中一个部分被保留用于验证，另一个部分用于测试，其余的数据（约75%）用于训练。 请注意，每当特定的524 kb序列窗口不在训练集中时，所有实验测定（RNA、DNase、CAGE、ATAC、ChIP）的所有覆盖测量都被保留（并且模型不会看到这些）。
Model
模型
Para_01
该模型基于Enformer网络架构，但引入了许多简化和改进以优化RNA-seq预测。 补充图14展示了完整的架构。 Enformer包含两个主要阶段。 首先，重复应用卷积块实现了序列长度的两倍减少，提取局部序列模式，直到序列中的每个位置代表128 bp。 其次，重复应用自注意力（或变压器）块使得每对序列位置之间能够进行长距离交互和信息交换。 Enformer接受一个196 kb的输入序列，并预测在128 bp分辨率下聚合的覆盖率数据。
Para_02
RNA-seq 是对转录 RNA 的基于碱基分辨率的读出。 我们认为增加序列长度和降低预测分辨率对于建模 RNA-seq 都很重要。 哺乳动物基因通常超过全长 >100 kb，如果一个基因的 5' 或 3' 端超出了训练序列窗口（使得其启动子和其他调控信号未被网络的接受域捕捉），可能会阻碍学习。 相反，哺乳动物外显子通常少于 128 bp，以如此粗糙的分辨率对这些外显子周围的覆盖模式进行建模可能会阻碍剪接位点的学习。 然而，计算限制使这些联合目标变得具有挑战性。 因此，我们旨在使用 524 kb 的输入序列，在 32 bp 的分辨率下进行预测。
Para_03
停止普通Enformer架构中的卷积和池化块在32 bp处将意味着自注意力块处理长度为16,384的序列。 这些块需要二次方的记忆复杂度，这超出了当代GPU/TPU硬件的能力，除非进行复杂的优化。 因此，我们选择保持自注意力块的128 bp分辨率。 为了在32 bp分辨率下进行预测，我们转而使用了来自图像分割和对象检测文献中的U-net上采样技术，这些技术解决了类似的问题，即将图像级别的内容确定下来并将其传递回像素分辨率的注释。 简而言之，由自注意力块在128 bp分辨率下预测的输出嵌入通过在每个位置复制嵌入向量的方式被上采样两次。 然后我们应用逐点卷积以匹配与原始卷积塔输出（在自注意力块之前的）在64 bp分辨率下的通道数量。 最后，我们将自注意力块的上采样特征图与卷积塔的中间特征图相加，并应用一个宽度为三的可分离卷积。 这一流程再次使用卷积塔在32 bp分辨率下的中间特征图重复一次。
Para_04
由于这种架构仍然非常计算密集，我们简化了几个Enformer组件。 首先，我们使用最大池化代替注意力池化，这虽然需要一个额外的卷积，但通常只会稍微提升性能。 其次，在初始卷积塔的每个块中，我们只应用一个宽度为五的卷积，放弃了Enformer中使用的带残差连接的第二个卷积。 第三，我们将自注意力块的数量从11减少到8，以减少内存使用。 第四，我们仅使用中心掩码相对位置嵌入，因为额外的距离函数对性能的影响微乎其微。
Training
训练
Para_01
我们在多任务设置中训练模型，以从一个物种预测所有实验的覆盖度，同时附加一个特定于物种的头部到共享的模型主体。 在训练过程中，我们通过动态切换相应物种特定的头部来交替使用人类和小鼠的训练批次。 为了避免在序列边界上预测不准确（由于非对称可见性），我们从每侧裁剪掉一部分，以便将损失计算集中在中间的196,608个碱基对上。 我们使用了泊松损失函数，但将其分解类似于BPnet，以分离幅度和形状项。 每个序列位置具有独立的泊松分布，在数学上等同于表示它们总和的一个泊松分布，然后使用多项分布将计数分配到各个序列位置。 因此，我们在观测和预测覆盖度之和上应用泊松损失，并在序列长度上对标准化的观测和预测覆盖度应用多项分布损失。 这种分解使我们可以更大幅度地加权多项分布形状损失（五倍），我们发现这提高了性能。
Para_02
使用 TensorFlow (版本 2.11)，在这个最大 524 KB 的序列上对该模型进行反向传播会使标准 NVIDIA A100 GPU 的 40 GB 内存达到上限。 每个模型实例使用 Adam 优化器进行训练，批量大小为两个，分布在两个 GPU 上大约 25 天，当验证集准确率趋于平稳时停止训练。
Para_03
我们训练了四个随机权重初始化和序列训练顺序的模型副本。 我们从这四个副本构建了一个集成预测器，其性能通常优于任何一个单独的模型。 请注意，在图1和图2的所有分析中，我们在保留的测试集片段上严格评估模型性能。 在后续分析（例如，图5中的变异效应预测）中，我们不对hg38的训练或测试分割进行区分。 这在技术上意味着集成应用于在训练期间看到的基因组位点。 我们主张这些仍然是无偏分析，因为评估是在未经过训练的域外测量上进行的（例如，精细映射QTL的替代等位基因及其估算效果不是训练数据的一部分）。
Model ablation experiments
模型消融实验
Para_01
Borzoi模型的实例在原始训练数据的较小子集上进行了训练，以评估各种数据模态对最终性能的贡献。 我们改变了模型是否仅使用人类实验数据或同时使用小鼠数据进行训练，是否除核心RNA-seq模态外还使用了其他检测（例如DNase-seq、ATAC-seq、ChIP-seq和CAGE）进行训练，以及模型是否使用U-net组件来提高输出分辨率。 由于组合数量众多，很难获得足够的NVIDIA A100 GPU，在合理的时间内允许它们作为全尺寸的Borzoi模型进行训练。 因此，我们将它们的大小减小到输入长度为393,192 bp（约3000万个可训练参数，每层四个自注意力头），这样我们就可以将它们与批次大小为二一起放在NVIDIA RTX 4090 GPU或NVIDIA TITAN RTX GPU上。 对于每个消融条件，我们训练了两个交叉验证的折叠，从每次折叠的八个基因组hg38或mm10分区中选择不同的保留验证和测试集。 对于基线条件（包含所有特征），我们训练了四个折叠。 训练持续了30至90天，具体取决于条件，并且在验证准确性饱和时停止。
Para_02
以下模型实例被训练了：['多物种'] 训练数据 - 人类（hg38）和小鼠（mm10）的CAGE、DNase-、ATAC-、ChIP- 和 RNA-seq。 架构变化 - 无（基线模型）。 ['多物种（无U-net）'] 训练数据 - 人类和小鼠的CAGE、DNase-、ATAC-、ChIP- 和 RNA-seq。 架构变化 - 移除了U-net。 在128 bp 输出分辨率下进行训练。 ['多物种（D/A/RNA）'] 训练数据 - 人类和小鼠的DNase-、ATAC- 和 RNA-seq。 架构变化 - 无。 ['多物种（RNA）'] 训练数据 - 人类和小鼠的RNA-seq。 架构变化 - 无。 ['人类'] 训练数据 - 人类的CAGE、DNase-、ATAC-、ChIP- 和 RNA-seq。 架构变化 - 无。 ['人类（D/A/RNA）'] 训练数据 - 人类的DNase-、ATAC- 和 RNA-seq。 架构变化 - 无。 ['人类（GTEx RNA）'] 训练数据 - 人类的GTEx RNA-seq。 架构变化 - 无。 ['K562'] 训练数据 - K562细胞中的CAGE、DNase-、ChIP- 和 RNA-seq。 架构变化 - 无。 ['K562（D/A/RNA）'] 训练数据 - K562细胞中的DNase- 和 RNA-seq。 架构变化 - 无。 ['K562（RNA）'] 训练数据 - K562细胞中的RNA-seq。 架构变化 - 无。
Enformer comparison
Enformer比较
Para_01
我们的研究目标是将这个建模框架扩展到新的数据（即RNA测序）上，并且不在重叠轨迹集上超过Enformer的表现，该重叠轨迹集包括CAGE、DNase、ATAC和ChIP测定。 几个建模决策使得Borzoi与Enformer之间的比较不完全准确。 首先，处理更大的序列需要重新处理基因组，因此Borzoi保留的测试集并不完全匹配Enformer的测试集。 其次，我们将数据聚合到了32 bp分辨率，而Enformer使用的是128 bp分辨率，这改变了bin值的分布。 第三，我们将CAGE数据集中的比对读取按链分开。 尽管进行了这些修改，我们还是检查了Borzoi与Enformer（v.3.0）在这组重叠数据集上的测试准确性，发现它们大致相似（扩展数据图1a–d）。
Tissue-specific expression, TSS and APA predictions
组织特异性表达，TSS和APA预测
Para_01
我们评估了从预测的 GTEx RNA-seq 覆盖轨迹中得出的三种不同的统计量，以量化（组织特异性）基因表达、替代 TSS 使用和 APA 异构体丰度（图 2）。 基因表达被量化为重叠外显子 bin 的预测覆盖度之和。 替代 TSS 使用通过计算 GENCODE (v.41) 中每个注释 TSS 下游九个 bin 中的最大覆盖度（因为外显子可能短于九个 bin），并计算每个基因最远端和最近端 TSS 的比值得出。 只有那些与 FANTOM5 中注释的 TSS 相距不超过 50 bp 的 TSS 被包括进来。 APA 位点使用通过计算最远端多聚腺苷酸化位点 (PAS) 上游四个 bin 和最近端 PAS 上游四个 bin 的平均覆盖度比值得出，基于 PolyADB44 中注释的多聚腺苷酸化位点。
Para_02
图2和扩展数据图3中可视化示例的选择方式如下：（1）差异表达示例是从基因中选择的，这些基因在目标组织的外显子聚集覆盖率与另外四个组织的平均覆盖率之间的测量倍数变化最大，基于GTEx RNA测序数据；（2）组织特异性TSS示例是从根据组织匹配的FANTOM5 CAGE数据测量的TSS使用差异最大的基因集中选择的；（3）组织特异性APA示例是从在目标组织中的覆盖率比值相对于另外四个组织的平均覆盖率比值测量倍数变化最大的基因中选择的。 为了降低因GTEx RNA测序数据中的3'偏差导致感知到的APA由基因驱动的风险，我们要求这些基因也在PolyASite 2.0数据库中细胞类型匹配的实验中表现出不同的远端多聚腺苷酸化。 所有示例基因都是从保留的测试集中挑选出来的，并使用四重复组合预测覆盖率。
Input sequence attribution
输入序列归因
错误！！！- 待补充
Para_02
作为预备知识，设 (\mathcal{M"}) 为 Borzoi 模型，x ∈ {0, 1}524,288×4 为输入序列的一热编码，(\boldsymbol{y"}=\mathcal{M"}(\boldsymbol{x"})∈(0,+∞]^{16,384×7,611}) 为（人类）覆盖度预测，(\mathcal{T"}={t_{0},\ldots ,t_{T"}}) 为我们要平均的覆盖轨道 y 中的 T 个索引集合（例如，结合所有与血液相关的轨道），并计算相应的归因分数。 请注意，Borzoi 的原始预测 y 是基于经过各种变换处理的训练数据，这些变换旨在稳定训练（如将数据立方根化，对超出目标值的残差进行额外的指数运算以及重新缩放）。这里，我们假设已经对 y 应用了逆变换，使得张量可以合理地认为反映了计数（还请注意，这些变换是可微的，这意味着梯度显著性可以通过逆操作传播）。
Para_03
以下是用于表达归因、polyadenylation 归因和剪接归因的三个不同汇总统计量的定义：
Log sum of exon coverage (expression attribution)
外显子覆盖的对数和（表达属性）
Para_01
汇总统计量 u 属于 实数集 R 通过聚合基因组中与目标基因外显子重叠的 32 个碱基对 bin 集合 B 等于 b0 到 bB（可选伪计数 C 属于 实数集 R）来计算：
Log ratio of PAS coverage (polyadenylation attribution)
PAS覆盖率的对数比率（多聚腺苷酸化归属）
Para_01
统计量 u 属于 R，通过将覆盖度在 bprox 重叠目标多聚腺苷酸化位点（PAS）之前的五个相邻bin中的覆盖度相加，然后除以匹配的一组在 bdist 之前的bin中的覆盖度（或在 bprox 立即下游，如果感兴趣基因不受多聚腺苷酸化位点后剪接影响）来计算：
Para_02
上述公式假设基因位于正（加）链上。如果基因位于负链上，则必须从bprox + 1到bprox + 5 + 1（以及从bdist + 1到bdist + 5 + 1）求和覆盖。 注意，如果基因位于负链上，覆盖率必须从bprox + 1到bprox + 5 + 1（以及从bdist + 1到bdist + 5 + 1）进行累加。
Log ratio of exon-to-intron coverage (splicing attribution)
外显子与内含子覆盖度的对数比（剪接归因）
错误！！！- 待补充
Para_02
上述定义的摘要统计量与以下归因方法结合使用：
Gradient × input (gradients)
梯度 × 输入 (梯度)
Para_01
给定汇总统计量 u(x)，属性评分 ({\boldsymbol{s"}}\in {{\mathbb{R"}}}^{524,288\times 4}) 通过将输入 x 的梯度与每个核苷酸位置上的平均值相减计算得出：
Para_02
当我们可视化s时，我们只提取位置i处对应于参考核苷酸j的分数（这可以通过与x相乘并在核苷酸上求和来轻松实现）：
ISM
ISM
Para_01
给定一个在x中计算ISM的起始位置pstart和结束位置pend，属性评分(\boldsymbol{s"}\in \mathbb{R"}^{524,288\times 4})的计算方法如下：创建一个新的张量(\widetilde{\boldsymbol{x"}}\in {0,1}^{({p"}{\text{end"}}-{p"}{\text{start"}})\times 4\times 524,288\times 4})，并使每个矩阵(\widetilde{\boldsymbol{x"}}_{u,v"})保存x的一个突变副本，在位置u处将参考核苷酸替换为核苷酸v。然后计算ISM评分为s。 令(\widetilde{\boldsymbol{x"}}_{u,v"})保存x的一个突变副本，在位置u处将参考核苷酸替换为核苷酸v。
Para_02
当我们可视化s时，我们对四个核苷酸的分数进行平均。
Window-shuffled ISM (ISM shuffle)
窗口洗牌ISM（ISM洗牌）
Para_01
给定起始位置和结束位置 pstart 和 pend，窗口大小 M 和重新洗牌次数 N，归因分数 ({\boldsymbol{s"}}\in {{\mathbb{R"}}}^{524,288\times 4}) 如下计算：创建张量 (\widetilde{{\boldsymbol{x"}}}\in {{0,1}}^{({p"}{\rm{end"}}-{p"}{\rm{start"}})\times N\times 524,288\times 4})，包含 (pend − pstart) × N 份输入模式 x 的副本。 对于每个矩阵 ({\widetilde{{\boldsymbol{x"}}}}_{u,v"})（其中 v 表示 N 个独立样本之一），对局部区域 [u − M/2, u + M/2 + 1] 进行双核苷酸洗牌或用均匀随机核苷酸替换该区域中的参考核苷酸。 当计算增强子显著性时进行双核苷酸洗牌（M = 7 和 N = 24，或大窗口尺寸时 N = 8），而用于启动子、剪接位点和 PAS 时使用均匀随机替换（M = 5 和 N = 24，或大窗口尺寸时 N = 8），因为显著特征通常是重复核苷酸的片段。 然后计算归因分数 s 为：
Para_02
当我们可视化s时，我们在N个样本上取平均分：
Tissue-specific motif discovery
组织特异性基序发现
Para_01
我们通过组合（1）选择大量（已测量的）高度组织特异性基因，（2）计算它们的梯度显著性并排除组织共享显著性，以及（3）使用TF-MoDISco（v.0.5.14.1）51和Tomtom MEME套件（v.5.5.2）52对显著性分数进行聚类和注释，可视化了通过GTEx轨道驱动RNA覆盖的学习组织特异性顺式调节基序。 我们首先下载了GTEx（v.8）的测量TPMs（GTEx_Analysis_2017-06-05_v8_RNASeQCv1.1.9_gene_median_tpm.gct.gz）。 我们通过添加一个大致为所有值第一百分位数的小伪TPM来启发式地清理数据（为了避免零），然后根据每个组织略大于第99百分位数的值进行裁剪（为了避免极大数值）。 然后，对于五个预期的GTEx组织（全血、肝脏、大脑-皮质、肌肉-骨骼和食道-肌层），我们计算了相对于其他四个组织平均TPM表达的感兴趣组织的TPM表达基因特异性对数倍变化。 对于每个组织，我们按此指标降序排列TPM矩阵，并选择了前1,000个差异表达最显著的基因，最终得到了总共5,000个基因。
Para_02
我们计算了与每个基因的聚集外显子覆盖率的对数相关的核苷酸级归因分数（输入梯度），涉及5000个基因中的每一个，并且对每个GTEx组织重复了梯度计算。 具体来说，我们将每个GTEx组织与通过recount3获得的两到三个RNA覆盖轨迹相匹配，这些轨迹是我们训练所用的（例如，对于大脑-皮层，我们计算了与三个GTEx大脑元组织轨迹相关的输入梯度显著性）。 针对所有四个模型副本，对前向互补和反向互补输入序列进行了梯度计算，并取平均值。
Para_03
上述计算梯度方法产生了所有5000个基因的五个独立的显著性评分集合（每个组织一套评分）。 接下来，我们通过切片提取出最初被选为在组织x中差异上调表达的1000个基因，并对组织x的残差梯度评分运行TF-MoDISco。 残差评分是通过从组织x的评分中减去其他四个组织评分的平均值来计算的，从而减弱了共享调控基序的显著性并突出组织x特有的基序。 此外，在运行MoDISco之前，我们首先通过计算每个位置上四种核苷酸的标准差来重新加权梯度，然后应用高斯滤波器（标准差=1280；截断=2）处理由此产生的标准差向量，并将梯度评分除以这个平滑后的向量。 这一操作导致下调了具有长连续大数值区域的调控区（通常是启动子区域）的权重，并上调了稀疏调控区（转录增强子）的权重。 为了提高计算效率，我们在调用MoDISco之前提取了中心化的131kb梯度评分（而不是完整的524kb）。 TF-MoDISco使用以下参数执行：'revcomp = true'，'trim_to_window_size = 24'，'initial_flank_to_add = 8'，'sliding_window_size = 18'，'flank_size = 8' 和 'max_seqlets_per_metacluster = 40,000'。 其他参数保持默认值。
Para_04
五种组织特异性MoDISco结果对象被过滤和合并如下：使用Tomtom MEME将每个MoDISco集群的位置权重矩阵与HOCOMOCO（v.11）53个基序进行匹配（每个位置权重矩阵通过信息含量阈值>0.1进行修剪）。只保留E值≤0.1的匹配。选择具有最低P值的匹配作为该集群的代表基序。通过匹配具有相同HOCOMOCO基序的集群并将序列片段坐标合并，将五个MoDISco对象合并，从而为每个假定的基序生成一个序列片段坐标列表。然后通过计算每个序列片段的输入门控梯度的平均重叠坐标来计算每个序列片段的标量组织特异性显著性得分。这些序列片段级别的梯度显著性分布用于评估每个基序的组织特异性。
Para_05
在预测的外显子覆盖度总和应用对数计算梯度之前添加伪计数来复制整个分析，结果几乎相同。 复制分析而不使用残差归因分数运行TF-MoDISco，而是将每个组织特异性覆盖轨迹的原始梯度作为TF-MoDISco的输入，同样产生了可忽略不计的差异。
Tissue-pooled splice motif discovery
组织混合剪接基序发现
Para_01
剪接调控基序是通过计算与剪接属性统计量（外显子对内含子覆盖率的对数比）相关的输入梯度生成的，对于Gasperini数据集中4,778个基因中的每个基因随机选择的一个外显子进行计算。 梯度是根据所有89个Borzoi的GTEx RNA-seq轨迹的平均预测覆盖率来计算的。 梯度在基因之间进行了标准化处理如下：首先我们计算了四个核苷酸的标准差，并且在每个基因的524,288个位置上找到了最大的标准差。 我们将4,778个最大偏差的低端部分在第25百分位处进行了截断（以避免放大幅度非常小的梯度），并将每个基因的梯度除以这个数值。 我们尝试改变百分位阈值（从1到100），结果对此参数具有鲁棒性（使用大致相同数量的支持序列识别出了相同的基序簇）。 最后，为了获得5'剪接基序，我们从每个梯度中提取了一个以剪接供体为中心的192 bp窗口。 为了获得3'剪接基序，我们从剪接受体周围提取了一个192 bp窗口。
Para_02
TF-MoDISco 在 4,778 × 192 × 4 个假设得分上进行了执行，使用了我们通过实验证明更适合退化RNA结合蛋白基序的自定义参数设置：'revcomp = false'，'trim_to_window_size = 8'，'initial_flank_to_add = 2'，'sliding_window_size = 6'，'flank_size = 2'，'max_seqlets_per_metacluster = 40,000'，'kmer_len = 5'，'num_gaps = 2' 和 'num_mismatches = 1'。
Tissue-pooled polyadenylation motif discovery
组织 pooled 的 polyadenylation 模序发现
Para_01
与可变多聚腺苷酸化位点相关的显著基序是在一个类似于用于发现剪接调控基序的过程中获得的。 我们计算了来自 Gasperini 数据集58的每个基因最远端可变多聚腺苷酸化位点相对于多聚腺苷酸化统计量（可变多聚腺苷酸化位点覆盖率的对数比）的组织合并梯度。 这些梯度被归一化为每种基因的最大标准差（截断后的）。 最后，使用一个以每个基因3'非翻译区中显著性模式为中心的192碱基对窗口来提取短梯度片段。 这些梯度片段被用作TF-MoDISco的假设评分，TF-MoDISco的执行采用了与用于剪接基序发现相同的自定义参数。
Attention matrix visualization
注意力矩阵可视化
Para_01
我们通过自注意力层的注意力评分矩阵直接可视化了Borzoi学习到的高阶结构和长程相互作用。 这类高阶结构的例子包括内含子和外显子区域、UTR、启动子和基因跨度。 长程相互作用描述了Borzoi学习到的这些结构之间的关系或依赖性，这将在注意力矩阵中表现为非对角线强度。 这样的例子包括一种现象，其中内含子关注它最近的外显子接合点，3' UTR关注它的多聚腺苷酸化信号（PAS），或者基因跨度关注启动子和转录增强子。 在探索了几个不同位点的预测注意力图谱后，我们注意到与GENCODE注释匹配的高阶结构通常出现在较后的自注意力层。 然而，为了减轻可能的实验特定偏差，并专注于通用知识，我们决定不使用最后两个注意力层，而是使用倒数第二的两个自注意力层进行所有分析。 我们进一步注意到，不同的注意力头往往捕获了大致相同的趋势，因此我们分析了所有八个注意力头的平均注意力值。错误！！！- 待补充
Fine-mapped eQTL classification and regression tasks
精细映射的eQTL分类和回归任务
Para_01
eQTL研究提供了有价值的数据，用于评估Borzoi是否能够识别驱动表达的正确核苷酸及其对特定替代等位基因的敏感性。 我们研究了GTEx（v.8）eQTL结果中的49种不同样本量的组织。 我们利用了之前发表的一篇文章中使用SuSiE生成的汇总统计和精细映射结果。 只有后验因果概率（PIP）≥0.9的精细映射因果eQTL被保留为阳性。 因为我们只关注单核苷酸变异，所以将插入和缺失（indels）排除在外，因为它们引入了由于预测边界移动而产生的技术方差，这是我们未来工作希望解决的问题。 为了可视化具有或不具有感兴趣次要等位基因的个体的测量RNA-seq覆盖轨迹，我们还利用了通过dbGAP获得的GTEx受试者的全基因组测序基因分型数据（http://www.ncbi.nlm.nih.gov/gap）。错误！！！- 待补充
Para_03
对于 L2 分数向量，我们计算预测值 y(ref) 和 y(alt) 在长度轴上每条轨迹的差异的 L2 范数。 在应用 L2 范数之前，我们首先对覆盖轨迹的箱型数据进行对数变换，以关注倍数变化而不是绝对变化。 最终的度量是通过上述方法计算得出的。
Para_04
L2得分提取了更多信息，在此任务上Borzoi的表现更好。 所有先前的Enformer工作都使用SUM得分，但我们观察到它也受益于L2，尽管程度不及Borzoi。错误！！！- 待补充
Para_06
为了完成第三项任务，我们评估了Borzoi识别由eQTL影响的基因的能力，这些基因是从局部基因集中选出的，目的是估计模型在更广泛的GWAS位点上优先正确基因的准确性。 我们从eQTL目录（第5版）下载了49个GTEx组织的精细映射eQTL可信集及其相关的e基因。 可信集文件是从ftp://ftp.ebi.ac.uk/pub/databases/spot/eQTL/credible_sets/下载的（例如，ftp://ftp.ebi.ac.uk/pub/databases/spot/eQTL/credible_sets/GTEx_ge_adipose_subcutaneous.purity_filtered.txt.gz）。
Para_07
注意：这些文件路径已经更改，但历史版本可以在https://github.com/eQTL-Catalogue/eQTL-Catalogue-resources/blob/00ea8a7abca895f26c3aee74ece1307dc5054ace/tabix/tabix_ftp_paths.tsv找到。若要使用最新文件路径表下载可信集，请使用列‘ftp_cs_path’（例如，对于皮下脂肪，下载文件ftp://ftp.ebi.ac.uk/pub/databases/spot/eQTL/susie/QTS000015/QTD000116/QTD000116.credible_sets.tsv.gz）。
Para_08
对于可信集合中的每个变异，我们预测了基因特异性的L2得分，该得分仅考虑与基因外显子重叠的序列位置，范围是在变异中心的360,448 bp序列窗口内的所有基因。 对于每个可信集合，我们通过平均变体的基因得分来计算每个周围基因的单个得分，该变体得分为后验因果概率加权。 对于每个GTEx组织，我们使用匹配的GTEx RNA-seq轨迹的模型预测来计算变体的L2得分。 我们只分析与蛋白质编码基因相关的可信集合。 由于上述插入缺失挑战，我们进一步移除了其中精细映射变体（PIP > 0.1）是插入或缺失的可信集合。 我们预测可信集合的目标基因为该集合中具有最高累积PIP加权L2得分的基因。 作为基线，我们将可信集合的目标基因预测为最近的基因。 我们定义‘最近的基因’为与可信集合具有最大PIP加权逆距离的基因。 最大化PIP加权逆距离的表现优于先前描述的最小化PIP加权距离的方法。 值得注意的是，单一远端可信集合变体可以增加最小平均距离统计量，从而导致错误的eGene预测，而最大化逆距离不会导致此问题。
Fine-mapped paQTL classification task
精细映射的paQTL分类任务
Para_01
我们使用通过 txrevise 处理从 eQTL 目录获得的精细映射的 3′ QTL（本文中称为多聚腺苷酸化 QTL）来评估 Borzoi 预测改变 mRNA 3′ 异构体相对丰度的遗传变异的能力。 这些精细映射结果的文件路径是从 https://github.com/eQTL-Catalogue/eQTL-Catalogue-resources/blob/master/tabix/tabix_ftp_paths.tsv 获取的。
Para_02
表格行通过研究 = ‘GTEx’ 和 quant_method = ‘txrev’ 进行了筛选。所得的汇总统计文件（例如，‘GTEx_txrev_adipose_subcutaneous.all.tsv.gz’）被转换为精细映射文件（‘GTEx_txrev_adipose_subcutaneous.purity_filtered.txt.gz’），并从 ftp://ftp.ebi.ac.uk/pub/databases/spot/eQTL/credible_sets/ 下载（例如 ftp://ftp.ebi.ac.uk/pub/databases/spot/eQTL/credible_sets/GTEx_txrev_adipose_subcutaneous.purity_filtered.txt.gz）。
Para_03
注意：这些文件路径已经更改，但可以找到历史版本的文件路径表在https://github.com/eQTL-Catalogue/eQTL-Catalogue-resources/blob/00ea8a7abca895f26c3aee74ece1307dc5054ace/tabix/tabix_ftp_paths.tsv。 要使用最新的文件路径表下载可信集，请使用列 'ftp_cs_path'（例如，对于皮下脂肪，下载文件ftp://ftp.ebi.ac.uk/pub/databases/spot/eQTL/susie/QTS000015/QTD000119/QTD000119.credible_sets.tsv.gz）。
Para_04
为了构建不是任何 txrevise 可信组一部分的 GTEx SNP 的阴性集合，我们从文件路径表中获取了 quant_method = ‘ge’ 的行，并从 ftp://ftp.ebi.ac.uk/pub/databases/spot/eQTL/sumstats/GTEx/ge/ 下载了完整的汇总统计文件（例如 ftp://ftp.ebi.ac.uk/pub/databases/spot/eQTL/sumstats/GTEx/ge/GTEx_ge_adipose_subcutaneous.all.tsv.gz）。 这些文件路径也发生了变化；要使用最新的文件路径表下载汇总统计文件，请使用列 ‘ftp_path’（例如，对于 adipose_subcutaneous，下载文件 ftp://ftp.ebi.ac.uk/pub/databases/spot/eQTL/sumstats/QTS000015/QTD000116/QTD000116.all.tsv.gz）。
Para_05
对于给定组织的精细映射因果paQTLs是从相应的精细映射文件（'XYZ.purity_filtered.txt.gz'）中通过过滤包含子串'.downstream.'的分子_traid_id行、SNP距离基因跨度（GENCODE v.41）不超过50 bp、与PolyADB（v.3）44中最近注释的3′UTR PAS的距离不超过10,000 bp且PIP大于等于0.9获得的。 有效的负样本从该组织的汇总统计文件（'XYZ.all.tsv.gz'）中获取，并使用与精细映射paQTLs相同的基因跨度和PAS距离过滤条件。 负样本SNP必须完全不在所有可信集合中或在所有GTEx组织中的PIP小于0.01。 最后，我们通过对每个精细映射的因果paQTL选择一个负样本SNP，要求它们与注释的PAS具有相同的距离，并且负样本SNP位于表达水平在（且低于）paQTL基因表达水平1.5倍之内的基因中（在同一组织内）。这导致保留了1,058个唯一的精细映射因果paQTLs。 为了高效地搜索满足这些要求的负样本SNP，采用了以下程序：
Para_06
将所有基因（GTEx v.8）的log2（TPM）值离散化并分桶，每个桶的大小为0.4（在对数2空间中）。 对于给定的目标基因（及其相应的log2（TPM）值），选取所有映射到同一桶中的候选基因。 扫描这一组基因子集，寻找包含距离匹配非因果SNP的基因。 如果桶中没有合适的候选基因（即没有基因具有距离匹配的非因果SNP），则从目标log2（TPM）值中减去0.15，并选取所有被分配到新桶中的候选基因（如果减去0.15后桶未改变，则跳转至步骤4）。 扫描这个新的桶，寻找合适的基因。 如果没有找到合适的基因，则重复步骤3，但这次向原始log2（TPM）值中加上0.15，而不是减去0.15。 扫描这个（可能）新的桶，寻找合适的基因。 如果仍然没有找到合适的基因，则退出并报错（无法匹配）。
Para_07
两个基因之间的最大对数2倍变化，如果它们仍然匹配，则为0.4 + 0.15 = 0.55（约1.464倍）。 根据这些参数设置，每个桶至少包含100个基因，并且我们从未因错误而退出第5步。
Para_08
由于精细映射的paQTL数量相对较少，我们决定将所有组织合并起来而不是分别针对每个组织进行基准测试。 鉴于许多阳性结果在不同组织之间是共享的（总共有1,058个独特的paQTL，每个至少出现在一个组织中），我们在组织合并后得到了大约2.5倍的独特阴性SNP数量。 因此，在基准测试中，我们进行了100次随机匹配的排列，每次从不同组织中选择一个有效的阴性SNP与每个相应的阳性SNP配对，并评估了每组1,058个阳性SNP和1,058个采样阴性SNP的性能。
Para_09
内含子paQTLs（和匹配的阴性样本）是从上述相同的文件中提取的，但它们必须位于内含子区域，并且比任何3′非翻译区多聚腺苷酸化位点更接近注释的内含子多聚腺苷酸化位点。阴性样本现在通过与最近的内含子PAS的距离进行匹配。总共保留了567个精细映射的因果内含子paQTL。 A total of 567 fine-mapped causal intronic paQTLs were retained.
Polyadenylation variant effect prediction
多聚腺苷酸化变异效应预测
Para_01
我们通过计算Borzoi预测的RNA覆盖轨迹的最大比值来计算多聚腺苷酸化中心变异效应评分，该比值是在同一基因的UTR内的任何已注释的3'切割位点之间覆盖倍数变化的最大值。 具体而言，我们将输入窗口的中心定位在SNP上，给定参考和替代等位基因序列x^(ref)和x^(alt)作为输入，预测覆盖轨迹y^(ref)=M(x^(ref)), y^(alt)=M(x^(alt))∈R^(16,384×7,611)，并计算统计量u(y^(ref), y^(alt))_t如下：
Para_02
上述方程中的K表示UTR内的总PAS数量。 有序集合({\mathcal{B"}}={b_{1},\ldots ,b_{K"}})是与K个PAS重叠的y中的bin索引。 基准测试中使用的最终分数是从Borzoi的所有89个GTEx覆盖率轨迹计算出的所有统计量的平均值。 该分数还通过对输入格式为前向补码和反向补码的两种模型复制进行平均得到。
Comparison to APARENT2 and Saluki
与APARENT2和Saluki的比较
Para_01
我们通过两种方式比较了Borzoi的分类性能与APARENT2（v.1.0.2）：首先，我们使用APARENT2对受变异影响的参考和替代PAS序列进行评分，并简单地将预测的对数几率比的绝对值作为变异效应评分。 其次，我们使用预测的几率比来缩放组织合并的参考PAS使用率（如PolyADB中报告的），并使用PAS使用率差异的绝对值作为最终的变异效应评分。后一种统计方法有效地降低了变异的影响程度，根据APARENT2的预测，这种变异具有较大的预测倍数变化，但根据测量结果，它发生在低频使用的PAS上（由于存在竞争性PAS）。
Para_02
当将性能与由APARENT2和Saluki（v.1.0.0）组成的集成模型进行比较时，在paQTL分类任务上，我们遵循了APARENT2论文中的方法。 简而言之，我们根据PolyADB的注释整理每个基因的多聚腺苷酸化信号（PAS）序列和相应的mRNA异构体（最多30个），并拟合了一个逻辑回归模型，用于预测组织汇集的远端异构体比例（如PolyADB中所报告）。 作为输入，该模型使用了APARENT2的PAS评分（最多30个标量）和Saluki的异构体评分（最多30个向量，提取自Saluki的倒数第二层的前四个主成分分析组件）。 利用这个校准后的集成模型，我们预测在诱导特定变异（可能影响多个PAS和异构体序列）时，基因的参考和替代远端比例。 我们从预测的远端比例估计最终的比值比，并使用该比值比基于测量到的参考远端比例重新计算替代远端比例。 最后，我们将替代远端比例从参考比例中减去，并使用该差值的绝对值作为最终的变异效应评分。
Fine-mapped sQTL classification task
精细映射sQTL分类任务
Para_01
精细映射的sQTLs和匹配的阴性样本是从eQTL目录中获得的，使用了与paQTL分类任务相同的汇总统计文件和精细映射文件。 通过筛选包含子字符串‘.contained.’的molecular_trait_id行来提取精细映射的因果sQTLs。 这些QTLs进一步筛选条件为PIPs ≥ 0.9，并且与注释剪接点（GENCODE v.41）的最大距离≤10,000 bp。 每个组织中构建了一组表达匹配和距离匹配的阴性样本，方法与paQTL任务相同，但不包括按最近距离到剪接点进行匹配。 我们总共保留了4,105个独特的精细映射的因果sQTL SNP。
Splicing variant effect prediction
剪接变异效应预测
Para_02
索引 bstart 和 bend 指向 y 中重叠于基因跨度起始和结束位置的区间。 相对较多的精细映射因果 sQTL 数量使得组织特异性基准比较成为可能。 为此，对于给定的 SNP 和 GTEx 组织，我们仅针对与该组织相对应的跟踪子集计算统计量的平均值。
Comparison to Pangolin
与穿山甲的比较
Para_01
我们使用预装的命令行工具通过Pangolin（版本1.0.1）来评分sQTL SNP。 为了便于比较，我们将程序修改为输出六位小数而不是两位。 我们使用以下命令来评分正向和负向vcf文件：pangolin -d 2,000 -m False < sqtl文件 > .vcf hg38.fa gencode41_basic_nort_protein.db < 输出目录 >。
Para_02
尽管该命令最多允许从注释的剪接连接处偏离2,000 bp，但Pangolin还会对变异位置处潜在的新剪接增益进行评分，这意味着该命令将为所有变异（即使它们与剪接位点的距离超过2,000 bp）产生变异效应评分。 我们解析了命令行输出，并将Pangolin输出的基因标识符与SNP发生的基因进行了匹配。最终的变异效应评分被计算为预测的最大增加和减少的绝对值之和。
Splice site identification task
剪接位点识别任务
Para_01
识别DNA序列中的剪接位点构成了成功解释剪接代码和优先考虑致病性剪接变异的基础15,16。 为了评估Borzoi识别剪接位点的能力，我们构建了一个类似的分类任务，并将其与Pangolin16进行了比较。 我们从recount3下载了所有GTEx样本的剪接连接计数，并选择了覆盖量超过对齐读取计数第50百分位的注释连接的正例。 我们过滤了这些集合，只保留那些落在Pangolin和Borzoi测试集交集中的实例。 对于每个正例，我们选择了一个匹配的负例位点，该位点具有相同的三核苷酸环境，距离在100bp到2000bp之间，并且没有证据表明自身是一个剪接连接。 对于Borzoi，我们将每个位点的得分定义为预测的外显子-内含子覆盖度比值的对数值，在对应GTEx组织的样本中取平均。 对于Pangolin，我们使用其预测的剪接位点概率作为得分，并在所有组织中取平均。
Classifying rare and common variation from gnomAD
对gnomAD中的罕见和常见变异进行分类
Para_01
我们从 GnomAD（v.3.1）数据库（https://gnomad.broadinstitute.org）中抽取了一组14,198个单例和14,198个匹配的常见变异体（等位基因频率> 5％），抽样限制在与ENCODE候选顺式调节元件重叠的区域。 为了控制序列突变性，我们排除了位于CpG岛和低复杂度区域内的变异体。 对于每个抽取的单例，我们抽取了一个负样本作为匹配的常见变异体，该变异体与单例具有相同的参考和替代等位基因。 我们还匹配了变异体的背景DNA环境，抽取了位于单例相同三核苷酸内的常见变异体。 最后，我们移除了与基因编码序列中的基因外显子重叠的变异体，仅关注调控变异体进行评估。 对于所有抽取的变异体，我们使用了它们从GnomAD（v.3.1）数据集中获得的CADD原始得分和CADD phred得分（v.1.6）。 我们训练了岭回归模型来区分常见变异体和单例，并使用十折交叉验证来评估模型。 基于CADD的模型使用CADD得分作为特征，而基于Borzoi的模型则使用所有RNA-seq轨迹上的L2得分作为特征，四个模型副本的平均值。 我们通过平均基于Borzoi的模型和基于CADD的模型预测的变异体排名得出了第三个（组合）模型。 为了进行第二次全基因组基准测试，我们从整个基因组中均匀抽取变异体，而不是将变异体抽取限制在ENCODE候选顺式调节元件内。 这导致变异集包含17,360个单例和17,360个匹配的常见变异体。
Predicting TRIP expression
预测TRIP表达
Para_01
我们从先前发表的一篇论文的补充材料中下载了TRIP插入坐标，并测量了七个不同启动子的表达水平。启动子序列列于表S1中，插入坐标（以及测量结果）列于该论文的数据S2中。 为了预测TRIP报告基因的活性，我们遍历了每个启动子序列和坐标，在插入坐标上将524 kb的输入窗口居中并对序列进行插入。 当从Borzoi的RNA-seq或CAGE预测中得出统计数据时，我们将整个TRIP报告基因插入到基因组位置（包括启动子序列、GFP编码区、PAS以及PiggyBac末端重复区域）。 相比之下，当从Borzoi的DNase或组蛋白修饰轨迹（例如H3K4me3）中得出统计数据时，我们只插入启动子，因为当插入整个报告基因时，这些预测变得稍微更差。 我们认为这种现象是由于报告基因两侧的PiggyBac可转座元件引起的，Borzoi本身不擅长预测这些元件，这可能是由于在原始训练数据处理过程中剪切了不可映射的区域。
Para_02
根据预测的覆盖度 ({\boldsymbol{y"}}={\mathcal{M"}}({\boldsymbol{x"}})\in {\mathbb{R"}}^{16,384\times T"})（例如，K562 DNase 轨道），我们通过平均覆盖轨道来计算一个标量预测 (u({\boldsymbol{x"}})\in {\mathbb{R"}})，在插入位点为中心的局部窗口大小为 W 的范围内聚合信号，并应用对数2变换：
Data availability
Para_01
处理后的Borzoi训练数据（包括独热编码序列和覆盖轨迹）可在'gs://borzoi-paper/data'（Google Cloud Storage）下载。 基因注释是从https://www.gencodegenes.org（版本41）获得的。 人类变异数据是从gnomAD（版本3.1）（https://gnomad.broadinstitute.org）获得的。 多聚腺苷酸化位点的注释是从PolyADB（版本3.2）（https://exon.apps.wistar.org/polya_db/v3）和PolyASite（版本2.0）（https://polyasite.unibas.ch/atlas）获得的。 CRISPRi数据是从Nasser等人（2021年）65和Gasperini等人（2019年）的数据从GEO访问号GSE120861获得的。 DNase-seq、ChIP-seq和RNA-seq数据是从ENCODE（https://www.encodeproject.org）下载和处理的；有关RNA-seq实验的详细信息和统计数据，请参见ENCODE门户。 GTEx个体的处理过的RNA-seq样本是从recount3（https://rna.recount.bio）下载的。 CAGE数据是从FANTOM5（https://fantom.gsc.riken.jp/5）下载的。 ATAC-seq数据是从CATlas（http://catlas.org/catlas_hub）下载的。 用于训练的所有实验，包括它们的独特标识符，在https://storage.googleapis.com/seqnn-share/borzoi/hg38/targets.txt中列出了人类样本，在https://storage.googleapis.com/seqnn-share/borzoi/mm10/targets.txt中列出了小鼠样本。 精细映射的eQTLs是从Wang等人（2021年）的补充材料中获得的。 精细映射的eQTL可信集和其他QTLs（sQTLs和paQTLs）是从eQTL目录（https://www.ebi.ac.uk/eqtl）下载的。 本研究中使用的正（精细映射的因果）和负eQTL、sQTL和paQTL集合可在'gs://borzoi-paper/qtl/'（Google Cloud Storage）下载。 TRIP数据是从Leemans等人（2019年）的补充材料中下载的。
Code availability
Para_01
用于训练RNA测序深度学习模型的代码库，包括使用该模型的示例代码以及变异评分的脚本，在 https://github.com/calico/borzoi 下以Apache 2.0开源许可证提供。 预训练的Borzoi模型权重可以通过GitHub获取。 一个单独的GitHub仓库（同样以Apache 2.0开源许可证授权）包含与手稿中分析和结果相关的代码，位于 https://github.com/calico/borzoi-paper 下。
