CN103546162B

CN103546162B - 基于非连续上下文建模和最大熵原则的基因压缩方法

Info

Publication number: CN103546162B
Application number: CN201310433453.9A
Authority: CN
Inventors: 熊红凯; 李平好
Original assignee: Shanghai Jiaotong University
Current assignee: Shanghai Jiaotong University
Priority date: 2013-09-22
Filing date: 2013-09-22
Publication date: 2016-08-17
Anticipated expiration: 2033-09-22
Also published as: CN103546162A

Abstract

本发明提供一种基于非连续上下文建模和最大熵原则的基因压缩方法，包括：在第一阶段，同时考虑了有参考序列和无参考序列两种实际情况，用字典方法表示基因序列内或间的重复序列，提高压缩效率；第二阶段，所述统计编码器由预测器和算术编码器构成，预测器对于非重复序列，采用在传统的连续型上下文模型的基础上引入非连续比特上下文组合式的建模方案，并利用基于最大熵原则产生的Logistic回归公式合并各个模型单独产生的预测概率，从而得到最终更准确的预测概率值送入算术编码器。该方法能显著提高压缩效率，实现高效存储。

Description

基于非连续上下文建模和最大熵原则的基因压缩方法

技术领域

本发明涉及一种面向超大规模基因序列的信息压缩方法，具体涉及的是一种将非连续上下文建模思想与传统的连续上下文建模思想相结合，并利用最大熵原则得到最终预测概率的统计式压缩方法。

背景技术

DNA是生物生存、延续和发展的重要物质基础，具有重大的科学价值和社会价值。目前，DNA的研究广泛应用于生物学、医学、遗传科学等诸多重要领域，如通过收集和保存DNA信息以保护濒临灭绝的生物物种、基于人类基因序列的信息预测以及找到基因变异规律以治疗癌症肿瘤等。为这些学科研究提供基础实验数据的各种DNA序列测定工程已成为各国重点发展的研究项目。随着这些测序项目的展开，每天都有海量的DNA序列数据产生，相关数据量呈指数方式增长，生物信息数据这种急速的积累增长在人类的科学研究历史中是空前的。存储和使用这些数据的成本已越来越面临着无法承担的规模，如何在有限的存储资源内有效储存急剧膨胀的DNA序列数据成为了计算机专家和生物学家面临的新课题，也是国内外诸多重大计划所面临的前进障碍。因此，采用更有效的压缩编码方式，用较小的存储空间存放较大的基因信息序列是必然的选择。

经过对现有技术的文献检索发现，2000年T Matsumoto和K Sadakane在GenomeInformatics上的“Biological sequence compression algorithms”里提出了CTW+LZ方法，将上下文树加权(Context tree weighting,CTW)方法和LZ压缩方法相结合，使用多个编码模型对DNA序列的不同片段上进行压缩；2002年X Chen和M Li等在Bioinformatics上的“DNACompress:fast and effective DNA sequence compression”里提出了DNACompress压缩方法，使用Pattern Hunter工具搜索DNA序列的重复与近似重复片段，提高了方法的整体速度；2005年，G Korodi和I Tabus在ACM Transactions on InformationSystems期刊上的“An Efficient Normalized Maximum Likelihood Algorithm for DNASequence Compression”中提出了GeNML方法，对具有不同数据特点的DNA片段使用不同的编码策略和概率模型进行压缩；同时还有针对自身片段重复的COMRAD(Shanika Kuruppu等发表于2012年的IEEE/ACM Transactions on Computational Biology andBioinformatics期刊的“Iterative Dictionary Construction for Compression ofLarge DNA Datasets”上)、DSRC(Sebastian Deorowicz等发表于2011年的Bioinformatics期刊上的“Compression of DNA sequence reads in FASTQ format”上)，Recoil(Vladimir Yanovsky等发表于2011年的Algorithms for Molecular Biology期刊的“ReCoil-an algorithm for compression of extremely large datasets of DNA data”上)等方法；另外利用参考序列结构组成相似的GRS(Congmao Wang等发表于2011年的Nucleic Acids Research期刊的“A novel compression tool for efficient storageof genome resequencing data”上)、GReEn(Armando J.Pinho等发表于2012年的NucleicAcids Research期刊的“GReEn:a tool for efficient compression of genomeresequencing data”上)、SLIMGENE(Christos Kozanitis等发表于2011年的Journal ofComputational Biology期刊的“Compressing genomic sequence fragments usingSlimGene”上)等方法也逐渐兴起。这些方法将序列视作由特殊字符构成的长字符串，从数据的构成特点出发进行整体处理，有效的提高了压缩效率和压缩时间。但总体而言DNA序列压缩技术仍处于起步阶段，组成基因片段的生物信息特征以及片段内部的细节重复特性并没有被充分发挥利用。在对序列进行统计建模求取概率进而使用算术编码的历程中，几乎所有的压缩方法都是逐一观测符号序列并顺序建模进行预测，然后通过贝叶斯平均法得到最终预测概率，这种过分单纯的预测模式对于获取以非传统规律排列的基因序列形式内的相互关联是十分不利的。

