GPN-Star: 利用系统发育信息指导的基因组语言模型预测跨进化时间尺度的功能约束

应该把图放到对应文字部分，其他保持不变：加州大学伯克利分校统计系、加州大学伯克利分校计算机科学系、杰克逊实验室

德国癌症研究中心（DKFZ）计算基因组学与系统遗传学部

创新基因组学研究所，加州大学伯克利分校，美国加利福尼亚州伯克利市

摘要

基因组语言模型（gLMs）已成为一种强大方法，可直接从DNA序列中学习全基因组范围内的功能约束。然而，源自自然语言处理的标准gLMs通常需极大模型规模与计算资源，却仍在预测任务中落后于经典进化模型。本文提出GPN-Star（融合物种树与比对表征的基因组预训练网络），这是一种具有生物学基础的gLM，其系统发育感知架构显式利用全基因组比对与物种树建模进化关系。该模型在涵盖脊椎动物、哺乳动物及灵长类进化时间尺度的比对数据上训练，在人类基因组编码区与非编码区的各类变异效应预测任务中均达到当前最优性能。跨时间尺度分析揭示：针对不同任务，建模较近或更深层进化各有优势。为展示其推动人类遗传学研究的潜力，我们证明GPN-Star在致病性变异与精细定位GWAS变异的优先排序上显著优于既有方法；实现复杂性状遗传力前所未有的富集；并提升罕见变异关联检验的统计效能。进一步拓展至人类以外物种，我们为五种模式生物——小鼠（Mus musculus）、鸡（Gallus gallus）、黑腹果蝇（Drosophila melanogaster）、秀丽隐杆线虫（Caenorhabditis elegans）和拟南芥（Arabidopsis thaliana）——训练了GPN-Star，验证了该框架的鲁棒性与泛化能力。综上，这些结果确立GPN-Star为一种可扩展、强大且灵活的新型基因组解读工具，能充分挖掘日益丰富的比较基因组学数据。

---

通讯作者联系邮箱：yss@berkeley.edu

---

引言

生物学的一项根本挑战在于理解遗传变异的功能意义。尽管基因组学实验技术取得巨大进展，确定哪些变异影响表型或导致疾病仍十分困难。进化为此问题提供了宝贵洞见：有害突变往往被自然选择清除，而可耐受或有利的改变则可能累积。利用数百万年来此类自然实验，可获得一份全基因组、体内水平的功能约束记录，为变异解读提供独特且高度信息丰富的视角。这种“向进化学习”的思想源远流长，最早可追溯至分子进化的数学模型及早期比较序列分析工作 ( \left\lbrack {1,2}\right\rbrack ) 。

基因组语言模型（gLMs）近期作为一种颇具前景的方法崭露头角，可直接从原始DNA序列中提取全基因组尺度的进化信息（参见[3]及其中引用文献）。这类大规模深度学习模型以无标签序列为训练数据，采用自监督目标进行训练——这一范式推动了人工智能与机器学习（ML）领域的近期突破，尤其在自然语言处理领域[4, 5]。通过学习从序列上下文预测被掩码的核苷酸，gLMs为每个位点分配似然值，该值定量反映各等位基因所受的功能约束，且无需依赖带标签变异体的监督信号。这些似然值已被证实可有效预测全基因组变异效应[3, 6]。然而，基于标准语言建模框架的gLMs在某些变异解读任务上仍逊于结构简单得多的经典系统发育模型——尤其是在人类等复杂真核生物基因组及增强子等远端调控元件中[7]——即使模型规模极为庞大( \left\lbrack {8,9}\right\rbrack )。

对进化数据建模的一种长期沿用方法是构建多序列比对（MSAs）。通过算法将不同生物序列中的同源位点对齐，MSAs可揭示位点特异性的保守模式，并有助于推断特定位置上各类变异的进化偏好。蛋白质MSAs已成为多个著名蛋白质序列深度学习模型的基础，包括AlphaFold[10]、MSA Transformer[11]与EVE[12]。最近，在观察到蛋白质语言模型持续扩大规模所带来的收益日趋递减后，一波基于MSAs（或更广义而言，基于同源搜索）的新方法应运而生[13–20]。将此概念从蛋白质拓展至基因组，数十至数百个物种完整基因组组装体构成的全基因组比对（WGAs）支持开展全基因组尺度的序列进化研究[21, 22]。在WGAs上拟合的经典参数化系统发育模型——如GERP[23]、PhastCons[24]与PhyloP[25]——长期以来一直是人类及其他物种变异解读的核心工具。鉴于这些资源的重要性，若干大型合作项目近期已投入大量资源构建高质量、大规模的WGA数据集[26–28]。

因此，将gLM框架与全基因组比对数据相结合，是一条极具前景的研究方向。我们近期通过GPN-MSA[29]验证了这一潜力：该模型是一种在90个脊椎动物物种WGA数据上训练的Transformer模型，在多项变异效应预测任务中均展现出优异性能。然而，其采用的是适用于未比对序列的通用架构，限制了模型充分挖掘比对数据中固有的丰富结构信息（尤其是物种间的进化关系）的能力。此外，实证表明，若在训练中纳入与人类亲缘关系过近的物种，会导致变异效应预测性能下降，原因可能是模型所学概率分布过度偏向这些物种的基因组。这促使我们从训练数据中排除绝大多数灵长类物种——这一权宜之计并不理想，因为灵长类为理解与人类密切相关的近期进化约束提供了宝贵洞见。最后，GPN-MSA高度适配于某一特定比对，因而难以灵活适应新数据。

本文提出GPN-Star（融合物种树与比对表征的基因组预训练网络），一种通用而强大的gLM框架；其创新性模型架构具备系统发育感知能力，从而更高效地利用多物种WGA信息。这一灵活框架使我们得以将GPN-Star应用于三套分别对应脊椎动物、哺乳动物与灵长类进化时间尺度的人类基因组WGA数据。GPN-Star在涵盖编码区与非编码区的广泛变异解读基准测试中均达到当前最优性能。其中，脊椎动物模型在相同比对数据上持续优于GPN-MSA；而哺乳动物与灵长类模型则在非编码变异效应预测方面尤为突出。我们进一步通过分析复杂性状的遗传力，评估其在人类遗传学中的实用性。值得注意的是，基于灵长类比对训练的GPN-Star模型，在百余种复杂性状中实现了前所未有的遗传力富集水平[30]。此外，我们发现：对约束预测最具信息量的进化时间尺度，与性状的有效多基因性之间存在显著关联。

为进一步展示本方法的普适性，我们在五种模式生物——小鼠（M. musculus）、果蝇（D. melanogaster）、鸡（G. gallus）、秀丽隐杆线虫（C. elegans）与拟南芥（A. thaliana）——上应用GPN-Star，仅做极少调优，即证明其在评估这些物种变异效应方面的有效性。

结果

基于比对与系统发育信息的基因组语言模型

图1：GPN-Star概览。（A）GPN-Star模型架构示意图。模型输入为全基因组比对窗口；目标序列与源序列均由此窗口构建。源序列经注意力池化压缩为支系层级嵌入，池化过程遵循物种树；目标序列则通过多层GPN-Star编码器模块编码，其中系统发育信息引导的交叉注意力模块依据物种树推断的物种间进化距离，整合源序列信息；最终，分类层将编码后的嵌入转换为目标序列各位点的核苷酸概率，用于计算训练损失或预测变异影响。模型完整描述见“方法”部分。（B）GPN-Star在人类基因组上的应用。分别基于脊椎动物、哺乳动物和灵长类比对数据训练了三个模型，以学习不同进化时间尺度上的功能约束。Mya：百万年前；cLLR：校准对数似然比（参见“方法”）。（C）脊椎动物、哺乳动物及灵长类GPN-Star模型在下游任务中的综合性能。对于ClinVar、COSMIC、OMIM、HGMD及GWAS精细定位数据集，性能指标为精确率-召回率曲线下面积（AUPRC）；对于ProteinGym，指标为各实验中Spearman相关系数的均值ρρ；对于S-LDSC，指标为遗传力富集度。所有性能指标经线性缩放至0–1区间：其中0定义为三个GPN-Star模型及在相同三个进化时间尺度上拟合的PhyloP与PhastCons中该任务最低指标值，1定义为最高指标值。

GPN-Star 通过利用多物种系统发育树和全基因组比对数据中蕴含的进化信号，学习目标基因组中遗传变异的功能约束。受经典进化模型启发，该模型旨在刻画每个基因组位点在不同物种间的演化模式。我们提出一种专用的 Transformer 架构，以捕捉高度表达性的、描述进化信息的函数，同时有效整合每个位点的基因组上下文（图1A）。

该模型采用仅编码器架构，并以掩码语言建模（MLM）为目标进行训练。每个输入包含一个跨越多个物种的全基因组比对（WGA）窗口及对应的物种树。基于该比对，我们构建两组序列：目标序列与源序列，每条序列对应一个物种。模型学习在自身序列上下文及源序列提供的进化上下文条件下，预测目标序列中被掩码的核苷酸。这一过程由若干编码器模块堆叠实现；每个模块包含：用于编码序列内上下文的序列级自注意力模块、用于编码源序列进化上下文的系统发育信息引导的交叉注意力模块，以及用于整合信息的前馈网络。模型核心的交叉注意力模块，通过注意力机制，依据物种树推导出的进化距离，自适应地加权各源序列的贡献。由此生成的目标序列编码嵌入，融合了丰富的序列内上下文与进化信息，并经由 MLM 头输出最终的单核苷酸概率预测，用于训练及下游变异评分。完整架构与训练细节见“方法”部分。

GPN-Star 框架在三个关键方面超越了 GPN-MSA：第一，GPN-MSA 训练时仅掩码人类基因组，而 GPN-Star 则利用比对中全部基因组，跨多个物种预测被掩码的核苷酸，从而大幅增加训练数据的数量与多样性；第二，GPN-Star 通过专用注意力模块显式整合物种间的系统发育关系，实现更准确、更具生物学基础的建模；第三，它可灵活适配任意组成与规模的比对，无需人工筛选（例如剔除亲缘关系过近的物种），而此类操作在 GPN-MSA 中是必需的 [29]。实证表明，我们在广泛评估中持续观察到显著的性能提升（补充图1）。

补充图1：与GPN-MSA在变异效应预测上的直接对比。基准测试与图2相同。为确保公平比较，我们展示了参数量相同（85M）且使用未校准原始分值的GPN-Star性能。

本研究另一方法学创新在于首次在基因组语言模型（gLMs）中纳入突变率变异校正。我们开发了一种高效校准流程，从原始模型输出中去除突变率影响，从而实现对全基因组功能约束更准确的量化（详见下文及“方法”部分）。

GPN-Star 框架具有通用性与灵活性，可适配任意物种的任意比对数据，仅需极少超参数调优即可实现稳健性能（“方法”部分）。我们首先将该框架应用于人类基因组，分别使用目前最大规模的脊椎动物、哺乳动物与灵长类全基因组比对（WGA）训练了三个独立的 GPN-Star 模型（图1B，“方法”部分）。本研究中，我们将这些模型分别记为 GPN-Star (V)、GPN-Star (M) 和 GPN-Star (P)。我们尝试过不同模型规模，但本文聚焦于最大规模模型——其含 2 亿参数，在 8 块 NVIDIA A100 GPU 上训练数日（补充表1

补充表 1：本研究中所有 GPN-Star 模型训练设置的汇总。

模型名称

基因组组装（大小）

训练染色体

验证染色体

预留测试染色体

全基因组比对（WGA）

比对物种

\( {\phi }_{\text{ clade }} \)

支系数量

上下文长度

模型规模

学习率

训练步数

训练时间（8×A100）

hg38-vertebrate-200m

hg38（3.1 Gbp）

1–20、X、Y

21

22

multiz100way

100 种脊椎动物

0.2

45

128 bp

200 M

5e-5

200 k

4.5 d

hg38-vertebrate-85m

85 M

1e-4

150 k

1.8 d

hg38-vertebrate-25m

25 M

1e-4

150 k

0.9 d

hg38-mammal-200m

cactus447way

447 种哺乳动物

0.05

145

256 bp

200 M

5e-5

200 k

6.6 d

hg38-mammal-85m

85 M

1e-4

150 k

2.9 d

hg38-mammal-25m

25 M

1e-4

150 k

1.2 d

hg38-primate-200m

243 种灵长类动物

21

200 M

5e-5

200 k

4.2 d

hg38-primate-85m

85 M

1e-4

150 k

1.6 d

hg38-primate-25m

25 M

1e-4

150 k

0.8 d

mm39-85m

mm39（2.7 Gbp）

1-17,

18

19

multiz35way

35 种脊椎动物

0.2

26

128 bp

85 M

1e-4

150 k

1.6 d

mm39-25m

X、Y

25 M

1e-4

150 k

0.7 d

galGal6-85m

galGal6（1.1 Gbp）

1-24,

25-28

30-33

multiz77way

77 种脊椎动物

0.2

29

128 bp

85 M

1e-4

150 k

1.7 d

galGal6-25m

Z、W

25 M

1e-4

150 k

0.7 d

dm6-85m

dm6（144 Mbp）

2L、3L、3R、X、Y

2R

4

multiz124way

124 种昆虫

0.2

93

128 bp

85 M

1e-4

\( {50}\mathrm{k} \)

0.8 d

dm6-25m

25 M

1e-4

\( {50}\mathrm{k} \)

0.3 d

ce11-85m

ce11（100 Mbp）

II、IV、V

I

III

multiz185way

135 种线虫

0.2

92

128 bp

85 M

1e-4

50 k

0.8 d

ce11-25m

X

25 M

1e-4

50 k

0.3 d

tair10-85m

tair10（120 Mbp）

1-3

4

5

multiz18way

18 种十字花科植物

0.05

18

256 bp

85 M

1e-4

50 k

0.5 d

tair10-25m

25 M

1e-4

50 k

0.2 d

），资源开销显著低于此前的 gLMs，例如 Nucleotide Transformer（在 128 块 A100 GPU 上训练一个月）[8] 与 Evo-2（在逾 2000 块 H200 GPU 上训练数月）[9]。

与此前在极宽泛进化时间尺度（如从原核生物至人类）数据上训练的 gLMs 不同，GPN-Star 聚焦于更窄、更贴近人类的近期系统发育距离（图1B）。如下文所示，建模更长的进化历史并非总能带来最优效果，具体取决于下游应用场景；相反，捕获近期进化约束对解读特定类别的遗传变异尤为有利（图1C）。

GPN-Star 对致病性编码区变异展现出强大预测能力

人类遗传学中一项关键挑战是理解遗传变异对疾病易感性的影响——这是提升遗传病诊断水平、识别药物靶点、实现精准医学承诺的关键一步。目前已开发多种计算模型，用于预测人类基因组中变异的致病性。本文系统评估了 GPN-Star 在一整套综合基准上的性能。

我们对编码区变异的评估聚焦于错义变异，这是最常见的一类变异。我们首先考虑在ClinVar（一个广泛使用的临床变异数据库，包含专家审定的致病性标签[31]）中对致病性与良性变异进行分类。在所有主流全基因组变异效应预测工具中——包括经典进化方法（PhyloP、PhastCons）、集成模型CADD，以及近期的基因组语言模型（Nucleotide Transformer 2.5B多物种模型、Evo-2 40B模型和GPN-MSA）——GPN-Star（V）在精确率-召回率曲线下面积（AUPRC）上表现最优，性能与蛋白质语言模型ESM-1b相当（图2A）。为便于可视化，图2中仅展示性能最优的GPN-Star模型；全部三个版本的结果见补充图2

补充图2：包含所有竞争模型的变异效应预测基准测试。基准测试与图2相同，但额外纳入GPN-Star、PhyloP和PhastCons（涵盖全部三种进化时间尺度）、Evo-2（两种参数规模）以及Roulette突变率估计。

。

随后，我们考察了来自COSMIC数据库[32]的体细胞癌症变异，构建了一个基准任务，用于区分肿瘤中高频出现的错义变异与一般人群中常见的错义变异（gnomAD[33]）。如图2B所示，GPN-Star（V）显著优于所有竞争模型，展现出对胚系变异之外致病性的强大预测能力。

除临床变异外，我们还在深度突变扫描（DMS）数据上评估了模型性能；DMS是一种功能实验，可测量特定蛋白质中所有可能错义突变的适应度效应。在ProteinGym[34]提供的31个人类DMS数据集中，GPN-Star（V）优于所有全基因组模型，尽管略逊于蛋白质特异性模型ESM-1b（图2D）。

我们进一步将GPN-Star与两种近期错义变异效应预测工具AlphaMissense[35]和PrimateAI-3D[36]进行比较。这些模型均以群体等位基因频率数据为监督信号进行训练——该信号对致病性预测高度有效，亦被用于ClinVar等数据库的标签定义[37]。有趣的是，尽管AlphaMissense与PrimateAI-3D在ProteinGym和ClinVar上均取得顶尖性能，但仅通过gnomAD等位基因频率对GPN-Star（V）预测结果进行简单的事后校正，即可使其在ClinVar上的性能超越PrimateAI-3D，并逼近AlphaMissense（补充图3

补充图3：与AlphaMissense及PrimateAI-3D在错义变异效应预测上的直接对比。该对比包含图2中的编码区基准测试（COSMIC数据集除外），因其阴性集完全位于这两个模型的训练数据内。部分对比的样本量较小，原因是这两种方法提供的预计算分值未能覆盖全部变异。在ClinVar基准测试中，“GPN-Star (V) + gnomAD AF”方法通过将gnomAD v3中等位基因频率( > {2e} - 4 )的变异赋以最良性分值，对GPN-Star (V)预测结果进行调整；该阈值依据AlphaMissense训练数据中良性变异的定义选定。

