CN113268459A - 基于fastq基因大数据的批量分布式压缩方法 - Google Patents

Abstract

本发明公开了基于FASTQ基因大数据的批量分布式压缩方法，包括：S100部署Spark集群、Hadoop集群和Java环境：S200、选取同种群物种、同一条染色体的FASTA文件，作为所有输入的基因序列的参考序列；S300为读取到的碱基序列与标识符序列创建Spark任务，S500driver本地程序从HDFS拉取压缩后的碱基序列与标识符序列、质量分数序列数据到本地，使用bsc压缩组件作最后的压缩运算，以获取压缩结果。该方法主要针对于FASTQ格式文件，主要用于解决基因数据量过大，减小基因数据存储和传输成本问题，并通过集群大大缩短大量数据的压缩时间。经验证性能将会较大幅度超越通用压缩器。

Description

基于FASTQ基因大数据的批量分布式压缩方法

技术领域

本发明涉及大数据领域与生物数据压缩技术领域，特别是涉及一种基于FASTQ基因大数据的批量分布式压缩方法。

背景技术

基因是DNA上有遗传效应的片断，和生命息息相关。基因数据因其重要的社会价值和科研价值受到国际社会的广泛重视。面对如此庞大并仍在不断增长的基因数据，以压缩算法对大量基因数据压缩可以减轻基因文件的存储和数据迁移压力，在基因研究和应用中显得尤为重要。FASTQ文件属于不同测序技术的直接产物，是一种存储了生物序列 (通常是核酸序列)以及相应的碱基测序质量评价的文本格式。这是序列分析中使用最广泛的格式，通常是定序器提供的格式。

基于该种格式的基因文件的压缩算法多运行于单机环境，但在面对大数据环境的情况下，单机算法也只有使用串行方式应对。中国专利申请CN110299187A-一种基于Hadoop的并行化基因数据压缩方法，虽然公开了能将大片段基于序列分成多个小片操作，并且可以针对多条基因序列分配多个计算节点对多个基因序列并行计算，以实现缩短基因压缩时间的目的。但是其在使用时，压缩时间仍然相对较长，难以达到预定要求。而中国专利CN111028897A-一种基于Hadoop的基因组索引构建的分布式并行计算方法，虽然公开了能够使用较少计算资源快速完成基因组索引构建，然后再进行数据清理以供给调用的分布式并行计算方法，在一定程度上提升了基因的压缩以及寻找时间。

上述两个现有技术，虽然其解决了部分串行模式资源占用率低的问题，但在面对TB 级的基因文件时，完成这些文件的压缩可能需要几十个小时的时间，耗时太久。基因序列数据过大，存储和传输成本高等问题。

发明内容

为了克服上述现有技术的不足，本发明提供了一种基于Spark平台，使用集群并行计算能力的优势，使用Spark平台框架规划运行基于参考序列的基因匹配压缩方法。将一个或多个基因文件在不同节点上分块并行读取并转化成二元组形式表示基因序列，并与参考序列匹配去除重复部分，保留不同的部分，提高压缩效率，通过多节点计算也可实现快速的压缩与高效稳定的存储，实现基因大数据并行压缩的方案，以解决基因序列数据过大、存储和传输成本高的问题的基于FASTQ基因大数据的批量分布式压缩方法。

本发明所采用的技术方案是：

基于FASTQ基因大数据的批量分布式压缩方法，其特征在于：包括如下步骤：

S100、部署Spark集群、Hadoop集群和Java环境：

将所有参与压缩的FASTQ文件上传至HDFS上；

在压缩环境中预先选取人类基因组序列hg38的参考序列信息和质量分数序列信息，并且放置到程序本地；

将Spark集群调整成Yarn standalone模式，在单独节点启动Driver；

S200、选取同种群物种、同一条染色体的FASTA文件，作为所有输入的基因序列的参考序列；从作为参考序列的FASTA文件应用哈希算法生成第一类参考文件，然后使用广播变量方法上传到每个计算节点的executor进程作为缓存，用于加速匹配进程；

S300、为读取到的碱基序列与标识符序列创建Spark任务，此Spark任务的输入信息包括FASTQ文件包含的碱基序列、标识符序列和参考序列，输出信息为包括<参考的序列号，匹配位置，匹配长度>三元组列表，然后转换成可以被二次识别的字符串格式输出至HDFS；并且：

S301、分片读取:

使用MapPartition算子对碱基序列与标识符序列进行分片，对每一分片逐行读取基因的ACGT字符信息、N字符信息、特殊字符信息和换行符信息，记录每一行的长度；

对于N字符信息，用<相对行首的偏移位置,出现长度>的二元组形式表示；特殊字符以<相对行首的偏移位置，特殊字符>的二元组表示；行尾的换行符使用整数表示其相对行首的偏移位置；

ACGT字符信息按照k-mer方式进行哈希编码，在内存中建立哈希索引数组和哈希冲突数组；根据步骤S200中的哈希算法，计算该ACGT信息在哈希索引数组与哈希冲突数组中对应的哈希值；再遍历参考序列，找到该哈希值对应的k-mer元素并记录其在参考序列中的位置号，以下简称参考位置号，然后将参考序列与待压缩序列的位置指针同时后移，连续比较直到ACGT信息不再匹配为止；

设置片段的匹配度门限，对低于该值或者违背规则的片段，以<哈希值，相对偏移量，字符串>三元组形式另外cache序列化到内存中，其中的相对偏移量表示该片段相对前一匹配片段的偏移量，字符串即该片段的文字表示；

对高于该值而且符合规则的片段，则生成一个格式为<参考位置号，匹配位置，匹配长度>的三元组，其中的参考位置号表示该片段在参考序列中对应的位置号，匹配位置表示该片段在待压缩FASTQ文件中的位置号，匹配长度即该片段的长度；一个待压缩 FASTQ文件的所有碱基序列信息被编程语言表示为一个Java ArrayList RDD返回，使用cache算子的MEMORY_SER方式缓存这些中间结果，以降低内存利用率；

S302、三元组哈希编码

