Nature | 基于四代家系参考的人类新生突变率

生信菜鸟团

发布于 2025-07-12 16:47:55

9280

文章被收录于专栏：生信菜鸟团生信菜鸟团

Basic Information

英文标题：Human de novo mutation rates from a four-generation pedigree reference
中文标题：基于四代家系参考的人类新生突变率
发表日期：23 April 2025
文章类型：Article
所属期刊：Nature
文章作者：David Porubsky | Evan E. Eichler
文章链接：https://www.nature.com/articles/s41586-025-08922-2

Abstract

Para_01

了解人类新发突变（DNM）率需要完整的序列信息。
在此，我们使用五种互补的短读长和长读长测序技术，对一个四代二十八人家系（CEPH 1463）中每个二倍体人类基因组超过95%的部分进行了分型和组装。
我们估计每次遗传过程中会出现98至206个DNM，其中包括74.5个新发单核苷酸变异、7.4个非串联重复插入缺失、65.3个来源于串联重复的新生插入缺失或结构变异，以及4.4个着丝粒DNM。
在男性个体中，我们发现每代平均有12.4个新发Y染色体事件。
短串联重复序列和可变数目串联重复序列是最容易发生突变的区域，其中有32个位点在几代之间表现出反复突变。
我们在多代样本中准确组装了288个着丝粒和六个Y染色体，并证明DNM率会因重复内容、长度和序列一致性而变化达一个数量级。
我们展示了所有类型生殖系DNM存在强烈的父系偏差（75-81%），但我们估计16%的新发单核苷酸变异起源于合子形成之后，没有父系偏差，包括早期生殖系嵌合突变。
我们将所有这些变异置于高分辨率重组图谱（约3.4 kb断点分辨率）背景下进行分析，并发现减数分裂交叉事件与新发结构变异之间没有相关性。
这些接近端粒到端粒的家族基因组提供了一个真实数据集，用于理解人类遗传变异背后最基本的过程。

Main

Para_01

人类基因组的端粒到端粒（T2T）组装增加了大约8%的高重复DNA区域，这些区域通常被排除在人类遗传变异研究之外，例如着丝粒、片段重复序列（SDs）以及近端着丝粒区域。
许多单倍型解析的人类二倍体基因组的长读长测序已经为突变机制的研究提供了新的见解，使得无论变异类型或复杂程度如何，都有可能被发现。
将父母基因组与其子代进行直接比较，能够增强识别新发突变的能力，而不是通过将测序数据比对到中间参考基因组（如GRCh38或T2T-CHM13）来实现

Para_02

本研究的目标是构建一个高质量的人类家系资源，实现染色体的完整组装和分型，并在代际间研究其遗传过程，以增进我们对重组和新生突变（DNM）过程的理解。
我们试图消除在发现过程中存在的三种选择偏倚，这些偏倚可能针对特定基因组区域、遗传变异类型以及参考基因组效应。
为实现这一目标，我们将重点放在了拥有四代、28名成员的CEPH 1463家系上。
在过去三十年中，该家系得到了深入研究，并且我们使用五种具有不同且互补性错误模式的测序技术对该家系成员进行了测序。
这个特殊的家系曾作为早期连锁图谱研究的基准模型，这些研究采用了短读长测序技术（SRS）并继续作为理解人类变异的重要参考。
该家系目前仍在用于研究人类变异，包括嵌合体现象的模式。

Para_03

正如最初的T2T基因组被用作理解基因组所有区域的参考，我们的目标是为遗传变异和新生变异都建立一个标准参考数据集。

Genome sequence and assembly

Para_01

我们为一个四代家庭（CEPH 1463 家系）中的大部分 28 名成员生成了 PacBio 高保真（HiFi）数据、超长读长 Oxford Nanopore Technologies（UL-ONT）数据、Strand-seq 数据、Illumina 数据以及 Element AVITI Biosciences（Element）的全基因组测序（WGS）数据（图 1 和补充表 1）。

Fig. 1: Sequencing the CEPH 1463 pedigree with five technologies.

- 图片说明

◉ 我们使用五种正交的下一代和长读长测序技术对四代系谱CEPH 1463中的28名成员进行了测序：对于G2至G4代，采用HiFi测序、Illumina和Element测序技术对其外周血进行了分析；而对于G1至G3代，则利用可用的永生化淋巴细胞系生成了UL-ONT和Strand-seq数据。该系谱数据集已被扩展，纳入了第四代以及G3代的配偶（编号为200080和200100）。◉ 该研究采用了多种先进的基因组测序平台，确保在不同世代之间获得全面且高精度的遗传信息，从而增强了对遗传变异传递模式的理解。

Para_02

为了发现变异，我们专注于从血液来源的DNA生成长读长PacBio、短读长Illumina和Element数据，以避免细胞系特异性的人为假象。
我们还利用相应的细胞系生成UL-ONT读长，用于构建接近端粒到端粒（near-T2T）的组装序列，以及生成Strand-seq数据来检测大的多态性倒位并评估组装准确性（方法和补充表2）。
简而言之，我们从多个正交测序平台生成了深度全基因组测序数据，主要关注前三个世代（G1–G3）（扩展数据图1a），并使用第四个世代（G4）来验证新生种系变异。
我们应用了两种混合基因组组装流程Verkko和hifiasm，用于为G1–G3生成高度连续且分相的基因组组装结果，而G4成员则仅使用HiFi数据进行组装（方法）。

Para_03

总之，Verkko 组装的连续性最强（AuN（类似于平均 contig 长度的度量）：102 Mb）（扩展数据图 1b、补充图 1 和 2 以及补充说明 8），我们估计在 G1–G3 中有 63.3%（504 条染色体中有 319 条）接近端到端完整（扩展数据图 1c）。
此外，在非近端着丝粒染色体中，有 42.3%（504 条中有 213 条）被单个 contig 覆盖，并且每端都具有规范的端粒重复序列（方法部分、扩展数据图 1d、补充图 3 和补充表 3）。
我们在 G1–G3 中测序并组装了 288 个着丝粒（占 44.7%，即 644 个中的 288 个），并且注意到不同的组装软件更倾向于组装不同的人类着丝粒（方法部分、扩展数据图 1e 和补充图 4）。
无论是序列质量（QV 范围为 47–58）还是分型准确性都很高（方法部分、补充图 5–12 和补充表 4）。

A multigenerational variant callset

Para_01

这一数据资源使我们能够追踪任意基因组片段及其相关变异在全部四代中的遗传情况（扩展数据图2a）。
我们在第二代和第三代中总共识别出595万个单核苷酸变异（SNVs）和插入缺失（indels），以及35,662个结构变异（SVs），这些变异均符合孟德尔遗传规律（方法，补充表5，补充图13及数据可用性）。
在595万个变异中，77%的小型变异得到了三种技术的支持，通过使用原始材料进行变异检测，有助于排除由细胞系人工产物引发的新生突变（补充说明1）。
与‘瓶中基因组’（GIAB）（2.51 Gb）或Illumina全基因组测序（WGS）（2.58 Gb）相比，长读长测序（LRS）让我们能够额外获取约2.6亿碱基对的人类基因组序列（共2.77 Gb），其中包含2.01亿碱基对是这两项研究中均未包含的。
在重复区域（SDs）及其相关基因中，我们观察到了一些最大的改进。
我们将85.5%的重复区域（覆盖度>95%，共计8,048个合并后的重复区域中有6,883个）分类为高置信度，而此前GIAB分析中仅有25.6%（共计8,048个中有2,060个）被分类为此等级，这标志着高度拷贝数可变区域的重大改进。
我们发现大部分已知拷贝数可变区域（>91%）在此家系中稳定传递，而剩余9%的区域通常被标记为可能组装错误（补充说明2）。
同样地，我们提供了移动元件插入的全面统计信息（方法，补充表6，补充图14和补充说明3），并识别出120个以孟德尔方式分离的倒位（其中21个存在歧义）（补充表7和补充图15-18）。
后者包括一个罕见的倒位（约703 kb），其覆盖了位于15号染色体q25.2区段的一个与疾病相关的拷贝数可变区域，以及一个位于16号染色体q11.2区段的倒转重复序列（约295 kb）（补充图19和补充图20）。

Sequence-resolved recombination map

Para_01

通过使用三种不同的方法13,22（方法和扩展数据图2b），我们鉴定了G3代（n = 8）相对于T2T-CHM13参考基因组的539个减数分裂断点，并且其中99.8%（538/539）得到了多于一种方法的支持（补充表8和补充图21）。
在初始分辨率达到约3.4 kb的基础上，通过父母与子代之间的直接基因组比较（方法和补充图22），我们将90.4%（487/539）的断点进一步精确定位到中位大小约为2.5 kb的区域。
值得注意的是，由于T2T-CHM13参考基因组中存在的参考偏差，其中有191个断点的实际大小有所增加（补充图23）。
我们将重组断点分为两类：一类是在父母单倍型之间具有非常陡峭过渡的断点，另一类是在两个父母单倍型上都存在较长同源区域的断点（扩展数据图2b和补充图24）。
此外，我们还鉴定出78个小的单倍型片段‘转换’事件（中位大小约为1 kb）23,24,25，这些事件可能与双交换重组或等位基因转化事件一致，但由于我们采用了严格的过滤标准，这个数字可能是低估的结果（方法、补充表9和补充图25）。
将重组定位分析扩展至G4代染色体后，我们又新增了964个断点，从而使22次遗传过程中共鉴定了1,503个减数分裂重组断点（补充图26）。
这些断点包括16个重组热点，其中11个与之前报道的重组率升高区域一致26（补充表8和补充图27）

Para_02

总体而言，15–20%的父系和母系同源染色体在没有可检测到的减数分裂断点的情况下被传递（即非重组染色体）（补充图28）。
我们观察到母系重组事件显著过剩（Wilcoxon符号秩检验，P = 6.4 × 10−5），预期的母系与父系断点比例为1.4（参考文献27）（补充图29）。
父系重组明显偏向人类染色体的末端区域，有55个父系重组事件位于距端粒2 Mb以内的区域，而女性个体中仅有1个这样的事件（扩展数据图2c和补充图30）。
在G2–G3代中，我们观察到随着父母年龄的增长，男性和女性生殖细胞中的交叉事件均有所减少（扩展数据图2d）。
我们在G1–G4代中使用具有对数连接函数的泊松广义线性模型（GLM）对此观察进行了建模，并继续观察到重组断点的数量随父母年龄增长及性别呈现出显著下降的趋势（父母年龄和性别的P值分别为7.17 × 10−3和1.22 × 10−9；AIC = 284.2）（补充图31）。
尽管目前尚无已知的生物学机制可以解释为何双亲生殖细胞中的重组事件均会减少，但这一观察结果与基于SRS数据进行的人群水平分析结果相悖。
我们认为该观察结果在更大数量的家庭被分析之前仍属初步发现。

De novo SNVs and small indels

Para_01

为了发现小型变异，我们检查了比对到T2T-CHM13的HiFi读段，然后使用正交的ONT和Illumina数据来确认某个变异确实存在于样本中且在父母中不存在（方法见下文）。
这种策略减少了平台偏差，但将新发突变（DNM）的发现限制在G2（n = 2）和G3（n = 8）代个体，因为G4代没有生成ONT数据。
我们的从头变异集合包括常染色体上的755个单核苷酸变异（SNVs）和73个插入缺失（indels）（图2a），以及X染色体上的27个SNVs和1个indel。
我们使用侧翼SNVs构建单倍型、确定变异的相位，并追溯突变来自父母配子还是早期胚胎。
如果一个突变满足以下两个标准之一，我们就判定该突变为体细胞发生，而且很可能是在胚胎早期发育阶段发生的：与父母单倍型不完全连锁（n = 122），或者如果无法进行相位分析，则其等位基因平衡在所有三个测序平台上均显著低于0.5（n = 7）（图2b）。
此外，我们通过追踪单倍型跨世代向后传播的情况，以及针对四个有后代测序数据的个体向前传播的情况，验证了每一个合子后突变（PZM）（补充说明4）。
这一结果进一步通过Element数据得到了确认（补充图32）。

Fig. 2: Summary of DNM rates.

- 图片说明

◉ a，CEPH 1463 家系中父母（G2）和八个子女中新发种系突变（PZMs）以及小片段插入缺失（<50 bp）的数量。由于第三代（G3）父母的测序深度更高且可以评估突变传递情况，因此仅在 G3 中展示了串联重复序列区域的新发突变（TR DNMs）（方法详见原文）。图中阴影柱状图表示被确认传递到下一代的单核苷酸变异（SNVs）数量。◉ b，种系 SNVs（n = 626）在不同测序平台上的平均等位基因平衡接近 0.50，而合子后 SNVs（n = 119）的平均等位基因平衡低于 0.25。箱型图显示了中位数（中心线）、四分位距（IQR，箱体边界），触须延伸至 25% − 1.5 × IQR 和 75% + 1.5 × IQR；异常值以点的形式显示。◉ c，对于种系中新发的 SNVs，观察到了显著的父亲年龄效应（每年增加 1.55 个新突变；双侧 t 检验，P = 0.013），但在 PZMs 中未观察到这种效应（P = 0.72）。对于 DNMs 或 PZMs，母亲年龄效应均不显著（DNMs：每年增加 0.20 个新突变，P = 0.54；PZMs：P = 0.74）。实线为使用线性模型拟合的回归线；周围阴影区域代表其 95% 置信区间。◉ d，根据基因组区域估计的 SNV 新发突变率表明，在大型重复区域（包括着丝粒和复制区段 SDs）中新发突变明显富集。基于组装的方法对着丝粒和 Y 染色体上的新发突变进行检测，结果显示卫星 DNA 区域中新发突变增加。使用双侧 t 检验判断这些区域与常染色体 DNMs 或 PZMs 突变率之间是否存在显著差异；*P < 0.05，**P < 0.001。各比较对应的 P 值分别为：0.0066（SDs 中基于比对的 DNMs）、0.049（SDs 中基于比对的 PZMs）、0.017（着丝粒中基于比对的 DNMs）、0.34（着丝粒中基于比对的 PZMs）、0.13（着丝粒侧翼区域中基于组装的 DNMs）、0.14（着丝粒 HORs 区域中基于组装的 DNMs）、0.59（Y 染色体常染色质区域中基于组装的 DNMs）和 0.00025（Yq12 区域中基于组装的 DNMs）。

Para_02

