基于泛基因组的序列实时比对方法

技术领域

本发明涉及一种序列比对方法，具体涉及一种面向边合成边测序类测序技术的基于泛基因组的序列实时精确比对方法，属于序列比对领域。

背景技术

DNA序列比对是将测序仪生成的数据在参考基因组上进行定位并比较差异的过程，是发现变异、进行进一步的生命科学以及精准医学分析之前的核心步骤。传统的DNA序列比对算法需要在测序仪生成全部数据，进行单碱基召回(base calling)生成fastq类型的文件之后进行，鉴于序列比对本身需要消耗大量时间，所以从测序开始到获取完整的序列比对结果就需要大量的时间。

高通量边合成边测序类测序技术(第二代测序技术，简写NGS)由于极大的降低了DNA等生命大分子测序的成本，因此，极大地促进了基因科学乃至生命科学的发展，在生命科学以及精准医学领域具备广阔的应用前景。但由于其具备廉价、高通量、测序速度快等特点，使得DNA序列比对步骤耗费大量的计算资源以及计算时间。

发明内容

本发明为了解决由于DNA序列比对是在序列测序之后，导致边合成边测序类测序技术在序列比对时会耗费大量的计算资源以及计算时间的问题，进而提出了一种基于泛基因组的序列实时比对方法。

本发明采取的技术方案是：

它包括以下步骤：

S1、获取测序序列、线性参考基因组和变异数据库，构建泛基因组，得到构建好的泛基因组；

S11、获取测序序列，测序序列包括单端测序序列和双端测序序列；

S12、假设测序序列的长度为L，对线性参考基因组中的每条染色体进行固定长度的窗口块划分，每个窗口块长度均为2L，相邻窗口块之间重叠L个碱基；

S13、假设线性参考基因组对应的变异数据库有h个样本，则对应有2h个单体型序列，按照S12的窗口块长度对2h个单体型序列进行划分，每个窗口块内均包含2h个局部单体型序列，对每个窗口块内完全一致的序列进行去重处理，按照窗口块的顺序得到处理后的每个新单体型序列，将所有新单体型序列合并作为泛基因组；

S2、将构建好的泛基因组转换为线性序列，并与线性参考基因组合并，得到新线性序列；

S3、将新线性序列的正链序列和逆反序列合并，构成新线性参考基因组，并利用BWT索引构建算法构建新线性参考基因组的BWT索引；

S4、删除S3构建的BWT索引的SA列表上无效或者重复的信息，得到新BWT索引；

S5、获取待测序序列，根据新BWT索引对待测序序列进行比对，得到待测序序列的比对结果；

S6、将待测序序列的比对结果进行存储；

S7、根据存储的待测序序列的比对结果，获取待测序序列在泛基因组上的比对位置；

S8、根据待测序序列在泛基因组上的比对位置，选取待测序序列最优的比对位置作为待测序序列的最终比对结果；

S9、利用基数排序算法对待测序序列的最终比对结果进行排序，完成待测序序列的比对。

进一步地，S2中将构建好的泛基因组转换为线性序列，并与线性参考基因组合并，得到新线性序列，具体过程为：

将构建好的泛基因组内相邻的单体型序列首尾相接，得到一条线性序列，将此线性序列与线性参考基因合并，得到一条新线性序列。

进一步地，S5中获取待测序序列，根据新BWT索引对待测序序列进行比对，得到待测序序列的比对结果，具体过程为：

S51、定义待测序序列比对时的输入格式，得到固定格式的待测序序列，具体过程为：

(A)头文件包含每条待测序序列所属的样本ID，样本ID的范围值为(0,255]，每个样本ID使用一个字节存储；

(B)待测序序列包含M个数据，M表示单端测序序列或双端测序序列中的长度；

S52、用S4得到的新BWT索引对固定格式的待测序序列进行比对，得到比对结果；

S53、根据比对结果处理比对错误的待测序序列，得到待测序序列的最终比对结果，具体过程为：

当某个碱基比对后的比对结果的BWT区间为空区间，即表示此碱基存在错误，则检测碱基的质量分数，若质量分数小于等于设定值，遍历其他三种碱基，得到三种碱基对应的BWT区间，若得到的三个BWT区间中存在唯一非空区间，将非空区间对应的碱基作为正确碱基，非空区间作为比对结果；若质量分数大于设定值或者不存在唯一非空区间，则比对终止，将本次比对结果与S52得到的比对结果进行合并，得到待测序序列的最终比对结果。

