高通量测序数据质量过滤方法和过滤装置与流程

本发明涉及生物和计算机领域，特别涉及通过大数据技术对高通量测序数据进行质量过滤的方法和装置。

背景技术：

基于高通量测序数据寻找和致病基因、癌症治疗、个性化用药相关的染色体突变位点为临床应用提供了不可估量的前景。由于测序技术的不断进步，获取到的数据越来越多，如何快速地处理不断增加的高通量测序数据已成为亟待解决的问题。

在获取到高通量测序数据后，需要对原始数据进行过滤，将低质量的数据过滤掉，保留高质量的数据进行下一步处理。目前常用的过滤软件Trime对1.1G X 2的高通量测序数据进行双端过滤需要九分钟以上，随着数据的不断增大，质量过滤需要的时间也越来越长。如何缩短质量过滤消耗的时间，让高质量的数据尽快地进入下游的分析环节，从而缩短科研人员、患者等待分析结果的时间，已成为亟待解决的问题。

技术实现要素：

有鉴于此，本发明基于分布式计算框架提供了一种对高通量测序数据进行质量过滤的方法和装置，能够将低质量的数据过滤掉，与以往的过滤方法相比，大大提升了处理速度。

本发明的实施例提供了一种对高通量测序数据进行质量过滤的方法，所述方法包括：

根据所述高通量测序数据为并行计算做准备；

通过并行计算过滤掉准备好的数据中质量不达标的数据。

优选地，所述根据所述高通量测序数据为并行计算做准备包括：

根据所述高通量测序数据中的质量数据确定质量转换方式；

对已确定质量转换方式的高通量测序数据进行切分；

生成对切分后的数据进行并行计算的执行实体。

优选地，所述对已确定质量转换方式的高通量测序数据进行切分包括：

将包含已确定质量转换方式的高通量测序数据的第一文件和第二文件分别转换为各自对应的第一RDD和第二RDD；

将第一RDD和第二RDD分别切分为各自对应的第一partition组和第二partition组；

根据第一文件和第二文件中对应的数据将第一RDD和第二RDD合并为第三RDD；

将第三RDD切分为第三partition组。

优选地，所述生成对切分后的数据进行并行计算的执行实体为：生成对第三partition组进行并行计算的执行实体task。

优选地，其特征在于，所述通过并行计算过滤掉准备好的数据中质量不达标的数据包括：

根据预定质量值阈值和质量值转换方式通过执行实体对所述高通量测序数据并行地进行过滤；

根据预定序列长度阈值通过执行实体对保留下来的高通量测序数据并行地进行过滤。

优选地，所述根据预定质量值阈值、所述质量值转换方式对所述高通量测序数据进行过滤包括：

如果所述高通量测序数据某一记录中的质量行中的某一位置上的质量值小于预定质量值阈值，则通过执行实体并行地过滤掉所述质量行该位置及以后的数据，以及同一记录中的序列行中的对应位置及以后的数据。

优选地，所述根据预定序列长度阈值对保留下来的高通量测序数据进行过滤包括：

在保留下来的高通量测序数据中，如果第一文件和第二文件相对应的两个记录中有任意一个记录中的序列行长度小于预定长度阈值，则通过执行实体并行地过滤掉第一文件和第二文件中相对应的这两个记录。

另一方面，本发明的实施例还提供了一种对高通量测序数据进行质量过滤的装置，所述装置包括：

并行准备模块，用于根据所述高通量测序数据为并行计算做准备；

质量过滤模块，用于通过并行计算过滤掉准备好的数据中质量不达标的数据。

优选地，所述并行准备模块包括：

质量转换方式确定单元：用于根据所述高通量测序数据中的质量数据确定质量转换方式；

数据切分单元：用于对已确定质量转换方式的高通量测序数据进行切分；

执行实体生成单元：用于生成对切分后的数据进行并行计算的执行实体。

优选地，所述数据切分单元具体用于：

