基于高通量测序的短序列生物数据压缩研究

孟 倩

(复旦大学计算机科学技术学院 上海 200433)

基于高通量测序的短序列生物数据压缩研究

孟 倩

(复旦大学计算机科学技术学院 上海 200433)

高通量测序技术(NGS)的发展带来了测序数据量的极速增长，给数据的存储和传输带来了极大的压力。数据压缩技术是解决这个问题的重要方法。传统的压缩方法并没有很好地利用数据本身的特性。因此，计算机学者们关注于NGS测序数据专用的压缩方法。全面总结针对高通量测序技术产生的Fastq和Fasta数据的压缩算法，介绍了Fastq和Fasta数据的特点，总结了目前常用的压缩方法。并通过不同物种、不同测序平台、不同规模的测序数据对多个具有代表性的压缩工具进行测试，比较它们的压缩性能并且验证相应的工具特点，为研究人员提供工具选择指导或改善工具性能提供帮助。最后总结阐述短序列数据压缩工具存在的问题和发展趋势。

数据压缩 短序列数据压缩 高通量测序

0 引 言

测序是对生命中的基本元素{ACGT}及其修饰和衍生物进行测定和解读。20世纪诞生的第一代测序技术Sanger测序技术能对几百到几千的DNA片段进行快速准确的读取，这直接促成了人类历史上一项伟大的计划—人类基因组计划(HPG)的诞生，它对人类30亿碱基进行测序，耗资30亿美金，前后耗时10余年，涉及6个国家，最终于2001年公布草图。2005年以后出现的高通量测序技术(High-Throughput Sequencing)，也叫下一代测序技术(Next-Generation Sequencing)或深度测序，是生命科学领域内的一项重要的技术变革，给个人全基因组测序带来了可能。主要以Illumina 公司的Solexa基因组分析仪，Roche公司的454测序仪和ABI公司的SOLiD测序仪为代表，又以Illumina应用最广，它可以同时处理几百万到几亿的基因数据，成本低，效率高，速度快。第三代测序相比于第二代测序在成本和效率上并没有显著的提高。因此，市场上还是以第二代测序仪为主。高通量测序的低成本和高效率，也在产生海量的测序数据。计算机是存储和处理这些数据的重要工具。但是数据量的增长速度还是远超过计算机的存储能力，这给数据的存储带来了极大的压力。同时，这些测序数据需要传输，但是带宽的增长速度也是远远无法满足需求。如何存储与传输大量的测序数据，已经成为了测序领域的重要问题。

数据压缩是有效解决测序数据带来的存储和传输问题的一种重要方法。通用的压缩方法如bzip、gzip、LZMA等对数据的压缩结果并不是特别理想，并没有充分利用生物数据本身的特点。因此，计算机学者们开始关注于NGS测序数据专用的压缩方法。

NGS测序数据主要包括原始测序数据Fastq数据，拼接数据Fasta和比对数据SAM/BAM数据，以及变异文件VCF。

本文较全面地对目前流行的NGS的短序列Fastq和Fasta格式的数据压缩方法进行分类总结，并通过真实的序列数据对几种压缩方法进行测试评估，比较其压缩性能。最后讨论目前的方法存在的问题和未来的发展趋势。

1 数据格式

1.1 Fastq格式

Fastq格式是NGS测序出来的原始数据格式。Fastq格式的序列一般有4行。第一行是元数据，由‘@’开始，后面跟着序列的描述信息，通常被认为是标志行。第二行是序列，由{ACGTN}组成，其中N代表不确定，表示该位置可能是ACG或T。第三行由‘+’开始，后面也可以跟着序列的描述信息，通常和@后面的内容一样。第四行是第二行序列的质量评价，因此字符数和第二行一样，表示对应碱基的错误概率的相关信息，越大代表错误概率越低，通常用ASCII表示。对Fastq格式的压缩包括对元数据、序列和质量分数三者的压缩。

1.2 Fasta格式

Fasta格式通常有两行。第一行由‘@’开始，后面跟着序列的描述信息，通常被认为是标志行。可能包括仪器名，测序平台，短读的标识等信息。第二行是序列，序列信息通常由ACGT四个碱基组成，有时会有N代表不确定该位置是ACG还是T。仅从数据内容的格式上而言，相比于Fastq格式数据，它缺少了Fastq格式中的第三行和第四行内容，因此相比于Fastq格式压缩，它的压缩算法主要是集中于对元数据和序列的压缩。