进一步地，S51中(B)的数据包括四行：

第一行记录数据的基本信息，以@开头，包括数据的名称和注释项,名称与注释项之间使用制表符分隔；

第二行记录待测序序列在某特性位置的碱基，所有碱基按照其所属待测序序列的ID顺次排列，定义固定长度的碱基添加一个换行符；

第三行以+开头，包括数据名称和注释项，名称与注释项之间使用制表符分隔；

第四行记录待测序序列在某特性位置的质量分数值，所有质量分数值按照其所属测序序列的ID顺次排列。

进一步地，S52中用S4得到的新BWT索引对固定格式的待测序序列进行比对，得到比对结果，具体过程为：

利用S4得到新BWT索引的OCC算法对每条固定格式的待测序序列的每个碱基进行处理，OCC算法的输入为当前碱基c和当前碱基c前面的碱基串对应的BWT区间OCC_PRE，OCC算法的输出为包含当前碱基c及其前面碱基串对应的BWT区间OCC_CUR，将BWT区间OCC_CUR进行压缩并存储，同时，将BWT区间OCC_CUR作为下一个碱基比对时的OCC_PRE，直至完成固定格式的待测序序列所有碱基的比对，得到所有BWT区间OCC_CUR，将BWT区间OCC_CUR作为比对结果。

进一步地，S6中将待测序序列的比对结果进行存储，具体过程为：

每条待测序序列的最终比对结果的BWT区间均占据10个字节，当待测序序列是“单端测序”时，“单端测序”的比对结果存储在硬盘；当待测序序列是“双端测序”时，第一端测序的比对结果存储在硬盘中，第二端测序的比对结果存储在内存中。

进一步地，S7中根据存储的待测序序列的比对结果获取待测序序列在泛基因组上的比对位置，具体过程为：

根据每条待测序序列的比对结果在新BWT索引的SA列表中的位置，查找得到待测序序列在泛基因组上的比对位置，所述查找为哈希法和二分查找法。

进一步地，S8中根据待测序序列在泛基因组上的比对位置，选取待测序序列最优的比对位置作为待测序序列的最终比对结果，具体过程为：

当待测序序列是“单端测序”时，对于每条待测序序列，随机选择待测序序列在泛基因组上的任一比对位置作为最终比对结果；

当待测序序列是“双端测序”时，对于互相配对的每两条待测序序列，从第一端测序序列的比对结果中选取某位置M1，从第二端测序序列的比对结果中选取某位置M2，当M1与M2满足配对成功的条件时，双端测序序列的M1和M2两个位置配对成功，称为位置对；若得到两对或两对以上配对成功的位置对，则随机选取其中一对位置对作为双端测序序列的最终配对结果；若得到一对配对成功的位置对，则将所述位置对作为最终配对结果；若没有得到配对成功的位置对，则按照两条单端测序序列进行处理。

进一步地，配对成功的条件为：

(A)双端测序序列的其中一端相对泛基因组索引正向匹配，另一端相对泛基因组索引逆反匹配；

(B)双端测序序列中正向匹配的一端比对位置所在的窗口块的ID值小于逆反匹配的一端比对位置所在的窗口块的ID值；

双端测序序列的插入长度与插入长度的平均值之差不超过标准差的三倍。

进一步地，S9中利用基数排序算法对待测序序列的最终比对结果进行排序，完成待测序序列的比对，具体过程为：

当S8比对失败时，待测序序列的比对结果为0条；当S8比对成功时，待测序序列的比对结果为1条48bit的比对结果，将48bit的比对结果中的前24bit的比对结果按照比对结果对应的窗口块ID进行基数排序，完成待测序序列的比对。

有益效果：