将包含已确定质量转换方式的高通量测序数据的第一文件和第二文件分别转换为各自对应的第一RDD和第二RDD；

将第一RDD和第二RDD分别切分为各自对应的第一partition组和第二partition组；

根据第一文件和第二文件中对应的数据将第一RDD和第二RDD合并为第三RDD；

将第三RDD切分为第三partition组。

优选地，所述执行实体生成单元具体用于：生成对第三partition组进行并行计算的执行实体task。

优选地，所述质量过滤模块具体用于：

根据预定质量值阈值和质量值转换方式通过执行实体对所述高通量测序数据并行地进行过滤；

根据预定序列长度阈值通过执行实体对保留下来的高通量测序数据并行地进行过滤。

优选地，所述质量过滤模块用于根据预定质量值阈值和质量值转换方式通过执行实体对所述高通量测序数据并行地进行过滤包括：

如果所述高通量测序数据某一记录中的质量行中的某一位置上的质量值小于预定质量值阈值，则通过执行实体并行地过滤掉所述质量行该位置及以后的数据，以及同一记录中的序列行中的对应位置及以后的数据。

优选地，所述质量过滤模块用于根据预定质量值阈值和质量值转换方式通过执行实体对所述高通量测序数据并行地进行过滤包括：

在保留下来的高通量测序数据中，如果第一文件和第二文件相对应的两个记录中有任意一个记录中的序列行长度小于预定长度阈值，则通过执行实体并行地过滤掉第一文件和第二文件中相对应的这两个记录。

本发明提供的技术方案采用大数据处理技术Spark集群开发了用于对高通量测序数据进行质量过滤的软件工具Sfastq_filter，极大地提高了对高通量测序数据进行质量过滤的速度：Sfastq_filter在12核18G的配置环境下双端过滤1.1GX2的数据，只需要三分钟，与传统的Trime软件相比处理速度快两倍以上。如果机器配置更高，数据量更大，那么对比效果也会更加明显。

附图说明

图1为本发明实施例一提供的对高通量测序数据进行质量过滤的方法的示意图；

图2为本发明实施例二提供的对高通量测序数据进行质量过滤的方法的示意图；

图3为本发明实施例三提供的对高通量测序数据进行质量过滤的方法的示意图；

图4为本发明实施例四至六提供的对高通量测序数据进行质量过滤的装置的示意图。

具体实施方式

为使本发明的目的、技术方案和优点更加清楚，以下将参照本发明实施例中的附图，通过实施方式清楚、完整地描述本发明的技术方案，显然，所描述的实施例是本发明一部分实施例，而不是全部的实施例。

在本发明的实施例中，相关名词解释如下：

高通量测序数据：用高通量测序方法获得的数据；

高通量测序：能够一次并行对几十万到几百万条DNA片段进行序列测定的测序方法；

第一阈值：在本发明的实施例中为58，在实际应用中可以根据具体情况取适当的值；

第二阈值：在本发明的实施例中为76，在实际应用中可以根据具体情况取适当的值；

第一碱基质量值转换方式：将质量行中每一位置上的质量数据(即ASCII码字符的ASCII码值)减去第一阈值，得到对应碱基的质量值；

第二碱基质量值转换方式：将质量行中每一位置上的质量数据(即ASCII码字符的ASCII码值)减去第二阈值，得到对应碱基的质量值。

第一文件、第二文件：作为Sfastq_filter软件的输入文件成对出现的高通量测序数据文件；其中第一文件和第二文件中的记录是一一对应的。

质量达标：如果第一文件和第二文件中一一对应的两个记录中的质量行某一位置及该位置以前的质量数据的质量值均大于等于预定质量值阈值，且这两个一一对应的记录中的序列行的碱基数均大于等于预定序列长度阈值，那么满足上述条件的高通量测序数据的记录称为质量达标。

实施例一

请参阅图1，在本发明的第一个实施例中，在根据测序后得到的高通量测序数据为并行计算做好准备后，对数据进行并行计算，过滤掉其中低质量的数据。