在这四个样本的62个PZM中，有64.5%（n = 40）被传递到下一代，相比之下，97.2%的生殖系SNV（249个中有242个）和100%的插入缺失（indel）被传递（扩展数据图3）。
我们发现有10个PZM未能通过基于单倍型的验证，最终得到119个PZM的集合，占全部常染色体SNV总数的16%（745个新生SNV）。
此前基于Illumina平台对该家族的研究识别出605个新生SNV，这些变异来源于生殖细胞（G2和G3）或合子后突变（仅G2），其中92.4%（n = 559）被包含在我们的最终变异集合中，而未包含在内的变异中除4个外均未能通过长读长数据验证。
我们首次在G3中鉴定出另外72个PZM，其中包括总共186个新的DNM，在生殖系SNV和插入缺失的检出率上分别提高了6.1%和21%。

Para_03

总体而言，81.4%的种系小DNMs起源于父系单倍型（父系：母系比例为4.38:1，Wilcoxon符号秩检验，P < 2 × 10⁻¹⁶）。
通过线性回归拟合发现，父亲年龄每增加一岁，种系DNMs数量平均增加1.55个（双侧t检验，P = 0.013）。
相比之下，PZMs在亲本来源上没有显著差异（父系：母系比例为1.38:1，Wilcoxon符号秩检验，P = 0.09），也没有观察到父母年龄的影响（图2c）。
尽管我们的样本量较小，不足以检测新发突变与合子后突变谱之间的显著差异（补充图33a）。
但我们确实观察到一种新的现象，即CpG>TpG类型的PZMs有所减少（χ²检验，P = 0.17），以及之前曾观察到的合子后T>A替换增加的趋势（χ²检验，P = 0.268）。

Para_04

我们成功检测了91.9%的常染色体基因组（2.66 Gb）（补充图33b和补充说明4）。
排除所有被归类为合子后变异的位点，我们发现父母生殖细胞系每代每个碱基对的新生单核苷酸变异（SNV）率为1.17 × 10−8（95%置信区间（CI）= 1.08 × 10−8–1.27 × 10−8）。
新生SNVs显著富集在重复序列中，在着丝粒区域最多可达到2.8倍（每代每碱基对的SNV率：95% CI = 1.79 × 10−8–5.51 × 10−8，双侧t检验，P = 0.017）。
在SD区域（segmental duplication，节段性重复区域），该富集程度为1.9倍（95% CI = 1.64 × 10−8–2.88 × 10−8 SNVs每碱基对每代，双侧t检验，P = 0.0066）（见图2d，补充图33c和补充表10）。
我们在整个常染色体范围内观察到较低的合子后突变（PZM）率，为2.04 × 10−9 SNVs每碱基对每代（95% CI = 1.68 × 10−9–2.47 × 10−9）。
然而，在SD区域中，PZMs的富集程度达到了3.9倍（95% CI = 4.84 × 10−9–1.25 × 10−8 SNVs每碱基对每代，双侧t检验，P = 0.049）。
在传递到下一代的PZMs中（共检测到4个样本中的33个PZMs），我们在SD区域观察到2.69倍的富集（95% CI = 1.15 × 10−9–1.08 × 10−8 SNVs每碱基对每代）。
但由于样本量较小，这一富集效应未达到统计学显著水平（双侧t检验，P = 0.218）。

De novo TRs

Para_01

在此，我们研究了串联重复序列（TRs），包括短串联重复序列（STRs，1–6 bp的基序）和可变数目串联重复序列（VNTRs，7–1,000 bp的基序）。
我们使用串联重复序列基因分型工具（TRGT）对家系中所有成员的HiFi数据中的782万个TR位点中的768万个成功进行了基因分型（方法）。
其中，有717万个（93.4%）位点在所有三联体中完全符合孟德尔遗传规律。
我们使用TRGT-denovo来识别候选新生突变（DNMs），这些位点在一个三联体的所有成员中至少被10条HiFi读长覆盖；平均而言，有688万个TR位点满足这一标准。
我们通过正交测序和传递分析对这些初步识别的DNMs进行了进一步验证（方法）。
由血液DNA生成的Element测序在同聚物区域后的错误率显著较低，因此我们测试了它是否能够更准确地测量同聚物和其他TR等位基因的长度。
我们在Element数据中观察到同聚物区域的影子峰（stutter）很低；在TRGT鉴定出的1000个纯合同聚物位点的随机样本中，平均有99.5%的Element读长能完美支持GRCh38中TRGT基因分型的等位基因大小，相比之下，Illumina测序读长的支持率为93.5%（补充图34和35）。

Para_02

我们使用Element数据进一步验证了TRGT-denovo调用的新发TR等位基因。
由于Element数据的读长较短，我们仅能评估613个新发STR等位基因中的80个（占13.1%）（每个样本平均评估10个STR）。
如果Element读数支持子代中TRGT等位基因长度，而在双亲中均不支持，则我们认为该DNM得到了验证（允许存在±1个碱基对的误差；见方法部分）。
在我们可以评估的80个新发STR中，有56个（70%）通过了我们严格的前后一致标准。
在均聚物区域的验证率（20个中有3个，占15%）低于非均聚物区域（60个中有53个，占88.3%）。
这表明我们对均聚物区域突变率的估计可能不够精确。
利用系谱信息，我们要求在两个具有测序子代的G3个体（NA12879和NA12886）中观察到的新发TR候选等位基因必须至少传递给后续一代（G4）中的一个子代。
在两个G3个体中观察到的128个新发TR等位基因中，有96个（75%）被传递至下一代。
这一比例显著低于新发SNV，反映出在准确表征新发TR方面仍存在的挑战。

Para_03

在完成元件和传递验证后，我们发现每个样本平均有65.3个串联重复从头突变（TR DNMs，包括短串联重复序列STRs、可变数目串联重复序列VNTRs和复杂位点）。
估计每个单倍体每代每个位点的TR DNM发生率为4.74 × 10⁻⁶（95%置信区间CI = 4.06 × 10⁻⁶–5.43 × 10⁻⁶），不同重复基序大小之间存在显著差异（图3a）。
总体而言，TR DNMs在每个样本中插入或删除了平均978个碱基对，每次事件平均为15.0个碱基对（补充表10）。
平均有54.9个突变为STR基序的扩增或收缩，2.6个影响VNTR基序，7.8个影响包含STR和VNTR基序的‘复杂’位点。
STR的突变率为每个单倍体每代每个位点5.50 × 10⁻⁶ DNMs（95% CI = 5.0 × 10⁻⁶–6.04 × 10⁻⁶）。
VNTR的突变率为0.83 × 10⁻⁶（95% CI = 0.51 × 10⁻⁶–1.27 × 10⁻⁶），主要由在SRS研究中无法评估的位点组成。
此前几项关于全基因组STR突变率的估算仅考虑了多态性STR位点；当我们仅分析CEPH 1463系谱中具有多态性的STR位点时。
我们发现每个单倍体每代每个位点有5.98 × 10⁻⁵个新的STR突变事件（95% CI = 5.43–6.57 × 10⁻⁵），这与之前估算的4.95 × 10⁻⁵–5.6 × 10⁻⁵大致一致。
总体而言，75.0%的相位明确的新发TR等位基因来源于父亲（图3b）。
二核苷酸基序的突变率高于同聚物，并且对于长度大于6个碱基对的基序，我们观察到随着基序大小的增加突变率也上升（图3a）。
正如先前研究所报道的，较大的TR位点（定义为参考基因组序列中TR位点的总长度）表现出更高的突变率（补充图36a）。
我们未观察到扩增或收缩的显著偏向性（双侧二项检验，P = 0.19）（补充图36b）。

Fig. 3: TR DNMs show motif-size-dependent mutation rates, paternal bias and are highly recurrent at specific loci.

- 图片说明

◉ a，TR DNM率（每单倍型每基因座每世代的突变数）根据T2T-CHM13参考基因组中每个TR基因座观察到的最小基序大小（n = 522），按照每个TR类别（STR、VNTR或复合型）进行展示（蓝色；左侧y轴）。每个个体中通过过滤标准的每个基序大小的平均基因座数量显示为灰色（右侧y轴）。◉ 误差条表示围绕平均突变率估计值的95%泊松置信区间（使用χ²分布计算）。◉ 突变率包括所有通过TRGT-denovo过滤标准和Element一致性分析的非重复性突变调用。◉ b，在G3代中对于可信定相的TR新发突变（DNMs）推断其亲本来源。斜纹填充表示至少传递给一个G4子代个体的情况（如存在）。◉ c，8号染色体上一个反复出现的VNTR位点的系谱概览：2376919–2377075（T2T-CHM13），其基序组成为GAGGCGCCAGGAGAGAGCGCT(n)ACGGG(n)。◉ 等位基因颜色表示由遗传向量确定的遗传模式，灰色表示数据不可用。◉ 符号表示相对于遗传自父母的等位基因的遗传类型：加号（+）表示新发扩增，减号（−）表示新发收缩，仅针对发生突变的等位基因显示；数字表示等位基因长度（以bp为单位）。◉ 在G3代8个个体中的7个中检测到新的TR等位基因，并传递给了4个G4代个体，其中两个在传递后进一步发生了扩增。◉ 一名G3代个体（200080）的配偶携带一种不同的TR等位基因，在后续传递中发生了新发收缩。◉ d，来自个体测序读段的垂直线代表c图中新发突变（DNM）的读段级别证据，按样本分别展示。◉ 如存在，同时展示了HiFi（顶部）和ONT（底部）测序读段。◉ 颜色与c图中的遗传模式一致；带有加号或减号标记的边框框突出显示了DNMs。

Para_04

我们识别出了一组TR位点，这些位点在系谱成员中反复发生突变。
我们确定了一个高可信度的包含32个位点的数据集（方法和补充表11）：其中五个位点显示代内重复出现（在至少两个G3个体中观察到DNMs），另外27个位点显示跨代重复出现（在至少两代中观察到DNMs）。
由于这五个代内出现的DNMs仅出现在单一世代中，它们可能代表亲代生殖系嵌合现象，而非重复发生的突变事件。
值得注意的是，在表现出重复现象的16个位点上，我们观察到了三种或更多不同的新生扩展或收缩事件（扩展数据表1）。
例如，我们重点介绍了一个跨代反复突变的TR位点，该位点出现了十个独特的新生扩展和/或收缩事件（图3c,d）。
所有等位基因的传递完全符合遗传向量（补充说明5），并且得到了HiFi和ONT读段的支持。

Centromere transmission and de novo SVs

Para_01

在288个完全组装的着丝粒中，我们评估了150次代际间的传递（图4a）。
通过比较父母与子女之间组装的着丝粒，我们鉴定出18个（占12%）经ONT和HiFi数据共同验证的新生结构变异（SV），插入和缺失的数量大致相当（图4b及方法部分）。
所有在这个研究中对子代进行了测序的新生SV（n=8）均确认被传递到了下一代（补充表10）。
我们发现，72.2%（18个中有13个）的SV映射到α-卫星高阶重复序列（HOR）阵列（扩展数据图4a），其余（18个中有5个，占27.8%）对应于各种着丝粒周围侧翼序列，但不包括单体α-卫星的侧翼区域。
所有α-卫星HOR新生SV事件都涉及每个着丝粒特异的基本α-卫星HOR模块的整数变化，其大小范围从680 bp（第9号染色体上的一个4聚体α-卫星HOR）到12,228 bp（第6号染色体上的四个18聚体α-卫星HOR）不等（图4c和扩展数据图4b）。
一次来自第9号染色体的传递在一个G2到G3代的传递过程中同时涉及了2,052 bp（六个二聚体α-卫星HOR单元）的增加和1,710 bp（一个4聚体α-卫星HOR和三个α-卫星二聚体单元）的丢失（图4d–f）。
第6号染色体的着丝粒具有最多的反复出现的结构变异事件，在三代中观察到了三次（图4a）。
第6号染色体的着丝粒拥有最多数量的近乎完全相同（>99.9%）的α-卫星HOR（扩展数据图4c）。

Fig. 4: De novo SVs among centromeres transmitted across generations.

- 图片说明

◉ a，正确组装的着丝粒数量（深灰色）以及传递到下一代的着丝粒数量（浅灰色）的总结。携带新生缺失、插入或两者兼有的已传递着丝粒用颜色标出。◉ b，在α卫星HOR阵列及其侧翼区域中新生结构变异的长度。◉ c，在G2-NA12878个体中6号染色体α卫星HOR阵列存在一个新生缺失，并在G3-NA12887个体中被遗传的例子。每个单倍型上的红色箭头表示α卫星HOR结构，单倍型之间的灰色区块表示共线区域。缺失区域由红色轮廓突出显示。Mat.表示母源；pat.表示父源。◉ d，G3-NA12885个体19号染色体α卫星HOR阵列中新发插入和缺失的一个例子。◉ e,f，对d图中插入（蓝色轮廓；e）和缺失（红色轮廓；f）的α卫星HOR结构的放大展示。每个单倍型顶部的彩色箭头表示α卫星HOR结构。◉ g，在G2-NA12877个体21号染色体着丝粒中存在两个新生缺失的例子。这些缺失位于着丝粒α卫星HOR阵列中的一个低甲基化区域，被称为CDR（着丝粒识别区域），被认为是动粒组装的位置。在CDR区域内三个α卫星HOR的缺失导致了G2-NA12877个体中CDR位置约260 kb的偏移。

Para_02

我们还评估了18个结构变异（SV）事件对与着丝粒凹陷区域（CDR）相关的低甲基化口袋的潜在影响——这是动粒附着位点的一个标志（方法部分）。
我们发现，有11个位于CDR区域之外的结构变异在世代之间对CDR中心点的变化影响较小（小于100 kb）（扩展数据图4d,e）。
而那些映射到CDR内部的结构变异则表现出更为显著的影响（平均偏移约260 kb），或者完全改变了CDR的分布（图4g和扩展数据图4f,g）。
尽管需要后续的CENP-A染色质免疫沉淀测序实验来确认动粒的实际结合位点，但这些发现表明结构突变可能会在改变动粒位置方面产生表观遗传学后果。

Para_03

最后，通过分析31个着丝粒的亲子传递（共计150.5 Mb），我们在着丝粒中识别出16个单核苷酸变异（SNV）的新生突变（DNMs），其中包括五个位于α-卫星高阶重复序列（HOR）阵列内的突变，对应的DNM率为每代每碱基对1.01 × 10−7（95%置信区间为5.75 × 10−8至1.63 × 10−7）。
该突变率与我们基于测序读段的作图方法所获得的结果相当，在超过10倍的序列数据中识别出了14个着丝粒SNV，得出的DNM率为每代每碱基对3.27 × 10−8（95%置信区间为1.79 × 10−8至5.51 × 10−8）（见图2d和补充表10）。
通过合并这两部分数据，我们估算出着丝粒的SNV DNM率显著提高，达到4.94 × 10−8（双侧t检验，P = 0.017）。
我们认为这一估计是保守的，因为所有突变事件都需要通过ONT和HiFi两种测序平台共同验证。