本发明利用线性参考基因组和变异数据库构建泛基因组，泛基因组包含了所有常见变异的基因组，在覆盖绝大多数变异或者变异组合的前提下，节约了对内存以及计算资源的占用，能够良好地运作于实时的序列比对中。将泛基因组与线性参考基因组合并组成新线性参考基因组，构建新线性参考基因组的BWT索引，并定义待测序序列在比对时的输入格式，利用BWT索引对定义格式后的待测序序列进行比对，得到待测序序列的比对结果，使得精确比对能够跳过测序序列错误而非变异的碱基完成后续的序列比对。根据比对结果处理比对错误的碱基，利用碱基对应的质量分数恢复正确的碱基，以减少测序错误对于精确比对的影响，当待测序序列是“单端测序”时，“单端测序”的比对结果存储在硬盘；当待测序序列是“双端测序”时，第一端测序的比对结果存储在硬盘中，第二端测序的比对结果存储在内存中。根据存储的待测序序列的比对结果获取每条待测序序列在泛基因组上的比对位置，选取每条待测序序列最优的比对位置作为每条待测序序列的最终比对结果，当待测序序列是“单端测序”时，对于每条待测序序列，随机选择待测序序列在泛基因组上的任一比对位置作为最终比对结果；当待测序序列是“双端测序”时，对于互相配对的每两条待测序序列，从第一端测序序列的比对结果中选取某位置M1，从第二端测序序列的比对结果中选取某位置M2，当M1与M2满足配对成功的条件时，双端测序序列的M1和M2两个位置配对成功，称为位置对；若得到两对或两对以上配对成功的位置对，则随机选取其中一对位置对作为双端测序序列的最终配对结果；若得到一对配对成功的位置对，则将此位置对作为最终配对结果；若没有得到配对成功的位置对，则按照两条单端测序序列进行处理。利用基数排序算法对待测序序列的最终比对结果进行排序，完成待测序序列的比对。

本发明使用依赖于BWT索引的测序序列的精确比对算法对所有测序序列进行序列比对。基于BWT的精确比对算法的特性不同于传统的每条测序序列为一个整体的数据结构，而采取所有测序序列的每一个特定位置的碱基为整体的数据结构，具体为：测序仪每合成并测序一个碱基，就进行一次针对此碱基的序列比对，当测序仪测序完成最后一个碱基之后，即能够完成对于此条测序序列的完全序列比对，实现了在测序时同时进行序列比对，使得在完成测序之后的极短时间内，就能够完成全部测序序列的比对，显著地压缩了从样本开始建库到获取变异分析结果的整体时间，对于一些实时性要求较高的任务有很大的帮助，例如医学检测。

本发明提出的基于泛基因组的序列比对方法相对于基于线性基因组的序列比对方法，能够在多变异、存在结构变异时更加准确且无偏向地进行序列比对，在面向边合成边测序类技术产生的测序序列时，不仅适应超高通量，显著节约计算资源，还节省了测序和比对的时间。

附图说明

图1是本发明的流程示意图；

图2是S12中窗口块划分的示意图；

图3是S13中窗口块划分的示意图；

具体实施方式

具体实施方式一：结合图1-图3说明本实施方式，本实施方式所述一种基于泛基因组的序列实时比对方法，它包括以下步骤：

S1、获取测序序列、线性参考基因组和变异数据库，构建泛基因组，得到构建好的泛基因组。具体过程为：

根据人体DNA序列线性参考基因组以及人群DNA序列变异数据库，并针对预先设定的长度为L的DNA测序序列构建DNA序列的泛基因组，L为单端测序序列的长度或双端测序序列中其中一端的长度，常见为150碱基。泛基因组包含了所有常见变异的基因组，初始使用的数据库越大，对于新的个体的变异覆盖度就越高。

S11、获取人体DNA的测序序列，测序序列包括单端测序序列和双端测序序列。

在序列比对时，单端测序序列或双端测序序列中其中一端均执行S1-S9。

S12、假设测序序列的长度为L，对线性参考基因组中的每条染色体进行固定长度的窗口块划分，每个窗口块长度均为2L，相邻窗口块之间重叠L个碱基。

假设人体DNA的单端测序序列的长度或双端测序序列中其中一端的长度均为L，长度为L表示测序序列包含L个碱基；以人体DNA序列线性参考基因组GRCH38为例，对每条染色体进行长度为2L的窗口块划分，划分后具体的窗口块(window block，WB)大小为：第一个窗口块的范围是[1,2L]，第二个窗口块的范围是[L+1,3L]，……，第n个窗口块的范围是[(n-1)L+1,nL]，如图2所示。

S13、假设线性参考基因组对应的变异数据库有h个样本，则对应有2h个单体型序列，按照S12的窗口块长度对2h个单体型序列进行划分，每个窗口块内均包含2h个局部单体型序列，对每个窗口块内完全一致的序列进行去重处理，按照窗口块的顺序得到处理后的每个新单体型序列，将所有新单体型序列合并作为泛基因组。