使用MapPartitionWithindex算法获取所有分片的序号，再用filter算子过滤出需要参与二次匹配的第一类参考文件并记录其序号；然后，对S301筛选出的一次匹配三元组结果进行哈希编码；

S303、三个数组发布

将S302得到的三元组哈希编码以及第一类参考文件序号信息通过first算子拉取到 driver本地，对其封装成广播变量，然后发布到每个进行二次匹配的计算节点；

S304、二次匹配

在计算节点侧，根据得到的三元组哈希编码以及第一类参考文件序号信息，使用MapPartition算子对每一个三元组运用第一哈希算法，依次与第一类参考文件进行匹配，直到选取出最长匹配的元组，记录其出现的位置、长度以及对应的参考序列文件号；将所有匹配后的信息序列化为字符串，并输出至HDFS；

S400、为读取到的质量分数序列创建Spark任务，此任务的输入信息为以String表示的质量分数序列，输出信息为字节块，然后输出到HDFS，具体步骤为：

S401、聚类门限计算

使用MapPartition算子，从质量分数序列文件中读取出十万行数据，对其中的每一行记录，使用k-mer读取方式对每一个k-mer计算其质量估算总值，其中：所述质量估算总值为将行内每个字符与7按位或，然后进行移位计算；如果出现两行的质量估算值相同，则判定为它们的k-mer质量相同，统计所有符合此条件的k-mer的出现频次，并按出现频次从大到小排序，选取出现频次前70％的k-mer，找到频次最高的k-mer的质量估算总值乘以一个划分类别的系数，定义为第一聚类门限，以作为聚类划分依据；

S402、频次极值和字符寻找

对质量分数序列，计算字符ASCII码分布频次的极值，将极值对应的字符记为第一字符，以在映射阶段作编码对照；

S403、字节块分类

以行为单位，从质量分数序列的每行中按k-mer计算每个k-mer的质量估算值，将每行计算的质量估算总值与第一聚类门限作比较，低于则进入缓冲类别1，高于则进入缓冲类别2；这两种缓冲类别将以字节块的形式存储；使用MapPartitionWithIndex获取文件号，并且使用reduceByKey算子，利用shuffle将具有相同缓冲类别与相同文件号的块数据划分到一个分区，以提升并行度；

S404、字符块映射

使用MapPartition进行字节的映射，将连续的字节或值较大的字符块映射成较小的字节块；对于连续的或者与第一字符相近范围内的字节块，一并映射成更小的字节块；处理完成后，上述字节块将输出到HDFS。

S500、driver本地程序从HDFS拉取压缩后的碱基序列与标识符序列、质量分数序列数据到本地，使用bsc压缩组件作最后的压缩运算，从而输出最终压缩结果。

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

本发明是给出一种基于Spark的参考序列匹配压缩基因序列的压缩方法，该方法主要针对于FASTQ格式文件，主要用于解决基因数据量过大，减小基因数据存储和传输成本问题，并通过集群大大缩短大量数据的压缩时间。经验证性能将会较大幅度超越通用压缩器。

第一，事先从待压缩基因序列中选取并通过k-mer构建Hash表编码参考序列，并将参考序列存储为索引文件。

第二，启动Hadoop与Spark集群，将待压缩文件传到HDFS上，将，配置Spark任务运行脚本确定运行参数。

第三，在Spark任务中分两次读取待压缩序列的碱基序列、标识符与质量分数。对这两部分数据分别设计Spark算子流来完成。对碱基序列使用两次匹配的压缩方式，一次匹配使用同类染色体FASTA文件作为参考序列，二次匹配使用部分一次匹配结果进行匹配，使用哈希编码方式实现重新寻址匹配。

第四，在另一Spark任务中对于质量分数使用聚类并分数映射的方式缩减体量。

第五，使用广播变量优化Executor的内存配置，对质量分数使用重分区的方式提升映射并行度以提升计算效率。

使用分布式计算方式可以在读取单个基因实现并行化并提高效率，还可以实现多条基因序列并行处理，以实现处理压缩大批量基因文件的加速。

综上所述，本发明的批量分布式压缩方法，使用Spark作为大数据平台来进行分布式运算，它能将批量输入的基因封装到多个计算槽中并行操作，并可以分配多个计算节点来对这些不同的基因序列并行计算。将Spark的批处理并行化优势与基因的压缩算法结合起来，实现一个新的压缩方法，可以充分利用网络上的服务器资源，使用多台节点的计算能力大大缩短大量基因压缩的时间，以达到提高效率的目的。

附图说明

图1为基于FASTQ基因大数据的批量分布式压缩方法流程图；

图2为步骤S300的算法流程图；

图3为步骤S400的算法流程图。

具体实施方式

下面详细描述本发明的实施例，所述实施例的示例在附图中示出，其中自始至终相同或类似的标号表示相同或类似的元件或具有相同或类似功能的元件。下面通过参考附图描述的实施例是示例性的，仅用于解释本发明，而不能理解为对本发明的限制。

在本发明的描述中，需要理解的是，术语“中心”、“上”、“下”、“前”、“后”、“左”、“右”等指示的方位或位置关系为基于附图所示的方位或位置关系，仅是为了便于描述本发明和简化描述，而不是指示或暗示所指的组合或元件必须具有特定的方位、以特定的方位构造和操作，因此不能理解为对本发明的限制。另外，本发明实施例的描述过程中，所有图中的“上”、“下”、“前”、“后”、“左”、“右”等器件位置关系，均以图1为标准。

如图1所示，基于FASTQ基因大数据的批量分布式压缩方法，其特征在于：包括如下步骤：

S100、部署Spark集群、Hadoop集群和Java环境：

将所有参与压缩的FASTQ文件上传至HDFS上；

在压缩环境中预先选取人类基因组序列hg38的参考序列信息和质量分数序列信息，并且放置到程序本地；

将Spark集群调整成Yarn standalone模式，在单独节点启动Driver；

S200、选取同种群物种、同一条染色体的FASTA文件，作为所有输入的基因序列的参考序列；从作为参考序列的FASTA文件应用哈希算法生成第一类参考文件，然后使用广播变量方法上传到每个计算节点的executor进程作为缓存，用于加速匹配进程；

