CN108197428B - 一种并行动态规划的下一代测序技术拷贝数变异检测方法 - Google Patents

Abstract

本发明公开了一种并行动态规划的下一代测序技术拷贝数变异检测方法，属于基因测序技术领域。本发明基于图形处理器并行动态规划的分割算法，使用该算法处理下一代测序NGS技术获得的数据，能够快速、准确地检测出拷贝数变异CNV，并取得碱基水平的分辨率。本发明方法可以有效应用于生殖健康诊断、孕妇的产前筛查、新生儿遗传病的基因诊断等临床应用，以及考古、生物、医学等科学研究，预防及减少患儿的出生，对于提高国民身体素质、推进科学研究工作具有重大意义。

Description

一种并行动态规划的下一代测序技术拷贝数变异检测方法

技术领域

本发明属于基因测序技术领域，具体涉及一种并行动态规划的下一代测序技术拷贝数变异检测方法。

背景技术

现有的拷贝数变异(CNV，Copy Number Variation)技术有荧光原位杂交(FISH，Fluorescent in situ Hybridization)和基于微阵列的比较基因组杂交技术(CGH，Comparative Genomic Hybridization)。

FISH是上世纪80年代由Langer-Safer等人开发，以检测和定位染色体上特定DNA区段的存在或者缺失，即CNV的有无。FISH使用与染色体的目标区域高度互补的荧光探针，因此可以结合到感兴趣的区域。在荧光显微镜观测下，可以找出荧光探针结合的目标区域。而且，当使用多种颜色的荧光探针时，可以同时染色更多的目标区域。通过分析颜色通道的不同组合，可以研究多个目标区域。这种技术被称为多色FISH，常被用来检测染色体之间区段重排引起的变异及CNV。但是由于它的分辨率受限于使用显微镜的性能，FISH检测的CNV分辨率通常较低，一般只能达到显微水平的5～10Mbp。

CGH最初是为了研究肿瘤而开发的另一种细胞遗传学技术。该技术对测试样品和对照样品的DNA进行差异化标记，并通过竞争杂交的方式与参考DNA进行结合。如果测试样品具有不平衡的结构变异如CNV，竞争杂交将是有差异的，导致这两种荧光探针密度之间的比率会偏离1，进而可以检测到拷贝数的增加或者缺失。结合微阵列技术开发的阵列CGH可以提升变异检测的分辨率，达到10～25kbp。以Affymetrix公司的SNP6.0阵列为例，该芯片可以检测人类全基因组中的186万个位点(91万个SNP位点以及95万个CNV位点)，等效的分辨率约为15kbp。

因此，应用上述技术检测CNV，主要的缺点是通量小、速度慢、需要设计探针、对小片段不灵敏等问题。

发明内容

为了克服上述现有技术存在的缺陷，本发明的目的在于提供一种并行动态规划的下一代测序(Next Generation Sequencing，NGS)技术拷贝数变异检测方法，使用该检测方法处理NGS技术获得的数据，能够快速、准确地检测出CNV，并取得碱基水平的分辨率。

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

本发明公开了一种并行动态规划的下一代测序技术拷贝数变异检测方法，包括以下步骤：

1)使用对齐映射工具将下一代测序平台得到的短序列数据对齐映射到参考基因组；

2)计算RD信号，RD信号是基因组上固定宽度窗口或滑动窗口内短序列的个数或者每个碱基位点的覆盖深度；

3)采用基于图形处理器并行动态规划的分割算法将RD信号分割成不同深度的区域；

4)通过统计假设检验来确定每个分割出的段的拷贝数，并合并拷贝数相同的相邻段；

5)检测拷贝数变异并输出记录。

优选地，步骤2)和步骤3)之间还包括对RD信号进行归一化处理以矫正偏差影响的操作。

进一步优选地，由于基因组上的不均一特性，导致RD会存在偏差，因此需要对RD进行归一化以矫正偏差的影响。常用的归一化如鸟嘌呤胞嘧啶(Guanine-Cytosine，GC)含量矫正，映射矫正等。对于单一样品的CNV检测，此步骤是必须的；但是对于病例-对照研究的双样品的CNV检测，此步骤不是必须的。

