CN103995988B

CN103995988B - 一种高通量dna测序质量分数无损压缩系统及压缩方法

Info

Publication number: CN103995988B
Application number: CN201410240933.8A
Authority: CN
Inventors: 周家锐; 华韵之; 纪震; 朱泽轩; 曾启明
Original assignee: Individual
Current assignee: Individual
Priority date: 2014-05-30
Filing date: 2014-05-30
Publication date: 2017-02-01
Anticipated expiration: 2034-05-30
Also published as: CN103995988A; WO2015180203A1

Abstract

本发明公开一种高通量DNA测序质量分数无损压缩系统及压缩方法，所述方法包括：A、预先基于文化基因算法构造质量分数压缩码本；B、接收输入的原始高通量DNA测序质量分数数据，针对数据中每个原始质量分数序列在质量分数压缩码本中搜索与其最相似的编码矢量；C、利用所搜索到的编码矢量对原始质量分数序列进行压缩。本发明通过对质量分数数据进行整体设计得到压缩码本，并借助文化基因算法优化，从而实现最佳压缩编码性能。其整体压缩率显著优于现有方法。另外，本发明的每个寻优个体表示单一编码矢量以及采用多模优化方式有效提升了码本设计效率。同时码本设计与压缩/解压缩过程相分离大大减少了运算时间。

Description

一种高通量DNA测序质量分数无损压缩系统及压缩方法

技术领域

本发明涉及数据压缩领域，尤其涉及一种高通量DNA测序质量分数无损压缩系统及压缩方法。

背景技术

DNA序列数据由针对DNA物质的测序技术采集获得，是遗传学、基因组学、生物信息学、医学等诸多领域的基础研究对象，具有重要科学价值与实际意义。随着新一代高通量测序技术 (Next-generation Sequencing, NGS) 日益成熟并大量使用，获取DNA数据所需时间有效降低，成本显著下降。但另一方面，其所产生的DNA数据量也在急剧增长，从而对现有存储与传输技术造成了巨大压力。亟须设计具针对性的压缩方法。

NGS高通量测序所获取的DNA数据常以FASTQ格式存储。与传统测序数据不同，FASTQ由一连串的短读记录 (Reads) 构成，每个短读包含三部分内容：(1). 元数据(Metadata)，用于描述短读名、测序平台等信息；(2). DNA碱基序列 (NucleotideSequence)，用于记录在当前短读中所获得的DNA片段；(3). 质量分数 (Quality Scores)，用于表示所对应DNA碱基序列中各符号测定的可信程度。在同一条短读记录内，其DNA碱基序列长度与质量分数序列长度是一致的。

现有的高通量测序数据压缩算法，一般只着眼于其短读中DNA碱基序列的压缩，而忽略了其它两个部分。对于元数据，因其整体相似度较高，仅使用差异编码即可获得较好的压缩结果。但对于质量分数部分，则需设计更具针对性的编码方法。其原因在于：(1). 质量分数与测序仪器、对应碱基序列等因素相关，其数据间差异度较高；(2). 与DNA碱基序列仅含A、T、G、C四种符号不同，质量分数往往包含数十种不同的字符，压缩难度更高；(3). 质量分数与DNA碱基序列长度相同，所占用的数据大小一致。

现有算法一般使用常见的熵编码方法，如哈夫曼编码 (Huffman Encoding)、游程编码 (Run-length Encoding, RLE) 等对高通量DNA测序质量分数进行无损压缩。而另外一些算法如Quip等，则使用高阶马尔科夫模型 (High-order Markov Chain) 对其进行预测编码。对于传统的熵编码压缩算法，由于其主要设计用于处理普通字符序列，并未考虑质量分数的独有数据特点，导致压缩性能不佳。在极端情况下，甚至出现编码后数据量反而有所增长的情况。而基于高阶马尔科夫模型的预测编码算法，一方面，其建模需统计全序列上各符号的出现频率，耗时较长。另一方面，预测模型所占存储体积较大，不适用于压缩较小的高通量测序数据。此外，模型的预测准确率与输入数据有着较大关联，对某些序列压缩率较低，算法鲁棒性能不佳。

因此，现有技术还有待于改进和发展。

发明内容

鉴于上述现有技术的不足，本发明的目的在于提供一种高通量DNA测序质量分数无损压缩系统及压缩方法，旨在解决目前高通量DNA测序数据压缩算法对质量分数数据针对性不强，压缩效果不理想的问题。

本发明的技术方案如下：

一种高通量DNA测序质量分数无损压缩方法，其中，所述方法包括以下步骤：

