利用泛基因组解析拷贝数变异实现重复基因的序列分型与平行同源体特异性表达差异研究

《Nature Genetics》：Genotyping sequence-resolved copy number variation using pangenomes reveals paralog-specific global diversity and expression divergence of duplicated genes

本研究提出创新性基因分型方法ctyper，通过将测序读段直接比对至泛基因组单倍型，实现拷贝数变异（CNV）和复杂基因的序列分辨率分型。该方法在3,351个CNV基因和212个临床相关基因中分别实现≥99.1%的拷贝数准确率和94.8%的相位变异捕获率，单CPU仅需1.5小时即可完成全基因组分型。研究揭示7.94%的平行同源体存在等位基因特异性表达差异，为脊髓性肌萎缩症（SMA）和代谢疾病机制研究提供新视角。

　　

Abstract

拷贝数变异（CNV）基因在进化与疾病中具有重要作用，但其序列变异仍是大规模研究中的盲区。研究团队开发了ctyper方法，利用泛基因组从二代测序样本中生成包含局部相位变异的等位基因特异性拷贝数。通过对3,351个CNV基因和212个具有临床挑战性的医学相关（CMR）基因进行基准测试，ctyper在CNV基因中捕获了96.5%的相位变异（正确率≥99.1%），在CMR基因中捕获了94.8%的相位变异。单CPU基因组分型仅需1.5小时。与已知表达数量性状位点（eQTL）变异相比，ctyper基因型对基因表达的预测能力提升4.81倍。等位基因特异性表达分析显示7.94%的平行同源体存在差异表达，4.68%存在组织特异性偏好。研究发现SMN1基因转换导致SMN2表达降低（可能影响脊髓性肌萎缩症），以及易位重复的AMY2B基因表达增加。该研究实现了CNV和CMR基因的生物样本库规模分型。

Main

基因分型方法概述

人类基因组频繁通过复制和缺失发生突变，导致拷贝数变异。约10%的蛋白质编码基因存在拷贝数变异，且在不同人群中呈现特异性分布。虽然全基因组范围内CNV频率较低，但被称为分段重复的长低拷贝重复区域在基因中富集，并且由于非等位基因同源重组成为 recurrent CNV 的催化剂。这些区域包括TBC1D3、NPIPA1（NPIP）和NBPF1（NBPF），与脑功能和适应性相关，且比非重复DNA经历更快的核苷酸替换。因此，仅考虑重复基因聚合拷贝数（aggreCN）的研究会遗漏非参考基因重复和多等位基因CNV之间的变异。

现有CNV检测工具依赖过量测序覆盖度，掩盖了拷贝间的变异。此外，NGS比对存在模糊性和偏倚，会遗漏 divergent 非参考重复。新型高质量泛基因组组装提供了来自多样化人群的序列分辨率CNV，包括非参考重复。虽然基于图的泛基因组减少了参考偏倚，但可能将相似序列（包括替代等位基因和功能不同的平行同源体）合并为共享序列，从而模糊了平行同源体特异性变异和序列多样性。随着泛基因组规模扩大，替代序列、基因转换和重排在图谱中的表示将更具挑战性，这促使需要新工具和表示方法来利用NGS和泛基因组分析重复基因。

ctyper方法通过将NGS读段与泛基因组单倍型进行比较，识别NGS样本与单倍型组装之间最相似的基因组片段并分配拷贝数状态。未合并的单倍型序列保留了局部相位变异，并捕获了复杂变异，包括结构变异和基因转换。研究聚焦于使用NGS和参考基因组难以分析的复杂重复基因。利用无比对技术和多项式时间分型模型，ctyper实现了高准确性和计算效率，足以满足未来生物样本库分析的需求。

方法概述

研究将变异表示为单倍型片段，其长度足够短以最小化重组干扰，通过血缘同一性（identity by descent）实现与NGS样本的精确共享，同时捕获包括相位小变异、结构变异和基因转换事件在内的结构信息。这些单倍型片段被标记为泛基因组衍生等位基因（PA），并检测与NGS样本共享的PA。PA边界经过安排以研究蛋白质编码基因的变异。每个PA包含相距<20 kb的连续外显子以及5 kb的侧翼序列，反映了短程转录因子的功能邻近性和群体水平的基因组连锁。PA通常范围在10到100 kb之间，对应于连锁不平衡（LD）区块的规模。

