组会论文 | preprint 构建泛基因组图

题目：Building pangenome graphs

DOI：10.1101/2023.04.05.535718

期刊：preprint

即使是同一个species，不同个体之间的基因组序列都是存在差异的，存在各种variation，比如inversion, SV, indel, SNP等等，那么用单一的线性的参考基因组是不能表征整个物种的遗传多样性。所以考虑整合多个基因组来作为一个参考系统，也就是泛基因组。现在最新研究是基于图形泛基因组泛基因组构建和分析。在图形泛基因组中，经常用variation graph构建，节点表示任意长度的共享序列，边表示节点之间的关系。有些图会把路径表示出来，路径表示一个基因组或者是单倍型。

图源：https://docs.google.com/presentation/d/13kOe67brYH4XcsnnuFajdN0g4_vtZCGCk2Xj7DoDclA/edit#slide=id.g28a5cec7b4e_0_34

引言

limitations:

1. 现有的构建泛基因组图的方法由于其基于参考（minigraph, minigraph-cactus）和树引导方法（PanGraph, Progressive Cactus）而存在偏差，这可能导致遗传变异的不完整和不稳定的表示。

2. 尽管已经提出了无偏泛基因组图构建的方法，比如基于variation graph的seqwish和基于compacted de bruijn graph的TwoPaCo，但这些方法仅限于构建图的步骤，而经验表明，需要专门的泛基因组比对和完善的技术来获得高质量的泛基因组构建。

为了解决这些问题，作者提出PanGenome Graph Builder (PGGB)

1. reference-free pipeline：整合多个软件的pipline，这些软件基本是由作者自己或者参与开发的。
2. Unbiased：所有基因组都被同等对待，PGGB不对系统发育关系、直系亲属群或进化史进行假设，没有输入顺序或系统发育依赖性，任何基因组都可以作为下游分析的参考框架。
3. base-level representation：它的输出呈现了泛基因组的碱基水平表示，包括从单核苷酸多态性（SNPs）到结构变异（SV）的所有规模的变异。

算法流程

图A展示了PGGB算法的整个流程。主要算法流程（红色）从FASTA比对（FASTA to alignment）、图构建（graph induction）和图平滑（graph smooth），然后使用GFAFIX进行归一化，生成最终图形（橙色）。可选的下游输出（蓝色），PGGB集成了一些常用步骤，使用ODGI，它生成基本的图形统计信息，如大小、节点计数和碱基含量；ODGI创建1D和2D可视化，可以应用到千兆规模的图，也就是Gbp。使用vg可以进行变异识别。并且可选择在MultiQC输出质量评估（黄色）。

一、 使用wfmash进行all vs all alignment

使用all vs all alignment避免基于参考和顺序产生的bias，图A中的这个表就是简要的互相比对。具体来说，首先是对序列切成小块，使用MashMap算法的变体根据所给的平均核苷酸相似度ANI找到大规模的同源序列。MashMap本身可以实现序列局部的jaccard距离估计，然后将其转换为Mash距离，也就是序列之间的相似性。然后对这些同源序列使用wflign，作者根据WFA（wavefront alignment）进行改编，WFA本身就是一个比对算法，利用序列的同源加速比对过程，替代以往的动态规划算法。

二、 使用seqwish进行图构建

图A中的alignment graph比较直观地说明是如何根据比对结果生成的，算法本身非常复杂。

三、 使用smoothxg进行图平滑（graph smooth）

这一步是文章详细介绍的

为什么要进行graph smooth？

1. all vs all结果会产生难以处理的复杂looping motifs。由于成对比对不是相互标准化的，这导致在低复杂度序列中对小的变异（如indel）的不同表示，从而产生难以处理的复杂looping motifs。
2. graph induction结果存在局部under-alignment。作者通过SEQWISH删除short matches，部分缓解了复杂looping motifs这个问题。这降低了复杂性，但同时也创建了一个局部“under-aligned”的图，并且不表示所有局部成对关系。

在图的所有部分应用local realignment kernel，偏序排序（partial order alignment, POA）。这样做的规模约为1kbp，这比基因组中发现的大多数非线性结构变异模式的大小要小。这使得PGGB图可以将复杂结构变异位点表示为重复序列的单个拷贝的简单环。