A、预先基于文化基因算法构造质量分数压缩码本；

B、接收输入的原始高通量DNA测序质量分数数据，针对数据中每个原始质量分数序列在质量分数压缩码本中搜索与其最相似的编码矢量；

C、利用所搜索到的编码矢量对相应的原始质量分数序列进行压缩。

所述的高通量DNA测序质量分数无损压缩方法，其中，所述原始高通量DNA测序质量分数数据为ASCII码编码的FASTQ格式。

所述的高通量DNA测序质量分数无损压缩方法，其中，所述步骤A具体为：

A1、设定质量分数压缩码本的大小M 以及编码矢量长度N，统计待输入的原始高通量DNA测序质量分数数据中的符号种类形成符号集合，并据此设置搜索范围；

A2、在搜索范围内随机构造M 个候选解长度为N 的寻优个体，形成进化种群，设定文化基因算法迭代次数为K，初始化迭代计数器k = 1；

A3、在每次迭代时，计算进化种群中每个寻优个体的适应度函数值；

A4、在计算所有寻优个体的适应度函数值后，使用适应度共享技术计算各寻优个体的共享适应度函数值；

A5、基于各寻优个体的共享适应度函数值，使用文化基因算法优化进化种群。

A6、更新迭代计数器k = k+1。若k<K，则返回步骤A3，否则执行步骤A7；

A7、将最终获得的进化种群中各寻优个体映射为各编码矢量，从而构成质量分数压缩码本。

所述的高通量DNA测序质量分数无损压缩方法，其中，所述步骤A3具体为：

A31、按照预定映射关系将寻优个体转换为编码矢量

A32、将编码矢量与原始高通量DNA测序质量分数数据中每个原始质量分数序列进行匹配，计算匹配编码后的数据体积；

A33、将该数据体积作为当前寻优个体的适应度函数值。

所述的高通量DNA测序质量分数无损压缩方法，其中，所述最相似的编码矢量指编辑距离最小的编码矢量。

所述的高通量DNA测序质量分数无损压缩方法，其中，所述步骤C具体为：

通过动态规划对原始质量分数序列和其最相似的编码矢量进行差异匹配，实现压缩编码形成压缩数据。

所述的高通量DNA测序质量分数无损压缩方法，其中，所述方法还包括：

D、利用所有原始质量分数序列的压缩数据构成数据集合，并将其与质量分数压缩码本作为系统输出结果。

所述的高通量DNA测序质量分数无损压缩方法，其中，解压缩时，将所述数据集合中的每个压缩数据根据所述质量分数压缩码本中的编码矢量恢复成原始质量分数序列，从而得到原始高通量DNA测序质量分数数据。

一种高通量DNA测序质量分数无损压缩系统，其中，所述系统包括：

压缩码本设计模块，用于基于文化基因算法构造质量分数压缩码本；

质量分数压缩模块，用于接收输入的原始高通量DNA测序质量分数数据，针对数据中每个原始质量分数序列在质量分数压缩码本中搜索与其最相似的编码矢量；并利用所搜索到的编码矢量对相应的原始质量分数序列进行压缩；

所述的高通量DNA测序质量分数无损压缩系统，其中，所述系统还包括：

质量分数解压模块，用于将所述数据集合中的每个压缩数据根据所述质量分数压缩码本中的编码矢量恢复成原始质量分数序列，从而得到原始高通量DNA测序质量分数数据。

有益效果：本发明提供一种高通量DNA测序质量分数无损压缩系统及压缩方法，本发明的压缩码本针对输入的NGS质量分数数据进行整体设计，并使用高效的文化基因算法予以优化。从而可获得最佳的压缩编码性能。使得本系统具有显著优于现有方法的整体压缩率，且在各数据文件上都保持了较好的鲁棒性能。另外，本发明在文化基因算法中，使用每个寻优个体表示单一的编码矢量，并以多模优化方式构造整个压缩码本。从而有效提升了码本设计效率。此外，码本设计过程与压缩、解压过程相分离，可使用离线构造的码本，压缩多个不同的质量分数数据文件，从而大幅度减少运算时间。

附图说明

图1为本发明具体实施例中高通量DNA测序质量分数无损压缩方法流程图。

图2为本发明基于码本的高通量DNA测序质量分数序列压缩过程示意图。

图3为本发明具体实施例中使用编码矢量对质量分数进行压缩编码的示意图。

图4为本发明基于文化基因算法构造质量分数压缩码本的算法示意图。

图5为图1中步骤S100的具体方法流程图。

图6为图5中步骤S130的具体方法流程图。