发明内容

本发明针对现有技术的不足，提出了一种更加有效的基于非连续上下文建模的基因压缩方法，该方法能显著提高压缩效率，实现高效存储。

本发明是通过以下技术方案实现的：

本发明所述的基于非连续上下文建模和最大熵原则的基因压缩方法，该方法包括两大阶段：

在第一阶段，对于无参考序列的压缩，通过在已压缩的子序列中向前搜索找到基因序列内部近距离和远距离的重复(包括直接重复和互补回文)，以重复发生的属性、位置和长度高效表示基因序列的重复子序列；对于有参考序列的压缩，通过在参考序列中向前向后双向搜索找到目标序列与参考序列中的重复部分，并通过重复发生的位置和长度高效表示重复子序列达到压缩目的。两种情况下剩余的非重复序列均将送入第二阶段的以比特位(Bit)为单位的统计编码器。

第二阶段，统计编码器由预测器和算术编码器构成。预测器对于待编码比特的预测准确性将直接影响到算术编码器的编码效率，预测的准确概率越高，算术编码器所用的编码比特越少，即压缩效率越高。为提高预测概率，与已有的编码不同，除了传统的连续上下文模型(模型上下文选择分别为当前字节的前1,2,4,6,8,10,11,12,13,14,16个字节)、动态马尔科夫模型和基因压缩中常用的重复匹配模型外，本发明根据实验测试又选择了13个非连续上下文模型，即选取在位置上相互互补且对测试实验产生很好压缩效果的11个利用非连续的比特位组合作为上下文来进行预测的模型(考虑到模型过多会增加算法复杂度，与传统的连续上下文模型相对应地，本发明确定非连续比特位组合的上下文模型个数也为11个，然后按照位置上相互互补的原则选取了二十对11个模型的组合，最后选取了一组对测试实验产生最好压缩率的11个上下文组合)、1个根据最近距离的重复匹配(超过3个以上比特位相同)的上下文进行预测的距离模型和1个根据前两个比特位的上下文进行预测的间接模型，以达到优化整体压缩性能的目的。

具体来说，上述模型选定以后，第二阶段的编码将经过以下过程：

①、首先根据当前需要编码(预测)的比特位y_j的上下文(对于第i个模型假设有m_ij个0比特，n_ij个1比特)计算上述所描述的每个模型对当前比特位为1的独立预测概率并求出每个概率的伸展形式

②、基于最大熵原则得到使得熵最大的合成概率为Logistic回归模型：其中M_c为模型的总个数，M_c＝11+13＝24，w_i是第i个模型此时的权值大小，w_i越大，表明该模型对最终的预测概率影响力越大，e是自然对数的底数；

③、将合成概率送入算术编码器输出；

④、根据牛顿法，以使编码消耗C＝-ln(1-|E|)最小(即预测误差E＝y_j-P(y|t)最小)为目标求出更新后的w_i，为：

w_{i} = w_{i} - \frac{\frac{\partial C}{\partial_{w_{i}}}}{\frac{\partial^{2} C}{\partial_{w_{i}^{2}}}} = w_{i} - \frac{1}{t_{i} (1 - y_{j} - P (y_{j} | t_{1}^{M_{c}}))};

⑤、使用新的权值w_i返回第①步，开始预测y_j的下个比特位y_j+1。

解压缩时，首先按照上述同样的①-⑤步恢复出原基因序列中非重复序列部分，然后根据压缩时第一阶段所产生的码字恢复重复序列，最后可得到完整无损的原基因序列。