为了实现计算效率并避免重复DNA中的比对模糊性，研究使用无比对方式比较NGS样本中测量的低拷贝k-mer（固定长度k的DNA片段；k=31）来对PA进行分型。对于每个基因，将所有相似的PA（包括直系同源基因、平行同源基因和同源假基因）分组，并构建一个用于分型的矩阵，其中包含所有分组PA的k-mer组成。矩阵的行对应单个PA，列对应分组PA独有的k-mer，单元格值代表每个PA中的k-mer多重性。

分型过程通过识别PA（行）及其拷贝数的组合进行，使它们的k-mer计数与NGS样本之间的最小平方距离最小化。样本k-mer计数被投影到每个k-mer矩阵的向量空间中，并基于PA序列的系统发育使用递归舍入分配整数拷贝数，从而产生PA特异性拷贝数（paCN）列表。

PA数据库构建

研究为先前报道的3,351个CNV基因构建了PA数据库，使用了来自人类泛基因组参考联盟（HPRC）、人类基因组结构变异联盟（HGSVC）、中国泛基因组联盟（CPC）的两个端粒到端粒组装的114个二倍体PacBio HiFi组装，以及GRCh38和CHM13。总共定义了1,408,209个PA，组织成3,307个矩阵。

由于有限的人类遗传多样性和短距离上更强的LD，PA通常高度相似或相同。为了降低维度并便于群体分析，研究使用它们的系统发育关系将相似PA合并为高度相似亚组（亚组），将其视为等效状态。总共定义了89,236个亚组，用于枚举所有PA，类似于人类白细胞抗原（HLA）命名法。

为了注释低频变异和参考基因组位置以确定直系同源或平行同源关系，研究将PA映射到GRCh38。总共确定了6,389个位点上的164,237个平行同源PA。与相应参考位点相似（≥80% k-mer相似性）的平行同源PA被标记为重复性，其余低相似性的平行同源PA被标记为分化性。总共 across 333个矩阵 identified 2,734个亚组中的10,792个分化平行同源基因。这些分化平行同源基因代表了难以进行典型参考分析的新序列。

虽然大多数重复是远端的，但6,673个PA反映了近端（<20 kb）重复，包括36个基因中的1,646个PA表现出“失控重复”，至少有三个近端重复。近端重复基因作为遗传单位被包含在与它们的直系同源基因相同的PA中。直系同源PA如果与参考基因属于同一亚组则被分类为参考等位基因，否则被分类为替代等位基因。所有PA都进行分型，无论平行同源-直系同源注释如何，因此 resulting 基因型包含群体和拷贝数变异。

Ctyper数据库捕获群体多样性

研究评估了PA是否捕获了其他CNV表示无法复制的基因组信息的独特方面，包括参考基因的拷贝数、单独特异核苷酸k-mer（SUNK）和大多倍型结构。研究发现PA提供了更高分辨率的变异（例如单核苷酸变异），因为94.7%的变异未在GRCh38序列中反映。此外，附近的SUNK标记和大倍型结构被发现是PA的不良代理，只有一小部分PA被发现与SUNK或更大倍型结构相关联。尽管维度大幅减少，亚组捕获了超过80%的总群体变异。最后，使用饱和分析，研究估计当前队列代表了非非洲人中98.7%的亚组和非洲人中94.9%的亚组，表明是一个近乎饱和的数据库。

NGS样本基因分型基准测试

研究对2,504名无关个体和641名来自1000基因组计划（1kGP）的后代进行了分型。大多数亚组（99.25%）显示哈迪-温伯格平衡，因此偏倚很小。有27个矩阵中超过15%的亚组处于不平衡状态，这些大多是短基因（中位数=4,564 bp），且具有少量低拷贝k-mer。基因型准确， trio 一致性平均F1分数为97.58%，而18个矩阵具有高不一致性（>15%），主要发生在亚端粒基因或性染色体上，这些区域组装质量较差。