图7为本发明具体实施例中高通量DNA测序质量分数序列压缩系统原理框图。

图8为本发明高通量DNA测序质量分数序列压缩系统的工作示意图。

具体实施方式

本发明提供一种高通量DNA测序质量分数无损压缩系统及压缩方法，为使本发明的目的、技术方案及效果更加清楚、明确，以下对本发明进一步详细说明。应当理解，此处所描述的具体实施例仅仅用以解释本发明，并不用于限定本发明。

本发明提供的高通量DNA测序质量分数无损压缩方法是基于码本(Codebook)的编码方法来压缩NGS质量分数，如图1所示，其包括以下步骤：

S100、预先基于文化基因算法构造质量分数压缩码本。具体是基于文化基因算法（Memetic Algorithms, MAs）的多模优化(Multimodal Optimization)技术来设计质量分数压缩码本。

S200、接收输入的原始高通量DNA测序质量分数数据，针对数据中每个原始质量分数序列在质量分数压缩码本中搜索与其最相似的编码矢量。其中，所述原始高通量DNA测序质量分数数据为ASCII码编码的FASTQ格式，所述的最相似的编码矢量指编辑距离（EditDistance）最小的编码矢量(Code Vector)。

S300、利用所搜索到的编码矢量对相应的原始质量分数序列进行压缩。

其中，所述步骤S300具体为：通过动态规划对原始质量分数序列和其最相似的编码矢量进行差异匹配，实现压缩编码形成压缩数据。

如图2所示的是本发明的基于码本的高通量DNA测序质量分数序列压缩过程，对于某一输入短读记录中的质量分数序列，本发明会搜索质量分数压缩码本中与之最相似的编码矢量，使用这一编码矢量的序号及其与原始质量分数序列间的符号差异作为压缩编码。

如图3所示的具体实例，对于输入短读记录中的质量分数序列Q = “CCCGFF”，在质量分数压缩码本={C ₁,C ₂,…,C _M}中存在与之最相似的编码矢量C _m = “CCGHFFC”。则此质量分数序列可被编码为{m, Q ^*}，其中Q ^*表示Q与C _m间的符号差异：

其中U表示符号相同 (Unchanged)，I表示插入 (Insertion，以“”标注)，D表示删除 (Deletion，以“−”标注)，S表示符号替换 (Substitution)。对于插入与替换情况，原质量分数符号也需一并记录 (例如第三个符号处插入的原质量分数“C”)。质量分数序列与编码矢量的这一差异匹配过程可通过动态规划 (Dynamic Programming, DP) 快速完成。

在FASTQ文件中，原始质量分数序列Q使用ASCII码 (8bits/字符) 存储每个符号，其数据大小为L _O = 8 × |Q|。当进行编码后，Q ^*中包含4种符号差异类型{U, I, D, S}，每个差异占用2bits存储空间。因此编码后的质量分数序列大小为：

；

其中M为压缩码本中编码矢量的总个数，T为编码后插入与替换情况下需存储的原质量分数符号个数。在绝大多数情况下，编码后的数据体积L _C将远小于原始数据体积L _O，从而获得压缩效果。而编码矢量C _m与Q越相似，则编码后|Q ^*|与T数值越小，数据体积L _C也将越小。亦即，质量分数压缩码本对质量分数数据的代表性越好，压缩率也将越高。

因此，质量分数压缩码本设计是本发明的另一关键，本发明基于文化基因算法构造质量分数压缩码本的过程如图4，其具体实施步骤如图5所示，即图5是所述步骤S100的具体展开，结合图4对图5中内容进行阐述：

S110、设定质量分数压缩码本的大小M以及编码矢量长度N，统计待输入的原始高通量DNA测序质量分数数据中的符号种类形成符号集合，例如集合S，并据此设置搜索范围，该搜索范围为R = [0, |S|]^N。例如，若输入的质量分数由“B”、“C”、“G”、“I”及“F”五种符号组成，则有S = [B, C, G, I, F]，R = [0, 5]^N。

S120、（优化前）在搜索范围R内随机构造M个候选解(Candidate Solution)长度为N的寻优个体，形成进化种群ps，设定文化基因算法迭代次数为K，初始化迭代计数器k = 1。

S130、在每次迭代时，计算进化种群ps中每个寻优个体的适应度函数值。

其计算过程如图6所示，所述步骤S130进一步包括：

S131、按照预定映射关系将寻优个体转换为编码矢量。例如设若输入的第m个体候选解为X _m = [x ₁, x ₂, ..., x _N]，其在各维度上均为R范围内的连续实数值。首先将X _m转换为离散符号序列编码矢量C _m = "s ₁ s ₂... s _N"，其中有映射关系：。