Y chromosome mutations

Para_01

在四代中有九名男性成员携带R1b1a-Z302 Y单倍群（图5a和补充表12）。
我们使用曾祖父（G1-NA12889；图1）的Y染色体组合作为参考，用于检测48.8 Mb的男性特异性Y染色体区域（MSY）中的新生突变（DNM）（方法和补充说明6）。
与基于HiFi读长的检测相比，基于从头组装的方法使可检测的碱基对数量增加了两倍以上，而新发现的新生单核苷酸变异（SNV）的数量则增加了七倍以上。
总共我们在5个G2–G3男性个体的MSY中识别出48个新生SNVs，每次Y染色体传递中发现的SNVs数量从7到11个不等（平均值为9.6，中位数为10）（补充表13）。
其中仅有2个SNVs位于Y染色体的常染色质区域，1个位于着丝粒周围区域，其余48个中的45个位于Yq12异染色质卫星区域（图5b）。
总体而言，我们估计整个MSY区域的新生SNV率为每代每碱基1.99 × 10−7突变（95%置信区间= 1.59 × 10−7–2.39 × 10−7）。
由于能够检测到Yq12卫星DNA，该估计值比先前报道的Y染色体常染色质区域的SNV率高出一个数量级（补充表13）。
我们注意到，在45个DNMs中有13个（29%）在Yq12区域的其他位置具有100%相同的匹配序列（但不在正交位置），这可能是DYZ1/DYZ2重复序列内的跨位点基因转换事件所致（方法）。
此外，我们总共识别出九个新生插入缺失（<50 bp，排除同聚物），每个样本中检测到1–3个插入缺失（平均每条Y染色体传递1.8个事件）以及五个新生结构变异（SVs，≥50 bp）（图5b和补充表13）。
这些SVs的大小范围从2,416到4,839 bp不等，每个都影响了整个DYZ2重复单元，平均每条Y染色体传递一个SV。
所有可验证的DNMs（SNVs，48个中的20个；插入缺失，9个中的6个；SVs，5个中的4个）均与预期的世代遗传结果一致（即从G2到G3–G4，以及从G3-NA12866到其三个G4男性后代）（图5b）。
总体而言，在Y染色体上识别的62个DNMs中，有82%（共51个）位于短读长无法可靠比对的区域（比对质量=0）（48个SNVs中的42个，9个插入缺失中的4个，5个SVs中的全部5个）

Fig. 5: Chromosome Y and an example of a de novo mobile element.

- 图片说明

◉ a，携带R1b1a-Z302 Y染色体的九个男性个体的系谱（左图）以及Y染色体组装的成对比较：来自HG00731（R1b1a-Z225）的密切相关Y染色体与三个世代中组装最连续的R1b1a-Z302 Y染色体进行比较。以每100 kb为一个区间，显示了不同样本之间的Y染色体序列一致性类别，并在系谱男性个体中鉴定出的质量控制通过的结构变异以蓝色和红色轮廓表示。◉ b，Y染色体从头突变（DNMs）汇总。顶部为G1-NA12889的Y染色体结构；在其下方展示了G1至G3各代Y染色体组装中鉴定出的所有DNMs；底部按突变类型和样本分别进行了展示。灰色部分展示了从G2传到G3至G4、以及从G3-NA12886传到其G4男性后代中的DNMs。◉ c，在G3-NA12887中发现的一个新生SVA插入事件。◉ d，G3-NA12887中新生SVA插入事件的HiFi读段支持情况。

De novo SVs

Para_01

我们总共在八名个体（G3）中验证了41个新生结构变异（SV），其中包括16个插入和25个缺失（方法），其中68%（41个中有28个）来源于父系生殖系，并且随着父亲年龄的增长，结构变异有增加的趋势（补充图37）。
几乎所有的结构变异（41个中有40个）对应于串联重复序列（TRs），包括着丝粒中的突变、Y染色体卫星序列和成簇的片段重复（SDs）（补充表10）。
我们估计每次遗传过程中大约有5个结构变异（95%置信区间=3–7），影响约4.4 kb的DNA（中位数，4,875 bp）。
如果我们排除映射到着丝粒和Y染色体的新生结构变异（n = 14），事件的中位大小会下降一个数量级（中位数，362 bp）。
非等位同源重组（NAHR）经常被用来解释串联重复序列的扩增和收缩机制42,43。
然而，我们发现27个常染色质区域的新生结构变异中没有一个与重组交叉事件重合（补充图37e）。
这表明减数分裂I期间同源染色体之间的非等位同源重组并非这些变异的主要发生机制，尽管我们不能排除其他与不涉及重组的双链断裂相关的机制。
我们鉴定出一个逆转录转座事件：一段完整的（3,407 bp长）SVA元件（SVAF亚家族）新生插入（G3-NA12887）44，其预测的供体位于上游约23 Mb处（图5c和补充图38）。
这个插入在一个亲本（G2-NA12878）中以较低频率存在（约占读段的11%），但在祖代传递的单倍型中不存在，这与该变异是在G2个体形成合子后出现的生殖系嵌合事件一致（图5d和补充图39）。

Discussion

Para_01

大多数DNM研究40,45,46,47,48,49基于来自大量三人组的SRS数据，得出每代人大约有60至70个DNM的结果；然而，这些研究通常排除了基因组中高度可变的区域。
我们基于多种平台和多代系谱的组装方法能够覆盖一些最重复的区域，例如着丝粒和Y染色体上的异染色质区域。
除了使用标准参考序列外，还利用父母本参考序列，并且具备在后续几代中确认遗传的能力，从而提高了检测的灵敏度和特异性。
在这个多代家系中，我们估计每次遗传事件中DNM的数量范围为98至206个（平均每代152个），并观察到显著的父系新生突变偏向（占70–80%）。
此外，随着父亲年龄的增长，这种突变数量会增加，这一现象不仅适用于单核苷酸变异（SNVs），也适用于插入缺失（indels）、结构变异（SVs），包括串联重复序列（TRs）。

Para_02

新生单核苷酸变异（de novo SNVs）的发生率根据基因组背景的不同存在超过一个数量级的差异，这与近期基于人群的人类分析结果7,50以及理论预测结果51一致。
串联重复区域（SD区域）的突变率增加了88%（分别为2.2 × 10−8（95%置信区间 = 1.64 × 10−8–2.88 × 10−8）和1.17 × 10−8）。
这种增加主要由具有>95%序列相似性的SD区域所驱动。
我们还观察到新生转换/颠换比例相较于整个基因组显著降低（卡方检验，P = 0.0109），这一结果与之前的预测结果7一致（见补充说明7）。
我们估计，Y染色体长臂1区2带（Yq12）异染色质区域41,52中的卫星DNA的突变率至少是常染色体常染色质区域的30倍（每代每碱基对突变率为3.86 × 10−7）。
该区域由数千个短串联重复序列组成（DYZ1/Hsat3A6 和 DYZ2/Hsat1B），这些序列组织成兆碱基大小的片段，并且具有>98%的序列一致性41,52。
此外，有29%的突变事件与Yq12区域中非直系同源位点相匹配，
这一现象与“位点间基因转换”（interlocus gene conversion）导致的超过20倍的突变率升高相符，可能源于姐妹染色单体交换事件的增加41。

Para_03

之前的研究预测，有6%-10%的DNMs并非来源于种系，而是在受精后某个时间点出现的，从而形成了嵌合变体。
这种区分是基于等位基因平衡阈值或跨越三代的与邻近SNV不完全连锁的依据。
长读长测序（LRS）通过将近乎所有新生单核苷酸变异（de novo SNV）分配到亲本单倍型来提高灵敏度，并通过其与该单倍型的不完全连锁来定义合子后变异（PZM）。
我们将16%的de novo SNVs归类为合子后起源（n = 119 PZMs/745 de novo SNVs）。
由于本研究中所有的测序数据均来源于血液，我们无法证明每个PZM都存在于多种组织中，但我们可以使用向下一代传递的结果作为替代指标，因为这表明突变也存在于生殖细胞中。
PZMs占传给下一代的所有SNVs的12%（n = 33 PZMs/275个传递的SNVs），比之前的估计有所增加。
人类胚胎的早期细胞分裂通常容易出错，快速的细胞分裂速率可能有助于高等位基因平衡（>25%）的PZMs占比显著。
这类事件以前会被归类为种系变异，但与PZM预期一致的是，我们发现这些DNMs没有父系偏向（图2c）。

Para_04

TRs 是我们基因组中最易突变的位点之一，这类新发突变事件的数量可与种系 SNV 相当，但每代影响的碱基对数量要高出一个数量级。
我们发现 TR 的新生突变率（DNM）随着重复次数和基序长度的增加而呈现出三倍的差异，这通常与突变率相关。
然而，我们在二核苷酸和更长的基序长度（>10 bp）之间观察到了明显的突变率低谷（图 3b），这可能反映了基于位点大小、基序长度和复杂性的不同突变机制。
例如，较大的 TR 基序可能更容易通过非等位同源重组（NAHR）、合成依赖链退火或位点间基因转换发生突变，而 STR 的突变事件则可能偏向于传统的复制滑动突变机制。
与一些早期对小卫星序列的全基因组分析结果一致，我们并未找到证据表明 TR 的变化是通过减数分裂过程中同源染色体之间的不等交换造成的，因为我们所有的 TR 新生 SV（n = 27）均未与重组断点重合。
在这方面特别值得关注的是发现了 32 个反复突变的 TR——这些位点很少在不稳定疾病等位基因之外的背景下被发现。
在其中五个这种反复突变的位点上，我们在单一代（G3）中发现了多个新生突变（DNMs）；这些 DNMs 可能是由于 G2 代父母生殖细胞嵌合现象或高度可变 TR 活性所致。
几乎所有这些高度反复的新发事件都产生了显著长于平均短读段长度的 TR 等位基因，只有使用长读长测序（LRS）才能检测到。
其中包括约 7% 的人类着丝粒长度变化，其插入和缺失均以主要 HOR 单位的整数倍形式发生。
新生结构变异（SV）的发生率从此前估计的每代 0.2–0.3 提高到了本研究报道的每代 3–4 个新生 SV。

Para_05

本研究存在几项局限性。
首先，即使使用Element数据，同聚物仍然具有挑战性，因为较长的等位基因和嵌入在更大重复序列中的模体仍无法通过短读长测序可靠地检测。
其次，由于端着丝粒区域的重复特性以及频繁发生的异常重组，我们无法对这些区域中的新发突变（DNMs）进行准确表征。
第三，我们将新发突变的发现仅限定在一个多代家庭的前三代，并使用第四代（G4）对已传递变异进行验证。
我们认识到，家系中的遗传变异依赖于遗传背景，因此需要更多的家系来建立可靠的突变率估计值，尤其是在基因组复杂区域。
在这方面，值得注意的是，目前已有相关工作正在进行，旨在对更多家系谱系进行表征。
尽管如此，本研究强调，仅依靠一种测序技术和一个单一的人类基因组参考序列，无法全面估计突变率。
此类多代资源将进一步优化新发突变率的估算，并为新的算法和测序技术提供另一个有价值的基准。

Methods

Ethics declarations

伦理声明

Human participants

人类参与者

Para_01

已获得CEPH/Utah个体的知情同意，且该研究已得到犹他大学机构审查委员会的批准（犹他大学IRB编号为IRB_00065564）。
其中包括23名家庭成员对研究成果数据发表的知情同意；其余5人提供了用于生物样本库保存的知情同意，并实施受控访问（数据可用性）。

Cell lines

细胞系

Para_01

CEPH 1463家族的14个成员（G1-GM12889、G1-GM12890、G1-GM12891、G1-GM12892、G2-GM12877、G2-GM12878、G3-GM12879、G3-GM12881、G3-GM12882、G3-GM12883、G3-GM12884、G3-GM12885、G3-GM12886 和 G3-GM12887）的细胞系来自科里尔医学研究所（CEPH细胞库）。
G3代配偶及G4代家庭成员（n=13）的细胞系为本实验室构建的EB病毒转化的淋巴母细胞系，包括：G3-200080-配偶、G4-200081、G4-200082、G4-200084、G4-200085、G4-200086、G4-200087、G3-200100-配偶、G4-200101、G4-200102、G4-200103、G4-200104 和 G4-200106。

Para_02

所有细胞系均通过DNA的全基因组测序及后续变异检测进行鉴定，以确认与个体预期性别相符，并与来自同一受试者血液DNA的测序结果一致。
此外，我们还探讨了所获得的测序数据是否符合亲子之间预期的遗传模式。
据我们所知，上述提及的细胞系均未进行支原体污染检测。

Sample and DNA preparation

样本与DNA制备

Para_01

G2和G3代的家庭成员被重新纳入研究，目的是更新知情同意书和健康史，并为其子女（G4代）和新加入的配偶父母（G3代）进行登记。
从全血中提取了来自G2和G3代的存档DNA。
新招募的家庭成员进行了知情同意程序，并采集血液以获取DNA和细胞系。
使用Flexigene系统（Qiagen 51206）从全血中提取DNA。
所有样本均广泛同意用于科研目的，这使得该数据集成为未来工具开发和基准研究的理想选择。

Sequence data generation

序列数据生成

Para_01

正交的短读长和长读长测序平台生成了如下测序数据：

Illumina data generation

Illumina 数据生成

Para_01

G1至G3的Illumina全基因组测序数据是如先前所述生成的。
G4以及G3配偶的Illumina全基因组测序数据由西北基因组中心生成，使用的是无PCR的TruSeq文库制备试剂盒，并在NovaSeq 6000平台上进行测序，双端读长为150 bp，测序深度约为30×。

PacBio HiFi sequencing

PacBio HiFi 测序

Para_01

