CN113270141B - 一种基因组拷贝数变异检测整合算法 - Google Patents

Abstract

本发明提出了一种基因组拷贝数变异检测整合算法，待测样本的测序序列进行数据筛选,保留高质量的测序序列,比对至相应参考基因组，使用固定长度的滑动窗口将比对后的高质量测序序列分成非重叠的片段，并计算每个窗口片段的原始Read数目的平均值作为该窗片段的深度信号；采用平均值校正法来纠正GC含量偏差，采用平滑分割算法将所有的经GC校正后的深度信号进行平滑，将相邻的窗口深度值一致的窗片段归并成大的片段；将一维空间中的平滑深度信号变换成二维平面，对二维平面建立高斯混合模型，采用步长搜索求解混合高斯模型的参数，对大于阈值概率的基因组片段的片段序列进行断点分析，计算相应序列处的拷贝数的增加或减少。

Description

一种基因组拷贝数变异检测整合算法

技术领域

本发明涉及生物信息学领域，具体涉及一种基因组拷贝数变异检测整合算法。

背景技术

近年来，高通量测序技术的迅猛发展给生命科学带来了巨大变革，测序通量不断提高，测序成本不断下降，如单个人类基因组的测序成本逐年降低，基因组数据规模不断增长。基因组数据大约每7个月就会增加一倍。现有服务器以及序列分析算法已经无法及时有效地处理如此大规模的数据，需要借助并行计算机的强大算力以及并行算法的有效支撑来实现对基因组大数据的有效处理。

测序仪产生Read后，首先进行质量控制，去除测序质量较差的数据；然后进行序列比对，找到Read最可能的位置，并输出比对文件；再基于比对文件检测基因组的变异情况，主要包括单核苷酸多态性(single nucleotide polymorphism，SNP)、插入删除变异(indel)、结构变异(structure variation，SV)和拷贝数变异(copy number variation，CNV)等；最后根据需要，进行特定的功能分析。序列比对和变异检测是基因组数据分析的基础环节，是后续功能性分析的基础，也是数据分析流程中最耗时的步骤。CNV是基因组结构变异的一种形式。CNV的狭义定义通常是指染色体中DNA片段的拷贝数变化。这种形式的基因组结构变异的类型和原因可以包括：缺失(末端缺失、间质缺失)；2.易位(相互易位、罗伯逊易位)；3.反转；4.环状染色体；5.双着丝粒染色体；6.CNV的更广泛的定义还包括例如结构变异，例如染色体非整倍性和部分非整倍性。目前可用的检测拷贝数变异的方法主要包括高分辨率染色体核型分析、FISH(荧光原位杂交)、阵列CGH(阵列比较基因组杂交)、MLPA(多重连接依赖性探针扩增)、PCR(聚合酶链反应)等，其中FISH检测被认为是遗传诊断的黄金标准，其可以有效地用于检测大多数已知的染色体缺失或重复。然而，这些方法通常具有低效率，特别是当用于全基因组扫描时，这可能消耗大量资源或可能无法检测未知CNV。

发明内容

为了解决上述技术问题，本发明提出了一种基因组拷贝数变异检测整合算法，包括如下步骤：

S1.对待测样本的测序序列进行数据筛选,去除含有接头序列的测序序列,保留高质量的测序序列，比对至相应参考基因组；

S2.使用固定长度的滑动窗口将比对后的高质量测序序列分成非重叠的片段，并计算每个片段的原始Read数目的平均值作为该片段的深度信号；

S3.采用平均值校正法来纠正GC含量偏差:

其中R_j表示第j个片段的原始Read深度信号值，R′_j表示第j个片段被校正后的深度信号值，R′_a表示基因组序列上所有片段的深度信号的平均值，R′_c表示所有与第j个片段具有相同GC含量的片段的深度信号的平均值；

S4.采用平滑分割算法将所有的GC含量偏差纠正后的深度信号进行平滑处理，将相邻的深度信号值一致的片段归并成大的片段；