S300、为读取到的碱基序列与标识符序列创建Spark任务，此Spark任务的输入信息包括FASTQ文件包含的碱基序列、标识符序列和参考序列，输出信息为包括<参考的序列号，匹配位置，匹配长度>三元组列表，然后转换成可以被二次识别的字符串格式输出至HDFS；并且从图2可以看出步骤S300在具体执行时，包括下列步骤：

S301、分片读取:

使用MapPartition算子对碱基序列与标识符序列进行分片，对每一分片逐行读取基因的ACGT字符信息、N字符信息、特殊字符信息和换行符信息，记录每一行的长度；

对于N字符信息，用<相对行首的偏移位置,出现长度>的二元组形式表示；特殊字符以<相对行首的偏移位置，特殊字符>的二元组表示；行尾的换行符使用整数表示其相对行首的偏移位置；

ACGT字符信息按照k-mer方式进行哈希编码，在内存中建立哈希索引数组和哈希冲突数组；根据步骤S200中的哈希算法，计算该ACGT信息在哈希索引数组与哈希冲突数组中对应的哈希值；再遍历参考序列，找到该哈希值对应的k-mer元素并记录其在参考序列中的位置号，以下简称参考位置号，然后将参考序列与待压缩序列的位置指针同时后移，连续比较直到ACGT信息不再匹配为止；

设置片段的匹配度门限，对低于该值或者违背规则的片段，以<哈希值，相对偏移量，字符串>三元组形式另外cache序列化到内存中，其中的相对偏移量表示该片段相对前一匹配片段的偏移量，字符串即该片段的文字表示；

对高于该值而且符合规则的片段，则生成一个格式为<参考位置号，匹配位置，匹配长度>的三元组，其中的参考位置号表示该片段在参考序列中对应的位置号，匹配位置表示该片段在待压缩FASTQ文件中的位置号，匹配长度即该片段的长度；一个待压缩 FASTQ文件的所有碱基序列信息被编程语言表示为一个Java ArrayList RDD返回，使用cache算子的MEMORY_SER方式缓存这些中间结果，以降低内存利用率；

S302、三元组哈希编码

使用MapPartitionWithindex算法获取所有分片的序号，再用filter算子过滤出需要参与二次匹配的第一类参考文件并记录其序号；然后，对S301筛选出的一次匹配三元组结果进行哈希编码；

具体规则需脱离权利要求，从属的实施方式)其编码规则具体为：在实体变量中的未匹配到的字符串，对其字符每个字符的ASCII码乘以92803，加上三元组实体中的位置值乘以69091，再加上匹配相对长度值乘以51787，以保证哈希值的足够分散；

S303、三个数组发布

将S302得到的三元组哈希编码以及第一类参考文件序号信息通过first算子拉取到driver本地，对其封装成广播变量，然后发布到每个进行二次匹配的计算节点；

S304、二次匹配

在计算节点侧，根据得到的三元组哈希编码以及第一类参考文件序号信息，使用MapPartition算子对每一个三元组运用第一哈希算法，依次与第一类参考文件进行匹配，直到选取出最长匹配的元组，记录其出现的位置、长度以及对应的参考序列文件号；将所有匹配后的信息序列化为字符串，并输出至HDFS；

S400、为读取到的质量分数序列创建Spark任务，此任务的输入信息为以String表示的质量分数序列，输出信息为字节块，然后输出到HDFS，从图3可以看出步骤S400 在具体执行时，包括下列步骤：

S401、聚类门限计算

使用MapPartition算子，从质量分数序列文件中读取出十万行数据，对其中的每一行记录，使用k-mer读取方式对每一个k-mer计算其质量估算总值，其中：所述质量估算总值为将行内每个字符与7按位或，然后进行移位计算；如果出现两行的质量估算值相同，则判定为它们的k-mer质量相同，统计所有符合此条件的k-mer的出现频次，并按出现频次从大到小排序，选取出现频次前70％的k-mer，找到频次最高的k-mer的质量估算总值乘以一个划分类别的系数，定义为第一聚类门限，以作为聚类划分依据；

S402、频次极值和字符寻找

对质量分数序列，计算字符ASCII码分布频次的极值，将极值对应的字符记为第一字符，以在映射阶段作编码对照；

S403、字节块分类

以行为单位，从质量分数序列的每行中按k-mer计算每个k-mer的质量估算值，将每行计算的质量估算总值与第一聚类门限作比较，低于则进入缓冲类别1，高于则进入缓冲类别2；这两种缓冲类别将以字节块的形式存储；使用MapPartitionWithIndex获取文件号，并且使用reduceByKey算子，利用shuffle将具有相同缓冲类别与相同文件号的块数据划分到一个分区，以提升并行度；

S404、字符块映射

使用MapPartition进行字节的映射，将连续的字节或值较大的字符块映射成较小的字节块；对于连续的或者与第一字符相近范围内的字节块，一并映射成更小的字节块；处理完成后，上述字节块将输出到HDFS。

S500、driver本地程序从HDFS拉取压缩后的碱基序列与标识符序列、质量分数序列数据到本地，使用bsc压缩组件作最后的压缩运算，从而输出最终压缩结果。

更佳的批量分布式压缩方法的实施方式为，步骤S301输出的序列化字符串支持使用 Hadoop指令或者访问HDFS可视化网页端下载访问。

更佳的批量分布式压缩方法的实施方式为，Spark使用Yarn任务参数与配置调优来优化任务中每个计算槽的性能。使用Yarnstandalone模式，能够降低有着高吞吐任务的Driver对于工作节点的内存占用。

实施例

更为具体地讲，FASTQ基因序列并行化压缩的方案，实现批量化的基因序列碱基序列压缩的并行化以及质量分数全过程的并行化，以提高基因文件压缩速率。该方案分为碱基序列与质量分数两个部分，分别对这两个部分构建Spark算子，在碱基序列压缩阶段，使用与参考序列多次匹配的方式消去重复，在质量分数阶段使用分数值的聚类并按位映射。将多条待压缩序列并行处理，实现了对基因文件分片并行化读取与压缩操作。最终以达到实现压缩与提高基因压缩速率的目的。以下部分内容仅仅用于解释说明本发明的具体实施方式，请结合上述具体发明内容进行理解：

