CN103546160B - 基于多参考序列的基因序列分级压缩方法 - Google Patents

Abstract

本发明提出了一种基于多参考序列的基因序列分级压缩方法，该方法首先将BAM格式文件转化成SAM格式的文件，SAM格式的基因序列由11个强制域和多个可选域构成，将原文件按域提取成12个独立文件，然后对12个文件进行并行压缩：（1）对‘Sequence’域采用基于多个参考序列逐步减半序列长度的分级压缩方法；（2）对于‘Quality Value’域采用k均值聚类结合上下文建模PPMVC压缩的方法；该方案相对于现有的同格式的压缩方案既提高或保证了压缩效率，又提供了压缩等级的多选择性，使得其更有适应性与扩展性。

Description

基于多参考序列的基因序列分级压缩方法

技术领域

本发明涉及一种面向超大规模SAM/BAM（两种序列比对格式标准）格式基因序列的信息压缩方法，具体是一种利用目标序列与相同物种的不同参考序列间的相似重复性，多次改变参考序列并减少短序列长度的分级基因压缩方法。

背景技术

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

迄今为止，所有的基因压缩研究主要围绕三种序列格式展开：（1）最初的DNA基集合形式的FASTA格式（一种基于文本用于表示核苷酸序列或氨基酸序列的格式）；（2）未比对短序列形式的FASTQ格式（一种基于文本用于表示核苷酸序列或氨基酸序列和相应质量信息的格式）；（3）比对后短序列形式的SAM/BAM格式。由于序列比对是进行序列分析和处理的首要步骤，且SAM/BAM格式包含了丰富完整的基因比对信息，因此近几年的基因分析与压缩研究均聚焦于SAM（由于BAM是SAM的二进制表示，因此压缩BAM格式时可先将其解压，然后还原为压缩SAM格式的序列）格式。2011年Muhammad等人在PLOSOne期刊的“Improvingtransmissionefficiencyoflargesequencealignment/map(SAM)files”上提出了专门针对于SAM格式和特征的压缩方法SAMZIP，利用基本压缩技巧，如游程编码（Run-LengthEncoding，RLE），差分（Delta）编码，哈弗曼（Huffman）编码和字典编码，对SAM格式中的每一列独立处理。同年，Kozanitis等人在JournalofComputationalBiology期刊的“CompressinggenomicsequencefragmentsusingSlimGene”上提出了压缩方法SLIMGENE，其中首次尝试了对SAM格式中的QualityValue（质量值）项进行保证后续处理的有损压缩。考虑到相同物种核苷酸差异性的微小不同（人类核苷酸差异性仅为0.1%左右），研究学者们开始将已知的参考序列引入SAM格式的基因压缩中。2011年，Fritz等人在Genomeresearch期刊的“EfficientstorageofhighthroughputDNAsequencingdatausingreference-basedcompression”上发表了一种基于参考序列的压缩方法CRAM，将SAM格式序列中的每一子序列（Read）与参考基因做比对，然后压缩其比对结果，用较少比特表示相应的差别。2012年Hach等人在Bioinformatics期刊的“SCALCE:boostingsequencecompressionalgorithmsusinglocallyconsistentencoding”上提出了另一种基于参考序列的压缩方法SCALCE，该方法基于短序列的重组以适应参考序列的局部特性，便于进一步比对压缩。以上所述的这些方法将SAM格式分离成多个独立项并行处理。在压缩DNA序列时将其视作由特殊字符（即‘A’，‘C’，‘T’，‘G’，‘N’）构成的长字符串，从数据的构成特点和自身冗余性出发进行整体处理，有效的提高了压缩效率和压缩时间。但总体而言SAM格式的大基因压缩技术仍处于起步阶段，已知公开的参考基因的信息并没有被充分发挥利用，目标信息中的准确比对率也仍待提高。

发明内容

针对现有技术中的缺陷，本发明的目的是提供一种更加有效的基于多参考序列的SAM/BAM格式基因序列分级压缩方法。