S101、根据高通量测序数据为并行计算做准备。

在包含高通量测序数据的FastQ文件中，每一个记录包括四行，分别为：

以“@”开头后面附加测序介绍信息的标识行；

由A、T、G、C四种碱基组成的序列行(测序仪无法识别的碱基用N表示)；

“+”行(或者“+”后面附带标识行中@后面的内容，但该内容一般被省略)；

由ASCII码字符组成的质量行(质量行和序列行长度相等且质量行中的ASCII码字符与序列行中的碱基一一对应，质量行中每一位置的质量数据为该位置碱基质量值的ASCII码字符表示方式)。

用于Sfastq_filter软件的输入文件是成对出现的，包括第一文件和第二文件两个FastQ文件，其中第一文件和第二文件中的记录是一一对应的。将包含原始高通量测序数据的第一文件和第二文件切分为多个数据块，为并行计算做准备。

S102、通过并行计算过滤掉准备好的数据中质量不达标的数据。

根据预定阈值对准备好的数据进行过滤，去掉其中低质量的数据。预定阈值包括预定质量值阈值和预定序列长度阈值。

根据预定质量值阈值同时过滤第一文件和第二文件中一一对应的两个记录。

在经过预定质量值阈值过滤后保留下来的数据中，根据预定序列长度阈值同时过滤第一文件和第二文件中一一对应的两个记录。

将经过预定质量值阈值和预定序列长度阈值过滤后保留下来的数据根据其原始数据来源分别输出到第一文件和第二文件各自对应的过滤结果文件中。

实施例二

请参阅图2，在本发明的第二个实施例中，使用Hadoop并行计算框架对高通量测序数据进行并行计算，从而过滤掉低质量的数据。

在本实施例中，相关名词解释如下：

Hadoop：由Apache基金会所开发的分布式并行计算框架。

HDFS(Hadoop Distributed File System)：由Hadoop实现的一个分布式文件系统。

S201、根据所述高通量测序数据中的碱基质量数据确定碱基质量值转换方式。

输入包含原始高通量测序数据的FastQ文件。用于Sfastq_filter软件的输入文件是成对出现的，包括第一文件和第二文件两个FastQ文件，其中第一文件和第二文件中的记录是一一对应的。根据读入的FastQ文件质量行中的碱基质量数据(即该位置碱基质量值的ASCII码字符表示方式)确定本文件对应的碱基质量值转换方式：

如果读取到大于第一阈值且小于等于第二阈值的碱基质量数据，则忽略此值，继续读入下一位置的碱基质量值；

如果读取到小于等于第一阈值的碱基质量数据，则确定本文件对应第一碱基质量值转换方式，选择碱基质量值转换方式的过程结束；

如果读取到大于第二阈值的碱基质量数据，则确定本文件对应第二碱基质量值转换方式，选择碱基质量值转换方式的过程结束。

S202、对已确定碱基质量值转换方式的高通量测序数据进行切分。

HDFS将输入的FastQ文件根据一定的规则切分成小数据块并保存。切分规则如下：例如输入文件为3G，在Hadoop中将一个数据块的大小设置为128M，那么输入文件总共将被切分为3*1024/128＝24块。

S203、生成包含切分后的数据的执行实体。

在hadoop集群中，参与并行计算的多台计算机并行地读取HDFS中的小数据块，并启动map job和reduce job，job为每一个小数据块生成一个map task。Map task的计算结果存储在中间结果文件中，中间结果文件保存在HDFS中。Reduce job 从HDFS中读取中间结果文件，并根据用户指定的数量生成多个reduce task。Map task和reduce task是并行计算的执行实体。

S204、通过并行计算过滤掉准备好的数据中质量不达标的数据。

Hadoop在参与并行计算的多台计算机上并行地运行map task和reduce task，先运行map task，然后运行reduce task。

在hadoop集群中，通过多个map task进行并行计算：

