超低深度全基因组测序数据中基因型填补工具的综合评估

摘要

在不同测序深度和样本量下，使用CKB和EAS参考面板对三种填补工具——STITCH、QUILT2和GLIMPSE2——进行了基准测试。研究结果为在超低深度测序数据中选择合适的基因型填补策略提供了实用框架，有助于开展大规模人群基因组学和复杂性状研究。

摘要

超低深度测序（ULDS）是一种适用于大规模基因组研究的经济高效策略，但其应用效果依赖于准确的基因型填补。本研究利用中国嘉道理生物样本库（CKB）和千人基因组计划（1KGP）东亚人群（EAS）参考面板，评估了三种填补工具——STITCH、QUILT2 和 GLIMPSE2——在不同测序深度和样本量下的表现。研究揭示了关键的性能差异：样本量敏感性方面，STITCH 的准确性随样本量增大显著提升，而 QUILT2 和 GLIMPSE2 对样本量的依赖性极小；参考面板优化方面，人群特异性的 CKB 面板显著提高了 QUILT2 和 GLIMPSE2 的准确性，但对主要依赖内部单倍型推断的 STITCH 影响甚微；深度阈值方面，所有工具在中等测序深度（≥ 0.5x）下均达到稳健的准确性，但在超低深度（≤ 0.1x）时，STITCH 表现明显不足。GLIMPSE2 结合 CKB 实现了最高的整体准确性，而 QUILT2 在精度与计算效率之间实现了良好平衡。对于无创产前检测（NIPT）数据，GLIMPSE2+CKB 仍能保持足够的准确性以支持下游分析。本研究提出一个决策框架，优先推荐使用人群匹配的参考面板和适应测序深度的工具，为在不同研究场景下优化 ULDS-WGS 提供了可操作的指导建议。这些发现弥合了方法学进展与实际应用之间的差距，使基因组研究能够在不牺牲数据质量的前提下实现经济高效的规模化扩展。

引言

超低深度测序（ULDS）是指测序深度低于1x的测序方法，因其成本低、基因组覆盖范围广以及适用于多种样本类型而受到广泛关注。ULDS已在无创产前检测（NIPT）¹、癌症监测²和染色体拷贝数变异（CNV）检测^3,4等临床应用中展现出重要价值。除临床诊断外，随着测序成本的持续降低和生物信息学的快速发展，ULDS在群体基因组学和复杂性状研究中正发挥着日益重要的作用。通过将ULDS数据与大规模群体单倍型参考面板相结合，基因型填补技术能够恢复个体水平的全基因组变异信息。因此，ULDS已成为传统单核苷酸多态性（SNP）芯片和高深度全基因组测序（WGS）的一种经济高效的替代方案，尤其适用于全基因组关联分析（GWAS）和群体结构分析等大规模研究⁵。

先前的研究已证明，利用无创产前检测（NIPT）测序数据开展多种遗传学研究是可行的，包括变异检测、群体历史重建、病毒感染模式推断以及全基因组关联分析（GWAS）⁶。

尽管具有这些优势，超低深度测序（ULDS）数据极度稀疏的特性仍带来了独特的挑战。在变异位点层面，许多位点完全未被检测到，或每个个体仅由单个等位基因代表，导致下游分析所需的数据质量不足。因此，基因型填补至关重要，它利用大型参考面板（例如 1000 Genomes⁷ 或人群特异性资源）中的单倍型结构，通过统计学方法推断缺失或不确定的基因型。先前的研究表明，从无创产前检测（NIPT）数据进行填补可实现高准确性，并在全基因组关联分析（GWAS）中识别性状相关变异时保持强大的统计功效⁸。通过使用 STITCH⁹ 算法，对中国20,900名孕妇队列的NIPT数据（平均深度约0.15x）成功进行了基因型填补，进而鉴定了与妊娠相关的位点。填补后的基因型与高深度全基因组测序（WGS）数据在GWAS结果中表现出高度一致性（Pearson R² > 0.8）¹⁰。