2 NGS数据的压缩技术

现有的NGS数据压缩算法主要可以分为两类：基于参考基因组压缩以及非基于参考基因组压缩。在对生物数据压缩的过程中，大多数的压缩工具采用了对数据的各个部分(元数据)分别压缩的方式。对每个部分可能采取不同的压缩策略和方法。常用的压缩方法可归类为字典编码、熵编码和变换变换。字典编码是利用字典把每个字符串编码为一个标识，字典保存字符串，可以是静态的，也可以是动态(自适应)的[1]，如增量编码、LZ等。熵编码在NGS中有很多运用，典型的熵编码有：算术编码、哈夫曼编码、游程长度编码、马尔可夫模型等。变换编码包括经典的BWT(Burrows-Wheeler Transformation)算法[2]、聚类算法等，通过变换编码，把相似字符聚集在一起，方便后续的游程长度编码或者MTF(Move To Front)进行压缩，达到较好的压缩率。

2.1 元数据

ID部分是数据的标志行。常常存储数据的编号，测序平台等信息。有些信息譬如说仪器名等信息会在每个短读中都大量出现，这就造成了大量冗余。考虑到这些冗余信息，常常把数据分为固定的部分和可变部分两部分进行压缩，采用增量编码的方式来进行，例如QUIP[3]、FQC[4]等都是采用的这种方法。QUIP是把数据和上一条数据进行比较，对完全相同的部分、有规则的数字符号和无规则的数字符号采用不同的压缩方案。FASTQZ[5]的快速模式也采用了相似的概念来进行压缩，存储差异信息，包括差异数据的位置和差异类型以及差异数据。

有些算法采用了K阶马尔科夫模型也叫有限上下文模型的方式来对该部分进行压缩。在使用马尔科夫模型的时候，要注意阶数的选择，减少内存消耗。阶数越多，对数据的预测准确率越高，但是对内存的消耗也会越大。例如如果字符集为256个字符，K阶马尔科夫模型的每一次预测都有256k个参数。因此k的选择要兼顾到内存和压缩率，需要在两者之间找到一个平衡点。例如FASTQZ的慢速方法就是采用了混合有限上下文模型，把4种上下文模型按照一定的权重比例组合起来来减少预测错误率，提高压缩率。MFCompress[6]算法利用了连续的元数据有极大的关联性的特点，借助于要压缩元素之前的位置和要压缩元素所在短读的上一条短读的相关位置的关联性，采用了有限上下文模型。同时考虑到压缩率和内存的平衡，在压缩前先对元数据的不同的元素个数进行了统计，根据不同元素个数的大小来确定阶数。准确地说是大于50时，选择3阶模型，小于时选择4阶模型。

BEETL[7]算法认为该部分信息只是作为一个短读的唯一标识，因此为了提高压缩率，获得更好的压缩效果，选择抛弃这部分信息，改为在解压缩的时候根据短读内容生成它的唯一标识。DSRC2[8]算法可由用户选择对元数据采用有损压缩还是无损压缩。

2.2 序 列

序列是Fastq数据中最核心的信息，不能接受任何的信息损失，因此序列的压缩都是无损压缩。对序列的压缩策略包括：

(1) 利用和质量分数的一对一关系，对序列进行处理，使得序列数据更容易压缩。例如Fastqz算法利用了序列中的N和质量分数之间的对应关系，N在质量分数中一般都为0。因此它在处理序列的时候，直接删除所有的N元素，解压缩的时候在质量分数中为0的地方在对应的序列中插入N。

(2) 利用同一条数据间的相关性，采用k阶马尔科夫模型。QUIP[3]算法直接采用了12阶马尔科夫模型，12阶马尔科夫模型需要32bit来存储中间参数。

(3) 用BWT或者聚类算法处理数据，使得相似数据聚集在一起，更容易被压缩。BEETL[7]算法是典型的基于BWT算法进行改进的针对高通量测序技术的压缩算法。

