CN117577184A - 一种可用于大规模基因组的多基因组比对方法 - Google Patents

Abstract

本方案公开了一种可用于大规模基因组的多基因组比对方法，属于计算机生物学技术领域，提出BWT‑FM‑LIS循环分治的方式检索全基因组间公共子串，尽可能多地检索出基因组间的公共子串和关键路径，缩短需动态规划对齐的差异子字符串长度，进而降低比对的复杂度以使本方案可用于大规模大基因组的比对，适用于大规模长序列或超长基因组数据，可以识别出基因组比对中多种结构变异。

Description

一种可用于大规模基因组的多基因组比对方法

技术领域

本方案属于计算机生物学技术领域，提出了一种可用于大规模基因组的多基因组比对方法。

背景技术

目前研究比较成熟的基因组序列比对工具主要为双序列比对，已被广泛地应用于基因组学研究。大多数工具采用的是“种子-扩展”方法，如BLAST，首先识别出序列间短且无插空的序列匹配(种子)，然后使用Smith-Waterman改进算法从种子序列两端进行扩展匹配，直到对齐的分数低于指定阈值。种子类型可按照是否要求精确、连续匹配或匹配长度是否固定进行分类，BLAT和STELLAR采用k-mers作为其精确匹配种子，LASTZ(BLASTZ的更新版本)使用间隔种子，LAST利用后缀数组可查找长度不等的自适应种子，MUMmer和CHAOS则分别使用后缀树和线程树数据结构快速查找所有具有一定最小长度的精确连续匹配种子。各类型种子中往往自适应、非精确或者间隔种子具有更高的灵敏度。扩展延伸之后，邻近且顺序、方向一致的局部对齐还可以被AXTCHAIN等链接程序链接起来，形成更大的对齐作为输出。

基因组的多序列比对工具开发相比于双序列比对较晚，巨大的时空复杂度也限制了多序列比对工具的应用。目前已开发的基因组多序列比对工具主要有分层和局部比对方法。分层比对代表工具有：M-GCAT、progressiveMauve、Mugsy、progressiveCactus、Parsnp和FAME；能完成渐进多序列全局比对的Parsnp、MAVID、MLAGAN、SeqAn::T-Coffee和PECAN；采用贪心多序列全局比对的DIALIGN；以及基于概率统计隐马尔可夫模型的FSA。

目前种间基因组多序列比对软件最强的为progressiveCactus(下简称Cactus)，可比对约600个羊膜动物基因组，但是其限制是需要科学的指导树作为先验知识输入，无法处理兼容windows中的字符。群体内基因组多序列比对软件最强的为Parsnp，可比对数千个微生物基因组，但是其缺点是具有木桶效应，即只能找到所有序列之间共有的区域，导致随着序列数量增多，其比对质量严重下降，而且其无法处理相似度较低或长度差异较大的数据。虽然Parsnp，Cactus已经分别在数万个线粒体基因组，数千个高相似性微生物基因组和约600个羊膜动物基因组上的成功应用。但是目前还尚不存在可以比对数千个哺乳动物基因组甚至更大规模的多基因组比对软件。

发明内容

本方案的目的是针对上述问题，提供一种可用于大规模基因组的多基因组比对方法及其系统，适用于大规模长序列或超长基因组数据，可以识别出基因组比对中多种结构变异。

一种可用于大规模基因组的多基因组比对方法，该方法包括：

S1.对待比对的基因组序列进行数据预处理并确定中心序列；

S2.为中心序列建立BWT索引得到BWT数据结构，并建立预测字典；

S3.使用预测字典预测每条非中心序列的方向；

S4.根据公共子串长度阈值为每条非中心序列在中心序列的BWT数据结构中进行正向的公共子串检索；这里的正向是相对中心序列而言，以中心序列的正向为正向，中心序列的反向为反向，例如一条非中心序列的方向是A-B，中心序列的方向为B-A，则对该非中心序列进行一次B-A的单向检索

S5.筛选最长公共子串组合为相应非中心序列的主链；

S6.过滤主链，得到差异子字符串对，判断差异子字符串的长度是否满足动态规划，若是，则继续S7，否则，剔除主链并降低公共子串长度阈值，重复S4-S6；非中心序列和中心序列主链对应的区间是完全一样的，是对齐的，无需继续比对，直接至最终结果即可，除了第一次循环，后续每次执行步骤S4，其检索的中心和非中心序列都是剔除了最新主链后的序列。

S7.使用动态规划方法对差异子字符串对进行精细比对；

S8.输出多序列比对结果。

本方案提出BWT-FM-LIS循环分治的方式检索全基因组间公共子串，尽可能多地检索出基因组间的公共子串和关键路径，缩短需动态规划对齐的差异子字符串长度，进而降低比对的复杂度以使本方案可用于大规模大基因组的比对。

在上述的可用于大规模基因组的多基因组比对方法中，步骤S1中，

数据预处理包括：采用二元组的形式记录简并碱基的位置，并将其与gap一起删除；

确定中心序列的方式为，将去掉简并碱基后的最长的序列确定为中心序列，或者由用户指定；

步骤S2中，通过对中心序列的正反链的Kmer分别建立两个哈希字典得到预测字典，预测字典包括中心序列的正向链Kmer字典和反向链Kmer字典；

步骤S2还包括，采用二进制压缩存储方式对中心序列进行存储；

步骤S3中，在中心序列的正向链Kmer字典和反向链Kmer字典中查找每个非中心序列与之共有的Kmer数量，并依此确定非中心序列的方向。

在上述的可用于大规模基因组的多基因组比对方法中，步骤S4中，正向的公共子串检索方式为：

S1.从非中心序列的起始位置B_index开始，初始len为0，不断对len加1在中心序列上查找全部的相同子串，直到没有公共子串，回退一步，得到当前步长len下的x个公共子串；

