基于信号谱差异的混合样本单核苷酸多态性的检测方法与流程

文档序号:12167854阅读:267来源:国知局
基于信号谱差异的混合样本单核苷酸多态性的检测方法与流程

本发明涉及一种单核苷酸多态性的检测方法,尤其涉及一种基于信号谱差异的混合样本单核苷酸多态性的检测方法。



背景技术:

单核苷酸多态性(Single nucleotide polymorphism,SNP)是群体基因组中单个核苷酸位置上发生的突变。目前,SNP是使用最广泛的分子标记之一。以人类基因组为例,单核苷酸多态性占遗传突变的绝大部分,而且可能是决定个体差异的主要原因。据统计,人类基因组中大约每1000-2000个核苷酸中就有一个SNP,SNP的高丰度、易于检测等特点决定了它非常适用于对复杂性状与疾病的遗传研究以及基于群体的基因识别的研究等方面。

SNP已经成为鉴定遗传因子的一种重要工具并在许多领域得到了实际应用,例如法医鉴定、临川诊断、动植物育种以及改善生物分子产量等。随着高通量测序技术的发展,对SNP进行大规模鉴定分析在许多物种中已经变得可行。然而,从高通量测序数据中检测SNP通常会受到测序错误以及测序片段比对错误的影响。因此,在实际应用之前,应该对发现的SNP进行进一步的验证以筛选出真正的SNP。

现有技术中有一种基于单核苷酸添加焦磷酸测序技术,从混合样本中检测SNP并计算其比例的方法。该方法基于正态性检验和动态规划算法,输入纯野生型样本和混合样本的焦磷酸测序信号谱,能够检测出SNP的位置、比例以及突变核苷酸。该方法的一个缺点是,在焦磷酸测序过程中,核苷酸的添加过程需要经过严格设计以保证野生型序列与突变型序列的延伸节奏不一致。



技术实现要素:

发明目的:针对以上问题,本发明提出一种基于信号谱差异的混合样本单核苷酸多态性的检测方法。

技术方案:为实现本发明的目的,本发明所采用的技术方案是:一种基于信号谱差异的混合样本单核苷酸多态性的检测方法,包括以下步骤:

第一步:提取野生型样本的DNA和混合样本的DNA,进行聚合酶链式反应,得到单链PCR产物,其中,混合样本包括野生型序列和突变型序列;

第二步:对野生型样本和混合样本的单链PCR产物进行三轮循环的实时合成测序,获得三组信号谱;

第三步:根据三轮测序实验获得的野生型样本与混合样本的信号谱,计算野生型样本与混合样本的信号谱之间的差异,利用枚举方法推断混合样本中可能存在的单核苷酸多态性位点;

第四步:综合分析三轮测序实验所检测出的可能的单核苷酸多态性位点,进行一致性判断,判断混合样本中是否存在单核苷酸多态性位点,如果混合样本中含有单核苷酸多态性位点,则归纳出单核苷酸多态性位点的位置、突变核苷酸及突变核苷酸的比例。

第三步中利用枚举方法推断混合样本中可能存在的单核苷酸多态性位点的步骤如下:

(1)定义每轮测序实验获得野生型样本的信号谱为W=(W1,W2,…,Wn),其中,Wi为野生型样本第i次添加两核苷酸时的信号强度;混合样本的测序信号谱为P=(P1,P2,…,Pn),其中,Pi为混合样本第i次添加两核苷酸时的信号强度;

(2)定义突变型序列的信号谱为M=(M1,M2,…,Mn),其中,突变序列的比例为a,则P=(1-a)W+aM;

(3)实际测序试验中,Wi、Pi服从正态分布:其中,CV为信号强度的变异系数,和为信号强度的理论值;

(4)定义差别谱D为混合样本的信号谱P与野生型样本的信号谱W之差,即D=P-W,则结合步骤(2)可得D=a(M-W),枚举野生型序列与突变型序列信号谱之间的差别,计算混合样本中突变序列的比例a;

(5)结合步骤(3)和(4)可得定义Sp来评估所枚举的野生型序列与突变型序列信号谱之间每种差别的可能性,且选择最大Sp值时的差别谱作为野生型序列与突变型序列信号谱之间的差别谱,其中:

(6)野生型序列与突变型序列信号谱之间存在五种差别类型:Type-0、Type-1、Type-2、Type-3和Type-4,假设野生型序列与突变型序列之间的不一致延伸发生在第j步循环测序反应过程,定义S为第j步循环测序反应中一致延伸的核苷酸个数,定义函数f(Wj)为根据信号强度Wj计算的延伸核苷酸的个数,可以得出S的表达式和单核苷酸多态性位点的位置L的表达式分别为:

其中,为突变序列的比例的平均值;

(7)结合步骤(6)检测出的单核苷酸多态性位点的位置L和加入的两核苷酸信息得到突变核苷酸的信息。

第四步中综合分析进行一致性判断的步骤如下:

(1)定义Sc来评价野生型样本的信号谱与混合样本的信号谱之间的一致性,同时定义Sm=Sc-Sp;当突变序列存在时,Sc的值大于Sp,突变序列的比例越高,Sm的值越大;其中:

(2)根据三轮测序实验结果,定义具有最低Sm值的一轮测序实验中的野生型序列与突变型序列的信号谱完全一致,即野生型序列与突变型序列信号谱之间的差别类型为Type-0;

(3)根据另两轮测序实验确定的核苷酸突变位置和突变核苷酸信息是否一致来确定单核苷酸多态性;如果两次推断出的单核苷酸多态性位点的位置一致,则判定混合样本中包含有突变序列,并以两次计算的突变序列比例的平均值作为突变比例的比;根据两次推断出的可能的突变核苷酸的集合,取其交集作为检测出的突变核苷酸;如果两次推断出的单核苷酸多态性位点的位置不一致,则判定混合样本中不包含单核苷酸多态性。

有益效果:本发明相比于现有技术,具有更高的准确度;两核苷酸的循环添加顺序不需要经过复杂的实验设计,实验更容易操作;当混合样本中不包含单核苷酸多态性时,具有更好的鲁棒性以及更低的假阳性发现率;允许野生型序列未知,并能够确定单核苷酸突变位点的位置、突变核苷酸及其比例;成本低廉。

附图说明

图1是本发明所述的检测方法的逻辑流程图;

图2是进行两核苷酸实时合成测序所获得的三组信号谱;

图3是五种信号谱差别类型;

图4是三轮测序中,野生型样本与混合样本的信号谱。

具体实施方式

下面结合附图和实施例对本发明的技术方案作进一步的说明。

如图1所示是本发明所述的基于信号谱差异的混合样本单核苷酸多态性的检测方法,首先提取野生型样本的DNA和混合样本的DNA,进行聚合酶链式反应,即PCR扩展,得到单链PCR产物,利用一种两核苷酸实时合成技术对野生型样本和混合样本分别进行独立测序获得各自的信号谱。以待测序列“GATCGGTTCACGTC”为例,对野生型样本和混合样本的单链PCR产物进行三轮循环的两核苷酸实时合成测序(A+G)/(C+G)、(A+C)/(G+T)、(A+T)/(C+G),获得三组信号谱,如图2所示。

在实际实验中,信号强度与每次延伸的核苷酸个数成正比,并且可以经过标准化转换为延伸的核苷酸的个数。为了展示方便,本发明中直接使用延伸的核苷酸的个数来表示信号强度。

记R(reference)为野生型序列在单核苷酸多态性位点处的核苷酸,V(variant)为突变序列在单核苷酸多态性位点处的核苷酸,B(before)为野生型序列在单核苷酸多态性位点前一位的核苷酸,A(after)为野生型序列在单核苷酸多态性位点后一位的核苷酸。理论上,野生型与单碱基突变型序列的测试信号谱之间的差别仅有五种模式,如图3所示。

在给定一种两核苷酸添加方案的前提下,以AT/CG为例,定义A与T为一对双胞胎核苷酸(twin-nucleotide),C与G一对双胞胎核苷酸,也就是说,突变为A与T之间的突变或者C与G之间的突变,相邻的双胞胎核苷酸将会在同一步测序中一起延伸。如果R和V为一对双胞核苷酸,野生型序列与突变型序列的延伸步骤将完全一致,即信号谱将完全一致,定义这种情况下野生型序列与突变型序列的信号谱之间的差别模式为Type-0。