与现有技术相比，本发明的有益效果是：

本发明所提出的基于非连续上下文建模和最大熵原则的基因压缩方法，提高了基因压缩的效率和完整性。本发明将非连续上下文建模思想与传统的连续上下文建模思想相结合，提高了预测模型的多样性与全面性，更适应于以非传统规律排列的基因序列；利用最大熵原则得到Logistic回归模型进行模型合成，克服了以往贝叶斯平均法过度依赖单一模型而忽略其他模型有效性的缺点，并使得合成概率更加准确；考虑了有参考序列和无参考序列两种情况，使该发明更具实际应用性。

附图说明

通过阅读参照以下附图对非限制性实施例所作的详细描述，本发明的其它特征、目的和优点将会变得更明显：

图1是本发明方法的压缩/解压流程图；

图2是本发明方法的第一阶段在无参考序列时的压缩方法示意图；

图3是本发明方法的第一阶段在有参考序列时的压缩方法示意图；

图4是本发明方法中第二阶段与连续上下文对应地所选择的非连续上下文示意图；

图5是本发明方法的第二阶段方法示意图。

具体实施方式

下面结合具体实施例对本发明进行详细说明。以下实施例将有助于本领域的技术人员进一步理解本发明，但不以任何形式限制本发明。应当指出的是，对本领域的普通技术人员来说，在不脱离本发明构思的前提下，还可以做出若干变形和改进。这些都属于本发明的保护范围。

本发明提供一种基于非连续上下文建模和最大熵原则的基因压缩方法，通过结合非连续上下文建模思想与传统的连续上下文建模思想提高统计模型全面性，并根据最大熵原则得到Logistic回归模型确定更加准确的预测结果送入算术编码器。另外还考虑了有参考序列和无参考序列两种情况，提高了方法的实用性。以下提供本发明的实施例。

如图1所示，本实施例的压缩过程包括如下步骤：

第一阶段：

步骤一，在基因序列内部(无参考序列时)或参考序列上(有参考序列时)搜索近距离和远距离的准确重复(包括直接重复和互补回文)，当重复长度超过设定值时以该重复发生的属性、位置和长度高效表示基因序列的重复子序列，高效且针对性强。具体方法见图2和图3；

第二阶段：

步骤二，首先选择好11个传统的连续上下文模型(模型上下文选择分别为当前字节的前1,2,4,6,8,10,11,12,13,14,16个字节)、1个动态马尔科夫模型、1个重复匹配模型、11个利用非连续的比特位组合作为上下文来进行预测的模型(具体见图4)、1个距离模型和1个间接模型，然后计算每个模型对于当前比特位的预测概率值；

步骤三，当所有上下文模型计算出预测概率后，根据Logistic回归模型得到最终预测概率值；

步骤四，将步骤四的预测概率值送入算术编码器(流程图见图5)，输出码流。

如图2和图3所示，本实施例压缩过程中的第一阶段具体实施包括如下细节：

1、对于无参考序列的情况，在自身序列当前编码位置之前的64K(即64*1024)窗口内搜寻准确重复(包括直接重复和互补回文)项，对最终找到的最长重复(且要求大于重复长度设定值)用四倍体<重复标志‘D’，重复类型‘r’，重复发生的偏移位置Offset，重复长度M>来替代重复的子序列，并用log-skewed编码方法和差分+Gamma编码方法分别对四倍体中的偏移位置Offset和重复长度M进一步压缩；

2、对于有参考序列的情况，从参考序列上次重复位置开始处的窗口内双向搜寻准确重复(仅包括直接重复)项，对最终找到的最长重复(且要求大于重复长度设定值)用三倍体<重复标志‘I’，在参考中的偏移位置Offset，重复长度M>来替代重复的子序列，并用log-skewed编码方法和差分+Gamma编码方法分别对三倍体中的重复长度M和偏移位置Offset进一步压缩。

如图4所示，本实施例的第二阶段与11个连续上下文相对应地选择了11个比特位非连续组合的上下文模型，选择的原则是尽可能与连续上下文相互互补并能产生很好的压缩效果，图中纵向代表模型代码，横向代表距离当前编码比特的位数，黑色代表该处被选为上下文，空白则代表未被选中。

实施效果