研究评估了高重复基因家族（例如淀粉酶、NBPF、GOLGA和TBC1D3）的拷贝数准确性和偏倚。将分型得出的拷贝数与相应组装体的拷贝数进行了比较，使用了包含这些样本的数据库。为了限制错误组装序列带来的复合误差，排除了具有低置信度序列的样本。对于每个样本，在所有对应组装体拷贝数高（>10）的矩阵上进行了基准测试。拷贝数高度相关（ρ=0.996，皮尔逊相关），偏倚很小，0.2%的缺失拷贝（假阴性）和2.4%的额外拷贝（假阳性），可能来自组装体中未组装的基因。当测试扩展到所有分型基因时，仍然保持高度一致性（ρ>=0.996，皮尔逊相关）。

研究评估了39个HPRC基准测试基因组的基因分型等位基因与真实基因组组装体之间的序列相似性。每个样本使用完整数据库（全集）或排除其对应PA的数据库（留一法）进行分型。将分型的PA与相应的组装体PA进行匹配，排除内含子和诱饵以及非重复碱基少于1 kb的序列，并测量了分型等位基因与指定查询之间的相似性。进行了类似的分析，将数据库中每个组装体PA的最远邻视为正确的分型位点。由于数据库抽样或错误组装的错配，留一法实验中2.9%的PA和全集实验中1.0%的PA未与真实拷贝配对进行评估。对于全集，配对的PA在非重复区域每10 kb有0.36个错配，93.0%在非重复区域没有错配。留一法测试在非重复区域每10 kb有2.7个错配，比最优解（最远邻）每10 kb多1.2个错配；57.3%的等位基因没有错配，77.0%映射到最优解。留一法结果与原始PA的相似度比最接近的GRCh38基因高96.5%，后者每10 kb有79.3个错配。

为了在错误组装情况下分离错误来源，研究直接将留一法分型结果与基因PA的端粒到端粒组装体进行比较。样本基因型有11,627个正确匹配的亚组，599个（4.8%）错误分型到其他亚组，131个来自组装体独有的亚组（1.1%；参考外），127个假阳性（0.5% F1）和93个假阴性（0.4% F1），总F1错误为6.7%，拷贝数一致性为99.1%。这与 trio 不一致性相比，错误分型增加了3%。

计算需求足以满足生物样本库分析。在30倍覆盖度下对3,351个基因进行分型的平均运行时间为80.2分钟（样本预处理每1倍覆盖度1.0分钟，每个基因分型0.9秒），使用单个核心，内存约20 GB，支持并行处理。

研究将HLA、KIR和CYP2D6基因型与位点特异性方法T1K和Aldy进行了比较。对于31个HLA基因，ctyper在所有四个HLA命名字段上的F1分数为98.9%（与全集分析相比），留一法分析为86.3%，而T1K为70.8%。对于蛋白质编码产物（前两个字段），ctyper全集分析达到99.98%（拷贝数F1正确性99.9%），留一法分析为96.5%（拷贝数F1正确性99.5%），T1K为97.2%。对于14个KIR基因，ctyper全集分析在所有字段上达到98.5%，留一法分析为70.6%，而T1K由于数据库有限为32.0%。对于蛋白质编码产物（前三位数字），ctyper全集分析达到99.2%（拷贝数F1正确性99.9%），留一法分析为88.8%（拷贝数F1正确性99.2%），而T1K为79.6%。基准测试CYP2D6星号注释，ctyper全集分析达到100.0%，留一法分析为83.2%，而使用Aldy为80.0%。ctyper与全集分析在SNP变异上完全一致，留一法分析为95.7%，而使用Aldy为85.2%。

最后，研究使用ctyper对273个CMR基因进行分型。非重复区域与全集分析相比平均每10 kb有0.29个错配，比将组装体与相应GRCh38序列进行比较（基线）少99.7%。使用留一法数据库的基因型每10 kb有4.9个错配，比基线少94.8%。包括重复屏蔽的低复杂度序列（例如可变数目串联重复），与全集分析相比每10 kb有10.5个错配（比基线少97.6%），留一法分析每10 kb有74.7个错配（比基线少82.7%）。