）。

GPN-Star是当前最先进的致病性非编码变异预测工具

接下来，我们考察非编码变异相关任务；此类变异已知对预测建模尤为具有挑战性。此处我们证明，GPN-Star是识别人类基因组中致病性非编码变异的有力工具。

除既往全基因组进化模型外，我们还纳入三种主流序列到功能模型——Enformer[38]、Borzoi[39]和AlphaGenome[40]；这些模型基于海量功能基因组学数据训练，已被广泛应用于非编码变异解读。我们在OMIM[41]与HGMD[42]数据库中的致病性非编码变异分类任务上评估各模型；这两个数据库均包含专家审定的人类疾病相关变异注释。GPN-Star（M）在两项基准测试中均取得最佳性能（图2E、F）。值得注意的是，序列到功能模型在此任务上的表现显著差于进化模型，这与近期研究发现一致[7, 43]。

鉴于启动子区域在转录起始与基因调控中的关键作用，已有多种专用模型被开发用于预测启动子变异的影响。因此，我们在OMIM中的启动子变异上评估GPN-Star，并将其性能与其他方法（包括三种启动子特异性模型：PromoterAI[44]、SpeciesLM[45]与GPN-Promoter[7]）进行比较。如图2H所示，GPN-Star（M）展现出优于所有竞争模型的预测性能；尤其相较于启动子特异性模型，其性能优势相当显著。

GPN-Star能够识别GWAS精细定位所得的潜在因果变异

全基因组关联研究（GWAS）在鉴定遗传性疾病易感性相关变异方面发挥了关键作用。为进一步评估GPN-Star的实用性，我们在UK Biobank[46]中65个性状的GWAS变异精细定位所得的潜在因果与非因果错义变异分类任务上检验其性能。在所有竞争模型中，GPN-Star（M）在这些精细定位的错义变异上取得了最高的预测性能（图2C）。值得注意的是，尽管AlphaMissense与PrimateAI-3D均利用了群体等位基因频率信息，但在本任务中仍被GPN-Star大幅超越（补充图3）。

随后，我们在来自英国生物银行（UK Biobank）[46]的83个性状的精细定位非编码GWAS变异上评估了各模型，再次检验其区分潜在致病性与非致病性变异的能力。GPN-Star（M）在此基准测试中持续优于所有其他模型（图2G）。序列到功能模型（Enformer、Borzoi和AlphaGenome）表现中等，而Evo-2预测能力相对有限，这与先前观察结果一致[7]。对于位于启动子区域的精细定位变异，GPN-Star（M）再次超越所有模型，包括启动子特异性模型PromoterAI、SpeciesLM和GPN-Promoter（图2H）。

不同变异效应预测任务所对应的相关进化时间尺度

上述结果表明，GPN-Star是一种强大且通用的全基因组变异解读框架。此外，我们的分析凸显了进化模型训练数据所表征的进化时间尺度的重要性。在各项基准测试中，我们观察到一致规律：编码区变异及效应量较大、频率较低的变异，通常更易被基于更深层进化时间尺度训练的模型准确预测；此类变异在ClinVar、COSMIC和ProteinGym基准中富集。相反，非编码区变异及效应量较小、频率较高的变异，总体上更易被基于较浅层进化时间尺度训练的模型准确预测，如OMIM、HGMD及精细定位GWAS基准所示。PhyloP与PhastCons评分在三个进化时间尺度上亦呈现类似趋势，尽管二者在各时间尺度上均持续落后于GPN-Star（图2I）。我们将在“讨论”部分进一步探讨该现象的潜在成因。

GPN-Star提升罕见变异关联检验效能

鉴于GPN-Star在预测致病性变异及精细定位变异方面表现出色，我们进一步探索其在罕见变异关联检验（RVAT）中的应用价值——这是统计遗传学中一项重要但具挑战性的任务。RVAT备受关注，原因在于大多数大效应变异往往较为罕见，因其受到负向选择压力[47]。为克服低频变异带来的统计功效不足问题，RVAT通常在基因水平开展，即对每个基因内的变异进行聚合[48]。此类检验流程常依赖反映功能重要性的变异注释，以在聚合过程中优先纳入关键变异。DeepRVAT[49]是近期提出的一种方法，采用深度集合网络整合此类变异注释用于罕见变异关联检验。通过将富有表达力的深度学习框架与一套强大的变异注释相结合，其在英国生物银行全外显子组测序（WES）数据的大规模评估中，相较既有方法展现出更高的统计功效与计算效率[49]。

我们采用DeepRVAT开展RVAT实验，在训练与关联检验阶段，将发表版本所用注释全面升级为全部三种GPN-Star预测结果。遵循文献[49]的基准测试流程，我们利用DeepRVAT对英国生物银行中34个数量性状开展RVAT实验，所用WES数据来自161,822名无亲缘关系的欧洲裔个体，并仅保留次要等位基因频率（MAF）( < {0.1}% )的变异（方法部分）。为控制随机模型初始化引入的差异，该流程使用不同随机种子重复运行三次。如图2J所示，将GPN-Star模型预测结果作为变异注释加入DeepRVAT后，在家族错误率( < {0.05} )下发现的基因数量增加（各次运行平均为402个基因，高于原方法的383个）。在显著发现中，亦有更多基因-表型关联在基于更大样本量的英国生物银行传统RVAT研究中得到复现（方法部分）( \left\lbrack {{50},{51}}\right\rbrack )（平均353个，高于338个），表明新增发现具有高度稳健性。值得注意的是，原始DeepRVAT注释中已包含AlphaMissense、PrimateAI等当前最先进的错义变异效应预测器，以及序列到功能模型DeepSEA。即便如此，GPN-Star仍提供了互补信息，并进一步提升了检验效能。

GPN-Star预测结果对复杂性状遗传力具有高度信息量

为进一步探究复杂性状全基因组关联研究（GWAS）中因果变异的预测，我们采用了分层连锁不平衡评分回归（S-LDSC）[52]。S-LDSC是一种原理严谨的方法，可在利用全部SNP（包括那些未被高置信度精细定位或未达全基因组显著性的SNP）信号的同时，评估某注释对复杂性状遗传力的信息贡献度。它亦是功能导向精细定位[53]和多基因风险评分[54]的重要基础。由于S-LDSC需对( \sim {10} )百万个变异打分，故仅能评估扩展性最强的模型（或已预先计算好分数的模型）。我们对模型分数进行二值化处理，以选取特定比例的常见变异；在主分析中，我们采用约束性最强的前( {0.1}% )个变异。我们针对每种模型分数分别运行S-LDSC，并同时校正96个基线特征，随后对106个独立性状[30]的结果进行荟萃分析。

自S-LDSC早期工作以来，保守性评分已被证实是对复杂性状遗传力富集程度最高的注释[52]。近期，灵长类特异性保守性——尤其是PhastCons（P）——已成为复杂性状遗传力富集的最先进指标[55, 56]。值得注意的是，GPN-Star（P）显著超越该基准，其次为GPN-Star（M）（图3A）。当考察遗传力系数( {\tau }^{ \star } )（图3A）时，这些提升尤为突出；该系数量化了某注释在校正基线特征后对遗传力的独特贡献（方法部分）。尽管单个性状的富集估计值噪声更大，GPN-Star（P）或GPN-Star（M）始终稳居榜首（补充图4

补充图4：基于20个遗传力最强性状的模型预测值所计算的遗传力富集分析。与图3一致，此处展示的是各模型定义的前0.1%约束性最强变异中的富集情况。

）。鉴于复杂性状遗传力富集长期作为人类遗传学中一项有意义的基准，这一进展尤为显著。

将分析范围分别限定于编码区与非编码区中排名最靠前的变异时，GPN-Star的性能优势依然存在（图3B、C），且在非编码区增益更强。此外，在一系列二值化阈值及全部三个进化时间尺度下，GPN-Star均持续优于既往保守性评分（图3D）。

进化时间尺度与复杂性状的遗传架构

我们进一步分析了基于不同进化时间尺度训练的GPN-Star模型如何对变异进行优先排序。在GPN-Star（P）所优先排序的常见变异中，39%与GPN-Star（M）及（V）共享；另有59%仅与GPN-Star（M）共享，36.6%则为GPN-Star（P）独有（补充图5

补充图5：三种模型各自定义的前0.1%约束性最强变异之间的重叠（仅常见变异）。

；所有S-LDSC变异见补充图6

补充图6：三种模型各自定义的前0.1%约束性最强变异之间的重叠（全部S-LDSC变异，含常见及低频变异）。

）。受以下观察启发——GPN-Star（M）在孟德尔性状上表现更优，而GPN-Star（P）在多基因性状上更胜一筹（图1C、图2I）——我们检验了二者在最新估算的性状多基因性谱系[57]上的相对遗传力富集表现。在27个具备多基因性估计值的性状中，我们发现（对数）有效多基因性与GPN-Star（P）和GPN-Star（M）间遗传力富集差异呈中等程度相关（Pearson’s ( r = {0.47} )）（图3E，补充图7

补充图7：GPN-Star (P) 与 (M) 富集差异随估计有效多基因性变化的趋势。( p )值为单侧。虚线为普通最小二乘拟合线。

）。例如，对于低密度脂蛋白胆固醇（LDL-C）（有效多基因性 = 89），GPN-Star（M）获得更高富集；而对于精神分裂症（有效多基因性 ( = {13},{069} )），GPN-Star（P）富集更高。这表明进化时间尺度与性状遗传架构之间存在实质性关联，并进一步佐证了灵长类特异性进化约束在人类复杂性状变异中所起的独特作用。

为更深入理解最优模型GPN-Star（P）所优先排序的变异类型，我们分析了Ensembl变异后果注释[58]以及ENCODE SCREEN候选顺式调控元件（cCREs，补充表2

补充表2：ENCODE SCREEN cCRE类别。

	TSS 距离	DNase/ATAC	H3K4me3	H3K27ac	CTCF	TF
启动子（PLS）	\( < {200} \)	✓	✓
近端增强子（pELS）	\( < {2000} \)	✓		✓
远端增强子（dELS）	\( > {2000} \)	✓		✓
染色质可及性伴 H3K4me3（CA-H3K4me3）		✓	✓	✘
染色质可及性伴 CTCF（CA-CTCF）		✓	✘	✘	✓
染色质可及性伴 TF（CA-TF）		✓	✘	✘	✘	✓
染色质可及性（CA）		✓	✘	✘	✘	✘
转录因子（TF）		✘	✘	✘	✘	✓

）[59]。常见变异的结果汇总于图3F，聚焦于各占优先排序变异( > 1% )的后果类别。常见变异的完整结果，以及包含低频变异的结果，分别见补充表3

补充表3：根据GPN-Star（P）筛选出的前( {0.1}% )个约束性最强的变异（仅常见变异）的功能后果分布。比值比及( p )值采用双侧Fisher精确检验计算。加粗项：在错误发现率（FDR）<5%下具有统计学显著性。

后果	计数	比例	比值比
错义突变	1983	0.33	170.27
dELS	1564	0.26	2.39
dELS侧翼区	468	0.08	0.35
长链非编码RNA	287	0.05	2.14
3'非翻译区	280	0.05	4.23
pELS	209	0.04	2.14
5'非翻译区	151	0.03	9.83
CA	135	0.02	0.92
内含子	99	0.02	0.07
PLS	96	0.02	7.18
CA侧翼区	92	0.02	0.21
同义突变	77	0.01	4.11
基因间区	67	0.01	0.05
pELS侧翼区	61	0.01	0.57
CA-CTCF	60	0.01	0.94
CA-H3K4me3	49	8e-03	1.07
剪接区	40	7e-03	7.52
CA-CTCF侧翼区	36	6e-03	0.19
剪接受体位点	28	5e-03	50.38
剪接受体区	27	5e-03	23.04
TF侧翼区	24	4e-03	0.14
剪接受体第5位碱基	22	4e-03	52.98
CA-TF	20	3e-03	1.30
剪接多嘧啶区	18	3e-03	4.19
CA-H3K4me3侧翼区	15	3e-03	0.14
剪接受体位点	11	2e-03	30.65
CA-TF侧翼区	10	2e-03	0.26
终止密码子获得	10	2e-03	41.91
转录因子结合位点	9	2e-03	0.15
起始密码子丢失	4	7e-04	68.72
下游基因	3	5e-04	0.15
成熟miRNA	3	5e-04	48.21
终止密码子丢失	2	3e-04	22.64
PLS侧翼区	1	2e-04	0.17

–4及补充图8

补充图8：按GPN-Star (P) 定义的前( {0.1}% )约束性最强变异中各类变异后果的比例，及其相对于99.9%约束性最弱变异（仅常见变异）的优势比。

–9。

错义变异占比最高（33%，富集170倍），其次为具有远端增强子样特征（dELS）区域的变异（26%，富集2.4倍）及其侧翼区域变异（8%，耗竭0.35倍）。远端增强子变异对复杂性状遗传力的重要性不容低估，因其仍是当前序列到功能模型[60]（包括最新的AlphaGenome[43]）以及无比对gLMs（包括规模最大的Evo-2[7]）的主要短板。

三种GPN-Star模型所优先识别的变异类型存在有趣差异。GPN-Star（P）与（M）之间仅存在细微差别，前者在错义变异中富集1.14倍（补充表5

补充表5：根据GPN-Star (P) 与 GPN-Star (M) 分析，在顶部常见变异位点中显著富集的功能后果。

后果	比值比	\( p \)-值	\( q \)-值
错义	1.14	7e-04	2e-02

）。不出所料，GPN-Star（P）与（V）之间的差异则显著更多（补充表6

补充表6：根据GPN-Star（P）与GPN-Star（V）分析，在顶部常见变异中显著富集的后果。

后果	比值比	\( p \)-值	\( q \)-值
同义突变	3.89	4e-09	5e-08
5′非翻译区	1.89	4e-06	2e-05
偏最小二乘法	1.76	1e-03	3e-03
错义突变	1.52	5e-25	2e-23
3′非翻译区	1.36	1e-03	4e-03
增强子样染色质开放区-侧翼区	0.62	7e-15	1e-13
CTCF结合区	0.62	1e-05	6e-05
CTCF结合区-侧翼区	0.56	8e-06	4e-05
CTCF结合区-CTCF结合位点-侧翼区	0.55	5e-03	1e-02
内含子	0.54	1e-06	9e-06
基因间区	0.49	2e-06	2e-05
CTCF结合区-H3K4me3修饰-侧翼区	0.34	2e-04	8e-04

）。GPN-Star（P）总体上在外显子变异（如错义变异富集1.52倍、同义变异富集3.89倍）中富集，而在非外显子变异（如基因间区变异）中耗竭。

组织特异性预测进一步提升遗传力富集

作为复杂性状遗传力的最终分析，我们评估了组织特异性评分的表现。Enformer[38]等序列到功能模型并不提供单一的变异效应评分，而是预测数千个代表不同实验与组织的功能基因组信号轨迹上的活性变化。在我们的比较中，我们利用了Fabiha等人[61]近期生成的Enformer评分：既涵盖跨所有轨迹的谨慎聚合（组织无关型），也涵盖九种特定组织内的聚合。鉴于组织无关型Enformer在全性状荟萃分析中富集程度较低（图3A），我们进一步考察了组织特异性Enformer评分的表现。我们仅在预期该组织相关的性状中对这些评分进行荟萃分析（补充表7

补充表7：我们预期它们相关的性状和组织。

	大脑	血液	肝脏	肠道	肾脏	脂肪	心脏	皮肤	肺
注意力缺陷多动障碍（ADHD）	1 0	0 0	0 0	0 0	0 0	0 0	0 0	0	0
腹主动脉瘤 寻常痤疮 急性阑尾炎	0 0	0 0	0 0	0 1	0 0	0 0	0 0		0 0 0
初产年龄 父母寿命	1 0	0 0	0 0	0 0	0 0	0 0	0 0		0 0
饮酒量／AUDIT评分	1	0	0	0	0	0	0		0
碱性磷酸酶	0	0	1	1	0	0	0		0
阿尔茨海默病 神经性厌食症	1	0	0	0	0 0	0 0	0 0		0 0
天冬氨酸氨基转移酶									1
心房颤动	0	0		0	0	0	1		0
心房颤动与扑动	0	0		0	0	0	1		0
自闭症谱系障碍 体重指数（BMI） 脱发 嗜碱性粒细胞百分比	1 1 0 0	0 0 0 1	0 0	0 0 0 0	0 0 0 0	0 1 0 0	0 0 0 0	0	0 0 0 0
双相情感障碍（所有病例） 乳腺癌（女性） 支气管炎、肺气肿、哮喘、鼻炎、湿疹、医生诊断的过敏 CHIP 慢性阻塞性肺疾病（COPD） COPD－FVC 大麻使用障碍	1 0 0 0 0 0 1	0 0 1 1 0 0 0	0 0 0 0 0 0	0 0 0 0 0 0 0	0 0 0 0 0 0 0	0 1 0 0 0 0 0	0 0 0 0 0 0 0	0	0 0 1 0 1 1 0
乳糜泻		1		1 0	0	0	0 0		0 0
冠状动脉疾病（Aragam）		0		0	0	0	0		0 0
新冠疫苗接种 肌酐		0		0	1	0	0		0
日间小睡	1	0 0	0	0 0	0 0	0	0 0		0 0
抑郁症 舒张压（血压）	0	0	0	0	1	0 0	1		0
憩室病与憩室炎	0 1	0 0	0 0	1 0	0 0	0 0	0 0		0 0
每周饮酒量 受教育年限	1	0	0	0	0	0	0		0
子宫内膜异位症 嗜酸性粒细胞百分比 食管癌 是否曾吸烟 胎儿出生体重	0 0 0 1 0	0 1 0 0 0	0 0	0 0 1 0 0	0 0 0 0 0	0 0 0 0 0	0 0 0 0 0		0 0 0 0 0
胃肠道疾病 一般风险耐受性 青光眼 痛风 高密度脂蛋白（HDL）		0 0			0	0	0 0 0 0		0 0 0 0 0
心力衰竭 跟骨T值 身高 甲状腺功能减退症（自报）		0 0 0 0	0	0 0 0 0	1 0 0 0	0 0 0 0	1 0 0 0		0 0 0 0
炎症性肠病（IBD） 胰岛素样生长因子1（IGF1）		1 0		1 0	0 0	0 0	0 0		0 0
腹股沟疝 失眠 智力		0 0 0		0 0 0	0 0 0	0 0 0	0 0 0		0 0 0
	0	0		0	0 0	0 0	0 0		0 1