S5.将一维空间中的平滑深度信号变换成二维平面，分别反映深度信号的幅度和位置空间；将片段的位置索引t_i视为第一维度，将幅度c_i视为第二维度，将一维深度信号C通过下式转换为二维平面C′，C′＝{(t_i，c_i)，1≤i≤m},其中m表示样本经平滑后的片段数目，c_i表明该样本平滑后第i个片段的深度信号值；

S6.对第一维度索引t_i进行校正，使其范围与第二维度的范围对齐，

S7.对二维平面C′建立如下高斯混合模型:

μ_k为均值矢量,

为协方差,π_k为参数，K的取值范围是2-6,第i个片段的二维Read深度信号(t_i，c_i)对应的拷贝数标签表示如下:

其中，

表示第k个高斯分量为生成(t_i，c_i)的权重；

完整的高斯混合模型参数

S8.采用步长搜索求解混合高斯模型的参数：先设定一个初始值λ₀，并以Δλ为固定步长动态的增长，进行模型调整，步长搜索结束后，最终选定参数λ’，分别计算计算每个基因组片段的高斯混合模型的概率值，将所述概率值和阈值概率进行比较，基因组片段的高斯混合模型的概率值小于阈值概率的，则该基因组片段发生拷贝数变异，否则则是正常基因组片段；

S9.对大于阈值概率的基因组片段的片段序列进行断点分析，鉴定位于断点区域之间的至少一个序列变异，基于至少一个序列变异计算相应碱基位置处的拷贝数的增加或减少，确定拷贝数变异的起始和终止位置。

进一步地，其特征在于：步骤S8中采用步长搜索求解混合高斯模型的参数的具体步骤为：首先假设高斯模型的参数是已知的，先设定一个参数初始值λ₀，利用参数初始值λ₀去估计每个高斯混合模型参数中的协方差

均值矢量μ_k、参数π_k和权重

以固定步长Δλ为1.05λ₀的倍数增长，进行固定步长的模型调整，选定最终的参数λ’后，将混合高斯模型固定为以λ’为参数的模型，分别计算每个基因组片段的高斯混合模型的概率值。

进一步地，步骤S9中，断点分析的具体步骤为：陆续扫描大于阈值概率的基因组片段的片段序列，记录每个窗口的断点bp1、……、bpf，多个断点将片段序列分成(f–1)个次级片段，计算每个次级片段窗口拷贝数的三均值，和拷贝数正常范围比较，三均值落在异常范围的次级片段为拷贝数变异区域，该拷贝数变异区域起始和终止的断点即为拷贝数变异的起始和终止位置。

进一步地，步骤S2中，平滑处理后的深度信号表示为:C＝{c₁，c₂…，c_m}，其中c_i表示平滑后的第i个片段的平均深度信号，m表示经平滑后的片段总数目。

进一步地，步骤S2中，固定长度为50bp、100bp、300bp、500bp、1000bp、1500bp。

进一步地，步骤S1中，设定碱基质量阈值，当一条测序序列有多个碱基质量值低于所述碱基质量阈值时，则过滤该测序序列，最后得到碱基质量均大于碱基质量阈值的高质量待测样本序列。

本发明还提出了一种计算机可读存储介质，所述计算机可读存储介质上存储有上述基因组拷贝数变异检测整合算法的步骤。

本发明与现有技术相比，具有如下优点与有益效果：

1、本发明通过待测样本的测序序列进行数据筛选,去除含有接头序列的测序序列,保留了高质量的测序序列,相比于现有技术具有更高的精确度和更好的稳定性。

2、本发明采用平滑分割算法将所有的经GC校正后的Read深度信号进行平滑，进一步地校正了偏差；通过采用平滑策略，可以将具有深度信号恒定相似性的局部相邻片段合并成一个大的片段，从而帮助减少一些系统噪声。