基于超低深度测序（ULDS）的分析的成功在很大程度上取决于填补（imputation）的准确性，而填补准确性又受到测序深度、参考面板的质量与人群匹配度、填补算法的性能、样本量以及等位基因频率谱的影响¹¹。在这些因素中，参考面板的选择是决定填补准确性的主要因素。常用的参考面板包括具有全球代表性的资源，如1000基因组计划（1KGP）⁷、TOPMed¹²和单倍型参考联盟（HRC）¹³，以及越来越多可用的人群或地区特异性面板，例如新加坡一万个基因组（SG10K）¹⁴和中国嘉道理生物样本库（CKB）¹⁵。影响填补性能的另一个关键因素是算法的选择。已有多种工具被开发出来，以应对低深度测序带来的独特挑战，显著推动了填补方法在大规模遗传学研究中的实际应用。尽管Beagle（v5+）¹⁶、Minimac4¹⁷和IMPUTE5¹¹等填补方法被广泛用于SNP芯片和中高深度全基因组测序（WGS）数据，但它们在ULDS条件下的表现往往不够理想。近年来，一些专用工具被开发出来以应对这些挑战。STITCH⁹可直接从低深度测序读段中推断单倍型，特别适用于大规模同质性队列。QUILT2¹⁸采用压缩单倍型库和局部似然模型，能够在使用大规模参考面板时实现高效的填补，并在产前基因组学中具有独特应用价值。GLIMPSE2¹⁹是原始GLIMPSE框架的扩展版本，在准确性和计算效率方面均实现了进一步提升。

尽管这些工具代表了重要的技术进步，但它们在不同实验设计（例如测序深度、队列规模和参考面板选择）下的相对性能尚未得到系统评估，导致研究人员在选择最合适的策略时缺乏明确指导。为填补这一空白，本研究对三种广泛使用的超低深度测序（ULDS）基因型填补工具——STITCH、QUILT2 和 GLIMPSE2——在多种测序深度和样本量条件下进行了系统性基准测试。评估采用了两个与中国人群高度相关的东亚参考面板。研究结果表明，在测序深度 ≥0.5x 时，ULDS 填补通常具有较高的可靠性；而测序深度 <0.1x 时，则需要显著更大的队列规模才能达到可接受的准确性。参考面板的选择应根据具体研究背景进行优化，使用与目标人群匹配的面板（如 CKB）可提高填补准确性。此外，这些方法可直接应用于大规模人群研究和无创产前检测（NIPT）中产生的超低深度测序数据。因此，本研究为基于 ULDS 的研究提供了工具选择的实用框架，为未来在群体遗传学和复杂性状分析中的应用提供了方法学指导。

方案

所有参与者在参与前均签署了书面知情同意书。涉及高深度全基因组测序（WGS）数据的研究已通过华大基因机构审查委员会（BGI-IRB 23058-T2）的审查和批准，并已获得中国人类遗传资源管理办公室关于人类遗传资源采集的批准（[2023] CJ0262）。涉及来自无创产前检测（NIPT）的超低深度测序（ULDS）数据的研究已获得武汉市儿童医院机构审查委员会（2021R062）、华大基因机构审查委员会（BGI-IRB 21088）的批准，并获得了中国人类遗传资源管理办公室的额外批准（[2021] CJ2002）。

注意：本研究包含两种类型的全基因组测序（WGS）数据。第一类为来自深圳自然人群队列中500名个体血液样本的高深度WGS数据（30x），用于构建高质量的基准数据集，并进行后续的降采样及准确性评估。第二类为来自武汉地区10,000名孕妇无创产前检测（NIPT）的超低深度测序（ULDS）数据。

1. 高深度全基因组测序数据

1. 在获得知情同意后，从普通人群队列中采集500份外周血样本（每份5 mL）。将样本储存于EDTA管中，并在2-8 °C条件下运输。
2. 在4 °C条件下，以1,600 x g 离心10分钟，分离血浆和白膜层。小心收集白膜层，并在-80 °C保存，直至进行DNA提取。
3. 使用基于磁珠的试剂盒，按照制造商说明书从白膜层中提取基因组DNA。
4. 采用荧光定量法测定DNA浓度，并通过琼脂糖凝胶电泳评估DNA完整性。选择总DNA产量≥1 µg、浓度≥12.5 ng/µL、片段长度>20 kb且无明显降解的样本用于文库构建。
5. 通过超声处理将80-200 ng高质量基因组DNA打断至平均大小为350-400 bp。
6. 在20 °C下进行末端修复30分钟，20 °C下接头连接15分钟，37 °C下环化30分钟，以构建无PCR扩增的文库。利用滚环扩增（RCA）生成DNA纳米球（DNBs）。在DNBSEQ测序平台上对双端文库（PE100，读长100 bp）进行测序，目标测序深度约为30x（平均每样本100 Gb）。将原始测序读段以FASTQ格式存储，用于后续分析。
   注意：所有来源于人类的样本均应在生物安全二级（BSL-2）实验室条件下操作。避免反复冻融，以防DNA降解。血液来源材料应作为生物危害废物处理；化学试剂应按照机构规定的危险废物指南进行处置。

