CN103336916A - 一种测序序列映射方法及系统 - Google Patents

Abstract

本发明公开了一种测序序列映射方法及系统。该方法包括对所述参考基因组进行预处理，以生成参考基因组压缩结构、地址索引结构和等分点索引结构；参考基因组压缩结构以压缩形式存储整个参考基因组，地址索引结构按照一定次序存储参考基因组中所有子序列的地址值，等分点索引结构用于存储一部分地址值在地址索引结构中所处的位置，用于加速实现测序序列的初步定位；基于参考基因组的特征、测序序列的整体信息、测序序列所属物种的多态发生情况，根据概率设计映射算法参数，以达到或折中对灵敏度、特异度、映射速度的要求；根据预处理得到的上述结构，通过子序列定位、基于自匹配函数的延拓、定量分析步骤将每一个测序序列映射至所述参考基因组。

Description

一种测序序列映射方法及系统

技术领域

本发明适用于基因工程技术领域，尤其涉及一种基于高通量测序技术的测序序列快速映射及有关定量分析的方法和系统。

背景技术

高通量DNA测序是实现个体化医疗和开展现代分子生物学研究的核心技术。在个体化医疗中，高通量DNA测序可以获得一个人的全基因组、表达组、以及各种调控分子的定性和定量信息，可以综合利用遗传序列中的多态和变异信息、功能性基因组学中的表达信息，从分子水平上实现疾病诊断，患病风险预测，从而更好地进行治疗或预防。特别地，我们可以根据个人的遗传序列和功能性基因组信息预测药物对于个体的影响程度，并基于此设计最佳的治疗方案。

除了人类健康，农业、环境、能源等对人类生活至关重要的发展都离不开我们对生物学在分子层面上的全面认识。而分子生物学研究的一个主要手段就是DNA测序。

在对基因组进行测序时，基因组被切割成很多小片段，通过复制、碱基辨识等步骤，我们可以获得这些短序列的碱基序列(测序序列)。然而在切割基因组后，我们无法知道各个测序序列的相对位置。如果没有参考基因组，就只能通过装配技术来得到所测的基因组。如果有一个已被测得的基因组作为参照，这就是一个相对容易的重测序问题。现在我们在生物学研究、个体化医疗中面临的测序问题，绝大部分是或可以近似转化为重测序问题。在重测序问题中，我们要寻找每一个测序序列在参考基因组上的位置或坐标，我们称之为测序序列映射。

下面列出DNA重测序可以测量的分子生物信息的若干主要类型。

DNA多态——通过将基因组测序序列与参考基因组序列进行比较探寻两个基因组之间存在的单点多态(SNP)以及插入/删失多态；

甲基化——通过将基因组测序序列映射至参考基因组序列，来探寻原基因组中的甲基化位点；

mRNA表达谱——通过将转录组的测序序列映射至参考基因组序列来测量不同种类的RNA的含量；

可变剪切——通过将转录组的测序序列映射至参考基因组序列来探测mRNA可变剪切的模式；

非编码RNA丰度——通过将转录组的测序序列映射至参考基因组序列来测量不同种类的非编码RNA的含量，如microRNA和lncRNAs；

CHIP-seq——通过染色质免疫共沉淀技术(ChIP)特异性地富集目的蛋白结合的DNA片段，将这些DNA片段的高通量测序数据映射至参考基因组序列上，从而获得全基因组范围内与组蛋白、转录因子等互作的DNA区段信息。

DNA重测序的用途广泛，不仅限于以上所列。相应地，测序序列映射也就成为个体化医疗以及分子生物学研究中必不可少的、日常的计算分析工作。虽然映射的概念很清楚，但是高通量的新一代测序技术可以在短时间内产生海量的测序序列，如何能够运用相对通用的计算机设备高速地完成映射工作，是一个非常有挑战的计算生物问题。

除了映射速度，正确评估一个映射方法和系统的映射率和准确率是测序数据下游分析的基石。映射率是指能够映射到参考基因组上的序列比例，而准确率是指可映射序列中的正确映射的比例。有时我们也用另一对指标灵敏度和特异度来代替映射率和准确率。在很多情况下，由于技术的局限性，基于特定的数据，我们往往不能同时提高灵敏度和特异度，这时，如何在灵敏度和特异度之间找到合适的平衡也是一个具有高度挑战的问题。

我们的发明的目的就是设计一个计算系统，具有以下功能：

1.它可以高速地映射重测序序列；

2.对于给定的数据特征如读长和质量值，和一组系统参数数值，它可以评估映射的灵敏度和特异度；

3.对于特定的生物和医学问题，基于给定的数据特征如读长和质量值，设计系统参数数值，以达到合理的映射的灵敏度和特异度。如有余地，还可以通过参数设计优化映射速度；

4.这个系统是模块化的。

5.对于不同的模块，用不同的硬件和软件实现。其中的硬件可以是通用的CPU，RAM，硬盘存储器的组合，如工作站、服务器，也可以包括GPU、FPGA(可重构计算单元)、DSP等。

测序序列映射是将测序平台给出的测序序列定位至参考基因组上这一过程。由于测序序列映射问题对新一代测序数据处理的重要性，人们构造了很多映射工具。这些映射工具都采用了多个子序列同时映射的方法，即使用多个子序列同时对测序序列进行初步定位，有些方法还允许用于定位的子序列内存在若干个不匹配以提高序列映射的准确率。但这样的方法往往需要占据很大的内存，且每一个用于定位的子序列的长度有限；另外，所允许的不匹配数越多，算法的复杂性越大。早先的新一代测序技术所产生的测序序列长度有限，因此这样的方法是合理的。但是随着新一代测序技术的发展，测序序列从最初的只有二十几个碱基到现在可以达到上百的长度，使用多个子序列映射方法，或者允许子序列内有不匹配会降低序列映射的速度。同时，为了找到插入/删失，现有的映射工具均使用动态规划算法(Smith-Waterman算法或其变种)，其时间复杂度和空间复杂度对于测序序列长度都是平方级的。

在具体实施上，现有的方法采用集群计算机进行映射，数据量大时需要几天的时间才能完成。如果能够采用单一服务器或工作站，以几倍的速度完成映射工作，将大大降低维护成本，并且提高操作的灵活性(在同样的时间内可以调整参数进行多次映射)。

此外，测序序列映射的定量分析是基于新一代测序技术科学研究必不可少的，而目前，除了MAQ这一序列映射工具提供了简单的定量分析外，其它映射工具都缺少参数设计准则以及映射率、准确率的评估。

以下为相关名词解释：

新一代测序技术：英文名称为Next Generation Sequencing，英文缩写为NGS，此技术利用若干生物化学技术，将被测基因组序列切割成多个子序列，并确定每个子序列碱基组成.其特点是能并行地产生大量的测序序列数据。

测序序列：由测序平台输出的某一生物体基因组子序列的测量结果，由若干个代表四种碱基的字符(A，C，G，T)组成，代表基因组的一个子序列。

bp：英语basepair的缩写。100bp的序列指一个DNA序列有100个碱基长。

碱基误读：由于测序技术的误差使得测序序列上某些碱基不同于真实碱基。

质量值：在测序过程中，测序序列的每一个碱基都有被误读的可能性，质量值是这一可能性的反应。测序序列的每一个碱基都对应一个这样的质量值，质量值越高，被误读的可能性越低。

参考基因组：被测序的生物体所属物种的一个已完成测序的基因组序列。

测序序列映射：寻找测序序列在参考基因组上所处的位置或坐标。

n位子序列：参考基因组或者测序序列上长度为n的子序列。

完全匹配前缀子序列：测序序列上的用于对其进行初步定位的子序列的前缀序列，该前缀序列与参考基因组上某些等长的子序列完全相同。

碱基替换：测序序列被映射到参考基因组上后某些碱基和参考基因组上对应的碱基不相同。

插入：指相对于参考基因组，测序序列的某两个相邻的碱基之间额外插入了一段碱基序列。

删失：指相对于参考基因组，测序序列丢失一段或若干段碱基序列。

碱基值：一个碱基序列被编码后所得的整数值，如用二进制码00，01，10，11分别编码A，C，G，T，则序列ATGCTA的碱基值为001110011100。

地址值：参考基因组子序列的起始位点，如果参考基因组为ATGTAGCTA，那么子序列ATGTA的地址值为0，子序列AGCTA的地址值为4，地址值为2的4-字词为GTAG，地址值为6的3-字词为CTA。

发明内容

本发明针对现有技术所存在的问题，基于以下(1)、(2)、(3)所述三个事实，提供了一种针对新一代测序技术的高通量数据映射系统，旨在提高序列映射的速度，并对序列映射的准确性予以定量分析：

(1)近年来新一代测序技术不断改进，测序平台输出的序列长度比以前有所增长，可以达到100bp或更长，且碱基误读概率也在下降。例如，Illumina测序技术进行测序时所发生的错误基本上都是碱基误读，出现插入/删失很少；

(2)对于如人类基因组这样的保守基因组，在一个长度为100bp的子序列中出现两个或以上的插入/删失的可能性小于0.01％，而出现的每一个插入/删失的平均长度只有3；

(3)科研用计算机的内存大幅度增加，使得以占用更多内存来换取序列映射速度的提高成为可能。