系统性红斑狼疮（SLE） 皮肤恶性肿瘤	0 0 0	1 0 0		0 0 0	1 0 0	0 0 0	0 0 0		0 0 0
母亲出生体重 降糖药物使用 初潮年龄	0 1	0 0		0 0	0 0	0 0	0 0		0 0
单核细胞百分比 晨型人 口腔牙齿问题 曾生育子女数	0 1 0 1 1	1 0 0 0 0	0 0 0 0	0 0 0 0	0 0 0 0	0 0 0 0	0 0 0 0 0		0 0 0 0 0
帕金森病 磷酸盐 身体活动量 血小板分布宽度	0 1 0	0 0 0	0 0 0	0 1 0 0	0 1 0 0	0 0 0 0	0 0 0		0 0 0
原发性胆汁性肝硬化									0
原发性开角型青光眼	0	0
	0	0		0	0	0	0		0
前列腺癌	0 0	0 0	0 0	0 0	0 1	0 0	0 0		0 0
红细胞计数 红细胞分布宽度	0	0	0	0	0	0	0		0
反应时间 网织红细胞计数	1 0	0 0	0 0	0 0	0 0	0 0	0 0		0 0
类风湿性关节炎（石垣）	0	1	0	0	0	0	0	0	0
类风湿性关节炎（塞瓦尔斯多蒂尔）	0	1	0	0	0	0	0	0	0
SARS-CoV-2 感染（C2）	0	0	0	1	0	0	0	0	1
精神分裂症 睡眠时长	1 1	0 0	0 0	0 0	0 0	0 0	0 0		0 0
打鼾 中风	1 0	0 0	0 1	0 0	0 0	0 0	0 1	0	0 0
结构连接全局指标（影像学）	1 0	0 0	0	0 0	0 0	0 0	0 0	0	0 0
晒伤 T1D——1型糖尿病	0	1	0	0	0	0	0	1	0
端粒长度 拇指骨关节炎	0 0	0 0		0	0	1	0 0		0 0
总胆红素	0	0		0	0	0	0		0
全髋关节置换术 总蛋白	0 0	0 0	0 1	0 0	0 1	1 0	0 0		0 0
抽动秽语综合征	1	0	0	0	0	0	0		0
2型糖尿病	0 0	0 0	1 0	1 0	0 0	1 0	0 0		0 0
静脉曲张手术 静脉曲张	0	0	0	0	0	0	0	0	0
静脉血栓栓塞症（VTE）	0	0	1	0	0	0	0		0
内脏脂肪组织体积 维生素D	0 0	0 0	0 1	0 0	0 1	1 0	0 0		0 0
腰臀比BMI比值 白细胞计数	0 0	0 1	0 0	0 0	0 0	1 0	0 0	0 0	0 0

）。不出所料，组织特异性Enformer评分始终优于组织无关型Enformer评分（图3G，补充图10

补充图10：组织特异性分析的遗传力富集度与标准化系数。

）。

受此观察启发，我们设计了一种简单方法将组织特异性整合进GPN-Star：仅考虑位于组织特异性基因附近的高分变异（灵感源自LDSC-SEG[62]）。该策略虽提升了部分组织（如脑与血液/免疫系统）的表现，却未能持续超越组织无关型GPN-Star注释的性能，可能因我们整合性状组织特异性的方式尚不理想（图3G，补充图10）。尽管如此，GPN-Star在全部九种组织中均优于Enformer。

聚焦具体性状示例：脑特异性GPN-Star对精神分裂症富集最高；血液/免疫特异性GPN-Star对狼疮表现最佳；而组织无关型与肝特异性GPN-Star对IGF1的预测结果相近（图3H）。这些发现凸显了在变异效应预测中纳入组织特异性信息的重要性。

GPN-Star学习基因组中的功能元件及其依赖关系

gLMs可通过预测被掩蔽核苷酸来学习强大的序列表征，因为某区域的局部核苷酸分布深受其功能影响。尽管先前研究已表明gLM嵌入可在无监督下区分基因组区域[6]，但基于比对的模型是否具备同等能力仍不明确。为此，我们可视化了来自基因区、cCREs及背景区域的基因组窗口嵌入（图4A）。此处重点展示含( {85}\mathrm{M} )参数的GPN-Star（M）模型，并将其余模型置于补充材料中。

所得嵌入在很大程度上按基因组区域分离（图4A，补充图11

补充图11：窗口嵌入的可视化（所有模型）。行表示在不同进化时间尺度上训练的模型（V：脊椎动物，M：哺乳动物，P：灵长类），列表示不同模型规模（参数数量）。

）。例如，编码序列（CDS）与增强子各自形成独立簇。另一方面，启动子与5' UTR窗口难以区分，可能因5' UTR嵌入同时编码转录与转录后信号。相比非保守窗口，保守窗口的嵌入在功能区域上表现出更强的聚类性（图4B）。这些结果表明，GPN-Star在做出预测时已感知基因组的核心功能元件。

位点无关模型（如phyloP）按定义无法捕捉核苷酸间的依赖关系。相比之下，无需比对的生成式语言模型（gLMs）已被证实能够学习已知互作元件之间的依赖关系，例如转录因子结合位点（TFBS）基序内部、或剪接供体与受体之间的依赖关系[45]。为进一步探究GPN-Star对基因组语法的理解，我们通过系统性地突变每个位置，并量化序列中其他位置对应的概率变化，分析了模型所学得的核苷酸依赖关系( \left\lbrack {{45},{63}}\right\rbrack )。在编码酪氨酸羟化酶的( {TH} )基因的启动子区及第一外显子中，该分析揭示了编码区内一个强相互作用模块，以及转录因子CREB结合位点处另一模块；已知该位点突变可导致酪氨酸羟化酶缺乏症和肌张力障碍[44, 64–67]（图4C，补充图12

补充图12：( {TH} )启动子区的核苷酸依赖性图谱（脊椎动物模型被排除，因该区域超出其上下文长度）。行表示在不同进化时间尺度上训练的模型（V：脊椎动物，M：哺乳动物，P：灵长类），列表示不同模型规模（参数数量）。

）。模型预测：CREB的功能依赖于TATA框，且有趣的是，还依赖于编码区。

我们在HBA1基因中也观察到跨外显子的预测依赖关系（补充图13

补充图13：HBA1第一个外显子周边区域的核苷酸依赖性图谱（脊椎动物模型被排除，因该区域超出其上下文长度）。行表示在不同进化时间尺度上训练的模型（V：脊椎动物，M：哺乳动物，P：灵长类），列表示不同模型规模（参数数量）。

）；该基因具有异常短小的内含子，完全落入模型上下文长度范围内。值得注意的是，剪接受体与供体区域之间的依赖关系尤为显著，与既往研究结果一致[45]。

接下来，我们考察了与家族性高胆固醇血症相关的LDLR启动子区；该区域含有多个已知的转录因子结合位点（TFBS），并已被大规模平行报告基因实验（MPRA）深入研究[55, 68, 69]。TFBS的位置可较好地从核苷酸依赖图谱的模块结构中预测得出（图4D，补充图14

补充图14：LDLR启动子区的核苷酸依赖性图谱（所有模型）。行表示在不同进化时间尺度上训练的模型（V：脊椎动物，M：哺乳动物，P：灵长类），列表示不同模型规模（参数数量）。

）[45]。此外，模型还预测出TFBS之间的依赖关系，包括SREBP2与SP1之间已知的经典互作[70]。

最后，我们分析了一个推测受灵长类特异性约束的开放染色质区域[56]。支持该假说的是，GPN-Star（P）模型在推定的TEAD4结合位点周围展现出最高的依赖强度（补充图15

补充图15：GRIA4附近一段潜在灵长类特异性约束区域的核苷酸依赖性图谱（所有模型）。行表示在不同进化时间尺度上训练的模型（V：脊椎动物，M：哺乳动物，P：灵长类），列表示不同模型规模（参数数量）。

）。

这些结果表明，GPN-Star可利用共进化信号学习具有生物学意义的核苷酸依赖关系，且该关系与已知功能依赖高度吻合。这代表了相较于传统保守性评分（如PhyloP和PhastCons）的一项显著进步。

GPN-Star 学习全基因组尺度的进化约束

为更直接评估模型预测与基因组中进化约束之间的关联，我们利用了gnomAD v3.1.2中的等位基因频率数据；该数据库整合了来自76,156名人类个体的全基因组测序样本[33]。人类群体中的等位基因频率是选择约束的重要指示指标：有害性越强的等位基因，因受到净化选择作用，其频率往往越低[6, 29, 33]。

在此评估中，我们重点将GPN-Star与同样基于多物种比对（WGA）数据学习进化约束的PhyloP和PhastCons进行比较。为评估各模型对等位基因频率与约束关系的刻画能力，我们分别采用各模型的脊椎动物、哺乳动物及灵长类版本，并获取所有gnomAD v3染色体22上的变异位点（该染色体为GPN-Star训练时预留的未使用染色体）的预测值。随后，我们比较了不同模型所划分的若干分位数区间内的平均等位基因频率。如图5A所示，在全部三个进化时间尺度下，GPN-Star较低分位数区间内的变异位点，其平均等位基因频率始终低于对应PhyloP和PhastCons区间，表明GPN-Star能更准确地刻画人类基因组中的选择约束。

接下来，我们开展了更定量的评估，重点关注模型得分分布中最具有有害性的尾部区域——这一区域在许多人遗传学应用中尤为重要。我们量化了各模型所预测的最受限尾部中罕见变异相对于常见变异的富集程度。此处，我们将罕见变异定义为单例变异（singletons），而将常见变异定义为等位基因频率为( > 5% )的变异。由于罕见变异平均而言更具危害性，因此约束预测更准确的模型应在最受限尾部中展现出更高的单例变异富集度。如图5B所示，全部三种GPN-Star模型在罕见变异富集度上均显著高于PhyloP、PhastCons或CADD。在GPN-Star模型中，脊椎动物模型整体富集度最强，并优于同样基于脊椎动物基因组训练的GPN-MSA。当按分子效应对变异进行分层时，GPN-Star在每一类中均实现了最高富集度（图5C）。值得注意的是，GPN-Star（V）在错义变异上表现最佳，而GPN-Star（M）则在同义变异和非编码变异类别中领先，这与既往基准测试中观察到的趋势一致。

控制突变率变异可提升变异效应预测性能

我们研究了依赖序列背景的突变率变异对模型预测的影响。由于GPN-Star是在自然界中观测到的基因组序列上训练的，其预测结果天然同时反映了突变与选择过程。正如预期，原始模型得分与Roulette[71]估算的突变率呈中等程度相关性（在22号染色体上，各模型的Spearman’s ( \rho = {0.31} - {0.34} )）。为分离选择效应的影响，我们设计了一种简单而有效的校准流程，以从模型得分中去除突变率变异的影响（方法部分）。该校准流程将模型得分与突变率估算值的相关性降至可忽略水平，低于PhyloP和PhastCons在相同进化时间尺度下的对应值，并在大多数下游基准测试中提升了性能（图5D，补充图16

补充图16：突变率校准。（A）GPN-Star预测值与Roulette突变率估计值( \left( {x\text{ axis }}\right) )及Gnocchi约束得分( \left( {y\text{ axis }}\right) )在校准前后的Spearman相关性。（B）–（I）GPN-Star预测值在校准前后于变异效应预测基准上的性能。

依据CC-BY-NC-ND 4.0国际许可协议提供。

）。

我们进一步通过将GPN-Star得分与Gnocchi（一种基于gnomAD v3全基因组估算、并显式纳入突变率变异影响的约束得分）进行比较，验证了该校准流程的有效性。经校准的GPN-Star得分与Gnocchi的一致性高于PhyloP或PhastCons（图5D）。

因此，本文呈现的所有分析均采用校准后的得分。作为最终合理性检验，我们确认Roulette突变率本身在其低突变率尾部对单例变异几乎无富集（图5B），且在致病性预测任务中表现接近基线水平（补充图2）。这证实了我们的罕见变异富集度与致病性基准测试主要评估的是选择约束，而非突变偏差。

这些发现凸显了在未来通用语言模型（gLMs）开发中处理突变率变异的关键重要性，尤其对于涉及功能约束预测的应用场景。

GPN-Star是一种通用且灵活的框架，可用于跨多样物种及比对数据学习进化约束

旨在理解遗传变异的重大实验工作主要集中于人类基因组[72, 73]。相比之下，研究人员在其他物种中大规模探究变异功能影响的手段十分有限。我们提出GPN-Star作为一种强大且通用的框架，用于解析多种物种中的遗传变异。GPN-Star仅需一个多物种全基因组比对（WGA）即可获得全基因组变异效应预测；而此类数据正日益广泛地覆盖众多物种[26–28]。

作为示范，我们在五种重要模式生物上训练了GPN-Star模型：小鼠（M. musculus）、鸡（G. gallus）、黑腹果蝇（D. melanogaster）、秀丽隐杆线虫（C. elegans）和拟南芥（A. thaliana）。在每种情况下，我们都证明所学习到的约束信息对解读该物种的遗传变异极具价值。针对这五个物种，我们分别收集了包含18至135个物种基因组比对的WGA数据集，并采用与人类模型相似的训练流程，仅作少量调整（补充表1）。

由于非人类物种中缺乏大规模、经人工审编且带有功能注释的变异数据集，我们采用群体遗传学数据进行评估。具体而言，我们收集了五个物种各自的群体水平变异数据（方法部分），并检验预测分布中最具有有害性尾部区域的罕见变异富集程度——该策略与用于人类模型的gnomAD基准分析一致。在全部五个物种中，GPN-Star得分均显著高于PhyloP和PhastCons（图6A），表明其对全基因组进化约束的预测更为准确。在基于分子效应的不同变异类别中，GPN-Star在几乎所有情况下均展现出最高的富集程度（补充图17

补充图17：按分子效应分层的非人源模型评估。（A）群体基因组数据库中罕见变异与常见变异在有害性评分尾部的富集情况（阈值设定为使各评分产生30个假阳性发现），与适配于相同多序列比对的PhyloP和PhastCons比较。（B）MMrdb致病变异与WMGP常见变异在( M ). musculus中的分类。（C）FlyBase致死变异与DEST常见变异在( D ). melanogaster中的分类。（D）( C ). elegans致死变异与CaeNDR常见变异的分类。

）。在拟南芥（A. thaliana）中，GPN-Star与其他模型之间的性能差距较其他物种更小，这可能源于比对规模较小，或源于植物全基因组比对（WGAs）已知质量较低[74]。

对于其中三个物种，我们成功收集到经人工审编的致病性变异，以开展进一步评估。针对( M ) . musculus模型，我们采用了MMRdb数据库中的致病性变异（方法部分）。在区分这些致病性变异与群体中常见变异方面，GPN-Star在全基因组范围（图6B）及各变异类别内（补充图17）均优于PhyloP和PhastCons。类似地，我们收集了黑腹果蝇（D. melanogaster）的实验验证致死变异（来自FlyBase [75]）以及秀丽隐杆线虫（C. elegans）的实验验证致死变异（来自一项实验研究[76]），结果表明GPN-Star在这两套数据集上均持续优于其他模型（图6C、D及补充图17）。

除变异优先级排序外，GPN-Star还可通过核苷酸依赖性分析，探索所学习的功能元件及其共进化关系。例如，我们在黑腹果蝇（D. melanogaster）中广为人知的MSE增强子位点处，展示了GPN-Star所学习到的转录因子结合位点（TFBS）间的核苷酸依赖性（图6E）。这些转录因子的协同作用驱动果蝇胚胎中精确的空间图式形成[77, 78]。这凸显了GPN-Star作为一种高效、低成本工具的价值：可用于跨基因组探究功能元件及其依赖关系，并可辅助和指导实验研究。

讨论

GPN-Star是一种新颖的基因组语言建模框架，其归纳偏置根植于进化生物学，具备强大表达能力。它在全人类基因组范围内预测有害变异时，始终优于现有方法；在精细定位的GWAS变异分类任务中亦表现最优，并提升了罕见变异关联检验中的变异优先级排序效果。在涵盖106种复杂性状的评估中，GPN-Star大幅推进了遗传力富集分析的前沿水平；我们的分析还为灵长类特异性约束及组织特异性的重要性提供了新见解。模型可解释性进一步揭示：GPN-Star无需任何监督信号，即可自主学习识别基因组中广泛的功能元件——从增强子到单个转录因子结合位点（TFBS），乃至其共进化依赖关系。GPN-Star是一个通用且可扩展的框架，仅需输入多序列比对的基因组序列，并消耗相对有限的GPU资源，即可训练出高效模型。以全基因组比对（WGA）为基础的经典模型（如PhastCons [24] 和 PhyloP [25]）二十余年来一直是生物学家与临床医生不可或缺的工具；我们的结果表明，GPN-Star可作为这些工具的有力补充与延伸，在多种物种及不同比对数据上均提供更优性能。