依据上述步骤，在无参考序列时，实验用基因序列分别是来源于美国生物科技信息中心(National Center for Biotechnology Information，(NCBI))的5个细菌序列和10个不同物种的DNA序列；在有参考序列时，实验用基因序列分别为首个韩国人基因的两个版本KOREF_20090131、KOREF_20090224和第一个亚洲人的二倍体基因YH。关键参数的设置为：

(1)无参考序列时的允许最小重复长度设置为25(默认值)；

(2)有参考序列时的允许最小重复长度设置为50(默认值)；

(3)有参考序列时的窗口大小通过先对目标序列与参考序列的差别度进行测试然后决定，即当差别率小于0.65％时默认窗口大小为[-12,650]；当差别率介于0.65％与5％之间时默认窗口大小设置为[-12,812]；否则，窗口大小设置为[-12,11560]。

当无参考序列压缩时，本实施例比较了采用本发明所述的方法以及Duc Cao等人提出的基于多个专家模型和贝叶斯平均法的XM200、XM500方法和Armando J.Pinho等人提出的基于常规连续上下文模型相组合的FCM-M、FCM-Mx等方法的性能：

对于来源于美国生物科技信息中心(National Center for BiotechnologyInformation，(NCBI))的5个细菌序列，本发明方法、Duc Cao的XM500以及Armando J.Pinho等人的FCM-Mx得到的平均压缩比特率分别为1.7204bpb(bits per base)、1.7870bpb以及1.7543bpb，分别相当于节省了3.73％和1.93％的比特率。对于来源于美国生物科技信息中心的10个不同物种的DNA序列，本发明方法、Duc Cao的XM200以及Armando J.Pinho等人的FCM-M得到的平均压缩比特率分别为1.8457、1.8587以及1.8753，分别相当于节省了0.7％和1.58％的比特率。

当有参考序列压缩时，本实施例比较了采用本发明所述的方法以及ArmandoJ.Pinho等人提出的基于多个参考专家模型的GReEn方法和Wang,Congmao等人提出的基于直接表示参考与目标序列差别的GRS等方法的性能：

当以KOREF_20090131为参考压缩KOREF_20090224时，本方法、Wang,Congmao等人的GRS以及Armando J.Pinho等人的GReEn得到的压缩倍数分别为186.1倍、158.2倍以及171.4倍，分别相当于提升了17.6％和8.6％。当以KOREF_20090224为参考压缩YH时，Wang,Congmao等人的GRS由于目标基因与参考基因的差别大于其所规定的设定值而无法得到最终结果。本发明方法和Armando J.Pinho等人的GReEn得到的压缩倍数分别为177.2倍和98.1倍，相当于提升了80.6％。

实验表明，无论是在无参考序列时较之于Duc Cao等人提出的基于多个专家模型和贝叶斯平均法的XM200、XM500方法和Armando J.Pinho等人提出的基于常规连续上下文模型相组合的FCM-M、FCM-Mx等方法，还是在有参考序列时较之于Armando J.Pinho等人提出的基于多个参考专家模型的GReEn方法和Wang,Congmao等人提出的基于直接表示参考与目标序列差别的GRS方法，本发明提出的基于非连续上下文建模和最大熵原则的基因压缩方法均能显著提高压缩效率，实现高效存储。

从上述实施例可以看出，本发明方案中同时考虑了有参考序列和无参考序列两种实际情况，提出了将压缩过程分为两个步骤，分别对基因内或基因间的重复序列和非重复序列单独处理，以提高本方案的压缩效率。在处理重复序列时，使用简洁的四倍体(无参考序列)或三倍体(有参考序列)表示法，并对表示后的分量进一步压缩。在处理非重复序列时，在传统的连续上下文模型中引入了非连续上下文比特组合建模思想，更全面地表示基因内部相互关系信息；并采用由最大熵原则得到的Logistic回归模型合并所有上下文模型的预测结果，克服了以往贝叶斯平均法过度依赖单一模型而忽略其他模型有效性的缺点，得到更加准确的预测概率值送入算术编码器。该方案相对于现有的无参考序列和有参考序列压缩方案均可明显提高压缩效率。

以上对本发明的具体实施例进行了描述。需要理解的是，本发明并不局限于上述特定实施方式，本领域技术人员可以在权利要求的范围内做出各种变形或修改，这并不影响本发明的实质内容。

Claims