1.Spark输入格式设计

在云存储环境下，基因文件的分布式存储往往面临着数据块的多重备份，这也会带来计算任务向数据存储迁移时的便捷。而面对大型文件，存储系统一般会将大文件切成许多分片(Hadoop默认为128M)来进行存储，并有存储管理节点对这些分片的索引及位置进行划分并管理。经由性能优化分析一节的介绍，需要对文件为最小的输入单位进行并行计算，并且大数据计算平台一般会对接来自于分布式存储系统(HDFS)的输入，那么划分输入数据的分块结构便显得尤为重要。

文件在HDFS中较大的序列文件会被切分成若干数据块保存，如果不进行任务前的分片容量调整，一个待压缩文件会被分割给不同的Spark任务分别处理，由于处理过程相互独立，最终导致压缩结果出错。故在划分处理文件块之前需要统计记录划分的位置。输入分片容量C_split调整范围应为：C_file≤C_split＜C_block+C_file，其中C_block 代表HDFS中一个数据块的大小，C_file代表待压缩文件大小。此处为了保证压缩输入文件的文件局部一致性，将C_split设置为输入文件的大小C_file，可知，输入分片的数目也应等于输入文件的数目。在读取文件的同时，本发明将使用k-mer算法对FASTQ 文件序列的碱基序列与质量分数进行读取，k-mer是指将一行读取到的字符序列(蛋白质、碱基或者质量因子)分成包含k个碱基的字符串，一般长短为m的序列可以分成m-k+1 个k-mers。使用k-mer读取的好处即为，利用到了基因序列的相似性，使其在匹配寻址或者划分阈值的过程中更加准确灵活。

为了区分碱基字符与质量分数，并提升内存使用的效率，我们将规划两支Spark任务来分别完成质量分数与碱基和标识符序列。

2.质量分数的压缩

由于FASTQ质量分数分布的毫无规律，所以除非是有损压缩，使用无损压缩处理时无法使用一般的文本建模算法对其进行统一建模与预测。一般采取的无损压缩模式均为对其根据质量分数的不同进行聚类，并根据不同的聚类区域进行映射编码，以起到缩减重复部分的作用。但为了提升其压缩鲁棒性与压缩速度，将其聚类集合减少为2个。

1)质量分数聚类阈值获取

若要对质量分数序列进行字符编码层面的聚类，首先要对这些字符编码进行统一划线。由于质量分数分布的零散性，对于不同输入文件的质量分数部分中的质量得分划线将会有所差异，为了减少后续压缩器的等待时间，将采用下面的划线方法。

由于k-mer包含质量得分的上下文信息并也可表示读取位置的变化，因此质量分数的权重得分将会使用k-mers来计算，以确保压缩率和鲁棒性。考虑到简单性和有效性，此处的质量得分线由k-mers表示。同时，当质量分数线共享高频k-mers时，它们往往趋于相似，反之也是如此。因此，应为高频k-mers分配较高的权重，而为低频k-mers 分配较低的权重。为了简化权重计算过程，我们从输入数据中采样出一个子集进行权重的计算，当然是以假设该子集和整个数据集遵循相同的分布为前提的。我们以读取数据集的前M行(默认为十万行)为选择采样子集的标准。之后按行来采集每个k-mer出现的信息和所有k-mer的出现总数。一个k-mer的权重被指定为它的出现次数与所有k-mer 的总出现次数之比。每个质量得分行的权重等于其所有k-mers权重之和与其中k-mers 的数量之比。那么一行k-mer的权重得分越高，即表示该行将更加频繁地出现相同的质量分数k-mer。获得M个采样质量得分线中的最大权重，并将其用于归一化所有其他质量得分线的权重。将样本中出现频次最高的质量分数k-mers的最高权重得分记为S，此处设置一个常数α(0<α<1)，α与S的乘积即为划分聚类的阈值。在后边的步骤中，接下来从文件中读取的质量得分将根据其不同的线宽分为两部分，这种粗略的聚类方式足以实现压缩速度和资源占用之间的更好平衡。

在计算采样部分我们还将计算一个字符ASCII码分数的分布频次的极值，将其记为 C，它通过样本体现了这个质量分数文件的字符分布情况。它也在采样计算这一阶段完成，但与聚类逻辑无关。这个测量的标准值将被用于后边的映射步骤，通过C作为基准来测定出现连续且频次高的质量分数。

2)使用字节缓冲池构建块索引

在上一步在质量序列划分聚类阈值测定后，根据质量序列分布在ASCII码中33～104 位的特性，通过上一步的聚类方法把读取的字符串进行聚类并转化成字节流暂存。构建一个复用的字节缓冲池存储两种类别的聚类字符数组，以降低内存占用率。这里我们使用Netty框架提供的Java NIO缓存池ByteBuf类作为暂存媒介，并对这些数据进行下一阶段的处理。下面将介绍划分缓冲块的方式。

在获得了数据流的聚类划分线αS之后，将对文件中的数据按行来进行块类型的区分与字节缓冲区的转化。对于缓冲区1，如果一条质量得分行的内容属于该缓冲区，那么该行应该是采样行，将其复制到缓冲区1中。同时，如果读取出新的一行质量分数，测定权重小于等于聚类划分线，则属于缓冲区1，将其复制到其中，并标记索引以跟踪其位置信息。对于缓冲区2，如果一条质量得分行的内容测定权重(权重大于聚类划分线)后属于该缓冲区，则将其复制到缓冲区2中。在填充缓冲区2时，在缓冲区1中添加前一行与后一行的索引，以便于在解压中合并缓冲区2中的数据。

3)质量分数映射

对于聚类出的两种数据流，上一步将其转化为字节流暂存。之后将其通过某种方法进行数字映射，以实现多个相匹配字段，并最后作为字节类型存储以压缩存储空间。具体地，通过一对一映射规则将每个质量得分值从33到104的质量得分k聚体无损地转换成从1到255的单个数字。对质量分数编码映射方式的特性如下：