3、本发明考虑到基因组位置的重要性，将平滑的深度信号与其对应的基因组位置结合起来，将一维空间中的平滑深度信号变换成二维平面，分别反映拷贝数的幅度和位置空间。依据此二维平面，可从水平和垂直两个角度分析深度信号数据。

4、基于断点分析，鉴定位于断点区域之间的至少一个序列变异，基于至少一个序列变异计算相应碱基位置处的拷贝数的增加或减少。通过该断点分析过程可以完成CNV检测，这可以准确地检测包括拷贝数变异的区域，该拷贝数变异包括较小的微缺失/微重复。

本文使用的术语解释如下：

拷贝数变异(CNV)：是指与正常样品中的相应核酸序列相比，待测试样品中核酸分子的至少一部分的拷贝数变化，其中所述部分具有大于1kb的长度。拷贝数变异的情况和原因可包括：缺失，诸如微缺失；插入，例如微插入、微复制、复制；倒位、转座和复杂的多位点变异。

测序：是指获得样品的核酸序列的信息的过程。测序可以通过各种方法进行，包括但不限于双脱氧链终止；优选地，高通量测序方法，包括但不限于下一代测序技术或单分子测序技术。测序深度越高，检测的灵敏度越高，即可以检测的缺失片段和重复片段的长度越小。

Read：是指具有一定长度(通常长于20bp)的核酸序列，例如由测序仪产生的序列的测序结果，其可以通过序列比对方法与参考序列的特定区域或位置比对。

索引：指具有特定长度并发挥标记功能的核酸序列。当待测试的DNA分子衍生自多个待测试的样品时，多个样品中的每一个可以添加有不同的索引，用于在测序期间区分多个样品。

GC含量偏差：批次之间或一个批次内存在一定的GC偏差，这可能导致拷贝数偏差呈现在基因组的具有高GC含量或低GC含量的区域中。用基于对照集的测序数据进行CG校正，以获得每个窗口中的校正的相对读段数，由此可以消除这种偏差，并且可以提高检测拷贝数变异的准确性。

附图说明

附图1为本发明的一种基因组拷贝数变异检测整合算法的流程图。

附图2为示意性的高斯混合模型分布图。

附图3为4例样本的基因拷贝数变异断点分析结果。

附图4为示例性的应用本发明的基因组拷贝数变异检测整合算法检测血液中游离DNA的拷贝数变异后的显示图。

具体实施方式

下面通过具体实施方式对本发明作进一步详细说明。但本领域技术人员将会理解，下列实施例仅用于说明本发明，而不应视为限定本发明的范围。实施例中未注明具体技术或条件者，按照本领域内的文献所描述的技术或条件或者按照产品说明书进行。

本文使用的词语“包括”、“包含”、“具有”或其任何其他变体意欲涵盖非排它性的包括。例如，包括列出要素的工艺、方法、物品或设备不必受限于那些要素，而是可以包括其他没有明确列出或属于这种工艺、方法、物品或设备固有的要素。

针对上述现有技术中算法的局限性，本发明提出了一种基因组变异检测数据整合算法。

对样本基因组进行测序，以获得基因组序列,本发明对样本的类型不受特别限制，可以是含有大量核酸的样本，如植物的器官，动物的组织、血液、尿液、唾液、羊水，也可以是含有微量核酸的样本，如肿瘤的单细胞、血液、尿液、唾液中游离的单细胞、游离的核酸、生殖细胞、胚胎发育过程中的单细胞、单细胞或只有少量细胞的微生物。

对待测样本的基因序列进行数据筛选，去除含有接头序列的测序序列，保留高质量的测序序列，如果一条测序序列有多个碱基质量值低于20，则过滤该序列，得到碱基质量均大于20的高质量待测样本序列。碱基质量值为20时，认为该碱基的准确率为99％。根据需要可以调整该碱基质量阈值。例如，选择碱基质量值为30，即保留准确率大于99.9％的测序序列为高质量序列。