1.一种基于非连续上下文建模和最大熵原则的基因压缩方法，其特征是，该方法包括两大阶段：

在第一阶段，对于无参考序列的压缩，通过在已压缩序列内向前搜索的方法找到基因序列内部近距离和远距离的重复，包括直接重复和互补回文，以重复发生的属性、位置和长度高效表示基因序列的重复子序列；对于有参考序列的压缩，通过在参考序列内向前向后双向搜索的方法找到目标序列与参考序列中的重复部分，并通过重复发生的位置和长度高效表示重复子序列达到压缩目的；上述两种情况下剩余的非重复序列均将送入第二阶段的以比特位为单位的统计编码器；

第二阶段，所述统计编码器由预测器和算术编码器构成，预测器对于非重复序列，采用在传统的连续型上下文模型的基础上引入非连续比特上下文组合式的建模方案，并利用基于最大熵原则产生的Logistic回归公式合并各个模型单独产生的预测概率，从而得到最终更准确的预测概率值送入算术编码器，输出码流。

2.根据权利要求1所述的基于非连续上下文建模和最大熵原则的基因压缩方法，其特征是，所述第一阶段中：将已压缩序列重复部分与非重复部分分开，对于能在基因序列内部和参考基因序列上找到重复的子序列分别采用四倍体<重复标志‘D’，重复类型‘r’，重复发生的偏移位置Offset，重复长度M>表示和三倍体<重复标志‘I’，在参考中的偏移位置Offset，重复长度M>表示，并对该四倍体表示或三倍体表示进一步压缩。

3.根据权利要求2所述的基于非连续上下文建模和最大熵原则的基因压缩方法，其特征是，所述对该四倍体表示或三倍体表示进一步压缩，具体为：对于无参考序列中的四倍体表示，用log-skewed编码方法和差分+Gamma编码方法分别对偏移位置Offset和重复长度M进一步压缩；对于有参考序列中的三倍体表示，用log-skewed编码方法和差分+Gamma编码方法分别对重复长度M和偏移位置Offset进一步压缩。

4.根据权利要求1所述的基于非连续上下文建模和最大熵原则的基因压缩方法，其特征是，所述的用在参考序列内向前向后双向搜索到的最长重复序列表示基因序列间的重复方案时，自适应地选取寻找重复的窗口大小；具体为：在已知参考序列的情况下，对目标序列进行压缩之前，首先计算两个序列的差别率，根据差别率的不同取值分配给其不同的滑动窗窗口大小，从而节省压缩中寻找重复步骤的时间，并提高方案的自适应性。

5.根据权利要求1-4任一项所述的基于非连续上下文建模和最大熵原则的基因压缩方法，其特征是，所述预测器引入非连续比特上下文组合式的建模方案，具体为：除了传统的连续上下文模型外，根据实验测试又选择了13个非连续上下文模型，即选取在位置上相互互补且对测试实验产生很好压缩效果的11个利用非连续的比特位组合作为上下文来进行预测的模型、1个距离模型和1个间接模型，以达到优化整体压缩性能的目的。

6.根据权利要求5所述的基于非连续上下文建模和最大熵原则的基因压缩方法，其特征是，所述利用基于最大熵原则产生的Logistic回归公式合并各个模型单独产生的预测概率，具体为：

①、首先根据当前需要编码的比特位y_j的上下文，对于第i个模型假设有m_ij个0比特，n_ij个1比特，计算每个模型的独立预测概率P_i，并求出每个概率的伸展形式

②、基于最大熵原则得到使得熵最大的合成概率为Logistic回归模型：其中M_c为模型的总个数，M_c＝11+13＝24，w_i是第i个模型此时的权值大小；e是自然对数的底数；

③、将合成概率送入算术编码器输出；

④、根据牛顿法，以使编码消耗C＝-ln(1-|E|)最小，即预测误差E＝y_i-P(y|t)最小为目标求出更新后的w_i，为：

w_{i} = w_{i} - \frac{\frac{\partial C}{\partial_{w_{i}}}}{\frac{\partial^{2} C}{\partial_{w_{i}^{2}}}} = w_{i} - \frac{1}{t_{i} (1 - y_{j} - P (y_{j} | t_{1}^{M_{c}}))};

⑤、使用新的w_i返回第①步，开始预测比特位y_j+1。