优选地，由于NGS技术的的特点是霰弹枪测序，因此随机打碎的片段重新对齐映射后所得到的RD正比于局部的拷贝数。大多数的基因组位点都处于正常区域，在RD上的一个高起的平台反映了复制增益状态，而一个低陷的盆地则反映了复制损失的状态。

步骤3)中，采用基于图形处理器并行动态规划的分割算法进行RD信号分割的操作，对于低于3倍覆盖度，首先应用式(1)计算：

其中，x表示RD信号，

和e_i:j分别表示RD信号从位置i到j的平均值和负对数似然函数误差；

对于高于6倍覆盖度，应用式(2)计算：

其中，

表示RD信号从位置i到j的方差；

对于在3倍～6倍覆盖度之间，则应用式(1)或式(2)均可；

得到e_i:j之后，使用如下式所示的并行动态规划的算法对RD进行分割：

其中，M为RD信号的长度，函数f_k(j)表示从位置1到j分割为k段的最优目标函数，常数c表示每个分段的惩罚系数，j_k表示第k个分段位置，g_k(j)是一个后向指针，存储优化过程中的临时分段位置，K是最大分段数；

最终分割的段与段之间的断点位置存储在j_k里；

具体计算时：

1)由于f_k(j)依赖于f_k-1(j)，但不依赖于其他的f_k(i),i≠j，使用图形处理器将f_k(j)计算并行化，使用图形处理器计算式(3)，即同时计算多个f_k(j)；

2)由于计算时间决定于K值的大小，K越大计算时间越长，则通过下式(5)计算出K的最大值：

其中，e₀是所有e_i:j的最小值；

3)为限制分割出的段的最小长度，即所有分割的段应大于长度L，则：

即将位置i和j小于距离L的负对数似然函数误差e置为正无穷。

优选地，步骤1)中，使用对齐映射工具将NGS测序平台得到的短序列数据对齐映射到参考基因组(如NCBI37/hg19)，例如bowtie，MAQ等。对齐映射工具将通常输出小片段的变异作为副产品；将映射索引、映射质量及原始测序数据信息存储在SAM格式的文件或者压缩的二进制BAM格式文件后输出，便于后续大片段变异的检测。

优选地，步骤2)中，通过对SAM格式或BMA格式文件的分析计算出RD信号。

优选地，步骤4)中，通过统计假设检验来分割出的每段的拷贝数。例如根据假设RD服从泊松分布或者负二项分布来估计每段的拷贝数；合并拷贝数相同的连续的段是因为由于步骤3)可能导致过度的分割，造成过多的小片段，因此需要通过这一步来拼接上被错误分割的段。

使用公式(5)的方法也能够有效抑制过度分割的问题。

优选地，步骤5)中，记录的拷贝数变异信息包括拷贝数增益或缺失、起始位点及长度以及拷贝数状态信息。

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

本发明公开的并行动态规划的下一代测序技术拷贝数变异检测方法，采用图形处理器并行动态规划的分割算法，使用该算法处理下一代测序NGS技术获得的数据，能够快速、准确地检测出拷贝数变异CNV，并取得碱基水平的分辨率。本发明方法可以有效应用于生殖健康诊断、孕妇的产前筛查、新生儿遗传病的基因诊断等临床应用，以及考古、生物、医学等科学研究，预防及减少患儿的出生，对于提高国民身体素质、推进科学研究工作具有重大意义。

附图说明

图1为本发明方法与循环二分分割(Circular Binary Segmentation，CBS)方法的分割结果比较(仿真数据)；

图2为本发明方法与CBS方法的检测结果比较(真实数据)；

图3为本发明方法对计算时间的改善。

具体实施方式

下面结合具体的实施例对本发明做进一步的详细说明，所述是对本发明的解释而不是限定。

本发明出现和使用的缩略语和关键术语定义如下：

CNV Copy Number Variation 拷贝数变异

