CN110021342B - 用于加速变异位点的识别的方法及系统 - Google Patents

Abstract

本发明提供一种用于加速变异位点的识别的方法，其将参考基因组的各染色体划分为一个或多个区间，基于各区间的测序片段密度来设置各区间的处理优先级，并根据各区间的处理优先级来处理各区间中的测序片段，获取各区间中与每个位点的堆叠相关的统计数据，从而判断各位点是否发生变异。该方法将参考基因组划分成多个区间，可以在多个核上同时并行处理各区间的相关数据，能有效利用现有的多核计算资源并缩短计算时间。

Description

用于加速变异位点的识别的方法及系统

技术领域

本发明涉及基因数据处理，尤其涉及对变异位点的识别进行加速的方法和系统。

背景技术

随着下一代测序技术(NGS)技术的快速发展，基因测序通量不断提升。在DNA测序中，DNA分子首先会被随机打断成很多个片段，然后将这些无序片段克隆，再经由测序仪进行测序，产生数百万很短的DNA reads(下文中也可称为测序片段)。之后将这些产生的DNAreads与参考基因组进行比对，找到这些DNA reads在参考基因组上的位置，接着识别变异位点SNP(Single Nucleotide Polymorphisms)，以发现基因组上单个核苷酸的变异。

现有NGS中识别变异位点的方法通常包括遍历参考序列上的各个位点，获得覆盖各位点位置的reads堆叠(pileup)，接着基于各位点位置的堆叠计算基因型概率等信息，然后基于这些信息进行统计分析来识别变异位点。

发明内容

但是，现有的识别变异位点的方法流程主要是单线程实现，要求顺序遍历基因参考序列上的每个位点，计算时间长而且对进行计算资源要求很高，不能有效地利用现有的多核计算资源。

因此，本发明的目的在于克服上述现有技术的缺陷，提供一种用于加速变异位点的识别的方法及系统。

本发明的目的是通过以下技术方案实现的：

一方面，本发明提供了一种用于加速变异位点的识别的方法，包括：

将参考基因组的各染色体划分为一个或多个区间；

基于各区间的测序片段密度来设置各区间的处理优先级；

由多个处理核根据各区间的处理优先级来处理各区间中的测序片段，以获取各区间中与每个位点的堆叠相关的统计数据；

基于与每个位点的堆叠相关的统计数据来判断该位点是否发生变异。

在上述方法中，所述区间可以是按照所述处理核的个数来划分的。

在上述方法中，所述区间的大小可以为所述参考基因组的各染色体中的最小长度。

在上述方法中，所述每个区间的测序片段密度可以为该区间中测序片段的个数与该区间的大小的比值。

在上述方法中，其中相邻两区间之间有重叠。

在上述方法中，所述重叠的长度至少为测序片段的长度。

在上述方法中，对于跨越两个区间的测序片段，可允许在这两个区间重复出现。

又一方面，本发明提供了一种用于加速变异位点的识别的系统，包括：

用于将参考基因组的各染色体划分为一个或多个区间的装置；

用于基于各区间的测序片段密度来设置各区间的处理优先级的装置；

用于由多个处理核根据各区间的处理优先级来处理各区间中的测序片段，以获取各区间中与每个位点的堆叠相关的统计数据的装置；

用于基于与每个位点的堆叠相关的统计数据来判断该位点是否发生变异的装置。

在上述系统中，所述每个区间的测序片段密度为该区间中测序片段的个数除以该区间的大小。

在上述系统中，相邻两区间之间有重叠。

与现有技术相比，本发明的优点在于：

将参考基因组的各个染色体划分成多个区间，基于各区间的reads密度安排各区间的处理顺序，在多个核上并行处理各区间的相关数据，能有效利用现有的多核计算资源并缩短计算时间。

附图说明

以下参照附图对本发明实施例作进一步说明，其中：

图1给出了位点的堆叠示例示意图；

图2为现有的识别变异位点的方法的流程示意图；