我们结果的核心洞见之一是：所学习的进化约束强度高度依赖于WGA训练数据所涵盖的进化时间尺度。总体而言，所观察到的模式与既有的基因组进化认知一致：编码区变异及更罕见、效应更大的变异，由更深层时间尺度训练所得模型评估得更为准确。这类变异多发生于高度受限的基因组位点，其演化速率缓慢，使得远缘物种的同源序列仍能保留具信息性的基因组语法。一个平行例证见于蛋白质语言模型——其得益于跨越原核生物至人类的极度多样化的训练序列[79]。

相比之下，非编码区变异及更常见、效应更小的变异，则在更浅层的时间尺度下建模效果更佳。这类变异通常出现在演化更快的位点，其序列背景变化更为迅速。例如，某些调控元件已被证实仅在灵长类内部存在约束[56]。综上，我们的发现强调：进化时间尺度是基因组语言模型（gLMs）开发与应用中一项关键考量因素。

GPN-Star 在模型规模和上下文尺寸均显著小于现有单序列gLMs的情况下，实现了最先进的性能。这种高效性很可能源于比对所提供的显式同源信息。值得注意的是，基于比对的蛋白质模型同样以少得多的参数量超越了单序列模型（例如，AlphaFold [10] 相较于 ESMFold [80]，或 MSA Transformer [11] 相较于 ESM-1b [79]）。尽管我们观察到增大模型尺寸带来适度的性能提升，但扩大上下文尺寸带来的改进却相对有限（补充图18

补充图18：模型规模与上下文规模的缩放实验。（A）参数量分别为( {25}\mathrm{M},{85}\mathrm{M} )和( {200}\mathrm{M} )的GPN-Star模型在各基准测试上的性能。（B）参数量为( {25}\mathrm{M} )、上下文长度分别为( {128}\mathrm{{bp}},{256}\mathrm{{bp}} )和( {512}\mathrm{{bp}} )的GPN-Star模型在各基准测试上的性能。所用基准测试与图2相同。除ProteinGym外，所有基准测试的性能指标均为精确率-召回率曲线下面积（AUPRC），误差线表示1000次自助重采样的标准误；对于ProteinGym，点及误差线分别表示各DMS实验中Spearman相关系数的均值与标准差。

bioRxiv预印本 doi: https://doi.org/10.1101/2025.09.21.677619 ；

）。一种合理的解释是：全基因组比对（WGAs）通常以参考物种为基准构建，且常由小而高度碎片化的共线性区块组成；因此，当上下文尺寸增大时，窗口内的核苷酸在非参考物种的实际基因组中可能并不连续，从而向模型引入误导性的上下文信息。如何在WGA数据中有效利用更大上下文尺寸，是未来研究的一个有前景方向。

尽管本研究聚焦于变异解读，gLMs还具有更广泛的应用潜力，包括基因组注释与功能元件发现 [3]，这些仍有待探索。我们的模型可解释性结果表明，GPN-Star在这些任务中具备应用潜力。在此场景下，扩展上下文尺寸可能尤为有用，值得进一步研究。

我们在致病性预测与复杂性状遗传力方面的结果，展现了GPN-Star在人类遗传学中的实用性。初步开展的罕见变异关联检验实验已显示出明显改进，但当前分析仅限于外显子组数据。我们预期其在非编码区将带来更大收益——该区域的变异优先排序历来更具挑战性。GPN-Star的预测结果还可增强其他下游应用，例如功能导向的精细定位、多基因风险评分及表型预测。我们正与领域专家合作，积极拓展这些方向。

除识别进化约束与有害变异外，另一重要未来方向是研究人类特异性适应。当前GPN-Star模型侧重于捕获跨物种的进化信号，缺乏解析人类群体内部选择过程所需的分辨率。我们预计这将是未来工作的一项关键方向，需依赖更高时空分辨率的进化数据进行学习，例如人类群体基因组与古人类基因组。

仅凭无标签DNA序列即可获取如此丰富的信息，实属惊人。功能基因组学与序列到功能模型为遗传变异提供了互补视角。尽管它们在致病性与遗传力分析中表现逊于我们的进化方法，但在研究基因调控机制、组织与细胞特异性以及基因表达等分子表型方面仍不可或缺。在我们的组织特异性遗传力分析中，采用GTEx来源的组织特异性基因构建了组织特异性GPN-Star注释，该注释在各类性状中实现了前所未有的富集效果。未来一项有前景的研究方向是整合进化与功能基因组学数据，开发更强大的多模态gLMs。

跨越生命之树的大规模基因组测序与比对工作正在加速 ( \lbrack {27},{28} ) , ( {81},{82}\rbrack ) 。GPN-Star恰逢其时，可充分受益于这一增长趋势：它能便捷地扩展至新比对数据，并有望从日益多样化的训练数据中进一步提升能力。在基因组数据迅猛扩张的时代，我们期待GPN-Star成为推动遗传变异认知进步的有力工具。

致谢

本研究部分得到美国国立卫生研究院（NIH）资助，项目编号为R35-GM134922和3P40-OD011102-24S1 7772。图1与图6中的动物图标源自Flaticon（ http://www.flaticon.com/ ）。

图1：GPN-Star概览。（A）GPN-Star模型架构示意图。模型输入为全基因组比对窗口；目标序列与源序列均由此窗口构建。源序列经注意力池化压缩为支系层级嵌入，池化过程遵循物种树；目标序列则通过多层GPN-Star编码器模块编码，其中系统发育信息引导的交叉注意力模块依据物种树推断的物种间进化距离，整合源序列信息；最终，分类层将编码后的嵌入转换为目标序列各位点的核苷酸概率，用于计算训练损失或预测变异影响。模型完整描述见“方法”部分。（B）GPN-Star在人类基因组上的应用。分别基于脊椎动物、哺乳动物和灵长类比对数据训练了三个模型，以学习不同进化时间尺度上的功能约束。Mya：百万年前；cLLR：校准对数似然比（参见“方法”）。（C）脊椎动物、哺乳动物及灵长类GPN-Star模型在下游任务中的综合性能。对于ClinVar、COSMIC、OMIM、HGMD及GWAS精细定位数据集，性能指标为精确率-召回率曲线下面积（AUPRC）；对于ProteinGym，指标为各实验中Spearman相关系数的均值( \rho )；对于S-LDSC，指标为遗传力富集度。所有性能指标经线性缩放至0–1区间：其中0定义为三个GPN-Star模型及在相同三个进化时间尺度上拟合的PhyloP与PhastCons中该任务最低指标值，1定义为最高指标值。

图2：GPN-Star在人类全基因组变异效应预测中的性能。（A）ClinVar致病性与良性错义变异的分类。（B）COSMIC高频错义变异（>0.1%）与gnomAD v3常见错义变异的分类。（C）UK Biobank中66个性状的GWAS精细定位结果：潜在因果错义变异（PIP>0.9）与非因果错义变异（PIP<0.01）的分类。（D）GPN-Star与ProteinGym中31个人类蛋白深度突变扫描适应度评分的Spearman相关性箱线图。（E）OMIM致病性与gnomAD v3常见非编码变异的分类。（F）HGMD致病性与gnomAD v3常见非编码变异的分类。（G）UK Biobank中84个性状的GWAS精细定位结果：潜在因果非编码变异（PIP>0.9）与非因果非编码变异（PIP<0.01）的分类。（H）OMIM致病性与gnomAD v3常见（( y )轴）及GWAS精细定位潜在因果与非因果( \left( {x\text{ axis }}\right) )启动子变异的联合分类。图（A）–（C）及（E）–（G）所用性能指标为精确率-召回率曲线下面积（AUPRC）；误差条表示1000次自助重采样的标准误。（I）雷达图对比GPN-Star、PhyloP与PhastCons在脊椎动物、哺乳动物及灵长类三个时间尺度上于七项致病性预测基准任务中的性能。性能指标经线性缩放至0–1区间：其中0定义为全部九个模型中该任务最低指标值，1定义为最高指标值。（J）基于UK Biobank中161,822名欧洲血统无关个体的外显子组测序（WES）数据，对34个定量性状开展罕见变异关联检验：比较原始DeepRVAT与整合三个GPN-Star预测（V、M与P）的DeepRVAT在显著基因数量（( x )轴）及与既往大样本研究重复发现的显著基因数量（( \left( {y\text{ axis }}\right) )）上的表现。

图3：人类复杂性状中的SNP遗传力分析。（A）不同模型对复杂性状遗传力的判别能力。我们将模型得分二值化，以选取约束性最强的前( {0.1}% )个常见变异位点。我们针对每种模型注释分别运行S-LDSC，同时控制96个基线特征，并对106个独立性状的结果进行荟萃分析。遗传力富集度指某注释所解释的遗传力比例除以其自身尺寸。条件效应( \left( {\tau }^{ \star }\right) )衡量该注释对遗传力的独特贡献（即未被现有注释所解释的部分）。误差线表示标准误。（B）将模型注释限定于编码区时的性能。（C）将模型注释限定于非编码区时的性能。（D）随着高约束SNP比例变化的性能比较。（E）GPN-Star（P）与（M）在27个性状中富集度差异随估计有效多基因性变化的趋势。p值为单侧检验。虚线为普通最小二乘拟合线。（F）GPN-Star（P）优先识别的常见变异类型。变异依据Ensembl功能后果与ENCODE SCREEN候选顺式调控元件联合注释。仅显示占比高于( 1% )的类型。比值比参照常见变异中底部( {99.9}% )部分。（G）GPN-Star（P）与Enformer在组织无关与组织特异性适配版本间，于分组性状上的性能比较。（H）GPN-Star（P）与Enformer在组织无关与组织特异性适配版本间，于单个性状上的性能比较。

图4：GPN-Star可解释性揭示功能元件及其依赖关系。（A）不同基因组窗口嵌入向量的可视化，按其注释区域着色。（B）不同基因组窗口嵌入向量的可视化，按保守性着色；保守性定义为窗口内PhastCons（P）值的( {75}^{\text{ th }} )百分位数。（C）( {TH} )启动子中的核苷酸依赖性图谱。热图中位置( i,j )处的强度对应于位置( i )核苷酸对位置( j )预测核苷酸概率的影响程度。（D）LDLR启动子中的核苷酸依赖性图谱。MPRA效应指平滑后的绝对log倍数变化。

bioRxiv预印本 doi: https://doi.org/10.1101/2025.09.21.677619 ；