2. 超低深度无创产前检测数据（约0.1x WGS）

1. 采集10,000份孕妇血液样本（每份5 mL），用于常规无创产前检测（NIPT）。使用EDTA管采集，并在2-8 °C条件下运输；采集后8小时内处理血浆。
2. 对于稳定游离DNA的采血管（K管或G管），使用温控运输设备在6-35 °C条件下运输，并在采集后96小时内按照制造商的标准操作程序进行处理。
3. 在4 °C条件下以1,600 × g离心10分钟，分离血浆。使用移液器小心吸取上层血浆层，避免扰动白膜层或细胞沉淀，并转移至新管中。将所得血浆再次于4 °C条件下以16,000 × g离心10分钟，以去除任何残留的细胞或碎片。小心将澄清的上清液（无细胞血浆）转移至新管中，用于DNA提取。
4. 使用核酸提取试剂盒从血浆中提取游离DNA（cfDNA）。在20 °C下进行末端修复30分钟，20 °C下接头连接15分钟，PCR扩增（12个循环：98 °C变性10秒，60 °C退火30秒，72 °C延伸30秒）。
5. 纯化PCR产物，并在37 °C下环化文库30分钟。通过滚环扩增（RCA）生成DNB。在BGISEQ-500平台上对单端文库（SE35，读长35 bp）进行测序。将原始测序数据以FASTQ格式存储。
   注意：在BSL-2条件下将血浆样本视为潜在感染性材料处理。尽量减少冻融循环以降低cfDNA降解。血浆废弃物及塑料耗材应作为生物危害物处理。

3. 数据预处理流程

1. 为系统评估不同测序深度下基因型填补工具的性能，对原始的高深度全基因组测序数据（30x）和超低深度无创产前检测（NIPT）数据（<0.1x）执行标准化的预处理流程，包括模拟降采样、质量控制、读段比对、重复序列去除以及碱基质量评分重校准（BQSR）。
   注意：从本步骤开始至填补准确性评估为止的流程构成了本研究的主要实验方案（图1）。具体代码见补充文件1。
2. 降采样
   1. 从原始的30x高深度测序样本中生成一系列降采样数据集。采用以下两种策略，以真实模拟NIPT数据的测序特征。
   2. 随机子采样：使用seqtk v1.5（https://github.com/lh3/seqtk）并设定固定随机种子为100，生成四个低深度水平的数据（0.05x、0.1x、0.5x和1.0x）。
   3. 模拟NIPT样读段结构：仅保留每对双端读段中的第一端读段（R1），并使用seqtk trimfq -L 35将所有保留的读段截断至35 bp，以符合超低深度NIPT测序典型的单端、短读段特征。
3. 质量控制
   1. 使用fastp v0.23.4²⁰处理所有原始FASTQ文件。采用以下参数：--qualified_quality_phred=5（碱基质量阈值），--unqualified_percent_limit=50（允许的低质量碱基最大百分比），--n_base_limit=10（每条读段中N碱基的最大数量），并使用自定义接头序列进行接头去除：--adapter_sequence=AAGTCGGAGGCCAAGCGGTCTTAG
      GAAGACAA（R1）和--adapter_sequence_r2=AAGTCGGATCGTAGCC
      ATGTCGTTCTGTGAG
      CCAAGGAGTTG（R2）。
   2. 禁用Poly-G尾部修剪（--disable_trim_poly_g），并为每个样本生成JSON和HTML格式的报告。