使用固定长度的滑动窗口将筛选后的高质量基因组序列分成非重叠的片段，固定长度可为50bp、100bp、300bp、500bp、1000bp、1500bp的长度。在一个示例中，300bp是优选的，500bp是更优选的。读段的长度可能由于不同的测序仪而具有大的差异。

计算每个片段的原始Read数目的平均值作为该片段的Read深度信号。对于GC含量偏差，采用平均值校正法来纠正GC含量偏差:

其中R_j表示第j个片段的原始Read深度信号值，R’_j表示第j个片段被校正后的深度信号值，R’_a表示基因组上所有片段的深度信号的平均值，R’_c表示所有与第j个片段具有相同GC含量的片段的深度信号的平均值。

为了进一步地校正偏差，采用平滑分割算法将所有的经GC校正后的深度信号进行平滑，从而将相邻的深度值一致的片段归并成大的片段。

平滑后的深度信号表示为:C＝{c₁，c₂…，c_m}，其中c_i表示平滑后的第i个片段的平均深度信号，m表示经平滑后的片段总数目。

通过采用平滑策略，可以将具有深度信号恒定相似性的局部相邻片段合并成一个大的片段，从而帮助减少一些系统噪声。

传统的基因变异检测算法通常依赖于相邻片段间深度值的差异，并对偏离整体水平的深度信号进行全局定位。虽然这种方法可以检测出基因变异数变化大的区域，但在其局部范围之外的发生小幅度基因变异区域却很难被检测到。

考虑到基因组位置的重要性，将平滑的深度信号与其对应的基因组位置结合起来，将一维空间中的平滑深度信号变换成二维平面，分别反映拷贝数的幅度和位置空间。

依据此二维平面，可从水平和垂直两个角度分析深度信号数据，鉴于深度信号C中基因组片段位置的连续依赖性，将片段的位置索引t_i视为第一维，将幅度c_i(即深度值)视为第二维。因此，将一维深度信号C通过下式转换为二维平面C’，

C’＝{(t_i，c_i)，1≤i≤m}，

其中m表示样本经平滑后的片段数目，c_i表明该样本平滑后第i个片段的读对深度信号值。

对第一维度索引t_i进行校正，使其范围与第二维度的范围对齐，如公式所示:

假定C′是由一个K个高斯分量模型混合生成的，且K的取值范围是2-6，对于一个基因片段，通过寻找生成该片段Read深度信号的最大概率值对应的高斯分量模型，来确定该片段的拷贝数状态。

对二维基因组片段的Read深度信号C′建立如下高斯混合模型:

μ_k为均值矢量,

为协方差，π_k为参数。那么第i个片段的二维Read深度信号(t_i，c_i)对应的拷贝数标签表示如下:

其中，

表示第k个高斯分量为生成(t_i，c_i)的权重，τ_i表示第i个片段的拷贝数标签。

完整的高斯混合模型参数由协方差

均值矢量μ_k和权重

组成，表示为:

由于Read深度信号有平滑的概率密度函数，所以有限数目的高斯密度函数足以对深度信号的密度函数形成一种平滑逼近。适当地选择高斯混合模型的权重，可以完成对一个概率密度函数的建模，减少需要估计的未知变量的数目。图2为示意性的高斯混合模型分布图。

本发明中用步长搜索算法对混合高斯模型的最优参数进行确定。采用步长搜索算法将求解分成两个步骤：首先假设高斯模型的参数是已知的，先设定一个初始值λ₀，利用参数去估计每个高斯混合模型参数中的协方差

均值矢量μ_k、参数π_k和权重

第二步，并以Δλ为固定步长动态的增长，进行模型调整，优选地，以固定步长Δλ为1.05λ₀的倍数增长，进行固定步长的模型调整，选择最优的参数λ′后，将混合高斯模型固定为以λ′为参数的模型，分别计算计算每个基因组片段的高斯分量模型的概率值。每个数据的概率值和阈值概率进行比较，小于阈值概率的认为是发生拷贝数变异，大于阈值概率的认为是正常区域。