更多的算法是对这几种策略进行部分修改后叠加。DSRC[9]以及DSRC2中的一种压缩策略也是采用了这种方法，把所有的非ACGT都转移到质量分数上，这样序列中只有ACGT四个字符。在接下来的压缩中，可以采用每个字符用2bit来表示，也可以采用哈夫曼编码，也可以采用算术编码和有限上下文模型。MFCompress[6]先对序列数据进行预处理，然后再采用混合有限上下文模型来进行压缩。对数据的预处理是将序列分为2个或者3个流程来单独处理。主流程的目的是使得该流程中只有0123四个元素，因此它把ACGT和acgt都替换成0123，所有非ACGT部分全部用0来替换。然后再用extrastream来标注所有的0代表的含义，因为序列中一般非ACGT部分很少，因此这样做可以有效减少压缩量。最后再用一个流程来区分大小写。

2.3 质量分数

QUIP采用了4阶马尔可夫模型，主要是因为质量分数中有四十多个不同的字符，当阶数越高，参数越多，而这限制了马尔科夫模型的阶数。KungFq[10]和FQZip[11]采用了采用了游程长度编码。LWFQZip[12]采用了限制游程长度编码，限制游程长度编码是对游程长度编码的一点改进。文献[13]在质量分数的压缩中，观察到“#”号出现以后一般后面就不会有别的质量分数。因此它去除了所有的“#”号以节省压缩数据量。同时它考虑到直接运用哈夫曼编码效率不高，因此把质量分数切割成k-mer大小的块，对每一块进行或者哈夫曼或者游程长度编码，具体由压缩率来决定。

NGC[14]算法是对匹配成功的质量分数和未匹配成功的质量分数采用了不同的压缩方案。FQC算法是采用了有损和无损压缩两个选择，无损压缩是对数据采用游程长度编码后再采用ppmd[15]压缩。在有损压缩中，是直接根据压缩率的需求度不同提供了三种压缩率不同信息损失率不同的压缩方案，第一级小于0 的用3来代替，大于3的奇数用偶数代替。第二级是小于3的用0来代替，大于3 的中每4个数用第2个数来代替。第三级是小于7的用0来代替，大于7的每8个数用第4个数来代替。但是它还提供了由有损压缩文件和补充文件还原出原数据的方式，把压缩过程中损失的信息单独记录，然后通过有损压缩解压后的文件和补充文件可以还原出原始文件，给后续数据传输带来了极大的方便。QualComp[16]算法是基于率失真理论，对质量分数用一个固定的失真指标量化到一个特定的比率中，以期望在对后续分析的影响和压缩率之间取得一个平衡。

有损压缩提高了压缩效率，但是丢失了部分信息，这可能给后续的数据分析带来挑战。目前这种压缩方法仍然面临着巨大的挑战。在近期涌现了很多优秀的质量分数有损压缩算法，如RQS[17]、QVZ[18]等。

2.4 数据转换技术

数据转换技术主要是指在数据压缩之前，对数据进行处理的技术。通过使用Burrows-WheelerTransformation(BWT)或者聚类算法等对数据进行处理，使得相似的聚集在一起，然后再进行压缩，以达到提高压缩率的结果。BEETL算法使用了BWT对输入的数据先进行了处理，使得文件更易于压缩，然后使用通用的压缩方法如7-zip，bizp2等对处理后的文件进一步压缩。SCALCE[19]算法通过局部一致性分析LocallyConsistentparsing(LCP)算法对有最大相同子串的数据识别后聚类，然后再用gzip等算法对数据进行后续压缩。

2.5 基于参考基因组压缩算法

(1)Cpv是指施工项目在考核期内，基于拟投入的安全成本，在确保计划安全保障水平的基础上安全成本的计划费用值。

同源基因组之间有很高的相似性，甚至可能高达99%。例如人类基因组，不同的人之间的基因相似度超过99.9%，而这些0.001 的不同就构成了不同的个体。相比于相同信息，差异信息是极少量的。如果已经有了同源基因组作为参考基因组，在已知参考基因组和差异信息的情况下，显然可以还原出原始数据。因此，在有参考基因组的情况下，对基因组的存储，就变成了对差异信息的存储。编码器只需要存储差异信息，这极大的提高了压缩率，节省了存储空间和传输时间。