本发明是通过以下技术方案实现的：该方法通过利用多个公开的参考基因序列，并将短序列长度逐步减半，多次比对目标序列以提高被压缩序列的比对准确率，进而提高压缩效率。另外，对于SAM格式中的QualityValue项，本发明采用了用户可指定压缩等级的k均值（k-means）聚类后再压缩的策略，提高了方法的广泛实用性。由于BAM格式是SAM格式的简单转化，因此下面仅讨论针对SAM格式的压缩技术，对于BAM格式，只需将其转换成SAM格式然后再按照SAM格式的压缩方法处理即可。

本发明所述的基因序列分级压缩方法，该方法首先将BAM格式文件转化成SAM格式的文件，SAM格式的基因序列是由比对工具输出的文本形式文件，它由11个强制域（Field）和一系列可选域（看作第12个域）构成，因此本发明首先将原文件按域提取成12个独立文件，然后对12个文件进行并行压缩：

（1）对‘Sequence’域采用基于多个参考序列逐步减半序列长度的分级压缩方法；

（2）对于‘QualityValue’域采用k均值聚类结合上下文建模PPMVC（部分串匹配的预测（PredictionwithPartialstringMatching）法的一种扩展）压缩的方法；

（3）对于剩下的十个域采用基于域内特征和域间相关性的压缩方法。

进一步的，所述对‘Sequence’域采用基于多个参考序列逐步减半序列长度的分级压缩方法，具体为：利用快速比对工具SOAP3将SAM/BAM格式基因序列文件的‘Sequence’域中的短序列分线程地与参考序列作比对，对于准确匹配序列高效压缩，对于非准确匹配和未匹配的短序列，将其序列长度减半，即一个序列分成长度相同的两个序列，并改变参考序列，再进行第二次比对，得到比对结果，如此重复三至四次结束，剩余的非准确匹配和未匹配的短序列进行PPMVC编码即可。这样多次比对后需要处理压缩的非准确匹配和未匹配序列变少。

优选的，所述对于准确匹配序列高效压缩，该压缩方法基于比对结果的特点，富有针对性。具体为：对于准确比对的子序列（Read），使用<Read编号，参考序列上重复发生的染色体号，参考序列上重复发生的偏移位置，重复类型>这四个量来替代目标序列上重复的子序列（Read），分别使用差分编码+Huffman编码、游程编码、差分编码+Huffman编码和游程编码来压缩这四个分量。

进一步的，所述的对于‘QualityValue’域采用k均值聚类结合上下文建模PPMVC压缩的方法，具体为：采用k均值聚类法将n个Q_ASCII值（代表对应基的比对质量分值的ASCII码值）聚成k类，使得聚类后每类内所有Qualityvalue的值与聚类前的值差值平方最小，然后采用基于上下文建模和统计的自适应文本压缩方法PPMVC压缩聚类后的‘QualityValue’文件。这样同等压缩条件下失真较小。

进一步的，所述的对于剩下的十个域采用基于域内特征和域间相关性的压缩方法，具体为：

对于‘QNAME’域，用‘0’表示之前未出现过的QNAME，用逐渐递增的数字编号与当前位置只差表示之前已经出现的QNAME，然后采用Huffman编码压缩这些非均匀分布的小型数值；

对于‘FLAG’域，用一个字节表示1～255之间的数值，用三个字节（即，0,x/256andx%256）表示其它数值，然后采用Huffman编码压缩变换后的数值；

对于‘RNAME’域，用相同的数字标记整个SAM文件中的相同的参考序列名字，记录下来所有参考序列，然后用游程编码进行压缩；

对于‘POS’域，采用差分编码+Huffman编码；

对于‘MAPQ’域，采用游程编码；

对于‘CIGAR’域，采用Lempel-Ziv-WelchLZW字典压缩方法；

对于‘MRNM’域，采用游程编码；

对于‘MPOS’域，结合‘MRNM’域的字符串，采用差分编码+Huffman编码；

对于‘TLEN’域，‘TLEN’域的值与‘MPOS’域减去‘POS’域的值的差(即，TLEN-(MPOS-POS))的绝对值服从于一个有限的字符集，因此，对于该域的压缩，本发明结合‘POS’，‘MPOS’和‘MRNM’三个域的信息采用Huffman编码压缩变换后的值；

对于‘OPTIONAL’域，使用bzip2压缩工具。

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

