基于群体遗传学的二倍体基因组浅测序无参考遗传多样性估计方法DipSkmer

《PLOS Computational Biology》：DipSkmer: Reference-free population genomics with diploid genome skims

　　

生态学家和保护生物学家依赖遗传多样性作为关键的基本生物多样性变量（essential biodiversity variable，EBV），用以追踪种群健康状况和动态变化，并利用群体规模突变参数θ作为多样性的核心指标。尽管全基因组测序（whole-genome sequencing，wgs）日益经济可行，但在高质量组装基因组能够代表全部生命多样性之前仍需相当长的时间；即便届时实现，持续监测仍需要在不同时间和空间进行重复采样。相比之下，基因组浅测序（genome skimming，即低覆盖度、短读长wgs）具有极高的成本效益，但由于覆盖度过低而难以进行组装和可靠的错误校正，分析颇具挑战。成熟方法如Mash基于k-mer集合的Jaccard相似性，利用草图技术（sketching techniques）估计成对基因组距离。部分方法如Skmer还额外对低覆盖度的影响进行建模。这些方法已成功应用于免组装物种鉴定和系统发育学；然而，其在群体遗传学中的应用仍受限。除其他局限外，这些方法隐含地将基因组视为单倍体，杂合性会混淆二倍体生物基因组距离的估计。在本研究中，研究人员通过若干技术进展解决了这一问题。首先，研究人员利用溯祖理论（coalescent theory）在近中性Kingman溯祖模型下，从数学上推导了两个二倍体样本的k-mer集合之间Jaccard指数随θ变化的关系。其次，研究人员对覆盖度和测序误差的影响进行建模，从而从Jaccard指数以及若干辅助变量（同样从基因组浅测序数据中估计）推导出θ的估计值。由此产生的方法DipSkmer能够对二倍体样本的覆盖度、测序误差和θ进行更准确的估计。对模拟数据集和实证数据集的分析均表明，对于二倍体样本和低遗传距离（例如<2%），DipSkmer的估计优于现有的免比对方法（如Mash和Skmer），并接近基于参考基因组和比对的工具ANGSD。

DipSkmer：面向二倍体基因组浅测序数据的无参考群体基因组学方法

一、研究背景与问题

遗传多样性作为基本生物多样性变量（essential biodiversity variable，EBV），是生态学家和保护生物学家监测物种健康状况与种群动态的核心指标。核苷酸多样性，即两个生物体基因组中差异位点所占比例，被广泛应用于系统发育学、微生物组分析、比较基因组学和群体遗传学。然而，全基因组测序（whole-genome sequencing，wgs）虽然日益经济，但要对地球上全部物种多样性进行高质量基因组组装仍需漫长时日；即便组装完成，跨时间和空间的持续监测仍要求反复高覆盖度测序，成本高昂。此外，地球上大量生物多样性尚未有基因组代表，保护生物学家无法等待此类采样完成。单一标记基因的条形码方案在物种水平分辨率有限。

基因组浅测序（genome skimming）作为一种低覆盖度、免组装的替代方案被提出并取得一定成功。基于k-mer草图技术的方法，如Mash，可通过最小哈希估计两个样本k-mer集合间的Jaccard指数（J），并转换为核苷酸距离；Skmer进一步建模低覆盖度和测序误差的影响；ReSkmer则考虑了重复序列的影响。然而，现有方法均隐含地将基因组视为单倍体，忽略二倍体样本内杂合性对距离估计的干扰。对于微生物等单倍体生物，这一简化尚可接受；但对于生态学关注的大多数二倍体物种，忽略样本内杂合性会在群体遗传分析中引入实质性偏差。因此，研究人员旨在建立一种能够显式建模二倍体杂合性、适用于低覆盖度基因组浅测序数据的遗传多样性估计方法。

二、研究内容与主要结论

研究人员基于Kingman溯祖模型（Kingman's coalescent model），从理论上推导了两个二倍体个体（n=4条单倍体染色体）k-mer集合间Jaccard指数与群体规模突变参数θ之间的关系。通过枚举四种单倍体在两种可能溯祖拓扑下的突变分布模式，研究人员推导出二倍体条件下期望交集（I）与并集（U）的表达式，进而得到从Jaccard指数估计θ的解析公式。随后，研究人员引入低覆盖度和测序误差模型，估计k-mer覆盖度（C）和测序错误率（ε），并利用k-mer频次谱（k-mer frequency spectrum）推导出二倍体样本的参数估计方程。最终形成的DipSkmer方法以两套亲缘关系较近的二倍体生物读段为输入，输出有效群体规模突变率θ。