对大于阈值概率的基因组片段的一组片段序列进行断点分析，由于拷贝数数目的变化反映在深度信号值的大小上，较大的深度信号值代表较多的拷贝数，较小的深度信号值代表较小的拷贝数，深度信号的走向反映拷贝数数目的变化，通过找出拷贝数数目有多变少或者由少变多的点，即断点。

当基于位于断点区域之间的碱基位置的深度计算的拷贝数可能性低于指定阈值时，可以执行断点分析，鉴定位于断点区域之间的至少一个序列变异，基于至少一个序列变异计算相应碱基位置处的拷贝数的增加或减少。

具体地断点分析过程为：陆续扫描大于阈值概率的基因组片段的片段序列。定义该组片段序列的第一个窗口为第1断点bp1，然后计算该组片段序列每个窗口及周围3个窗口的平均值。逐一计算每个窗口，当出现至少连续2个平均值落在异常范围时，记录该窗口为第2断点bp2，继续扫描，直到出现至少连续2个平均值回到正常范围时，记录该窗口为第3断点bp3，这样每遇到正常和异常转换的窗口，记录一个断点bpi，直到该片段序列的最后一个窗口，记录为bpf。断点bp1到断点bpf将组片段序列分成(f–1)个次级片段，计算每个次级片段窗口拷贝数的三均值，和拷贝数正常范围比较，三均值落在异常范围的次级片段即为精确的拷贝数变异区域，其中三均值为该次级片段的拷贝数，该次级片段起始和终止的断点即为拷贝数变异的起始和终止位置。通过该断点分析过程可以完成CNV检测，这可以准确地检测包括拷贝数变异的区域，该拷贝数变异包括较小的微缺失/微重复。

图3为采用的4例样本的基因拷贝数变异断点分析结果，对该4例检测结果进行了验证，结果显示通本发明分析的断点位置与验证结果平均误差<5bp，其中3个断点(共4个断点)位置完全相同。

图4为示例性的应用本发明的基因组拷贝数变异检测整合算法检测血液中游离DNA的拷贝数变异后的显示图。每个圆圈表示一个CNV片段，其尺寸与实际CNV片段的尺寸成比例。坐标平面旁边是对应第一维度索引ti和深度值维度的加权直方图，权重为每个CNV片段所含的杂合SNP位点的数量。

本发明另一方面将上述提出的基因组拷贝数变异检测整合算法存储在计算机可读存储介质，所述计算机可读存储介质上存储有机器可执行指令，所述机器可执行指令在被执行时使机器执行根据本发明提出的基因组变异检测数据整合算法的步骤。

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

用于执行本公开操作的机器可执行指令可以是汇编指令、指令集架构(ISA)指令、机器指令、机器相关指令、微代码、固件指令、状态设置数据、或者以一种或多种编程语言的任意组合编写的源代码或目标代码，该编程语言包括面向对象的编程语言—诸如Smalltalk、C++等，以及常规的过程式编程语言—诸如“C”语言或类似的编程语言。机器可执行指令可以完全地在用户计算机上执行、部分地在用户计算机上执行、作为一个独立的软件包执行、部分在用户计算机上部分在远程计算机上执行、或者完全在远程计算机或服务器上执行。

以上所述仅为本发明的较佳实施例而已，并不用以限制本发明，凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims (4)

1.一种基因组拷贝数变异检测整合算法，其特征在于，包括如下步骤：

S1.对待测样本的测序序列进行数据筛选,去除含有接头序列的测序序列, 保留高质量的测序序列，比对至相应参考基因组；

S2.使用固定长度的滑动窗口将比对后的高质量测序序列分成非重叠的片段，并计算每个片段的原始Read数目的平均值作为该片段的深度信号；所述固定长度为50bp、100bp、300bp、500bp、1000bp、1500bp；

