CN105760706B - 一种二代测序数据的压缩方法 - Google Patents

Abstract

本发明公开了一种二代测序数据的压缩方法。方法包括：对每个样本的二代测序数据，根据第一预设长度进行分割，生成BSSL初始文件；根据第二预设长度构建切割标签文件；根据切割标签文件对BSSL初始文件进行处理，获得BSSL中间文件；将BSSL中间文件合并得到BSSL最终文件；统计BSSL最终文件中种子序列的频次分布结果，根据结果得到种子文件；结合测序数据的格式特征，确定压缩规则，基于种子文件，对每个样本的二代测序数据进行压缩。本发明通过将二代测序数据进行分割以并行处理，提高处理速度，同时结合种子序列选择并得到种子文件，根据二代测序数据的格式特征和种子文件对二代测序数据进行压缩，从而可将二代测序数据的存储空间大大缩小。

Description

一种二代测序数据的压缩方法

技术领域

本发明涉及生物信息和数据压缩技术领域，尤其涉及一种二代测序数据的压缩方法。

背景技术

DNA（脱氧核糖核酸）是细胞中用以长期存储生物遗传指令信息的双螺旋长链状聚合物。它是由四种碱基：腺嘌呤（A）、胸腺嘧啶（T）、鸟嘌呤（G）和胞嘧啶（C）两两配对而成的碱基对序列。

随着人类基因组计划等大型国际合作研究项目的实施，一系列包括基因组学、转录组学、RNA（核糖核酸）组学和蛋白质组学等在内的研究导致了海量数据的产生，进而对数据存储和传输提出了更多的挑战。数据压缩有利于节约存储空间，提高数据交换及网络传输效率，这一点对于海量生物信息数据来说同样重要。

到目前为止，已有针对生物数据的压缩而开发的算法和工具，如Gencompress、DNAcompress、RNACompress等。但是，这些算法和工具的压缩比例不高。

发明内容

本发明提供一种二代测序数据的压缩方法，包括：

初文件生成步骤：对每个样本的二代测序数据，根据第一预设长度进行分割，生成固定长度碱基序列BSSL初始文件；

标签文件生成步骤：根据第二预设长度构建切割标签文件；

中间文件生成步骤：根据所述切割标签文件对所述BSSL初始文件进行处理，获得以切割标签为界的BSSL中间文件；

终文件生成步骤：将得到的所有所述BSSL中间文件合并，得到BSSL最终文件；

种子文件生成步骤：对所述BSSL最终文件，统计种子序列的频次分布结果，根据所述频次分布结果对所述BSSL最终文件进行选择，得到种子文件；

压缩步骤：结合二代测序数据的格式特征，确定压缩规则，基于所述种子文件并根据所述压缩规则，对所述每个样本的二代测序数据进行压缩。

进一步地，所述初文件生成步骤包括：根据第一预设长度，将每个样本的二代测序数据分割成多个小文件；对每个所述小文件，按预设序列固定长度生成BSSL初始文件。

进一步地，所述标签文件生成步骤包括：根据第二预设长度，编写切割标签生成程序，所述切割标签是长度固定的短的碱基序列的全集合；运行所述切割标签生成程序，得到切割标签文件。

进一步地，所述中间文件生成步骤包括：根据所述切割标签文件，编写中间文件生成程序；运行所述中间文件生成程序，得到和切割标签相关的多个文件夹；将所述多个文件夹中的文件合并，排序成以各个文件夹名称命名的多个文件，所述多个文件为所述BSSL中间文件。

进一步地，所述终文件生成步骤包括：将所述BSSL中间文件按次序合并，得到一个所述BSSL最终文件，其中所述BSSL最终文件中的每条序列为BSSL条目。

进一步地，所述种子文件生成步骤包括：统计所述BSSL最终文件中种子序列的出现次数；根据统计结果从所述BSSL最终文件中选出预设数目个所述BSSL条目，选择结果构成所述种子文件。