PacBio HiFi 数据的生成遵循制造商的推荐方案。
简而言之，使用 Monarch HMW DNA 提取试剂盒（New England Biolabs，T3050L）从血液样本或培养的淋巴母细胞中提取 DNA。
在所有步骤中，使用 Qubit dsDNA HS（Thermo Fisher Scientific，Q32854）在 DS-11 FX（Denovix）上进行定量，并使用 FEMTO Pulse（Agilent，M5330AA 和 FP-1002-0275）检查片段大小分布。
HMW DNA 使用 Megaruptor 3（Diagenode，B06010003 & E07010003）系统进行打断，根据初始质量检测结果选择参数 28/30、28/31 或 27/29，以获得约 22 kb 的目标峰值大小。
打断后的 DNA 使用 SMRTbell 制备试剂盒 3.0（PacBio，102-182-700）生成 PacBio HiFi 文库。
根据实验方案，使用稀释的 AMPure PB 磁珠进行片段大小选择，或者使用 Pippin HT（Sage Science，HTP0001 和 HPE7510）进行高通量截留，截留范围基于打断大小设定为 10–17 kb。
文库在 Sequel II 平台上使用 Sequel II 测序化学试剂 3.2（PacBio，102-333-300）和 SMRT Cells 8M（PacBio，101-389-001），预延伸 2 小时，测序运行时间 30 小时，在 SMRT Link v.11.0 或 11.1 上完成；或者在 Revio 平台上使用 Revio SMRT Cell（PacBio，102-202-200）和 Revio 聚合酶试剂盒 v1（PacBio，102-817-600），预延伸 2 小时，测序运行时间 24 小时，在 SMRT Link v.12.0 上完成。

ONT sequencing

ONT测序

Para_01

为了生成超过100 kb的超长测序读段（UL），我们使用了牛津纳米孔技术（ONT）进行测序。
根据先前发表的方法，我们从淋巴母细胞系中提取了超高分子量的基因组DNA。
简而言之，将3–5 × 10⁷个细胞在含有10 mM Tris-Cl（pH 8.0）、0.1 M EDTA（pH 8.0）、0.5%（w/v）SDS和20 mg/ml RNase A的缓冲液中于37°C孵育1小时以裂解细胞。
然后加入终浓度为200 μg/ml的蛋白酶K，并在50°C下继续孵育2小时。
通过两轮25:24:1的酚-氯仿-异戊醇抽提后，再进行乙醇沉淀来纯化DNA。
将沉淀的DNA溶于含0.02% Triton X-100的10 mM Tris（pH 8.0）溶液中，并在4°C条件下溶解两天。

Para_02

文库是使用超长DNA测序试剂盒（牛津纳米孔科技公司，SQK-ULK001）构建的，并对生产商的实验方案进行了修改：按照方案所述，将约40μg DNA与FRA酶和FDB缓冲液混合，在室温下孵育5分钟，随后在75°C下热灭活5分钟。
将RAP酶与DNA溶液混合后，在室温下孵育1小时，然后进行纯化步骤。
纯化采用Nanobind UL文库制备试剂盒（Circulomics，NB-900-601-01），并用450μl EB洗脱。
之后，将75μl文库加载到预处理过的FLO-PRO002 R9.4.1测序芯片上，使用PromethION平台进行测序（使用的MinKNOW软件版本为v.21.02.17–23.04.5），并在测序开始后的24小时和48小时分别进行两次核酸酶清洗和重新加载。
所有G1至G3代ONT碱基识别均使用Guppy软件（版本6.3.7）完成。

Element (AVITI) sequencing

元素（AVITI）测序

Para_01

Element WGS数据是根据制造商的推荐生成的。
简而言之，DNA如上所述从全血中提取。
使用机械剪切制备了无PCR的文库，得到约350 bp的片段，并使用Element Elevate文库制备试剂盒（Element Biosciences，型号830-00008）。
通过定量PCR对线性文库进行定量，并在AVITI 2 × 150 bp流动池（Element Biosciences，尚未商业化）上进行测序。
使用Bases2Fastq软件（Element Biosciences）生成解复用的FASTQ文件。

Strand-seq library preparation

链特异性测序文库制备

Para_01

使用改良版的已建立的OP-Strand-seq方案制备单细胞Strand-seq文库，改良内容如下。
简而言之，将来自G1–3阶段的EBV细胞在BrdU存在的情况下培养24小时，并如先前所述使用荧光激活细胞分选技术对处于细胞周期G1期且含有BrdU的细胞核进行分选。
接下来，将单个细胞核分配到开放式的72×72孔纳米孔板的各个孔中，并用热不稳定性蛋白酶处理。
然后用限制性内切酶AluI和HpyCH4V（NEB）代替微球菌核酸酶（MNase）消化DNA。
随后，对片段进行加A尾处理，连接叉状接头，经过紫外照射处理后，使用带有索引序列的引物进行PCR扩增。
与MNase产生的DNA末端不同，使用限制性内切酶能够产生短而可重复、平端的DNA片段（超过90%小于1 kb），这些片段在连接接头前不需要末端修复。
省略末端修复步骤使得可以在分配单个细胞核之前就加入索引引物。
预先点样并干燥的引物能够在不影响后续文库制备步骤的前提下存活下来，直到PCR扩增阶段才发挥作用。
相比于在文库制备过程中在阵列之间转移索引引物的方法，提前点样索引引物更为可靠。
将Strand-seq文库合并后用AMPure XP磁珠纯化，并在NextSeq 550或AVITI（Element Biosciences）系统上进行PE75测序之前，通过凝胶纯化300至700 bp之间的文库片段。
补充图40展示了使用限制性内切酶制备的Strand-seq文库的示例。

Strand-seq data post-processing

链读长测序数据后处理

Para_01

使用BWA66（v.0.7.17-r1188）将解复用的FASTQ文件比对到GRCh38和T2T-CHM13参考基因组上（补充表14），以进行标准文库筛选。
通过SAMtools67（v.1.10）按照基因组位置对已比对的读段进行排序，并使用sambamba68（v.1.0）标记重复读段。
利用ASHLEYS69（v.0.2.0）预筛选通过质量过滤的文库。
我们还对手动评估了这些被选中的Strand-seq文库，并进一步排除了覆盖度不均匀，或者存在过多‘背景读段’（对于预期仅继承Crick链或Watson链的染色体，映射到相反方向的读段）的文库，如前所述。
这样做的目的是确保准确检测倒位并进行分型。

Strand-seq inversion detection

链读序反转检测

Para_01

如前所述，通过将Strand-seq读段的方向映射到参考基因组来检测G1–G3的多态性倒位71,72。
对于每个样本，我们选择了60多个Strand-seq文库（范围为62–90），每个文库中位数约为274,000条读段，且每条文库中 mapping quality ≥10的读段覆盖了大约0.67%的基因组（T2T-CHM13）（补充图41）。
然后我们在选定的Strand-seq文库中运行breakpointR73（v.1.15.1），以检测链状态变化的断点73。
我们使用这些结果，根据之前描述的方法71，利用breakpointR中的‘synchronizeReadDir’功能生成样本特异性的复合文件。
同样地，我们在这种复合文件上再次运行breakpointR，以检测Strand-seq读段以反向方向映射、提示存在倒位的区域。
最后，我们通过在UCSC基因组浏览器74中检查Strand-seq读段的映射情况，手动评估每一个报告的倒置区域，并去除任何低置信度的调用。
我们还使用Strand-seq数据对所有倒位进行了分型，并根据单倍型一致性将分型结果与分型的基因组组装结果同步。
最后，我们评估了所检测到并完全分型的倒位的孟德尔一致性。
我们将那些至少有一半的G3样本被Strand-seq完全分型，并且与可能遗传自G2父母的等位基因一致的位点标记为符合孟德尔遗传规律的位点（补充表7）。

Generation of phased genome assemblies

分阶段基因组组装的生成

Para_01

使用两种不同的算法生成了分阶段的基因组组装，分别是Verkko（v.1.3.1和v.1.4.1）和带有ONT支持的hifiasm（UL）（v.0.19.5）。
由于Verkko和hifiasm算法在积极开发中，因此使用了两个不同版本生成组装结果。
G2–G3的分阶段组装是通过结合HiFi和ONT测序数据，并利用父母本的Illumina k-mer进行分型而生成的。
为了生成G1的分阶段基因组组装，我们仍然使用HiFi和ONT测序数据结合Verkko流程，并使用Strand-seq对组装图进行分型。
最后，G4样本仅使用HiFi测序数据并通过hifiasm（v.0.19.5）进行组装。

Para_02

请注意，使用 Verkko 进行基于三人组的分型时，单倍型1被指定为母本，单倍型2被指定为父本。
相比之下，对于 hifiasm 组装，我们报告的单倍型标签进行了调换，即单倍型1为父本，单倍型2为母本，以符合 HPRC 对 hifiasm 组装的标准。

Evaluation of phased genome assemblies

分阶段基因组组装的评估

Para_01

为了评估每个分型组装结果的碱基对和结构准确性，我们使用了多种组装评估工具，以及正交数据集，如 PacBio HiFi、ONT、Strand-seq、Illumina 和 Element 数据。
已知的组装问题列于补充表4中。
我们注意到，在基于组装的变异检测结果中，我们修正了四处单倍型切换错误，以避免后续分析中出现偏差。
本文中使用的组装质量术语详见补充说明8。

Strand-seq validation

Strand-seq 验证

Para_01

我们使用Strand-seq数据来评估每个分型组装的方向性和结构性准确性。
首先，我们使用BWA66（版本0.7.17-r1188）将每个样本选定的Strand-seq文库比对到分型从头组装结果上。
接下来，我们以比对生成的BAM文件作为输入运行breakpointR73（版本1.15.1）。
然后，我们使用breakpointR中的createCompositeFiles函数创建方向性组合文件。
之后通过runBreakpointR函数在这些组合文件上运行breakpointR。
这为我们提供了在任意给定样本中所有单细胞Strand-seq文库间发生链状态变化的区域。
许多这样的区域表明是真实的杂合倒位。
然而，相对于周围区域，Strand-seq读段以相反方向映射的区域可能是由取向错误导致的。
此外，在多个文库中链状态反复改变的位置可能是组装连接错误的标志。
因此，我们对这些区域进行了更详细的检查，以排除任何大型结构组装不一致的问题。

Read to assembly alignment

读取到组装的比对

Para_01

为了评估从头组装的准确性，我们使用Winnowmap76（版本2.03）将样本特异性的PacBio HiFi读段比对到相应的分型基因组组装结果中，所使用的参数如下：

Para_02

-I 10G -Y -ax map-pb --MD --cs -L --eqx

Flagger validation

标记验证

Para_01

使用 Flagger9 通过将 HiFi 读段比对到组装结果以及将组装结果比对到参考基因组来检测错误组装。
根据读段比对的分歧程度和特定的参考偏倚区域来标记可疑区域。
使用了一个参考特异性的 BED 文件（chm13v2.0.sd.bed），设置最大读段分歧度为 2%，并指定了参考偏倚区块。
分析这些被标记的区域以识别坍塌区域、假重复区域、错误区域以及具有预期读段覆盖度的正确组装的单倍区块。

Para_02

我们使用了 Flagger v.0.3.3（https://github.com/mobinasri/flagger）来运行 flagger_end_to_end WDL。

Para_03

所需的输入包括以下内容：

[ol]- (1) Read-to-contig alignments—Winnowmap alignments of all HiFi reads to the assembly (hap1, hap2 and unassigned.fasta) - (2) A combined assembly fasta file with hap1, hap2 and unassigned contigs - (3) BAM alignments of assembly to the CHM13v2.0 reference

Para_04

使用 GitHub 上提供的一个流程（https://github.com/mrvollger/asm-to-reference-alignment），将组装得到的 hap1、hap2 和未分配的 fasta 文件与 CHM13v2.0 进行了比对。

NucFreq validation

NucFreq 验证

Para_01

使用 NucFreq77（v.0.1）计算了通过 Winnowmap76（v.2.03）比对的 HiFi 读段的核苷酸频率。这用于识别发生坍缩的区域，即次高核苷酸计数超过 5 的位置；以及错装的区域，即所有核苷酸计数均为零的位置。

Para_02

NucFreq分析流程可在GitHub上获取（https://github.com/mrvollger/NucFreq）。

Assembly base-pair quality

组装碱基对质量

Para_01

为了评估基因组组装的准确性，我们使用了一个流程，该流程使用 Meryl78（v.1.0）通过以下命令从 Illumina 读段中计数长度为 21 的 k-mers：

Para_02

meryl k=21 count {input.fastq"} output {output.meryl"}

Para_03

我们随后使用了 Merqury78（v.1.1），该工具将测序读取中的 k-mer 与组装基因组中的 k-mer 进行比较，并标记那些仅在组装结果中唯一存在的 k-mer 的位置。
这些唯一的 k-mer 表明可能存在碱基对错误。
Merqury 接着基于来自 Meryl 的 k-mer 计数所估计的 k-mer 存活率来计算质量值，从而提供一种量化指标以评估基因组组装的完整性和正确性。

Gene completeness validation

基因完整性验证

Para_01

为了评估我们组装结果中单拷贝基因的完整性，我们使用了Compleasm79（v.0.2.4）。
更多详细信息请访问GitHub（https://github.com/huangnengCSU/compleasm）。

Para_02

我们使用以下参数运行了 compleasm：

Para_03

compleasm.py run -a {assembly.fasta"} -o results/{sample.id"} -t {threads"} -l {params.lineage"} -m {params.mode"} -L {params.mb_downloads"}

Para_04