如果R和V不是双胞胎核苷酸,同时B和A也不是双胞胎核苷酸,那么在测序过程中,野生型序列与突变型序列的延伸步骤几乎完全一致。然而在B与A分别延伸的两步循环测序过程中,野生型序列与突变型序列所对应的信号强却存在着差别。具体来说,如果R和B是双胞胎核苷酸时,此时V和A也将成为双胞胎核苷酸,在这种情况下,在B延伸的那一步测序过程中,野生型序列所延伸的核苷酸将比突变型序列所延伸的核苷酸多一个,而在A延伸那一步测序过程中,野生型序列所延伸的核苷酸将比突变型序列所延伸的核苷酸少一个。对应的,定义这种情况下野生型与突变型序列的信号谱之间的差别模式为Type-1。反过来,如果R和A是双胞胎核苷酸时,此时V和B也将成为双胞胎核苷酸,定义这种情况下野生型与突变型序列的信号谱之间的差别模式为Type-2。

如果R和V不是双胞胎核苷酸,而B和A是双胞胎核苷酸,那么在测序过程中,在单核苷酸多态性位点之后,与野生型序列的延伸过程相比,突变型序列的延伸将会领先两步或者滞后两步。具体来说,当B和V是双胞胎核苷酸时,此时B和R将不是双胞胎核苷酸,因此B、V和A能够在同一步测序中延伸,而B、R和A却需要在三步循环测序过程才能延伸完成,所以突变型序列的延伸将会领先两步。对应的,定义这种情况下野生型序列与突变型序列的信号谱之间的差别模式为Type-3。反过来,当B和R是双胞胎核苷酸时,此时B和V将不是双胞胎核苷酸,突变型序列的延伸将会滞后两步,定义这种情况下野生型与突变型序列的信号谱之间的差别模式为Type-4。在Type-3型和Type-4型这两种情况下,野生型与突变型序列的信号谱可能在突变位置之后的多步循环测序中存在着差别。

根据三轮测序实验获得的野生型样本与混合样本的信号谱,针对每种两核苷酸添加方式,分别计算野生型样本与混合样本的信号谱之间的差异,使用枚举方法推断混合样本中可能的单核苷酸多态性位点的位置、比例和突变核苷酸。假设两核苷酸测序结果中,野生型样本的信号谱为W=(W1,W2,…,Wn),其中Wi为第i次添加两核苷酸时的信号强度。类似的,定义突变型序列的信号谱为M=(M1,M2,…,Mn),混合样本的测序信号谱为P=(P1,P2,…,Pn)。在实际实验过程中,需要使用延伸核苷酸个数为1的那一步循环测序所对应的峰值对信号谱进行标准化,以消除野生型样本与混合样本DNA浓度不一致所带来的影响。我们假设W和P是标准化之后的信号谱。混合样本由野生型序列和突变序列组成,其中突变序列的比例为a。理论上满足:

P=(1-a)W+aM 式1

研究显示,Wi的值实际上服从正态分布,如式2所示,类似的,混合样本的信号谱中,第i次添加两核苷酸时的信号强度Pi同样服从正态分布,如式3所示。

其中,和为信号强度的理论值,即理论上延伸的核苷酸个数。CV为信号强度的变异系数,能够直接反应焦磷酸测序实验的精准度。CV的具体值为标准差除以均值。在实际中,CV的值取决于测序实验的准确度。

在给定野生型样本的信号谱W和混合样本的信号谱P的前提下,定义差别谱D为:

D=P-W 式4

根据式1和式4可以推断出:

D=a(M-W) 式5

同样地,在第i次循环测序中,Di的值同样服从一个正态分布,基于式2-式5,可以推断出,Di服从如下分布:

在给定野生型样本信号谱W与混合样本信号谱P的前提下,使用枚举方法来推测突变发生的位置。具体来说,根据野生型序列与突变型序列信号谱之间的五种差别类型,假设野生型序列与突变型序列中不一致的延伸可能开始于任意一步循环测序中,枚举出野生型序列与突变型序列信号谱之间所有可能的差别谱。因为D是观察值并且是已知的,所以根据式5能够计算出混合样本中突变序列的比例。由于野生型序列与突变型序列信号谱在至少两步测序过程中存在着差别,所以会获得多个突变比例的计算值,取其平均值为最终的突变序列比例的计算值,并记为

定义一个得分Sp来评估所枚举的野生型序列与突变型序列信号谱之间每种差别的可能性,如式7所示。

其中,n为循环测序的步数。简而言之,Sp为差别谱中的每个元素对应的标准分(也被称为Z-score)的平均值。所以,Sp能够反应在假设的野生型序列与突变型序列信号谱的差别的前提下,观察值D与理论值P-W之间的偏差程度。对应地,较低的Sp值表示观察值D与理论值P-W之间的偏差较小,即所假设的野生型序列与突变型序列信号谱之间的差别谱具有更高的概率。