首先根据第一文件和第二文件已确定的质量值转换方式对每个记录的质量行中的每一个位置上的碱基计算其碱基质量值。

然后根据预定质量值阈值和预定序列长度阈值对准备好的数据进行过滤，去掉其中的低质量数据，过滤过程如下：

同时遍历第一文件和第二文件中一一对应的两个记录中的质量行，如果质量行中某一位置上碱基的质量值小于预定质量值阈值时，则过滤掉该位置及其以后的所有数据；同时过滤掉相应记录的序列行在对应位置及其以后的所有数据；将标识行和“+”行原样输出。

在经过预定质量值阈值过滤之后保留的数据中，同时获取第一文件和第二文件中一一对应的两个记录中的序列行中的碱基数，如果这一对记录中的任何一个记录的碱基数小于预定序列长度阈值时，则将这一对记录全部过滤掉。

Map task的中间计算结果存储在HDFS上的中间结果文件中，reduce job读取到中间结果文件后，根据用户指定的数量生成多个reduce task。

在hadoop集群中，多个reduce task将经过预定质量值阈值和预定序列长度阈值过滤后保留下来的数据根据其原始数据来源分别输出到第一文件和第二文件各自对应的最终结果文件中。

实施例三

请参阅图3，在本发明的第三个实施例中，使用spark并行计算框架通过并行计算过滤掉高通量测序数据中质量不达标的数据。

S301、根据所述高通量测序数据中的碱基质量数据确定碱基质量值转换方式。

输入包含原始高通量测序数据的FastQ文件。用于Sfastq_filter软件的输入文件是成对出现的，包括第一文件和第二文件两个FastQ文件，其中第一文件和第二文件中的记录是一一对应的。在FastQ文件中，每一个记录包括四行，其中序列行由A、T、G、C四种碱基组成(测序仪无法识别的碱基用N表示)；质量行中每一位置上的质量数据为该位置的碱基质量值的ASCII码表示方式，质量行和序列行长度相等。

根据读入的FastQ文件质量行中的碱基质量值确定本文件对应的碱基质量值转换方式：

如果读取到大于第一阈值且小于等于第二阈值的碱基质量数据，则忽略此值，继续读入下一位置的碱基质量值；

如果读取到小于等于第一阈值的碱基质量数据，则确定本文件对应第一碱基质量值转换方式，选择碱基质量值转换方式的过程结束；

如果读取到大于第二阈值的碱基质量数据，则确定本文件对应第二碱基质量值转换方式，选择碱基质量值转换方式的过程结束。

S302、将包含已确定质量转换方式的高通量测序数据的第一文件和第二文件分别转换为各自对应的第一RDD和第二RDD。

在本实施例中：

Spark：是UC Berkeley AMPLab开发的一种计算框架。

RDD是指弹性分布式数据集(Resilient Distributed Datasets)，它是可容错的并行数据结构，使用户能够显式地在内存中保存中间的运算结果，通过控制RDD的分区来优化数据的布局，并使用丰富的转换算子进行操作。

在读取输入的FastQ文件时，spark将第一文件转换为第一RDD，将第二文件转换为第二RDD。

S303、将第一RDD和第二RDD分别切分为各自对应的第一partition组和第二partition组。

在本实施例中，partition是指spark在计算过程中，生成的数据在计算空间内的最小单元。

在生成RDD时，用户可以根据实际需要指定将RDD切分为partition的数量。例如输入文件为3G，设置将RDD切分为24个partition，那么每一个partition所占存储空间为3*1024/24＝128M。实际生成的partition的数量最少为(该文件所占存储空间/128M)，如果指定的partition数量少于(该文件所占存储空间/128M)，则实际将生成(该文件所占存储空间/128M)个partition。Spark根据用户指定的partition数量将RDD切分成partition，每个RDD对应的所有partition即为一个partition组。

在本实施例中，spark将第一RDD和第二RDD分别切分为各自对应的第一partition组和第二partition组。

S304、根据第一文件和第二文件中对应的数据将第一RDD和第二RDD合并为第三RDD。