参考基因组的选择可以是同源基因组，也可以是由自身拼接得到，基因组拼接技术也有很多，文献[20]是对基因组拼接技术的综述。

参考基因组压缩需要记录比对信息和差异信息，包括匹配位置、匹配长度、匹配类型、差异位置、差异类型和差异内容。差异类型可分为插入、删除和替代这三类。然后对这些记录信息压缩，可以采用哈夫曼编码、算术编码或者游程长度编码等。

迄今为止，有很多的高通量测序压缩算法采用了基于参考基因组压缩的方式，如SlimGene[21]，QUIP，NGC，SAMZIP等。QUIP除了给出了基于参考基因组压缩的方案，还给出了基于拼接的压缩方式，是对deBruijin采用了dlCBF(d-leftcountingBloomFilter)改进的拼接。基于参考基因组的拼接压缩理论上而言也是基于参考基因组压缩的方式。

2.6 并行处理

在对大数据进行压缩和解压缩的时候，时间成为了一个不可忽略的评判因素。MFCompress为了减少数据压缩和解压缩时间，把数据分割成几块来并行处理，一般而言压缩是双线程处理。但是解压缩的时候，线程数量可以由用户自己来决定。它把数据并行处理后再合并处理，达到了提高压缩时间的目的，同时也能有效地利用了计算机资源。

3 压缩算法性能评估

在这一节中，我们选择了4个数据用于测试现有的主要的fastq和fasta数据压缩算法的压缩率和压缩时间。我们选择不同物种(人类，小鼠，微生物)，来自于不同的数据的数据平台(Illumina,454)和不同的大小(从217MB到12GB)。所有数据来自于http://sra.dnanexus.com/，通过sratoolkithttp://trace.ncbi.nlm.nih.gov/Traces/sra/sra.cgi?view=software工具转换成fastq格式。数据如表1所示。

表1 测试数据集

测试环境采用64位Ubuntu14.04操作系统，硬件平台为4核8线程的3.5GHzIntel(R)Core(TM)i7-4770KCPU，16GB内存。在测试中尽量不运行别的程序。

首先我们对FASTQZ、QUIP、DSRC2、SCALCE和BEETL方法进行了压缩。表2展示的是对fastq格式数据压缩的结果，未做标注的都是采用默认压缩模式。SRR001540是由变长的454 数据组成，BEETL不能压缩变长的数据。也就是不能压缩SRR001540。但是我们对数据进行了格式化处理，所有短读统一长度为200。SOLiD数据格式和Illumina数据格式不一样，它含有的是Colorspace而并不是碱基ACGT。很多算法都不支持SOLiD格式数据，例如QUIP方法，这里我们没有选择SOLiD格式数据做测试数据。所有的压缩算法都支持Illumina平台测序的fastq格式数据。表中，我们对DSRC2 的三种运行模式都跑了一遍，三种模式分别对应速度最快、居中以及压缩率最好。从表中可以看出来，在FASTQZ、SCALCE、QUIP、DSRC2，BEETL这几个压缩方法里面，BEETL的压缩效率相比于别的几个压缩算法，压缩率有个很明显的不正常，这主要是因为BEETL采用的是ASCII码格式存储。BEETL算法提供了两种存储格式，另一种存储格式的压缩率很好，大约是ASCII码格式的1/8，但是解压缩需要还原成ASCII码格式后才可以。同时BEETL压缩解压缩以后，我们发现它是元数据有损的无损压缩。除了元数据部分以外，序列和质量分数都是无损。其他的几种算法都是无损压缩，经过压缩和解压缩还原以后的数据和原数据完全一致。同时，在实际运行中，发现BEETL对内存消耗很大。SCALCE算法结果也对整个数据进行了修改，元数据部分被改变了，同时解压后的文件是以序列和对应的质量分数为一组的乱序排列，尽管数据本身除了元数据部分以外并没有丢失。几种算法中FASTQZ压缩效果最好，可能是因为FASTQZ中采用了很高效的混合有限上下文模式，并且在采用这个模式之前还对数据进行了处理，把所有的”N” 都传递给质量分数进行压缩，减少了对数据N的压缩，而这有效地提高了压缩效率。

表2 FASTQ格式数据压缩算法的压缩率结果图

