CN110021349B - 基因数据的编码方法 - Google Patents

Abstract

本发明提供一种基因数据的编码方法。该方法包括：对基因数据中AGCT字符及其原质量分数进行编码；基于预定的编码转换表对基因数据中除AGCT字符之外的特殊字符及其原质量分数进行编码，以获得编码之后的碱基序列和质量分数序列。根据本发明的编码方法，能够对基因数据中的碱基序列和质量分数序列进行压缩，从而节约存储空间并加快基因数据的处理速度。

Description

基因数据的编码方法

技术领域

本发明涉及基因测序数据分析技术领域，尤其涉及一种基因数据的编码方法。

背景技术

随着第二代基因测序技术的发展，测序通量直线上升，以人为例，完成一次基因组测序可产生数以百万计数目的测序序列，对大量的基因数据分析逐渐成为整个流程中的瓶颈。目前，人们通常使用分布式技术来提高对基因数据分析的处理速度，例如，将当前主流的大数据技术Spark云平台应用于基因数据分析的并行和加速中。

然而，采用Spark计算时，绝大部分数据都位于内存中，对于集群的整体内存和单节点的需求都较高。在发生内存不足的情况下，将数据写入磁盘时，会严重拖慢系统的整体性能。而在目前的需求下，基因数据的总量非常大，原始输入数据总量可以达到数百GB，产生的中间文件体积可以达到几TB。因此，在用Spark进行基因数据分析处理时，往往面临集群内存不足的情况。此外，海量的数据总量容易导致Spark在进行Shuffle操作(例如排序、聚合等)和Cache操作(例如写回)时引发庞大的网络和磁盘开销。

现有的基因数据压缩算法，通常需要预先统计数据集上的一些特征，以便于达到更高的压缩效率。这在Spark中又带来了额外的遍历、计算和Shuffle开销。

因此，需要进一步考虑针对基因数据的处理方法，以解决分析基因数据时负载过重的问题，从而使得更适合于利用Spark平台分析处理基因数据。

发明内容

因此，本发明的目的在于克服上述现有技术的缺陷，提供一种基因数据的编码方法，其包括通过对碱基序列和质量分数序列进行编码，以实现数据压缩。

根据本发明的第一方面，提供了一种基因数据的编码方法。该方法包括：

步骤1：对基因数据中AGCT字符及其原质量分数进行编码；

步骤2：基于预定的编码转换表对基因数据中除AGCT字符之外的特殊字符及其原质量分数进行编码，以获得编码之后的碱基序列和质量分数序列。

在本发明的基因数据的编码方法中，步骤2包括：

步骤21：将所述特殊字符的原质量分数转化为有效质量分数；

步骤22：基于所述预定的编码转换表确定所述特殊字符及其有效质量分数所对应的编码转换之后的编码，其中，所述预定的编码转换表存储的是所述特殊字符和其有效质量分数的组合在编码转换之后与AGCT字符和质量分数的组合的对应关系。

在本发明的基因数据的编码方法中，所述特殊字符的原质量分数的范围是[33-126]，在编码转换之后，所述特殊字符的质量分数的范围是[0-32，127]。

在本发明的基因数据的编码方法中，所述有效质量分数的范围为[0-10]。

在本发明的基因数据的编码方法中，在步骤1和步骤2中，以基因数据中的一条基因测序序列为存取单元进行编码。

在本发明的基因数据的编码方法中，所述特殊字符包括R、Y、K、M、S、W、B、D、H、V、N和占位符。

在本发明的基因数据的编码方法中，在步骤2之后还包括：

步骤3：将所述编码之后的质量分数序列存储为差值形式的序列；

步骤4：对所述差值形式的序列基于哈夫曼编码方式进行编码压缩。

根据本发明的第二方面，提供了一种基因数据的处理方法。该处理方法包括：

根据本发明的基因数据的编码方法对基因数据编码，以获得编码之后的碱基序列和质量分数序列；在spark环境中，使用预定的哈夫曼树对所述编码之后的质量分数序列进行压缩。