(a)在聚类阈值计算步骤曾介绍过，在采样集合中获取出现次数最高的质量得分值并将其表示为C。一旦获得C，则映射规则是固定的。我们选择出现次数最多的质量得分作为转换间隔的边界的原因在于，它与其他质量得分值组合形成k-mer的可能性最高。k-mer满足转化条件的条件越多，可以减少的含量就越大；

(b)质量得分的分组和映射。对于每个质量得分，k-mer分组和映射的优先级从左到右

增加。

其中详细映射方法如下步骤所示：

①开始，定义一个比较阈值变量l，

②打开一个ByteBuf缓冲块，从中取出一个字节的数据d；

③如果这个字节d＝C或者满足d>＝l并且d<＝l+7，则再取出一个字节d2(第一个字节将标记为d1)，进入步骤③。如果步骤③的条件不满足，则会映射这个单独的字节d。将会产生两种映射结果：如果d＝10，那么将映射为0；如果d满足l+7<d<l，则将其质量分数平移缩小，映射为d-32。

④如果连续两个字节(d1,d2)在满足d1＝C并且d2＝C，或者d1，d2均满足l+4<＝d<＝ l+7的条件下，继续从缓冲区中取出新的字节d3，进入步骤④。如果步骤④的条件不满足，那么将单独映射(d1,d2)这两个字节。首先如果d1＝d2＝C，则d1 被映射为199+1，对d2使用②中的映射方法；如不满足以上条件，如果d2满足l+4<＝d2<＝l+7则将它们映射成73+(d1-l)+8×(d2-l)；如果不满足则将d1输出为d2-32，对d2使用②中的映射方法输出。

⑤如果连续三个字节(d1,d2,d3)满足d1＝d2＝d3＝C则取出新的字节d4进入步骤⑤。如果步骤⑤的条件不满足，则映射这3个字节。如果d1＝d2＝d3＝C，则将d1，d2 映射为2+199，对d3使用①中的映射方法；如不满足以上条件，如果d3满足 l+4<＝d<＝l+7，则将这三个字节映射为137+(d1-l3)+4×(d2-l3)+16 *(d3-l3)；如果不满足则将d1，d2映射为73+(d1-l)+8×(d2-l)，对d3使用①中的映射方法输出。

⑥对于连续n(n>＝3)个的质量分数d＝C，映射这n个字节为n+199，仅限n<＝55。如果n>＝55则从此截断，对后继部分依旧按n’＝n’+199映射为一个新的字节。

这种数字映射方式可以保证每一个字节将会被尽可能地变成连续或者有规律的字节，映射后的字符将会分布在某个特定的区间内。便于通用压缩器对其进行统一建模并去重。

4)Spark并行化算法构建

这一部分介绍质量序列的并行化计算任务流。Spark的输入格式依赖于HDFS与Hadoop的Inputformat，我们按照默认格式String进行输入，输入的数据是按照文件为单元划分分区。

在Spark程序流程中，使用各种算子组织对数据流的计算。首先对第一个RDD对象进行MapPartitionWithIndex算子的处理，该算子的功能是在使用MapPartition处理每一个分区的同时，并且在每一个分区的迭代器中加入一个分区号变成二元组形式的数据。在这一个算子中，我们对每一个分区的文件依次进行读取、采样与划分每一个质量分数文件中的质量分数聚类，这一步骤将在每个RDD的分区中生成许多基于字节与划分类型的缓冲块。之后，算子将会给予每一个分区一个唯一的分区号，将在后面阶段用于区分不同分区的缓冲块。我们将在这个分区号的基础上，增加每个块聚类的结果，将在下一阶段划分至不同的分区。

下一步骤将使用reduceByKey算子，根据上一步骤划分的二元组的键数据进行排列重组，相同的Key将被划分至一个分区。由于聚类的数量一般为两个，所以此处将会把分区数依据划线阈值αS提升至最多两倍，那么并行度也将随分区数目的提升而提升。此处定义分区器来实现在shuffle write中，相同key的聚类。

之后再使用MapPartition算子进行下一操作，此时分区数目已经改变，同一个文件的同一种聚类缓冲块被划分至一个分区内。Partition算子与Map算子最大的区别是，MapPartition算子一次性将全部分区的数据代理发送给Executor进行计算，Executor 将在内存中缓存全部的数据，而Map算子是逐条发送。使用Partition算子减少了数据发送给Executor的次数，可以确保执行的效率。在该算子中，执行字节的映射算法，将放缩之后的字节刷新ByteBuf缓冲区中的内容，最后将从ByteBuf中依次提取出缓冲区里的每一个字节，存放到一个Byte类型数组中。分区依旧定义为String类型，包含了文件分区与聚类分区两层的分区索引信息。以此形式通过Spark的saveAsObjectFile 存储至HDFS中。这种存储对象的方式将在按照格式输出对象内容的同时，还会在文件的头部留存一行Hadoop序列化存储对象类型信息。Byte类型即为ByteWritable类型， String类型将被转化为Text类型。在解压缩的时候，程序将调用Hadoop的解析类来解析这个特殊类型。

至此质量分数序列已经完全压缩完毕，将会保存至一个单独的HDFS目录中，每一个分区的文件将会按照默认方式单独保存为一个文件。

3.碱基序列与标识符序列的压缩

本方法使用k-mer算法对基因序列进行读取，当k取合适值时，DNA序列中k-mer 频数分布包含了基因组的全部信息，从而构成序列的一个等价表示。使用k-mer匹配会更加精准快速，Hash表的创建不再是将庞大的序列全部转化成整数数组，而是基于k-mer 逐段创建。这样会节省复杂的Hash计算，简化Hash表创建过程，节省Hash表创建的内存消耗。

1)序列读取

由于匹配算法需要使用两个碱基序列文件进行交互，所以需要先提取它们的关键信息。由于这些基因文件被存储于HDFS中，这是一个分布式的抽象存储系统，在处理输入文件时需要考虑输入分片大小的设置与输入文件的缓存配置。

首先从文件中读取参考序列，经分析可知，文件中ACGT相似度较高，其余字符如N、Y、R等出现频率低，分布十分散乱，所以只对其记录位置与编号，不作消重处理。所以提取参考序列过程只需提取碱基字符序列信息。将参考序列的小写字符序列转化成大写字符，通过分配一个连续空间的数组G_reference储存所有的序列字符，存储的字符只有大写的A、C、G、T。