-l 灵长类 -m busco -L {params.mb_downloads"}

Para_05

并使用以下命令下载：compleasm_kit/compleasm.py download primates。

Assembly to reference alignment

将装配体对齐到参考

Para_01

所有从头组装的数据均使用minimap2（版本2.24）比对到GRCh38以及完整的人类参考基因组T2T-CHM13（版本2）上，使用的命令如下：

Para_02

minimap2 -K 8G -t {线程数"} -ax asm20 \

Para_03

--secondary=no --eqx -s 25000 \
让我们一步一步地思考。
记住，直接给我整理的结果，不要有除了结果的其他任何语句。一般情况下，待整理原本有多少句话，输出的 json 就有多少个 item。

Para_04

input.ref"} {input.query"} \samtools view -F 4 -b - > {output.bam"

Para_05

该参考序列比对的完整流程可在GitHub上获取：(https://github.com/mrvollger/asm-to-reference-alignment)。

Para_06

我们还使用 rustybam（v.0.1.33）（https://github.com/mrvollger/rustybam）中的 trim-paf 功能生成了这些比对结果的修剪版本，以去除在高度相同的 SD 区域中出现的冗余比对。
通过这种方式，我们旨在减少单个连续序列在这些重复区域上多次比对所产生的影响。

Definition of stable diploid regions

稳定二倍体区域的定义

Para_01

为了进行这项分析，我们使用了比对到参考基因组的组装结果（见‘组装与参考基因组比对’部分），其以PAF文件格式呈现。
我们采用了rustybam工具中的trim-paf功能所输出的经过修剪的PAF文件。
稳定的二倍体区域被定义为：在单倍型1和单倍型2中，分型后的基因组组装结果分别只对应一个连续序列（contig）的比对，并且这些区域被标记为‘2n’区域。
若某一区域中每个单倍型存在两个或两个以上的比对结果，则被标记为‘multi’比对。
最后，那些在一个单倍型中仅出现单一连续序列比对的区域被标记为‘1n’区域。
这些报告是通过使用‘getPloidy’ R函数生成的（代码可用性见后文）。

Detection and analysis of meiotic recombination breakpoints

减数分裂重组断点的检测与分析

Para_01

我们使用三种正交方法构建了G2和G3个体的高分辨率重组图谱，这些方法要么基于不同的测序技术，要么应用了不同的数据分析检测算法。
第一种方法是基于Strand-seq数据，利用R软件包StrandPhaseR81（版本0.99）提取染色体长度的单倍型。
第二种方法是通过家系谱系中遗传变异的孟德尔一致性推导出的遗传向量来实现。
最后一种方法是采用基于三人组（trio）的相位解析基因组组装，并随后使用PAV和Dipcall进行小变异检测，以更精确地确定减数分裂的断裂点。

Para_02

相比之下，G4个体的重组图谱是通过结合G3配偶的Strand-seq数据以及基于组装的G4样本变异集合（Dipcall）构建的。
由于使用了基于Strand-seq且变体密度较低的G4配偶单倍型，所报告的重组断点分辨率低于G2和G3样本（补充表8）。

Detection of recombination breakpoints using circular binary segmentation

使用圆形二进制分割检测重组断点

Para_01

为了使用环状二进制分割法绘制减数分裂重组断点，我们使用了两个不同的数据集。
第一个数据集代表由基于Strand-seq的相位分析（SSQ）报告的已分型的小变异（SNVs和indels）。
另一个数据集基于三联体分型组装中报告的小变异，这些变异常通过PAV（v.2.3.4）或Dipcall（v.0.3）进行检测。
通过这种方法，我们旨在检测重组断点，即某个子代单倍型从匹配一个亲本的H1切换到H2，或者相反的位置。
为了检测这些位置，我们首先通过计算子代等位基因与每个亲本中的纯合变体的一致性水平，确定子代中哪一个同源染色体来自特定的亲本。
接下来，我们将每个子代的同源染色体与对应亲本的两个同源染色体进行比较，并根据是否匹配H1或H2分别编码为0或1。
然后，我们在这样的二值向量上应用了一个环状二进制分割算法，具体实现是R语言fastseg包（版本1.46.0）中的fastseg函数，参数设置如下：fastseg(binary.vector, minSeg=20, segMedianT=c(0.8, 0.2))。
在Strand-seq单倍型稀疏的情况下，我们将fastseg参数minSeg设为20；而在基于组装的密集单倍型情况下，为了达到可比较的检测重组断点的灵敏度，我们对基于Dipcall和PAV的变异调用分别采用了更大的窗口大小400和500。
随后将分割均值≤0.25的区域标记为H1，而分割均值≥0.75的区域则被指定为H2。
介于这两个值之间的分割均值区域被认为是模糊的，因此被排除。
此外，我们过滤掉了短于500kb的区域，并合并了连续且分配相同单倍型的区域（代码可用性）。

Detection of meiotic recombination breakpoints using inheritance vectors

使用遗传向量检测减数分裂重组断点

Para_01

通过来自G1、G2和G3家系成员的HiFi测序数据进行DeepVariant变异检测（见‘基于读长的变异检测’部分），使我们能够识别G3中杂合位点的来源单倍型，并在染色体上不同位点间来源单倍型发生变化时推断出重组事件的发生。
首先应用深度过滤器去除每个样本中预期覆盖度分布之外的变异，从而初步确定了遗传向量的大致轮廓；
然后使用自定义脚本勾勒出遗传模式，要求至少需要有10个SNV支持特定的单倍型，并通过手动校正以删除生物学上不太可能的单倍型区块，或在存在支持证据的情况下添加额外的单倍型区块，并进一步优化单倍型的坐标信息。
缺失的重组事件通过出现违反家系遗传规律的变异区块来识别，并与基于组装的重组检测结果中的位置相对应。
我们开发了一个隐马尔可夫模型框架，利用Viterbi算法从SNV位点识别最可能的遗传向量序列。
有关转移概率和发射概率的详细信息，请参见参考文献18以及相关的GitHub仓库（https://github.com/Platinum-Pedigree-Consortium/Platinum-Pedigree-Inheritance）。

Para_02

转移矩阵定义了给定遗传状态转移（重组）的概率。
发射矩阵定义了在特定位点上的变异调用准确描述遗传状态的概率。
通过优化转移矩阵和发射矩阵中的数值，以重现先前已识别的遗传向量，并正确识别缺失的向量。
Viterbi算法识别出539次重组，母系重组率为每兆碱基1.29厘摩，父系重组率为每兆碱基0.99厘摩。
在系谱中观察到母系偏差，其中57%的重组来源于G3的母系

Merging of meiotic recombination maps

合并减数分裂重组图谱

Para_01

通过不同的正交技术和算法报告的减数分裂重组断点（参见‘使用圆形二元分割检测减数分裂重组断点’和‘使用遗传向量检测减数分裂重组断点’章节）分别在G2和G3样本中进行了合并。
我们从G3重组图谱开始，其中使用基于遗传的图谱作为参考，并随后在基于PAV、Dipcall和Strand-seq（SSQ）分型变异体报告的重组图谱中查找每个参考断点的支持证据。
如果对于给定样本和同源染色体，另一种技术报告的断点与参考断点的距离不超过1 Mb，则认为该参考断点得到了支持。
任何距离更远的重组断点都被报告为唯一断点。
我们也对G2重组图谱重复了这一过程。
然而，在G2重组图谱的情况下，我们使用基于PAV的图谱作为参考。
这是因为基于遗传的方法需要三个世代来绘制G3中的重组断点。
此外，我们还报告了一列名为‘best.range’的内容，它表示跨越所有正交重组图谱中与给定参考断点直接重叠的最窄断点范围。
最后，我们报告了一列‘min.range’，它表示对于任何给定断点，在所有正交数据集中覆盖度最高的范围。
合并后的重组断点详见补充表8。

Meiotic recombination breakpoint enrichment

减数分裂重组断点富集

Para_01

我们检测了在G2–G4中发现的所有（n = 1,503）个重组断点相对于T2T-CHM13的富集情况，以确定它们是否根据亲本同源染色体的来源而聚集在染色体末端。
为此，我们分别统计了每条染色体末端最后5%区域内母系和父系重组断点的数量。
然后，我们在每条染色体上将检测到的重组断点随机打乱1,000次，并重新进行计数。
在置换分析中，我们使用了R软件包regioneR84（v.1.32.0）及其permTest函数，并采用了以下参数：

Para_02

permTest( A=breakpoints, B=chrEnds.regions, randomize.function=circularRandomizeRegions, evaluate.function=numOverlaps, genome=genome, ntimes=1000, allow.overlaps=FALSE, per.chromosome=TRUE, mask=region.mask, count.once=FALSE)

Refinement of meiotic recombination breakpoints using MSA

使用多重序列比对（MSA）对减数分裂重组断点进行精细分析

Para_01

到目前为止，所有的减数分裂重组断点都是通过与单一的线性参考序列（GRCh38 或 T2T-CHM13）进行比对而确定的。
为了消除与单一参考基因组比对可能引入的任何潜在偏差，我们尝试直接将每个子代所继承的同源染色体与提供该同源染色体的亲本单倍型进行比较，并据此优化检测到的重组断点。
我们首先从G3中选择‘best.range’列（补充表8），得到一组合并后的T2T-CHM13参考断点。
然后，对于每个断点，我们在断点边界两侧各设置一个750 kb的‘查找’区域，并使用SVbyEye85（v.0.99.0）中的subsetPafAlignments函数，将分相组装与参考序列（T2T-CHM13）的PAF比对结果截取为指定区域的子集。
接下来，我们从分相组装中提取特定区域的FASTA序列。
我们分别对子代继承的同源染色体（已重组）以及提供该同源染色体的对应亲本单倍型进行了上述操作。

Para_02

接下来，我们使用 R 软件包 DECIPHER（版本 2.28.0；函数为 AlignSeqs）为三个序列（子代继承的同源序列、亲本同源序列 1 和亲本同源序列 2）创建了多序列比对（MSA）。
由于使用 DECIPHER 对差异较大的序列进行最优比对所需的计算时间显著增加，因此跳过了长度差异超过 100 kb 或核苷酸数目差异超过 10,000 个碱基的序列的比对过程。
在构建 MSA 后，我们选取至少存在一个错配的位置，并排除那些两个亲本单倍型携带相同等位基因的位点。
重组断点是指子代继承的同源序列部分匹配来自亲本同源序列 1 的等位基因，部分匹配来自亲本同源序列 2 的等位基因的区域。
因此，我们跳过了子代等位基因与某一个亲本同源序列相似度超过 99% 的 MSA 的进一步分析。
如果通过该过滤条件，我们则使用自定义的 R 函数 getAlleleChangepoints（代码可用性说明）来检测子代所继承的单倍型从匹配亲本单倍型 1 的等位基因转换为匹配亲本单倍型 2 的等位基因的转换点。
此类特定于 MSA 的转换点随后被报告为一个新的区间，其中可能发生了重组断点。
最后，我们尝试通过从断点边界提取 1 kb 长的 k-mer 并使用 R 软件包 Biostrings（版本 2.70.2）中的函数 'matchPattern' 将这些 k-mer 与参考序列（按染色体）进行匹配（允许最多有 10 个错配），从而报告此类 MSA 特定断点的参考坐标。
经过细化的重组断点列表详见补充表 8。

Detection of allelic gene conversion using phased genome assemblies

使用分阶段基因组组装检测等位基因转换

Para_01

我们试图使用此前定义的该家族重组图谱来检测父母等位基因遗传中的较小局部变化。
我们对所有G3样本（n = 8）与G2代父母进行了这一分析。
为此，我们遍历每个孩子同源染色体（在每个样本中），并将其与该孩子同源染色体所继承自的双亲同源染色体进行比较。
我们通过比较孩子与其相应父母之间来自分相基因组组装的SNV和indel变异结果来进行此分析。
为了仅考虑可靠的变异，我们只保留那些至少被两个基于测序读长的变异检测工具支持的变异（分别为DeepVariant-HiFi、Clair3-ONT或dragen-Illumina变异集合）。
此外，我们仅保留那些在父母中为杂合且在孩子中也被检出的变异位点。
经过如此严格的变异过滤后，我们一次滑动两个连续的孩子变异位点，并将它们分别与相应来源父母的单体型1和单体型2进行比对。
对于这种相似性计算，我们使用了自定义的R函数getHaplotypeSimilarity（代码可用性部分）。
然后，针对由重组断点定义的每个单型段，我们报告那些至少有两个连续变异与重组图谱定义的预期父母同源染色体相反的父系单体型匹配的区域。
此外，我们将间隔≤5 kb的连续区域进行合并。
对于疑似基因转换事件列表，我们仅保留未被Flagger标记为问题区域的区域。
我们还去除了距离先前定义的重组事件≤100 kb的区域，以及与着丝粒卫星区域和高度相同的片段重复序列（≥99%相同）重叠的事件。
最后，我们通过可视化检查分相的HiFi读长对疑似等位基因转换事件列表进行了评估。

Read-based variant calling

基于读取的变异检测

Para_01

PacBio HiFi 数据使用了可在 GitHub 上获取的人类全基因组测序工作流（human-WGS-WDL）进行处理（版本v1.0.3，地址：https://github.com/PacificBiosciences/HiFi-human-WGS-WDL/releases/tag/v1.0.3）。
该流程对数据进行比对、分型并检测小变异（使用 DeepVariant v1.6.0）以及结构变异（使用 PBSV v2.9.0；地址：https://github.com/PacificBiosciences/pbsv）。
我们使用经过比对并带有单倍型标签的 HiFi BAM 文件进行所有下游的 PacBio 分析。

Clair3

Clair3

Para_01

Clair3（版本1.0.7）的变异检测是分别基于PacBio HiFi和ONT（ont_guppy5）数据的默认模型比对结果进行的，且启用了染色体分型和gVCF生成功能。
变异检测在每条染色体上单独进行，并合并生成一个VCF文件。
随后，将gVCF文件通过自定义配置文件输入到GLNexus中进行进一步处理。

Para_02

PacBio HiFi

Para_03

run_clair3.sh --bam_fn={input.bam"} --sample_name={sample"} --ref_fn={input.ref"} --threads=8 --platform=hifi --model_path=/path/to/models/hifi --output={output.dir"} --ctg_name={contig"} --enable_phasing --gvcf

Para_04

Para_05

run_clair3.sh --bam_fn={input.bam"} --sample_name={sample"} --ref_fn={input.ref"} --threads=8 --platform=ont --model_path=/path/to/models/ont_guppy5 --output={output.dir"} --ctg_name={contig"} --enable_phasing –gvcf

Para_06

ONT测序数据用于Clair3变异检测时，采用minimap2（v.2.21）进行比对，使用的参数为：-L --MD --secondary=no --eqx -x map-ont。

Generation of truth set of genetic variation using inheritance vectors

使用遗传向量生成真实遗传变异集合

Para_01

我们采用了一个先前已建立的框架来定义真实的基因变异。
与基于家系的过滤方法相比，我们的分析使用全部四个等位基因来检测基因型错误，而在一个三人组（trio）中只能传输并观察到两个等位基因。
通过将第三代的基因型模式与第一代的单倍型（A,B,C,D）进行比对，我们可以检验从第二代到第三代的等位基因传递是否正确。
我们构建了一张覆盖第三代的单倍型图谱（遗传向量），据此可以对变异调用进行判定。
为了检验家系一致性，我们编写了代码，利用遗传向量作为预期的单倍型，并测试查询VCF文件中的可能基因型组合。
利用单倍型结构，我们将符合家系遗传规律的变异进行单倍型推断（phasing）。
这些功能被实现为一个独立的二进制工具，该工具需要输入遗传向量和标准格式的VCF文件，例如：

Para_02

concordance -i ceph.grch38.hifi.g3.csv –father NA12877 –mother NA12878 –vcf input.vcf –prefix pedigree_filtered > info.stdout

Para_03

构建小变异真实集的谱系过滤及其他步骤可在GitHub上获取（https://github.com/Platinum-Pedigree-Consortium/Platinum-Pedigree-Inheritance）。

Detection of small de novo variants

小型新生变体的检测

Para_01

根据之前列出的参数，我们使用GATK HaplotypeCaller（版本4.3.0.0）和DeepVariant（版本1.4.0）对与T2T-CHM13比对的HiFi数据进行变异检测，并简单地识别了每个G2和G3样本特有的变异。
我们将SNV和indel的检测结果分开，并应用了基本的质量过滤标准，例如删除在1kb窗口内出现三个或更多SNV的簇。
我们将这种由次级检测方法生成的变异检测集合合并（https://github.com/Platinum-Pedigree-Consortium/Platinum-Pedigree-Inheritance/blob/main/analyses/Denovo.md），并对所有检测结果进行了以下验证流程。

Para_02

我们通过检查HiFi、ONT和Illumina读数中的SNVs和indels来对它们进行验证，并排除那些未达到比对质量阈值（长读长为59，短读长为0）的读数。
在变异位点处，分别统计具有高碱基质量（>20）和低碱基质量（<20）的读数。
我们保留了在儿童样本中至少出现在两种测序数据中的变异，并且这些变异应不存在于高质量的父母读数中。
对于SNV检测，我们接下来检查了系谱中每个样本的HiFi数据。
如果某个SNV在不属于新发样本直接后代的所有家庭成员中均未出现，则判定该SNV为真实的新生突变。
最后，我们检查了每个变异的等位基因平衡，确定哪些变异位于TRs中，并跨所有测序平台重新评估了父母读数数据，去除了那些测序数据噪声较大或支持替代等位基因的低质量父母读数超过两个的变异（补充说明9）。

DNM phasing and postzygotic assignment

DNM 分型与合子后分配

Para_01

为了确定新发SNVs的亲本来源，我们重新检查了包含新发等位基因的长读段数据。
首先，我们使用最初的GATK变异检测结果，在DNM周围80 kb的窗口内识别具有信息量的位点。
我们选择那些可以将一个等位基因唯一地分配给一个亲本的单核苷酸多态性（SNP）位点（例如，某个位点在父亲中是纯合参考型，而在母亲中是杂合型）。
对于每一个DNM，我们评估了比对到该新发等位基因位点上的每一条ONT和HiFi读段，并根据是否存在有信息量的SNP将其分配到父源或母源单倍型。
分配过程按照之前描述的方法计算遗传得分来进行。
那些仅被分配到母源或父源单倍型的DNM被认为是成功定相的，而那些分配到矛盾单倍型的DNM则被排除在最终的变异集合之外。
如果未定相的变异在不同测序平台间的等位基因比例没有显著差异（通过χ2检验），并且其综合等位基因比例与0.5存在显著差异，则被判定为合子后起源（n = 7）。

Para_02

一旦我们将每一条读段分配到亲本单倍型后，就统计了母本和父本读段中具有参考等位基因或替代等位基因的数量。
如果某个新生突变（DNM）出现在某一亲本单倍型的所有读段中，则我们判定该突变为种系起源（germline）。
相反，如果某个新生突变仅出现在某一亲本单倍型的部分读段中，则我们判定该突变为合子后起源（postzygotic）。

Sex chromosome DNM calling and validation

性染色体DNM检测与验证

Para_01

为了识别X染色体上的新生突变（DNMs），我们采用了与常染色体变异相同的策略，只有一个例外：我们仅使用了由GATK生成的变异检测结果。
对于男性个体，我们以单倍体型模式重新运行GATK，使其在X染色体上仅识别一种基因型。

Para_02

为了识别Y染色体上的新生突变（DNMs），我们将男性样本的HiFi、ONT和Illumina数据比对到G1-NA12889染色体Y组装序列上，然后在比对后的HiFi数据上使用GATK以单倍体模式进行变异检测。
我们将每个男性个体直接与其父亲进行比较，筛选出仅存在于子代的变异。
我们通过检查父亲的HiFi、ONT和Illumina数据来验证SNVs和indels，并排除在亲本读段中出现的任何变异，采用与常染色体变异相同的判断逻辑。

Callable genome and mutation rate calculations

可调用的基因组和突变率计算

Para_01

为了确定我们在基因组中能够识别新发变异的位置，我们对每个三人组的HiFi数据进行了评估。
我们首先使用GATK HaplotypeCaller90（版本4.3.0.0）并选择‘ERC BP_RESOLUTION’选项，在基因组的每个位点生成基因型调用。
仅考虑父母双方都被基因型为纯合参考型（0/0）的位点，因为具有父母替代等位基因的位点被排除在我们新发变异发现流程之外。
随后，我们检查了样本及其父母的HiFi读段，仅限于比对质量至少为59的主要比对结果。
对于子代，我们只考虑来源于血液的HiFi读段，而对于父母则同时考虑血液和细胞系数据。
我们在基因组的每个位点统计了最小碱基质量得分达到20的读段数量，并将这些信息与我们的变异调用结果结合。
如果一个位点的双亲和子代各自至少拥有一条高质量读段且该位点碱基质量高，则该位点被认为是可以被调用的。
在常染色体上，我们观察到平均有2.67 Gb的可及序列（总共2.90 Gb，标准差=24.9 Mb）。
对于女性子代，X染色体的可调用性以同样的方式确定；而对于男性子代，在分析X染色体时我们仅考虑母亲的HiFi数据，在分析Y染色体时则仅考虑父亲的HiFi数据。
此外，男性的性染色体并不局限于双亲都被基因型为参考型的位点——允许每个父母携带替代等位基因。

Para_02

我们通过将每个样本的种系常染色体新生突变（DNMs）数量除以被判定为可分析的碱基对数目的两倍，计算出了该样本的种系常染色体突变率。
对于PZMs，我们使用了相同的分母进行计算。
在女性个体中，可分析性染色体的数量被定义为X染色体上可分析碱基数目的两倍；
而在男性个体中，则被定义为X染色体和Y染色体上可分析碱基总数之和。
针对每一个特定特征的突变率（例如SDs），我们将样本的新生单核苷酸变异（SNVs）以及该样本的可分析区域与相关特征的坐标进行了交集处理。
接着，我们通过将区域内SNVs的数量除以可以可靠进行比对的、可分析的基因组序列总量来计算突变率。

Analysis of STRs and VNTRs

STRs 和 VNTRs 的分析

Para_01

鉴于检测STRs（1–6 bp重复单元）和VNTRs（≥7 bp重复单元）中的突变存在挑战，我们采用了一种靶向的HiFi基因分型策略，并通过遗传传递分析和正交测序进行验证。

Defining the TR catalogues

定义TR目录

Para_01

使用了命令 trf-mod -s 20 -l 160 {reference.fasta"}，从而得到最小参考位点大小为 10 bp，基序大小范围为 1 至 2,000 bp（https://github.com/lh3/TRF-mod）。
将距离在 50 bp 内的位点合并，然后丢弃任何大于 10,000 bp 的位点。
利用 tr-solve（https://github.com/trgt-paper/tr-solve）对剩余位点进行注释，以解析复合位点的结构。
仅考虑了染色体 1–22、X 和 Y 上的串联重复序列注释（数据可用性）。

TR genotyping with TRGT

使用TRGT进行TR基因分型

Para_01

TRGT32 是一个利用 PacBio HiFi 测序读长进行 TR 等位基因分型的软件工具（https://github.com/PacificBiosciences/trgt）。
在提供比对好的 HiFi 测序读长（BAM 格式）以及一个列出了多个 TR 位点的基因组位置和基序结构的文件的前提下，TRGT 将会返回一个包含每个 TR 位点推断基因型的 VCF 文件。
在本分析中，我们使用上述定义的 TR 目录，对 CEPH 1463 家系中的每一个成员运行了 TRGT（版本号：v.0.7.0-493ef25）。
TRGT 的运行采用了默认参数。

Para_02

trgt --threads 32 --genome {in_reference"} --repeats {in_bed"} --reads {in_bam"} --output-prefix {out_prefix"} --karyotype {karyotype"}` bcftools sort -m 3072M -Ob -o {out_prefix"}.sorted.vcf.gz {out_prefix"}.vcf.gz
bcftools index --threads 4 {out_prefix"}.sorted.vcf.gz
samtools sort -@ 8 -o {out_prefix"}.spanning.sorted.bam {out_prefix"}.spanning.bam
samtools index -@ 8 {out_prefix"}.spanning.sorted.bam

Measuring concordant inheritance of TRs

测量TR的协调遗传性

Para_01

为了确定TRs的符合孟德尔遗传规律的情况，我们计算了所有可能组合中曼哈顿距离的最小值，这些组合包括先证者等位基因长度（AL）与父母双方AL值的所有可能情况。
如果从所有计算出的距离中发现最小的曼哈顿距离为0，则认为该位点是符合孟德尔遗传规律的，这意味着先证者的AL值的某种组合能够完美匹配父母的AL值。
相反，如果最小的曼哈顿距离大于0，这表明先证者AL值的所有组合都与父母AL值存在某种程度的偏差，我们将该位点视为不符合孟德尔遗传规律，并将其记录为潜在的孟德尔遗传错误。
对于每个TR位点，我们分别统计了符合遗传规律的三人组数量、出现孟德尔遗传错误的三人组数量以及因缺失值而无法完全进行基因分型的三人组数量。
在计算符合率时，排除了任何存在缺失基因型数据的位点；然而，在以下分析中新发现变异时，仍然考虑了完整的个体三人组。

Calling de novo TRs

调用从头TRs

Para_01

我们选择对G3进行从头TR检测，这出于几个原因。
首先，它们的G2代父母（NA12877和NA12878）分别被测序至99倍和109倍的HiFi深度，这使得相较于测序深度较低的样本，发生父母等位基因扩增丢失的可能性要低得多。
其次，G1代DNA来源于细胞系，这增加了在G2代中检测DNMs时产生假象的风险。
最后，在本研究中拥有已测序子代的G3代两个个体中的DNMs可以通过遗传得到进一步验证。

Para_02

我们使用了TRGT-denovo33（版本0.1.3），这是TRGT的配套工具，旨在利用HiFi测序数据对家庭三重组中的TR新生突变（DNMs）进行深入分析（https://github.com/PacificBiosciences/trgt-denovo）。
TRGT-denovo 使用 TRGT 生成的共识等位基因序列和基因分型数据，并结合预测这些等位基因序列时所使用的跨域 HiFi 读段提供的额外证据。
简而言之，TRGT-denovo 从每个家庭成员（母亲、父亲和孩子）中提取并划分跨域读段到它们最可能的等位基因。
将父母的跨域读段重新比对到孩子两个共识等位基因序列中的每一个，并为每个读段计算比对得分（该得分总结了父母读段与共识等位基因序列之间的差异）。
在每个TR位点上，分别将孩子的两个等位基因视为潜在的新生候选等位基因。
对于孩子的每个等位基因，TRGT-denovo 报告是否存在支持新生事件的证据，这些证据包括：denovo_coverage（支持孩子中独特 AL 的读段数量，且该 AL 在父母读段中缺失）；
overlap_coverage（父母中支持某个 AL 的读段数量，该 AL 与潜在的新生等位基因高度相似）；以及新生事件的大小（表示为具有 de novo coverage 的读段比对得分相对于最近的父母等位基因的绝对平均差异）。

Calculating the size of a de novo TR expansion or contraction

计算一个新生TR扩增或收缩的大小

Para_01

我们根据最可能发生过收缩或扩展事件的亲本TR等位基因，测量了新生TR等位基因的大小。
如果TRGT-denovo在特定位点报告了一个新生的扩展或收缩事件，我们会采用以下方法来计算该事件的大小。

Para_02

根据 TRGT 为三人组（trio）中的每个成员报告的 ALs，我们计算了子代中新生的 TR 等位基因与子代双亲中所有四个 TR 等位基因之间的大小差异（我们称之为‘diff’）。
例如，如果 TRGT 报告的父亲的 ALs 是 100、100，母亲是 50、150，孩子的 ALs 是 200、100，并且报告称孩子中的长度为 200 的等位基因为新生等位基因，则父亲的 diffs 会是 100、100，而母亲的是 150、50。
如果我们能够将新生 TR 等位基因定位于某一特定父母来源，则只需在该父母的 ALs 中找出最小的 diff，并将其视为可能的扩增或收缩大小。
否则，我们假设所有父母 ALs 中最小的 diff 表示可能的新生等位基因的大小。

De novo filtering

从头过滤

Para_01

我们对候选的TR-DNMs（由TRGT-denovo识别得到）应用了一系列过滤器，以去除可能的假阳性结果。
对于在每个子代中发现的每一个新生等位基因，我们要求满足以下条件（补充说明9和10）：

[ul]- HiFi sequencing depth in the child, mother, and father ≥10 reads. - The candidate de novo AL in the child must be unique: as in ref. 37, we removed candidate de novo TR alleles if (1) the child’s de novo AL matched one of the father’s ALs and the child’s non-de novo AL matched one of the mother’s ALs or (2) the child’s de novo AL matched one of the mother’s ALs and the child’s non-de novo AL matched one of the father’s ALs. - The candidate de novo allele must represent an expansion or contraction with respect to the parental allele. - At least two HiFi reads supporting the candidate de novo allele (denovo_coverage ≥ 2) in the child, and at least 20% of total reads supporting the candidate de novo allele (child_ratio ≥ 0.2). - Fewer than 5% of parental reads likely supporting the candidate de novo AL in the child.

Para_02

为了计算某个特定个体中的串联重复序列（TR）新生突变（DNM）率，我们首先计算了使用TRGT基因分型的约780万个TR位点中，在目标个体及其父母三人组的每个成员中至少被10条HiFi测序读长覆盖的总TR位点数量。
接着，我们将新发现的TR等位基因总数除以可检测的位点总数，得到一个每代每个位点上的总体DNM率。
最后，我们将该突变率再除以2，以获得每个位点、每个单倍型、每代的突变率。
如图3a所示，我们还根据位点内的最小基序大小估计了DNM率的变化情况。
例如，具有AT(n)AGA(n)T(n)基序结构的位点其最小基序大小为1。
我们统计了最小基序大小为N的位点上发生的TR DNM数量，并将该数量除以通过过滤阈值的最小基序大小为N的TR位点总数。
然后我们将这个比率再除以2，以获得每个位点、每个单倍型、每代的突变率。
在计算短串联重复序列（STR）、可变数目串联重复序列（VNTR）和复杂突变率时，我们将STR位点定义为所有组成基序长度在1至6个碱基对之间的位点；
将VNTR位点定义为所有基序长度大于6个碱基对的位点；
将复杂位点定义为同时包含STR（1–6 bp）和VNTR（≥7 bp）基序的位点。
例如，A(n)位点和AT(n)AGA(n)T(n)位点都会被归类为STR，因为它们都仅包含STR基序。

Para_03

以往的研究通常在目标群体中具有多态性的位点上测量短串联重复序列（STR）的突变率。
为了生成与这些先前研究更为一致的突变率估计值，我们还计算了CEPH 1463家系中表现出多态性的STR位点数量。
如果在一个特定的TR位点上，CEPH 1463个体间至少观察到两个不同的等位长度（AL），则该位点被定义为具有多态性。
我们注意到，这种对多态性STR的定义既受研究群体规模的影响，也受用于基因分型的测序技术的影响。
如之前研究所讨论的，多态性位点的数量与群体的大小成正比。
此外，如果我们观察到群体中某个位点存在多个不同的AL，则将该位点定义为多态性位点；然而，在这种情况下，如果HiFi测序读长在这些位点上产生了大量影子峰（stutter），导致不同个体间的STR AL估计值出现变化，则可能会错误地将某些位点归类为多态性位点。
总计有1,096,430个STR位点在该群体中表现出多态性。
为了计算每个G3个体中的突变率，我们采用了上述描述的相同覆盖度质量阈值。

Phasing of TRs

TR 的阶段划分

Para_01

由TRGT基因分型的STR通过HiPhase92（版本1.0.0-f1bc7a8）进行了单倍型推断。
我们按照HiPhase的指南，将来自DeepVariant、PBSV和TRGT的相关VCF文件输入HiPhase，对小变异、结构变异（SV）和串联重复序列（TR）进行联合单倍型推断，从而为每个分析样本生成了三个经过单倍型推断的VCF文件。
我们还通过--global-realignment-cputime参数启用了全局重新比对，以提高等位基因分配的准确性。
请注意，HiPhase在单倍型推断过程中会特别排除完全位于基因分型STR区域内的变异。
这样设计的原因是，STR区域通常包含大量较小的变异。

Para_02

hiphase --threads 32 --io-threads 4 --sample-name {sample_id"} --vcf {in_vcf_deepvariant"} --vcf {in_vcf_pbsv"} --vcf {in_vcf_trgt"} --output-vcf {out_vcf_deepvariant"} --output-vcf {out_vcf_pbsv"} --output-vcf {out_vcf_trgt"} --bam {in_bam"} --reference {in_reference"} --summary-file {out_summary"} --blocks-file {out_blocks"} --global-realignment-cputime 300

Parent-of-origin determination

亲本来源鉴定

Para_01

我们使用HiPhase推断的单倍型基因型来确定新生串联重复（TR）扩增和收缩的可能亲本来源。
对于我们在一个孩子中观察到的每个已分型的新生等位基因，我们检查该孩子父母在新生等位基因±500 kb范围内的所有信息性SNV。
我们使用以下标准定义信息性位点：这些位点必须是双等位基因SNV；母亲、父亲和孩子的总测序深度至少为10个读段；母亲、父亲和孩子的Phred标度基因型质量至少为20；孩子的基因型必须是杂合的；并且父母的基因型不能在状态上完全相同。
利用孩子的单倍型分型SNV VCF文件，我们随后确定在信息性位点处孩子携带的REF或ALT等位基因是从母亲还是父亲遗传而来。
例如，如果母亲的基因型是0/0，父亲的基因型是0/1（请注意，父母的基因型不需要进行单倍型分型），而孩子的基因型是1|0，则我们知道孩子的第一个单倍型来自父亲，第二个单倍型来自母亲。
我们对±500 kb区间内的所有信息性位点重复这一过程。
然后我们寻找N个信息性位点，它们满足以下条件：（1）距离新生TR等位基因最近（可以是上游或下游方向）；（2）支持孩子中一致的遗传模式（即都支持孩子两个单倍型来自同一个亲本来源）；以及（3）全部位于同一个HiPhase单倍型区块内（根据HiPhase输出VCF中的PS标签定义）。
最后，我们使用HiPhase生成的分型TR VCF文件来检查新生等位基因是否被分型到孩子的第一个或第二个单倍型中。
我们进一步确认新生等位基因与上述识别的信息性位点具有相同的PS标签，并利用这N个信息性位点来确定新生等位基因所在的单倍型很可能源自母亲还是父亲。

Measuring concordance with orthogonal sequencing technology

使用正交测序技术评估一致性

Para_01

在每一个候选的新生TR等位基因上，我们计算了由TRGT估算出的新生AL与通过Element、ONT或HiFi读段支持的AL之间的一致性。我们将一致性分析限制在仅含单次扩增或收缩的常染色体TR位点（即，我们未分析包含多个独特扩增和/或收缩事件的‘复杂’TR位点）。

Para_02

TRGT 为三联体中的每个成员在常染色体 TR 位点报告两个等位长度（AL）估计值，而 TRGT-denovo 将这两个 AL 中的一个分配为子代中的新生 AL。
在每个 TR 位点，我们计算参考基因组中该位点的长度（以碱基对为单位）与特定个体中两个 AL 中的每一个之间的差异。
我们将 TRGT AL 与参考位点大小之间的差异称为相对 AL。
然后，我们在 BAM 文件中查询包含 Element、Illumina、ONT 或 PacBio HiFi 读段的每个 TR 位点的数据。
使用 pysam 库（https://github.com/pysam-developers/pysam），我们遍历了所有完全跨越 TR 位点且比对质量为 60 的读段。
为了估计相对于参考基因组某个读段中 TR 扩展/收缩的 AL，我们统计了所有与 TR 位点重叠的 CIGAR 操作所关联的核苷酸数量。
例如，一个 Element 读段可能具有如下 CIGAR 字符串：100M2D10M6I32M。
对于与 TR 位点重叠的每个 CIGAR 操作，我们通过 OP * BP 增加计数器，其中 OP 对于‘匹配’操作为 0，对于‘插入’操作为 1，对于‘缺失’操作为 -1，BP 表示给定 CIGAR 操作所涉及的碱基对数目。
因此，在每个 TR 位点，我们在三联体的每个成员中生成了一个净 CIGAR 操作的分布。

Para_03

我们使用这些净CIGAR操作来验证每个子代中的候选新生TR等位基因。
对于每一个新生TR等位基因，我们计算了子代中支持该新生AL（允许Element读段支持新生AL ± 1 bp）的Element读段数量。
接着，我们又计算了该子代父母中支持该新生AL的Element读段数量（同样允许存在一个碱基的误差）。
如果至少有一个Element读段支持子代中的该新生TR AL，并且双亲中没有任何一个Element读段支持该新生TR AL，我们就认为该新生TR得到了验证。

Validating recurrent TR DNMs

验证重复出现的TR DNM

Para_01

为了整理出一份可靠的候选重复新生TR等位基因列表，我们首先整理了一份TR位点列表，这些位点上CEPH 1463家系中的两个或更多个体（在G2、G3或G4代中）存在新生TR等位基因的证据。
对于每一个候选位点，我们进一步要求CEPH 1463系谱中的所有成员都在该位点上进行了TR等位基因的基因分型，并且在该位点至少有10条比对的HiFi读段。
这些过滤条件得到了一个包含49个候选位点的列表，在这些位点上我们观察到了代内或代间重复出现的证据。
我们使用整合基因组浏览器（IGV）以及定制的HiFi CIGAR操作图，对每个位点上的HiFi读段证据进行了目视检查，以确定候选的新生TR等位基因是否合理。

Detection and filtering of de novo SVs

新发结构变异的检测与过滤

Para_01

我们尝试从三个不同的来源获得推定的从头SV（结构变异）。
第一种方法是基于基于读长的调用集报告从头SV（PBSV（v.2.9.0），Sniffles（v.0.12.0），Sawfish（v.2.2））。
第二种方法是通过单倍型基因组组装中鉴定出的变异来报告推定的从头SV。
最后一种方法利用基于单倍型基因组组装构建的泛基因组图谱来报告从头SV。

Assembly-based detection of de novo SVs

基于组装的新生结构变异检测

[ol]- (1) SVPOP8 (v.3.4.0) (https://github.com/EichlerLab/svpop) was used to produce a merged PAV callset across all samples. It merges a single source (single SV caller) across multiple samples. The merge definition used was: nr::ro:szro:exact:match. The samples were provided in this order (G1–G2–G3): NA12889, NA12890, NA12891, NA12892, NA12877, NA12878, NA12879, NA12881, NA12882, NA12883, NA12884, NA12885, NA12886, NA12887. - (2) For each sample in G3, we selected variants unique to that sample alone. - (3) To compare variant calls against the previous generation, SVPOP was used again to do a PBSV/PAV intersection. This involved intersecting the PAV calls for G3 with the PBSV calls for G2, comparing each sample in G3 against each sample in G2. - (4) The callable BED files from PAV, intersections with G2’s PBSV calls, and the list of putative de novo calls went into our validation pipeline. - (5) The pipeline (1) checks if the putative de novo variant was called by PBSV in either parent. (2) Checks if the putative de novo variant is seen in HiFi reads in either parent by running subseq (https://github.com/EichlerLab/subseq). (3) Checks if the variant was in a callable region in either parent. (4) Performs an MSA using DECIPHER of the two haplotypes of the sample, and both parents, in the location of the SV with 1,000 bp flank on either side.

Pangenome graph detection of de novo SVs

泛基因组图谱检测新生结构变异

Para_01

通过使用WFMASH（v.0.13.1-251f4e1）泛基因组比对工具将Verkko组装结果与GRCh38、T2T-CHM13和HG002（v.1.0.1）人类参考基因组进行比对，从而按染色体对组装结果进行划分。
在每组连续序列（contigs）上，我们应用PGGB（v0.6.0-87510bc）构建染色体水平的无偏泛基因组变异图谱，使用的参数为：-s 20k -p 95 -k 47 -V chm13:100000, grch38:100000。
接着我们使用Variation Graph Toolkit（v.1.40.0）相对于T2T-CHM13和GRCh38参考基因组从图谱中检测变异。
然后应用VCFLUB工具（v.0.1.0-26a1f0c）对变异进行分解，保留那些锚定在作为参考的基因组上的顶层气泡中的变异。
再通过VCFWAVE（v.1.0.3）使不同样本间的结构变异（SV）表示趋于一致。
随后，原始VCF文件被用作基于家系过滤潜在新生结构变异（de novo SV）的输入数据。

De novo SV filtering in SV callsets (PGGB, PAV, PBSV, Sniffles, Sawfish)

在SV调用集中进行从头SV过滤（PGGB、PAV、PBSV、Sniffles、Sawfish）

Para_01

使用BCFtools（v.1.17）+fill-tags对新生SV进行过滤，随后在所有样本都有基因型调用的位点上，对联合调用的VCF文件中单例衍生等位基因进行过滤。
通过纳入所有G2/G3家族成员（而不仅限于三联体），我们提高了新生SV的特异性。
我们使用的命令行如下：

Para_02

bcftools view -i 'INFO/AC = 1'VCF FILE"}bcftools +fill-tags -- -t 'all,F_MISSING' | bcftools view -i 'F_MISSING = 0.0' --max-alleles 2 | bcftools view --samples {SAMPLE"| bcftools +fill-tags | bcftools view -i 'INFO/AC = 1' | bcftools view -i '(ILEN < -49 || ILEN > 49)' | bcftools view -i 'QUAL > 49' | vcf2tsv

Para_03

所有在基因组各区域收集的候选新发结构变异（SV）均进一步通过分型基因组组装和长读长比对进行了评估。
更多详细信息见补充说明10。

Extracting donor site of de novo SVA insertion

提取新生SVA插入的供体位点

Para_01

我们首先从NA12887的新生Verkko组装中提取了一个插入的SVA元件（母系单倍型，单倍型1）。
接下来，我们使用minimap2（版本2.24）将这段约3.4 kb长的DNA序列分别与母系和父系的Verkko组装进行比对，所用参数如下所示：

Para_02

minimap2 -x asm20 -c --eqx --secondary=yes {assembly.fasta"} {sva.fasta"} > {output.paf"}

Para_03

根据这些参数，我们报告了该DNA片段的所有位置。我们将一个推定的供体位点定义为母系单倍型中与SVA新生插入近乎完全匹配的比对位置。

Analysis of centromeric regions

着丝粒区域的分析

Para_01

为了从每个基因组组装中完全且准确地识别着丝粒，我们首先使用minimap2（参考文献80）及以下参数：-a --eqx -x asm20 -s 5000 -I 10G -t {threads"}，将通过Verkko16或hifiasm (UL)17生成的基因组组装结果与T2T-CHM13参考基因组进行比对。
然后，我们将全基因组比对结果过滤，仅保留那些比对到T2T-CHM13参考基因组中着丝粒区域的连续序列（contigs）。
我们检查了这些包含着丝粒的连续序列是否跨越了整个着丝粒区域，方法是查看它们是否包含了紧邻着丝粒区域的p臂和q臂序列。
接着，我们通过使用pbmm2（v.1.1.0；https://github.com/PacificBiosciences/pbmm2）及以下命令：align --log-level DEBUG --preset SUBREAD --min-length 5000 -j {threads"}，将来自相同来源基因组的原始PacBio HiFi数据比对到每个全基因组组装结果上，以验证着丝粒区域的组装情况。
随后，我们使用NucFreq77（v.0.1）评估这些组装结果在着丝粒区域内的读段覆盖深度是否均匀。
我们还使用minimap2（v.2.28）将来自同一来源基因组、长度超过30kb的原始ONT数据比对到各个全基因组组装结果上，并利用IGV浏览器93评估这些组装结果在着丝粒区域的读段覆盖深度是否均匀。

Para_02

为了识别每个着丝粒区域内的新生结构变异（SVs）和单核苷酸变异（SNVs），我们首先使用minimap2将每个子代的基因组组装结果与相应亲代的基因组组装结果进行比对，所使用的参数为：-a --eqx -x asm20 -s 5000 -I 10G -t {threads"}。
随后，我们利用生成的PAF文件，通过SVbyEye85（版本0.99.0）识别新生的结构变异和单核苷酸变异，并将结果筛选限定在那些被完全且准确组装的着丝粒区域内。
我们使用NucFreq、Flagger9以及原始的ONT数据对每一个结构变异和单核苷酸变异的检出结果进行了验证，以确保支持每一个变异调用的数据基础可靠。
更多详细信息请参见补充说明9和10。

Analysis of telomeric regions

端粒区域分析

Para_01

我们使用 Tandem Repeats Finder（TRF）91 处理了所有的 G1、G2 和 G3 组装序列，以确定每个组装片段的末端区域是否存在经典的端粒重复序列（p 臂为 CCCTAA；q 臂为 TTAGGG）。
运行 TRF（v.4.09.1）时使用的参数为 '2 7 7 80 10 50 10 -d -h-ngs'，此参数设置适用于较年轻（在此语境中即未退化）的重复序列，该配置已在 RepeatMasker（v.4.1.6）中实现。
随后，使用 minimap2（参考文献 80）（v.2.24）将组装得到的片段比对到 T2T-CHM13 参考基因组上，采用 asm20 预设参数，以确定每个序列的身份信息（即某个特定片段是代表整个参考染色体还是其部分区域，以及是否需要取反向互补链以符合标准表示方式）。
在确认身份后，从外向内爬取 TRF 注释信息（在 p 臂上从 5′ 端开始，在 q 臂上从 3′ 端开始，并根据 minimap2 报告的反向互补情况相应调整），直到遇到经典重复序列为至。
同时保留了非经典散在重复序列的出现情况。

Para_02

此外，将 PacBio HiFi 读长比对到连续序列（contig）上，以评估每个组装结果的每个区域获得了多少条 HiFi 读长的支持（覆盖深度）；对于支持的 HiFi 读长少于五条的远端区域则进行了屏蔽处理。
在所有 G1、G2 和 G3 样本中的非近端染色体末端中，Verkko 组装结果中有 74.2%（所有参与者和单倍型中共有 1,204 个可能的末端，其中 893 个满足条件）终止于一个经典的端粒重复序列（要么跨越连续序列的最起始或终止位置，要么紧邻因覆盖度低而被屏蔽的区域），这些重复序列的中位长度为 5,608 bp（补充表 3）。
此外，在 T2T-CHM13 染色体中，对于其 p 臂和 q 臂端粒末端均被成功恢复的染色体，其中有 64.6%（共 342 条染色体中的 221 条）均由一条从 p 端粒延伸至 q 端粒的单一组装连续序列所代表。

Para_03

G4 hifiasm组装采用了相同的方式进行处理；然而，仅恢复了56.8%的端粒区域（共602个可能中的342个）（补充图3），且经典重复序列的中位长度为4,674 bp（补充表格3；与G1–G3相同），此外连续性明显更差：只有一个染色体（个体G4-200101单倍型1中的第9号染色体）被一个单一的连续序列（h1tg000017l）明确覆盖。

CpG methylation analysis

CpG甲基化分析

Para_01

为了确定每个着丝粒的CpG甲基化状态，我们首先使用Guppy（https://community.nanoporetech.com；v.6.5.7）对原始ONT数据进行碱基识别，采用sup-prom模型和dna_r9.4.1_450bps_modbases_5hmc_5mc_cg_sup_prom.cfg配置文件。
接下来，我们使用minimap2（参考文献80）（v.2.28）将每个样本的ONT数据比对到相应的基因组组装上，使用的参数为：-ax lr:hq -y -t 4 -I 8g。
随后，我们使用modbam2bed（https://github.com/epi2me-labs/modbam2bed）并使用以下参数将生成的BAM文件转换为bedMethyl文件：-e -m 5mC --cpg -t {threads"} {input.bam"} > {output.bed"}。
接下来，我们使用以下命令将bedMethyl文件转换为bedGraph：awk ‘BEGINOFS=“ ”}; {print 2, 11}’ {input.bed"}grep -v “nan” | sort -k1,1 -k2,2n > {output.bedgraph"。
然后使用bedGraphToBigWig（https://www.encodeproject.org/software/bedgraphtobigwig/）将bedGraph转换为bigwig，并使用Integrative Genomics Viewer（版本2.16.0）可视化bigwig文件。
为了确定每个着丝粒中低甲基化区域（称为CDR2，参考文献39）的大小，我们使用了CDR-Finder（https://github.com/arozanski97/CDR-Finder）。
该工具首先将bedGraph按5 kb窗口分箱，计算包含α-卫星序列（由RepeatMasker，参考文献99（v.4.1.0）确定）的窗口内的CpG甲基化频率中位数。
选择甲基化频率低于该区域中位数的窗口，将连续窗口合并为更大的窗口，筛选出长度大于50 kb的合并窗口，并报告这些窗口的位置。

Y-chromosomal analysis

Y染色体分析

Construction and dating of Y phylogeny

Y染色体系统发育树的构建与年代测定

Para_01

总共58个样本的Y染色体系统发育树的构建和定年，包括本研究中的14个家系男性以及此前已发表基于长读长的Y染色体组装数据的44个个体，其过程如前所述52。
简而言之，使用此前定义为适用于短读长测序（SRS）的约10.4 Mb Y染色体序列，从14名家系男性的Illumina高覆盖度数据中识别出所有位点。
使用BCFtools101,102（版本1.16），设置最低碱基质量20、比对质量20以及倍性为1进行分析。
移除了在插入缺失（indel）呼叫位点前后5 bp范围内的单核苷酸变异（SNV），并过滤掉所有的插入缺失位点。
随后对所有位点进行进一步过滤，要求最小读深为3，并且支持所检测基因型的读段比例需≥85%。
使用BCFtools将所得VCF文件与参考文献52中44个个体经过类似过滤的VCF文件合并，
然后利用VCFtools103（版本0.1.16）去除缺失率≥5%的位点，即在58个样本中超过3个样本缺失的位点。
经过过滤后，总共保留了10,404,104个位点，其中包括13,443个多态位点。

Para_02

每个样本的Y单倍群组均如先前所述进行预测104，并与国际基因谱系学会命名法（ISOGG；https://isogg.org；版本15.73）相对应。
使用BEAST105（版本1.10.4）中实现的一种基于共祖（coalescence）的方法来估计内部节点的年龄。
采用RAxML106（版本8.2.10）并结合GTRGAMMA替换模型构建用于BEAST分析的起始最大似然系统发育树。
马尔可夫链蒙特卡洛抽样基于2亿次迭代，每1000次迭代记录一次，前10%的迭代结果被作为预烧期（burn-in）丢弃。
采用恒定大小的共祖树先验模型、考虑位点异质性（gamma分布）的GTR替换模型以及严格时钟模型（其替代速率基于95%置信区间呈正态分布：每年每碱基对为0.76 × 10−9（95% CI = 0.67 × 10−9–0.86 × 10−9）单核苷酸突变）107。
利用Tree-Annotator软件（版本1.10.4）生成汇总树，并通过FigTree软件（版本1.4.4）进行可视化。

Identification of sex-chromosome contigs

性染色体重叠群的鉴定

Para_01

对Y染色体上的DNMs进行了详细分析，研究对象为7名男性个体（R1b1a-Z302 Y单倍群，G1-NA12889、G2-NA12877、G3-NA12882、G3-NA12883、G3-NA12884和G3-NA12886），这些个体生成了分相的Verkko组装序列。
包含X染色体和Y染色体序列的重叠群（contigs）是从全基因组组装数据中鉴定并提取出来的，具体方法如前所述。
此外，通过使用minimap2（参考文献80，版本2.26）将T2T-CHM13参考基因组中的相应序列比对到G1祖母NA12890和G2母亲NA12878的基因组组装结果中，从而鉴定了这些组装中的假常染色体区域。

Annotation of Y-chromosomal subregions

Y染色体亚区段的注释

Para_01

Verkko组装的Y染色体亚区注释是使用GRCh38和T2T-CHM13 Y参考序列进行的，方法如先前所述。
为了进行Y染色体亚区注释，使用RepeatMasker（版本4.1.2-p1）及其默认参数识别着丝粒α卫星重复序列。
Yq12重复序列的注释是使用HMMER（版本3.3.2dev）完成的，所用的序列为已发表的DYZ1、DYZ2、DYZ18、2k7bp和3k1bp序列。
随后通过人工检查重复单元的方向及其相互之间的距离对结果进行了验证。
为比较Y染色体序列生成点阵图（dot plots），我们使用了Gepard（版本2.0）软件。

Detection and validation of DNMs

DNMs的检测与验证

Para_01

人类Y染色体在重复区域的大小和组成方面存在广泛变异，其中包括T2T-CHM13 Y（单倍群J1a-L816）以及家系中七位男性携带的R1b1a-Z302单倍群Y染色体（补充说明6）。
因此，G1代祖父NA12889的Y染色体组装被用作检测DNM（新发突变）的参考序列。
这些DNM是通过对五名G2代（NA12877）和G3代（NA12882、NA12883、NA12884和NA12886）男性的Y染色体组装数据使用Dipcall软件（版本0.3）调用获得的，参数采用适用于男性样本的默认设置。
变异仅从MSY（男性特异性区域）中识别得出，也就是说，伪常染色体区域未包含在本次分析中。
所有识别出的变异都经过如下过滤：任何在参考或查询组装中被Flagger或NucFreq标记的区域中的变异都被排除在外。

Para_02

对于SNVs，最终过滤后的变异位点均得到了100% HiFi读长的支持（即在子代中没有任何读长支持参考等位基因，或在父亲中支持替代等位基因），并且检查了比对到参考序列和各个个体组装结果中的ONT读长以确认其支持情况。

Para_03

对于小片段插入缺失（≤50 bp），分析中排除了同聚物区域，其余的插入缺失通过读取数据（HiFi、ONT、Illumina）进行验证，具体方法如下。
使用subseq工具（https://github.com/EichlerLab/subseq）从所有样本中提取覆盖插入缺失位点及其上下游各150 bp序列的个体读段，这些读段被映射到参考基因组（G1 NA12889 Y组装结果）。
随后使用MAFFT110,111（v.7.508）对提取的读段进行比对，采用默认参数。
所有比对结果均经过人工检查，若HiFi数据中存在两个或以上支持参考等位基因的读段以及一个或以上支持替代等位基因的读段，则该插入缺失会被排除。
所有最终的SNV和插入缺失结果，在能够唯一映射到目标区域的情况下，都得到了Illumina和Element读段数据的支持。

Para_04

对于所有的结构变异（SV）检测，我们可视化了参考等位基因和替代等位基因的HiFi读长深度，并去除了那些读长深度变化显著且与单核苷酸变异（SNV）簇重合的区域中的SV，这些SNV中超过10%的读长支持替代等位基因。同时，我们检查了HiFi和ONT读长在参考基因组和个体组装结果中的比对支持情况。

Para_05

对于所有的变异位点，均确认了其在世代间传递的预期结果具有一致性。
此外，还利用三个G4代雄性个体（编号200101、200102和200105）的HiFi数据，验证了所发现变异位点的支持情况。

Y-chromosomal DNM rate calculation

Y染色体DNM率计算

Para_01

基于每个男性个体Y染色体组装的可及区域（即被Flagger和/或NucFreq标记的任何区域均已被去除），计算了这五个男性个体各自的基于组装的DNM率。

Mobile element analysis

移动元件分析

Para_01

使用 xTea（版本 0.1.9）对 PacBio HiFi 读段进行移动元件分析。
通过 IGV 可视化 xTea 鉴定到的潜在非参考移动元件插入（MEI），以确保这些插入在测序读段中可识别，并判断其中是否存在新生事件。
使用 BEDTools，我们将非参考插入与来自 T2T-CHM13 的内含子、外显子、5'-UTR 和 3'-UTR 进行比对交叉分析。
为了鉴定非参考 LINE-1 插入的潜在来源元件，我们使用 BLAT 在 T2T-CHM13 参考基因组中寻找最佳匹配的插入序列。
如果参考基因组中存在多个得分相同的匹配项，则不判定某个来源元件。
代表已知 Alu、L1 和 SVA 子类的 MEI 序列来自之前的研究、Dfam 数据库和 UCSC 基因组浏览器。
每个 MEI 类别的参考序列和新发现序列被合并到类别特定的文件中。
将序列方向统一为正链方向。
去除高度截断的序列。
使用 MUSCLE（版本 3.8.31）比对 MEI 序列。
采用 Kimura 两参数法计算 MEI 序列间的成对距离，并将其转换为相关系数。
通过对成对相关系数矩阵进行特征值分解获得主成分。
绘制前三个主成分，以可视化非参考 MEI 与已知 MEI 亚家族序列之间的关系。

Reporting summary

报告摘要

Para_01

有关研究设计的进一步信息，请参阅与本文关联的《Nature Portfolio》报告摘要。

Data availability

Para_01

家族中28位成员的所有基础数据均可作为AWS开放数据计划、欧洲核苷酸档案库（ENA）或dbGaP的一部分获取。对于其中23位家庭成员（G1-GM12889、G1-GM12890、G1-GM12891、G1-GM12892、G2-GM12877、G2-GM12878、G3-GM12879、G3-GM12881、G3-GM12882、G3-GM12885、G3-GM12886、G3-200080-配偶、G4-200081、G4-200082、G4-200084、G4-200085、G4-200086、G4-200087、G3-200100-配偶、G4-200101、G4-200102、G4-200104 和 G4-200106），他们已同意其数据像千人基因组计划样本一样公开访问，以支持新技术开发、人类变异研究、DNA生物学研究以及健康与疾病研究，这些数据可通过AWS开放数据计划（s3://platinum-pedigree-data/）以及欧洲核苷酸档案库（BioProject: PRJEB86317）获得。有关如何访问这些数据的具体细节可在GitHub上查阅（https://github.com/Platinum-Pedigree-Consortium/Platinum-Pedigree-Datasets）。
另外五位未同意开放访问的家庭成员（G3-NA12883、G3-NA12884、G3-NA12887、G4-200103 和 G4-200105）的比对测序数据和组装数据可在dbGaP中获取（phs003793.v1.p1；Platinum系谱联盟LRS）。这些数据也包括整个家庭（28个成员）的变异调用。
TR目录可在Zenodo上获取（https://doi.org/10.5281/zenodo.13178746）。
一个密切相关的R1b单倍群样本HG00731的Y染色体组装数据从人类基因组结构变异联盟IGSR网站下载（https://ftp.1000genomes.ebi.ac.uk/vol1/ftp/data_collections/HGSVC3/working/20230927_verkko_batch2/assemblies/HG00731/）。
本研究所使用的参考基因组及其注释见补充表14。