为此，本发明提供了一种测序序列映射方法，对于所获取参考基因组和至少一个测序序列进行操作；其中，参考基因组为已完成测序的基因组序列；操作包括以下步骤：

步骤1、对所述参考基因组进行预处理，以生成参考基因组压缩结构、地址索引结构和等分点索引结构；所述参考基因组压缩结构以压缩形式存储整个参考基因组，所述地址索引结构按照一定次序存储所述参考基因组中所有子序列的地址值，所述等分点索引结构用于存储一部分地址值在地址索引结构中所处的位置，用于加速实现测序序列的初步定位；

步骤2、基于参考基因组的特征、测序序列的整体信息、测序序列所属物种的多态发生情况，根据概率计算，设计映射算法的参数，以达到或折中对灵敏度和特异度、映射速度的要求；

步骤3、根据经过预处理后得到的上述三种结构，通过子序列定位、基于自匹配函数的延拓、定量分析步骤将每一个测序序列映射至所述参考基因组上；

步骤4、输出每个测序序列的映射信息。

本发明还提供了一种测序序列映射系统，该系统对于所获取参考基因组和至少一个测序序列进行操作；其中，参考基因组为已完成测序的基因组序列；该系统包括：

参考基因组预处理模块，其用于对所述参考基因组进行预处理，以生成参考基因组压缩结构、地址索引结构和等分点索引结构；所述参考基因组压缩结构以压缩形式存储整个参考基因组，所述地址索引结构按照一定次序存储所述参考基因组中所有子序列的地址值，所述等分点索引结构用于存储一部分地址值在地址索引结构中所处的位置，用于加速实现测序序列的初步定位；

参数设计模块，其用于基于参考基因组的特征、测序序列的整体信息、测序序列所属物种的多态发生情况，根据概率计算，设计映射算法的参数，以达到或折中对灵敏度和特异度、映射速度的要求；

映射模块，其用于根据经过预处理后得到的上述三种结构，通过子序列定位、基于自匹配函数的延拓、定量分析步骤将每一个测序序列映射至所述参考基因组上；

结果输出模块，其用于输出每一个测序序列的映射信息。

本发明提出的上述方案经过序列映射以及有关定量分析评估的处理，输出每一个测序序列是否被成功映射、被成功映射的测序序列在参考基因组上所处位置，以及测序序列和参考基因组的比对情况等信息。在对测序序列进行映射操作之前，本发明先对所输入的参考基因组进行预处理，之后利用预处理后所得的各个结构，采用同样的方法并行地实现对每一个测序序列的映射和定量分析。

附图说明

图1为本发明中测序序列映射方法流程图；

图2为本发明中将人体基因组全体32位子序列中的A、C、G、T分别编码为00、01、10、11后得到的碱基值在排序后与均匀分布的关系示意图；图中纵轴为每个碱基值，横轴为均匀分布的等分点数值；

图3为本发明中对参考基因组进行预处理得到地址索引结构的操作示意图；

图4为本发明中选用16位无符号二进制数存储20位二进制地址值的方法示意图；

图5为本发明中选用16位无符号二进制数存储28位二进制地址值的示意图；

图6为本发明中等分点索引结构构造示意图；

图7为本发明中映射一个测序序列的方法流程图；

图8为本发明中采用等分点索引结构缩小查找范围的方法示意图；

图9为本发明中采用二分法对测序序列上的子序列进行初步定位的实施过程示意图；

图10为本发明中匹配向量的定义示意图；

图11为本发明中使用自匹配函数方法探测插入/删失的长度示意图；

图12为本发明中探测插入/删失的起始位置示意图；

图13为本发明中探测插入/删失的长度以及起始位置的比选方案示意图。

具体实施方式

为使本发明的目的、技术方案和优点更加清楚明白，以下结合具体实施例，并参照附图，对本发明作进一步的详细说明。

本发明提供了一种基于高通量测序技术的测序序列快速映射方法。

图1示出了本发明提出的基于高通量测序技术的测序序列快速映射方法流程图。

该方法的输入包括参考基因组和由测序平台给出的一个或多个测序序列数据集。所述参考基因组和测序序列由代表四种碱基的字符(A，C，G，T)组成；参考基因组可以是任意一个已完成测序的物种的基因组序列；所述测序序列与参考基因组应属于同一物种或相近物种。所述测序序列可以包含表示测序序列中每个碱基被误识别的可能性的质量值信息，所输入的测序序列的长度可以不同。如图1所示，该方法包括以下步骤：

步骤1、对所述参考基因组进行预处理的步骤。

首先用二进制码00、01、10、11的某个排列来编码四种碱基A、C、G、T。编码规则可以任意选择，也可以通过更细致的分析进行选择。然后建立参考基因组压缩结构、地址索引结构、等分点索引结构三种结构。其中，所述参考基因组压缩结构用于存储参考基因组碱基序列的二进制码，每个字节存储四个碱基，用于查找参考基因组上从给定位置开始的序列片段或其碱基值。所述地址索引结构按照一定的顺序存储参考基因组中每个n位子序列对应的地址值，所述一定的顺序为所述每个n位子序列碱基值的大小顺序，其如下生成：将参考基因组中全部n位子序列的碱基值从小到大排序后，在一个向量中依次存储它们在参考基因组中的地址值。所述等分点索引结构用于存储一部分地址值在地址索引结构中所处的位置，用于加速实现测序序列在参考基因组中的初步定位。对于给定的参考基因组，预处理只需做一次，结果保留在存储系统中，可以对不同的测序序列数据集反复使用。

步骤2、参数设计步骤。

基于基因组的特征、测序序列的整体信息、测序序列所属物种基因组的多态发生情况，根据概率计算，设计映射算法的参数，以达到或折中对灵敏度、特异度、映射速度的要求。其中，所述测序序列的整体信息包括全体测序序列的长度分布情况和总体质量值分布情况；所述测序序列所属物种基因组的多态发生情况包括SNP、插入/删失发生的频率和插入/删失的长度分布情况。在所有参数中，完全匹配前缀子序列长度下界k₀是主要参数，它依赖于参考基因组的长度、碱基频率和多态率、以及测序序列长度和质量值。步骤1中构造地址索引结构时所用的参数n要大于k₀。

步骤3、映射步骤。

基于经过预处理后得到的三种结构，或者只基于参考基因组压缩结构和地址索引结构两种结构，对给定的每一个测序序列进行映射，即寻找所述每一个测序序列在参考基因组上所处的位置；具体步骤包括“子序列定位”、“基于自匹配函数的延拓”、“定量分析”三个子步骤。

子序列定位子步骤：按照设定的次序逐个扫描测序序列的n位子序列，直到找到某个n位子序列，它有一个长度大于k₀的前缀与参考基因组的某一子序列完全匹配，用这一个参考基因组的子序列的地址作为测序序列的初步定位。

基于自匹配函数的延拓子步骤：利用或不用测序质量值，使用本发明中新提出的自匹配函数方法将测序序列与参考基因组之间的差异分为三种情况：(1)只有碱基替换；(2)在完全匹配前缀子序列的一侧只有一个插入/删失；(3)其它更复杂的情况。对(2)，进一步用自匹配方法确定插入/删失的类型、位置和长度。对(3)，用比对算法确定差异的细节。特别地，采用局部-全局比对算法确定测序序列上被延拓的片段(全局)和参考基因组相应DNA片段(局部)之间的差异细节。在平均意义下，这个方法对于测序序列长度具有线性的时间复杂度。

定量分析子步骤：对于测序序列上完全匹配前缀子序列以外部分，根据所属生物物种基因组的多态率和每个碱基辨识的质量值，将每一个碱基与参考基因组上相应碱基的匹配状态以一个“软分值”表示，并利用概率分布来评估，以决定是否需要选择新的n位子序列重新映射测序序列。

步骤4、反馈步骤。

在对所有测序序列完成一轮映射后，对未能映射成功的测序序列，或者对敏感、特异度较差的测序序列，通过调整映射参数重新进行映射；其中所述映射参数包括步骤2中所设计的全部或部分参数。

步骤5、输出步骤。

输出每个测序序列的映射信息，将映射信息写入具有标准格式的文件。

上述测序序列映射方法适用于生物分析和个体化医疗中DNA遗传信息和RNA功能信息的定性或定量测量，这包括DNA多态、甲基化、信使RNA表达谱、非编码RNA丰度、CHIP-seq、可变剪切等信息。

以下详细阐述本发明提出的方法和原理。首先详细阐述映射非甲基化测序序列的方法和原理，对于甲基化测序序列的映射，只需对个别步骤进行修改，其它步骤与非甲基化测序序列的映射方法相同，具体的区别将在后面阐述。

本发明提出的上述方法的步骤1中预处理阶段通过对参考基因组进行操作，生成映射所用到的三个结构：参考基因组压缩结构、地址索引结构、等分点索引结构，步骤1具体包括以下内容：

步骤11、构造参考基因组压缩结构。参考基因组压缩结构以压缩的形式存储参考基因组，用于查寻参考基因组中从给定地址值开始的子序列或者其碱基值，其构造包含以下几个步骤：

步骤111：设定每个被编码的子序列长度n，n为正整数，优选地将其范围设为5≤n≤32)，对于碱基字符集合{A，C，G，T}与二进制码集合{00，01，10，11}之间的每一种映射方式，进行以下操作：