图3为根据本发明实施例的用于加速变异位点的识别的方法的流程示意图；

图4为根据本发明实施例的参考基因划分示意图；

图5为利用本发明实施例的方法的变异位点的识别流程示意图。

具体实施方式

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

在NGS技术中，通常采用samtools和bcftools这两个工具来进行变异位点的识别，下面也以这两个工具为例来说明变异位点的识别过程。其中主要包括两个关键的部分：samtools mpileup和bcftools call，其中samtools mpileup主要负责计算每个位点的基因型概率(genotype likelihood)等辅助信息，而bcftools call根据这些信息并基于统计模型的方法来判定该位点是否为SNP。在变异位点识别中很关键的一个概念是堆叠(pileup)，变异位点的判定的依据是基于覆盖该位点的堆叠得到的统计数据。图1给出了一个位点的堆叠示例示意图。堆叠指的是多个reads的堆叠，如图1所示，标号为i、j、k、l、m、n、l的6条read都覆盖了位点s，那么这六条read构成了位点s的堆叠，该位点s的深度为覆盖该位点s的reads的总数，即该位点s的深度为6。图2给出了使用samtools mpileup和bcftoolscall进行变异位点识别的流程示意。如图2所示，首先samtools mpileup顺序遍历参考序列上的每个位点，获取该位点的堆叠，基于该堆叠计算基因型概率等信息然后以vcf文件格式写入管道中，之后bcftools call从管道中读取每个位点的信息，先通过Is ref过程判断其是否和参考序列上该位点相同，如果是则直接丢弃掉，否则再进行调用变异检测(callvcf)过程发现变异位点。从图2可以看出，该识别变异位点的过程是单线程实现，计算时间长而且对进行计算资源要求很高，因此硬件成本也较高，无法有效利用现有的多核计算资源。并且利用samtools mpileup和bcftools call在进行变异位点发现过程中，参考序列上位点大概有30亿，即使每个位点10字节数据，大概有30G(30*10^8*10/10^9＝30G)的数据，管道读写数据量大，非常耗时。

针对目前samtools mpileup和bcftools识别变异位点方法中的问题，在本发明的实施例中提出了利用多核装置加速变异位点的识别的方法。图3给出了根据本发明一个实施例的用于加速变异位点的识别的方法的流程示意图。该方法主要包括下列步骤：

在步骤S1)，将参考基因组的各个染色体划分成一个或多个区间，这样每个区间对应的reads数据可以在不同的核上并行处理。区间的数量和/或大小可以根据实际情况来进行设定。例如，可以根据执行变异位点识别方法的计算装置的处理核的个数来划分，例如，如果多核计算装置的处理器有4个核，则可以将每个染色体划分为4个区间。在这种情况下，由于不同染色体长度不同，因此不同染色体对应的区间的大小也不同。又例如，考虑到各个染色体长度差异较大，可以根据最小染色体的长度设置区间大小，这样所有的染色体都按照相同大小的区间划分。在划分区间之后，可以在多核装置的多个处理核上同时处理不同区间内的reads数据。

通常可以给在步骤S1)划分的各个区间设置相同的优先级，每个核随机地调度待处理的区间数据进行处理。但实际上各个read在参考基因组上的位置分布式不均匀的，这会使得有些区间对应的reads很多，其每个位点的深度过大，需要很长的处理时间，而有些区间对应的reads较少，只需要很少处理时间，这样可能会引起多核装置中各个核上处理任务的不均衡，出现长尾任务现象，从而拖延整个变异位点识别流程的处理时间。

因此，优选地，在步骤S2)可以基于各区间的reads密度来设置各区间的处理优先级。每个区间的read密度为该区间中read的个数与该区间的大小的比值，即read密度＝区间内read数目/区间大小。read密度越大，处理优先级越高。也就是优先针对read密度大的区间进行处理，这样能有效地减少长尾的发生，缩短流程的处理时间。

