技术领域
[0001] 本
发明是给出一种面向多条FASTA格式基因序列的无损压缩方法，尤其涉及
生物信息学基因测序
数据压缩领域。
背景技术
[0002] 自从2001年发布第一条人类基因图谱以来，随着高通量测序技术的出现，基因测序成本不断降低，测序速度不断提高。基因技术在生物研究，
疾病预防，精准医疗等方面发挥着越来越重要的作用。众多国家和组织纷纷启动基因工程计划，其中著名的有1000Genome Project，它测试了1092个人的全序列基因。接着，人们又不断宣布开展更为宏大的测试计划，比如UK10K project,the Personal Genomes Project,and the Million Veteran Project(MVP)，这些计划测序的基因数分别为10K,100K,1M。除此之外，类似的针对其他种群的基因测序计划也有很多，比如最著名的是1001Genome Project，大概1000多种拟南芥基因被测序。另外，还有theEarth Biogenome Project(EBP),计划利用10年的时间测序世界上所有生物基因序列。伴随着基因科学发展的是庞大的基因数据给存储和传输带来了巨大的压
力。最新论文显示，基因数据存储和传输的代价已经超过了测序的代价。如何以更高的效率来压缩基因数据，降低存储和数据迁移成本，在基因研究和应用中有着十分重要的作用。
[0003] 目前已经提出了很多种FASTA格式专用基因压缩方法，但更多的是单序列压缩方法，比如GRS,GReEn,RLZ,RLZ-OPT,ERGC,NRGC,HiRGC,SCCG等。单序列压缩方法只适用于一次压缩一条基因序列，不支持多条FASTA格式基因序列批量压缩。即如果存在m条待压缩序列，则单序列压缩方法需要执行m次。这种压缩方式不能充分利用待压缩序列间的相似性提升压缩率。而目前基因应用中，同时压缩多条序列甚至成千上万条基因序列用于基因数据的存储与迁移越来越普遍。面向多条基因序列的压缩方法能够一次压缩批量基因序列于一个压缩文件，而且在待压缩序列间能够实施二次匹配，充分利用待压缩序列间的相似性提高压缩率。但是截止目前在已发表的国内外论文中，宣称可支持多条FASTA格式基因序列批量压缩方法只有FRESCO,COGI,GDC2三种。FRESCO是2013年由Sebastian Wandelt和Ulf Leser在“IEEE/ACM Transactions on Computational Biology and Bioinformatics”提出的第一种可用于多条基因序列同时压缩的方法。但FRESCO是
有损压缩，它只支持{A,C,G,T,N}五种字符。而实际上，基因序列除了ACGT四种基本
碱基字符，尚未探测出的碱基位N之外，还包含基因序列标识信息，用来标识序列长度，测序日期等；大量其他特殊字符，比如R代表该碱基位可能是G或A；Y代表该碱基位可能是T或C等。另外，还常用小写字符来表示基因序列内部的重复
片段，因此基因序列内部还包含丰富的大小写字符信息。所以FRESCO在压缩之前需要对序列进行预处理，将所有非{A,C,G,T,N}字符之外的其他字符丢弃掉，这样会丢失部分信息，而FASTA格式基因序列任何基因信息的丢失都可能给生信人员带来不可预测的后果。COGI是2015年5月由Xiaojing Xie等在“IEEE/ACM Transactions on Computational Biology and Bioinformatics”提出的基因序列压缩方法，COGI的压缩原理也可以用于多条基因序列同时压缩。但首先COGI也是有损压缩，它只支持{a,c,g,t,n}五种字符。其次COGI在多条基因序列同时压缩时，压缩率会随着基因数量的增多而下降。因为COGI的实现原理决定了随着待压缩序列数量增多，经过异或与序列矩阵化后的二维比特矩阵会非常复杂，矩阵分割编码方法在这样的矩阵上压缩效果会急剧下降。在COGI论文中的实验部分，以KOREF_20090224为参考基因压缩YH基因，平均压缩率为244。然而，当YH,KOREF_20090224和KOREF_20090131同时压缩时，平均压缩率只有11.7。所以，COGI虽然能够实现多序列同时压缩，但压缩效果并不理想。GDC2是Sebastian Deorowicz等人于2015年6月在“Scientific Reports”提出的支持多条序列批量压缩的FASTA格式基因序列压缩方法。GDC2的优点在于首次实现了多条基因序列的无损压缩，不足的地方在于GDC2方法采用所有种类字符参与哈希值计算的线性探测构建哈希表的方法，哈希索引构建效率和匹配阶段的查询效率都很低，导致GDC2压缩速度非常慢。GDC2在二次匹配后的四元组表示法导致其基因压缩率也不高。而且GDC2内存消耗很大，在32GB RAM配置的
服务器上同时压缩超过110个人类基因时，会报内存分配失败的错误，导致压缩失败，因此GDC2的易用性很差。
[0004] 目前已公开的面向多条FASTA格式基因序列的压缩方法
专利有2019年7月5号公开的公开号为CN109979537A的专利“一种面向多条序列的基因序列数据压缩方法”，该方法通过将所有待压缩序列与参考序列按字符对齐然后异或的方式，将待压缩序列与参考序列相同的碱基字符变成’0’字符。然后将所有处理过的待压缩序列矩阵化后通过矩阵分割编码方法表示矩阵中的非’0’字符实现多条基因的压缩存储。该方法采用二维比特图模型来实现基因序列的压缩，与本发明采用信息提取与多次匹配实现多条基因序列无损压缩完全不同。该方法对待压缩序列与参考序列的相似性要求非常高。在待压缩序列与参考序列相似性不够时，异或操作会产生大量的非‘0’字符，影响压缩率，而且该方法也是一种有损压缩方法。所以，总体而言，目前国内外还没有一种有效地面向多条基因序列的无损压缩方法。
发明内容
[0005] 为解决上述技术问题，本发明所采用的技术方案如下：提供一种面向多条FASTA格式基因序列的无损压缩方法，而且压缩率、压缩时间和内存消耗可以由压缩者根据
硬件资源和待压缩基因序列数据量进行动态地调节。
[0006] 将大规模基因序列的压缩依次分为三个主要步骤：序列信息提取、序列信息匹配和序列信息编码。序列信息提取阶段提取所有基因序列的基因文件标识信息、基因文件行宽、碱基小写字符的长度和
位置、ACGT基本碱基信息、ACGT基本碱基字符之外的其他特殊字符信息。这些信息包括了所有待压缩序列的信息，所以在恢复阶段能够根据这些信息无损地恢复出原始序列。序列信息匹配阶段对所有待压缩序列第一阶段提取的小写字符信息和ACGT基本碱基信息与参考序列对应信息进行匹配，小写字符信息为一次匹配，基本碱基信息分为首次匹配和二次匹配。首次匹配目的在于去除待压缩序列相较于参考序列的冗余信息；二次匹配为所有待压缩序列基本碱基信息首次匹配结果的再次匹配，目的是去除待压缩序列间基本碱基信息的冗余信息。两次匹配过程结束后，综合运用游程编码，差分编码、静态
熵编码和PPMD编码，输出压缩文件。详细步骤如下：
[0007] 1)序列信息提取。序列信息提取采用
逐行扫描基因序列文件的方式，首先提取文件中的第一行信息，即是基因序列文件的标识信息，因为所有FASTA格式基因文件第一行均是标识信息，标识基因序列
染色体号，测序时间等信息。从第二行开始即为碱基字符信息，采用逐个字符检测方式，如果该字符C∈{A,C,G,T}中的一个，则将它存入一个新建的ACGT序列，我们称之为基本碱基序列 序列；如果字符是{a,c,g,t}中的一个，则首先将它转化为大写字符，然后加入 序列，并同时记录小写字符的位置和连续的小写字符的长度，以(position,length)二元组形式加入小写字符向量Vlow，这样在恢复阶段可以根据Vlow将所有大写字符恢复成小写字符。如果是N字符，则记录它的位置和连续N字符的长度，加入N字符向量VN，因为N字符一般是连续出现的；如果是n字符，则同样首先将它转
化成大写N，然后将它们的小写信息加入向量Vlow，N字符信息加入向量VN。如果是特殊字符，则将该特殊字符信息记录到特殊字符信息数组Vspe。最后，记录下该文件一行的字符个数即行宽，因为FASTA基因文件都是等距换行的，所以只要记录一个行宽即可。以此方式扫描到文件最后一行，序列信息提取即结束。序列信息提取结束后，将输出6个信息：基因序列标识信息identifier，基本碱基序列 序列，小写字符信息Vlow，N字符信息VN，特殊字符信息Vspe，行宽lineWidth。
[0008] 2)小写字符信息匹配。小写字符信息匹配是依次对第一阶段提取出的待压缩序列小写字符向量 与参考序列小写字符信息向量 进行匹配，目的是去除待压缩序列向量相较于参考序列向量 的冗余信息。小写字符信息匹配采取顺序搜索的方式，获取所有待压缩序列 与参考序列 的相同元素，然后将相同元素用参考序列 下标表示，以达到减小小写字符信息存储体积的目的。
[0009] 3)基本碱基序列首次匹配。基本碱基序列只包含A,C,G,T四种碱基字符。将复杂的FASTA文件转化成基本碱基序列进行匹配，可以大大降低匹配的复杂度，提升匹配的成功率和效率。首次匹配以参考文件 序列为参考，查找待压缩文件 序列与参考文件 序列匹配片段。当有多个匹配片段时，则选取最长匹配片段。然后将待压缩文件 序列中被匹配片段表示成(position,length)二元组序列，同时组合相邻的未匹配片段mismatched string，最终将待压缩文件 序列分解成(position,length,mismatched)三元组序列。
[0010] 4)基本碱基序列二次匹配。二次匹配是对所有待压缩序列在首次匹配时输出的三元组序列的再次匹配。二次匹配的元素是(position,length,mismatched)三元组，而不是首次匹配中的A,C,G,T基本碱基字符。而且首次匹配后，待压缩序列输出的三元组序列已远小于其基本碱基序列 序列。所以为了提高压缩率，二次匹配选择多条待压缩序列首次匹配结果作为二次匹配的参考序列，也不至于导致内存消耗过大。当存在多条参考序列时，待压缩序列需要在所有参考序列中找到最长匹配作为匹配结果。二次匹配结果中匹配成功以(sequence_id,position,length)三元组表示，未匹配成功则保留首次匹配结果(position,length,mismatched)三元组形式。在解压缩时通过检查三元组最后一个元素是数值还是碱基字符串来区分是首次匹配结果还是二次匹配结果。因为二次匹配结果三元组的最后一个元素length是数值，而首次匹配结果三元组最后一个元素mismatched是碱基字符串。
[0011] 参考序列数量的不同会导致压缩率不同，同时也会影响压缩时间和内存消耗。为了能让该压缩方法满足不同的压缩率和压缩时间要求，以及能够运行于不同配置的机器上。本方法定义了一个二次匹配百分比，即二次匹配中用于作为参考序列的待压缩序列数量占所有待压缩序列数量的百分比。该百分比可以由压缩者自行定义。
[0012] 5)匹配结果编码。对待压缩序列信息提取阶段提取出的序列标识信息identifier，行宽lineWidth，特殊字符信息Vspe，N字符信息VN；小写字符信息匹配步骤产生的匹配结果；及基本碱基序列二次匹配后产生的匹配结果进行编码存储。具体为：
[0013] 对于标识信息identifier组成的数组及行宽信息lineWidth组成的数组，采用游程编码，连续相同的标识信息及行宽信息只记录连续出现的次数即可；
[0014] 对于特殊字符信息Vspe，采用Huffman编码。
[0015] 对于N字符信息VN，采用增量编码。
[0016] 对于小写字符信息匹配步骤产生的匹配结果，采用游程编码与增量编码；
[0017] 对于基本碱基序列二次匹配后产生的匹配结果采用增量编码；
[0018] 最后对所有信息采用PPMD编码存储。
[0019] 6)解压缩。因为PPMD编码，游程编码、Huffman编码、增量编码都是可逆的，所以可以依次恢复出待压缩序列标识信息和行宽信息，特殊字符信息，N字符信息。对于小写字符信息的解压过程为，首先按照压缩阶段相同的方法提取出参考序列的小写字符向量 然后依据 与匹配结果恢复各待压缩序列小写字符信息 对于基本碱基序列的解压缩过程为，首先从待压缩文件中恢复出第一条三元组序列，然后以第一条三元组序列为参考序列恢复第二条三元组序列，再以第一条和第二条三元组序列为参考序列恢复第三条三元组序列，依此类推，直至恢复出所有的三元组序列，二次匹配结果恢复完成。然后所有恢复出来的三元组序列再结合参考序列基本碱基序列恢复出原始的基本碱基序列，直至所有的基本碱基序列恢复。最后一步综合基本碱基序列、小写字符信息、行宽、N字符信息、特殊字符信息，标识信息恢复出基因文件。压缩阶段的所有步骤都是可逆的，所以能够实现基因序列的无损压缩与恢复。
[0020] 本发明与
现有技术相比，其有益效果是：
[0021] 本发明所提出的基于信息提取与多次匹配的面向多条基因序列的无损压缩方法，相比于当前存在的同类方法，提高了基因压缩率和压缩速度，减小了系统内存消耗。在信息提取阶段，通过对参考序列和待压缩序列进行信息提取，将待压缩文件分成相似性各不相同的六部分。区别于部分同类方法的直接丢弃部分信息比如小写字符信息和特殊字符信息而导致无法做到无损恢复，本压缩方法能够实现无损恢复。在信息匹配阶段，对提取后的基本碱基序列进行匹配，因为基本碱基序列只包含ACGT四种基本碱基字符，所以降低了匹配的复杂度，提高了匹配的效率。区别于部分同类方法未经过信息提取的原始碱基序列所有字符参与匹配，匹配复杂度高，压缩效率低，本发明提升了压缩效率。本方法另一个特点是对信息提取之后的相似性最大的基本碱基序列B序列采用二次匹配方法，而且对二次匹配结果全部采用三元组形式表示，通过首次匹配结果的三元组与二次匹配结果的三元组最后一个元素不同数值类型来识别，区别于部分同类方法通过标识符来标识首次匹配结果与二次匹配结果的四元组表示方法，进一步提升了压缩率。图4和图5实验结果表明，与当前顶级的FASTA格式基因压缩方法对比，在110条基因序列同时压缩时，本发明不管从压缩率，还是从压缩速度和内存消耗，在不同长度的染色体上均远远优于同类优秀方法。
附图说明
[0022] 图1是本发明方法的压缩流程示意图
[0023] 图2是本发明在首次匹配与二次匹配中的索引构建与匹配查找示意图[0024] 图3是本发明的基本碱基序列首次匹配与二次匹配效果示意图
[0025] 图4是本发明在的实际压缩率和压缩速度对比图
[0026] 图5是本发明在实际内存消耗对比图
具体实施方式
[0027] 如图所示，一种面向多条FASTA格式基因序列的无损压缩方法依次分为以下步骤：
[0028] 1)对参考序列进行信息提取。参考序列只提取基本碱基字符序列和小写字符的位置信息，其他字符(换行符，N字符，特殊字符)信息丢弃。因为本方法并不压缩参考序列。将提取出的基本碱基字符序列存储为数组 将小写字符信息存储为向量Vlow，且它们都暂时存储于内存中，供第二阶段匹配。
[0029] 2)对待压缩序列进行信息提取。待压缩序列除了提取基本碱基字符序列和小写字符的位置信息，其他字符(换行符，N字符，特殊字符)也需要分别提取。提取出的小写字符信息存储为向量 基本碱基字符序列存储为数组 基因序列标识信息存储为字符串identifier1，N字符信息存储为向量VN，特殊字符信息Vspe，行宽lineWidth1。将identifier1、lineWidth1都暂时存储于内存中， 通过步骤3进行压缩， 通过步骤5进行压缩，identifier1、lineWidth1待所有待压缩序列信息提取完成后统一经过游程编码进行压缩。VN直接写入压缩文件中，Vspe通过静态熵编码后写入压缩文件中。
[0030] 3)小写字符信息匹配。将步骤2中提取出的 和步骤1中的Vlow进行相似匹配。对中的元素如果在Vlow中存在，只存储其在Vlow中的位置即可。对于连续出现在Vlow中的元素，记录其在Vlow中的起始位置和连续出现的长度。这样可以实现 压缩，减小小写字符信息存储体积。
[0031] 4)为参考序列基本碱基序列构建哈希索引表。对 按照图2的方法建立哈希索引表。创建 数组，数组元素全部初始化为-1。创建 数组，用于存储 表中的冲突元素。依次计算 中所有k-mer(即连续k个碱基字符)的hash值。假设第i个k-mer计算出的哈希值记录为vi，按照 方法将k-mer位置信息存储在 表中，将冲突k-mer位置存储在 表中。依此，通过 表和 表可以存储 中所有k-mer信息，在匹配时，也可以遍历 中所有k-mer信息。
[0032] 5)将步骤2中提取出的待压缩基本碱基序列与参考序列基本碱基序列进行匹配。从待压缩基本碱基序列第一个碱基字符开始，依次计算k-mer哈希值，然后查找步骤3中构建的 表中该哈希值下元素是否为-1。如果为-1，表明该k-mer在参考基本碱基序列中不存在，则待压缩序列下该k-mer的首个字符作为未匹配字符记录在mismatched字符串中。如果元素值不为-1，表明该k-mer在参考序列中存在，则遍历 表和 表查找最长匹配，记录下该最长匹配在参考序列中的位置和长度二元组(position,length)，同时组合相邻的未匹配片段mismatched字符串，将待压缩序列被匹配片段表示成(position,length，mismatched)三元组。待压缩序列跳过被匹配片段继续匹配，直至所有碱基字符记录在三元组中，待压缩基本碱基序列分解完成。将匹配结果三元组序列同时写入待压缩文件中和存储在内存中。如果待压缩序列只有一条，压缩结束；否则继续进行步骤6。
[0033] 6)对剩下的所有待压缩序列依次按照步骤2进行序列信息提取。第i个待压缩序列提取出的小写字符信息存储为向量 基本碱基字符序列存储为数组 基因序列标识信息存储为字符串identifieri，行宽存储为整数lineWidthi且它们都暂时存储于内存中。N字符信息以向量 形式直接存储于压缩文件中，特殊字符信息 经静态熵编码后存储于压缩文件中。对于内存中小写字符向量 按照步骤3匹配并将匹配结果写入压缩文件。对于内存中基本碱基字符序列，按照步骤4匹配，将匹配结果仅暂时存储于内存中，通过步骤8进行二次匹配。identifieri、lineWidthi待所有待压缩序列信息提取完成后统一经过游程编码进行压缩。
[0034] 7)为待压缩序列首次匹配结果构建哈希索引。二次匹配在所有暂时存储在内存中的待压缩序列首次匹配结果间进行，所以需要为待压缩序列首次匹配结果构建哈希索引。二次匹配哈希索引表创建与首次匹配哈希索引表创建相同，如图2所示。为每个用于二次匹配的参考序列创建1个 数组和 数组。 数组元素全部初化为-1。L数组，用于存储 表中的冲突元素。按照公式 计算每个
首次匹配结果(position,length,mismatched)三元组的哈希值.按照
的方法将三元组的位置依次存储在 数组中，冲突元素存储在
数组中。至此，用于二次匹配的所有参考三元组序列哈希索引构建完成。
[0035] 8)对第i个待压缩序列的基本碱基字符序列数组 进行二次匹配。基于步骤7构建的哈希索引，查找所有 在 至 中的最长匹配，以(sequence_id,position,length)三元组表示并依次写入压缩文件中。未匹配三元组，依次直接写入压缩文件中。至此，完成所有待压缩序列基本碱基序列的二次匹配和存储。本发明的基本碱基序列首次匹配与二次匹配效果示意图如图3所示。
[0036] 9)设置二次匹配百分比，控制二次匹配中作为参考序列的待压缩序列数量。在默认情况下，本方法二次匹配百分比为10％，即10％的待压缩序列的基本碱基序列 存储在内存中，并构建相应的哈希索引表供步骤8二次匹配使用。当压缩者自定义二次百分比为p时，假设待压缩序列数目为m，则第i(i<＝mp)个待压缩序列基本碱基序列 在 至中查找最长匹配；第i(i>mp)个待压缩序列基本碱基序列 在 至 中查找最长匹配。
[0037] 如上所述，尽管参照特定的优选
实施例已经表示和表述了本发明，但其不得解释为对本发明自身的限制。在不脱离所附
权利要求定义的本发明的精神和范围前提下，可对其在形式上和细节上作出各种变化。