进一步地，所述压缩步骤包括：根据二代测序数据的格式特征，编写压缩程序，其中，对于二代测序数据的第一行，压缩时丢弃“*”开始后的所有信息，只保留序列的标识符；对于二代测序数据的第二行，丢弃非确定碱基出现在非末尾位置的读段序列，利用ASCII码进行编码；对于测序数据的第三行，不进行压缩处理；对于测序数据的第四行，将碱基质量值大于第一预设阈值的质量值全部处理为等于第一预设阈值，将碱基质量值小于第二预设阈值的质量值全部处理为第二预设阈值减去第一预设阈值，并将具有相同字符多于预设数量个的字符修改为该字符及其相同字符的个数的形式。

其中，所述利用ASCII码进行编码包括：利用ASCII码对所述种子文件进行编码，并将所述每个样本的二代测序数据中的第二行替换为编码后的种子文件，如果存在不能采用编码后的种子文件替换的序列，则将该序列中的每三个碱基替换为一个ASCII码。

本发明的有益效果是：通过将二代测序数据进行分割，并行处理分割后得到的文件，提高处理速度，同时结合种子序列选择并得到种子文件，然后根据二代测序数据的格式特征和种子文件对二代测序数据进行压缩，从而可将二代测序数据的存储空间大大缩小。

附图说明

为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍。

图1为fastq格式的一种示例的示意图；

图2为本发明一种实施例的二代测序数据压缩方法的流程示意图；

图3至图8分别为采用ASCII码进行编码的示意图。

具体实施方式

fastq格式是一种用于存储生物序列（主要是核酸序列）及其测序质量分数的基于文本的格式，目前已成为高通量测序仪输出序列的标准格式。一个fastq文件中每个序列通常由四行组成，如图1的示例所示：第一行（line 1）以“@”开头，后面接序列的标识符和可选的描述信息；第二行（line 2）是测序的序列；第三行（line 3）以“+”开头，后面可接相同于第一行的序列的描述信息；第四行（line 4）是每个碱基的测序质量分数，质量分数用相应数值的ASCII码表示。显然，fastq数据的reads（读段，指测序出来的一条条碱基序列）主要由A、T、G、C四种碱基构成（部分存在没有测定且尚不明确的碱基N）。对于一条reads，以固定的碱基长度对其做分割，由该固定碱基长度组成的碱基序列（BSSL，base sequence withset length，简称固定长度碱基序列）便称为一个BSSL条目。由于高深度基因组fastq数据相当庞大，存在几亿至几十亿条reads，对于某些BSSL条目，在整个基因组fastq数据中，重复出现次数可高达几十万至几百万次。而对不同的碱基长度值，各BSSL条目的重复频次也会有所差异。

本发明利用fastq数据中某些BSSL条目的高频次重复出现的特性，以千人基因组fastq数据为示例进行BSSL分析，寻找最佳碱基长度值，构建合适的BSSL条目高重复频次表（种子文件，seed），建立相关的压缩算法，将可使fastq存储空间大大缩小。具体地，提出了一种二代测序数据的压缩方法，如图2所示，包括步骤S11~S16。

步骤S11用于生成初文件，即，对每个样本的二代测序数据，根据第一预设长度进行分割，生成BSSL初始文件。在该步骤S11中，首先根据第一预设长度，将每个样本的二代测序数据分割成多个小文件，然后对每个小文件，按预设序列固定长度生成BSSL初始文件。第一预设长度可以根据经验值来设定，也可以通过若干次试验确定出，例如，首先设定为具有一定梯度的预设长度，比如10、20、30、40、50等，然后比较最终生成的BSSL文件（即步骤S14中的BSSL最终文件），如果BSSL文件大小越小，说明对应的第一预设长度越优，由此可以确定出第一预设长度。

步骤S12用于生成标签文件，即根据第二预设长度构建切割标签文件。在步骤S12中，根据第二预设长度，编写切割标签生成程序，该切割标签是指长度固定的短的碱基序列的全集合，然后运行切割标签生成程序，得到切割标签文件。切割标签文件是由碱基A、T、C、G构成，如果其长度为5，则一共有4^5=1024个标签；切割标签文件的作用是用来有序分割BSSL初始文件。