在优选的实施例中，为了改善结果的准确性，允许通过添加重叠的方式来统计各区间的reads。如图4所示，在将参考基因组的各个染色体划分成一个或多个区间时，可能会导致某些reads横跨两个区间，为了保持结果更精确，可以对这些区间两端进行适度延长，即和相邻的区间有重叠部分，也就是各区间可以分别向前后扩展所设置的重叠(overlap)的长度。这样设置的重叠的长度至少等于或大于测序片段的长度。对于跨越两个区间的测序片段，允许在这两个区间重复出现或重复统计。

继续参考图3，在步骤S3)由多核装置的多个核按照各区间的处理优先级来调度和处理各个区间的数据，不同的区间数据可以在不同的核上同时进行处理。可以对于每个区间，依次遍历该区间参考序列上的各个位点，获得覆盖各位点位置的reads堆叠(pileup)，接着计算与各位点的堆叠相关的统计数据例如，可以对于每个区间，依次遍历该区间参考序列上的各个位点，获得覆盖各位点位置的reads堆叠(pileup)，接着计算与各位点的堆叠相关的统计数据，例如基因型概率等信息。

在步骤S4)基于这些统计数据识别变异位点。图5给出了利用本发明实施例的方法的变异位点的识别流程的一个示例。如图5所示，对于每个区间，依次遍历该区间参考序列上的各个位点，获得覆盖各位点位置的reads堆叠(pileup)，接着计算与各位点的堆叠相关的统计数据，例如基因型概率等信息，并将这些统计数据保存在vcf文件中。多个核的并行处理的结果会得到多个vcf文件，但这些vcf是无序的，可按照参考基因组上区间的顺序对这些vcf文件排序并合并成一个文件。然后利用bcftools的变异检测(call vcf)函数处理合并后的vcf文件来识别其中的变异位点。在另外的实施例中，也可以不进行各个核产生的vcf文件的合并，而是由每个核在处理各区间时直接给出识别结果。即针对每个区间来执行步骤S3和S4，直接识别出该区间的变异位点。这样，各区间的变异位点识别是在多个核上并行处理的，大大缩短了处理时间。

虽然本发明已经通过优选实施例进行了描述，然而本发明并非局限于这里所描述的实施例，在不脱离本发明范围的情况下还包括所做出的各种改变以及变化。

Claims (10)

1.一种用于加速变异位点的识别的方法，该方法包括：

将参考基因组的各染色体划分为一个或多个区间；

基于各区间的测序片段密度来设置各区间的处理优先级；

由多个处理核根据各区间的处理优先级来处理各区间中的测序片段，以获取各区间中与每个位点的堆叠相关的统计数据；

基于与每个位点的堆叠相关的统计数据来判断该位点是否发生变异。

2.根据权利要求1所述的方法，其中所述区间是按照所述处理核的个数来划分的。

3.根据权利要求1所述的方法，其中所述区间的大小为所述参考基因组的各染色体中的最小长度。

4.根据权利要求1所述的方法，其中每个所述区间的测序片段密度为该区间中测序片段的个数与该区间的大小的比值。

5.根据权利要求1所述的方法，其中相邻两区间之间有重叠。

6.根据权利要求5所述的方法，其中，所述相邻两区间之间的重叠的长度至少为测序片段的长度。

7.根据权利要求1所述的方法，其中对于跨越两个区间的测序片段，允许在这两个区间重复出现。

8.一种用于加速变异位点的识别的系统，该系统包括：

用于将参考基因组的各染色体划分为一个或多个区间的装置；

用于基于各区间的测序片段密度来设置各区间的处理优先级的装置；

用于由多个处理核根据各区间的处理优先级来处理各区间中的测序片段，以获取各区间中与每个位点的堆叠相关的统计数据的装置；

用于基于与每个位点的堆叠相关的统计数据来判断该位点是否发生变异的装置。

9.根据权利要求8所述的系统，其中每个所述区间的测序片段密度为该区间中测序片段的个数与该区间的大小的比值。

10.根据权利要求8所述的系统，其中相邻两区间之间有重叠。

Images