图5：GPN-Star评分反映人类基因组的进化约束。（A）基于gnomAD染色体22上双等位基因位点，按GPN-Star、PhyloP及PhastCons（分别基于脊椎动物、哺乳动物和灵长类比对）划分的十分位数区间( \left( {\left\lbrack {0,{10}^{4}}\right) ,\left( {{10}^{4},{10}{3}}\right\rbrack ,\ldots ,\left( {{10}^{1},1}\right\rbrack }\right) )的平均MAF。（B和C）罕见变异（单例）与常见变异（MAF>5%）在有害性评分尾部的富集情况（阈值设定为各评分产生30个假阳性发现）。（B）显示全基因组富集结果，（C）按变异分子后果分层显示富集结果。罕见变异经下采样，使其在各类别中数量与常见变异一致。（D）GPN-Star与PhyloP、PhastCons在三个进化时间尺度上，分别与染色体( {22}\left( {x\text{ axis }}\right) )上Roulette突变率估计值及Gnocchi约束估计值（( y )轴）相关性的性能比较。

图6：GPN-Star 在非人类物种中的应用。(A) 在五个非人类物种中评估 GPN-Star 模型对罕见变异与常见变异富集能力（方法中定义的阈值），所用数据为对应物种群体基因组数据库中致病性评分尾部的变异（阈值设定为使各评分产生30个假阳性发现），并与在同一多序列比对上拟合的 PhyloP 和 PhastCons 进行比较。(B) 在 ( M ) . musculus 中对 MMrdb 致病变异与 WMGP 常见变异的分类。(C) 在 D. melanogaster 中对 FlyBase 致死变异与 DEST 常见变异的分类。(D) 在 C. elegans 中对线虫致死变异与 CaeNDR 常见变异的分类。(B)-(D) 中采用的性能指标为精确率-召回率曲线下面积（AUPRC）。误差线表示1000次自助重采样的标准误。(E) ( D ) . melanogaster 模型在 MSE 增强子内某位点的核苷酸依赖性图谱，该位点已标注已知转录因子结合位点（TFBS）。

参考文献

1. Jukes, T. H. & Cantor, C. 载于《哺乳动物蛋白质代谢》21–132页（Academic Press，纽约，1969年）。

2. Dayhoff, M. O., Eck, R. V. & Park, C. M. 载于《蛋白质序列与结构图谱》（主编：Dayhoff, M. O.）89–99页（国家生物医学研究基金会，华盛顿特区，1972年）。

3. Benegas, G., Ye, C., Albors, C., Li, J. C. & Song, Y. S. 基因组语言模型：机遇与挑战。《遗传学趋势》41, 286–302 (2025)。

4. Achiam, J. 等。GPT-4 技术报告。arXiv 预印本 arXiv:2303.08774。https://arxiv.org/abs/2303.08774 (2023)。

5. Vaswani, A. 等。“注意力机制即全部”载于《神经信息处理系统进展》（编者：Guyon, I. 等）第30卷（Curran Associates, Inc., 2017）。

6. Benegas, G., Batra, S. S. & Song, Y. S. DNA语言模型是全基因组变异效应的强有力预测器。《美国国家科学院院刊》120, e2311219120 (2023)。

7. Benegas, G., Eraslan, G. & Song, Y. S. 人类遗传学中因果调控变异预测的DNA序列模型基准测试。bioRxiv 预印本。https://doi.org/10.1101/2025.02.11.637758 (2025)。

8. Dalla-Torre, H. 等。核苷酸Transformer：构建并评估面向人类基因组的鲁棒基础模型。《自然·方法》22, 287–297 (2024)。

9. Brixi, G. 等。Evo 2：跨生命所有领域的基因组建模与设计。bioRxiv 预印本。https://doi.org/10.1101/2025.02.18.638918 (2025)。

10. Jumper, J. 等。基于AlphaFold的高度准确蛋白质结构预测。《自然》596, 583–589 (2021)。

11. Rao, R. M. 等。“MSA Transformer”载于《国际机器学习会议论文集》(2021)，8844–8856页。

12. Frazer, J. 等。利用深度生成模型基于进化数据预测疾病变异。《自然》599, 91–95 (2021)。

13. Truong Jr, T. & Bepler, T. PoET：一种将蛋白家族建模为“序列之序列”的生成模型，载于《神经信息处理系统进展》（编者：Oh, A. 等）第36卷（Curran Associates, Inc., 2023），77379–77415页。

14. Ma, C. 等。“用于蛋白质表征学习的检索式序列增强”，载于《2024年自然语言处理经验方法会议论文集》（编者：Al-Onaizan, Y., Bansal, M. & Chen, Y.-N.）（计算语言学协会，美国佛罗里达州迈阿密，2024年11月），1738–1767页。 https://aclanthology.org/2024.emnlp-main.104/ 。

15. Alamdari, S. 等. 利用进化扩散生成蛋白质：序列即全部所需。bioRxiv，2023-09（2023）。

16. Li, P., Cheng, X., Song, L. & Xing, E. 面向蛋白质结构预测的检索增强型蛋白质语言模型。bioRxiv，2024-12（2024）。

17. Sgarbossa, D. & Bitbol, A.-F. RAG-ESM：通过序列检索提升预训练蛋白质语言模型性能。PRX Life 3, 033013（2025）。

18. Yang, K. K. 等. Dayhoff 图谱：扩展序列多样性以提升蛋白质生成效果。bioRxiv，2025-07（2025）。

19. Akiyama, Y., Zhang, Z., Mirdita, M., Steinegger, M. & Ovchinnikov, S. 借助MSA Pairformer降低蛋白质语言建模规模。bioRxiv，2025-08（2025）。

20. Truong Jr, T. F. & Bepler, T. 利用多模态检索增强基础模型理解蛋白质功能。arXiv 预印本 arXiv:2508.04724（2025）。

21. Blanchette, M. 等. 使用线程块集比对器（Threaded Blockset Aligner）比对多条基因组序列。Genome Research 14, 708–715（2004）。

22. Armstrong, J. 等. Progressive Cactus：面向“千基因组时代”的多基因组比对工具。Nature 587, 246–251（2020）。

23. Davydov, E. V. 等. 利用GERP++识别大量受选择约束的人类基因组区域。PLoS 计算生物学 6, e1001025（2010）。

24. Siepel, A. 等. 脊椎动物、昆虫、线虫及酵母基因组中的进化保守元件。Genome Research 15, 1034–1050（2005）。

25. Pollard, K. S., Hubisz, M. J., Rosenbloom, K. R. & Siepel, A. 在哺乳动物系统发育树中检测非中性替换率。Genome Research 20, 110–121（2010）。

26. Christmas, M. J. 等. 跨数百种胎盘哺乳动物的进化约束与创新。Science 380, eabn3943（2023）。

27. Rhie, A. 等. 迈向所有脊椎动物物种完整无误的基因组组装。Nature 592, 737–746（2021）。

28. Stiller, J. 等. 基于科级基因组揭示鸟类演化的复杂性。Nature 629, 851–860（2024）。

29. Benegas, G., Albors, C., Aw, A. J., Ye, C. & Song, Y. S. 一种基于多物种比对的DNA语言模型可预测全基因组变异效应。Nature Biotechnology。https://doi.org/10.1038/s41587-024-02511-w（2025）。

30. Gazal, S. S-LDSC 参考文件（Zenodo，2024年1月）。 https://doi.org/10.5281/zenodo.10515792 。

31. Landrum, M. J. 等. ClinVar：序列变异与人类表型关系的公共档案. 《核酸研究》42卷，D980–D985页（2014年）。

32. Tate, J. G. 等. COSMIC：癌症体细胞突变目录. 《核酸研究》47卷，D941–D947页（2019年）。

33. Karczewski, K. J. 等. 基于141,456名人类个体变异量化得出的突变约束谱. 《自然》581卷，434–443页（2020年）。

34. Notin, P. 等. ProteinGym：用于蛋白质适应性预测与设计的大规模基准数据集，《神经信息处理系统进展》（Oh, A. 等编）第36卷（Curran Associates, Inc., 2023年），64331–64379页。 https://proceedings.neurips.cc/paper_files/paper/2023/file/cac723e5ff29f65e3fcbb0739ae91bee-Paper-Datasets_and_Benchmarks.pdf 。

35. Cheng, J. 等. 利用AlphaMissense实现全蛋白质组范围错义变异效应的精准预测. 《科学》381卷，eadg7492（2023年）。

36. Gao, H. 等. 人类与灵长类中可耐受遗传变异的图谱. 《科学》380卷，eabn8153（2023年）。

37. Ghosh, R. 等. 关于良性独立ACMG/AMP判别标准的更新建议. 《人类突变》39卷，1525–1530页（2018年）。

38. Avsec, Z. 等. 通过整合长程互作实现基于序列的有效基因表达预测. 《自然·方法》18卷，1196–1203页（2021年）。

39. Linder, J., Srivastava, D., Yuan, H., Agarwal, V. & Kelley, D. R. 以DNA序列为输入预测RNA-seq覆盖度：一种统一的基因调控模型. 《自然·遗传学》. https://doi.org/10.1038/s41588-024-02053-6（2025年）。

40. Avsec, Z. 等. AlphaGenome：以统一DNA序列模型推动调控变异效应预测. bioRxiv, 2025-06（2025年）。

41. Amberger, J. S., Bocchini, C. A., Schiettecatte, F., Scott, A. F. & Hamosh, A. OMIM.org：在线孟德尔遗传病数据库（OMIM®），一个人类基因与遗传疾病在线目录. 《核酸研究》43卷，D789–D798页（2015年）。

42. Stenson, P. D. 等. 人类基因突变数据库（HGMD®）：优化其在临床诊断或研究场景中的应用. 《人类遗传学》139卷，1197–1207页（2020年）。

43. Avsec, Z. 等. AlphaGenome：以统一DNA序列模型推动调控变异效应预测. bioRxiv. 预印本：https://www.biorxiv.org/content/early/2025/07/11/2025.06.25.661532.full.pdf. https://www.biorxiv.org/content/early/2025/07/11/2025.06.25.661532（2025年）。

44. Jaganathan, K. 等. 利用深度学习预测影响基因表达的启动子突变. 《科学》，eads7373（2025年）。

45. Tomaz da Silva, P. 等. DNA语言模型的核苷酸依赖性分析揭示基因组功能元件. bioRxiv预印本，2024-07. https://www.biorxiv.org/content/10.1101/2024.07.27.605418v1（2024年）。

46. Kanai, M. 等. 来自跨多元人群复杂性状精细定位的洞见. medRxiv，2021-09（2021年）。

47. Bomba, L., Walter, K. & Soranzo, N. 罕见及低频遗传变异在常见疾病中的影响。《基因组生物学》18, 77 (2017)。

48. Lee, S., Abecasis, G. R., Boehnke, M. & Lin, X. 罕见变异关联分析：研究设计与统计检验。《美国人类遗传学杂志》95, 5–23 (2014)。

49. Clarke, B. 等。利用深度集合网络整合变异注释可提升罕见变异关联检验效能。《自然·遗传学》，1–10 (2024)。

50. Backman, J. D. 等。对英国生物银行454,787名参与者的外显子组测序与分析。《自然》599, 628–634 (2021)。

51. Karczewski, K. J. 等。在394,841份英国生物银行外显子组数据中，系统开展单变异及基因水平的数千种表型关联检验。《细胞·基因组学》2 (2022)。

52. Finucane, H. K. 等。利用全基因组关联汇总统计，按功能注释划分遗传力。《自然·遗传学》47, 1228–1235 (2015)。

53. Weissbrod, O. 等。基于功能信息的精细定位与复杂性状遗传力的多基因定位。《自然·遗传学》52, 1355–1363 (2020)。

54. Márquez-Luna, C. 等。纳入功能先验信息可提升英国生物银行与23andMe数据集中多基因预测的准确性。《自然·通讯》12, 6052 (2021)。

55. Sullivan, P. F. 等。利用碱基对尺度的哺乳动物约束信息解析遗传变异与人类疾病的关系。《科学》380, eabn2937 (2023)。

56. Kuderna, L. F. 等。在239种灵长类基因组中鉴定受限序列元件。《自然》625, 735–742 (2024)。

57. O'Connor, L. J. & Sella, G. 衡量与估计性状多基因性的原理性方法。bioRxiv, 2025-07 (2025)。

58. Dyer, S. C. 等。Ensembl 2025。《核酸研究》53, D948–D957 (2025)。

59. Moore, J. E. 等。扩展版顺式调控元件候选者登记库，用于转录调控研究。bioRxiv (2024)。

60. Karollus, A., Mauermeier, T. & Gagneur, J. 当前基于序列的模型虽能捕捉启动子中的基因表达决定因子，却大多忽略远端增强子。《基因组生物学》24, 56 (2023)。

61. Fabiha, T. 等。一种共识性变异-功能评分，用于按功能优先排序致病相关变异。bioRxiv, 2024-11 (2024)。

62. Finucane, H. K. 等。特异性表达基因的遗传力富集可识别疾病相关组织与细胞类型。《自然·遗传学》50, 621–629 (2018)。

63. 张, Z. 等. 蛋白质语言模型学习相互作用序列基序的进化统计特征。《美国国家科学院院刊》121, e2406285121 (2024)。

64. Verbeek, M. M. 等. 酪氨酸羟化酶基因环磷酸腺苷反应元件中的突变。《神经病学年鉴》62, 422–426 (2007)。

65. 董, H.-Y., 冯, J.-Y., 岳, X.-J., 单, L. & 贾, F.-Y. 酪氨酸羟化酶缺乏所致多巴反应性肌张力障碍：三例报告及文献综述。《医学》99, e21753 (2020)。

66. Ribasés, M. 等. 一例多巴反应性脑病患者中酪氨酸羟化酶基因启动子区纯合突变：临床、生化及遗传学分析。《分子遗传学与代谢》92, 274–277 (2007)。

67. Stamelou, M. 等. 酪氨酸羟化酶缺乏所致肌阵挛-肌张力障碍综合征。《神经病学》79, 435–441 (2012)。

68. Kircher, M. 等. 在单碱基对分辨率下对二十种疾病相关调控元件进行饱和诱变。《自然·通讯》10 (2019)。

69. Khamis, A. 等. 家族性高胆固醇血症患者LDLR 5′UTR及启动子区四种变异的功能分析。《欧洲人类遗传学杂志》23, 790–795 (2015)。

70. Bennett, M. K., Ngo, T. T., Athanikar, J. N., Rosenfeld, J. M. & Osborne, T. F. 固醇调节元件结合蛋白与Sp1对低密度脂蛋白受体基因启动子的协同激活作用被阴阳1蛋白特异性破坏。《生物化学杂志》274, 13025–13032 (1999)。

71. Seplyarskiy, V. 等. 一种单碱基对分辨率突变率模型识别出DNA聚合酶III转录的诱变效应。《自然·遗传学》，1–8 (2023)。

72. ENCODE联盟等. 人类基因组DNA功能元件综合百科全书。《自然》489, 57–74 (2012)。

73. Lonsdale, J. 等. 基因型-组织表达（GTEx）计划。《自然·遗传学》45, 580–585 (2013)。

74. 宋, B., Buckler, E. S. & Stitzer, M. C. 挖掘植物多样性亟需新型全基因组比对工具。《植物科学趋势》29, 355–369 (2024)。

75. Öztürk-Çolak, A. 等. FlyBase：果蝇基因与基因组数据库更新。《遗传学》227, iyad211 (2024)。

76. 秦, Z. 等. 秀丽隐杆线虫必需基因的基因组鉴定与功能表征。《G3：基因、基因组、遗传学》8, 981–997 (2018)。

77. Small, S., Blair, A. & Levine, M. 果蝇胚胎中even-skipped第2条带的调控。《EMBO杂志》11, 4047–4057 (1992)。

78. Bothma, J. P. 等. eve第2条带表达的动态调控揭示活体果蝇胚胎中的转录爆发现象。《美国国家科学院院刊》111, 10598–10603 (2014)。

79. Rives, A. 等。生物学结构与功能源于将无监督学习扩展至2.5亿条蛋白质序列。《美国国家科学院院刊》118卷，e2016239118（2021年）。

80. Lin, Z. 等。基于语言模型的进化尺度原子级蛋白质结构预测。《科学》379卷，1123–1130页（2023年）。

81. 生命之树项目联盟，D. T. 局部测序，全局思考：达尔文生命之树项目。《美国国家科学院院刊》119卷，e2115642118（2022年）。

82. Lewin, H. A. 等。地球生物基因组计划2020：启动时钟2022。

83. Miles, A. 等。zarr-developers/zarr-python：v3.0.7 版本 v3.0.7。2025年4月。 https://doi.org/10.5281/zenodo.15255977 。

84. Kent, W. J. 等。UCSC人类基因组浏览器。《基因组研究》12卷，996–1006页（2002年）。

85. 用于科学发现与保护的比较基因组学多用途工具。《自然》587卷，240–245页（2020年）。

86. Albors, C., Li, J. C., Benegas, G., Ye, C. & Song, Y. S. 基于系统发育的基因组语言建模方法，载于《计算分子生物学研究》（主编：Sankararaman, S.）（Springer Nature Switzerland, Cham, 2025），99–117页。ISBN: 978-3-031-90252-9。

87. Tian, F., Yang, D.-C., Meng, Y.-Q., Jin, J. & Gao, G. PlantRegMap：绘制植物功能调控图谱。《核酸研究》48卷，D1104–D1113页（2020）。

88. Hubisz, M. J., Pollard, K. S. & Siepel, A. PHAST与RPHAST：基于时空模型的系统发育分析。《生物信息学简报》12卷，41–51页（2011）。

89. Su, J. 等。Roformer：采用旋转位置嵌入增强的Transformer。arXiv预印本 arXiv:2104.09864（2021）。

90. Li, S. 等。相对位置的功能插值提升长上下文Transformer性能，载于《第十二届国际表征学习会议》（2024）。 https://openreview.net/forum?id=rR03qFesqk 。

91. Aggarwala, V. & Voight, B. F. 一种扩展的序列上下文模型可广泛解释人类基因组中多态性水平的变异。《自然·遗传学》48卷，349–355页（2016）。

92. Tarailo-Graovac, M. & Chen, N. 利用RepeatMasker识别基因组序列中的重复元件。《当代生物信息学方案》25卷，4.10.1–4.10.14页。电子版： https://currentprotocols.onlinelibrary.wiley.com/doi/pdf/10.1002/0471250953.bi0410s25 。https://currentprotocols.onlinelibrary.wiley.com/doi/abs/10.1002/0471250953.bi0410s25（2009）。

93. Open2C 等。Bioframe：在Pandas数据框中对基因组区间执行操作。《生物信息学》40卷，btae088（2024）。

94. McInnes, L., Healy, J. & Melville, J. UMAP：用于降维的一致流形逼近与投影。arXiv预印本 arXiv:1802.03426（2018）。

95. Schubach, M., Maass, T., Nazaretyan, L., Röner, S. & Kircher, M. CADD v1.7：融合蛋白质语言模型、调控CNN及其他核苷酸层级评分，提升全基因组变异预测能力。《核酸研究》52卷，D1143–D1154页（2024）。

96. Chen, K. M., Wong, A. K., Troyanskaya, O. G. & Zhou, J. 基于序列的全基因组调控活性全局图谱，助力解析人类遗传学。《自然·遗传学》54卷，940–949页（2022）。

97. Bycroft, C. 等。具备深度表型与基因组数据的英国生物样本库资源。《自然》562卷，203–209页（2018）。

98. Chen, S. 等。基于76,156个人类基因组变异构建的全基因组突变约束图谱。bioRxiv，2022-03（2022）。

99. Brandes, N., Goldman, G., Wang, C. H., Ye, C. J. & Ntranos, V. 利用深度蛋白质语言模型进行全基因组疾病变异效应预测。《自然·遗传学》。ISSN: 1546-1718。https://doi.org/10.1038/s41588-023-01465-0（2023年8月）。

100. Halligan, D. L. 等。编码蛋白与调控区域变异对鼠科啮齿类适应性分子进化的贡献。《PLoS遗传学》9卷，e1003995（2013）。

101. Davies, R. W. 影响野生小鼠遗传变异的因素。博士论文（牛津大学，2015）。

102. Harr, B. 等. 家鼠（Mus musculus）及其近缘种斯氏小鼠（Mus spretus）野生种群的基因组资源。Scientific data 3, 1–14 (2016)。

103. Phifer-Rixey, M. 等. 家鼠环境适应的基因组基础。PLoS Genetics 14, e1007672 (2018)。

104. Payseur, B. A. & Jing, P. 戈夫岛巨型小鼠正向选择的基因组靶点。Molecular Biology and Evolution 38, 911–926 (2021)。

105. Fujiwara, K. 等. 全基因组分析揭示欧亚大陆家鼠（Mus musculus）种群结构。Genome biology and evolution 14, evac068 (2022)。

106. Morgan, A. P. 等. 不同地理尺度下野生家鼠（Mus musculus）的种群结构与近交程度。Heredity 129, 183–194 (2022)。

107. Lawal, R. A. 等. 基于全基因组测序对两种野生家鼠亚种的分类学评估。Scientific reports 12, 20866 (2022)。

108. Chen, S., Zhou, Y., Chen, Y. & Gu, J. fastp：一款超快速一体化FASTQ预处理器。Bioinformatics 34, i884–i890 (2018)。

109. Hickey, G. 等. 利用Minigraph-Cactus基于基因组比对构建泛基因组图谱。Nature biotechnology 42, 663–673 (2024)。

110. Ferraj, A. 等. 解析多种小鼠基因组中的结构变异，揭示转座元件引发的染色质重塑。Cell Genomics 3 (2023)。

111. Sirén, J. 等. 泛基因组学实现对5202个高度多样化基因组中已知结构变异的基因分型。Science 374, abg8871 (2021)。

112. Garrison, E. 等. 变异图谱工具包通过在参考基因组中表征遗传变异提升序列比对效果。Nature biotechnology 36, 875–879 (2018)。

113. Poplin, R. 等. 基于深度神经网络的通用SNP及小片段插入/缺失变异检出工具。Nature biotechnology 36, 983–987 (2018)。

114. Lin, M. F. 等. GLnexus：面向大规模队列测序的联合变异检出工具。BioRxiv, 343970 (2018)。

115. Fairfield, H. 等. 利用全外显子组测序在小鼠中发现突变。Genome biology 12, 1–12 (2011)。

116. Fairfield, H. 等. 外显子组测序在91株孟德尔遗传病小鼠品系中鉴定致病突变。Genome research 25, 948–957 (2015)。

117. Perez, G. 等. UCSC基因组浏览器数据库：2025年更新。Nucleic Acids Research 53, D1243–D1249 (2025)。

118. Nunez, J. C. 等. 利用扩展版 DEST 2.0 基因组资源揭示黑腹果蝇结构化群体中的全球适应痕迹。bioRxiv，2024-11（2024）。

119. Crombie, T. A. 等. CaeNDR：秀丽隐杆线虫自然多样性资源。《核酸研究》52 卷，D850–D858（2024）。

120. Fu, W. 等. Galbase：整合鸡多组学数据的综合数据库。《BMC 基因组学》23 卷，364（2022）。

121. Alonso-Blanco, C. 等. 1,135 个基因组揭示拟南芥全球多态性模式。《细胞》166 卷，481–491（2016）。

122. Gazal, S. 等. 连锁不平衡依赖的人类复杂性状遗传架构揭示负向选择作用。《自然·遗传学》49 卷，1421–1427（2017）。

123. 千人基因组计划联盟等. 人类遗传变异的全球参考图谱。《自然》526 卷，68–74（2015）。

方法

训练数据

通用工作流程。当前训练工作流程需满足以下条件：

• 多物种全基因组比对（例如，使用 MULTIZ [21] 或 Cactus [22] 获得）。

• 系统发育树。

• PhastCons [24] 和 PhyloP [25] 保守性评分。

MULTIZ 比对（MAF 格式）首先通过 maf2fasta [21] 转换为 FASTA 格式。与 GPN-MSA [29] 一致，我们剔除目标序列中含空位的列。最后，将比对结果转换为 Zarr 格式 [83]，以支持对基因组窗口的快速多线程随机访问。

与 Benegas 等 [29] 类似，我们进一步依据保守性从目标基因组中筛选训练区域。对每个模型，以模型上下文长度为滑动窗口，在未屏蔽的基因组区域上以半窗长为步长移动；定义窗口保守性为该窗口内 PhastCons 评分的第 75 百分位数；保留保守性最高的前 5% 窗口，以及剩余窗口中随机 ( {0.1}% ) 作为模型训练数据。窗口按染色体划分为训练集与验证集，详见补充表 1。对所选窗口，正链与反向互补链序列均纳入训练数据。

数据来源。MULTIZ 比对及关联数据在 UCSC 基因组浏览器 [84] 上可直接获取，涵盖本文涉及的多个物种（

补充表8：源自 UCSC 基因组浏览器的比对结果。

	参考基因组	比对
智人（H. sapiens）	hg38	multiz100way（100种脊椎动物）
小家鼠（M. musculus）	mm39	multiz35way（35种脊椎动物）
红原鸡（G. gallus）	galGal6	multiz77way（77种脊椎动物）
黑腹果蝇（D. melanogaster）	dm6	multiz124way（124种昆虫）
秀丽隐杆线虫（C. elegans）	ce11	multiz135way（135种线虫）

补充表 8）。人类与灵长类及哺乳类的比对 [56, 85] 采用 Cactus [22] 完成，并由 Albors 等 [86] 进行处理。拟南芥（TAIR10 参考基因组）与另外 17 种芸薹目物种的 MULTIZ 比对及保守性评分，自 PlantRegMap [87] 下载。

芸薹目比对的额外系统发育树构建步骤。由于芸薹目物种尚无公开可用的系统发育树，我们依据原始文献 [87] 描述的方法，利用 PHAST 软件包 [88] 基于四倍简并位点构建该树。我们采用 Ensembl [58] 第 60 版注释。由于流水线中 msa_view 工具每次仅支持单条染色体，我们任意选取第 1 号染色体进行处理。

---

msa_view ./maf/tair10_multiz18way/1.maf --4d --features chr1.gff > 4d-codons.ss

msa_view 4d-codons.ss --in-format SS --out-format SS --tuple-size 1 > 4d-sites.ss

phyloFit --tree "(Carica_papaya, (Tarenaya_hassleriana, (Aethionema_arabicum,

    ((Arabis_alpina,(Eutrema_salsugineum,(Thellungiella_parvula,(Sisymbrium_irio,

        (((Brassica_oleracea,Brassica_napus),Brassica_rapa),Raphanus_sativus)))),

        ((Boechera_stricta, (Camelina_sativa, (Capsella_grandiflora, Capsella_rubella)))),

    (Arabidopsis_thaliana, (Arabidopsis_halleri, Arabidopsis_lyrata)))))) "