SNP Single Nucleotide Polymorphism 单核苷酸多态性

NGS Next Generation Sequencing 下一代测序

PC Parallel Computing 并行计算

DP Dynamic Programming 动态规划

CPU Central Processing Unit 中央处理器

GPU Graphics Processing Unit 图形处理器

FISH Fluorescent in situ Hybridization 荧光原位杂交

CGH Comparative Genomic Hybridization 比较基因组杂交技术

RD Read Depth 读深

GC Guanine-cytosine content 鸟嘌呤胞嘧啶

CBS Circular Binary Segmentation 循环二分分割

本发明公开的并行动态规划的下一代测序技术拷贝数变异检测方法的处理流程主要包括以下七个步骤：

1)使用对齐映射工具将NGS测序平台得到的短序列数据对齐映射到参考基因组(如NCBI37/hg19)，例如bowtie，MAQ等。这些对齐映射工具通常输出SNP和indel等小片段的变异作为副产品。映射索引、映射质量、原始测序数据等信息存储在SAM格式的文件，或者压缩的二进制BAM格式文件作为输出。

2)计算RD信号。

RD是基因组上固定宽度窗口或滑动窗口内短序列的个数或者每个碱基位点的覆盖深度，通过SAM/BMA文件的分析可以计算出RD信号。

3)对RD信号进行归一化。

由于基因组上的不均一特性，导致RD会存在偏差，因此需要对RD进行归一化以矫正偏差的影响。常用的归一化如GC含量矫正，映射度矫正等。

因此，对于单一样品的CNV检测，此步骤是必须的；但是对于病例-对照研究的双样品的CNV检测，因为这些偏差会相互抵消，所以此步骤不是必须的。

4)将RD信号根据覆盖深度进行分割。

由于NGS技术的特点是霰弹枪测序，因此随机打碎的短片段重新对齐映射后所得到的RD信号正比于局部的拷贝数。大多数的基因组位点都处于正常区域，在RD上的一个高起的平台反映了拷贝数增益状态，而一个低陷的盆地则反映了拷贝数缺失的状态。

5)通过统计假设检验来确定每个分割出的段的拷贝数。

例如可以根据假设RD信号服从泊松分布或者负二项分布来估计每段的拷贝数。

6)合并拷贝数相同的相邻段。由于步骤4)可能导致过度的分割，造成过多的小片段，因此需要通过这一步来拼接上被错误分割的段。

7)检测CNV并输出记录。记录的信息包括CNV的类型(增益或缺失)，起始位点及长度，拷贝数状态等信息。

优选地，步骤4)如何精确的分割RD信号是核心问题。针对这个问题，本发明公开了一种基于GPU并行计算的DP算法。该方法的技术关键点是：

其中x表示步骤4)中的RD信号，长度为M，

和e_i:j分别表示RD信号从位置i到j的平均值、方差和负对数似然函数误差，函数f_k(j)表示RD信号从位置1到j分割为k段的最优目标函数，常数c表示每个分段的惩罚系数，j_k表示第k个分段位置，g_k(j)是一个后向指针，存储优化过程中的临时分段位置，K是最大分段数。

对于低覆盖度(低于3倍覆盖)应使用(1)式，高覆盖度(高于6倍覆盖)应使用(2)式，之间的覆盖度取两者皆可。

其中核心技术为：

1)由于f_k(j)依赖于f_k-1(j)，但不依赖于其他的f_k(i),i≠j，因此我们使用GPU将f_k(j)的计算并行化，使用GPU计算公式(3)，即同时计算多个f_k(j)。

2)由于计算时间决定于K值的大小，K越大计算时间越长，因此可以通过下边公式计算出K的最大值，从而避免不必要的计算，进而提升计算速度。

其中e₀是所有e_i:j的最小值。

3)为了限制分割出的段的最小长度，即所有分割的段应大于长度L，可以通过下面的方法实现：

即将位置i和j小于距离L的负对数似然函数误差e置为正无穷。

采用本发明方法与CBS方法进行比较：