4. 比对与重复序列去除
   1. 使用BWA v0.7.16a-r1181²²将高质量读段比对至人类参考基因组GRCh38（hg38）²¹。
   2. 采用aln算法进行比对（-e 10 -t 4 -i 5 -q 0），随后使用samse进行包含读段组信息的单端比对。
   3. 将生成的SAM文件转换为BAM格式，进行排序（samtools sort -@ 8），并使用SAMtools v1.3²³（samtools rmdup）去除重复序列。对所有BAM文件建立索引。
5. 碱基质量评分重校准（BQSR）
   1. 使用GATK v4.0.4.0²⁴执行BQSR。使用三个高置信度变异数据集训练重校准模型：dbSNP build 146²⁵、Mills和1000G金标准indel数据集²¹，以及GATK资源包²⁴中针对GRCh38的已知indel文件。所用资源包文件参考GATK官方示例（https://github.com/gatk-workflows/gatk4-data-processing/blob/master/processing-for-variant-discovery-gatk4.hg38.wgs.inputs.json）。共下载三个文件及其对应的索引文件。
   2. 依次运行BaseRecalibrator和ApplyBQSR以生成重校准后的BAM文件。使用SAMtools v1.3对所有BAM文件建立索引。
      注意：所有模拟数据集均经过相同的预处理步骤——降采样、质量控制、比对、重复序列去除和BQSR——以确保后续填补性能评估的一致性和可比性。

4. 基因型填补

1. 数据准备
   1. 填补数据集设置：根据下述方法构建多个评估数据集，以系统性地在不同测序深度和样本量条件下对基因型填补工具进行基准测试。基于上述不同的样本量和测序深度，共形成九种组合。输入文件包括经过质量控制后上述九个子集的测序数据 BAM 文件列表（bamlist.txt），存储于相应的 bamlist.txt 文件中。其他输入文件包括人类参考基因组（GRCh3821）以及来自 1000 Genomes Project⁷ 的遗传图谱。
      1. 稀释后的高深度全基因组测序数据：从 500 名以 30x 深度测序的个体中随机选取两个子集（200 和 500 个样本）。将每个子集进一步稀释至四个深度（1x、0.5x、0.1x 和 0.05x），生成八种实验条件。
      2. 基于无创产前检测（NIPT）的超低深度测序（ULDS）数据集：将 10,000 个超低深度 NIPT 样本（平均深度为 0.102x，图 2）与 50 个高深度样本稀释至 0.1x 后的数据合并。
   2. 分析区域设定：将所有分析限制在染色体 1 的区域 chr1:150,500,000-160,500,000（10 Mb）内，并为填补设置 500 kb 的缓冲区，以确保不同工具之间的直接可比性。
   3. 参考面板选择：使用两个参考面板（表 1）：CKB 面板，基于中国人群数据构建；1KGP-EAS 面板，源自 1000 Genomes Project 的东亚人群子集。
      注：CKB 参考面板¹⁵ 利用中国嘉道理生物样本库（China Kadoorie Biobank）这一大型前瞻性队列研究中 9,964 名中国成年人的高深度（约 15x）全基因组测序数据构建而成。这些样本来源于自然人群，表型偏倚最小，具有同质的汉族血统和一致的人群结构，特别适用于中国人群的基因型填补。Yu 等人¹⁵ 的研究表明，在一项真实表型的身高全基因组关联分析（GWAS）中，使用 CKB 面板进行填补使检测到的 SNP 数量增加了两倍，全基因组显著变异的数量增加了一倍。1000 Genomes Project（1KGP）^7,26 是目前最广泛使用的基因组参考资源，其东亚（EAS）第 3 阶段子集包含 585 名个体。该子集涵盖五个东亚人群，测序深度约为 30x，包括北京汉族（CHB）、南方汉族（CHS）、西双版纳傣族（CDX）、越南胡志明市京族（KHV）以及东京日本人（JPT）。