研究将HLA和CMR基因的分型与使用泛基因组的当代方法Locityper进行了比较，使用留一法分析。对于HLA，Locityper在预测所有四个命名字段上的F1分数为87.9%（对比ctyper为86.3%），而ctyper在蛋白质编码变异的前两个字段上表现稍好（96.5%对比94.0%），尽管ctyper由于无比对齐分型速度提高了约218倍。当分析CMR基因型时，在可比区域中，比Locityper基因型每10 kb少19.8个错配。

全球人群中CNV的序列水平多样性

研究使用主成分（PC）分析（PCA）检查了2,504名无关1kGP样本、879名基因型-组织表达（GTEx）样本和105个二倍体组装体（排除HGSVC due to质量过滤）中PA基因型的群体结构，排除了罕见亚组（<0.05等位基因频率）并将拷贝数限制为十个以平衡PC的权重。所有数据按群体而非来源聚类，表明分型与组装体之间或跨NGS队列的偏倚很小。PC1中权重最高的前0.1%亚组的平均aggreCN方差为26.33，显著大于整体的4.00（P值=1.11×10^-16，F检验）。类似地，PC2和PC3的平均aggreCN方差为19.73和7.20，表明CNV与序列变异的关联较弱。此外，PC1是唯一将所有样本聚类到相同符号且地理中心远离0的PC，表明如果將样本视为paCNs向量，则它对应于模数方差（因此是aggreCN）。同时，PC2和PC3类似于基于全球样本SNP数据的PCA图，表明它们与CNV基因的序列多样性相关。非洲人群中的重复总数增加，反映在PC1的顺序中。

研究检查了ctyper基因型以衡量重复显示群体特异性的程度。使用F统计量（F检验的推广，适应多于两种基因型）来测试跨大陆群体分布的差异。总共4.4%（223 of 5,065）的重复亚组显示群体特异性（F统计量>0.2）。F统计量最高（0.48）的PA亚组包含HERC2P9的重复，这是一个已知的分化基因。此外，注释为SMN1重复的SMN2转化拷贝在非洲人群中富集（F统计量=0.43）。

研究然后测量了重复基因与其参考拷贝的分化程度，指示近期或古代重复，并提供来自缺失平行同源基因的参考偏倚的度量。为每个矩阵的序列构建了多序列比对（MSA），并测量了非重复序列中的所有 pairwise 差异。确定了平行同源序列分化相对于直系同源分化的平均相对平行同源分化（RPD）。还使用平均绝对误差（MAE）测量了拷贝数多样性，指示群体中的CNV水平。基于RPD，使用基于密度的噪声应用空间聚类，在0.71和3.2处识别出两个峰值，MAE中心在0.18和0.93，分别对应具有罕见和近期CNV的基因以及更多分化和常见CNV的基因。后者反映了不同结构单倍型上的CNV，无法使用单个参考基因组进行分析。例如，AMY1A由于截断重复而具有3.10的高RPD。这些结果与人类进化中古代爆发重复一致。

研究然后使用ctyper基因型调查不同CNV位点的重组。使用无关1kGP样本确定了989个在GRCh38中间隔小于100 kb的相邻亚组之间的多等位基因LD（mLD），并报告了每个矩阵内的平均mLD。拷贝数MAE与mLD之间存在更强的负等级相关（ρ=-0.24，P值=3.4×10^-15，斯皮尔曼等级），高于mLD与位点长度之间的等级相关（ρ=-0.21，P值=1.5×10^-11，斯皮尔曼等级），表明在频繁CNV的基因中单倍型连锁减少。最低mLD（0.013）发现在FAM90中，这是一个具有频繁重复和重排的基因。最高mLD（mLD>0.7）的29个位点在性染色体中富集（n=19）。此外，HLA-B和HLA-DRB的mLD>0.7且仅有缺失CNV。

eQTL分析