由于第一文件和第二文件中的记录是一一对应的，在将此二者转换为第一RDD和第二RDD后，第一RDD和第二RDD中的记录也是一一对应的。而在通过并行计算对数据进行质量过滤时，需要同时遍历第一RDD和第二RDD中对应的记录，所以为了提高并行计算的速度，在进行并行计算以前，根据其中一一对应的记录，将第一RDD和第二RDD合并为第三RDD，用于后续的并行计算。

S305、将第三RDD切分为第三partition组。

在生成第三RDD时，可以根据实际需要指定将第三RDD切分为partition的数量，该数值独立于第一RDD和第二RDD分别对应的partition组中包含的partition数量。例如，原来第一RDD对应的第一partition组中有四个partition，第二RDD对应的第二partition组中有六个partition，那么在将第一RDD和第二RDD合并生成第三RDD后，可以将第三RDD切分为五个partition，即第三partition组中有五个partition。

S306、生成对第三partition组进行并行计算的执行实体task。

在本实施例中：

Job是指包在spark中由多个stage组成的并行计算，对RDD执行action操作后会生成job；

Stage是指在spark中，一个job会根据处理过程需要而分成不同的阶段即stage，stage由多个task组成；

Task是指被送到为某个应用启动的executor进程的工作单元。

在spark中对RDD进行action操作时生成DAG Scheduler(有向无环图调度器)，从而启动一个job。对一个job内的操作，根据处理过程需要分成不同的stage，并在每一个stage内产生一系列的task。通常一个RDD内的task数量与partition的数量相同。后续多个执行实体task将在多台计算机上对不同的partition执行并行计算过程。

在本实施例中，spark首先生成与第三RDD对应的job，再根据job的处理需要生成stage，并在每个stage中生成对第三partition组进行并行计算的多个执行实体task。后面将在spark集群的多台计算机上通过多个执行实体task并行地对每个partition进行后续计算。

S307、根据预定质量值阈值和质量值转换方式通过执行实体对所述高通量测序数据并行地进行过滤。

首先，计算数据的质量行中每一位置的碱基质量值：

如果本次并行计算数据对应的是第一碱基质量值转换方式，则将质量行中每一位置上的ASCII码字符的ASCII码值减去第一阈值即为对应的碱基质量值；

如果本次并行计算数据对应的是第二碱基质量值转换方式，则将质量行中每一位置上的ASCII码字符的ASCII码值减去第二阈值，即为对应的碱基质量值。

然后，同时遍历第一文件和第二文件中一一对应的两个记录中的质量行，如果质量行中某一位置上碱基的质量值小于预定质量值阈值时，则过滤掉该位置及其以后的所有数据；同时过滤掉相应记录的序列行在对应位置及其以后的所有数据；将标识行和“+”行原样输出。

S308、根据预定序列长度阈值对保留下来的高通量测序数据并行地进行过滤。

在经过预定质量值阈值过滤之后保留的数据中，同时获取第一文件和第二文件中一一对应的两个记录中的序列行中的碱基数，如果这一对记录中的任何一个记录的碱基数小于预定序列长度阈值时，则将这一对记录全部过滤掉。

将经过预定质量值阈值过滤和预定序列长度阈值过滤后保留下来的数据根据其原始数据来源分别输出到第一文件和第二文件各自对应的过滤结果文件中。

实施例四

如图4所示，本发明的第四个实施例提供了一种高通量测序数据的质量过滤装置，所述装置包括：

并行准备模块410，用于根据高通量测序数据为并行计算做准备；

并行计算模块420，用于通过并行计算过滤掉准备好的数据中质量不达标的数据。

在本实施例中，并行准备模块410根据高通量测序数据为并行计算做准备，并行计算模块420根据预定质量值阈值和预定序列长度阈值通过并行计算过滤掉准备好的数据中质量不达标的数据。

实施例五

如图4所示，本发明的第五个实施例提供了一种高通量测序数据的质量过滤装置，所述装置包括：