步骤S13用于生成中间文件，即根据切割标签文件对步骤S11得到的BSSL初始文件进行处理，获得以切割标签为界的BSSL中间文件，即该文件以标签为序列的起点和终点。在步骤S13中，根据步骤S12得到的切割标签文件，编写中间文件生成程序，然后运行该中间文件生成程序，得到和切割标签相关的多个文件夹，将这些多个文件夹中的文件合并，排序成以各个文件夹名称命名的多个文件，该多个文件即为BSSL中间文件。

步骤S14用于生成终文件，即将得到的所有所述BSSL中间文件合并从而得到BSSL最终文件。在步骤S14中，将步骤S13得到的BSSL中间文件按次序合并，得到一个BSSL最终文件，其中该BSSL最终文件中的每条序列称为BSSL条目。

步骤S15用于生成种子文件，即对步骤S14得到的BSSL最终文件，统计种子序列的频次分布，根据统计的频次分布结果对BSSL最终文件进行选择，选出预设数目个BSSL条目，选择结果构成种子文件。种子文件是按照序列的频次由高到低进行排序，取排在最前的预设数目个种子序列作为种子文件。

步骤S16用于压缩，即结合二代测序数据的格式特征，确定压缩规则，基于步骤S15得到的种子文件以及确定出的压缩规则，对每个样本的二代测序数据进行压缩。在步骤S16中，根据二代测序数据的格式特征，编写压缩程序，其中，对于二代测序数据的第一行，压缩时丢弃“*”开始后的所有信息，只保留序列的标识符；对于二代测序数据的第二行，丢弃非确定碱基出现在非末尾位置的读段序列，利用ASCII码进行编码；对于测序数据的第三行，不进行压缩处理；对于测序数据的第四行，将碱基质量值大于第一预设阈值的质量值全部处理为等于第一预设阈值，将碱基质量值小于第二预设阈值的质量值全部处理为第二预设阈值减去第一预设阈值，并将具有相同字符多于预设数量个的字符修改为该字符及其相同字符的个数的形式。其中，在利用ASCII码进行编码时，首先利用ASCII码对种子文件进行编码，并将每个样本的二代测序数据中的第二行替换为编码后的种子文件，如果存在不能采用编码后的种子文件替换的序列，则将该序列中的每三个碱基替换为一个ASCII码。

下面通过具体实施方式结合附图对本发明作进一步详细说明。

实施例1

本实施例以千人基因组的数据为例进行说明，其中样本NA12345是千人基因组的一千多个样本中的一个样本数据，这里为方便描述，以NA12345作为其中的一个例子来举例说明，该样本数据的二代测序数据存储为fastq格式，对应的文件名为example.fastq。以下利用前述步骤S11~S16对千人基因组的二代测序数据进行压缩。

本实施例中，步骤S11生成BSSL初文件。具体如下。

步骤S11中，首先用split命令将example.fastq分割成80000000行（即前述的第一预设长度，当然，该第一预设长度也可以是其它数值）的多个小文件；系统可以将自动为结果小文件命名，如第一个文件将被命名为exampleaa.fastq。split命令是linux系统下用于切割文件的一个命令。当然也可以采用脚本语言例如Perl来编写程序，以将大文件example.fastq分割成多个小文件。将大文件分割成小文件的目的是为了并行处理这些小文件，从而可以缩短程序运行时间，提高效率。