(1)对参考基因组的每一个n位子序列，将其中表示碱基的字符A，C，G，T按照给定的映射方式逐一映射为其所对应的二进制码，得到一个长度为2n的二进制码串，该二进制码串表示一个二进制整数值，该二进制整数值即为相应的n位子序列的碱基值；其中，对于参考基因组上的每个碱基，都要取出以其为起始点的n位子序列，且只要有一个碱基不同，就属于不同的n位子序列；

所述n位子序列为参考基因组中长度为n的子序列；所述给定的映射方式为碱基字符集合{A，C，G，T}与二进制集合{00，01，10，11}之间的24种一一映射法则之一；本发明中通过该24种映射方式将参考基因组中的每个子序列进行映射，并得到每一种映射方式下的n位子序列的碱基值序列；

(2)将全体碱基值从小到大排序。

步骤112：将每一个经排序的碱基值序列与[0，4ⁿ]上的均匀分布进行比较，以最接近均匀分布的碱基值序列所采用的映射方式作为编码规则；即在步骤111中得到的所有经过排序的碱基值序列中，找到最符合[0，4ⁿ]上的均匀分布的碱基值序列，并将该碱基值序列对应的映射方式作为最终的编码规则；

图2示出了对人体基因组，取n为32，将A、C、G、T分别编码为00、01、10、11后得到的碱基值在排序后与均匀分布的关系。图中纵轴为每个碱基值，横轴为均匀分布的等分点数值。

步骤113：依照步骤112所确定的编码规则，利用二进制码将参考基因组序列按照从左到右，或者从右到左的方向逐位存入到一个字节类型的向量中，每个字节可以存四个碱基；

步骤12、构造地址索引结构。

地址索引结构按一定次序存储参考基因组中所有n位子序列的地址值，用于将测序序列的子序列初步定位至参考基因组，其中所述n位子序列为参考基因组中长度为n的子序列，n为构造参考基因组压缩结构时设定的被编码的子序列长度。其构造包含以下步骤：

步骤121：对于111(1)中按照步骤112中确定的编码规则编码所得的碱基值序列，将每一个碱基值与该n位子序列在参考基因组中的地址值相关联，该n位子序列的碱基值与地址值共同组成一关联单元；

步骤122：将全体关联单元按照碱基值从小到大排序，碱基值相同的关联单元，按照地址值从小到大排序；

步骤123：根据参考基因组的长度，选用一个或多个无符号的16位，32位，或者64位整数存储排序后的关联单元中的全体地址值。这些无符号整数所占用的内存空间形成地址索引结构。

图3示出了本发明中对参考基因组进行预处理得到地址索引结构所进行的操作。如图3所示，以n＝8为例，示出了子序列截取、碱基编码以及地址索引结构的构造等操作。其中，“索引值”一列记录每一个地址值在地址索引结构中的位置，编码规则为：将A、C、G、T分别编码为00、01、10、11；

步骤122中，也可以将全体关联单元按照碱基值从大到小排序，碱基值相同的关联单元，按照地址值从大到小排序；后文所提到的地址索引结构都是按照从小到大排序的方法生成的，如果按照从大到小排序的方法生成地址索引结构，映射方法类似。

步骤123中，优选地，可以根据下面的规则存储地址值，其中b为不小于log₂N的最小整数，N为参考基因组的长度：

对于下面所述(1)至(3)种情况，选用16位无符号二进制数进行存储；对(4)和(5)种情况，选用32位无符号二进制数进行存储：

(1)b≤16时，将每个地址值顺次赋值给相应的每一个16位无符号二进制数；

(2)16＜b≤24时，令每c+1个二进制数存储c个地址值，其中前c个二进制数依次存储每一个地址值的前16位，第c+1个二进制数存储全部c个地址值的后b-16位；

为不大于

的最大整数；

图4示出了b＝20时，选用16位无符号二进制数存储地址值的方法。如图4所示，左边为需要进行存储的地址值，而右边为使用16位无符号二进制数进行存储的情况：前4个16位无符号二进制数分别存储4个地址值的前16位，最后一个16位无符号二进制数存储4个地址值的后4位；最后一个16位无符号二进制数可以和前4个16为无符号二进制数分开存储。

(3)24＜b≤28时，令

每

个二进制数存储c个地址值，其中前c个二进制数依次存储每一个地址值的前16位，此后的

个二进制数存储每两个地址值分别从第17位开始的8位，最后一个二进制数存储全部c个地址值的后b-24位；

为不小于的最小整数；

图5示出了b＝28时，选用16位无符号二进制数存储地址值的示意图，其中每7个16位无符号二进制数存储4个地址值。

(4)28＜b≤31时，令

每c个二进制数存储c+1个地址值，前c个地址值顺次存储在c个二进制数的前b位，第c+1个地址存储在所剩下的(32-b)c个bit位中；

(5)b＝32时，则使用32位无符号二进制数存储地址值，将地址值顺次赋值给每一个a位无符号二进制数；

(6)b＞32时，则将地址值分组存储，每一组最多存储2³²个整数，对每组内的地址值，选取适当位数的无符号整数将其除以2³²后所得的余数按照(1)至(5)所描述的方法进行存储。

对于上述第(4)种情况，可以使用无符号的32位整数存储地址值；对于(6)，可以使用无符号的64位整数存储地址值。

为了清楚地介绍映射步骤，引入一些必要的记号。将地址索引结构中所存储的第j个地址值记为a[j]，在实施时，可根据步骤123的存储方式，利用字位操作从地址索引结构中获取a[j]的值。将参考基因组上地址值为a[j]，即从第a[j]个碱基开始的长度为n的碱基序列记为S(a[j]，n)，将S(a[j]，n)称为地址值a[j]所对应的n位子序列，或者称为参考基因组上从地址值a[j]开始的n位子序列。将S(a[j]，n)按照步骤112中确定的编码规则编码之后所得的碱基值记为e(S(a[j]，n))。在实施时，获取a[j]的值后，可根据步骤113的存储方式，利用字位操作从参考基因组压缩结构中获取碱基值e(S(a[j]，n))的值。

步骤13、构造等分点索引结构。等分点索引结构用于加快映射子序列的速度，在一些情况下，比如按照步骤112中确定的编码方式编码并排序后所得到的碱基值序列与均匀分布偏离很大，在映射时也可以不用这个结构。其构造包含以下步骤：

步骤131：建立一维数组按从小到大的顺序存储区间[0，4ⁿ-1]的均匀分布等分点z_i＝i×2^2n-c，i＝0，1，2，...，2^c-1，其中n是步骤111中设定的被编码的子序列长度，c(1≤c≤2n)为预先设定的整数值，使得在定位每一个n位子序列时地址索引结构的访问次数平均可以减少大约c次；如果内存足够大，而且

则可以优选地将c设为不大于

的最大整数，或不小于

的最小整数，其中N为参考基因组的长度。

步骤132：对于步骤131中得到的每一个等分点z_i，利用二分法在地址索引结构中，找到两个相邻的地址值a[t_i]和a[t_i+1]，使得等分点z_i的大小介于碱基值e(S(a[t_i]，n))和e(S(a[t_i+1])，n)之间，即e(S(a[t_i]，n))＜z_i＜e(S(a[t_i+1])，n)，记录地址值的索引值t_i，使之与等分点z_i相对应，即将索引值t_i记录为等分点索引结构的第i个值；如果地址索引结构中存在一个或若干个地址值，它们所对应的n位子序列的碱基值等于z_i，则记这些地址值的索引值中最小的一个为t_i，使之与等分点z_i相对应，即将索引值t_i记录为等分点索引结构的第i个值；

步骤133：将步骤132中针对每一个等分点z_i所记录的数值t_i存储在一个数组中：q[i]，i＝0，1，2，...，2^c，在映射测序序列之前，先将数组q[i]置于内存中。

图6示出了本发明中构造等分点索引结构的示意图。如图6所示，等分点z_i的大小介于碱基值e(S(a[t_i]，n))和e(S(a[t_i+1]，n))之间，将t_i存储为等分点索引结构的第i个值。

下面详细阐述映射方法的步骤2中的算法参数设计。

算法中的参数包括：

●基因组预处理中的步骤111中被编码的子序列长度n；

●测序序列的每一个n位子序列与参考基因组的最大完全匹配前缀长度的下限值k₀，这是关键参数，映射的灵敏度和特异度主要由它决定，k₀要小于n。

●所允许的参考基因组上可以被完全匹配前缀子序列定位的位置数目的最大值G，这是处理基因组上同源序列所需要的一个参数，由该物种基因组特征决定。对人类基因组，G可以选择3、5等等。

●每个测序序列的子序列搜寻数目上限值U，这个参数越大，算法越接近全部子序列搜寻方案，但也越慢。

●可以被探测的插入/删失的长度的最大值B，这个由该物种的群体基因组学特征决定。对人类基因组，可以选择5等。

●定量分析中的检验水平α，映射的灵敏度和特异度部分由它决定。

需要设计的参数主要是k₀、α、U。此外，算法还需要预先得到参考基因组和测序数据集的以下信息。

●该参考基因组的总长度N。