大规模人群变异数据库如，国际千人基因组计划的数据库，其中包含的变异数据是线性参考基因组中部分基因数据的变异内容，假设所述数据库有h个样本，则对应有2h个单体型序列(如，人体有两套染色体，单体型序列指其中的一套染色体)，每个单体型序列均由一条局部碱基序列组成，则2h个单体型序列共有2h条碱基序列，如图3所示，每一条线代表一条单体型序列，每一个窗口块内均包含2h个局部单体型序列，不同的局部单体型序列均具有特定的碱基序列，则不同的单体型序列在同一个窗口块内的碱基序列可能一致，也可能不一致，所以需要对同一个窗口块内完全一致的局部碱基序列(局部单体型序列)进行去重，所有单体型序列都保留一个。如，在一个窗口块内，单体型1和单体型2分别有一组碱基序列(DNA序列)：碱基序列1和碱基序列2，如果碱基序列1和碱基序列2重复，就删除其中一个，直到此窗口块内2h个单体型序列没有重复为止，则此窗口块内剩余的局部单体型序列互异，如此得到每个新单体型序列，将所有新单体型序列合并组成单体型簇，按照窗口块的顺序将单体型簇排序，将排序后的单体型簇作为泛基因组。这些新单体型序列的局部碱基序列(局部单体型序列)每个都至少包含一个变异。

通过以上过程，面向边合成边测序方法的定长人体DNA测序序列的、特殊的泛基因组序列就构建完成。该泛基因组具有如下特点：

(A)如果一条待比对的测序序列包含与大规模人群变异数据库中的某单体型一致的变异序列或变异组合，则其能够在该泛基因组上找到完全的精确比对位置。

(B)除非待比对的测序序列发生全新变异(Novel mutation)，否则局部变异或者变异组合均能在大规模人群变异数据库中找到。

(C)人群的基因组在局部范围内通常仅仅包含几种特定的单体型，每种代表一种特定变异或者变异组合。因此泛基因组有极高的压缩率。

因此，所述泛基因组结构在覆盖绝大多数变异或者变异组合的前提下，节约了对内存以及计算资源的占用，能够良好地运作于实时的序列比对中。且此泛基因组基于大规模人群变异数据库(国际千人基因组计划的数据库)构建，占用内存空间小，并完全适用于短测序片段的比对。

S2、将构建好的泛基因组转换为线性序列，并与线性参考基因组合并，得到新线性序列，具体过程为：

将构建好的泛基因组内相邻的单体型序列首尾相接，得到一条线性序列，将此线性序列与线性参考基因组(如GRCH38)合并，得到更大的一条新线性序列，新线性序列也被称为群体基因组的线性展开(Linear expansion，LE)。

S3、将新线性序列LE的正链序列和逆反(或者称碱基互补)序列合并，构成长度为原始序列二倍的新线性参考基因组(Amplification linear expansion,记为LE_A)，并利用经典的BWT索引构建算法构建新线性参考基因组LE_A的BWT索引。

S4、删除S3构建的BWT索引的SA列表上无效或者重复的信息，得到新BWT索引。

BWT索引的SA列表记录着BWT索引的每一个碱基所对应的新线性参考基因组LE_A的位置，SA列表的作用是将BWT序列的某个位置在新线性参考基因组LE_A上进行定位。由于存储所有碱基的原始位置(在新线性参考基因组LE_A上的位置)所需要的内存空间过于庞大，所以经典的索引算法(如BWA)会对SA列表进行采样存储，通常32个位置或碱基存储为一个，再使用OCC查找(通常为平均32次)找到某特定BWT序列上的碱基所对应的位置信息，但是在实时系统上，该步骤的运算消耗过大，影响了序列比对的实时性，同时，SA上的一些信息是无意义或者冗余的，例如，一些比对到两条不同单体型交界地带的测序序列没有任何意义，另外，相邻的SA可能代表重复的信息，这些信息需要去重，删除无效或者重复的信息，则会使SA列表被压缩。无效定位有两种，一是测序序列只能定位到某个单体型序列内，而不能跨越不同的单体型片段；二是窗口块WB内有局部重复，定位到其中一个即可，局部重复存在的原因：人类变异是稀疏的，任何碱基上变异发生的概率大约千分之一。