S2.判断len是否大于或等于公共子串长度阈值，若是，则记录该len下找到的公共子串信息，并使B_index＝B_index+len，回到步骤S1重复检索；

若否，则使B_index＝B_index+1，回到步骤S1重复检索；

S3.上述过程持续检索至非中心序列的末尾。

ACGTAAACGT

ACGTGACGA

以上序列对为例，上方的为中心序列，下方的为非中心序列，该示例中，从非中心序列的A开始，初始步长len为0，有4个公共子串，然后加1为1，在中心序列中找AC，有2个公共子串，然后再加1为2，找ACG，有2个公共子串，再加1为3，找ACG，有2个公共子串，再加1为4，找ACGT，有2个公共子串，再加1为5，找ACGT，无法找到公共子串；

然后可以得到4时为两个公共子串

假设阈值为3，那么len＝4大于阈值3，则记录len＝4的两个公共子串信息。然后使起始位置＝B_index+4开始，len初始为0，继续上述步骤，即从非中心序列的G开始，初始步长为0，无法找到公共子串。

ACGT

CCGT

对于以上示例，C走一步以后为CC，中心序列没有公共子串，假设阈值大于1，那么此时len小于阈值，从第二个C开始，len＝0时，找到一个公共子串，len＝1时，找到一个公共子串CG，len＝2时，找到一个公共子串CGT，len＝3时，没有公共子串，所以len＝2，有一个公共子串CGT。

在上述的可用于大规模基因组的多基因组比对方法中，被记录的所述公共子串信息根据检索结果为一对一模式或一对多模式；

步骤S5中寻找主链的方式包括：

S51.首先将针对相应非中心序列的检索出来的全部公共子串按照相对位置最近原则筛选为一对一模式；

S52.查找重叠公共子串和交叉公共子串，针对每组重叠公共子串和每组交叉公共子串，保留非中心序列中对应的最长的同源区间，舍弃非中心序列中对应的其余同源区间，以筛选得到最长的同源区间组合，最终得到最长的不重叠、不交叉的公共子串组合就是主链。

公共子串有些是一对一的，有些是一对多的，这里统一筛选为一对一模式，即非中心序列上的一段区间对应中心序列上的一段区间，主链是由公共子串组合起来的，公共子串是指两条序列上完全一样的小片段，所以每条非中心序列的主链，中心序列都有对应的一条，这里按照相对位置原则筛选一对一模式，一方面能够从中心序列中选出对应的主链，另一方面能够尽可能地减少两序列的偏移，使后续比对更整齐。此外，本方案还利用寻找主链剩下的公共子串寻找结构变异链，利用全基因组的公共子串来识别规则的结构变异链，然后填充链条孔洞来得到结构变异区域比对，实现了结构变异的科学提取方法。

在上述的可用于大规模基因组的多基因组比对方法中，步骤S7中，采用Kband动态规划对差异子字符串对进行精细比对：

先初始化状态表；

然后利用仿射罚分计算出当前行得分，依据当前行和前一行得分的关系，填充路径转移方向到状态表；

最后由状态表回溯出最优路径。

在上述的可用于大规模基因组的多基因组比对方法中，Kband动态规划中，

采用仿射罚分对首次和继发插空进行区分，使首次插空的扣分大于继发插空；

通过限制动态规划范围节省时空成本：将动态规划限制在对角线附近，只对对角线附近的区域进行计算；条带的基线由两序列长度差确定，以所述条带为基准，右上扩，左下扩，从1开始以2的幂次变大，直到变化后一次的得分不比前一次的得分高；

通过路径存储简化数据结构：在计算的过程中单独设置一个存储路径的信息表，回溯时直接根据路径信息表反向查找最佳路径；计算过程中，存储当前三行状态，然后计算一行覆盖一行，全程记录转移到各状态的路径转移方向；

开辟全局空间替代反复申请和释放：一次性申请可能会利用到的最大内存，将重复申请释放内存改为重复读写。

在上述的可用于大规模基因组的多基因组比对方法中，步骤S7中，通过星比对的方式对差异子字符串对进行精细比对，得到每条序列的插空信息；

且本方法步骤S7中，对精细比对后的结果进行数据后处理：

将预处理去掉的简并碱基整合进最后的比对；

简并碱基插入过程为，首先计算插入后的长度，然后开辟相应大小的向量(vector)，并按顺序将原始串信息和插入信息写入该向量中，从而实现插入操作。

在上述的可用于大规模基因组的多基因组比对方法中，当步骤S8被要求输出结构变异比对结果时(用户可根据需要选择输出多序列比对结果(FSATA模式)或同时得到两个结果(MAF模式))，步骤S4中对每条非中心序列进行双向的公共子串检索，以找到可能发生的倒位结构变异；

步骤S5将得到由正向的最长公共子串组合构成的正向主链和由反向的最长公共子串组合构成的反向主链，并将正向主链和反向主链合并得到相应非中心序列的最终主链，合并过程包括：

找到正向主链和反向主链；将反向主链映射至正向；合并正向主链和映射后的反向主链。

在上述的可用于大规模基因组的多基因组比对方法中，双向检索中，正向检索方式为：

从非中心序列的起始位置起，经过步长len，若在中心序列上找到了x个公共子串，则判断len的长度是否大于公共子串长度阈值；

若是，则将检索结果记录为的公共子串信息，然后从步长+len的位置继续重复上述检索；

若否，则从步长+1的位置继续重复上述检索；

上述过程持续检索至非中心序列的末尾；

反向检索方式为：

从非中心序列的起始位置起，经过步长len，若在中心序列上找到了x个公共子串，则判断len的长度是否大于公共子串长度阈值；

若是，则将检索结果记录为的公共子串信息，然后从步长+len的位置继续重复上述检索；

若否，则从步长+100的位置继续重复上述检索；

上述过程持续检索至非中心序列的末尾；

在上述的可用于大规模基因组的多基因组比对方法中，本方法中，本方法还包括：