2. 填补工具
   1. 评估三种填补算法，选择依据为其不同的建模策略及对超低深度测序（ULDS）数据的适用性。
      1. STITCH：STITCH（v1.6.6）是一种无需参考的基于单倍型的填补算法，可选择性地整合外部参考单倍型。输入文件包括 BAM 列表和人类参考基因组（GRCh38）。在进行基于参考的填补时，需准备参考面板文件（hap/legend/pos）。包含以下关键参数：method=diploid，buffer=500 kb，K=10 个祖先单倍型，nGen=4x 样本量/K（遵循 STITCH 文档建议）。生成包含所有个体每个 SNP 基因型剂量的输出文件。
         注：根据 STITCH 官方文档，K 表示模型中祖先单倍型的数量。较大的 K 值可提高大样本和高覆盖度下的填补准确性，但也会增加计算时间，且在低覆盖度下准确性可能下降。
      2. QUILT2：QUILT2 采用一种针对 ULDS 数据优化的贝叶斯参考引导方法。使用提供的 prepare_reference 脚本进行参考面板准备，指定遗传图谱和区域坐标。以与 STITCH 相同的缓冲区大小（500 kb）和 nGen 设置运行二倍体模式填补，以确保可比性。
      3. GLIMPSE2：GLIMPSE2 是一种基于隐马尔可夫模型（HMM）的参考驱动填补工具，专为大规模和极低深度测序数据集设计。使用 GLIMPSE2_phase_static 执行填补，指定输入 BAM 列表、人类参考 VCF 面板、遗传图谱、输入区域 = chr1:150,000,000-161,000,000，输出区域 = chr1:150,500,000-160,500,000（以维持 500 kb 缓冲区）。此处输入的 BAM 列表文件需包含两列：第一列为 BAM 文件路径，第二列为样本名称。若未提供第二列，则每个 BAM 文件名将作为输出 VCF 文件中的样本名称。

5. 填补准确度评估

1. 真实集定义
   1. 选择50个测序深度为30x的个体作为真实数据集。将这些样本纳入实验中的下采样条件中，以确保可比性。
   2. 使用以下流程对真实集进行变异检测：采用SOAPnuke²⁷进行质控，BWA用于序列比对，Picard用于标记重复序列，GATK v4.0.4.0的BQSR和HaplotypeCaller进行变异检测，DPGT（https://github.com/BGI-flexlab/DPGT）进行联合检测，BCFtools v1.11²³进行变异质量过滤。
   3. 仅保留高置信度的PASS变异，生成基准VCF文件。评估范围限制在chr1:150,500,000-160,500,000区域，与填补后的数据集保持一致。
2. 数据标准化与过滤
   1. 使用PLINK2.0²⁸处理所有工具生成的填补VCF文件。提取剂量数据并转换为pgen格式。
   2. 对SNP水平进行质量控制，应用以下过滤条件：次要等位基因频率（MAF）≥ 0.05（--maf 0.05），哈迪-温伯格平衡（HWE）p值 ≥ 1e-6（--hwe 1e-6），仅保留双等位基因SNP（--max-alleles 2）。
   3. 将通过质控的变异导出为traw格式，用于后续比较。
3. 准确性评估指标
   1. 将每个SNP的填补剂量与真实剂量进行比较。逐个SNP计算皮尔逊相关系数（R），仅保留两个数据集中共有的位点，并计算所有评估SNP的平方相关系数均值（R²），以量化每种条件下的整体填补准确性。
   2. 使用该准确性指标来反映填补基因型剂量估计值与真实基因型之间的一致性。

结果

样本量对填补准确度的影响
将样本量从 N = 200 增加到 N = 500 提高了 STITCH 的填补准确度，尤其是在低测序深度条件下。例如，在 CKB 参考面板 1x 测序深度下，STITCH 的 R²> 从 0.882（N=200）提升至 0.916（N=500），提高了 3.4% (图 3；补充文件 2)。同样，在 0.5x 测序深度下，其准确度从 0.800 提高到 0.868（ΔR²> = 8.5%）。相比之下，QUILT2 和 GLIMPSE2 对样本量变化的敏感性极低，在所有测试条件下 R²> 波动均小于 0.5%。例如，QUILT2 在 CKB 面板 1x 测序深度下表现稳定（N = 200 时 R²> = 0.970，N = 500 时为 0.971）。这表明 STITCH 更能从更大的样本量中获益，这与其基于隐马尔可夫模型（HMM）的单...

讨论