由于上述问题，本发明在依赖泛基因组进行DNA序列比对的情况下提出了一种方法，对SA列表进行压缩，这使SA列表不进行采样也能维持较小的内存占用，避免了对于SA列表进行采样而导致的定位过慢的问题，并预先进行了去除无意义数据以及去重步骤。另外，SA列表会被进一步计算成为DNA序列比对的结果存储在内存中，DNA序列比对的时候可以直接读取结果，从而提高实时性。

本发明改进了BWT查找算法中的SA定位这一个步骤。每条测序序列的所有比对结果将会被预先计算好，并存储在内存中，最后一个碱基经过DNA序列比对之后直接获取所有在泛基因组上的比对位置。NGS的测序数据的长度是固定且有限的(通常为150碱基)。人体的变异虽然无穷无尽，但是常见变异(覆盖人99％以上的变异)比较固定。因此，人体(或任何特定物种)的二代测序结果是可以有限枚举的，大约数千亿个不同结果。这为结果的预先计算创造了条件。

S5、获取待测序序列，根据新BWT索引对待测序序列进行比对，得到待测序序列的比对结果，具体过程为：

在实时性的序列比对算法中，测序序列比对算法的输入和中间数据的存储与传统的算法有很大的不同。以T7测序仪为例，测序仪同时产生20G条(二百亿条)的测序片段，平均每五分钟为这20G条测序片段分别生成一个新的碱基。实时比对算法不能直接处理一条完整的测序序列，而是需要处理大量不同的测序片段某个相同位置的碱基。为了提高处理效率，所有测序序列的比对中间数据需要存储在内存中，而不能存储在硬盘上。为了解决上述问题，面向高通量实时系统，本发明设计了一种新的输入数据格式以及BWT比对过程中的中间数据压缩存储格式。

S51、定义待测序序列比对时的输入格式，得到固定格式的待测序序列，具体过程为：

通常情形下，测序仪每次测得N条测序序列的同一个位置上的碱基，连续测固定次数M，其结果为：测序仪一次运行产生的测序结果是一个N×M的碱基矩阵，传统测序序列指代矩阵的每一列，该数据经过矩阵转置，获得N条长度均为M的测序序列。但是本发明在实时比对算法中，为了加快计算速度，避免转置过程，直接使用N×M的碱基矩阵的每一行结果直接作为比对时的输入。即本发明设置了如下格式的测序序列作为DNA序列比对时的输入：

(A)头文件包含每条待测序序列所属的样本ID，样本ID的范围值为(0,255]，每个样本ID使用一个字节存储。

(B)待测序序列包含M个数据，M表示单端测序序列或双端测序序列(pair end，简称PE)的长度，双端测序序列中即为数据对的数量，以PE150为例，PE为双端测序序列的缩写，PE150为双端测序序列两边各测序150个碱基，共计150*2＝300个碱基，此时M＝300。该文件类似FASTQ格式。

每个数据包括四行：

第一行记录数据的基本信息，以@开头，包括数据的名称和注释项,名称与注释项之间使用制表符分隔。

第二行记录待测序序列在某特性位置的碱基，所有碱基按照其所属待测序序列的ID顺次排列，其顺次需要与头文件中的顺次相同。为了避免单行过长，定义固定长度的碱基添加一个换行符，所述固定长度一般为70bp。

第三行以+开头，具体内容同第一行，包括数据名称和注释项，名称与注释项之间使用制表符分隔。

第四行记录待测序序列在某特性位置的质量分数值，所有质量分数值按照其所属测序序列的ID顺次排列，即质量分数值的排序顺序与第二行的排序顺序相同，每个位置的质量分数值与第二行每个位置的碱基相对应。

S52、用S4得到的新BWT索引对固定格式的待测序序列进行比对，得到比对结果，具体过程为：

利用S4得到新BWT索引的OCC算法对每条固定格式的待测序序列的每个碱基进行处理，OCC算法的输入为当前碱基c和当前碱基c前面的碱基串对应的BWT区间OCC_PRE，OCC算法的输出为包含当前碱基c及其前面碱基串对应的BWT区间OCC_CUR，将BWT区间OCC_CUR进行压缩并存储，同时，将BWT区间OCC_CUR作为下一个碱基比对时的OCC_PRE，直至完成固定格式的待测序序列所有碱基的比对，得到所有BWT区间OCC_CUR，将BWT区间OCC_CUR作为比对结果，比对结果实质上是比对位置列表。