本发明所提出的基于多参考序列的SAM/BAM格式基因序列分级压缩方法，提高了基因压缩的效率和完整性。本发明将多个公开的基因序列作为参考结合使用，充分利用了相同物种间的基因相似性；将非准确匹配的子序列进行截断式再次比对，克服了以往方法无针对性地对准确匹配和非准确匹配序列统一处理的缺点，提高了准确比对率，进而节省了压缩比特；对‘QualityValue’域采用用户可指定压缩级的k均值聚类法，在提高压缩效率的同时保证了质量值的准确度；对剩余域既考虑了域间独立性也考虑了域内特征和部分域间的分布相关性，针对性的进行变换然后再编码，挖掘出了SAM序列格式中的潜在信息，进一步提高了压缩效率。

附图说明

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

图1是本发明实施例的压缩流程图；

图2是本发明实施例的‘Sequence’域经过不同参考序列不同子序列长度多次比对的示意图；

图3是本发明实施例的‘Sequence’域在某一实例中的压缩效果图。

具体实施方式

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

本发明首先将BAM格式文件转化成SAM格式的文件，由于SAM格式的基因序列是由比对工具输出的文本形式文件，它由11个强制域（Field）和一系列可选域（看作第12个域）构成，因此本发明首先将原文件按域提取成12个独立文件（这12个独立文件对应的是上述11个强制域和第12个域，本实施例中即以下所描述的各个域），然后对12个文件进行并行压缩。对‘Sequence’域的处理为本发明的第一部分；对‘QualityValue’域的处理为本发明的第二部分；对剩余10个域的压缩为本发明的第三部分。其中，‘Sequence’域和‘QualityValue’域占据了整个基因序列的50%左右内容且不易压缩，是本发明的处理重点。每个部分的编码经过如下：

1、‘Sequence’域的多参考分级式压缩方案

对由‘A’，‘C’，‘T’，‘G’，‘N’五个基构成的‘Sequence’域，将经过以下编码过程：

（1）通过快速比对工具SOAP3将每一Read分别与指定的参考序列进行比对；

（2）对于准确比对的子序列（Read），使用<Read编号，参考序列上重复发生的染色体号，参考序列上重复发生的偏移位置，重复类型>这四个量来替代目标序列上重复的子序列（Read）；

（3）分别使用差分编码+哈弗曼编码、游程编码、差分编码+哈弗曼编码和游程编码压缩<Read编号，参考序列上重复发生的染色体号，参考序列上重复发生的偏移位置，重复类型>中的四个分量；

（4）对于未准确比对或未匹配的子序列（Read），将其长度减半，并选择其他参考序列（可选）作为参考，然后重复（1）-（4），共重复三次。

2、‘QualityValue’域的用户可指定编码级的k均值聚类法压缩方案

对由多个（一般情况下有51个）不同字符构成的‘QualityValue’域，每个字符的ASCII值与当前位置所对应的基的误差概率的关系为：

$Q_{\mathrm{ASCII}}=Q_{\mathrm{phred}}+33,P_e=10^{-\frac{Q_{\mathrm{ASCII}}-33}{10}},Q_{\mathrm{phred}}=-10{\log }_{10}{P}_e$

其中P_e是当前位置基的测序错误率，Q_phred是其转化为整数后的值，Q_ASCII是可以在文本中显示出来的质量值。

采用k均值聚类法将n个Q_ASCII值聚成k类，使得聚类后每类内所有Qualityvalue的值与聚类前的值差值平方最小，即

$\arg \underset{S}{\min }\underset{i=1}{\overset{k}{\mathrm{\&Sigma;}}}\underset{q_j\&Element;S_i}{\mathrm{\&Sigma;}}{\left|\right|10^{-\frac{q_j-33}{10}}-{10}^{-\frac{u_i-33}{10}}\left|\right|}^2\&CenterDot;{\mathrm{num}}_j$

以保证聚类后的误差概率P_e的变化最小。其中u_i是类别S_i的均值，num_j是整个序列中q_j出现的次数。

最后，采用基于上下文建模和统计的自适应文本压缩方法PPMVC压缩聚类后的‘QualityValue’文件。

3、剩余域的基于域内特征和域间联系的压缩方案

根据剩余域域内和域间的数值关联及特征性，分别对每一域采用如下编码方法：