●该参考基因组A，C，G，T四种碱基的组成频率，记它们的平方和为ε，其含义为：以这四种碱基出现的频率为概率分布，独立地做两次有放回的随机抽样试验，所抽到的两个碱基相同的概率。

●该物种基因组中多态的频率γ。

●测序序列的长度L，或所有测序序列的长度分布。

●测序序列每个碱基的平均误读概率β，由测序序列中碱基辨识的质量值得到。

选择k₀、α，是为了满足对映射的灵敏度和特异度的要求。这里灵敏度定义为某个测序序列被映射到正确地方(测序序列所来的地方)的概率，特异度定义为该测序序列没有被映射到参考基因组上其它地方(除重复区域和高度保守的同源区域)的概率。下面阐述选择k₀、α的步骤：

(1)定义以下函数

其中L代表测序序列的长度，k代表完全匹配前缀子序列的长度，θ为碱基匹配概率，为max{s；L-ks≥g}与(g+1)的最小值，g为不匹配个数，w是长度至少为k的完全匹配前缀子序列的个数，s是w所能取到的最大值。

(2)将θ的值取为1-β-γ，计算灵敏度如下：Ψ(L，k，1-β-γ)×(1-α)；选取k^*和α使得灵敏度大于给定的值，如99％、95％、90％、80％等等。

(3)将θ的值取为ε，计算特异度如下：1-min{NΨ(L，k，ε)，1}；选取k^**使得NΨ(L，k，ε)小于给定的值，如1、0.5、0.1等等。

(4)若k^*小于k^**，需要调整α以及对灵敏度和特异度的要求，使得k^*大于k^**。k^*大于k^**时，k₀可选为闭区间[k^**，k^*]中的任意数值。

以上是基于概率模型选择k₀、α的方法。在此基础上，可以通过对测序序列数据集随机抽样，进一步优化参数。具体的操作为：(1)取一组k₀、α和充分大的U；(2)随机抽取一定量的测序序列；(3)对每一个被抽取的测序序列，按照给定参数的算法映射每一个被抽取的测序序列，如映射成功，记录已被搜寻的子序列数目；(4)根据(3)中已被搜寻的子序列数目的分布情况，调整子序列搜寻数目上限值U；(5)根据映射率，对k₀、α做必要调整。

本发明提出的上述方法的步骤3中，利用预处理得到的三个结构，或者只利用参考基因组压缩结构和地址索引结构，按照子序列定位、基于自匹配函数的延拓、定量分析的步骤将每一个测序序列映射至参考基因组。

图7示出了本发明中映射一个测序序列的方法流程图。如图7所示，该方法具体包括：

步骤31、子序列定位步骤。

在被扫描的子序列数目不大于U的情况下，按照设定的次序逐个扫描测序序列的n位子序列，直到找到某个n位子序列，它有一个长度大于k₀的前缀与参考基因组的某一子序列完全匹配。用这一个参考基因组的子序列的地址作为测序序列的初步定位，该子序列所处的位置称为参考基因组中对应于完全匹配前缀子序列的位置，其中完全匹配前缀子序列是指与所述参考基因组的子序列完全匹配的n位子序列的前缀。如果参考基因组上的长度大于k₀、且与n位子序列的前缀完全匹配的子序列数目小于G，则进入步骤32，否则继续扫描下一个n位子序列；如果被扫描的子序列数目大于U，则判定映射失败，结束步骤3的操作。

步骤32、基于自匹配函数的延拓步骤。

利用或不用测序序列的质量值，使用自匹配函数方法将测序序列上位于完全匹配前缀子序列左侧或右侧的片段与参考基因组中对应于完全匹配前缀子序列的位置的左侧或右侧片段之间的差异分为三种情况：(1)只有碱基替换，即两个片段之间只存在碱基不匹配；(2)完全匹配前缀子序列的一侧只有一个插入/删失，即从测序序列片段上的某一个位置开始，测序序列片段相对于参考基因组片段增加或减少了一个或连续若干个碱基；(3)其它更复杂的情况。对(1)，已经完成了差异的寻找；对(2)，进一步用自匹配函数方法确定插入/删失的类型、位置和长度；对(3)，用比对算法确定差异的细节。特别地，采用局部-全局比对算法确定测序序列延拓部分(全局)和参考基因组相应DNA片段(局部)之间的差异细节。在平均意义下，这个方法对于测序序列长度具有线性的时间复杂度。

步骤33、定量分析步骤。

对于测序序列上完全匹配前缀子序列以外的部分，根据所属生物物种基因组的多态率和每个碱基辨识的质量值，将每一个碱基与参考基因组上相应碱基的匹配状态以一个“软分值”表示，并利用概率分布来评估，以判断是否接受这些碱基的匹配状态；如果接受，则判定映射成功；如果拒绝，则返回步骤31。

上述步骤31具体包括以下步骤：

步骤311：将测序序列上被扫描的n位子序列按照步骤112中确定的编码规则进行编码，得到该n位子序列的碱基值；

步骤312：进行二分法查找初始化，即根据所述n位子序列的碱基值得到第V个等分点z_V，V为不大于n位子序列的碱基值除以2^c的商的最大整数，并得到等分点索引结构的第V个地址索引值t_V，然后将下限起始位置lower设置为t_V，上限起始位置upper设置为t_V+1+1；

图8示出了上述采用等分点索引结构缩小查找范围的方法示意图。

步骤313：将n位子序列的完全匹配前缀长度k设置为0；

步骤314：如果upper-lower＞1，则进入步骤315，否则进入步骤318；

步骤315：记f为S(a[lower+k]，n-k)和S(a[upper+k]，n-k)的最大公共前缀的长度，将k重新设置为k+f；

步骤316：记在参考基因组上查寻碱基值e(S(a[middle+k]，n-k))；

图9示出了上述二分法的实施过程。如图9所示，其中以n＝10为例，如果碱基值e(S(a[upper]，10))和e(S(a[lower]，10))这两个二进制数的前8位一样，则说明它们所对应的10位子序列的前4位是完全相同的，因此只需比较e(S(a[middle+4]，6))和测序序列上的10位子序列的6位后缀的碱基值。

步骤317：比较所述n位子序列的n-k位后缀的碱基值与步骤316中得到的碱基值e(S(a[middle+k]，n-k))，如果n位子序列的n-k位后缀的碱基值大于或等于e(S(a[middle+k]，n-k))，则将lower的值重新设置成middle并进入步骤314，否则将upper的值重新设置为middle并进入步骤314；

步骤318：比较子序列S(a[upper+k]，n-k)与n位子序列的n-k位后缀，记h_u为两者所共有的最长前缀长度，记k_u为k+h_u；比较子序列S(a[lower+k]，n-k)与n位子序列的n-k后缀，记h_l为两者所共有的最长前缀长度，记k_l为k+h_l；

步骤319：取完全匹配前缀子序列为n位子序列的k_l位前缀与k_u位前缀中较长的一个，如果完全匹配前缀子序列的长度小于k₀，则在测序序列上扫描下一个n位子序列并返回步骤311重新进行上述操作；否则利用地址索引结构和参考基因组压缩结构获取参考基因组上可以被完全匹配前缀子序列定位的位置数目，即参考基因组上与完全匹配前缀子序列相同的子序列的数目，如果该数目大于G，则在测序序列上扫描下一个n位子序列并返回步骤311重新进行上述操作，如果该数目小于或等于G，则结束操作。

下面以寻找完全匹配前缀子序列右侧的序列与参考基因组的差异为例详细阐述步骤32的操作过程，对于完全匹配前缀子序列左侧的序列，可以类似地进行操作。

在描述这部分之前，先对所需要的符号和变量进行介绍：

记测序序列上完全匹配前缀子序列右侧的片段为X，其长度记为l_X；在子序列定位阶段，可以在参考基因组上找到一个子序列，其与完全匹配前缀子序列完全一样，记这一子序列右侧的一个长度大于l_X的片段为Y，其长度记为l_Y。记X[i]为X的第i个碱基，记Y[i]为Y的第i个碱基。i＞0时，记X{i}为X的l_X-i位后缀，记Y{i}为Y的l_Y-i位后缀。以向量V(X，Y)表示X和Y左端对齐时每一个位置的匹配信息，V(X，Y)的长度为l_X的长度与l_Y的最小值。记X[i]的质量值为QA_X[i]。V(X，Y)的第i个元素的定义方式如下：

将质量值QA_X[i]转化为碱基误读概率β_X[i]；如果X[i]和Y[i]相同，则将V(X，Y)的第i个元素设置为β_X[i]，否则设置为1-β_X[i]。

在上述定义中，将质量值QA_X[i]转化为碱基误读概率β_X[i]时，可以使用公式

也可以在换算之前对质量值进行修正以得到更准确的误读概率。如果测序序列没有质量值信息或者质量值的可信度不高，则可以采用如下定义方式：如果X[i]和Y[i]相同，则将V(X，Y)的第i个元素设置为0，否则设置为1。

对于X和Y，以M(i)表示它们的匹配向量，其定义如下：

M(0)＝V(X，Y)；

M(i)＝V(X，Y{i})，若i＞0；

M(i)＝V(X{-i}，Y)，若i＜0；