利用S4得到新BWT索引的OCC算法对固定格式的待测序序列的每个碱基(对应矩阵中的每一列)进行处理，所有的列都将被处理。OCC算法的输入为当前碱基c和所对应的列的BWT区间(记为OCC_PRE)，OCC算法的输出为经过处理得到的新的BWT区间(记为OCC_CUR)。计算OCC之后，将OCC_CUR存储在原来的位置上，覆盖掉原本的OCC_PRE。此时即完成在特定行，特定列的序列比对过程。之后，等待测序仪完成下一个轮次的测序。对于每一个处理到的列而言，上一轮次的OCC_CUR将作为这一个轮次比对时的OCC_PRE。每一个轮次都依照上述方案进行序列比对(即OCC)，直至完成所有的轮次(即矩阵的所有的行)，最后一个轮次生成的所有列的OCC_CUR即为比对的结果。比对结果实质上是所有列的BWT区间的列表。

BWT区间是精确比对过程中的唯一的中间数据，其需要常驻内存，每个BWT区间占用内存空间两个64位值，开始位置以及结束位置各占据8字节，共计使用16字节。在具体的实践中，BWT区间的总长度不会超过1024G，因此，BWT区间的位置值(开始位置与结束位置)可选取范围被局限在(0-2^

S53、根据比对结果处理比对错误的固定格式的待测序序列，得到固定格式的待测序序列的最终比对结果，具体过程为：

当某个碱基比对后的比对结果的BWT区间为空区间，即表示此碱基存在错误，则检测碱基的质量分数，若质量分数小于等于设定值，遍历其他三种碱基，得到三种碱基对应的BWT区间，若得到的三个BWT区间中存在唯一非空区间，将非空区间对应的碱基作为正确碱基，非空区间作为比对结果；若质量分数大于设定值或者不存在唯一非空区间，则比对终止将本次比对结果与S52得到的比对结果进行合并，得到待测序序列的最终比对结果。

测序序列错误是由测序仪产生的，测序仪依赖各种生化反应对DNA序列进行测定，因此错误在所难免。测序序列错误通常伴随着碱基质量分数低下，利用这一点尽量地减少测序错误对于精确比对的影响。如果某个碱基比对后的比对结果的BWT区间(中间数据)为空区间，即区间长度为0，且所述碱基的质量分数小于等于设定值(默认值15)，则使用如下步骤恢复正确的碱基：遍历其他的三种碱基(DNA碱基共四种，分别简称ACGT)，并获取相应碱基所对应的BWT区间，如果这三个BWT区间中存在唯一的非空区间，则以非空区间对应的碱基作为正确碱基；如果此处质量分数大于特定值(默认15)，或者不存在唯一非空区间，则在该列上比对终止，并将该列标记为比对失败。后续的其他轮次中，该列将不再进行处理。

本发明改进了基于BWT的精确比对算法，使得精确比对能够跳过测序序列错误(而非变异)的碱基，完成后续的序列比对。

S6、将待测序序列的比对结果进行存储，具体过程为：

每条待测序序列的最终比对结果的BWT区间均占据10个字节，当待测序序列是“单端测序”时，“单端测序”的比对结果存储在硬盘；当待测序序列是“双端测序”时，“第一端测序”的比对结果存储在硬盘中，“第二端测序”的比对结果存储在内存中。

每条待测序序列(矩阵的每一列)的最后一个碱基(矩阵最后一行)比对后，取得的相应BWT区间是该测序序列比对的最终比对结果。如果待测序序列是单端测序序列，则单端测序序列的比对结果存储在硬盘上，对于每条单端测序序列，其比对结果的BWT区间将会被转化成为最终的比对位置信息。如果待测序序列是双端测序序列，则第一端测序序列的比对结果被存储在硬盘上，第二端测序序列的比对结果被存储在内存里。即再具体实现时，第一端测序序列的比对结果将会从硬盘中读取出来，第二端测序序列的比对结果将会从内存中读取出来，根据第一端测序序列的比对结果读取的顺序与第二端测序序列配对，再分别获取两端的最终比对位置。

本发明提出了适用于实时分析的测序仪下机数据结构，不同于传统的每条测序序列为一个整体的数据结构(以矩阵中列为基础单位)，而采取所有测序序列的每一个特定位置的碱基为整体的数据结构(以矩阵中行为基础单位)。

S7、根据存储的待测序序列的比对结果(即矩阵每列对应的最终OCC_CUR)，获取待测序序列在泛基因组上的比对位置，具体过程为：

根据每条待测序序列的比对结果在新BWT索引的SA列表中的位置，查找得到待测序序列在泛基因组上的比对位置，所述查找为哈希法和二分查找法。

对于每个特定碱基的BWT精确比对结果产生的BWT区间，其对应的最终比对位置集合已经通过对SA列表的预处理与压缩步骤进行了预先计算。根据每条待测序序列所有碱基比对完成后得到的BWT区间，在SA列表中通过哈希法和二分查找法得到待测序序列在泛基因组上的最终比对结果。每个最终比对结果均为48bit，其中，26bit用于存储比对结果所处的窗口块ID；另外的22bit用于存储具体的单体型以及待测序序列在所述单体型内的位置。

S8、根据待测序序列在泛基因组上的比对位置，选取待测序序列最优的比对位置作为待测序序列的最终比对结果，具体过程为：

当待测序序列是“单端测序”时，对于每条待测序序列，随机选择待测序序列在泛基因组上的任一比对位置作为最终比对结果；

当待测序序列是“双端测序”时，对于互相配对的每两条待测序序列，从第一端测序序列的比对结果中选取某位置M1，从第二端测序序列的比对结果中选取某位置M2，当M1与M2满足配对成功的条件时，双端测序序列的M1和M2两个位置配对成功，若得到两对或两对以上配对成功的位置对(如M1与M2为一组位置对)，则随机选取其中一对位置对作为双端测序序列的最终配对结果；若得到一对配对成功的位置对，则将所述位置对作为最终配对结果；若没有得到配对成功的位置对，则按照两条单端测序序列进行处理。

所述配对成功的条件为：

(C)双端测序序列的其中一端相对泛基因组索引正向匹配，另一端相对泛基因组索引逆反匹配；

(D)双端测序序列中正向匹配的一端比对位置所在的窗口块WB的ID值小于逆反匹配的一端比对位置所在的窗口块WB的ID值；

(E)双端测序序列的插入长度与插入长度的平均值之差不超过标准差的三倍。

由于线性参考基因组存在重复序列，因此一条待测序序列可能精确地比对到线性参考基因组的不同位置上。对于单端测序序列而言，不同的线性参考基因组的位置会被随机选择，每个位置被选取的概率相同。对于双端测序序列而言，从第一端的结果列表中选取某位置M1，从第二端的结果列表中选取某位置M2，当M1与M2满足如上条件的时候为完美配对，在选取或得到最终的配对结果；否则，为非成功配对。

选取某位置M1或M2时，具体为：假设第一端的结果列表有M行，第二端的结果列表有N行，则一共选择M*N次，每种组合各试一次。

序列配对(pairing)：双端测序结果的两端来源于同一个小片段的两端，这提供了一些额外的信息，这些信息可以辅助双端测序序列的两端都找到最好的唯一结果。优先选取配对成功(即配对的两条序列的“比对结果“与“两条序列分别来源于小片段的两端这一事实”一致)的比对结果。

S9、利用基数排序算法对待测序序列的最终比对结果进行排序，完成待测序序列的比对，具体过程为：

当S8比对失败时，待测序序列的比对结果为0条；当S8比对成功时，待测序序列的比对结果为1条48bit的比对结果，将48bit的比对结果中的前24bit的比对结果按照比对结果对应的窗口块ID进行基数排序，完成待测序序列的比对。

每条单端测序序列或者双端测序序列的某一端在完成S8定位后，会生成0条(比对失败时)或者1条(比对成功时)48bit的比对结果，本发明将48bit的比对结果中的前24bit的比对结果按照比对结果对应的窗口块ID进行基数排序，为了减少基数排序的时间消耗，仅仅对于48bit中的前24bit进行排序。此时得到的比对结果是有序的，可以直接进行后续的变异检测步骤。本发明对于测序序列的比对结果进行标准化定长存储，并通过基数排序算法直接按照位置排好序，并后接变异检测方法直接输出变异。