下一步将按Reads从文件中，按照分段方式提取待压缩碱基序列，使用上文所述的分步读取方式读取碱基序列。由于FASTQ文件体积庞大，并且许多Reads会有交错重复，为了节省内存，此处设置一个文件块大小的阈值用以分段读取文件，即对每个数据块的碱基序列按照下面方式读取：

(a)将该文件块的所有碱基字符存储在数组G_target中；

(b)读取N字符信息、换行符信息与其他字符信息。N字符信息同样使用<p_n,l_n>的向量形式存储在数组G_n；由于不同FASTA文件每行长度不定，而且同一文件中也可能有行不满的情况，使用换行符的位置来记录每行的长度，将这些长度保存到一个数组 G_line中；对于剩下的未知碱基，会使用一些排列无规律的特定字符对其表示，这些字符信息需要使用<p_s,c_info>的形式存储在数组G_special中，其中p_s为特殊字符的位置，依旧采取相对位置来表示，c_info表示其字符。该步骤将对待压缩序列扫描一遍，故时间复杂度为O(n)。

这一过程结束后，下面会直接使用此步骤提取的块信息进行匹配去重。

2)初次匹配压缩

首先需要构建初次匹配用于索引的哈希表。使用K-mer读取方法再次扫描参考基因序列G_reference，各位碱基字符按照A＝0，C＝1，G＝2，T＝3的编码方式进行整形数字编码后，依据哈希值计算公式创建散列表，使用散列表的哈希值来搜索匹配拥有相同的哈希值的G_target中的字符片段。根据生物学的匹配意义，使用K-mer按位读取方法可以使匹配过程无字符遗漏，在确定K-mer读取步长k后，每次读取一个K-mer都会包含一串长度为k的数字序列S_j，其中j＝0,1,2,…,k-1。根据S_j计算哈希值v，v计算方法如公式1所示。

碱基序列信息经第二步都已存入各自的数组中，接下来需要对其进行匹配去重操作。这一阶段主要使用待压缩序列字符数组G_target与参考序列数组G_reference匹配，并直接保存其他字符信息。

对参考序列构建哈希表后，对G_target中的数据进行哈希编码得到值v，由G_hash[v]中的值回溯G_loc数组元素的位置p，这表示从G_target[p]到G_target[p+k]与G_reference[G_hash[v]]至G_reference[G_hash[v]+k]，两段K-mer完全相同。在这次匹配中将会对待压缩序列进行重新寻址，以契合参考序列的不同位置。在此之后参考序列与待压缩序列继续向后逐位推进，继续寻找相同的碱基字符直到无匹配字符，记录这段匹配长度l_match。由于在G_hash中同值的元素可能有很多，如果G_loc相应位置不为-1，则会继续向前回溯并匹配，直到匹配至-1为止。为了保证匹配的长度，设置阈值m作为匹配长度l_match的最小值，若l_match小于m，则该处的G_loc位置不做考虑。最后取得一个匹配长度的最大值l_(match_max)，代表匹配到了最长序列。使用 l_match，表示l_relative超出m部分的长度，公式表示为l_relative＝l_(match_max)-m。在这其中m会影响匹配片段的长度，将会影响最终结果的压缩率。未匹配到的字符存入 G_mismatch中，将其由字符型转化为整形存储，以便下一阶段的数据流动与转存。存储片段后，接下来将从最大片段之后(位置为〖k+l〗_(match_max))继续计算哈希值并向后匹配，直至扫描完毕全部碱基序列。最终使用<p_relative，l_relative，G_mismatch> 的三元组形式保存一个输入文件k-mer在参考序列中匹配片段的信息。

在读取阶段还产生N字符信息数组G_n、换行符信息数组G_line、特殊字符信息数组G_special。这些数组中存储的信息不会在参考序列中提取。由于换行符信息数组 G_line储存的是基因文件一行的长度，而且每行长度几乎相同，所以数组G_line中会有很多冗余信息。使用游程编码对其进行处理，以<b_position,b_num>二元组形式覆盖G_line，其中b_position为相同行宽的起始位置，b_num为从b_position开始相同行宽的行数。这一过程结束后，下一步会把这些数字信息放入列表中等待输出。

3)二次匹配压缩

二次匹配压缩中，需要使用部分第一次匹配结果作为参考序列的输入集。这一步骤与初次匹配的操作类似，此处依然会依照一个足以解决这些输入集产生哈希碰撞问题的编码方式，并将二次参考序列输入集的数据依次编码。

此外，在二次匹配中将采用多参考序列顺序匹配的方式，需要在所有参考序列中搜索最长的匹配。为了尽可能确保不同的散列值代表不同的匹配实体，我们在计算中涉及到位置、长度和每个未匹配到的碱基字符，并使用大素数作为取余编码器。其编码方式为：在实体变量中的未匹配到的字符串，对其字符每个字符的ASCII码乘以92803，加上三元组实体中的位置值乘以69091，再加上匹配相对长度值乘以51787。这样的编码可以保证哈希值的足够分散。

在采样一次匹配结果作为参考序列组的时候，令输入的数目为N，当N越大时，二次匹配的压缩效果将会越好，但中间结果产生的内存占用也将越大。二次匹配建立哈希表的数据结构与初次匹配类似，也将建立两个数据结构，一个作为桶来存储哈希值，另一个通过链式方式处理冲突并将索引存放至数组中。

在建立哈希表之后，初次匹配的匹配结果将继续被进行二次匹配。二次匹配依然需要计算输入元素的哈希值，但与第一级匹配不同的是，用于匹配的元素不仅是四个基本字符{A、C、G、T}编码形成的整形数字，而是匹配的初次压缩结果集实体<p_relative， l_relative，G_mismatch>。