参见图1，展示了采用本发明的方法与现有技术CBS方法的分割结果。这里使用了仿真的数据。可以看到，在位置400处RD有一段幅度约为60，宽度为10的脉冲信号，对应的是一个CNV增益。本发明可以找到这个CNV，而作为对照的CBS则无法找到该CNV。

参见图2，展示了采用本发明的方法与CBS方法的检测结果。这里使用了真实的NGS数据，该数据来自于一名52岁的患有乳腺癌的白人女性。可以看到，在17号染色体的59.8到59.9Mbp位点有一个CNV缺失，对应的是乳腺癌相互作用蛋白1(Breast CancerInteracting Protein 1，BRIP1)。通过本发明方法可以找到这个CNV，而作为对照的CBS方法则无法找到该CNV。

图3展示了采用本发明方法对计算时间的改善结果。可以看到，对于相同维度的问题，相比于传统的CPU计算，本发明方法的计算速度提高了约10倍。

Claims (5)

1.一种并行动态规划的下一代测序技术拷贝数变异检测方法，其特征在于，包括以下步骤：

1)使用对齐映射工具将下一代测序平台得到的短序列数据对齐映射到参考基因组；

2)计算RD信号，RD信号是基因组上固定宽度窗口或滑动窗口内短序列的个数或者每个碱基位点的覆盖深度；

3)采用基于图形处理器并行动态规划的分割算法将RD信号分割成不同深度的区域，具体包括如下操作，对于低于3倍覆盖度，首先应用式(1)计算：

其中，x表示RD信号，

和e_i:j分别表示RD信号从位置i到j的平均值和负对数似然函数误差；

对于高于6倍覆盖度，应用式(2)计算：

其中，

表示RD信号从位置i到j的方差；

对于在3倍～6倍覆盖度之间，则应用式(1)或式(2)均可；

得到e_i:j之后，使用如下式所示的并行动态规划的算法对RD进行分割：

其中，M为RD信号的长度，函数f_k(j)表示从位置1到j分割为k段的最优目标函数，常数c表示每个分段的惩罚系数，j_k表示第k个分段位置，g_k(j)是一个后向指针，存储优化过程中的临时分段位置，K是最大分段数；

最终分割的段与段之间的断点位置存储在j_k里；

具体计算时：

1)由于f_k(j)依赖于f_k-1(j)，但不依赖于其他的f_k(i),i≠j，使用图形处理器将f_k(j)计算并行化，使用图形处理器计算式(3)，即同时计算多个f_k(j)；

2)由于计算时间决定于K值的大小，K越大计算时间越长，则通过下式(5)计算出K的最大值：

其中，e₀是所有e_i:j的最小值；

3)为限制分割出的段的最小长度，即所有分割的段应大于长度L，则：

即将位置i和j小于距离L的负对数似然函数误差e置为正无穷；

4)通过统计假设检验来确定每个分割出的段的拷贝数，并合并拷贝数相同的相邻段；

5)检测拷贝数变异并输出记录。

2.根据权利要求1所述的并行动态规划的下一代测序技术拷贝数变异检测方法，其特征在于，步骤2)和步骤3)之间还包括对RD信号进行归一化处理以矫正偏差影响的操作。

3.根据权利要求1所述的并行动态规划的下一代测序技术拷贝数变异检测方法，其特征在于，步骤1)中，对齐映射工具将通常输出小片段的变异作为副产品；将映射索引、映射质量及原始测序数据信息存储在SAM格式的文件或者压缩的二进制BAM格式文件后输出，便于后续大片段变异的检测。

4.根据权利要求3所述的并行动态规划的下一代测序技术拷贝数变异检测方法，其特征在于，步骤2)中，通过对SAM格式或BMA格式文件的分析计算出RD信号。

5.根据权利要求1所述的并行动态规划的下一代测序技术拷贝数变异检测方法，其特征在于，步骤5)中，记录的拷贝数变异信息包括拷贝数增益或缺失、起始位点及长度以及拷贝数状态信息。

Images