由于野生型序列和真实的突变比例a是未知的,导致和也都是未知的,根据式4和式5,可以推断出P-W等于a(M-W),在给定假设的野生型序列与突变型序列信号谱之间的差别谱,即M-W的前提下,可以使用来近似类似的,使用来近似此外,CV是一个取决于测序实验准确度的常量,所以修改Sp为式8,该值仍然能够反映观察值D与理论值P-W之间的偏差程度。

在给定野生型样本信号谱W和混合样本信号谱P的前提下,基于野生型序列与突变型序列信号谱之间的五种差别类型,枚举出野生型序列与突变型序列信号谱之间所有可能的差别谱,并计算突变型序列的比例。同时,对每种假设的差别谱进行打分,获得Sp值,该值能够间接反映各假设的差别谱的可能性。最终,选择具有最大Sp值的差别谱作为野生型序列与突变型序列信号谱之间可能的差别谱。

给定野生型序列与突变型序列信号谱之间差别谱,可以推断出单核苷酸多态性位点的位置为测序过程中野生型序列与突变型序列之间不一致延伸所开始的地方。假设不一致延伸发生在第j步循环测序反应过程中。显而易见,在第j-1步循环测序反应中所一致延伸的核苷酸个数能够直接根据野生序列的信号谱推断出来。而在第j步循环测序反应中所一致延伸的核苷酸个数,记为S,为和之间的较小值。理论上,如果差别谱属于Type-1型或Type-4型,S等于而如果差别谱属于Type-2型或Type-3型,S等于

对于Type-1型差别谱,的值等于而对于Type-4型差别谱,的值则可以通过式9计算得到。

因此,可以得出S的表达式,如式10,而最终计算的单核苷酸多态性位点的位置L的表达式,如式11。

由于野生型序列式未知的,所以在每步循环测试反应中真实延伸的核苷酸个数也是未知的。在实际的实验中,真实延伸的核苷酸个数是与信号强度成正相关但并非严格相等,而且每个对应的信号强度是处于一定的取值范围内的。反过来,给定信号强度,可以根据其值所处的区间范围直接推断出真实延伸的核苷酸个数。在这里,使用函数f(Wi)来标记根据信号强度Wi所推断而获得的核苷酸个数。

因此,式10可以修改为:

在式13中,使用来近似突变比例a的值,并使用Pj-Wj来近似的值。综合,最终推测出单核苷酸多态性位点的位置L如式14。

同时,推断出突变型序列与野生型序列在单核苷酸多态性位点处的可能核苷酸。假设不一致延伸发生在第j步循环测序反应过程中。对于Type-2和Type-3型差别谱,野生型序列在单核苷酸多态性位点处的可能核苷酸即为在第j+1步循环测序过程中添加的两种核苷酸,而突变序列在单核苷酸多态性位点处的可能核苷酸即为在第j步循环测序过程中添加的两种核苷酸。与此相反的是,对于Type-1和Type-4型差别谱,野生型序列在单核苷酸多态性位点处的可能核苷酸即为在第j步循环测序过程中添加的两种核苷酸,而突变序列在单核苷酸多态性位点处的可能核苷酸即为在第j+1步循环测序过程中添加的两种核苷酸。

综合分析三轮测序实验所检测出的可能的单核苷酸多态性位点,进行一致性判断,如果混合样本中含有单核苷酸多态性位点,归纳出其位置、突变核苷酸及其比例。

两核苷酸实时测序技术中,野生型样本和混合样本的信号谱均有三组,分别对应三种不同的两核苷酸组合添加方案。理论上,在三轮测序实验中,必定有一轮测序中,参考核苷酸R与突变核苷酸V属于一对双胞胎核苷酸,即野生型样本和混合样本的信号谱是完全一致的。因此给定三轮测序实验结果,首先要鉴定出在哪一轮测序实验中,野生型样本的信号谱与混合样本的信号谱是完全一致的。

当核苷酸R与核苷酸V属于一对双胞胎核苷酸时,野生型样本的信号谱与混合样本的信号谱是完全一致的,即差别谱中的每个元素的值均为0,与式8类似,定义一个打分Sc,如式15,来评价野生型样本的信号谱与混合样本的信号谱之间的一致性。

其中,n为循环测序的步数。具体而言,Sc是Sp的一个特殊情况,即假设的野生型序列信号谱与突变型序列信号谱的差别谱中所有的元素的值均为0。