将寻找正向主链和反向主链过滤掉的公共子串作为寻找其他结构变异链的数据来源；

在所述数据来源分别进行两次寻找工作，分别找到正向的结构变异链和与倒位组合的结构变异链；

针对覆盖重叠的结构变异链进行选择保留处理；

步骤S7中，使用动态规划方法对结构变异链的公共子串之间区域进行比对填充，得到整个对齐的结构变异比对，并写出至maf文件；

多序列比对结果被写出至fasta文件；

maf文件用于存储结构变异比对结果，fasta文件用于存储多序列比对结果；

步骤S8中，将maf文件和fasta文件合并输出最终的多基因组比对结果：

先将maf文件中的block依据中心序列index排序并且筛选block中的比对质量超过原始fasta文件中比对的序列；

然后依据block中每条序列的起点终点作为切割位点，将fasta比对矩阵进行列切割；

最后依据星比对将maf文件的block和fasta文件切割后对应的block进行整合，同时删除结构变异在原始fasta文件中的比对；

合并之前先对maf文件进行排序，并将排序结果保存，然后，逐个Block读入并进行合并操作，并将合并后的结果写出；且合并过程采用逐层合并的方式以使每个层级内部的合并均相互独立。

上述方法使本方案提供的算法能够实现结构变异maf文件的输出，并且在最后针对结构变异maf文件和多序列比对结果fasta文件进行合并，输出最终的整合了结构变异的最终MAF格式的多基因组比对结果。

本方案的优点在于：

(1)本方案创新了一种循环分治法检索全基因组间公共子串的方法，通过组合了线性时空复杂度的BWT、数据压缩、FM-index和最长上升子序列(LIS)等算法的BWT-FM-LIS循环分治法，尽可能多地检索出基因组间的公共子串和关键路径，缩短需动态规划对齐的差异子字符串长度，进而降低比对的复杂度以使本方案可用于大规模大基因组的比对；

(2)本方案在大规模基因组的多序列比对中引入仿射罚分、路径存储、Kband限制等方法，用于针对优化高相似性基因组的动态规划，在提高比对结果生物学意义的同时，降低动态规划过程中的计算时间和存储空间，在大规模人全基因组的多序列比对中，时空复杂度占比最高的是对差异子字符串的动态规划，该改进可大大降低比对的计算时间和存储空间；

(3)本方案利用全基因组的公共子串来识别规则的结构变异链，然后填充链条孔洞来得到结构变异区域比对，实现了结构变异的科学提取方法；

(4)本方案提出了一种结构变异与MSA结果整合的方法，采用比multiz更高效的创新型星比对策略来整合结构变异到多序列比对结果；

(5)基于本方案实现的软件是首个可以完成数千人类一号染色体多基因组比对的软件，大大提升了多基因组比对可处理的数据规模的上限，可以完成目前所有软件都无法完成的大规模基因组比对任务。而且其他软件可以完成的小规模基因组比对任务上，本方案也在表现出足够优秀质量的同时，大大降低了时空成本的消耗。

附图说明

图1为本发明实施的整体流程图；

图2为本发明流程中的具体算法流程汇总，其中，

A表示将一条DNA序列压缩成原来空间1/4的示例；

B表示在公共子串中寻找主链的示例过程；

C表示合并正逆主链的流程；

D是寻找严格结构变异链的示例；

E是仿射罚分的示意图；

F是Kband动态规划的一个具体示例；

G表示回溯路径表的压缩存储；

H表示一个3条序列删N插N的具体示例；

I表示结构变异的maf结果文件整合到多序列比对的FASTA结果文件的具体流程；

图3所示为Cactus、Parsnp和Halign-G在18个人的24组染色体数据集上的比较分析结果，其中，

A是三个软件M-score的分布情况；

B是三个软件在实验中的时间和内存消耗情况；

图4所示Cactus、Parsnp、Mugsy、Kalign、MLAGAN、PMauve和Halign-3共七个软件与Halign-G的FSATA和MAF两个模式在五个数据集上的对比分析，其中，

A是九个软件M-score的分布情况；

B是九个软件在实验中的时间和内存消耗情况；

图5所示为Cactus、Parsnp、Mugsy、Pmauve与Halign-G在9组模拟数据集识别到的结构变异情况对比。

具体实施方式

下面结合附图和具体实施方式对本方案做进一步详细的说明。

本实施例给出了一种可用于大规模基因组的多基因组比对方法，主要在于提出了循环分治法检索全基因组间公共子串的方法，针对优化高相似性基因组的动态规划方法、结构变异的科学提取方法以及结构变异与MSA结果整合的方法。本实施例以人类1号染色体基因组为例对本方案进行说明。

总体流程图见图1，以下为具体流程步骤：

S1.数据预处理与确定中心序列

本方案采用FASTA文件格式作为输入数据的格式，输入数据中可包含各种简并碱基和gap字符，输入形式可以为直接输入一个包含多条序列的FASTA文件，也可以为指定一个文件夹，本方案软件将遍历该文件夹中的所有文件进行处理。本方案采用二元组<index，number>的形式记录简并碱基的位置，并将其与gap一起删除。根据序列特征确定当前序列是RNA还是DNA序列。此外，本方案将剩余的所有核酸字符转换为大写形式写入到中间文件中，以便后续需要时在外存调用而不是常驻内存，以方便处理超大规模的数据。去掉简并碱基后的最长的那条序列所包含的有效信息最多，故选择该条序列作为中心序列。当然，实际应用中，中心序列也可由用户指定。

图1所示，原始序列中，最长的DNA序列为被确定的中心序列，在其他序列中，上面几条为正向的非中心序列，最下面一条为逆向的非中心序列，黑色五角星代表简并碱基，虚线框中的深色箭头环绕的流程就是BWT-FM-LIS循环分治。

S2.为中心序列建立BWT索引和预测字典

为中心序列建立BWT索引