S132、将编码矢量与原始高通量DNA测序质量分数数据中每个原始质量分数序列进行匹配，计算匹配编码后的数据体积。即将C _m与原始质量分数序列集合 = {Q ₁, Q ₂,…, Q _P}中的每个序列进行匹配，计算编码后的数据体积总和为：

；其中P为输入的原始高通量DNA测序质量分数序列的总数，L _C(C _m, Q _p)表示编码矢量C _m在质量分数序列Q _p上匹配编码后的体积。其计算可采取如图2所示的方法进行。

S133、将该数据体积作为当前寻优个体的适应度函数值。即设。适应度函数值越小，则表示当前个体对输入质量分数序列集合的代表性越好，则其所构成的码本可获得更佳的整体压缩率。

S140、在计算所有寻优个体的适应度函数值后，使用适应度共享（FitnessSharing）技术计算各寻优个体的共享适应度函数值。

，其中：

其中参数ε为小生境半径 (Niching Radius)，参数α用于控制共享适应度函数的形态，距离d _i,j计算公式如下：

其中dist(X _i, X _j)表示寻优个体X _i与X _j间的马氏距离 (Manhattan Distance)。若两个个体位于过分相近的寻优空间范围，则其共享适应度函数值将显著变差，导致个体被驱散至不同的搜索区域。通过使用适应度共享（Fitness Sharing）技术，可保证优化完成后，质量分数压缩码本中各编码矢量间的冗余度最小。

S150、基于各寻优个体的共享适应度函数值，使用文化基因算法优化进化种群。

其中文化基因算法常用的为差分进化 (Differential Evolution, DE) 与Davies, Swann, and Campey with Gram-Schmidt Orthogonalization (DSCG) 优化方法的混合算法。

S160、更新迭代计数器k = k+1。若k<K，则返回步骤S130，否则执行步骤S170。

S170、将最终获得的进化种群中各寻优个体映射为各编码矢量，从而构成质量分数压缩码本。

将最终获得的进化种群ps中各寻优个体X _m，通过与图3中步骤S131的相同方法映射为编码矢量C _m，从而构成压缩码本= {C ₁, C ₂, …, C _M}输出。

进一步地，所述的高通量DNA测序质量分数无损压缩方法还包括：

利用所有原始质量分数序列的压缩数据构成数据集合，并将其与质量分数压缩码本作为系统输出结果。

另外，解压缩时，将所述数据集合中的每个压缩数据根据所述质量分数压缩码本中的编码矢量恢复成原始质量分数序列，从而得到原始高通量DNA测序质量分数数据。

如图7所示的高通量DNA测序质量分数无损压缩系统，其中，所述系统包括：

压缩码本设计模块100，用于基于文化基因算法构造质量分数压缩码本；

质量分数压缩模块200，用于接收输入的原始高通量DNA测序质量分数数据，针对数据中每个原始质量分数序列在质量分数压缩码本中搜索与其最相似的编码矢量；并利用所搜索到的编码矢量对相应的原始质量分数序列进行压缩。即主要用于根据压缩码本设计模块设计的压缩码本，对输入质量分数数据进行无损压缩编码。

另外，所述系统还包括：

质量分数解压模块300，用于将所述数据集合中的每个压缩数据根据所述质量分数压缩码本中的编码矢量恢复成原始质量分数序列，从而得到原始高通量DNA测序质量分数数据。即用于对压缩后的数据文件进行解压恢复操作。

该高通量DNA测序质量分数无损压缩系统的大致工作过程如图8所示。S1、数据输入。

S2、输入的是否为原始质量分数序列集，即判断输入数据是否为原始质量分数序列集，若是，则执行步骤S3，若否，则输出给质量分数解压模块300执行步骤S5。

S3、输入数据是否包含压缩码本，若是，则将数据输出给质量分数压缩码本200执行步骤S4，若否，则将数据输出给压缩码本设计模块100完成压缩码本设计，并在之后将数据和设计的压缩码本输出给质量分数压缩模块200执行步骤S4。

S4、输出压缩码本与压缩后质量分数数据。即经质量分数压缩模块处理得到压缩后质量分数数据，将其与压缩码本一起输出。

S5、输出解压恢复的原始质量分数数据集。经质量分数解压模块对压缩数据的解压处理得到原始质量分数数据集。