其中，|i|＜min(l_X，l_Y)，M(0)表示片段X和Y中每一个碱基对的不匹配程度；若i＞0，则M(i)表示片段X与片段Y的l_Y-i位后缀之间的不匹配程度；若i＜0，则M(i)表示片段X的l_X-|i|位后缀与片段Y之间的不匹配程度。

图10示出了所述匹配向量的一个具体实例，构造匹配向量时所采用的向量V(X，Y)，V(X，Y{1})，V(X{1}，Y)采用简单的定义，其元素只包含0和1。

对于X和Y，定义自匹配向量w(i)如下：w(i)的第j个元素是M(0)的第j个元素和M(i)的第j个元素中较小的一个，w(i)的长度为M(0)的长度与M(i)的长度的最小值。其中|i|＜min(l_X，l_Y)。记AMF(i)为自匹配向量w(i)中各分量之和。

步骤32具体包含以下步骤：

步骤321：若片段X和Y的匹配向量M(0)的各分量之和小于h(0≤h≤l_X)，则判定X和Y之间只存在碱基替换，结束操作，否则进入步骤322；其中h为预先设定的值，一个优选的设定方法为：先设置一个检验水平α₁，取h为二项分布B(l_X，β+γ)的上α₁分位点，或者取h为参数为l_X(β+γ)的泊松分布的上α₁分位点，其中β为X的每一个碱基的碱基误读概率的平均值；

步骤322：设定自匹配函数值的阈值m(0≤m≤l_X)，顺次比较AMF(1)，AMF(-1)，...，AMF(B)，AMF(-B)，B为预设的可以被探测的插入或删失的最大长度值，其数值是根据测序序列所属物种不同个体的基因组之间存在的插入/删失长度的分布设定的；若存在某个整数j使得AMF(j)＜m，则判定X相对于Y只存在一个插入/删失，停止比较，进入步骤323，否则判定X与Y之间存在更复杂的差异，进入步骤324；所述阈值m的一个优选的设定方法为：先设置一个检验水平α₂，取m为二项分布B(l_X，β+γ)的上α₂分位点，或者取m为参数为l_X(β+γ)的泊松分布的上α₂分位点，其中β为X的每一个碱基的碱基误读概率的平均值；

图11示出了步骤322一个具体实例。如图11所示，其中测序序列片段对于参考基因组片段有一个长度为2的删失，此例中阈值h取为5，m取为4。

步骤323：如果j为正数，则判定X相对于Y存在一个长度为|j|的删失，否则判定X相对于Y存在一个长度为|j|的插入；利用匹配向量M(0)和M(j)探测插入或者删失的起始位点，结束操作；

步骤324：利用局部-全局比对算法寻找X与Y之间的细致差异；

上述步骤323所述的利用匹配向量M(0)和M(j)探测插入或者删失的起始位点：将M(0)的l_M(j)位前缀中每一个位置的数值减去M(j)相同位置的数值，其中l_M(j)为向量M(j)的长度；对于所得到的差向量中的每一个位置，计算从差向量的第一项至这一位置的所有项之和，将这些和中数值最小的一个所对应的差向量中的位置判定为插入或删失在碱基片段X上的起始位置。具体步骤如下：

步骤3231：将插入/删失的起始点P设置为0，设置变量MINI、DIAG、T，并将它们设置为0；

步骤3232：当T＜l_M(j)时，进入步骤3233，否则进入步骤3236；其中l_M(j)表示M(j)的长度；

步骤3233：将DIAG重新设置为DIAG+M(0)[T]-M(j)[T]，其中M(0)[T]，M(j)[T]分别为M(0)和M(j)的第T个元素；

步骤3234：若DIAG小于MINI，则将MINI重新设置为DIAG，并将P重新设置为T+1；

步骤3235：将T重新设置为T+1，返回步骤3232；

步骤3236：如果P为0，则判定插入/删失的起始点位于X[0]之前，否则判定插入/删失的起始点位于X的第P-1个碱基之后，即插入/删失起始位点之前有且仅有从X[0]至X[P-1]这P个碱基。

图12示出了本发明中探测插入/删失起始位点的一个具体实例。如图12所示，经过比较M(0)与M(2)，可以得出长度为2的删失的起始位点位于X[6]与X[7]之间；

上述步骤324中针对X与Y之间更加复杂的差异情况，利用局部-全局比对算法寻找X与Y的差异，即对两个序列进行比对。由于X与Y左端对齐但X的长度大于Y的长度，该算法在两条序列的左端采用全局比对算法，在两条序列右端采用局部比对算法，即在对得分矩阵(下文中提到的二维数组F)进行初始化时，将第一行初始化为插入/删失的罚分，而将第一列初始化为0。另外，该算法可以对不同种类的碱基替换采用不同的罚分。步骤324具体包含以下步骤：

步骤3241：建立二维数组F[l_Y+1][l_X+1]；将F[0][q]设置为δq，将F[p][0]设置为0，其中，l_X为X的长度，l_Y为Y的长度，q≤l_X，p≤l_Y，δ为插入/删失的罚分，δ＜0；

步骤3242：对p≥1且q≥1，按照以下方式设定二维数组的其它元素：

F[p][q]＝max{F[p-1][q-1]+θ(X[l_X-q]，Y[l_Y-p])，F[p-1][q]+δ，F[p][q-1]+δ}；

其中θ(X[l_X-q]，Y[l_Y-p])为X[l_X-q]与Y[l_Y-p]两种碱基匹配情况的罚分函数；可以根据质量值信息以及不同碱基组合设定不同罚分，如在对甲基化数据做比对时，可以将θ(T，C)设为接近0的数，而将θ(C，T)设为接近-1的数；

步骤3243：将起始位置p重新设置为l_Y，将起始位置q重新设置为l_X；

步骤3244：若q＞0，则进入步骤3245，否则结束操作；

步骤3245：若p＞0，且F[p][q]＝F[p-1][q-1]+θ(X[l_X-q]，Y[l_Y-p])，则进入步骤3246；若p＞0且F[p][q]＝F[p-1][q]+δ，则进入步骤3247；若p≥0，且F[p][q]＝F[p][q-1]+δ，则进入步骤3248；

步骤3246：如果X[l_X-q]和Y[l_Y-p]相同，则判定X[l_X-q]与Y[l_Y-p]匹配，否则判定两者不匹配；将p重新设置为p-1，将q重新设置为q-1，返回步骤3244；

步骤3247：判定Y[l_Y-p]这一碱基在X中被删失，将p重新设置为p-1，返回步骤3244；

步骤3248：判定X[l_X-q]相对于Y为一个插入碱基，将q重新设置为q-1，返回步骤3244。

作为一个比选方案，步骤32还可以按如下步骤实现：

步骤32B1：按照和步骤321相同的方法判定X与Y之间是否只存在碱基替换，若是则结束操作，否则进入步骤32B2；

步骤32B2：逐次对假设的事件“X相对于Y只存在一个长度为1的删失，X相对于Y只存在一个长度为1的插入，…，X相对于Y只存在一个长度为B的删失，X相对于Y只存在一个长度为B的插入”进行检验，B为预设的可以被探测的插入或删失的最大长度值；若其中的一个事件被接受，则结束操作，若所有的事件都被拒绝，则进入步骤32B3；具体地，在对其中的一个事件进行检验时，可以先采用和步骤323相同的方法进行操作，找到插入/删失的起始位点，得到这一事件成立时X与Y之间的所有碱基替换位点；若X[i]为一个碱基替换位点，则计算软分值

其中β_i为根据质量值换算的X[i]的误读概率；对所有软分值求和，设定检验水平α₃，若软分值之和大于二项分布B(l_X，γ)或者参数为l_Xγ的泊松分布的上α₃分位点，则拒绝该事件，否则接受该事件；在进行假设检验时，也可以用碱基替换的个数代替软分值；

图13示出了上述步骤的一个实施例，该实施例使用碱基替换个数代替软分值。如图13所示，测序序列片段相对于参考基因组片段存在一个长度为2的删失，对应于“长度为1的删失”、“长度为1的插入”这两个事件的碱基替换总数较大，两个事件被拒绝，而对应于事件“长度为2的删失”的碱基替换总数为0，该事件被接受；

步骤32B3：按照和步骤324相同的方法，采用局部-全局比对算法寻找X与Y的细致差异；

以上所述为寻找完全匹配前缀子序列右侧的片段与参考基因组的差异的方法，对于完全匹配前缀子序列左侧的片段和参考基因组上对应的片段，可以使用类似的方法进行操作，也可以将两者完全翻转，采用和步骤32完全一样的方法寻找差异。

在分别找到完全匹配前缀子序列两侧的片段与参考基因组的差异后，将其结合得到整个测序序列与参考基因组之间存在的差异。

步骤33具体操作包括以下步骤：

步骤331：取定测序序列所属生物物种基因组的多态率γ，根据测序序列上完全匹配前缀子序列以外的所有碱基与参考基因组上相应碱基的匹配状态计算“软分值”，具体地，若测序序列中的第i个碱基为碱基替换，则对这一碱基计算“软分值”

其中β_i为根据质量值换算的这一碱基的误读概率，其数值为

QA_i为测序序列上第i个碱基的质量值；也可以在换算之前对质量值进行修正以得到更准确的误读概率；

步骤332：将步骤331中计算得到的所有比例值