在搜索过程中，将相邻两个待压缩匹配实体以上文的方式计算哈希值，并逐个遍历每个参考序列的哈希表，找出是否存在具有相同哈希值的实体。但是，与第一级匹配不同，相同的哈希值并不意味着相同的实体，因为在第二级匹配中，匹配的实体和散列值之间的映射不是一对一的映射。因此，当发现相同的散列值时，有必要验证匹配的实体是否相同。如果是，则继续搜索连续的匹配实体，直到找到不同的匹配实体，然后记录匹配段的序列id、位置和长度。

对于前N条作为参考的实体而言，它们也需要被二次压缩，将被一个单独的处理方式处理：第n条作为参考的实例对象只会与1，2……，n-1条参考序列进行匹配(n<＝N)，以便于后期解压缩中参考序列的依次还原。而后边的所有序列将只会与这N条参考序列进行匹配。在遍历所有的参考序列后，取最长的匹配段作为最终的匹配段进行替换，并存储为三元组(序列号、位置、长度)。对于匹配的实体，如果在遍历所有参考序列后没有找到相同的连续匹配实体，则直接将该压缩实体标记为未匹配，并执行下一步对于三元组向量存储。

4)标识符序列的压缩

标识符行是通过测序仪器十分系统地生成出的，以Illumina测序技术为例，会在一个Read的第一行生成富有规律的有关该read测序信息的标识符序列。针对这些标识符序列采取以下几步编码策略：

(a)游程编码：使用<position,length>二元组压缩一组连续记录中具有相同Token 值的集合，以表示相同Token出现的位置与长度。对于整个数据集中具有恒定值的集合，我们也使用相同的压缩技术。

(b)垂直编码：使用垂直编码压缩连续的整数Token值单调递增或递减的差值，以缩短较长整数数列的数字位数。

(c)不压缩具有不属于上述任何类别的Token的集合。

5)Spark并行化流程

首先需要在任务提交端准备一个FASTA格式参考序列，它将被程序的主函数读取并封装成参考类型对象。为了使得所有工作节点都能共享这个分布式环境下的全局变量，我们使用广播变量的方式给每个节点Executor进程分发该变量并缓存一份副本。如果不使用广播变量，由于匹配算法的输入是参考与待压缩两条序列，将导致参考序列的 HDFS路径被不同计算进程重复读取，而且每次读取都有可能产生磁盘读取与网络传输，在读取大批量基因文件的情况下会产生大量时间损耗；而且这种重复读取的文件在执行过程中会缓存在每个Task中，由于Executor与Task的比例是一对多，那么将产生大量不必要的内存消耗。针对这一问题，在构建应用程序时采取Spark的广播技术，其目的是在每个Spark任务运行之前，将driver客户端程序中构建的广播变量或对象通过网络发送到每个Executor进程中，缓存在Spark的缓存内存中作为该Executor内部共享文件。缓存的对象个数为Executor启动进程的数目。这样网络传输的次数仅为节点总数，在实际处理中，即使是最大的1号染色体，在分布式多副本存储环境下，缓存速度也没有令人失望。

由于并行化算法的复杂性，对碱基序列的并行压缩算法将分两阶段介绍。与质量序列一致，以字符串格式从HDFS中读取文件输入分片，加载至每一个分区，并由每一个计算任务负责接收输入分片按行转化成的字符串迭代器。此时将使用MapPartition算子，并行地对每一个分区的数据完成由字符串至可操作的Read类型转化，此处将在读取迭代器内容的同时，转化标识符序列与N、特殊字符信息。在接收到主程序(运行在 driver端)封装的，同样转换成参考文件类型的广播变量对象之后，将在每一个计算节点的内存空间复制并执行初次匹配操作。完成初次匹配步骤之后，将形成一个二元对象组的数据结构，其中包括其他字符(标识符、N字符与特殊字符)组成的字符串列表，与初次匹配的结果实体列表。

接下来需要使用WithIndex算子对所有的分区加入分区号，以分区号作为索引，用于在下一阶段通过筛选算子选择二次匹配采样的参考序列数目。一般采样序列的筛选方式是选取输入列表中的前N个序列(N为上文中介绍的选取参数)，由开头所述的FIFO 任务调度过程，可以使用RDD的分区号作为判断序列输入的顺序的标志。

接下来进入匹配流程的第二阶段，这一阶段的RDD将会产生计算分支，也将划分出两个阶段。为了使得上一阶段的结果RDD不会丢失，也不会随着任务阶段的消失而重新计算。所以我们在这一阶段通过缓存方式暂存该RDD的数据，使用persist算子，将缓存等级调整为MEMORY-SER，使用kyro序列化方式将内存中数据压缩后缓存在Executor 中，既保证了存取速度，由不会占用太多的内存空间。

第二阶段首先对输入RDD使用Filter算子通过分区号的判断(x是否<＝N)过滤出二次采样的序列数目。在处理后将会过滤出N个分区的数据，下一步需要将这些分区合并，以打包成一个文件构建二阶段的广播变量。使用RePartition算子，利用shuffle 方式修改分区数量为1。这一算子之后将对这整个分区的数据进行计算，构建二次匹配的参考序列哈希表与匹配向量组。在按照上文介绍的方式构建结束后，构建一个包括哈希桶、哈希冲突位置与二次参考的实体列表三种类型的三元组。在Driver端通过first 方法读取该单一分区的值，并在提交端构建广播变量，这一阶段将会产生一定量的网络传输。

再次获取一阶段的RDD结果数据，通过MapPartition方式进行二次压缩的匹配过程。在该过程中缓存在内存中的数据将被取出，读取并交互广播变量，匹配后将向量间加入空格，全部转化成String格式，通过TextFile方法输出至HDFS。

下一阶段将会把HDFS中质量分数与碱基序列两部分的压缩结果通过Driver端的FileSystem客户端对象读取至本地，并使用bsc压缩工具进行最后的压缩。

4.通用压缩器压缩

BSC是一种基于Burrows-Wheeler变换(BWT)的通用压缩方法，可以实现文件高压缩比的压缩。BSC利用多线程并行计算，计算效率高。在该算法中，我们在算法的最后部分使用BSC来进一步提升压缩效率。将输出的压缩结果文件集使用apache的tar方式打包后，使用BSC压缩器作最后的压缩手段。该压缩器将会使用混合建模算法对质量序列的打包文件与碱基序列的匹配结果文件实现进一步建模去重，对于有序的数据块，压缩效率将会进一步提升。