接着，编写BSSL初文件生成程序。在本实施例中，采用脚本语言Perl进行程序编写，其它实施例中本领域技术人员可以采用自己熟悉的编程语言进行编写。本实施例的后续步骤比如S12~S16中都是采用脚本语言Perl进行程序编写，同样地，其它实施例中本领域技术人员可以采用自己熟悉的编程语言进行编写。对于本实施例，编写的程序为例如scriptpath /kmer_n2.pl。这是用来预处理待压缩文件（如example.fastq），kmer_n2.pl是用来将example.fastq文件转换为BSSL初始文件的程序，当然kmer_n2.pl会处理样本数据中的所有fastq格式的小文件，而不仅仅是exampleaa.fastq。然后运行程序命令，例如perlkmer_n2.pl –i exampleaa.fastq –n NA12345 –k 17 –d ./kmer/primary_BSSL。显然程序基本运行模式为perl kmer_n2.pl -i fast_file_name -n sample_ID -k seq_length-d out_dir1。其中fast_file_name为待处理的fastq格式的小文件，sample_ID表示样本数据如NA12345，-k选项为BSSL初始文件的序列固定长度，这里举例为17，但可以修改为其它数值。-i、-n、-d等选项都是对应的程序kmer_n2.pl自带的，可以参考相关资料得到，在此不作详述。

执行步骤S11后，将得到包括NA12345的所有千人基因组的样本的BSSL初始文件，为方便下文描述，这里假设其中一个文件名称为primary_01.bssl。

本实施例中，步骤S12构建BSSL文件切割标签list。具体如下。

步骤S12中，首先编写切割标签生成程序，例如编写为： scriptpath/New_idea_1/bin/generate_stop_character.pl。

然后运行程序命令：perl generate_stop_character.pl > stop_c-kmer17.list，获得切割标签文件 stop_c-kmer17.list。该切割标签是长度固定的短的碱基序列的全集合，本例中切割标签长度设置为5（即第二预设长度），当然也可以考虑设置其他的标签长度。generate_stop_character.pl是形成一个分割标签文件，后续将进一步按照这个分割标签对步骤S11得到的BSSL初文件进行分割。

本实施例中，步骤S13用于获得大量以list为界的BSSL中间文件。具体如下。

步骤S13中，首先编写BSSL中间文件生成程序，例如scriptpath /New_idea_all/bin/classify_by_sc.pl。然后运行程序命令：perl classify_by_sc.pl -i primary_01.bssl -o out_dir2 -l stop_c-kmer17.list。执行程序命令结束后，对于各样本而言，在其结果文件夹out_dir2中将存在1024个文件夹。这里，由于采用的标签是由ATCG四个碱基进行自由组合形成的，4的5次幂即是1024。接着，再编写程序，用以分别将1024个文件夹中各文件合并，并排序成一个以各文件夹名称命名的文件，例如scriptpath /New_idea_all/bin/kmer_merge_part.pl，然后再运行该程序命令，例如perl kmer_merge_part.pl -i out_dir2 -o out_dir3 -l stop_c-kmer17.list。该程序的作用是将1024个文件夹变成1024个文件，即将每个文件夹里的文件（按照切割标签形成的文件夹）进行合并，合并后的文件以对应的文件夹名来命名。由此，得到了BSSL中间文件。

本实施例中，步骤S14用于得到BSSL最终文件。在该步骤中，将步骤S13中获得的BSSL中间文件（此处即1024个文件）按次序合并便得到BSSL文件，文件名称为例如example_all.bssl。其中BSSL文件中每条序列（seq）称为一个BSSL条目。

本实施例中，步骤S15用于获得种子文件（seed）。

步骤S15中，首先统计BSSL最终文件中BSSL条目的频次（频次1）的频次（频次2）分布，将结果文件按频次1为键，按数字由大到小排序；然后依据该频次分布结果选出频次1最大的658个BSSL条目，即为种子文件，即一个种子文件由658（即前述的预设数目，该数目的来源将在下文中提及）个BSSL条目组成。这里的“结果文件”是指BSSL最终文件（即步骤S14提到的将1024个BSSL中间文件合并而得）。频次1是种子序列，频次2对对应的种子序列的出现次数。种子序列（seed sequence）的概念类似本领域中提及的种子序列，就是一些比较可信或者比较有代表性的序列集合，在这些序列的基础上，进行一些序列比对，可以获取更多的目标序列。例如在重复序列检索中，有些seed sequeces可以是基因组中重复次数很高的短片段，对这些短片段与基因组其他序列做blast等比对，获得一些比这些短片段长度更长，但是重复次数较少的序列，实现重复序列检索的目的。