BWT索引建立过程采用当前最好的后缀数组构造器—可并行的LibDivSufSort，它可以在线性时间复杂度内构建出后缀数组。构造器是一个功能函数，BWT索引也可以看作是后缀数组，建立BWT索引就是通过构造器的功能函数求出这个后缀数组。通过为中心序列建立BWT索引得到中心序列的BWT数据结构。

为中心序列建立预测字典

在多基因组比对中通常会处理来自不同物种或个体的基因组序列，这些序列可能存在不同的方向。本方案使用MurmurHash3技术对中心序列的正反链的Kmer分别建立两个哈希字典—正向链Kmer字典和反向链Kmer字典，后续将使用哈希码的键在两个字典中进行快速的搜索。

中心序列压缩存贮

为了在比对过程中节约内存，本方案采用了二进制压缩存储方式对中心序列进行存储。这种压缩存储方法将每个字节(byte)用来存储4个碱基，其中一个碱基占用2位。通过将多个碱基紧凑地存储在一个字节中(见图2中A)，可以显著减少内存的使用量(降低到原来的四分之一)，提升可比对数据的规模。

S3.检索公共子串

kmer预测方向

在中心序列的正向链和反向链的两个Kmer字典中查找每个非中心序列与之共有的Kmer数量。通过比较这两个数值的大小确定该非中心序列是正向链还是反向链(相对于中心序列而言，与中心序列同向为正向链，反向则为反向链，一条非中心序列若与正向链Kmer字典的Kmer数量多于与反向链Kmer字典的公共Kmer数量，则确定为正向链，反之则确定为反向链)。

在进行超长序列预测时，因为长序列包含的Kmer数量极多，即使只对该序列的前半部分Kmer在两个字典中进行检索预测错误的概率也极低，故只对前半部分进行检索以加速预测过程。而对于短序列，由于短序列包含的Kmer数量较少，有效信息较少，偶然性较大，预测错误的概率相对于长序列来说更高，故遍历整个序列来统计与中心序列两个方向上相同Kmer的数量。同时，因为短序列通常是基因级别的序列数据，其方向基本上都是与中心序列的方向相同，所以本方案在对同向计数上加一定的权重，更偏向于中心序列相同。

单向检索公共子串

对于每条非中心序列，本方案按顺序使用FM-index在中心序列的BWT数据结构上进行检索，以找到所有与该非中心序列匹配的公共子串。初始时，由于检索距离较近，因此能够找到许多公共子串。随着步长逐渐增加，直到没有公共子串被找到为止，然后回退一步停止。在这个过程中，可以确定该序列从起始位置B_index开始，经过步长len后，在中心序列上找到了x(x>0)个公共子串。同时，还可以获得这些公共子串在中心序列上的起始位置(即对应的后缀数组SA的值)A_index_1、A_index_2…A_index_x。接下来，判断len的长度是否大于公共子串长度阈值threshold_BWT。如果是，那么就将这些信息记录为公共子串信息，记录的公共子串信息根据检索结果为一对一模式或一对多模式，一对多的模式形如[B_index,A_index_1,len]、[B_index,A_index_2,len]…[B_index,A_index_x]的形式。然后，继续下一步的检索，从B_index+len的位置开始。如果len长度不满足threshold_BWT，那么就重新从B_index+1的位置开始检索。这个过程持续到该非中心序列检索到末尾为止，则该序列在该方向上的公共子串检索完成。threshold_BWT这个阈值在每次循环迭代中提前设定，并且一次比一次小，而这个阈值越大x就会越少，阈值越小x就会越多。

双向检索公共子串

对于fasta格式输出，只进行上方所描述的一次正向公共子串检索。而对于maf(multiple alignment format)格式输出，本方案则在正反两个方向上都进行检索，以找到可能发生的倒位结构变异。与正向检索不同的是，在反向检索失败时，不是从B_index+1开始下一步检索，而从B_index+100开始下一步检索。倒位结构变异在整个序列中并不常见，因此在逆向检索中找到达到threshold_BWT的公共子串的概率较低，导致检索失败的可能性较高，如果单步移动，检索速度非常慢。每次失败后增加100个位置进行跳跃移动，能够加快检索速度，快速跳过当前概率极低的区域，提高检索效率。

S4.筛选最长公共子串组合及寻找SV链

查找出的公共子串的覆盖率会非常高，位置偏移小且相互之间有重叠。而且其中包含着由于人种差异以及癌症等疾病引起大片段的复制/丢失、重排和倒位的结构变异，因此寻找主链和结构变异链是重中之重。

寻找主链

针对小规模基因组数据，因为在上方寻找到的公共子串数量不多，在使用本方案方法时，可以采取最精准但是最耗时的动态规划算法来寻找主链。针对大规模基因组数据，公共子串数量庞大，时间复杂度O(n2)的动态规划算法运行极其缓慢。所以本方案采取以下步骤来对其进行筛选：

1.首先将检索出来的全部一对多的同源区间，按照相对位置最近原则筛选成一对一的模式，以尽可能地减少两序列的偏移，使后续比对更整齐；

2.接下来解决重叠问题，将公共子串同源区间看作一个个区间，按区间右端排序后，利用贪心算法筛选出最长非重叠子区间问题。

3.最后利用最长上升子序列LIS算法解决交叉排列问题，以筛选得到最终同源区间组合。

如图2中B给出了在公共子串中寻找主链的示例过程，图中给出了样本序列(即非中心序列)相对参考序列(中心序列)检索出X1-X5五个同源区间，X1对应Y1、Y2，X2对应Y3，X3对应Y8，X4对应Y4、Y5、Y6，X5对应Y7，其中X1和X4为一对多的同源区间，故先通过相对位置筛选X1和X4在多个同源区域中的一个，使其成为一对一的模式，X1选Y1，X4选Y5。然后解决重叠问题，发现X4(Y5)和X5(Y7)存在重叠，而保留X4整体具有区间长度最长，所以舍弃X5(Y7)代表的同源区域。最后解决交叉排列问题，发现X3(Y8)和X4(Y5)存在交叉，而保留X4整体具有区间长度最长，所以舍弃X3(Y8)代表的同源区域，至此，得到了X1(Y1)、X2(Y3)、X4(Y5)组合而成的主链。