表3是针对Fasta格式进行压缩的压缩率实验结果。BEETL因为同样的原因所以压缩率看起来比较差。在剩下来的两种方法中，我们可以明显看出MFCompress的压缩效果要好于Deliminate方法。可能是因为MFCompress方法本身也是在应用混合的有限上下文压缩之前先对数据进行了处理，有效减少了要压缩的数据量。而无论是混合有限上下文模型还是数据预处理，都对压缩率的提高带来了极大的帮助。MFCompress压缩的三种模式的压缩率不同也证明了阶数越高，压缩率越好。

表3 FASTA格式数据压缩的压缩率结果

我们还选取DSRC2来比较有损压缩和无损压缩的区别。如表4所示。相比于无损压缩，有损压缩的压缩率表现会好很多。

表4 DSRC2的有损和无损压缩结果

在某些场景下，比如数据只是存储在本地，同时需要不断地压缩和解压缩的时候，压缩时间比压缩率更重要。在压缩率上，表5和表6展示了各个压缩方案的压缩时间结果图。时间是采用Linux命令time直接得到得的real时间，单位为秒。

表5 FASTQ格式数据压缩时间

表6 FASTA格式数据压缩时间

从表5 中我们可以看到DSRC2的压缩耗时明显小于别的几个压缩算法，相比于通用压缩算法gzip，DSRC2所耗时间也并不逊色。表6中可以看到这几种压缩算法的耗时都很长。尤其是BEETL方法，可能是因为BEETL是基于BWT算法改进的压缩方法，BWT 算法耗时很长，内存消耗也很大。

总体上而言，如果注重压缩率，不太重视压缩和解压缩时间，比如数据需要进行传输，或者数据压缩解压缩的次数比较少，针对fastq格式压缩，可以考虑FASTQZ方法或者BEETL方法，具体的可以根据需求来选择。如果注重压缩时间，可以采用DSRC2方法。在对FASTA格式进行压缩，可以考虑MFCompress的不同参数来选择偏向压缩时间还是压缩率。或者也可以采用BEETL方法。

4 结 语

针对高通量测序技术产生的fastq格式短序列生物数据，计算机学者们开发了很多相应的压缩工具并取得了一定的成果。本文对现有的压缩方法进行了介绍并比较了它们的压缩性能。生物学者可以根据各自的需求选择相应的压缩工具。但是可以看出现有的压缩工具仍然存在着问题，计算机资源和数据存储量之间的不匹配问题也仍然存在。下面是对这些问题的阐述以及可能的解决途径：

(1) 压缩和解压缩时间问题。目前的压缩工具在压缩/解压缩时间上耗时过大，并行处理可以很好地提高压缩速率。但是仍然需要计算机学者们提供一个速度快，压缩率高的压缩算法。

(2) 发展性问题。如同sanger测序发展到高通量测序技术，测序技术在不断更新，第三代测序采用了基于纳米孔的单分子测序方法，这将会产生更优质的短片段序列，压缩算法也需要不断适应生物测序技术的发展。

(3) 计算机资源消耗问题。数据的增大对计算资源带来了挑战，在数据压缩状态进行数据的查询、提取和检索分析将极大的缓解对资源的压力。

测序技术推动着生物医学领域的发展，是生命科学领域内一项重要的技术变革。生物数据压缩的发展能很好地解决计算机资源和生物数据量不匹配的问题，减少存储和传输的消耗，推动生物医学领域的进一步发展。但是现在面临的挑战仍然很多，仍需要更好的短序列压缩方法，来解决这些问题。

[1] Salomon D.Data Compression:The Complete Reference[M].Springer-Verlag New York,Inc.,2000.

[2] Burrows M.A block-sorting lossless data compression algorithm[J].Technical Report Digital Src Research Report,1994,57(4):425.

[3] Jones D C,Ruzzo W L,Xinxia P,et al.Compression of next-generation sequencing reads aided by highly efficient de novo assembly[J].Nucleic Acids Research,2012,40(22):2364-2367.

[4] Anirban Dutta,Mohammed Monzoorul Haque,Tungadri Bose,et al.FQC:A novel approach for efficient compression,archival,and dissemination of fastq datasets[J].Journal of Bioinformatics & Computational Biology,2015,13(3).