S3.采用平均值校正法来纠正GC含量偏差：

其中

表示第j个片段的原始Read深度信号值，

表示第j个片段被校正后的深度信号值，

表示基因组序列上所有片段的深度信号的平均值，

表示所有与第j个片段具有相同GC含量的片段的深度信号的平均值；

S4. 采用平滑分割算法将所有的GC含量偏差纠正后的深度信号进行平滑处理，将相邻的深度信号值一致的片段归并成大的片段；

S5. 将一维空间中的平滑深度信号变换成二维平面，分别反映深度信号的幅度和位置空间；将片段的位置索引

视为第一维度，将幅度

视为第二维度，将一维深度信号C通过下式转换为二维平面

，

，其中m表示样本经平滑后的片段数目，

表明该样本平滑后第i个片段的深度信号值；

S6. 对第一维度索引

进行校正，使其范围与第二维度的范围对齐，

；

S7. 对二维平面

建立如下高斯混合模型：

；

为均值矢量,

为协方差,

为参数，K的取值范围是2-6,第i个片段的二维Read深度信号

对应的拷贝数标签表示如下：

；

其中，

表示第k个高斯分量为生成

的权重；

完整的高斯混合模型参数λ={

；

；

；

}，i=1，2，…，m；

S8. 采用步长搜索求解混合高斯模型的参数：首先假设高斯模型的参数是已知的，先设定一个参数初始值λ₀，利用参数初始值λ₀去估计每个高斯混合模型参数中的协方差

、均值矢量

、参数

和权重

；以固定步长Δλ为1.05λ₀的倍数增长，进行固定步长的模型调整，选定最终的参数

后，将混合高斯模型固定为以

为参数的模型，分别计算每个基因组片段的高斯混合模型的概率值，将所述概率值和阈值概率进行比较，基因组片段的高斯混合模型的概率值小于阈值概率的，则该基因组片段发生拷贝数变异，否则是正常基因组片段；

S9. 对大于阈值概率的基因组片段的片段序列进行断点分析，鉴定位于断点区域之间的至少一个序列变异，基于至少一个序列变异计算相应碱基位置处的拷贝数的增加或减少，确定拷贝数变异的起始和终止位置；

所述断点分析的具体步骤为：陆续扫描大于阈值概率的基因组片段的片段序列，定义该组片段序列的第一个窗口为第1断点bp1，然后计算该组片段序列每个窗口及周围3个窗口的平均值；逐一计算每个窗口，当出现至少连续2个平均值落在异常范围时，记录该窗口为第2断点bp2，继续扫描，直到出现至少连续2个平均值回到正常范围时，记录该窗口为第3断点bp3，每遇到正常和异常转换的窗口，记录该窗口为断点bpi，直到该组片段序列的最后一个窗口，记录为bpf；断点bp1到断点bpf将该组片段序列分成(f-1)个次级片段，计算每个次级片段窗口拷贝数的三均值，和拷贝数正常范围比较，三均值落在异常范围的次级片段为精确的拷贝数变异区域，其中三均值为该次级片段的拷贝数，该次级片段起始和终止的断点为拷贝数变异的起始和终止位置。

2.根据权利要求1所述的基因组拷贝数变异检测整合算法，其特征在于：步骤S2中，平滑处理后的深度信号表示为：C={c₁，c₂…，c_m}，其中，m表示经平滑后的片段总数目。

3.根据权利要求1所述的基因组拷贝数变异检测整合算法，其特征在于：步骤S1中，设定碱基质量阈值，当一条测序序列有多个碱基质量值低于所述碱基质量阈值时，则过滤该测序序列，最后得到碱基质量均大于碱基质量阈值的高质量待测样本序列。

4.一种计算机可读存储介质，所述计算机可读存储介质上存储有如权利要求1-3中任意一项所述的基因组拷贝数变异检测整合算法的步骤。

Images