在本方案的MAF模式下，寻找公共子串需要在正向和反向两个方向上分别进行。在得到两个方向上的公共子串后，寻找主链的工作也需要在两个方向上分别进行，从而得到正向和反向上的最长公共子串组合。其中，与中心序列同向的最长公共子串组合最长且包含最全的信息，而反向最长公共子串组合实质上是由少数倒位结构变异所形成的。反向最长公共子串组合会被分成多个倒位簇，每个倒位簇对应着一个倒位结构变异。将这些倒位簇与正向主链进行合并。此外，在发生重叠时，本方案裁剪相应的倒位簇，最终得到真正的全局主链。一旦倒位簇整合到全局主链后，倒位簇对应的正向非同源区间则不需要进行动态规划或者分治比对，从而减少了比对低质量区域的开销。

如图2中C表示了合并正逆主链的流程。第一步先找到了实线箭头所表示的正向主链，第二步找到了虚线箭头所表示的反向主链，第三步将反向主链映射到正向，第四步合并正逆主链，图中实线箭头和虚线箭头组成了最终的主链，其上的区间都不必再进行比对，其间隙才需进一步动态规划比对。

寻找SV链

除了反向主链中的倒位簇是一种SV链以外，本方案将为寻找正逆主链而过滤掉的公共子串作为寻找其他SV链的数据来源。对这部分数据分别进行两次寻找SV链的工作，分别找到正向的结构变异和与倒位组合的结构变异。优选地，本方案将用几个阈值来限制疑似链中的公共子串，以找到这些由紧凑平行的公共子串构成的SV链。如图2中D，其具体限制规则是：Alen>thresh_SV；Blen>thresh_SV；disA<thresh_dis；disB<thresh_dis；all_len/min(ALen,BLen)>thresh_RATIO；abs(disA-disB)/min(disA,disB)<thresh_sub_ratio；abs(ALen-BLen)/min(ALen,BLen)<thresh_SUB_RATIO，thresh_SV是由用户指定的结构变异的最短长度，thresh_dis是相邻公共子串的最远距离，thresh_RATIO是公共子串总长占链长度的最小比值，thresh_SUB_RATIO是两条序列在该链中总长度的最大相对误差，thresh_sub_ratio是链中相邻两个公共子串之间距离的最大相对误差。

在得到一些紧凑平行规整的结构变异链后，若存在覆盖重叠，则进一步对这些链进行筛选。就双基因组比对而言，在比对结果文件中，针对某个局部对齐结果(block)是否保留，判断该block中的ref区域和sample区域是否均被其他block覆盖过，如果是则不保留，反之则保留。最后利用Kband动态规划对链中公共子串之间区域的比对填充，得到整个对齐的结构变异比对，以block的形式写出到maf文件中。

S5.非同源区间的BWT递归分治

BWT/FM-index可以快速将中心序列进行转换压缩和检索公共子串，而通过贪心和求解LIS问题可筛选出匹配中心序列最优的公共子串组合。因此以上方法就组成了一个完整的求解双序列间最优公共子串组合的链条：BWT-FM-LIS。

针对像人类这样的大型基因组序列，本方案进行BWT-FM-LIS循环分治，通过在检索环节逐步降低公共子串长度阈值(默认全局threshold_BWT＝15，局部循环分治中threshold_bwt＝5)，来提高BWT-FM-LIS的灵活性，使其可应用于不同长度/相似度的双序列间求解最优公共子串组合。对于长度超长的基因组序列(如人类1号染色体)，可以采用较长的阈值，快速将248Mb的序列拆分，过滤掉主链后，得到较短但相似性稍低的差异子字符串对。然后再采用较短的公共子串长度阈值进一步分解差异子字符串对。如此循环迭代，直至差异子字符串对被分割成长度适于动态规划的片段(默认threshold_DC＝10000)。最后使用Kband动态规划算法进行精细比对。

S6.非同源区间的Kband动态规划

经BWT-FM-LIS循环分治后，基因组中非同源区间的比对就化简成许多短序列间的动态规划。本方案从打分细则、限制范围、数据存储和底层系统内存特点四个方面来优化动态规划。

仿射罚分提升比对质量

在基因组的各种变异中，发生频率通常为替换>>缺失>插入>>大片段的结构变异，其中首次发生插入或缺失后，在此基础上继续发生插入或缺失就相对容易。为比对结果更贴近自然规律，采用仿射罚分对首次和继发插空进行区分：首次插空的扣分大于继发插空，使得比对结果中位置靠近的插空能够尽可能地接邻，进而帮助后续得到更准确的变异及生物学意义，如图2中E，-σ代表首次gap罚分，-e代表后续(非首次gap)罚分，+s就是正常的比对计分由评分矩阵决定。一般都是首次gap罚分重一些，后续gap罚分轻(gap罚分指的A/C/G/T和-进行比对，-代表gap)。比对积分就只由评分矩阵决定4[ACGT]*4[ACGT]。

限制动态规划范围节省时空成本

针对高相似性序列通过限制动态规划的范围来将时空复杂度由O(n2)降为线性的O(kn)，其中n代表序列的长度。本方案使用Kband方法将动态规划限制在对角线附近，其主要利用的是两条高相似性序列得到最佳比对的回溯路径一般都在对角线附近。那么就不需要对整个矩阵进行填充和计算，只对对角线附近的区域进行计算即可，这个区域便称为Kband。条带的基线由两序列长度差确定，是靠近对角线的两条平行线之间的区域，然后以此构成的条带为基准，右上扩，左下扩，从1开始以2的幂次变大，直到变化后一次的得分不比前一次的得分高，即可结束当前的动态规划。如图2中F是Kband动态规划的一个具体实例，两条序列长度差diff为1，kband的宽度为diff+2*k+1，所以中间浅灰色条带(位于两条深灰色之间)k＝0，浅灰色+深灰色条带区域k＝1，k从0到1其最大分数和回溯路径没有发生变化，所以k＝1就停止扩张。