[5] Bonfield J K,Mahoney M V.Compression of FASTQ and SAM Format Sequencing Data[J].Plos One,2013,8(3):1453-1456.

[6] Pinho A J,Diogo P.MFCompress:a compression tool for FASTA and multi-FASTA data[J].Bioinformatics,2014,30(1):117-118.

[7] Cox A J,Bauer M J,Jakobi T,et al.Large-scale compression of genomic sequence databases with the Burrows-Wheeler transform[J].Bioinformatics,2012,28(11):1415-9.

[9] Deorowicz S,Grabowski S.Compression of DNA sequence reads in FASTQ format[J].Bioinformatics,2011,27(6):860-862.

[10] Grassi E,Di Gregorio F,Molineris I.KungFQ:A simple and powerful approach to compress fastq files[J].IEEE/ACM Transactions on Computational Biology and Bioinformatics (TCBB),2012,9(6):1837-1842.

[11] Zhang Y,Li L,Xiao J,et al.FQZip:Lossless Reference-Based Compression of Next Generation Sequencing Data in FASTQ Format[C]//Proceedings of the 18th Asia Pacific Symposium on Intelligent and Evolutionary Systems-Volume 2.Springer International Publishing,2015:127-135.

[12] Zhang Y,Li L,Yang Y,et al.Light-weight reference-based compression of FASTQ data[J].Bmc Bioinformatics,2015,16(1):1-8.

[13] Nicolae M,Pathak S,Rajasekaran S.LFQC:a lossless compression algorithm for FASTQ files[J].Bioinformatics,2014,31(20).

[14] Popitsch N,von Haeseler A.NGC:lossless and lossy compression of aligned high-throughput sequencing data[J].Nucleic acids research,2013,41(1):e27-e27.

[15] Shkarin D.PPM:One step to practicality[C]//Data Compression Conference,2002.Proceedings.DCC 2002.IEEE,2002:202-211.

[16] Ochoa I,Asnani H,Bharadia D,et al.QualComp:a new lossy compressor for quality scores based on rate distortion theory[J].Bmc Bioinformatics,2013,14(1):1-16.

[17] Wan R,Anh V N,Asai K.Transformations for the compression of FASTQ quality scores of next-generation sequencing data[J].Bioinformatics,2012,28(5):628-635.

[18] Malysa G,Hernaez M,Ochoa I,et al.QVZ:lossy compression of quality values[J].Bioinformatics,2015,31(19):3122.

[19] Hach F,Numanagic I,Alkan C,et al.SCALCE:boosting sequence compression algorithms using locally consistent encoding[J].Bioinformatics,2012,28(23):3051-3057.

[20] 曾培龙.基于reads引导的基因组序列拼接[D].哈尔滨工业大学,2012.

[21] Kozanitis C,Saunders C,Kruglyak S,et al.Compressing Genomic Sequence Fragments Using SlimGene[J].Journal of Computational Biology,2011,18(3):310-324.

RESEARCH ON DATA COMPRESSION OF SHORT-SEQUENCE BIOLOGICAL DATA BASED ON NEXT-GENERATION SEQUENCING

Meng Qian

(SchoolofComputerScience,FudanUniversty,Shanghai200433,China)

Due to the development of next-generation sequencing technology (NGS), the rapid growth of sequential data has brought a heavy pressure to data storage and transmission. The data compression technique is an important method to solve this problem, but traditional compression methods do not exploit the characteristics of the data well. Therefore, scholars begin to focus on the compression algorithm which is the special one for NGS data. In this paper, we present a comprehensive summary of compression algorithms for the Fastq and Fasta data obtained from NGS. We introduce the features of Fastq and Fasta, and summarize the commonly used methods of sequential data compression. Then we evaluate these representative compression tools through tests on several data sets from various scales, species and sequencing platforms, in order to compare the compression performance and validate the characteristics so that they can support researchers as a guide for algorithm selection and improvement. Finally, some problems and the trends of short-sequence data compression algorithms are also proposed in this paper.

Data compression Short-sequence data compression Next-generation sequencing

2016-03-05。孟倩，硕士生，主研领域：生物信息学。

TP391

A

10.3969/j.issn.1000-386x.2017.04.005