在本实施例中，步骤S16用于制定压缩解压方案并编写程序对fastq文件进行压缩。

在步骤S16中是结合二代测序数据的格式特征来确定压缩规则，具体如下。

对于二代测序数据的第一行（如图1的line 1）：压缩时丢弃 “*”开始后的所有信息；同时，由于各reads的该行信息呈现@SRR12345.1、@SRR12345.2、 @SRR12345.3 的规律，故压缩时仅在压缩文件头保留@SRR12345样信息；压缩后文件内部不再保留line1。

对于二代测序数据的第二行（如图1的line 2）：首先，若reads中间出现N，该reads将被丢弃，仅末尾出现N的reads将被保留；接着，充分利用ASCII码，本实施例选定的ASC码为33~126，共94个ASC码/字符；其中AAA→TTT(4^3)占据chr(33)→chr(96)共64个字符，[ATGC][ATGC]N (4^2)占据chr(97)--->chr(112)共16个字符，[ATGC]NN (4^1)占据chr(113)→chr(116)共4个字符；N、NN、NNN分别用chr(117)、 chr(118)、 chr(119)来表示；余下7个字符（chr(120)→chr(126)）将与chr(3)→chr(126)联用来表示种子（seed）文件的kmer条目。这里对此解释下。由于一共有A、T、C、G四种碱基，当AAA→TTT时，长度为3，一共需要4的3次方，即64个ASCII码来替换；当这三个碱基中有一个为N（即未知序列，又称未确定碱基）时，还剩下两个碱基，一共需要4的2次方来表示；当有两个N时，同理需要4个ASCII码来替换。由此得到了如图3至图8所示的表格，其中，图3至图6表示碱基对应的ASCII编码，图8表示种子文件，图7表示利用ASCII码对图8的种子文件进行转换后的表示形式。当kmer条目为双字符编码时，受ASC码个数限制，seed文件kmer条目数最多为7*94，即658个，即前述的种子文件中涉及的BSSL条目个数。

对于line 2，在压缩时，优先编码与seed条目（即种子文件中的BSSL条目）匹配的序列块，然后再将不能匹配序列按照上述规则做3to1编码，不足3个碱基的，以N补全，若是序列以N或NN结尾，则无需补全。这里的规则是指，将待压缩文件（fastq）的序列利用种子文件的ASCII码来替换，假如种子文件不存在时，用3to1编码规则来替换。3to1规则是指将原始文件中的3个碱基用一个ASCII码来替换，从而实现压缩，显然压缩比为3。

对于二代测序数据的第三行（如图1的line 3）：保持原形，不做处理。

对于二代测序数据的第四行（如图1的line 4）：根据fastq序列本身的特点，质量值压缩规定为：若碱基质量值>=20，则全部视为质量值==20处理；若质量值<=-7,则全按质量值==-27处理；再将超过2个字符相同的字符（如BBBB）压缩成“B4”形式。这里20和-7表示碱基对应的测序质量值，该值越低，说明测序结果越不可信。

可见，本发明实施例基于千人基因组的fastq数据，结合其格式特征，通过设定不同的碱基长度值来切割其碱基序列，统计出每个BSSL的出现频次，找到最适的碱基长度值，建立seed文件；基于该seed文件，通过一定的算法处理，使得能够最大化压缩fastq数据中的碱基序列；同时，对fastq数据的序列标签和测序质量值进行有选择的无损和有损压缩处理，使得压缩比最大化。

与现有的常见压缩软件相比，本发明实施例基于fastq数据固有特性，与目前的常见压缩软件相比，压缩比显著提高，同时，经测试，本发明实施例可与已知常见的压缩软件结合使用，使得压缩比大大提高。随着生物数据的爆炸性增长，存储空间成本巨大，而其中fastq数据占得比例很高，通过对fastq数据进行压缩后存储，能大大降低存储空间，提高存储效率。