与现有技术相比，本发明的优点在于：本发明通过对碱基序列的特殊字符进行编码，能够压缩碱基序列；进一步地，通过对碱基序列中的质量分数序列进行哈夫曼编码等来压缩质量分数序列；在Spark中采用本发明的压缩方法对基因数据进行分析时，能够以较低的计算代价对单条测序序列实现压缩，解决了Spark处理基因数据时内存负载过重的问题，降低了Shuffle和cache操作对磁盘和网络的开销，同时不妨碍对基因数据的细粒度并行处理。

附图说明

以下附图仅对本发明作示意性的说明和解释，并不用于限定本发明的范围，其中：

图1示出了根据本发明一个实施例的基因数据的编码方法的流程图。

具体实施方式

为了使本发明的目的、技术方案、设计方法及优点更加清楚明了，以下结合附图通过具体实施例对本发明进一步详细说明。应当理解，此处所描述的具体实施例仅用以解释本发明，并不用于限定本发明。

图1示出了根据本发明一个实施例的基因数据的编码压缩方法的流程图。具体包括以下步骤：

1)S110，对碱基序列中的典型字符和其对应的质量分数进行编码

目前，基因测试技术产生的基因数据格式包括FASTQ、SAM和VCF等，以FASTQ格式为例，其每条记录的主要字段包括测序序列(read)的描述行、碱基序列(sequence)和质量分数(quality)序列。对于SAM格式，其每条记录的主要字段包括测序序列的名字、flag(正负链的标志位)、cigar(与参考基因的比对结果)、碱基序列、质量分数序列和一些属性标签等。在FASTQ、SAM等这些基因数据格式中，碱基序列和质量分数序列约占整个基因数据记录的80％至90％。因此，对碱基序列和质量分数序列的编码压缩能够有效的减小基因数据量。

在一个实施例中，对碱基序列中的典型字符A、G、C、T采用2-bit的编码方式，例如，A用00表示，G用01表示，C用10表示，T用11表示，而典型字符的质量分数序列的范围对应ASCII码[33-126]。

2)S120，对碱基序列中的特殊字符和其对应的质量分数进行编码

对于碱基序列字段来说，其典型地由AGCT四种字符组成，但同时也可能包括基因测序时无法确定的特殊字符，例如，N、R、Y、K等。简单的2-bit编码无法完整地表达这些特殊字符的信息。

为了仍可用2bit来实现对特殊字符和其质量分数的编码，在本发明中，增加质量分数的表示范围[0-32]和[127]，将特殊字符的质量分数转化为质量分数字符的ASCII码值减33，例如，如果某个特殊字符的质量分数的ASCII码值为35，则转化后的质量分数为2，在本文中，为了清楚起见，将经过这种转化之后的特殊字符对应的质量分数定义为有效质量分数。

将质量分数的表示范围扩展到0-127之后，可以用34*4(即136)种组合来表示特殊字符和其有效质量分数的组合。碱基序列中可能出现的特殊字符共有12种，例如，R、Y、K、M、S、W、B、D、H、V、N和占位符“-”，因此，可以表示的特殊字符的最大有效质量分数范围是0-10，即对于特殊字符和其有效质量分数有12*11＝132中组合。鉴于该类特殊字符本身的有效质量分数通常很低，而且该类碱基的出现本身就代表着可信度低，无法对测序结果进行支持，因此，对于这些特殊字符，当有效质量分数大于10时，可将其强制转换为10，这种转化不会对最终的基因测序结果产生明显影响。应理解的是，本领域的技术人员根据需要将有效质量分数的范围设定为其它值，例如[0-8]，[0-6]等。

为了进一步提高基因数据的处理速度，可将所有特殊字符和其有效质量分数组合的编码方式以预定的表格进行存储，以备后续使用。表1示出了根据本发明一个实施例的12种特殊字符的编码转换表(不考虑典型字符的正常编码部分)。其中，各单元格为特殊字符及其对应的有效质量分数的组合，行标为转换后的碱基字符，列标为转换后的质量分数。以表1中N1为例，其表示原碱基为N，该碱基对应的有效质量分数为1时，转换后的编码为A，质量分数为1。应注意的是，在表1中特殊字符中大于10的有效质量分数按照10进行处理；V11-V24仅为填补表中空白，并无特殊意义；表1仅用于示意，本领域的技术人员可基于本发明的思想进行多种变型，例如，调整特殊字符的顺序、缩减有效质量分数的范围等。

表1