相加，记为sum；如果l≥10且γ≤0.1，则取p_value为sum关于参数为lγ的泊松分布的上侧分位数，否则取p_value为sum关于二项分布B(l，γ)的上侧分位数；或者取p_value为关于标准正态分布N(0，1)的上侧分位数；同时取power为sum关于二项分布B(l，1-ε)的上侧分位数；如果p_value≥α，则判定接受测序序列上完全匹配前缀子序列以外的所有碱基与参考基因组上相应碱基的匹配状态，否则拒绝测序序列上完全匹配前缀子序列以外的所有碱基与参考基因组上相应碱基的匹配状态，其中，l为测序序列长度减去完全匹配前缀子序列长度所得的差，α为步骤2中取定的检验水平。

由于测序序列可以来自参考基因组的正向序列与反向互补序列中的任何一条，如果测序序列映射失败，还应该按照相同的方法将测序序列的反向互补序列映射至参考基因组。所述反向互补序列为：将正向序列的每一个碱基替换为与其互补配对的碱基，再进行左右翻转后得到的序列。

在对测序序列完成一轮映射后，根据映射速度、被成功映射的测序序列所占的比例，以及被成功映射的测序序列的灵敏度和特异度分布情况调整参数，利用和步骤3相同的方法对映射失败的测序序列，或者对灵敏度、特异度较差的测序序列重新进行映射。

步骤5为输出步骤，输出每一个测序序列的映射信息，所有输出的信息写入具有标准格式的文件，对于每个测序序列，输出的映射信息包含以下内容：

a)该测序序列在数据集中的编号，表示为整数或字符串；

b)表示测序序列是否被成功映射，表示为整数或字符串，映射失败的信息由步骤31得到，映射成功的信息由步骤33得到；

c)测序序列所被映射到的参考基因组名称，表示为字符串；

d)测序序列在参考基因组上所处的位置，由步骤32得到，表示为整数；

e)测序序列与参考基因组之间存在的碱基替换、插入、删失，由步骤32得到，表示为一个字符串；

f)测序序列的碱基序列，表示为字符串；

g)测序序列每个碱基的质量值信息，表示为字符串；

h)完全匹配前缀子序列长度k；

i)映射的灵敏度Ψ(L，k，1-β-γ)×(1-p_value)，由步骤33得到；

j)映射的特异度1-min{NΨ(L，k，ε)×(1-power)，1}，由步骤33得到；

与步骤2相比，此处灵敏度和特异度表达式中的k为步骤31结束时所找到的完全匹配前缀子序列的长度，且在灵敏度的表达式中以步骤332中得到的p_value代替检验水平α，在特异度的表达式中将NΨ(L，k，ε)替换为NΨ(L，k，ε)×(1-power)。所述标准格式为：对于映射成功的测序序列，依次按照a)至j)的顺序输出映射信息；对于映射失败的测序序列，只输出a)，b)，f)，g)。

对于甲基化测序序列的映射，在与上述方法相同的框架下，需要做相应的变化。

1.在预处理部分，参考基因组压缩结构的生成方式不变。但是生成地址索引结构的步骤121中，可以采用与非甲基化情况不同的编码方式：

编码A：仍然用两个二进制码代表一个碱基，但不区分C和T。比如，用00、01代表A、G，用10代表C或T；

编码B：用Huffman码来表示，比如，用00、01代表A、G，用1代表C或T。这时，碱基的编码字长是不等长的，为1或2。

在生成地址索引结构时，对于非甲基化情况，考虑全部n位子序列的碱基值，它们以2n位二进制码表示。在甲基化情况中，仍然以2n位二进制码表示从基因组的每个位置开始的子序列。如果采用编码A，它们仍然对应于n位子序列，如果采用编码B，它们所对应的子序列长度可以不同，其长度可以为闭区间[n，2n]中的任意整数，这依赖于C/T的含量。其它步骤如排序、存储地址值的操作不变。但在甲基化情况中，需对参考基因组的正向序列和反向互补序列都生成地址索引结构。

2.在参数设计部分，按照A、G、C/T三个字符和它们的频率计算特异度。

3.在子序列定位阶段，步骤311要做相应的变动，仍然以2n位二进制码表示从测序序列的每个位置开始的子序列，但编码方式要与生成地址索引结构时的相同。如果采用编码A，则该2n位二进制码所表示的子序列长度为n，如果采用编码B，则该2n为二进制码所表示的子序列长度可以为闭区间[n，2n]中的任意整数。

同时，在根据地址值从参考基因组压缩结构中获取碱基值时，要将从参考基因组压缩结构中取出的参考基因组上的子序列所对应的二进制码进行转化，得到按照生成地址索引结构时所使用的编码方式编码所得的该子序列的碱基值。

4.在延拓部分，先用自匹配方法将测序序列与参考基因组之间的差异分为三种情况：(1)只有碱基替换；(2)在无错误匹配子序列的一侧只有一个插入/删失；(3)其它更复杂的情况；在操作时不区分C和T，即将C和T当作同一种碱基。对(2)，进一步用自匹配方法的确定插入/删失的类型、位置和长度，此时区分C和T，将它们当作不同碱基。对(3)，用比对算法确定差异的细节，此时区分C和T并采用适用于甲基化数据的分值系统。特别地，采用局部-全局比对算法确定测序序列延拓部分(全局)和参考基因组相应DNA片段(局部)之间的差异细节。

5.如果测序序列映射至参考基因组的正向序列失败，还需将测序序列映射至参考基因组的反向互补序列，此时使用反向互补序列的地址索引结构。所述反向互补序列的获取方法为：对参考基因组正向序列的每一个碱基替换为与其互补配对的碱基，将所得到的序列左右翻转后得到的序列。

在系统实现上，采用软件和硬件的组合形式。这一系统主要分为预处理、子序列定位和基于自匹配函数的延拓三个环节。预处理部分对参考基因组的给定版本只需做一次，结果保留在存储系统中，可以对不同的测序数据集反复使用。采用快速排序或者堆排序算法将2n位二进制码排序。

计算量最大的子序列定位和延拓环节可由以下方法实现。

1.基于C++等编程语言的软件实现。

2.在映射部分采用多个CPU的并行运算。

3.子序列定位环节需要用到地址索引结构、参考基因组压缩结构、等分点索引结构，而在具体的操作中，只是从这些结构中读取数据，并不改写数据，这三个公共存储空间结构可以用高性能的存储硬件实现。

4.当参考基因组上的子序列被取出后，基于自匹配函数的延拓阶段要将其与测序序列进行比对，而比对过程中不需要从公共存储空间中读取数据，而且自匹配方法和局部-全局比对算法只是在两个字符串间进行比较，程序非常简单，在映射的延拓部分可采用GPU技术实现高性能并行。

5.类似于Smith-Waterman算法的FPGA生物计算系统(Timelogic，DeCypher)，在核心算法部分采用FPGA或可编程的DSP进行固化。

上述方法可以以软件和硬件的组合形式实现。这包括，在映射部分采用多个CPU的并行运算实现，在映射的延拓部分采用GPU技术实现，以及在核心算法部分采用DSP和FPGA技术。

本发明还公开了一种测序序列映射系统，其对于所获取参考基因组和至少一个测序序列进行操作；其中，参考基因组为已完成测序的基因组序列；该系统包括：

参考基因组预处理模块，其用于对所述参考基因组进行预处理，以生成参考基因组压缩结构、地址索引结构和等分点索引结构；所述参考基因组压缩结构以压缩形式存储整个参考基因组，所述地址索引结构按照一定次序存储所述参考基因组中所有子序列的地址值，所述等分点索引结构用于存储一部分地址值在地址索引结构中所处的位置，用于加速实现测序序列的初步定位；

参数设计模块，其用于基于参考基因组的特征、测序序列的整体信息、测序序列所属物种的多态发生情况，根据概率计算，设计映射算法的参数，以达到或折中对灵敏度和特异度、映射速度的要求；

映射模块，其用于根据经过预处理后得到的上述三种结构，通过子序列定位、基于自匹配函数的延拓、定量分析步骤将每一个测序序列映射至所述参考基因组上；

结果输出模块，其用于输出每一个测序序列的映射信息，将所有输出的信息写入具有标准格式的文件。

上述各个模块均可以通过软件、硬件或者软硬件结合的方式实现。

本发明提出的上述方法和系统，具有以下优点：

1.在预处理阶段构造地址索引结构时，存储排序后的所有参考基因组中n位子序列的地址值，而不是直接存储参考基因组中n位子序列，这样可以减小存储空间；

2.利用排序后的碱基值与均匀分布的差异确定编码规则，可以采用等分点索引结构加速二分法查找；

3.在预处理阶段可以根据参考基因组的长度选用不同长度的无符号整数实现地址索引结构的灵活存储；

4.在定位一个子序列时利用等分点索引结构初步缩小查找范围，二分法查找过程中不再考察已知的上下界的公共前缀长度，加快查找速度；

5.定位一个子序列时可以找到不同长度的完全匹配前缀子序列，相比基于哈希表的查找法更加灵活；

6.延拓阶段自匹配函数法的时间复杂度为线性，若自匹配法失败，局部-全局比对法可以找到更准确的比对情况；