为了研究paCNs对表达的影响，研究使用遗传欧洲疾病变异（GEUVADIS）和GTEx队列进行了eQTL分析。有4,512个基因可以在RNA-seq比对中唯一映射。另外44个基因，如SMN1、SMN2、AMY1A、AMY1B和AMY1C，具有无法区分的转录产物，通过 pooling 所有拷贝进行分析。基于外显子序列将PA分配给这些转录本，并与paCNs进行了关联分析。

将paCNs合并为aggreCNs后，5.5%（178 of 3,224）的转录本显示显著性（校正P=1.6×10^-5，皮尔逊相关），如先前观察到的。相比之下，当用个体paCNs更新aggeCNs并对表达进行多变量线性回归时，27.6%（890 of 3,224）的转录本的拟合有显著改善（校正P=1.6×10^-5，单尾F检验）。为了测试拟合是否由同一参考基因的不同等位基因的非均匀表达解释，研究使用线性混合模型（LMM）将总表达回归到个体亚组，并估计等位基因特异性表达，然后将这些值与分配给同一参考基因的同一矩阵的其他亚组进行比较。对于可解矩阵内且样本数超过十个的亚组，研究发现7.94%（150 of 1,890）的平行同源基因和3.28%（546 of 16,628）的直系同源基因具有显著不同的表达水平（校正样本量=平行同源基因数+直系同源基因数，校正P=2.7×10^-6，χ²检验）。总体而言，发现平行同源基因表达减少，与重复基因的先前发现一致。

研究比较了GTEx样本中57个组织的表达，以测试平行同源基因的偏好表达。在2,820个平行同源基因中有132个（4.68%）和19,197个直系同源基因中有225个（1.17%）存在替代组织特异性（校正P=6.4×10^-8，两个χ²检验的并集）。

此外，研究使用方差分析（ANOVA）估计了paCNs在GEUVADIS表达数据中解释的表达方差比例（R²），并将其与基于已知SNP、 indel 和eQTL结构变异（SV）的模型进行了比较。正如预期，高度 granular 的paCNs解释了最多的方差：平均为10.3%（包括基线为14.3%）。相比之下，58.0%的转录本是具有已知eQTL变异的基因，这些变异解释了2.14%的有效方差（考虑实验噪声为1.60%，与先前估计的1.97%一致）。平均而言，aggreCNs解释了1.98%的方差，亚组信息解释了8.58%的方差。当结合paCNs和已知eQTL位点时，解释了10.4%（包括基线为19.0%）的有效方差。

研究检查了SMN和AMY2B基因作为案例研究， due to 它们在疾病和进化中的重要性。SMN基因被分类为SMN1、SMN2和SMN-转化。研究发现所有SMN1和SMN2转录本的总表达没有显著差异（0.281±0.008对比0.309±0.009；P=0.078，χ²检验）。然而，发现SMN-转化与SMN1和SMN2之间存在显著差异（0.226±0.012对比0.294±0.002；P=1.75×10^-7，χ²检验），SMN-转化表达减少23.0%。相比之下，尽管总体表达较低，但SMN-转化的有效外显子7剪接表达是SMN2的5.93倍（P=2.2×10^-16，χ²检验），表明SMN-转化具有全功能剪接但总体表达较低。研究了AMY2B重复的表达，包括易位到其他AMY基因近端的等位基因，例如图2a中包含AMY1和AMY2B的PA。使用概率估计表达残差（PEER）校正的GTEx胰腺数据，发现易位的AMY2B基因表达显著高于其他重复（1.384±0.233对比-0.275±0.183，P=7.87×10^-9，χ²检验）。

Discussion

新的泛基因组为研究复杂遗传变异（例如CNV、 recurrent SV、易位和基因转换）提供了机遇和挑战：它们揭示了复杂变异的景观，但需要新的工具进行表示和分析。研究将基因组变异表示为PA：捕获基因组结构信息和相位变异的单倍型片段。为了支持大型NGS队列分析，研究开发了无比对齐分型工具ctyper，用于使用NGS对PA进行分型，提供等位基因特异性序列信息和拷贝数。分型基于一个新的数学模型，该模型将NP难问题放松为更有效的多项式半解析解，具有稳健的基因型和拷贝数估计。尽管这里的分析侧重于CNV基因，但ctyper适用于对复杂遗传变异和局部相位进行全基因组分型。