路径存储简化数据结构

动态规划的流程是先计算状态表，再根据状态表回溯出最优路径，然而加入仿射罚分就需要3个状态表去存储动态规划的打分。本方案在计算的过程中单独设置出一个存储路径的信息表，回溯时可以直接根据这个路径信息表反向查找最佳路径。据此计算过程中就不用存储整个状态表的打分，而只需要存储当前三行状态，然后计算一行覆盖一行，全程记录转移到各状态的路径方向即可。之前三个仿射罚分状态表用int存储，每个状态的int占4个字节，3个状态就是12个字节。而路径信息表每三个状态只需要1个字节，用其中的六个位就可以记录得到当前状态的路径转移方向(见图2中G)，这意味着此处内存需求将减少十倍左右。

开辟全局空间替代反复申请和释放

Kband算法已经处于程序较低层了，会被频繁调用；而在每个Kband方法内动态规划矩阵还要随着k的增大，频繁得申请释放。通过调研底层操作系统内存管理，发现对内存的频繁申请释放会使内存操作效率降低，致使内存虽然delete/free了，但是却没有释放掉，一直堆积，使内存虚假的增高。所以本方案一次性申请可能会利用到的最大内存，将重复申请释放内存改为重复读写，放弃Kband的kn空间复杂度选择了n2，使用Kband的kn时间复杂度来优化时间，达到使内存峰值降低到至理想状态的目的。

S7.数据后处理

星比对整合

在人类基因组大规模多序列比对中，以人的参考基因组为中心序列，正好与星比对策略的思想相契合。每条序列只需要跟参考基因组比对一次，然后按照“一旦为空，始终为空”的原则综合每组双序列比对的插空信息就可完成多序列比对。这是因为实际记录的只是插空信息(5号索引，插3个gap)就记录成<5，3>，因为插空不多所以这种记录方式不耗内存，而不是整条序列ACGT---ACGT，这样就很耗内存。该过程就是把原本双序列对齐的插空信息更新成在多序列比对可以使序列对齐的一些<>，最后跟N的存储方式一样，所以才有下面这个模块，将N和gap一起插回真实序列。

简并碱基N的整合

本方案通过星比对已经获得了每条序列的插空信息，但是当前的对齐是基于数据预处理去除简并碱基以后的(下面用N代指所有简并碱基)序列的。所以还要将预处理去掉的N再整合进最后的比对。有三类index对应于下面三类字符串：

原串：带N的原始串；

输入串：去掉N，仅有ACGT4种字符的串；

Align串：将“输入串”加上一些gap对齐后的串；(通过星比对整合的对齐)

1.程序主体不包含N，前期读入数据时去掉N，但是记录成(index,num)这样的二元组向量，后期得到gap的信息也是(index,num)这样的二元组向量。N和gap的index都是相对于“输入串”的。

2.整合两个插入向量为(index,n_num,g_num)这样的三元组向量，依据“输入串”索引倒序插入(因为先在后方插入不会影响到前方的索引，而先在前方进行插入，后续的索引会随之变化，插入就会紊乱出错)，但是保证以g_num为上限，即如果同一位点4个gap，6个N，只能插入4个N，剩下两个N后续考虑，因为要保证“Align串”的对齐不变。

3.对于多出的N和由于某条串多出N而使其它串多出的gap(同星比对原理)，再依据“Align串”的索引倒序插入。

4.最后处理特殊情况，假如3条序列比对，比对后某个位置本没有gap，但是有两条串有N，所以必然导致除自身外其余串该位点要插入gap，所以会导致多出一列全是gap，后来处理了一下对于“Align串”的插入向量，避免了这里的多余插入。如图2中H，表示了一个3条序列的具体删N，插N的具体示例，插回第1条序列的N导致了其余序列两个浅灰色的gap(“-”)，插回第2条序列的N导致了其余序列两个深色的gap，最后导致了有一列完全是gap构成的，这是不合理的，需要删除。

插入模块

在本方案中，对原序列的插入操作并没有直接进行实际的插入，而是将插入操作转换成了对插N和插gap的二元组向量的操作。这样的设计可将费时的实际插入操作延迟到最后进行，避免实际插入后因为某些原因未成功比对而导致时间的浪费。

虽然链表是对于插入操作最高效的数据结构，但是对于如此大规模的序列和插入信息来说，一般很难向操作系统申请到如此大的内存空间。同时，链表的空间利用效率也很低。在32位编译模式下，链表的每个节点由一个占用1字节的unsignedchar的数据块和一个4字节大小的指针构成；而在64位编译模式下，指针的大小为8字节。因此，链表的内存利用率仅为20％或11％。此外，链表还需要频繁进行内存申请和释放操作，导致底层操作系统对内存操作的效率降低。

为了解决这个问题，本方案采用了不同的插入思路：首先计算插入后的长度，然后开辟相应大小的向量(vector)，并按顺序将原始串信息和插入信息写入该向量中，从而实现插入操作。这样，每条序列的插入操作时间复杂度仅为线性遍历的O(n)复杂度。

S8.结构变异与MSA的合并

FSATA是本方案的输入格式，也是其中一种输出格式。MAF格式是本方案的另一种输出格式，与fasta格式相比，可以存储倒位，易位，重复等结构变异的对齐。本方案基于multiz合并两个maf文件的功能对其进行优化，使其能够被用于maf与fasta的合并。