--msa-format SS --out-root nonconserved-4d --EM --precision MED 4d-sites.ss

---

模型架构

GPN-Star 模型采用一种新颖的基于 Transformer 的架构，作用于多物种全基因组比对的分块数据。我们将一个输入比对块记为 ( \mathbf{X} \in {\mathbb{R}}^{{N}{\mathcal{S}} \times L} )，其中 ( {N}{\mathcal{S}} ) 表示物种数，( L ) 表示该比对块的长度。每个条目（即 token）记为 ( {x}_{j}^{i} )，以整数值表示，编码第 ( i ) 个物种在位置 ( j ) 处的核苷酸身份。

目标序列与源序列。GPN-Star 的核心目标是学习一种核苷酸的条件分布，以建模比对中不同物种序列间的进化过程。具体而言，对于待预测分布的某序列中的每个 token（称为目标序列），GPN-Star 学习：

[\mathbb{P}\left( {{x}{j}^{\left( i\right) } \mid {x}^{\left( i\right) },{\left{ {\mathbf{x}}^{\left( {i}^{\prime }\right) }\right} }{{i}^{\prime } \in \mathcal{S}},{\left{ {\phi }^{\prime }}\right} }_{{i}^{\prime } \in \mathcal{S}}}\right) ,]

其中 ( {x}{j}^{\left( i\right) } ) 是物种 ( i,{x}^{\left( i\right) } ) 的目标序列中位置 ( j ) 处的 token；( i,{\left{ {\mathbf{x}}^{\left( {i}^{\prime }\right) }\right} }{{i}^{\prime } \in \mathcal{S}} ) 是该序列中其余 tokens；( \mathcal{S} ) 是提供进化背景的一组物种的序列（称为源序列）；( {\left{ {\phi }^{\prime }}\right} }_{{i}^{\prime } \in \mathcal{S}} ) 是目标序列与各源序列之间的进化距离。

因此，GPN-Star 的第一步是从每个比对块中分别整理出两组不同的序列：

• 目标序列 ( {\left{ {\mathbf{x}}^{\left( t\right) }\right} }{t \in \mathcal{T}} )，其中 ( {\mathbf{x}}^{\left( t\right) } \in {\mathbb{R}}^{L},\left| \mathcal{T}\right| = {N}{\mathcal{T}} )：用于损失计算与推断的序列。在训练与推断阶段，从中选取 ( {N}{\mathcal{T}} ) 个物种的子集。首个目标物种被指定为主要目标物种，记为 ( {t} )，全基因组比对依其拼接，从而提供连续的基因组坐标以锚定位置信息；它也是使用训练后模型进行推断的主要物种。

• 源序列 ( {\left{ {\mathbf{x}}^{\left( s\right) }\right} }{s \in \mathcal{S}} )，其中 ( {\mathbf{x}}^{\left( s\right) } \in {\mathbb{R}}^{L},\left| \mathcal{S}\right| = {N}{\mathcal{S}} )：比对中所有物种的序列。它为目标序列的分布推断提供进化背景。

将物种归类为分支支系。为高效表示物种间的系统发育关系，我们基于成对系统发育距离将它们划分为离散的分支支系。具体而言，物种被划分为若干支系，使得属于不同支系的任意两个物种之间的成对系统发育距离大于预定义的、模型特定的阈值( {\phi }{\text{ clade }} )，这实质上即单链接聚类。该准则既保证了支系内部的系统发育一致性，又确保了不同支系间物种在进化上足够疏远。支系数量记为( {N}{\mathcal{C}} )。

目标序列与源序列的嵌入。对于目标序列，每个输入词元直接嵌入为( H )维向量：

[{\mathbf{e}}{j}^{\left( t\right) } = {f}\left( {x}_{j}^{\left( t\right) }\right) ,]

其中( {f}{e} : {\mathbb{R}}^{V} \rightarrow {\mathbb{R}}^{H} )为可学习的嵌入函数，( V )为核苷酸词表大小；( {x}^{\left( t\right) } )为目标序列( t )在位置( j )处的独热编码核苷酸。由此得到目标序列嵌入( {\left{ {\mathbf{E}}^{\left( t\right) }\right} }_{t \in \mathcal{T}} )，其中( {\mathbf{E}}^{\left( t\right) } \in {\mathbb{R}}^{L \times H} )。

对于源序列的嵌入，每个支系内各物种的序列通过注意力池化模块，在每个比对列上压缩为单一代表性嵌入。对每个支系( c )，各成员物种在列( i )处的词元在每个注意力头中按如下方式组合：

[\text{ Attn_pool }\left( {{\mathbf{x}}{j}^{\left( s\right) },s \in c}\right) = \mathop{\sum }\limits{{s \in c}}\operatorname{softmax}\left( {{\mathbf{x}}{j}^{\left( s\right) }{\mathbf{W}} + b\left( {\phi }{c,s}\right) }\right) {\mathbf{x}}^{\left( s\right) }{\mathbf{W}}_{\mathbf{v}}\text{ , }]

其中( {\mathbf{W}}{\alpha } \in {\mathbb{R}}^{V \times 1},{\mathbf{W}} \in {\mathbb{R}}^{V \times D} )为可学习权重矩阵，( b : \mathbb{R} \rightarrow \mathbb{R} )为进化距离编码函数（详见“系统发育信息引导的交叉注意力”一节），( {\phi }{c,s} )为物种( s )与支系( c.{x}^{\left( s\right) } )中所有物种最近共同祖先（MRCA）之间的系统发育距离；( j )为源序列( t )在位置( j )处的独热编码核苷酸。来自( A )个注意力头的输出嵌入经拼接后，再通过一个线性层变换为( H )维嵌入，即该模块的最终输出。该注意力池化过程在聚合进化信号的同时降低了下游计算复杂度。若某支系仅含单一物种，则该注意力池化模块退化为单一嵌入模块，其结构与目标序列嵌入模块类似。由此得到支系级源序列嵌入( {\left{ {\mathbf{E}}^{\left( c\right) }\right} }_{c \in \mathcal{C}} )，其中( {\mathbf{E}}^{\left( c\right) } \in {\mathbb{R}}^{L \times H} )。

GPN-Star编码器。GPN-Star采用纯编码器架构，其编码器由( K )个编码器块堆叠而成。每个块依次包含：

1. 序列级自注意力模块。

2. 系统发育信息引导的交叉注意力模块。

3. 前馈网络。

序列嵌入依次通过各编码器块。第一与第二注意力模块各自具有( A )个注意力头。各模块均采用标准残差连接，且嵌入在送入每个模块前经LayerNorm归一化，并全程应用Dropout（概率0.1）。第三模块与标准Transformer相同：它是一个双层全连接网络，接收来自( A )个注意力头的拼接嵌入，输出维度为( H )的嵌入，中间层维度为( F )，并使用GELU激活函数。下文将更详细地描述前两个专用注意力模块。

序列级自注意力。该模块沿基因组捕获局部核苷酸上下文，并应用于各目标序列。由于比对中各物种基因组存在大量重排与片段化，位置上下文在主要目标物种中定义最为可靠。因此，自注意力中的查询（query）与键（key）仅由主要目标物种的嵌入计算得出，随后将注意力分数应用于所有目标物种嵌入所得的值（value）张量。对于每个目标序列嵌入，单个注意力头的输出为：

[{\operatorname{Attn}}{\mathrm{{seq}}}\left( {\mathbf{E}}^{\left( t\right) }\right) = \operatorname{softmax}\left( \frac{\left( {{\mathbf{E}}^{\left( {t}\right) }{\mathbf{W}}{Q}^{\mathrm{{seq}}}}\right) {\left( {\mathbf{E}}^{\left( {t}\right) }{\mathbf{W}}{K}^{\mathrm{{seq}}}\right) }^{\top }}{\sqrt{D}}\right) {\mathbf{E}}^{\left( t\right) }{\mathbf{W}}^{\mathrm{{seq}}},]

其中 ( {\mathbf{W}}{Q}^{\text{ seq }},{\mathbf{W}}^{\text{ seq }},{\mathbf{W}}{V}^{\text{ seq }} \in {\mathbb{R}}^{H \times D} ) 为该注意力头可学习的键、查询与值权重矩阵，( {\mathbf{E}}^{\left( t\right) } \in {\mathbb{R}}^{L \times H} ) 为前一层中目标物种 ( t ) 的序列嵌入，( {\mathbf{E}}^{\left( {t}\right) } \in {\mathbb{R}}^{L \times H} ) 为前一层中主要目标物种 ( {t}{0} ) 的序列嵌入。我们在注意力模块中采用旋转位置编码（RoPE）[89]，以主要目标物种基因组为参考嵌入基因组位置；为便于描述，上述公式中省略了该编码。此外，跨 ( {N}{\mathcal{T}} ) 条序列共享单一注意力图，显著降低了内存占用，从 ( O\left( {{N}_{\mathcal{T}}{L}^{2}}\right) ) 减少至 ( O\left( {L}^{2}\right) )，这与MSA Transformer [11] 中引入的另一种方法类似。

系统发育信息引导的交叉注意力。该模块显式建模从源序列到目标序列的进化过程，同时将系统发育信息直接整合进注意力计算。交叉注意力在目标序列嵌入（用于计算查询张量）与聚合后的源支系嵌入（用于计算键与值张量）之间进行。单个注意力头的输出为：

[{\operatorname{Attn}}{\mathrm{{phy}}}\left( {{\mathbf{E}}^{\left( \mathcal{T}\right) },{\mathbf{E}}{j}^{\left( \mathcal{C}\right) }}\right) = \operatorname{softmax}\left( {\frac{\left( {{\mathbf{E}}^{\left( \mathcal{T}\right) }{\mathbf{W}}{Q}^{\mathrm{{phy}}}}\right) {\left( {\mathbf{E}}^{\left( \mathcal{C}\right) }{\mathbf{W}}{K}^{\mathrm{{phy}}}\right) }^{\top }}{\sqrt{D}} + b\left( \Phi \right) }\right) {\mathbf{E}}^{\left( C\right) }{\mathbf{W}}_{V}^{\mathrm{{phy}}},]

其中 ( {\mathbf{W}}{Q}^{\text{ phy }},{\mathbf{W}}^{\text{ phy }},{\mathbf{W}}{V}^{\text{ phy }} \in {\mathbb{R}}^{H \times D} ) 为该注意力头可学习的键、查询与值权重矩阵，( {\mathbf{E}}^{\left( \mathcal{T}\right) } \in {\mathbb{R}}^{{N}{\mathcal{T}} \times H} ) 为前一模块中第 ( j ) 列对齐位置的目标嵌入，( {\mathbf{E}}^{\left( \mathcal{C}\right) } \in {\mathbb{R}}^{{N}{\mathcal{C}} \times H} ) 为第 ( j ) 列处的支系级源嵌入，( \Phi \in {\mathbb{R}}^{{N}{\mathcal{T}} \times {N}_{\mathcal{C}}} ) 为各目标物种与各支系内源物种之间的平均系统发育距离。为编码各目标物种与源物种间的进化距离概念，我们将传统注意力模块中的位置编码替换为进化距离编码，即上式中的偏置项 ( b\left( \Phi \right) )。我们适配了“相对位置编码的功能插值法”（FIRE）[90] 来编码进化距离：

[b\left( \phi \right) = {f}{\theta }\left( \frac{g\left( \phi \right) }{g\left( {\phi }\right) }\right) ,]

其中 ( g : x \mapsto \log \left( {{cx} + 1}\right) ) 和 ( c \in {\mathbb{R}}^{ + } ) 是可学习的缩放因子，( {\phi }{\max } ) 是比对中所有物种对之间的最大距离，( {f} : \mathbb{R} \rightarrow \mathbb{R} ) 是一个双层感知机 ( {f}{\theta }\left( x\right) = ; {\mathbf{v}}^{\top }\sigma \left( {{\mathbf{v}}{1}x}\right) )，其中 ( \theta = \left{ {{\mathbf{v}},{\mathbf{v}}_{2} \in {\mathbb{R}}^{R}}\right} ) 为可学习参数，( \sigma ) 为 SiLU 激活函数。该进化距离编码模块在同层的所有注意力头之间共享。

选择该距离编码方案主要有两个原因：第一，与序列上位置间的距离不同，物种间的系统发育距离遵循树状分支路径，并非沿单一线性路径可加；因此，相对位置编码方案通过仅依据成对距离向注意力分数添加偏置项，提供了极大便利；第二，我们希望注意力分数能以灵活方式依赖于物种间的进化距离，而所采用的 FIRE 方法正通过使用双层感知机学习从原始进化距离到注意力计算中偏置项的灵活映射，实现了这一目标。

如前所述，同一分支（clade）内的物种在进化上关系密切，序列相似性极高。因此，为使模型能从有意义的进化过程中学习，我们在该交叉注意力模块中应用了分支内注意力掩码，以防止目标物种的 token 注意到其自身所属分支的嵌入表示。

输出层。在堆叠的 ( K ) 个 GPN-Star 编码器块之后，输出层是一个标准的掩码语言建模头，用于将隐状态嵌入转换为每个 token 的核苷酸概率：

[\widehat{p}\left( {{x}{j}^{\left( t\right) } \mid {\mathbf{E}}^{\left( t\right) }}\right) = \operatorname{softmax}\left( {{\mathbf{W}}{o}{\mathbf{E}}^{\left( t\right) } + {\mathbf{b}}_{o}}\right) ,]

其中 ( {W}{o} \in {\mathbb{R}}^{H \times V} ) 和 ( {b} \in {\mathbb{R}}^{V} ) 分别为该模块中可学习的权重矩阵和偏置项，( {e}_{j}^{\left( t\right) } ) 是第 ( t ) 个物种在位置 ( j ) 处、来自最终编码器块的嵌入。预测得到的核苷酸概率随后用于计算训练损失或进行推理。

训练设置

模型基于掩码语言建模目标进行训练。对每个训练样本，我们按如下方式从比对区块构建目标序列：主目标物种将出现在每个样本中；随后，我们通过分层随机抽样选取 19 个目标物种以校正系统发育偏差：首先，从全部分支中抽取 19 个分支，若分支总数少于 19 则允许重复抽样，否则不重复抽样；其次，从每个被选中的分支中各抽取一个物种。由此，每个训练样本共含 20 个目标物种。至于源序列，则如前所述使用全部物种。

接下来我们对输入序列施加掩码。在每个训练样本中，我们在每个目标序列中随机选取 15% 的非空位（non-gap）位置，并确保同一分支内各物种在相同位置被选取。我们将所选 token 的索引集合记为 ( \mathcal{M} )。随后，其中 90% 的被选 token 随机替换为 [MASK] token，剩余 ( {10}% ) 则保持不变，从而得到一组掩码后目标序列 ( {\left{ {\widetilde{\mathbf{x}}}^{\left( t\right) }\right} }{t \in \mathcal{T}} )。源序列也以同样方式掩码，即同一分支内的所有序列，在目标序列与源序列间于相同位置施加掩码。掩码后的源序列表示为 ( {\left{ {\widetilde{\mathbf{x}}}^{\left( s\right) }\right} }{s \in \mathcal{S}} )。

掩码序列被用作输入。模型在数据增强后的原始目标序列的选定token上，以加权交叉熵损失进行训练( {\widehat{x}}_{j}^{\left( t\right) } )（详见下文）：

[{\mathcal{L}}{\mathrm{{MLM}}} = - \frac{1}{\left| \mathcal{M}\right| }\mathop{\sum }\limits\right) \in \mathcal{M}}}{w}{j}\log \widehat{p}\left( {{\acute{x}}^{\left( t\right) } \mid {\left{ {\widetilde{\mathbf{x}}}^{\left( t\right) }\right} }{t \in \mathcal{T}},{\left{ {\widetilde{\mathbf{x}}}^{\left( s\right) }\right} }{s \in \mathcal{S}}}\right) ,]