本研究系统评估了三种广泛使用的基因型填补工具在超低深度测序（ULDS）中的表现，以高深度全基因组测序（WGS）作为金标准。该研究的关键方法学优势在于采用了统一的预处理流程——包括序列比对、质量控制和碱基质量评分重校准——从而最大限度地减少批次效应，并确保不同工具和条件之间的可比性。通过对深度测序样本进行降采样，在受控条件下模拟出超低深度数据，从而为性能评估提供了客观的基准框架。将分析限定于1号染色体上一个定义明确的10 Mb基因组区间，既保证了计算的可行性，又保留了足够的变异位点密度。对STITCH、QUILT2和GLIMPSE2的比较评估揭示了不同填补策略各自的优势与局限性。

测序深度和样本量对填补准确性具有可预测的影响。当深度低于约0.1x时，准确性急剧下降，反映出稀疏的读段数据无法提供足够的信息以可靠地重建单倍型。更大的样本量通过更准确地推断群体单倍型结构，缓解了这一限制，这与理论预期一致：更多的样本可提高填补性能²⁹。结果表明，即使在较小的样本量下，0.5x及以上的测序深度仍能对常见变异保持可接受的准确性；而低于0.1x的测序...

材料

本文使用的材料清单

姓名	公司	目录编号	评论
数据
10,000 例 NIPT 低深度样本	本文		用于填补基准测试的超低深度全基因组测序数据。
500 例高深度 WGS 样本	本文		30× 高深度 WGS，作为金标准/真实数据集。
参考面板
1KGP-EAS 参考面板	千人基因组计划（东亚）		千人基因组计划中用于东亚血统特异性基因型填补的子集。
CKB 参考面板	中国嘉道理生物样本库		用于基因型填补的定制人群特异性面板。
软件与算法
BCFtools v1.11	GitHub（samtools/bcftools）		用于合并和排序染色体水平结果，以及变异位点过滤。
GATK 4.0.4.0 工具集中的 BQSR	布罗德研究所		用于碱基质量评分重校准（BQSR）。
BWA-MEM .7.16a-r1181	Heng Li / GitHub		用于将原始测序读段比对至 GRCh38 参考基因组。
DPGT（分布式群体遗传学工具）	BGI		一种分布式群体遗传学分析工具，支持对数百万 WGS 样本进行联合检测。可访问 [GitHub - BGI-flexlab/DPGT](https://github.com/BGI-flexlab/DPGT) 获取。
fastp.0.23.4	开源软件（Chen 等，2018）		用于质量控制和接头序列修剪。
GLIMPSE2	牛津大学		适用于低覆盖度 WGS 的快速基因型分型与填补工具。
分析所用原始代码	本文		补充文件 1 分析所用原始代码
Picard 工具包	布罗德研究所		用于标记重复序列和文件格式转换。
Plink 2.0	C. Chang, S. Purcell / 布罗德研究所		用于基因型格式转换和关联分析。
Python 3.8	Python 软件基金会		用于脚本编写、自动化处理和数据分析。
QUILT2	牛津大数据研究所		基于隐马尔可夫模型（HMM）并使用外部参考面板进行基因型填补。
R 4.1.3	R 基金会		用于运行 STITCH、QUILT2 以及绘图和统计分析。
SAMtools v1.3	GitHub（samtools/samtools）		用于操作 SAM/BAM 文件。
Seqtk-1.5	GitHub（lh3/seqtk）		用于处理 FASTA/Q 格式序列的工具包。可访问 [GitHub - lh3/seqtk](https://github.com/lh3/seqtk) 获取。
SOAPnuke	BGI		用于高通量测序数据的质量控制和过滤。
STITCH v1.6.6	牛津大学		专为超低覆盖度测序优化的基因型填补工具。
tabix	GitHub（samtools/tabix）		用于索引和查询经 bgzip 压缩的 VCF 文件。
其他材料
GATK 捆绑文件	GATK		可访问 [https://github.com/gatk-workflows/gatk4-data-processing/blob/master/processing-for-variant-discovery-gatk4.hg38.wgs.inputs.json] 获取。
1000G 的遗传图谱（GRCh38）	牛津大学 / 千人基因组计划		基因型分型/填补工具所必需的输入文件。
GRCh38	基因组参考联盟		用于测序读段比对和变异检测。