ctyper基因型的使用扩大了NGS研究的范围，可以分析难以映射的CMR和CNV基因的变异。例如，关于CNV反映两种变异模式（高度相似（近期）和低一致性（古代和多态性）重复）的发现是基于1kGP基因型而非组装注释。此外，ctyper基因型产生了平行同源基因的组织特异性表达以及不同重复结构对表达的相对贡献，例如SMN基因。

研究调查了PA上ANOVA的显著改进，其基因型是多等位基因的，反映了与常规已知双等位基因eQTL变异不同的变异组合。首先，与PA相比，每个基因的eQTL变异要么非常少，要么非常多，表明存在LD，如精细映射所解决，并且在后者情况下，增加了多重检验负担。由于LD的间接关联也解释了为什么在具有更多CNV的基因中，传统eQTL变异解释的方差比例更大（例如HPR基因）。然而，随着CNV频率的增加，eQTL变异解释的方差增加（t=3.80，P值=1.6×10^-4，皮尔逊相关），而eQTL变异数量减少（t=-4.79，P值=2.1×10^-6，皮尔逊相关），表明更大的效应（如CNV）可能掩盖了不在LD中的其他eQTL变异的发现（总方差的增加降低了使用类高斯模型的关联分析中的显著性）。由于PA包含了LD，它们将较少受到此类基于LD的问题的影响。此外，基因表达可能不是变异的线性加性效应。例如，尽管SMN-转化包含来自SMN1和SMN2的变异，但其总体表达低于两者。通过这种方式，使用具有连锁变异的遗传模型（如PA）在线性加性模型的基础上改进了对基因表达的预测。由于这些限制也适用于非CNV基因，PA的概念可能在未来关联分析中具有更广泛的潜力。

由于样本量有限，研究的关联基于亚组而非个体PA。不同的队列规模可能需要不同水平的亚组粒度。当前的粒度是为生物样本库队列设计的；较小队列的研究可能需要定义更广泛的亚组。例如，SMN-转化的三种亚型在eQTL分析中显示出 little 差异，并在案例研究中合并，但更大的队列研究可能会发现它们的差异。分型的粒度还由PA的长度定义；较短的PA更准确地反映小变异，而较长的序列保留更多结构信息，可能在低重组区域（如HLA-DRB）中更可取。

Methods

构建泛基因组衍生等位基因数据库

从泛基因组组装中搜索和提取初始目标基因

泛基因组队列由来自HPRC（92个单倍型，排除由于大量 flagged 区域的HG02080）、CPC（114个单倍型）、HGSVC（18个单倍型；仅使用PacBio HiFi组装）、两个端粒到端粒二倍体组装（4个单倍型）和参考基因组（GRCh38包括替代位点和T2T-CHM13）的组装组成。使用的基因坐标来自基于GRCh38参考基因组的GENCODE版本39。

研究为在HPRC和CPC研究中发现具有拷贝数变异的3,203个基因构建了数据库。基因最初被组织成“查询集”，其中每个查询集包含具有功能或相似序列的基因，包括假基因和同一基因家族内具有远缘同源性的基因。查询集最初基于共享名称前缀的基因定义，并用于定位泛基因组中的所有相似序列。

高效的映射方法错过了包含k-mer匹配的序列的比对，当未包含在数据库中时降低了分型准确性，包括小假基因和分化平行同源基因。为了解决这个问题，研究开发了一种以k-mer簇为中心的敏感且高效的扫描方案，以检测泛基因组中目标基因的所有相似序列。

由k-mer定义的热点通常包括由查询集中多个基因映射的位点和串联重复基因。为了解释这种冗余，研究合并了相距小于10 kb的比对（加上5 kb侧翼序列，这合并了20 kb内的基因）。为了避免分型可能被重组分割的较长位点，研究在超过20 kb的内含子中点划分了位点。为确保最小位点长度，调整了侧翼序列以实现最小15 kb的长度。这些方法旨在标准化每个序列的大小以近似LD区块的大小。由