multiz本身并没有排序功能，且每次都会将整个文件读入内存进行链表合并，对此，本方案在合并之前先对maf文件进行排序，并将排序结果保存，然后，逐个Block读入并进行multiz的合并操作，并将合并后的结果写出，避免了将整个maf文件常驻内存中，更适合处理大规模基因组数据。

此外，multiz只能以单线程方式处理两个文件的合并，对此，本方案采用了逐层合并的方法，每个层级内部的合并是相互独立的，因此可以采用多线程并行处理，适应本方案需要处理多个MAF文件进行多序列/基因组比对。

针对最后一步合并结构变异maf文件和多序列比对结果fasta文件，本方案开发了效率远超multiz的无丢失信息的整合模块(如图2中I)：首先将sv的MAF中的block依据中心序列index排序并且筛选block中的比对质量超过原始fasta中比对的序列，然后依据block中每条序列的起点终点作为切割位点，将fasta比对矩阵进行列切割，最后依据星比对将结构变异sv的block和fasta切割后对应的block进行整合，同时删除结构变异在原始fasta中的比对，这样就得到整合了结构变异的最终maf格式的多基因组比对结果。

为了体现本方案相对现有技术方案的优势，对本方案与现有技术的一些方案进行了对比分析，将基于本方案的软件称为Halign-G。

图3所示为Cactus、Parsnp和Halign-G在18个人的24组染色体数据集上的比较分析。A表示了三个软件M-score的分布情况，上方子图为0.9-1范围的扩大展示，每个竖线框中的三个箱线依次是Cactus、Parsnp、Halign-G的分布情况，为了更清晰的表示，分别用a表示Halign-G的分布情况，b表示Parsnp的分布情况，c表示Cactus的分布情况。可以看到，Parsnp和Cactus在13，14，15，22数据集上没有结果。Cactus的M-score在5，7，8，16，17五组数据集上超过了Halign-G，剩余15组则是Halign-G最高。B是三个软件在实验中的时间和内存消耗情况，气泡的高度表示时间长度，气泡的size表示内存的size，Halign-G颜色最深，其次是Cactus，最浅的是Parsnp，Halign-G均位于每一排的最左边，耗时最短，且圈最小内存消耗最少。可以看到，在时空成本上Halign-G的消耗最低。

图4所示Cactus、Parsnp、Mugsy、Kalign、MLAGAN、PMauve和Halign-3共七个软件与Halign-G的FSATA和MAF两个模式在五个数据集上的对比分析。A表示了9个软件M-score的分布情况。B展示了9个软件在实验中的时间和内存消耗情况，气泡的高度表示时间长度，气泡的size表示内存的size。对于Mitochondrial-genome、Streptococcus-pneumoniae和Escherichia-coil数据集，Halign-G的MAF模式在平均M-score上表现最好。值得注意的是，比对质量相对较低的FASTA模式在多序列比对层面也超过了Halign-3的M-score。此外，Halign-G内存是全部五组数据中所有软件中最低。图中分布情况和气泡本以颜色区分各软件结果，为了避免灰度后影响理解，现用x表示Halign-G的MAF结果，用y表示Halign-G的FSATA结果，其余为对照组的结果，对理解影响不大，故不单独标记。

图5所示为Cactus、Parsnp、Mugsy、Pmauve与Halign-G在9组模拟数据集识别到的结构变异情况对比。图中每一个五边形代表一个结构变异的标准答案位点，根据五种软件分成5个小三角形，每个小三角形被点亮则表示对应软件识别到该位点，五边形中，Cactus对应于右上方的小三角形，Parsnp对应于右下方的小三角形，Halign-G对应于正下方的小三角形，Mugsy对应于左上方的小三角形，Pmauve对应于左下方的小三角形。

实验结果显示，在共计108个结构变异位点中，Pmauve未找到任何位点，Cactus找到了21个位点，Parsnp找到了35个位点，Mugsy找到了52个位点，而Halign-G找到了107个位点，仅在第七组模拟数据中有一个位点未找到。这充分证明了Halign-G在识别结构变异方面的强大能力。

本文中所描述的具体实施例仅仅是对本方案精神作举例说明。本方案所属技术领域的技术人员可以对所描述的具体实施例做各种各样的修改或补充或采用类似的方式替代，但并不会偏离本方案的精神或者超越所附权利要求书所定义的范围。

Claims (10)

1.一种可用于大规模基因组的多基因组比对方法，其特征在于，该方法包括：

S1.对待比对的基因组序列进行数据预处理并确定中心序列；

S2.为中心序列建立BWT索引得到BWT数据结构，并建立预测字典；

S3.使用预测字典预测每条非中心序列的方向；

S4.根据公共子串长度阈值为每条非中心序列在中心序列的BWT数据结构中进行正向的公共子串检索；

S5.筛选最长公共子串组合为相应非中心序列的主链；

S6.过滤主链，判断差异子字符串的长度是否满足动态规划要求，若是，则继续S7，否则，降低公共子串长度阈值，重复S4-S6以进一步分解差异子字符串；

S7.使用动态规划方法对差异子字符串对进行精细比对；

S8.输出多序列比对结果。

2.根据权利要求1所述的可用于大规模基因组的多基因组比对方法，其特征在于，步骤S1中，

数据预处理包括：采用二元组的形式记录简并碱基的位置，并将其与gap一起删除；

确定中心序列的方式为，将去掉简并碱基后的最长的序列确定为中心序列，或者由用户指定；

步骤S2中，通过对中心序列的正反链的Kmer分别建立两个哈希字典得到预测字典，预测字典包括中心序列的正向链Kmer字典和反向链Kmer字典；

步骤S2还包括，采用二进制压缩存储方式对中心序列进行存储；

步骤S3中，在中心序列的正向链Kmer字典和反向链Kmer字典中查找每个非中心序列与之共有的Kmer数量，并依此确定非中心序列的方向。

3.根据权利要求1所述的可用于大规模基因组的多基因组比对方法，其特征在于，步骤S4中，正向的公共子串检索方式为：