在模拟数据中，通过msprime模拟不同有效群体大小的二倍体群体，并利用ART模拟不同覆盖度（0.5×至8×）的读段，研究人员比较了DipSkmer与Mash、Skmer、ReSkmer的估计精度。结果表明，DipSkmer在大多数覆盖度和遗传距离条件下均显著优于其他免比对方法，整体平均绝对误差（mean absolute error，MAE）为25%，远优于ReSkmer（67%）、Skmer（86%）和Mash（619%）。在低覆盖度（0.5×和1×）条件下，DipSkmer的优势尤为突出，而Skmer在θ较低时经常输出零值。在群体分化模拟中，DipSkmer对Weir和Cockerham F_ST和D_xy的估计在多数条件下最为准确，尤其在低覆盖度下优势明显；但在接近物种形成级别的高分化水平（1个溯祖单位）下，其准确性下降。

在实证数据中，研究人员使用包含6个物种、12个种群的公开数据集，将DipSkmer的θ估计与基于参考基因组比对的ANGSD结果进行比较。DipSkmer在所有覆盖度下均与ANGSD高度相关（平均R=0.92），优于Skmer（R=0.9）和ReSkmer-noref（R=0.9），远优于Mash（R=0.43）。在绝对数值上，DipSkmer与ANGSD的MAE为28%，优于ReSkmer-noref（41%）、Skmer（83%）和Mash（738%）。在Sillago sinica三种群的种群结构分析中，基于DipSkmer距离的主坐标分析（PCoA）投影与ANGSD结果高度一致（t₀=0.821），且其F_ST估计与ANGSD最为接近。此外，DipSkmer计算速度快，每个样本平均仅需5分钟，与Skmer相当。

三、主要技术方法

研究人员首先利用Kingman溯祖模型推导二倍体条件下k-mer集合Jaccard指数与θ的解析关系，枚举四种单倍体在平衡与非平衡两种拓扑下的突变分布模式。其次，建立低覆盖度和测序误差模型，通过k-mer频次谱估计覆盖度C和错误率ε。第三，采用迭代策略：先用样本内k-mer频次谱估计θ，再据此改进覆盖度估计，交替迭代三次以获得更优的θ估计。模拟数据使用msprime生成二倍体群体基因组，ART模拟读段；实证数据来自已发表的6物种12种群基因组浅测序数据，并亚采样至目标覆盖度。

四、研究结果

4.1 二倍体改变了Jaccard指数与遗传多样性的关系。 研究人员通过溯祖模型推导出二倍体条件下Jaccard指数与θ的解析关系，发现二倍体估计值高于单倍体估计值，且仅在J>1/3时有定义。这是因为在无限位点假设下，即使遗传距离较高，每个个体两个等位基因仍从同一群体等位基因池中采样，保持一定相似性。

4.2 模拟结果。 在模拟数据中，DipSkmer在大多数条件下MAE为25%，显著优于ReSkmer（67%）、Skmer（86%）和Mash（619%）。在0.5×和1×极低覆盖度下，DipSkmer远优于其他方法；在2×覆盖度下达到最佳精度（17% MAE）。在重复序列含量较高的基因组中，DipSkmer仍保持竞争力，但在极高重复含量且低θ条件下，ReSkmer-ref可能更优。在群体分化模拟中，DipSkmer对F_ST和D_xy的估计在低至中度分化水平下最为准确，但在接近物种形成级别的高分化下准确性下降。迭代策略可将θ估计的平均误差从46%降至29%。

4.3 生物学结果。 在实证数据中，DipSkmer的θ估计与ANGSD高度相关（平均R=0.92），且MAE为28%，优于其他免比对方法。在Sillago sinica种群结构分析中，DipSkmer的PCoA投影与ANGSD高度一致，F_ST估计与ANGSD最为接近。在所有覆盖度下，DipSkmer均能稳定提供与ANGSD高度一致的估计，而Mash在低覆盖度下几乎无相关性。

五、讨论与结论

DipSkmer通过建模二倍体杂合性对k-mer集合Jaccard指数的影响，实现了从低覆盖度基因组浅测序数据中准确估计群体规模突变参数θ。该方法在模拟和实证数据中均优于现有免比对、无参考方法，且与基于参考基因组的高覆盖度方法ANGSD高度一致。DipSkmer消除了对组装基因组的需求，为研究较少的物种提供了低成本群体监测的可能性。值得注意的是，DipSkmer在2×覆盖度下反而比4×更准确，研究人员建议在实际覆盖度为3–5×时将读段降采样至2×。尽管基于Kingman溯祖模型和随机交配假设，DipSkmer对中等程度的群体结构具有经验稳健性，但应限于种内或极近缘物种间的比较。未来研究可尝试将重复序列建模与杂合性建模结合，并扩展到多倍体场景。研究结论表明，DipSkmer通过遗传多样性估计实现了可及、低成本的群体监测，为生物多样性保护提供了有力工具。