上述系统工作流程可进一步表述为：对于输入的原始高通量DNA测序质量分数序列集合，首先使用码本设计模块建立压缩码本，其过程如图4所示。而后，对于中的每个序列Q _p，选择与之最相似的编码矢量C _p∈对其进行压缩：

其编码方法如图2所示，从而形成压缩数据{m _p, Q _p ^*}。重复此过程直至中所有质量分数序列都已被压缩编码，从而构成压缩后的数据集合 = {{m ₁, Q ₁ ^*}, {m ₂,Q ₂ ^*}, …, {m _P, Q _P ^*}}。最后，将与作为系统的输出结果。

在进行解压缩时，将中的每个编码数据{m _p, Q _p ^*}，根据输入码本中的编码矢量C _p恢复其原始质量分数序列Q _p，从而还原出原质量分数序列集合即可。

本发明提供一种高通量DNA测序质量分数无损压缩系统及压缩方法，本发明的压缩码本针对输入的NGS质量分数数据进行整体设计，并使用高效的文化基因算法予以优化。从而可获得最佳的压缩编码性能。使得本系统具有显著优于现有方法的整体压缩率，且在各数据文件上都保持了较好的鲁棒性能。另外，本发明在文化基因算法中，使用每个寻优个体表示单一的编码矢量，并以多模优化方式构造整个压缩码本。从而有效提升了码本设计效率。此外，码本设计过程与压缩、解压过程相分离，可使用离线构造的码本，压缩多个不同的质量分数数据文件，从而大幅度减少运算时间。

应当理解的是，本发明的应用不限于上述的举例，对本领域普通技术人员来说，可以根据上述说明加以改进或变换，所有这些改进和变换都应属于本发明所附权利要求的保护范围。

Claims

1.一种高通量DNA测序质量分数无损压缩方法，其特征在于，所述方法包括以下步骤：

A、预先基于文化基因算法构造质量分数压缩码本；

C、利用所搜索到的编码矢量对相应的原始质量分数序列进行压缩；

所述步骤A具体为：

A1、设定质量分数压缩码本的大小M 以及编码矢量长度N ，统计待输入的原始高通量DNA测序质量分数数据中的符号种类形成符号集合，并据此设置搜索范围；

A2、在搜索范围内随机构造M 个候选解长度为N 的寻优个体，形成进化种群，设定文化基因算法迭代次数为K ，初始化迭代计数器k = 1；

A5、基于各寻优个体的共享适应度函数值，使用文化基因算法优化进化种群；

A6、更新迭代计数器k = k+1，若k<K ，则返回步骤A3，否则执行步骤A7；

2.根据权利要求1所述的高通量DNA测序质量分数无损压缩方法，其特征在于，所述原始高通量DNA测序质量分数数据为ASCII码编码的FASTQ格式。

3.根据权利要求1所述的高通量DNA测序质量分数无损压缩方法，其特征在于，所述步骤A3具体为：

A31、按照预定映射关系将寻优个体转换为编码矢量

A33、将该数据体积作为当前寻优个体的适应度函数值。

4.根据权利要求3所述的高通量DNA测序质量分数无损压缩方法，其特征在于，所述最相似的编码矢量指编辑距离最小的编码矢量。

5.根据权利要求4所述的高通量DNA测序质量分数无损压缩方法，其特征在于，所述步骤C具体为：

6.根据权利要求5所述的高通量DNA测序质量分数无损压缩方法，其特征在于，所述方法还包括：

7.根据权利要求6所述的高通量DNA测序质量分数无损压缩方法，其特征在于，解压缩时，将所述数据集合中的每个压缩数据根据所述质量分数压缩码本中的编码矢量恢复成原始质量分数序列，从而得到原始高通量DNA测序质量分数数据。

8.一种高通量DNA测序质量分数无损压缩系统，其特征在于，所述系统包括：

A6、更新迭代计数器k = k+1，若k<K，则返回步骤A3，否则执行步骤A7；

A7、将最终获得的进化种群中各寻优个体映射为各编码矢量，从而构成质量分数压缩码本;

质量分数压缩模块，用于接收输入的原始高通量DNA测序质量分数数据，针对数据中每个原始质量分数序列在质量分数压缩码本中搜索与其最相似的编码矢量；并利用所搜索到的编码矢量对相应的原始质量分数序列进行压缩。

9.根据权利要求8所述的高通量DNA测序质量分数无损压缩系统，其特征在于，所述系统还包括：

质量分数解压模块，用于将数据集合中的每个压缩数据根据所述质量分数压缩码本中的编码矢量恢复成原始质量分数序列，从而得到原始高通量DNA测序质量分数数据。