本领域技术人员可以理解，上述实施方式中各种方法的全部或部分步骤可以通过程序来指令相关硬件完成，该程序可以存储于一计算机可读存储介质中，存储介质可以包括：只读存储器、随机存储器、磁盘或光盘等。

以上内容是结合具体的实施方式对本发明所作的进一步详细说明，不能认定本发明的具体实施只局限于这些说明。对于本发明所属技术领域的普通技术人员来说，在不脱离本发明构思的前提下，还可以做出若干简单推演或替换。

Claims (8)

1.一种二代测序数据的压缩方法，其特征在于，包括：

初文件生成步骤：对每个样本的二代测序数据，根据第一预设长度进行分割，生成固定长度碱基序列初始文件；

标签文件生成步骤：根据第二预设长度构建切割标签文件；

中间文件生成步骤：根据所述切割标签文件对所述固定长度碱基序列初始文件进行处理，获得以切割标签为界的固定长度碱基序列中间文件；

终文件生成步骤：将得到的所有所述固定长度碱基序列中间文件合并，得到固定长度碱基序列最终文件；

种子文件生成步骤：对所述固定长度碱基序列最终文件，统计种子序列的频次分布结果，根据所述频次分布结果对所述固定长度碱基序列最终文件进行选择，得到种子文件；

压缩步骤：结合二代测序数据的格式特征，确定压缩规则，基于所述种子文件并根据所述压缩规则，对所述每个样本的二代测序数据进行压缩。

2.如权利要求1所述的方法，其特征在于，所述初文件生成步骤包括：

根据所述第一预设长度，将每个样本的二代测序数据分割成多个小文件；

对每个所述小文件，按预设序列固定长度生成所述固定长度碱基序列初始文件。

3.如权利要求1所述的方法，其特征在于，所述标签文件生成步骤包括：

根据所述第二预设长度，编写切割标签生成程序，所述切割标签是长度固定的短的碱基序列的全集合；

运行所述切割标签生成程序，得到所述切割标签文件。

4.如权利要求1所述的方法，其特征在于，所述中间文件生成步骤包括：

根据所述切割标签文件，编写中间文件生成程序；

运行所述中间文件生成程序，得到和切割标签相关的多个文件夹；

将所述多个文件夹中的文件合并，排序成以各个文件夹名称命名的多个文件，所述多个文件为所述固定长度碱基序列中间文件。

5.如权利要求1所述的方法，其特征在于，所述终文件生成步骤包括：

将所述固定长度碱基序列中间文件按次序合并，得到一个所述固定长度碱基序列最终文件，其中所述固定长度碱基序列最终文件中的每条序列为固定长度碱基序列条目。

6.如权利要求1所述的方法，其特征在于，所述种子文件生成步骤包括：

统计所述固定长度碱基序列最终文件中种子序列的频次；

根据统计结果从所述固定长度碱基序列最终文件中选出预设数目个所述固定长度碱基序列条目，选择结果构成所述种子文件。

7.如权利要求1所述的方法，其特征在于，所述压缩步骤包括：

根据二代测序数据的格式特征，编写压缩程序，其中，

对于二代测序数据的第一行，压缩时丢弃“*”开始后的所有信息，只保留序列的标识符；

对于二代测序数据的第二行，丢弃非确定碱基出现在非末尾位置的读段序列，利用ASCII码进行编码；

对于测序数据的第三行，不进行压缩处理；

对于测序数据的第四行，将碱基质量值大于第一预设阈值的质量值全部处理为等于第一预设阈值，将碱基质量值小于第二预设阈值的质量值全部处理为第二预设阈值减去第一预设阈值，并将具有相同字符多于预设数量个的字符修改为该字符及其相同字符的个数的形式。

8.如权利要求7所述的方法，其特征在于，所述利用ASCII码进行编码包括：

利用ASCII码对所述种子文件进行编码，并将所述每个样本的二代测序数据中的第二行替换为编码后的种子文件，如果存在不能采用编码后的种子文件替换的序列，则将该序列中的每三个碱基替换为一个ASCII码。