采用与Benegas等人相同的损失加权方案[29]。其目的是降低重复区域的权重、提高保守元件的权重，从而在训练过程中减轻对中性或非功能性区域错误预测的惩罚。PhyloP和PhastCons被用作衡量保守性的指标以计算权重。对于PhastCons，我们构建了一个平滑版本PhastCons( {}_{M} )，该版本取局部7 bp窗口内的最大PhastCons得分。某位点的损失权重定义为：

[w \propto \left( {{0.1} \times \mathbb{1}{ \text{ repeat }} + \mathbb{1}{ \neg \text{ repeat }} }\right) \times \max \left( {\text{ phyloP },1}\right) \times \left( {\text{ phastCons }{}_{M} + {0.1}}\right) .]

所有目标序列中的位点均相对于主目标序列进行定义。未显式施加按物种划分的损失加权。然而，由于我们仅在非空缺（non-gap）位点上计算损失，因此距离主目标物种更远的物种会受到隐式降权——因为随着物种与主目标物种的分化程度增加，比对中空缺token的比例通常会上升。

我们还采用了Benegas等人[29]的数据增强技术。在计算损失前，原始目标序列中每个位点以一定概率( q )被随机核苷酸替换：

[q = {0.5} \times \mathbb{1}\left{ {\text{ phastCons }{}{M} < {p}{\text{ neutral }}} .}\right.]

通过该过程，由PhastCons定义的高置信度中性位点将在不同训练样本中被随机扰动，从而使模型从更多样化且合理的序列中学习。阈值( {p}_{\text{ neutral }} )在除hg38灵长类模型外的所有模型中均设为0.1；而在hg38灵长类模型中，所定义中性位点的比例显著小于其他模型，因此我们将该阈值提高至0.2。

所有模型均在由8块NVIDIA A100 GPU组成的集群上训练。对每个比对，我们选取一小部分染色体作为验证集，并依据预设训练步数内最低的验证损失选择最优模型。不同模型的完整训练配置见附表1。不同规模模型的架构超参数列于附表9。

相同规模的所有预训练模型，其超参数完全一致，仅两个例外：上下文长度（context size）与界定演化支（clade）的系统发育距离阈值( {\phi }{\text{ clade }} )。针对这两个超参数，我们依据比对特征为模型设置了两组取值：本质上，在物种间分化足够大（最大系统发育距离大于1）的比对中，我们采用更大的clade界定阈值( {\phi }{\text{ clade }} )（0.2）和更小的上下文长度（128），因为此类比对高度碎片化；而在所有物种彼此相对接近（最大系统发育距离低于1）的比对中，我们采用更小的( {\phi }_{\text{ clade }} )（0.05）以建模更精细的演化信号，并采用更大的上下文长度（256），因为基因组间重排程度较弱。采用第一组设置的模型包括：hg38-脊椎动物、mm39、galGal6、dm6和ce11；采用第二组设置的模型包括：hg38-哺乳动物、hg38-灵长类和tair10。

推理

使用训练好的模型进行推理时，通常仅以一个目标序列（即主目标序列）作为输入，而所有序列均作为源序列。在计算某一位点的熵或对数似然比（LLR）时，我们选取以该位点为中心、长度等于模型训练上下文大小的比对窗口。输入的目标序列中，中心token被掩码，随后模型预测该位点的核苷酸概率。我们按如下方式计算某变异的对数似然比（LLR）：

[\mathrm{{LLR}} = \log \frac{p\left( \text{ alt alelle }\right) }{p\left( \text{ ref allele }\right) },]

熵如下所示：

[H\left( p\right) = - \mathop{\sum }\limits_{{x \in { A,C,G,T} }}p\left( x\right) \log p\left( x\right)]

沿用既往方法学[7]，我们在多数致病性预测及罕见变异关联检验分析中采用对数似然比（LLR）作为评分函数。在GWAS精细定位数据集中，我们使用LLR的绝对值，因相关研究依据效应大小的绝对值对变异进行优先级排序。在复杂性状遗传力分析中，我们最初亦采用LLR绝对值，但后续发现熵更具信息量——这一结果令人意外，因其丢失了等位基因特异性。一个熵比LLR更具信息量的实例是rs149514089，该位点为平均红细胞体积与血红蛋白的潜在因果变异，位于一远端候选增强子附近[59]（补充图19

补充图19：一个熵值预测能力优于对数似然比（LLR）的潜在因果变异示例（chr6-16167286-C-T）。预测熵值极低，但参考等位基因（C）与替代等位基因（T）的概率均不高。

）。在此变异中，参考等位基因与替代等位基因在多序列比对中的频率均较低，故LLR接近于0。

为获取目标序列基因组窗口的嵌入表示，将所有物种的比对序列用作源序列。此时输入目标序列为未掩蔽状态，并提取最后一层编码器模块的隐藏嵌入作为输出。

对于变异评分及嵌入表示，我们取正向链与反向链预测结果的平均值。

突变率校准

本模型基于自然界中观测到的基因组序列进行训练。这些序列受高度复杂的进化过程中的多种力量塑造。除自然选择外，突变率亦是决定观测序列的关键因素。为更准确地利用本模型估算全基因组选择约束，我们设计了一种简洁而有效的校准方法，以消除模型预测中突变率变异的影响。其基本思想与传统系统发育模型（如PhyloP和PhastCons）类似：先在高置信度中性位点上估计中性模型评分，再以此调整模型预测时的评分。

具体流程如下：首先，收集一组高置信度中性位点。对于脊椎动物基因组，我们通过将中等远缘物种基因组中的重复序列“提升映射”（lift over）至目标物种基因组，并与目标基因组内重复序列取交集，从而识别目标物种的祖先重复序列。具体而言，人类模型采用小鼠重复序列提升映射；小鼠模型采用人类重复序列提升映射；鸡模型则采用斑胸草雀重复序列提升映射。其次，进一步筛选最可靠的中性位点：选取祖先重复序列内PhyloP评分为-0.1至0.1、且PhastCons评分为0的位点。对于非脊椎动物物种，由于祖先重复序列提供的位点过少，不足以可靠估计中性评分，故跳过第一步，直接采用由PhyloP（-0.05至0.05）与PhastCons（0）定义的全基因组中性位点。第二步，在中性位点上生成模型评分（包括LLR与熵）。为获得描述位点特异性突变率的基线模型评分，我们将中性位点按其中心五核苷酸序列上下文分箱，因已知中心k-mer上下文是突变率的主要决定因素[91]。最后，计算每类五核苷酸上下文对应变异的模型评分均值，即由该五核苷酸上下文及突变类型所确定的特定突变率分箱之基线模型评分。我们称其为中性评分（即( {\mathrm{{LLR}}}{\text{ neutral }} )与( {H}{\text{ neutral }} )）。随后，将对应五核苷酸分箱的中性评分用于校正模型预测，公式如下：

[\text{ calibrated LLR } = \text{ LLR } - {\text{ LLR }}_{\text{ neutral }}\left( {x,\text{ alt }}\right)]

[\text{ calibrated entropy } = \frac{H}{{H}_{\text{ neutral }}\left( x\right) }\text{ , }]

式中( x )为目标位点的中心五核苷酸上下文。

如补充图16与补充图20

补充图20：非人源模型在校准突变率前后的评估结果。（A）–（E）群体基因组数据库中罕见变异与常见变异在有害性评分尾部的富集情况（阈值设定为使各评分产生30个假阳性发现），与适配于相同多序列比对的PhyloP和PhastCons比较。（F）MMrdb致病变异与WMGP常见变异在M. musculus中的分类。（G）FlyBase致死变异与DEST常见变异在D. melanogaster中的分类。（H）C. elegans致死变异与CaeNDR常见变异的分类。

所示，该突变率校准流程将GPN-Star评分与Roulette突变率估计值之间的相关性降至可忽略水平，提升了约束力估计精度，并在不同程度上改善了大多数基准测试的表现。

序列嵌入可视化

我们旨在可视化人类常染色体上一组相对均衡的功能性与非功能性基因组区域的嵌入表示。我们选择100 bp的分辨率，该值接近人类编码外显子长度的中位数。分析中排除了所有重复序列，因其可能占据大量表征空间[6]，并妨碍其余更富信息的基因组区域的可视化。我们从UCSC基因组浏览器[84]下载了RepeatMasker[92]重复序列注释。基因注释采用Ensembl第113版[58]，cCRE目录采用ENCODE SCREEN v4[59]。对每种基因组元件类型（如lncRNA或启动子），我们使用Bioframe[93]过滤掉与其他元件类型重叠的区域。我们将不与外显子或cCRE重叠的区域定义为背景区域。对于cCRE，我们聚焦于启动子和远端增强子。当灵长类PhastCons[55]在( {75}^{\text{ th }} )百分位上的值为100%时，我们将对应窗口标记为“保守”。对每种功能元件类型，我们分别采样最多10,000个“保守”区域和10,000个“非保守”区域；并随机采样20,000个背景区域。

100 bp窗口的模型嵌入通过添加最多256 bp侧翼序列获得，计算每个位置的嵌入，并对中心100个位置取平均。正链与负链嵌入取平均后进行标准化。UMAP[94]采用默认参数运行。

核苷酸依赖性分析

本方法基于文献[45]。具体而言，我们计算位置( i )与( j )之间的依赖性，公式如下：

[{e}{i,j} = {\begin{Vmatrix}\left{ \log \left\lbrack \frac{\mathbb{P}\left( {{n} = k \mid {n}{i} = {k}{\text{ alt }}}\right) }{\mathbb{P}\left( {{n}{j} = k \mid {n} = {k}{\text{ ref }}}\right) }\right\rbrack \mid k,{k}{\text{ alt }} \in { \mathrm{A},\mathrm{C},\mathrm{G},\mathrm{T}} ,{k}{\text{ alt }} \neq {k}{\text{ ref }}\right} \end{Vmatrix}}_{\infty }]

其中( {n}{i} )为位置( i )处的核苷酸。计算该得分时未掩蔽任何核苷酸。此外，我们通过对( {e} )与( {e}_{j,i} )取平均实现对称化。

LDLR MPRA数据[68]以预处理形式自文献[95]下载。我们应用论文中所述最严格的筛选条件：变异需至少包含10个条形码，且( p )-值满足( < {10}^{-5} )。为获取每个位置的效应大小度量，我们对不同等位基因的效应大小取平均；对无显著等位基因的位置，其效应大小定义为0。最终，效应大小经平滑处理：对连续5个位置取平均。

评估基准

智人（Homo sapiens）

数据集

• ClinVar[31]：下载20220924版本。目标是区分被标注为“致病性”与“良性”的错义变异，评价指标为受试者工作特征曲线下面积（AUROC）。

• COSMIC[32]：下载经文献[29]处理的数据。目标是区分COSMIC高频错义变异与gnomAD[33]常见错义变异，评价指标为精确率-召回率曲线下面积（AUPRC）。

• ProteinGym[34]：下载ProteinGym v0.1中的31个人类蛋白质数据。我们依据Benegas等[29]所述流程处理该数据集。

• OMIM[41]：下载自TraitGym[7]获取的经处理孟德尔性状数据集。目标是区分致病性与常见调控变异，评价指标为AUPRC。

• HGMD[42]：下载自文献[96]整理的数据，并参照TraitGym[7]中孟德尔性状数据集的方式进行处理。目标是区分致病性与常见调控变异，评价指标为AUPRC。

• GWAS精细定位[46, 97]：针对非编码变异，我们从TraitGym[7]下载经处理的复杂性状数据集；对编码变异，则自同一来源[46, 97]处理，方式类似TraitGym[7]，但仅匹配次要等位基因频率（MAF）。目标是区分高与低后验包含概率（PIP）的变异，评价指标为AUPRC。

• gnomAD [33]：gnomAD v3.1.2 等位基因频率数据从其官网 https://gnomad.broadinstitute.org/data 下载。我们遵循 GPN-MSA 论文 [29] 的评估协议：在有害性评分尾部（阈值设定为允许少量假阳性发现，此处设为 30）中，富集稀有变异（单例变异）与常见变异（MAF > 5%）的比例。我们将样本量较小的若干变异后果类别合并如下：“splice-region” 合并 Ensembl 分类中的 splice_donor、splice_acceptor、splice_donor_5th_base、splice_donor_region、splice_region 和 splice_polypyrimidine_trac；“start-or-stop” 合并 Ensembl 分类中的 start_lost、stop_gained 和 stop_lost。为减少待评分变异数量，我们在每个类别中对稀有变异进行下采样，使其数量与对应类别中常见变异数量一致。

我们还使用了基于 gnomAD v3 数据推导出的 Gnocchi [98] 约束得分。经质量控制筛选后的 ( 1\mathrm{;{kb}} ) 窗口的预计算全基因组 z 得分，从同一网站下载。得分范围为 -10 至 10。我们将窗口按长度为 1 的得分区间下采样至每区间 100 个窗口，共得 2,000 个窗口。在分析中，我们比较每个窗口内所有可能变异的 GPN-Star、PhyloP 和 PhastCons 得分均值与对应 Gnocchi z 得分。

对比模型

• AlphaGenome [40]：通过官方 API 获取得分，遵循 TraitGym [7] 中的 Enformer/Borzoi 协议，但将跨 track 的 L2 聚合替换为最大值操作（如文献 [40] 所推荐）。

• Borzoi [39]：下载 OMIM 和 GWAS 精细定位数据集 [7] 的预计算得分，并对 HGMD 采用类似方法。

• CADD [95]：下载预计算的 CADD v1.7 原始得分。

• ESM-1b [79]：由参考文献 [99] 预计算的对数似然比（LLR）得分，以处理后的形式从参考文献 [29] 下载。

• Enformer [38]：下载 OMIM 和 GWAS 精细定位数据集 [7] 的预计算得分，并对 HGMD 采用类似方法。对于 LDSC，我们从参考文献 [61] 下载预计算得分，具体包括组织无关特征 Enformer.Enformer_all_all 和组织特异性特征 Enformer.Enformer__all。

• Evo-2 [9]：按原文方法，使用 40B 和 7B 模型、8192 bp 上下文长度计算 LLR。

• GPN-MSA [29]：从 https://huggingface.co/datasets/songlab/gpn-msa-hg38-scores 下载预计算 LLR。

• Nucleotide Transformer [8]：使用 2.5B 多物种模型计算 LLR。

• PhastCons（V）：基于 100 种脊椎动物，从 https://hgdownload.soe.ucsc.edu/goldenPath/hg38/phastCons100way/hg38.phastCons100way.bw 下载。

• PhastCons（M）：基于 470 种哺乳动物，从 https://hgdownload.soe.ucsc.edu/goldenPath/hg38/phastCons470way/hg38.phastCons470way.bw 下载。

• PhastCons（P）[55]：基于 43 种灵长类动物，从 https://cg1.gi.ucsc.edu/data/cactus/zoonomia-2021-track-hub/hg38/phyloPPrimates.bigWig 下载。

• PhyloP（V）：基于 100 种脊椎动物，从 https://hgdownload.soe.ucsc.edu/goldenPath/hg38/phyloP100way/hg38.phyloP100way.bw 下载。

• PhyloP（M）[55]：基于 447 种哺乳动物，从 https://hgdownload.soe.ucsc.edu/goldenPath/hg38/phyloP447way/hg38.phyloP447way.bw 下载。

• PhyloP（P）：基于243种灵长类物种，数据下载自 https://hgdownload.soe.ucsc.edu/goldenPath/hg38/phyloP447way/hg38.phyloP447wayPrimates.bw

• AlphaMissense [35]：下载了针对hg38基因组变异的预计算得分。

• PrimateAI-3D [36]：通过学术许可从Illumina获取得分。

• SpeciesLM [45]：使用在动物基因组上训练的模型计算LLR。

• PromoterAI [44]：通过学术许可从Illumina获取得分。我们首先取得分的绝对值（以捕捉表达量变化的幅度，而非仅过表达或低表达方向）。若某变异被预测影响多个基因，则取所有基因中表达量变化的最大值。

• GPN-Promoter [7]：使用来自 https://github.com/songlab-cal/gpn 的预训练模型计算LLR。

• Roulette [71]：全基因组突变率预估数据从GitHub仓库 https://github.com/vseplyarskiy/Roulette 下载。

对所有我们计算的得分，均取正向链与反向链预测结果的平均值。

小家鼠（Mus musculus）

数据集

• 野生小鼠基因组计划（WMGP）：群体等位基因频率数据源自一组已发表的WMGP数据集 [100–107]（ https://andrewparkermorgan.github.io/wmgp/ ）。原始数据使用fastp [108] 进行接头序列去除和质量修剪。所得reads比对至一个由GRCm39小鼠参考基因组与协作杂交作图群体（Collaborative Cross mapping panel）七种额外奠基者基因组的PacBio连续长读长（CLR）组装、以及BALB/cByJ、BALB/cJ、C3H/HeJ、C3H/HeOuJ、C57BL/6NJ、DBA/2J和PWD/PhJ品系的PacBio CLR组装共同构建的minigraph-cactus [109]泛基因组，比对工具为VG giraffe [111]。所得比对结果利用VG [112] 投影至GRCm39坐标空间并输出为BAM格式。全基因组变异使用DeepVariant [113]（WGS模型）对每个样本进行识别；各单一样本gVCF文件经合并后，采用GLnexus [114] 进行联合变异检出。在富集分析中，我们将单倍型（singletons）定义为稀有变异，等位基因频率高于0.2的变异定义为常见变异。分子效应注释流程与gnomAD数据集一致。为减少打分变异数量，我们也采用相同流程对稀有变异进行亚抽样。