本发明的实施例公布的是较佳的实施例，但并不局限于此，本领域的普通技术人员，极易根据上述实施例，领会本发明的精神，并做出不同的引申和变化，但只要不脱离本发明的精神，都在本发明的保护范围内。

Claims (3)

1.基于FASTQ基因大数据的批量分布式压缩方法，其特征在于：包括如下步骤：

S100、部署Spark集群、Hadoop集群和Java环境：

将所有参与压缩的FASTQ文件上传至HDFS上；

在压缩环境中预先选取人类基因组序列hg38的参考序列信息和质量分数序列信息，并且放置到程序本地；

将Spark集群调整成Yarn standalone模式，在单独节点启动Driver；

S200、选取同种群物种、同一条染色体的FASTA文件，作为所有输入的基因序列的参考序列；从作为参考序列的FASTA文件应用哈希算法生成第一类参考文件，然后使用广播变量方法上传到每个计算节点的executor进程作为缓存，用于加速匹配进程；

S300、为读取到的碱基序列与标识符序列创建Spark任务，此Spark任务的输入信息包括FASTQ文件包含的碱基序列、标识符序列和参考序列，输出信息为包括<参考的序列号，匹配位置，匹配长度>三元组列表，然后转换成可以被二次识别的字符串格式输出至HDFS；并且：

S301、分片读取:

使用MapPartition算子对碱基序列与标识符序列进行分片，对每一分片逐行读取基因的ACGT字符信息、N字符信息、特殊字符信息和换行符信息，记录每一行的长度；

对于N字符信息，用<相对行首的偏移位置,出现长度>的二元组形式表示；特殊字符以<相对行首的偏移位置，特殊字符>的二元组表示；行尾的换行符使用整数表示其相对行首的偏移位置；

ACGT字符信息按照k-mer方式进行哈希编码，在内存中建立哈希索引数组和哈希冲突数组；根据步骤S200中的哈希算法，计算该ACGT信息在哈希索引数组与哈希冲突数组中对应的哈希值；再遍历参考序列，找到该哈希值对应的k-mer元素并记录其在参考序列中的位置号，以下简称参考位置号，然后将参考序列与待压缩序列的位置指针同时后移，连续比较直到ACGT信息不再匹配为止；

设置片段的匹配度门限，对低于该值或者违背规则的片段，以<哈希值，相对偏移量，字符串>三元组形式另外cache序列化到内存中，其中的相对偏移量表示该片段相对前一匹配片段的偏移量，字符串即该片段的文字表示；

对高于该值而且符合规则的片段，则生成一个格式为<参考位置号，匹配位置，匹配长度>的三元组，其中的参考位置号表示该片段在参考序列中对应的位置号，匹配位置表示该片段在待压缩FASTQ文件中的位置号，匹配长度即该片段的长度；一个待压缩FASTQ文件的所有碱基序列信息被编程语言表示为一个Java ArrayList RDD返回，使用cache算子的MEMORY_SER方式缓存这些中间结果；

S302、三元组哈希编码

使用MapPartitionWithindex算法获取所有分片的序号，再用filter算子过滤出需要参与二次匹配的第一类参考文件并记录其序号；然后，对S301筛选出的一次匹配三元组结果进行哈希编码；

S303、三个数组发布

将S302得到的三元组哈希编码以及第一类参考文件序号信息通过first算子拉取到driver本地，对其封装成广播变量，然后发布到每个进行二次匹配的计算节点；

S304、二次匹配

在计算节点侧，根据得到的三元组哈希编码以及第一类参考文件序号信息，使用MapPartition算子对每一个三元组运用第一哈希算法，依次与第一类参考文件进行匹配，直到选取出最长匹配的元组，记录其出现的位置、长度以及对应的参考序列文件号；将所有匹配后的信息序列化为字符串，并输出至HDFS；

S400、为读取到的质量分数序列创建Spark任务，此任务的输入信息为以String表示的质量分数序列，输出信息为字节块，然后输出到HDFS，具体步骤为：

S401、聚类门限计算

使用MapPartition算子，从质量分数序列文件中读取出十万行数据，对其中的每一行记录，使用k-mer读取方式对每一个k-mer计算其质量估算总值，其中：所述质量估算总值为将行内每个字符与7按位或，然后进行移位计算；如果出现两行的质量估算值相同，则判定为它们的k-mer质量相同，统计所有符合此条件的k-mer的出现频次，并按出现频次从大到小排序，选取出现频次前70％的k-mer，找到频次最高的k-mer的质量估算总值乘以一个划分类别的系数，定义为第一聚类门限，以作为聚类划分依据；

S402、频次极值和字符寻找

对质量分数序列，计算字符ASCII码分布频次的极值，将极值对应的字符记为第一字符，以在映射阶段作编码对照；

S403、字节块分类

以行为单位，从质量分数序列的每行中按k-mer计算每个k-mer的质量估算值，将每行计算的质量估算总值与第一聚类门限作比较，低于则进入缓冲类别1，高于则进入缓冲类别2；这两种缓冲类别将以字节块的形式存储；使用MapPartitionWithIndex获取文件号，并且使用reduceByKey算子，利用shuffle将具有相同缓冲类别与相同文件号的块数据划分到一个分区，以提升并行度；

S404、字符块映射

使用MapPartition进行字节的映射，将连续的字节或值较大的字符块映射成较小的字节块；对于连续的或者与第一字符相近范围内的字节块，一并映射成更小的字节块；处理完成后，上述字节块将输出到HDFS。

S500、driver本地程序从HDFS拉取压缩后的碱基序列与标识符序列、质量分数序列数据到本地，使用bsc压缩组件作最后的压缩运算，从而输出最终压缩结果。

2.根据权利要求1所述的批量分布式压缩方法，其特征在于，步骤S301输出的序列化字符串支持使用Hadoop指令或者访问HDFS可视化网页端下载访问。

3.根据权利要求1或2所述的批量分布式压缩方法，其特征在于，Spark使用Yarn任务参数与配置调优来优化任务中每个计算槽的性能。

Images