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

摘要

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

说明书

技术领域

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

背景技术

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

NGS高通量测序所获取的DNA数据常以FASTQ各式存储。与传统测 序数据不同，FASTQ由一连串的短读记录(Reads)构成，每个短读包含三 部分内容：(1).元数据(Metadata)，用于描述短读名、测序平台等信息；(2). DNA碱基序列(Nucleotide Sequence)，用于记录在当前短读中所获得的 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 格式，所述的最相似的编码矢量指编辑距离(Edit Distance)最小的编码矢 量(Code Vector)。

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

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

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

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

Q＝C C C  G - F F -

C_m＝C C ∧  G H F F C

Q^*＝U U (I,"C") U D U U D

其中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与原始质量分 数序列集合中的每个序列进行匹配，计算编码后的数 据体积总和为：

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

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

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

f_s(X_m)＝f(X_m)×τ_i，其中：

${\tau }_i=\underset{j\in \left|\mathrm{ps}\right|,j\ne i}{\Sigma }{(1-\frac{d_{i,j}}{ϵ})}^{\alpha }$

其中参数ε为小生境半径(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，从而构成压缩码本输 出。

进一步地，所述的高通量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，选择与之最相似的编码矢量对其进行压缩：

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

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

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

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