• 小鼠突变资源数据库（MMRdb）：致病性变异整理自最新版MMRdb。该库收录了杰克逊实验室（Jackson Laboratory）培育的300多种实验小鼠品系中自发疾病表型（符合孟德尔遗传）的变异检出结果。根据Fairfield等人[115]先前基于外显子组与全基因组数据集分析提出的、用于识别潜在孟德尔疾病相关变异的一套标准，包括：变异在数据库中为稀有变异 ( \left( { < 3% }\right) )；样本中观察到的等位基因数与预期基因型相符（杂合子为1个等位基因，纯合子为2个）；且变异染色体位置与现有连锁分析数据一致。部分变异源自早期版本MMRdb [116]，其余则来自MMRdb v.2.0.0（可访问 https://mmrdb.jax.org/ ）。该数据集中所有变异均以GRCm38参考基因组坐标记录。为适配本模型所用GRCm39参考坐标空间，我们使用UCSC liftOver工具及对应chain文件 [117] 进行转换。我们过滤掉非SNV变异及参考等位基因与GRCm39参考基因组不一致的变异。仅保留以下变异类型（因数量充足）：错义、同义、无义、3' UTR、5' UTR及剪接区变异。该人工审编数据集共包含470个致病性或疑似致病性变异。其中89个致病性变异经PCR验证（Fairfield等2015年报道58个；MMRdb v.2.0.0中报道31个），其余381个被归类为疑似致病性变异。我们在基准测试中，使用WMGP对应类别中的常见变异作为阴性对照。

对比模型

• PhastCons：基于35种脊椎动物，数据下载自 https://hgdownload.soe.ucsc.edu/goldenPath/mm39/phastCons35way/mm39.phastCons35way.bw 。

• PhyloP：基于35种脊椎动物，数据下载自 https://hgdownload.soe.ucsc.edu/goldenPath/mm39/phyloP35way/mm39.phyloP35way.bw 。

黑腹果蝇

数据集

• 果蝇时空进化（DEST）[118]：等位基因频率数据下载自 https://berglandlab.pods.uvarc.io/vcf/dest.all.PoolSNP.001.50.24Aug2024.ann.vcf.gz 。在我们的富集分析中，我们将等位基因频率低于0.002的变异定义为稀有变异，将等位基因频率高于0.2的变异定义为常见变异（该数据集中无单例变异）。我们采用与gnomAD数据集相同的流程对变异的分子效应进行注释。

• FlyBase [75]：从FlyBase数据库 https://flybase.org/（FB2025_02版，2025年4月17日发布）收集经实验验证的致死突变。我们使用关键词“lethal”、“point_mutation”和“na_change”检索数据库，获取导致致死表型的已注释SNV。我们发现绝大多数( \left( { > {99}% }\right) )变异属于错义、无义或剪接变异；因此仅保留这三类。在将参考等位基因与dm6参考基因组比对后，共获得2,929个致死变异。我们在基准测试中，以DEST数据集中对应类别中的常见变异作为阴性对照。

对比模型

• PhastCons：基于124种昆虫，数据下载自 https://hgdownload.soe.ucsc.edu/goldenPath/dm6/phastCons124way/dm6.phastCons124way.bw 。

• PhyloP：基于124种昆虫，数据下载自 https://hgdownload.soe.ucsc.edu/goldenPath/dm6/phyloP124way/dm6.phyloP124way.bw 。

秀丽隐杆线虫

数据集

• 秀丽隐杆线虫自然多样性资源（CaeNDR）[119]：等位基因频率数据下载自 https://storage.googleapis.com/caendr-site-public-bucket/dataset_release/c_elegans/20231213/variation/WI.20231213.hard-filter.isotype.vcf.gz 。在我们的富集分析中，我们将等位基因计数为2（即一个纯合品系）的变异定义为稀有变异，将等位基因频率高于0.2的变异定义为常见变异。我们采用与gnomAD数据集相同的流程对变异的分子效应进行注释。

• 实验验证的致死变异：源自Qin等[76]研究的补充数据。我们使用该数据集中的72个SNV。我们在基准测试中，以CaeNDR数据集中对应类别中的常见变异作为阴性对照。

对比模型

• PhastCons：基于135种线虫，数据下载自 https://hgdownload.soe.ucsc.edu/goldenPath/ce11/phastCons135way/ce11.phastCons135way.bw 。

• PhyloP：基于135种线虫，数据下载自 https://hgdownload.soe.ucsc.edu/goldenPath/ce11/phyloP135way/ce11.phyloP135way.bw 。

红原鸡

数据集

• Galbase [120]：从 http://animal.omics.pro/code/source/download/Chicken/variation/GRCg6a_SNPs.anno.tab.gz 下载等位基因频率数据。在我们的富集分析中，我们将等位基因频率为 0.001（该数据集中最小值）的变异定义为稀有变异，将等位基因频率高于 0.2 的变异定义为常见变异。我们采用与 gnomAD 数据集相同的流程注释变异的分子效应。

对比模型

• PhastCons：基于 77 种脊椎动物，从 https://hgdownload.soe.ucsc.edu/goldenPath/galGal6/phastCons77way/galGal6.phastCons77way.bw 下载。

• PhyloP：基于 77 种脊椎动物，从 https://hgdownload.soe.ucsc.edu/goldenPath/galGal6/phyloP77way/galGal6.phyloP77way.bw 下载。

拟南芥（Arabidopsis thaliana）

数据集

• 1001 基因组计划 [121]：从 https://huggingface.co/datasets/gonzalobenegas/processed-data-arabidopsis/resolve/main/variants/all/variants.parquet 下载 GPN 论文 [6] 中处理过的等位基因频率数据。在我们的富集分析中，我们将单倍型（singletons）定义为稀有变异，将等位基因频率高于 0.2 的变异定义为常见变异。我们采用与 gnomAD 数据集相同的流程注释变异的分子效应。

对比模型

• PhastCons：基于 18 种十字花科植物，从 PlantRegMap [87] 下载： http://plantregmap.gao-lab.org/download_ftp.php?filepath=08-download/Arabidopsis_thaliana/sequence_conservation/Ath_PhastCons.bedGraph.gz 。

• PhyloP：基于 18 种十字花科植物，从 PlantRegMap [87] 下载： http://plantregmap.gao-lab.org/download_ftp.php?filepath=08-download/Arabidopsis_thaliana/sequence_conservation/Ath_PhyloP.bedGraph.gz 。

遗传力富集分析

分层连锁不平衡评分回归（S-LDSC）背景

以下是对 S-LDSC [52] 的简要介绍。

符号说明。我们首先引入部分符号：

• ( N ) — GWAS 样本量（个体数）。

• ( M ) - 分析的SNP数量。

• ( \mathbf{X} \in {\mathbb{R}}^{N \times M} ) - 标准化基因型矩阵。

• ( \mathbf{y} \in {\mathbb{R}}^{N} ) - 标准化表型。

• ( \mathbf{\beta } \in {\mathbb{R}}^{M} ) - 真实加性效应（其中( {\beta }_{j} )表示SNP ( j ) 的效应）。

• ( {a}_{jc} ) - 注释( c )在SNP ( j ) 处的取值（二元型、连续型或概率型）。

• ( {\tau }_{c} ) - 每单位注释( c )对遗传力的单SNP贡献（在给定所有其他注释条件下的偏效应）。

• ( \varepsilon ) - 环境残差向量，其( \mathcal{N}\left( {0,{\sigma }_{e}^{2}}\right) )个分量相互独立。

模型：在加性多基因假设下，表型可表示为：

[\mathbf{y} = \mathbf{X}\mathbf{\beta } + \varepsilon]

我们假设每个SNP效应均独立采样自零均值高斯分布，其方差随功能注释线性变化：

[{\beta }{j} \sim \mathcal{N}\left( {0,\mathop{\sum }\limits{{c = 1}}^{C}{a}{jc}{\tau }{c}}\right) .]

总SNP遗传力计算如下：

[{h}{g}^{2} = \mathop{\sum }\limits{{j = 1}}^{M}\operatorname{Var}\left( {\beta }{j}\right) = \mathop{\sum }\limits{{c = 1}}^{C}{\tau }{c}\mathop{\sum }\limits{{j = 1}}^{M}{a}_{jc}.]

推断：参数估计值( {\widehat{\tau }}_{c} )通过将全基因组关联分析的卡方统计量对分层LD评分进行（广义）最小二乘回归获得。主要统计难点在于连锁不平衡所导致的SNP水平统计量之间存在强协方差；S-LDSC通过将LD压缩为每个SNP的评分，并在基因组窗口上采用分块Jackknife法计算稳健标准误，从而规避该问题。

输出：对每个注释( c )，S-LDSC报告两个汇总其与遗传力关系的估计统计量：

• 系数( {\widehat{\tau }}_{c} )：在调整所有其他注释的前提下，估计的单位注释增加所对应的每SNP遗传力贡献；该系数的无量纲版本[122]表示注释每增加一个标准差所引起的每SNP遗传力的比例变化。

[{\widehat{\tau }}{c}^{ \star } = \frac{{\widehat{\tau }}\operatorname{sd}\left( {a}{\cdot c}\right) }{{\widehat{h}}^{2}/M},]

以平均每SNP遗传力为基准进行标准化，便于跨性状比较。

• 富集度（仅适用于二元注释）：

[{\operatorname{Enrichment}}{c} = \frac{{\widehat{h}}^{{2}/{\widehat{h}}_{g}}{2}}{\mathop{\sum }\limits_{{j = 1}}^{M}{a}_{jc}/M},]

其中( {\widehat{h}}{c}^{2} = \mathop{\sum }\limits{jc}\widehat{\operatorname{Var}}\left( {\beta }\right) )为携带该注释的SNP所贡献的遗传力（包括通过与其他注释重叠所介导的贡献）。富集度将遗传力占比与SNP占比进行比较；当( {a}_{jc} \in { 0,1} )时该指标有定义。

S-LDSC分析

参考文件、基线注释及汇总统计均下载自文献[122]；我们采用baselineLD_v2.2作为基线注释。在具备汇总统计的107个独立性状中，我们保留了其中106个无使用限制的性状。为评估某一模型，我们运行S-LDSC，将该模型评分与基线注释一并纳入。运行S-LDSC需计算约( {10}\mathrm{M} )个变异位点的评分，因此该分析仅适用于可扩展性最强的模型。具体而言，S-LDSC使用来自1000基因组计划第三阶段[123]中479名欧洲个体的9997231个参考变异位点，其等位基因计数为( \geq 5\left( {\mathrm{{MAF}} \geq {0.52}% }\right) )；其中5961159个为常见变异（MAF ( \geq ) 5%），在遗传力估计中起核心作用（详见文献[52]）。我们依据分位数阈值对评分进行二值化处理。由于富集度仅在常见变异上计算，故仅使用这些变异确定阈值，以确保所有模型注释包含相同比例的变异位点。我们对各性状的富集度与系数进行随机效应荟萃分析；显著性检验采用单侧Wald检验。组织特异性分析复用LDSC-SEG[62]中的注释（gs://broad-alkesgroup-public-requester-pays/LDSCORE/LDSC_SEG_ldscores/Multi_tissue_gene_expr_1000Gv3_ldscores.tgz）。本研究中，各GTEx[73]组织特异性区域的定义方法为：首先识别该组织中前( {10}% )个差异高表达最显著的基因，再在这些基因上下游各延伸100 kb形成窗口。为更好匹配Enformer的组织聚合策略[61]，我们对部分组织取其窗口并集进行合并（补充表10）。最终，组织特异性GPN-Star评分定义为仅选取组织特异性区域内排名最靠前的变异位点。在比较组织无关型与组织特异性GPN-Star注释时，我们确保二者包含相同数量的变异位点。需注意的是，PhastCons（P）所用比对物种数（43种灵长类）远少于结果中呈现的GPN-Star（P）（243种灵长类）；我们通过仅使用两套物种集合的交集（36种）重新训练GPN-Star（P），验证了其优势并非主要源于物种数量增加（补充图21

补充图21：使用不同物种数量的灵长类全基因组比对（WGA）训练的灵长类模型比较。基准测试包括图2中的变异效应预测任务（（A）–（G））以及图3中的复杂性状遗传力分析（（H）–（I））。模型名称括号内的数字表示比对中包含的灵长类物种数：P243代表完整cactus447way比对中的243个物种；P43代表初始Zoonomia发布中cactus241way比对的43个物种；P36代表上述两个集合交集中的36个物种。

）。

基于DeepRVAT的罕见变异关联检验

DeepRVAT在已发表注释集基础上增强GPN-Star后的基准测试，遵循文献[49]中图4a所用的相同基准测试流程；此处简要说明细节。

UK Biobank 全外显子组测序（WES）数据按文献[49]中“UKBB 200k无关欧洲血统数据集”的流程进行处理与质量控制。该数据集源自 UK Biobank 20 万例中期 WES 数据，仅纳入具有欧洲遗传血统的个体，并排除三代及以内亲缘关系者，最终获得 161,822 名个体的队列。此限制旨在基准测试中消除人群分层带来的混杂效应。测序数据处理、质量控制及注释均遵循文献[49]，训练所用变异保留其等位基因频率（MAF）( < 1% )，关联检验所用变异保留其 MAF ( < {0.1}% )。

DeepRVAT 模型（软件包 v1.1）分别在有无 GPN-Star 变异注释条件下，采用三种不同随机种子的随机参数初始化方式训练，并使用与文献[49]相同的 21 个训练表型。负荷检验按文献[49]图 4a 所述方法进行，使用 REGENIE 对相同 34 个定量表型（其中 21 个为训练所用，13 个未用于训练）计算 DeepRVAT 基因得分。

DeepRVAT 发现指经 Bonferroni 校正后 ( p ) 值 ( < {0.05} ) 的结果。复制验证通过将显著的基因-表型关联与两项基于更大规模 UK Biobank WES 队列的研究[50, 51]（样本量分别为 394,841 和 454,787 人）所报道的发现进行比对完成。

数据可用性

预训练模型、训练数据集及基准测试数据集已在 Hugging Face 平台公开（ https://huggingface.co/collections/songlab/gpn-star-68c0c055acc2ee51d5c4f129 ）。人类基因组及五种模式生物全基因组所有可能单核苷酸变异的预测结果将在本文发表后公开提供。

代码可用性

用于模型训练、推理及复现主要分析的代码已在 GitHub 公开（ https://github.com/songlab-cal/gpn ）。

补充表格

补充图9：按GPN-Star (P) 定义的前( {0.1}% )约束性最强变异中各类变异后果的比例，及其相对于99.9%约束性最弱变异（全部S-LDSC变异，含常见及低频变异）的优势比。

补充表9：不同模型规模下GPN-Star架构的超参数。

模型尺寸	25M	85M	200M
编码器块数量（K）	8	12	16
隐藏层维度（H）	512	768	1024
中间层尺寸（F）	2048	3072	4096
自注意力头数量（A）	4	6	8
交叉注意力头数量（A）	4	6	8
注意力头维度（D）	64	64	64
Dropout 率	0.1	0.1	0.1

	GTEx 组织
脑	脑黑质 脑脊髓（颈段 C-1） 脑杏仁核 脑前扣带回皮层（BA24） 脑海马体 脑下丘脑 脑壳核（基底神经节） 脑小脑半球 脑额叶皮层（BA9） 脑伏隔核（基底神经节） 脑皮层 脑尾状核（基底神经节） 脑小脑
血液/免疫	脾脏 EBV转化淋巴细胞、全血细胞
肝脏	肝脏
肠道	小肠末端回肠、结肠乙状部 食管、胃食管交界处、胃 结肠横部 食管肌层 食管黏膜
肾脏	肾皮质
脂肪	内脏脂肪（大网膜） 皮下脂肪
心脏	心脏右心耳 心脏左心室
皮肤	非日晒暴露皮肤（耻骨上） 转化成纤维细胞 日晒暴露皮肤（小腿下部）
肺	肺

补充表4：根据GPN-Star（P）确定的前( {0.1}% )个约束性最强变异体（所有S-LDSC变异体，含常见及低频变异）的后果分布。比值比及( p )值采用双侧Fisher精确检验计算。加粗：在错误发现率（FDR）< 5%下显著。

后果	计数	比例	比值比
错义突变	6494	0.40	193.30
dELS	3734	0.23	2.02
dELS侧翼区	1092	0.07	0.30
3’非翻译区	798	0.05	4.25
非编码RNA	654	0.04	1.76
pELS	545	0.03	2.01
5’非翻译区	423	0.03	9.52
CA	322	0.02	0.82
PLS	255	0.02	6.72
内含子	235	0.01	0.06
CA侧翼区	188	0.01	0.16
同义突变	164	0.01	2.97
pELS侧翼区	143	9e-03	0.50
基因间区	138	9e-03	0.04
CA-CTCF	115	7e-03	0.67
CA-H3K4me3	106	7e-03	0.85
剪接区	100	6e-03	6.63
剪接受体位点	94	6e-03	58.71
CA-CTCF侧翼区	79	5e-03	0.16
剪接受体区	60	4e-03	18.68
终止密码子获得	51	3e-03	56.16
CA-H3K4me3侧翼区	51	3e-03	0.18
剪接受体第5位碱基	44	3e-03	38.80
转录因子侧翼区	42	3e-03	0.09
剪接多嘧啶区	38	2e-03	3.17
CA-TF	35	2e-03	0.85
剪接受体位点	30	2e-03	30.80
CA-TF侧翼区	19	1e-03	0.18
转录因子	17	1e-03	0.10
起始密码子丢失	16	1e-03	84.74
PLS侧翼区	7	4e-04	0.44
成熟miRNA	7	4e-04	38.37
终止密码子丢失	6	4e-04	26.35
下游基因	4	2e-04	0.07
上游基因	2	1e-04	0.06