S1.从非中心序列的起始位置B_index开始，初始len为0，不断对len加1在中心序列上查找全部的相同子串，直到没有公共子串，回退一步，得到当前步长len下的x个公共子串；

S2.判断len是否大于或等于公共子串长度阈值，若是，则记录该len下找到的公共子串信息，并使B_index＝B_index+len，回到步骤S1重复检索；

若否，则使B_index＝B_index+1，回到步骤S1重复检索；

S3.上述过程持续检索至非中心序列的末尾。

4.根据权利要求3所述的可用于大规模基因组的多基因组比对方法，其特征在于，被记录的所述公共子串信息根据检索结果为一对一模式或一对多模式；

步骤S5中寻找主链的方式包括：

S51.首先将针对相应非中心序列的检索出来的全部公共子串按照相对位置最近原则筛选为一对一模式；

S52.查找重叠公共子串和交叉公共子串，针对每组重叠公共子串和每组交叉公共子串，保留非中心序列中对应的最长的同源区间，舍弃非中心序列中对应的其余同源区间，以筛选得到最长的同源区间组合。

5.根据权利要求1所述的可用于大规模基因组的多基因组比对方法，其特征在于，步骤S7中，采用Kband动态规划对差异子字符串对进行精细比对：

先初始化状态表；

然后利用仿射罚分计算出当前行得分，依据当前行和前一行得分的关系，填充路径转移方向到状态表；

最后由状态表回溯出最优路径。

6.根据权利要求5所述的可用于大规模基因组的多基因组比对方法，其特征在于，Kband动态规划中，

采用仿射罚分对首次和继发插空进行区分，使首次插空的扣分大于继发插空；

通过限制动态规划范围节省时空成本：将动态规划限制在对角线附近，只对对角线附近的区域进行计算；条带的基线由两序列长度差确定，以所述条带为基准，右上扩，左下扩，从1开始以2的幂次变大，直到变化后一次的得分不比前一次的得分高；

通过路径存储简化数据结构：在计算的过程中单独设置一个存储路径的信息表，回溯时直接根据路径信息表反向查找最佳路径；计算过程中，存储当前三行状态，然后计算一行覆盖一行，全程记录转移到各状态的路径转移方向；

开辟全局空间替代反复申请和释放：一次性申请可能会利用到的最大内存，将重复申请释放内存改为重复读写。

7.根据权利要求2所述的可用于大规模基因组的多基因组比对方法，其特征在于，步骤S7中，通过星比对的方式对差异子字符串对进行精细比对，得到每条序列的插空信息；

且本方法步骤S7中，对精细比对后的结果进行数据后处理：

将预处理去掉的简并碱基整合进最后的比对；

简并碱基插入过程为，首先计算插入后的长度，然后开辟相应大小的向量，并按顺序将原始串信息和插入信息写入该向量中，从而实现插入操作。

8.根据权利要求1-7任意一项所述的可用于大规模基因组的多基因组比对方法，其特征在于，当步骤S8被要求输出结构变异比对结果时，步骤S4中对每条非中心序列进行双向的公共子串检索，以找到可能发生的倒位结构变异；

步骤S5将得到由正向的最长公共子串组合构成的正向主链和由反向的最长公共子串组合构成的反向主链，并将正向主链和反向主链合并得到相应非中心序列的最终主链，合并过程包括：

找到正向主链和反向主链；将反向主链映射至正向；合并正向主链和映射后的反向主链。

9.根据权利要求8所述的可用于大规模基因组的多基因组比对方法，其特征在于，双向检索中，正向检索方式为：

S1.从非中心序列的起始位置B_index开始，初始len为0，不断对len加1在中心序列上查找全部的相同子串，直到没有公共子串，回退一步，得到当前步长len下的x个公共子串；

S2.判断len是否大于或等于公共子串长度阈值，若是，则记录该len下找到的公共子串信息，并使B_index＝B_index+len，回到步骤S1重复检索；

若否，则使B_index＝B_index+1，回到步骤S1重复检索；

S3.上述过程持续检索至非中心序列的末尾。

反向检索方式为：

S1.从非中心序列的起始位置B_index开始，初始len为0，不断对len加1在中心序列上查找全部的相同子串，直到没有公共子串，回退一步，得到当前步长len下的x个公共子串；

S2.判断len是否大于或等于公共子串长度阈值，若是，则记录该len下找到的公共子串信息，并使B_index＝B_index+len，回到步骤S1重复检索；

若否，则使B_index＝B_index+100，回到步骤S1重复检索；

S3.上述过程持续检索至非中心序列的末尾。

10.根据权利要求9所述的可用于大规模基因组的多基因组比对方法，其特征在于，本方法中，本方法还包括：

将寻找正向主链和反向主链过滤掉的公共子串作为寻找其他结构变异链的数据来源；

在所述数据来源分别进行两次寻找工作，分别找到正向的结构变异链和与倒位组合的结构变异链；

针对覆盖重叠的结构变异链进行选择保留处理；

步骤S7中，使用动态规划方法对结构变异链的公共子串之间区域进行比对填充，得到整个对齐的结构变异比对，并写出至maf文件；

多序列比对结果被写出至fasta文件；

步骤S8中，将maf文件和fasta文件合并输出最终的多基因组比对结果：

先将maf文件中的block依据中心序列index排序并且筛选block中的比对质量超过原始fasta文件中比对的序列；

然后依据block中每条序列的起点终点作为切割位点，将fasta比对矩阵进行列切割；

最后依据星比对将maf文件的block和fasta文件切割后对应的block进行整合，同时删除结构变异在原始fasta文件中的比对；

合并之前先对maf文件进行排序，并将排序结果保存，然后，逐个Block读入并进行合并操作，并将合并后的结果写出；且合并过程采用逐层合并的方式以使每个层级内部的合并均相互独立。