7.子序列定位、基于自匹配函数的延拓阶段均有相对应的定量分析，对结果进行评估，保证映射的敏感性(测序序列可以被映射到参考基因组上)和特异性(测序序列不会被映射至错误的地方)，定量分析可以指导参数设计，使得操作更加灵活；

8.预处理和映射速度快于现有工具，映射率普遍比现有工具高，适用于普通的工作站。

9.对于甲基化数据，通过Huffman编码或其它编码来代表A、G、C、T。在映射甲基化数据时，对自匹配算法、局部-全局比对算法等处对碱基种类之间的罚分做相应调整。

10.方法可以根据不同的需求和所具备的硬件设备，以软件和硬件的组合形式实现系统。这包括，在映射部分采用多个CPU的并行运算实现，在映射的延拓部分采用GPU技术实现，以及在核心算法部分采用DSP和FPGA技术。

11.从表现上看，本方法的预处理和映射速度快于现有工具，映射率普遍比现有工具高，适用于普通的工作站。

本发明提出的上述方案适用于生物分析和个体化医疗中DNA遗传信息和RNA功能信息的定性或定量测量，包括DNA多态、甲基化、信使RNA表达谱、非编码RNA丰度、CHIP-seq、可变剪切信息的定性或定量测量。

以上所述的具体实施例，对本发明的目的、技术方案和有益效果进行了进一步详细说明，应理解的是，以上所述仅为本发明的具体实施例而已，并不用于限制本发明，凡在本发明的精神和原则之内，所做的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims (18)

1.一种测序序列映射方法，对于所获取参考基因组和至少一个测序序列进行操作；其中，参考基因组为已完成测序的基因组序列；操作包括以下步骤：

步骤1、对所述参考基因组进行预处理，以生成参考基因组压缩结构、地址索引结构和等分点索引结构；所述参考基因组压缩结构以压缩形式存储整个参考基因组，所述地址索引结构按照一定次序存储所述参考基因组中所有子序列的地址值，所述等分点索引结构用于存储一部分地址值在地址索引结构中所处的位置，用于加速实现测序序列的初步定位；

步骤2、基于参考基因组的特征、测序序列的整体信息、测序序列所属物种的多态发生情况，根据概率计算，设计映射算法的参数，以达到或折中对灵敏度、特异度、映射速度的要求；

步骤3、根据经过预处理后得到的上述三种结构，通过子序列定位、基于自匹配函数的延拓、定量分析步骤将每一个测序序列映射至所述参考基因组上；

步骤4、输出每个测序序列的映射信息。

2.如权利要求1所述的测序序列映射方法，其特征在于，步骤1中在生成参考基因组压缩结构、地址索引结构时，将表示参考基因组的碱基字符集合{A，C，G，T}按照一定的二进制编码规则映射至二进制表示集合{00，01，10，11}中的二进制码；其中，对于非甲基化数据的映射，所述编码规则指：

所述碱基字符集合{A，C，G，T}与二进制表示集合{00，01，10，11}中元素之间一一对应的映射方式。

3.如权利要求2所述的测序序列映射方法，其特征在于，所述编码规则通过如下步骤确定：

步骤11、在碱基字符集合{A，C，G，T}与二进制表示集合{00，01，10，11}之间的每一种映射方式下，将参考基因组中每一个n位子序列映射成二进制数，该二进制数即为该n位子序列的碱基值，n为预先设定的子序列长度；

步骤12、将所述每一种映射方式下的n位子序列的碱基值序列，按顺序排序，将最接近均匀分布的碱基值序列对应的映射方式作为所述编码规则。

4.如权利要求1所述的测序序列映射方法，其特征在于，步骤1中参考基因组压缩结构如下生成：

依照二进制编码规则，利用二进制码将参考基因组序列按照从左到右，或者从右到左的方向逐位存入到一个字节类型的向量中，每个字节存储四个碱基。

5.如权利要求1所述的测序序列映射方法，其特征在于，步骤1中地址索引结构如下生成：

对于参考基因组中的每一个n位子序列按照二进制编码规则获得其碱基值，并按顺序对每一个n位子序列的碱基值进行排序，排序后的每一个碱基值所对应的n位子序列在参考基因组中的地址值形成的序列为所述地址索引结构；其中，n为预设的参考基因组中子序列的长度。

6.如权利要求1所述的测序序列映射方法，其特征在于，步骤1中等分点索引结构如下生成：

对于区间[0，4ⁿ-1]中的等分点z_i＝i×2^2n-c，i＝0，1，2...，2^c-1，比较所述等分点z_i与地址索引结构中每个地址值对应的n位子序列的碱基值，找到一个地址值或者两个相邻地址值，使得该等分点z_i等于该一个地址值所对应的n位子序列的碱基值或者位于该两个相邻地址值所对应的n位子序列的碱基值之间，并将该一个地址值在地址索引结构中的索引值或该两个相邻地址值在地址索引结构中的较小索引值存储为该等分点索引结构的第i个值，最终得到等分点索引结构。

其中c为预先设定的整数值，其取值范围1≤c≤2n，n为预设的参考基因组中子序列的长度。

7.如权利要求1所述的测序序列映射方法，其特征在于，步骤2中所述参数包括完全匹配前缀子序列长度下界k₀和检验水平α，步骤1中构造地址索引结构所用的参考基因组中子序列的长度n大于k₀，k₀和检验水平α的选择满足对映射的灵敏度和特异度的要求，参数选择依赖于参考基因组的长度、多态率和碱基频率，测序序列长度和质量值，通过概率计算进行选择。

8.如权利要求1所述的测序序列映射方法，其特征在于，步骤3中将测序序列映射至参考基因组具体包括：

步骤31、对于一个测序序列，利用步骤1所生成的三种结构、或者只利用参考基因组压缩结构和地址索引结构两种结构，在测序序列的每一个n位子序列中寻找完全匹配前缀子序列，直到找到一个长度大于完全匹配前缀子序列长度下界k₀的完全匹配前缀子序列，进入步骤32，否则判定映射失败；其中所述完全匹配前缀子序列为所述n位子序列中与所述参考基因组中的相应子序列完全匹配的前缀序列；

步骤32、对于寻找到的完全匹配前缀子序列，利用自匹配函数法将测序序列中所述完全匹配前缀子序列左右两侧和在参考基因组中对应于所述完全匹配前缀子序列的位置左右两侧的碱基片段之间的差异分为三类：(1)只存在碱基替换，(2)完全匹配前缀子序列一侧只存在一个插入/删失，(3)其他更复杂的情况；对(2)，进一步用自匹配方法确定插入/删失的类型、位置和长度；对(3)，用局部-全局比对算法确定差异的细节；

步骤33、对于测序序列上完全匹配前缀子序列以外的部分，根据所属生物物种基因组的多态率和每个碱基辨识的质量值，利用二项分布、泊松分布或者正态分布将所有碱基与参考基因组上相应碱基的匹配状态进行评估，以判断是否接受这些碱基的匹配状态，如果接受，则判定映射成功，如果拒绝，则返回步骤31选择新的完全匹配前缀子序列。

9.如权利要求8所述的测序序列映射方法，其特征在于，步骤31具体包括：

步骤311、按照一定的编码规则得到测序序列上的n位子序列的碱基值；

步骤312、通过计算得到区间[0，4ⁿ-1]等分为2^c份后的两个相邻的等分点，使得测序序列上的n位子序列按照一定的编码规则编码后得到的碱基值大于或等于较小的等分点，同时严格小于较大的等分点；

步骤313、根据所述两个相邻的等分点在等分点索引结构中获取它们对应的地址索引值，以其中较小的地址索引值作为映射范围的下限起始位置，以其中较大的地址索引值加1后所得的值作为映射范围的上限起始位置；将完全匹配前缀子序列的长度k设置为0；

步骤314、获取参考基因组上从下限起始位置所对应的地址值开始的n位子序列的n-k位后缀和参考基因组上从上限起始位置所对应的地址值开始的n位子序列的n-k位后缀的最大公共前缀长度f，并更新完全匹配前缀子序列的长度k，使得k＝k+f；

步骤315、比较参考基因组上从中间起始位置所对应的地址值开始的n位子序列的n-k位后缀的碱基值与测序序列上的n位子序列的n-k位后缀的碱基值，所述中间起始位置为所述下限起始位置和上限起始位置的中间位置；若测序序列上的n位子序列的n-k位后缀的碱基值大于或等于参考基因组上从所述中间起始位置所对应的地址值开始的n位子序列的n-k位后缀的碱基值，则下限起始位置更新为所述中间起始位置，否则上限起始位置更新为所述中间起始位置；

步骤316、重复执行步骤314至315，直到上限起始位置和下限起始位置之差小于或等于1，然后比较参考基因组上从下限起始位置所对应的地址值开始的n位子序列的n-k位后缀与测序序列上的n位子序列的n-k位后缀，获得两者最大公共前缀，将该前缀的长度与k相加得到第一长度；并比较参考基因组上从上限起始位置所对应的地址值开始的n位子序列的n-k位后缀与测序序列上的n位子序列的n-k位后缀，获得两者最大公共前缀，将该前缀的长度与k相加得到第二长度；

步骤317、若所述第一长度和第二长度中的较大值大于k₀，则所述较大值对应的测序序列上n位子序列的前缀序列即为得到的完全匹配前缀子序列；