除了碱基序列采用2-bit编码之外，还需使用一个或多个字节来存储编码前碱基序列的长度。表2示出了对包含特殊字符的碱基序列进行编码的示意。具体地，将特殊字符N转换为A，将N对应的原质量分数字符#转换为字符STX(其转换过程为：字符#对应的ASCII值为35，其对应的有效质量分数为2(即35-33＝2)，根据ASCII表，2对应的字符是STX)，将转换后的碱基序列GGTTACCTA采用2bit编码后得到二进制序列，其中，字节00001001表示原序列长度为9个字符。通过这种编码方式，碱基序列的压缩率可以达到接近4倍。

表2

在通过上述方式对碱基序列进行转换之后，质量分数序列的取值范围扩展到了[0-127]。在一个优选实施例中，将质量分数序列转换为存储质量分数之差的序列(delta序列)，将质量序列的取值范围扩展为[-127-127]，参见表3，这主要是基于以下特征：一条测序序列的质量分数可以分为多段，在每一段内，某个位置的质量分数与前一个位置相近，因此，不同样本相邻质量分数之差相比质量分数本身分布更为集中且更容易预测，并且绝大多数相邻的质量分数差值都位于[0-10]之间。

表3

3)S130，将编码后的质量分数序列进行压缩

通过将质量分数序列存储为差值的方式，虽然不同的数据集在整体质量分数上存在差别，但是相邻质量分数差值却基本服从同样的分布，可对差值的质量分数序列进行压缩，并在压缩序列的最后加入一个EOF的符号，表示压缩结束。

例如，使用哈夫曼(Huffman)编码方式进行压缩，哈夫曼编码的基本思想是以根据待压缩数据中各符号出现的频率，采用可变长度的编码方法，对出现频率越高的符号采用越短的编码长度，从而最终获得更短的总编码长度。例如，对于表3中的差值序列[67、0、0、1、‐64、68、0、0、0]，由于0出现的频率最高，可将其编码为01(二进制)，而68出现的概率相对较低，则可以将其编码为0010，在实际应用中，通常通过构建哈夫曼树利用程序代码来实现对质量分数的编码，为了节省在整个数据集上进行字符统计的开销，对于所有数据集都可以使用预先设定好的Huffman树进行压缩。构建哈夫曼树的过程属于现有技术，在此不再赘述。

应理解的是，除了哈夫曼(Huffman)编码方式之外，还可以采用其它的方式实现数据压缩的目的，例如Shannon-Famo编码、游程编码(RLE：run-length encoding)或RICE编码等。

本发明的对特殊字符进行编码的方式以及对质量分数序列进行压缩的方式可以应用于任何基因数据分析平台，以实现节约存储空间目的，此外，为了使本发明更适合于分布式数据平台，上述的编码方式可以每条测序序列的记录为一个存取单元进行编码或压缩，以不妨碍例如用于去除重复数据的Mark Duplicate等算法的细粒度并行，从而达到加速数据处理的目的。

例如，利用本发明的编码压缩方法在Spark云平台进行计算的过程如下：

a)在Spark环境中，对Huffman树在多个executor之间进行广播。这个Huffman树可以是静态的。

b)根据碱基序列，对质量分数序列进行转码。

当碱基是A、G、C、T时，不需要转码；当碱基是N之类的特殊字符时，则与当前的质量分数一起，根据表1的转化表将其进行转化，得到的碱基序列全部由A、G、C、T四种碱基构成，质量分数由0-127的数组构成。

c)对碱基序列进行2-bit编码。

例如，使用一个字节存储编码后的二进制序列，根据需要，可调整为两个字节或更多字节，然后每四个碱基为一个字节进行编码，得到编码后的碱基序列。

d)使用Huffman树对转码后的质量分数进行压缩，得到压缩后的质量分数序列。

对应地，解压步骤如下：

a)在Spark环境中，对Huffman树在多个进程执行单元(executor)之间进行广播。

b)使用Huffman对压缩后的质量分数序列进行解压，得到转码后的质量分数序列。

c)对压缩后的碱基序列进行解压，得到转码后的碱基序列。

d)根据映射表，恢复碱基序列和质量分数序列中特殊字符碱基和质量分数。