具体来说，smoothxg的算法流程：

1. 预处理。通过Path-Guided Stochastic Gradient Descent (PG-SGD) 路径引导的随机梯度下降算法优化图一维节点顺序，可以最大程度匹配嵌入路径（也就是基因组的Path）中的核苷酸位置。梳理（grooming)步骤确保对于每个节点，节点方向遵循访问该节点的所有路径步骤的一致节点（consensus nodes）方向。至此，一维拓扑排序完成了整个排序步骤。chop步骤，图形被截断（chopped），使每个节点包含相对较少的核苷酸（SMOOTHXG默认值：100），从而保留节点拓扑结构和顺序。这将在下一节中描述的块创建过程中为准备图形提供更精确的切割点。
2. 创建块。smoothable blocks是通过按照先前计算的顺序在迭代所有节点创建的。如果一个节点的添加满足以下条件，不超过：1.块的总路径长度限制 2.块的最大边跳跃限制 3. 最大块长度。在可能存在的可变串联重复次数（VNTR）边界处断开块，以确保每个块内的路径范围不超过下一节中描述的POA步骤的最大序列输入大小。
3. 平滑每个块。对于每个块，填充将每个块向左和向右扩展。这改进了每个块边界处的局部对齐。POA应用于每个块。可选地，该步骤为每个块生成MAF格式的多序列比对。填充将被删除，块将被保存，以便稍后集成到完整的图中。
4. 把每个块放回图。smoothed blocks按照其初始输入顺序排列在一起，形成最终的泛基因组图。最后一步，图不会再被截断，保留图中可能的最大节点长度。

整个SMOOTHXG步骤需要进行多次迭代（默认为3次），以限制在块边界附近可能发生的边缘效应，逐步完善。

四、 应用GFAFIX来压缩冗余节点，并使用ODGI对生成的图进行排序

图源：https://github.com/marschall-lab/GFAffix

结果

一、交叉验证

因为缺乏包含所有变异的真实数据，所以对真实数据的质量评估变得非常困难。这里把两个工具的结果进行对比，粗略评估。一个是MUMMER4（NUCMER），广泛使用的双序列比对工具，进行all vs all比对得到潜在的SNPs，将所有的SNPs信息进行合并。另外PGGB构建得到图文件，通过vg deconstruct进行变异识别，并过滤得到潜在的SNPs。

这些结果显示，在所有测试环境中，交叉验证F得分>90%，表明该方法的性能与现有标准相当。然而，虽然MUMMER4仅提供与目标参考的成对比较，但PGGB在基因组之间产生了完全的all vs all比较，从而产生了全新的生物信息学分析模式。

第一列是泛基因组的名字以及单倍型数量。第二列是图的大小。第三列是压缩比，泛基因组的长度除以图的大小。第四列是运行时间。第五列是最大内存使用量。第六列SNV为NUCMER比对得到的所有SNV（single-nucleotide variant）的数量，作为基准与PGGB比较。第七列是得到的F1分数。

二、PGGB提供的实用通用性

1. 许多通常基于相对于单一参考基因组的变异（例如SNPs）极化的下游应用，可以使用PGGB和相关方法构建的变异图的空间中直接实现。
2. 在变异图中，节点是等位基因，而基因组可以简单地理解为等位基因计数的载体。
3. 基于该向量空间的方法使我们能够在没有参考偏差的情况下同时考虑下游分析中的所有类别的变异。

作为证明，基于共享碱基对的jaccard距离构建的系统发育树与基于人工构建的SNPs集合来构建的系统发育树相同，如图D。另外图B展示了T2T灵长类6号染色体的一维泛基因组对齐。与T2T-CHM13方向相同的为黑色，反向互补为红色。MHC、p臂、着丝粒和q臂的T2T-CHM3注释已作为路径注入图中。图C展示了二维结构，成环原因是亚端粒的结构变异。

总结

1. 一种新的、模块化的、概念上简单的方法，用于在泛基因组学和比较基因组学环境中理解许多完整基因组之间的序列关系。
2. 为基因组图构建技术提供了一个通用框架。
3. 通过简化变异图的构建，PGGB为多样化的下游群体和进化遗传方法打开了大门，这些方法可以同时考虑所有类别的序列变异。