QNAME:用‘0’表示之前未出现过的QNAME，用逐渐递增的数字编号与当前位置只差表示之前已经出现的QNAME，然后采用哈弗曼编码压缩这些非均匀分布的小型数值。

FLAG:用一个字节表示1～255之间的数值，用三个字节（即，0,x/256andx%256）表示其它数值，然后采用哈弗曼编码压缩变换后的数值。

RNAME:用相同的数字标记整个SAM文件中的相同的参考序列名字，记录下来所有参考序列，然后用游程编码进行压缩。

POS:差分编码+哈弗曼编码。

MAPQ:游程编码。

CIGAR:Lempel-Ziv-Welch(LZW)字典压缩方法。

MRNM:游程编码。

MPOS:结合‘MRNM’域的字符串，采用差分编码+哈弗曼编码。

TLEN:因为(TLEN-(MPOS-POS))的绝对值服从于一个有限的字符集，因此，对于该域的压缩，结合‘POS’,‘MPOS’and‘MRNM’三个域的信息采用哈弗曼编码压缩变换后的值。

OPTIONAL:bzip2压缩工具。

解压缩时，首先按照上述同样的步骤并行恢复出原基因序列中每一域的部分，然后将所有域合并得到完整无损的SAM或BAM基因文件。

以下提供本发明的一个具体应用实例：

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

步骤一，利用samtools工具将SAM/BAM文件提取成12个独立文件，其中每个域为一单独文件，所有的可选域（即OptionalFields）提取成一个文件，然后送入每个域单独的编码器；

步骤二，对有规律性的十个域采用上述‘剩余域的基于域内特征和域间联系的压缩方案’中描述的方法进行压缩；

步骤三，对重点的‘Sequence’域采用图1所示的分级压缩结构，其目的是为了提高短序列的准确匹配率。如图2所示，以单一序列为参考进行一次比对时的匹配率远远低于多次比对时的匹配率，更低于以多个序列为参考时多次比对的匹配率。图中的hg19和HuRef是两个具有代表性的参考序列，EMRs（ExactMappedReads）,IMR（InexactMappedReads）,UMR（UnmappedReads）分别指准确匹配序列、非准确匹配序列、未匹配序列；

步骤四，对重点的‘QualityValue’域先采用k均值聚类法，然后采用基于上下文建模和统计的自适应文本压缩方法PPMVC（命令行参数为‘e–o4–r1’）压缩聚类后的‘QualityValue’文件；

步骤五，将12个域的压缩文件打包成一个完整的压缩文件，作为最终的压缩结果输出。

如图3所示，本实施例压缩过程中的步骤三具体实施包括如下细节：

1、首先利用快速比对工具SOAP3将SAM/BAM文件的‘Sequence’域中的上百万个短序列分线程地与参考序列作比对，得到以<Read编号，参考序列上重复发生的染色体号，参考序列上重复发生的偏移位置，重复类型>标注的比对结果；

2、对于非准确匹配和未匹配的短序列，将其序列长度减半，即一个序列分成长度相当的两个序列，并改变参考序列，再进行第二次比对，得到比对结果；

3、返回第一步，如此重复三到四次直到此时的序列长度低于15，剩下的非准确匹配和未匹配的短序列数已非常少，如图3所示，经过分别以hg19.fa,HuRef.fa,hg19.fa和hg19.fa为参考序列的四次比对后，剩余的非准确匹配和未匹配的短序列比例仅为0.01%，对它们再进行PPMVC（命令行参数为‘e-o8’）编码即可。

实施效果：

依据上述步骤，选取了多组不同种类的下一代测序数据（NextGenerationSequencing,NGS）：7组来自千人基因组项目（1000GenomesProject），一组来自于ChIP-Seq的mouse数据，一组来自于RNA-Seq的E.coli数据，还有一组来自于癌症基因集（theCancerGenomeAtlas,TCGA）的基因数据。本实施例比较了采用本发明所述的方法，MarkusHsi-YangFritz等人提出的基于目标序列与参考序列间的区别编码的CRAM方法和JamesK.Bonfield等人提出的基于序列中每个base的位置坐标建模的Samcomp方法的性能：