当假设突变不会导致野生型序列和突变型序列具有不同的信号谱时,Sc可以反应观测到的差别谱D与理论上的差别谱P-W的偏差程度,而当假设突变会导致野生型序列和突变型序列具有不同的信号谱时,Sp可以反应观测到的差别谱D与理论上的差别谱P-W的偏差程度。最终,定义Sm,如式16,来评估混合样本中存在着突变型序列的可能性。

Sm=Sc-Sp 式16

当突变序列存在时,Sc的值会比Sp更大。而且突变序列的比例越高,Sm的值将会更大。直观上,Sm能够反映突变序列的存在对野生型样本信号谱与混合样本信号谱之间的一致性的影响。因此推测具有最低Sm值的一轮测序实验中,野生型序列与突变型序列的信号谱完全一致。

在三轮测序完成之后,针对每轮测序实验的结果,分别利用枚举方法检测可能的单核苷酸多态性位点,同时计算对应的Sm。认定具有最低Sm的那一轮测序实验中,野生型样本和混合样本的信号谱是完全一致的。那么,在该轮测序实验中,参考核苷酸R与突变核苷酸V为一对双胞胎核苷酸。

根据另外两轮测序实验中所检测出来可能的单核苷酸多态性位点的位置,进行一致性判断。如果两次推断出的单核苷酸多态性位点的位置一致,则认为混合样本中包含有突变序列,并以两次计算的突变序列比例的平均值作为突变比例的最终计算值,根据两次推断出的可能的突变核苷酸的集合,取其交集作为最终检测出的突变核苷酸。反之,如果两次推断出的单核苷酸多态性位点的位置不一致,则认为混合样本中不包含单核苷酸多态性。

下面以一个实际的测试实验为例进行说明。

基于信号谱差异推断拟定的混合样本中的单核苷酸多态性位点,该混合样本含野生型序列“CGACCAGCT”和突变型序列“CGATCAGCT”。

(1)将100%野生型序列和由90%的野生型序列与10%的突变型序列组成的混合样本,独立经过三轮两核苷酸实时合成测序,分别模拟获得野生型样本与混合样本的测序信号谱,结果如表1及图4所示。设定测序过程中,信号变异系数CV值为0.001。

(2)每轮两核苷酸实时合成测序中每步循环测序反应信号强度的结果如表1所示。

表1

(3)根据式4、5、8、13、14、15、16。针对每轮两核苷酸测序结果,首先计算野生型样本与混合样本的信号谱之间的差异,在此基础之上,采用枚举方法推断可能的单核苷酸多态性位点的位置、突变序列的比例和突变核苷酸,同时计算Sm值,结果如表2所示。

表2

(4)根据步骤(3)中表2的结果,第二轮AG/CT两核苷酸测序的Sm值在三轮测序中最低,可以确定这一轮测序实验中,野生型序列与突变型序列信号谱之间的差别类型为Type-0。在第二轮测序中,添加的两核苷酸分别为AG与CT,据此可以推断突变为A与G之间的突变或者C与T之间的突变。值得注意的是,根据第二轮测序结果,本发明方法推断混合样本中可能存在一个单核苷酸多态性位点,其位于第3位核苷酸处,对应的突变序列比例为0.0043,显而易见,这是由于测序过程中引入的随机误差所导致的。

根据第一轮和第三轮两核苷酸测序结果独立推断出的可能的单核苷酸多态性位点的位置均为第4位核苷酸,因此可以判定混合样本中包含有突变序列、并确定单核苷酸多态性位点发生在第4位核苷酸处。根据第一轮和第三轮测序计算混合样本中的突变序列的比例分别为0.1006和0.0999,因此最终计算出混合样本中突变序列的比例为0.10025。

根据第一轮测序结果,确定野生型序列在单核苷酸多态性位点处核苷酸可能为C或G,而根据第三轮测序结果,本发明方法检测出野生型序列在单核苷酸多态性位点处核苷酸可能为A或C,因此可以归纳出野生型序列在单核苷酸多态性位点处核苷酸为C。

根据第一轮测序结果,确定突变序列在单核苷酸多态性位点处核苷酸可能为A或T,而根据第三轮测序结果,本发明方法检测出突变序列在单核苷酸多态性位点处核苷酸可能为G或T,因此可以归纳出突变序列在单核苷酸多态性位点处核苷酸为T。

当前第1页1 2 3 
网友询问留言 已有0条留言
  • 还没有人留言评论。精彩留言会获得点赞!
1