综上所述，本发明提出的针对碱基序列的编码方式使得包含特殊字符的碱基序列仍然可用2-bit的编码方式；此外，通过采用静态的huffman树对质量分数序列进行压缩，避免了在数据集上进行统计。本发明在压缩/解压速度足够快的前提下，压缩效率也比较高，并且每条测序序列的记录能够独立压缩/解压，不妨碍Mark Duplicate等算法的细粒度并行。

在本文的描述中，提及的原质量分数、原碱基等是指基因测序技术产生的在经过编码、转换、压缩等操作之前的数据。

需要说明的是，虽然上文按照特定顺序描述了各个步骤，但是并不意味着必须按照上述特定顺序来执行各个步骤，实际上，这些步骤中的一些可以并发执行，甚至改变顺序，只要能够实现所需要的功能即可。

本发明可以是系统、方法和/或计算机程序产品。计算机程序产品可以包括计算机可读存储介质，其上载有用于使处理器实现本发明的各个方面的计算机可读程序指令。

计算机可读存储介质可以是保持和存储由指令执行设备使用的指令的有形设备。计算机可读存储介质例如可以包括但不限于电存储设备、磁存储设备、光存储设备、电磁存储设备、半导体存储设备或者上述的任意合适的组合。计算机可读存储介质的更具体的例子(非穷举的列表)包括：便携式计算机盘、硬盘、随机存取存储器(RAM)、只读存储器(ROM)、可擦式可编程只读存储器(EPROM或闪存)、静态随机存取存储器(SRAM)、便携式压缩盘只读存储器(CD-ROM)、数字多功能盘(DVD)、记忆棒、软盘、机械编码设备、例如其上存储有指令的打孔卡或凹槽内凸起结构、以及上述的任意合适的组合。

以上已经描述了本发明的各实施例，上述说明是示例性的，并非穷尽性的，并且也不限于所披露的各实施例。在不偏离所说明的各实施例的范围和精神的情况下，对于本技术领域的普通技术人员来说许多修改和变更都是显而易见的。本文中所用术语的选择，旨在最好地解释各实施例的原理、实际应用或对市场中的技术改进，或者使本技术领域的其它普通技术人员能理解本文披露的各实施例。

Claims (9)

1.一种基因数据的编码方法，包括：

步骤1：对基因数据中AGCT字符及其原质量分数进行编码；

步骤2：基于预定的编码转换表对基因数据中除AGCT字符之外的特殊字符及其原质量分数进行编码，以获得编码之后的碱基序列和质量分数序列，包括：

步骤21：将所述特殊字符的原质量分数转化为有效质量分数；

步骤22：基于所述预定的编码转换表确定所述特殊字符及其有效质量分数所对应的编码转换之后的编码，其中，所述预定的编码转换表存储的是所述特殊字符和其有效质量分数的组合在编码转换之后与AGCT字符和质量分数的组合的对应关系。

2.根据权利要求1所述的方法，其中，所述特殊字符的原质量分数的范围是[33-126]，在编码转换之后，所述特殊字符的质量分数的范围是[0-32，127]。

3.根据权利要求1所述的方法，其中，所述有效质量分数的范围为[0-10]。

4.根据权利要求1所述的方法，其中，在步骤1和步骤2中，以基因数据中的一条基因测序序列为存取单元进行编码。

5.根据权利要求1所述的方法，其中，所述特殊字符包括R、Y、K、M、S、W、B、D、H、V、N和占位符。

6.根据权利要求1所述的方法，在步骤2之后还包括：

步骤3：将所述编码之后的质量分数序列存储为差值形式的序列；

步骤4：对所述差值形式的序列基于哈夫曼编码方式进行编码压缩。

7.一种基因数据的处理方法，包括：

根据权利要求1至6任一项所述的方法对基因数据进行编码，以获得编码之后的碱基序列和质量分数序列；

在spark环境中，使用预定的哈夫曼树对所述编码之后的质量分数序列进行压缩。

8.一种计算机可读存储介质，其上存储有计算机程序，其中，该程序被处理器执行时实现根据权利要求1至7中任一项所述方法的步骤。

9.一种计算机设备，包括存储器和处理器，在所述存储器上存储有能够在处理器上运行的计算机程序，其特征在于，所述处理器执行所述程序时实现权利要求1至7中任一项所述的方法的步骤。

Images