对于来自不同项目组织的NGS基因数据，本发明所提出的方法均取得了明显优于CRAM和与Samcomp相当的压缩率（压缩后的文件大小/压缩前文件大小）。相对于CRAM的无损压缩方法，本发明的无损模式在BAM文件上所产生的压缩率降低了6%-20%，达到了0.5-0.65的压缩率，节省了35%-50%的存储空间。在压缩时间相当的情况下，其解压时间小于CRAM。相对于Samcomp，本发明在保证压缩效率相当的情况下，提供了‘QualityValue’压缩等级的多选择性，使得其更有适应性与扩展性。

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

Claims (5)

1.一种基于多参考序列的基因序列分级压缩方法，其特征是，首先将BAM格式文件转化成SAM格式的文件，SAM格式的基因序列由11个强制域和多个可选域构成，将可选域作为第12个域，原文件按域提取成12个独立文件，然后对12个文件进行并行压缩：

(1)对‘Sequence’域采用基于多个参考序列逐步减半序列长度的分级压缩方法；

(2)对于‘QualityValue’域采用k均值聚类结合上下文建模PPMVC压缩方法；

(3)对于剩下的十个域采用基于域内特征和域间相关性的压缩方法。

2.根据权利要求1所述的基于多参考序列的基因序列分级压缩方法，其特征是，所述对‘Sequence’域采用基于多个参考序列逐步减半序列长度的分级压缩方法，具体为：利用快速比对工具SOAP3将SAM文件的‘Sequence’域中的短序列分线程地与参考序列作比对，对于准确匹配序列高效压缩，对于非准确匹配和未匹配的短序列，将其序列长度减半，即一个序列分成长度相同的两个序列，并改变参考序列，再进行第二次比对，得到比对结果，如此重复三至四次结束，剩余的非准确匹配和未匹配的短序列采用上下文建模PPMVC压缩方法进行压缩。

3.根据权利要求2所述的基于多参考序列的基因序列分级压缩方法，其特征是，所述的对于准确匹配序列高效压缩，具体为：对于准确比对的子序列Read，使用<Read编号，参考序列上重复发生的染色体号，参考序列上重复发生的偏移位置，重复类型>这四个量来替代目标序列上重复的子序列，分别使用差分编码+哈弗曼编码、游程编码、差分编码+哈弗曼编码和游程编码来压缩这四个分量。

4.根据权利要求1-3任一项所述的基于多参考序列的基因序列分级压缩方法，其特征是，所述的对于‘QualityValue’域采用k均值聚类结合上下文建模PPMVC压缩方法，具体为：采用k均值聚类法将n个Q_ASCII值聚成k类，使得聚类后每类内所有Qualityvalue的值与聚类前的值差值平方最小，然后采用基于上下文建模PPMVC压缩方法处理聚类后的‘QualityValue’文件。

5.根据权利要求1-3任一项所述的基于多参考序列的基因序列分级压缩方法，其特征是，所述的对于剩下的十个域采用基于域内特征和域间相关性的压缩方法，具体为：

对于‘QNAME’域，用‘0’表示之前未出现过的QNAME，用逐渐递增的数字编号与当前位置只差表示之前已经出现的QNAME，然后采用哈弗曼编码压缩这些非均匀分布的小型数值；

对于‘FLAG’域，用一个字节表示1～255之间的数值，用三个字节即0,x/256和x％256表示其它数值，然后采用哈弗曼编码压缩变换后的数值；

对于‘RNAME’域，用相同的数字标记整个SAM文件中的相同的参考序列名字，记录下来所有参考序列，然后用游程编码进行压缩；

对于‘POS’域，采用差分编码+哈弗曼编码；

对于‘MAPQ’域，采用游程编码；

对于‘CIGAR’域，采用LZW字典压缩方法；

对于‘MRNM’域，采用游程编码；

对于‘MPOS’域，结合‘MRNM’域的字符串，采用差分编码+哈弗曼编码；

对于‘TLEN’域，‘TLEN’域的值与‘MPOS’域减去‘POS’域的值的差即TLEN-(MPOS-POS))的绝对值服从于一个有限的字符集，对于该域的压缩，结合‘POS’,‘MPOS’和‘MRNM’三个域的信息采用Huffman编码压缩变换后的值；

对于‘OPTIONAL’域，使用bzip2压缩工具。