并行准备模块410，用于根据高通量测序数据为并行计算做准备；

质量过滤模块420，用于通过并行计算过滤掉准备好的数据中质量不达标的数据。

所述并行准备模块包括：

质量转换方式确定单元4101：用于根据所述高通量测序数据中的质量值确定质量转换方式；

数据切分单元4102：用于对已确定质量转换方式的高通量测序数据进行切分；

执行实体生成单元4103：用于生成对切分后的数据进行并行计算的执行实体。

在本实施例中，并行准备模块410中的质量转换方式确定单元4101根据所述高通量测序数据中的碱基质量值确定碱基质量值转换方式。并行准备模块410中的数据切分单元4102将已确定碱基质量值转换方式的高通量测序数据进行切分。并行准备模块410中的执行实体生成单元4103生成包含切分后的数据的执行实体map task和reduce task。质量过滤模块420通过并行计算过滤掉准备好的数据中质量不达标的数据。

实施例六

如图4所示，本发明的第六个实施例提供了一种高通量测序数据的质量过滤装置，所述装置包括：

并行准备模块410，用于根据高通量测序数据为并行计算做准备；

质量过滤模块420，具体用于：

如果所述高通量测序数据某一记录中的质量行中的某一位置上的质量值小于预定质量值阈值，则通过执行实体并行地过滤掉所述质量行该位置及以后的数据，以及同一记录中的序列行中的对应位置及以后的数据；

在保留下来的高通量测序数据中，如果第一文件和第二文件相对应的两个记录中有任意一个记录中的序列行长度小于预定长度阈值，则通过执行实体并行地过滤掉第一文件和第二文件中相对应的这两个记录。

所述并行准备模块包括：

质量转换方式确定单元4101：用于根据所述高通量测序数据中的质量值确定质量转换方式；

数据切分单元4102：具体用于：

将包含已确定质量转换方式的高通量测序数据的第一文件和第二文件分别转换为各自对应的第一RDD和第二RDD；

将第一RDD和第二RDD分别切分为各自对应的第一partition组和第二partition组；

根据第一文件和第二文件中对应的数据将第一RDD和第二RDD合并为第三RDD；

将第三RDD切分为第三partition组。

执行实体生成单元4103：生成对第三partition组进行并行计算的执行实体task。

在本实施例中，并行准备模块410中的质量转换方式确定单元4101根据所述高通量测序数据中的碱基质量值确定碱基质量值转换方式。并行准备模块410中的数据切分单元4102将包含已确定质量转换方式的高通量测序数据的第一文件和第二文件分别转换为各自对应的第一RDD和第二RDD,将第一RDD和第二RDD分别切分为各自对应的第一partition组和第二partition组,根据第一文件和第二文件中对应的数据将第一RDD和第二RDD合并为第三RDD,将第三RDD切分为第三partition组。并行准备模块410中的执行实体生成单元4103生成对第三partition组进行并行计算的执行实体task。如果所述高通量测序数据某一记录中的质量行中的某一位置上的质量值小于预定质量值阈值，则质量过滤模块420通过执行实体并行地过滤掉所述质量行该位置及以后的数据以及同一记录的序列行中的对应位置及以后的数据；在保留下来的高通量测序数据中，如果第一文件和第二文件相对应的两个记录中有任意一个记录中的序列行长度小于预定长度阈值，则质量过滤模块420通过执行实体并行地过滤掉第一文件和第二文件中相对应的这两个记录。

本领域普通技术人员可以理解，实现上述本发明实施例中的高通量测序数据统计方法和统计装置可以通过程序指令相关的硬件来完成，所述的程序可以存储于可读取存储介质中，该程序在执行时执行上述方法中的对应步骤。

以上所述仅是本发明的优选实施方式，应当指出，对于本技术领域的普通技术人员来说，在不脱离本发明原来的前提下，还可以做出若干改进和润饰，这些改进和润饰也应视为本发明的保护范围。