步骤318、若参考基因组上可以被完全匹配前缀子序列定位的位置数目超过了预定值G，则转步骤311进行下一个测序序列上的n位子序列的定位。

10.如权利要求8所述的测序序列映射方法，其特征在于，步骤32中获取测序序列中所述完全匹配前缀子序列和在参考基因组中对应于所述完全匹配前缀子序列的位置的左侧碱基片段与右侧碱基片段之间的差异过程类似，其中获取右侧碱基片段之间的差异过程具体包括：

步骤321、获取所述测序序列中位于完全匹配前缀子序列右侧的片段X，并获取参考基因组中对应于所述完全匹配前缀子序列的位置右侧的片段Y，使得Y的长度l_Y大于X的长度l_X；

步骤322、若表示片段X和Y的匹配程度的匹配向量M(0)中所有元素之和小于预定值h，则判定片段X和Y之间的差异只包含碱基替换，结束操作，否则转步骤323；其中M(0)中的每个元素表示片段X和Y相应位置处的碱基的匹配程度，0≤h≤l_X，l_X为片段X的长度；

步骤323、顺次比较自匹配函数值的预设阈值m与AMF(1)，AMF(-1)，...，AMF(B)，AMF(-B)，B为预设可以被探测插入或删失的最大长度值；若存在某个整数j使得AMF(j)＜m，则进入步骤324，否则转步骤325；

步骤324、若j为正数，则判定片段X和Y之间的差异为片段X相对于片段Y存在一个长度为|j|的删失，否则判定片段X和Y之间的差异为片段X相对于Y存在一个长度为|j|的插入；探测插入或删失的起始位点，结束操作；

步骤325、利用局部-全局比对算法寻找X与Y之间的其它差异；

其中，M(0)＝V(X，Y)；若i＞0，则M(i)＝V(X，Y{i})，Y{i}为Y的l_Y-i位后缀；若i＜0，则M(i)＝V(X{-i}，Y)，X{-i}为X的l_X-|i|位后缀；V(X，Y)的长度为l_X和l_Y的最小值，V(X，Y)的元素具有两种定义：(一)若片段X和Y的第i个碱基相同，则V(X，Y)的第i个元素为β_i，否则为1-β_i；其中β_i是由片段X的第i个碱基的质量值转换得到的碱基误读概率；(二)若片段X和Y的第i个碱基相同，则V(X，Y)的第i个元素为0，否则为1；

AMF(i)为自匹配向量w(i)中所有项之和；而自匹配向量w(i)如下定义：w(i)中的元素值是M(0)和M(i)中相应位置处元素值的最小值，w(i)的长度为M(0)的长度与M(i)的长度的最小值，i为绝对值小于片段M(0)和M(i)长度的最小值的整数。

11.如权利要求10所述的测序序列映射方法，其特征在于，步骤324具体包括：

对于步骤323中得到的M(j)，将M(0)的l_M(j)位前缀中每一个位置的数值减去M(j)相同位置的数值，其中l_M(j)为向量M(j)的长度；对于所得到的差向量中的每一个位置，计算从差向量的第一项至这一位置的所有项之和，将这些和中数值最小的一个所对应的差向量中的位置判定为插入或删失在碱基片段X上的起始位置。

12.如权利要求10所述的测序序列映射方法，其特征在于，步骤325具体包括：

步骤3251：建立二维数组F[l_Y+1][l_X+1]；将F[0][q]设置为δq，将F[p][0]设置为0，其中，l_X为X的长度，l_Y为Y的长度，q≤l_X，p≤l_Y，δ为预设的插入/删失的罚分，δ＜0；

步骤3252：对p≥1且q≥1，按照以下方式设定二维数组的其它元素：

F[p][q]＝max{F[p-1][q-1]+θ(X[l_X-q]，Y[l_Y-p])，F[p-1][q]+δ，F[p][q-1]+δ}；

其中θ(X[l_X-q]，Y[l_Y-p])为X[l_X-q]与Y[l_Y-p]两种碱基匹配情况的罚分函数；

步骤3253：将p重新设置为l_Y，将q重新设置为l_X；

步骤3254：若q＞0，则进入步骤3255，否则结束操作；

步骤3255：若p＞0，且F[p][q]＝F[p-1][q-1]+θ(X[l_X-q]，Y[l_Y-p])，则进入步骤3256；若p＞0且F[p][q]＝F[p-1][q]+δ，则进入步骤3257；若p≥0，且F[p][q]＝F[p][q-1]+δ，则进入步骤3258；

步骤3256：如果X[l_X-q]和Y[l_Y-p]相同，则判定X[l_X-q]与Y[l_Y-p]匹配，否则判定两者不匹配；将p重新设置为p-1，将q重新设置为q-1，返回步骤3254；

步骤3257：判定Y[l_Y-p]这一碱基在X中被删失，将p重新设置为p-1，返回步骤3254；

步骤3258：判定X[l_X-q]相对于Y为一个插入碱基，将q重新设置为q-1，返回步骤3254。

13.如权利要求8所述的测序序列映射方法，其特征在于，步骤32通过以下步骤实现：

步骤321、判断片段X和Y之间是否只存在碱基替换，若是则结束操作，否则转步骤322；

步骤322、对长度不超过B的每一种插入或删失情况，探测其起始位点，得到比对结果，并评估比对结果，若该比对结果被肯定，则结束操作，否则对下一种插入或删失情况进行比对；其中B为预设可以被探测插入或删失的最大长度值；

步骤323、若所有插入或删失情况的比对结果均被否定，则使用局部-全局比对法对X与Y进行比对。

14.如权利要求8所述的测序序列映射方法，其特征在于，步骤33具体包括：

步骤331：获取测序序列所属生物物种基因组的多态率γ，若测序序列的第i个碱基与参考基因组上相应的碱基不匹配，则对测序序列的第i个碱基计算比例值

其中β_i为由该碱基的质量值转化而得的碱基误读概率；

步骤332：取sum为所有比例值之和，计算sum关于二项分布B(l，γ)，或参数为lγ的泊松分布的上侧分位数；或者计算

关于标准正态分布N(0，1)的上侧分位数；如果上侧分位数不小于α，则接受测序序列上位于当前完全匹配前缀子序列两侧的所有碱基与参考基因组上相应碱基的匹配状态，否则拒绝测序序列上位于当前完全匹配前缀子序列两侧的所有碱基与参考基因组上相应碱基的匹配状态；其中，l为测序序列长度减去完全匹配前缀子序列长度所得的差，检验水平α在步骤2中设定。

15.如权利要求1所述的测序序列映射方法，其特征在于，步骤4之前还包括：对未能映射成功的测序序列，或者对灵敏度、特异度较差的测序序列，通过调整映射参数重新进行映射；且步骤4所输出的映射信息具体包括：测序序列在数据集中的编号、测序序列是否被成功映射、测序序列被映射到的参考基因组名称、测序序列在参考基因组上所处的位置、测序序列与参考基因组之间存在的差异、测序序列的碱基序列、测序序列每个碱基的质量值信息、完全匹配前缀子序列长度、映射的灵敏度、映射的特异度。

16.如权利要求1所述的测序序列映射方法，其特征在于，对于甲基化测序序列数据的映射，生成参考基因组压缩结构时所用的编码方式与非甲基化情况相同，在生成地址索引结构时，采用如下的编码方式：

编码A：用两个二进制码代表一个碱基，但不区分C和T；

编码B：用Huffman码来表示，C或T用一位相同的二进制码表示，A、G分别用两位二进制码表示，其中A和G所用的编码的第二位不同但第一位相同，且与表示C或T的一位二进制码不同；

对于甲基化测序序列的映射，需对参考基因组的正向序列和反向互补序列都生成地址索引结构；且在步骤3中映射甲基化测序序列时，将甲基化测序序列映射至参考基因组的正向序列和反向互补序列上。

17.如权利要求1所述的测序序列映射方法，其特征在于，该方法适用于生物分析和个体化医疗中DNA遗传信息和RNA功能信息的定性或定量测量，包括DNA多态、甲基化、信使RNA表达谱、非编码RNA丰度、CHIP-seq、可变剪切信息的定性或定量测量。

18.一种测序序列映射系统，该系统对于所获取参考基因组和至少一个测序序列进行操作；其中，参考基因组为已完成测序的基因组序列；该系统包括：

参考基因组预处理模块，其用于对所述参考基因组进行预处理，以生成参考基因组压缩结构、地址索引结构和等分点索引结构；所述参考基因组压缩结构以压缩形式存储整个参考基因组，所述地址索引结构按照一定次序存储所述参考基因组中所有子序列的地址值，所述等分点索引结构用于存储一部分地址值在地址索引结构中所处的位置，用于加速实现测序序列的初步定位；

参数设计模块，其用于基于参考基因组的特征、测序序列的整体信息、测序序列所属物种的多态发生情况，根据概率计算，设计映射算法的参数，以达到或折中对灵敏度和特异度、映射速度的要求；

映射模块，其用于根据经过预处理后得到的上述三种结构，通过子序列定位、基于自匹配函数的延拓、定量分析步骤将每一个测序序列映射至所述参考基因组上；

结果输出模块，其用于输出每一个测序序列的映射信息。

Images