一代基因测序中自动识别杂合突变的方法与流程

本发明涉及一种基因测序结果分析方法，尤其是一种可自动识别其中杂合突变的方法，具体的说是一代基因测序中自动识别杂合突变的方法。

背景技术：

现有的一代基因测序（sanger测序法）结果缺乏成熟的杂合突变自动化识别技术，均须通过人工识别，导致杂合突变识别的效率及准确率低下，而人力及时间成本则非常高昂.严重影响了基因测序的发展，亟须改进。

技术实现要素：

本发明的目的是针对现有技术的不足，提供一种一代基因测序中自动识别杂合突变的方法，大大提升杂合突变的识别效率及准确率，降低人工及时间成本，具有广泛的应用前景

本发明的技术方案是：

一种一代基因测序中自动识别杂合突变的方法，包括以下步骤：

s1将一代基因测序结果中各个碱基的检出信号数据转换为横坐标为正实数点位，纵坐标为对应碱基检出信号强度值的坐标序列；其中，所述正实数点位的单位为bp；

s2基于所述各个碱基的坐标序列，计算各个碱基在各整数点位左右各0.5bp区间中的碱基检出信号强度点曲线与横坐标轴所围图形的面积，作为各整数点位上各个碱基检出信号强度包络面积；

s3以各整数点位上各碱基检出信号的最大包络面积作为分类依据，将各整数点位分组，使得每个分组内各碱基检出信号的最大包络面积数值相接近；

s4依据杂合突变基因中出现频率的统计特征及发生杂合突变的整数点位左右各0.5bp区间中的图形特征，在所述的每个整数点位分组内，各自确定用于识别杂合突变的面积差阈值和面积比阈值；

s5使用噪声识别算法标识出被测序列中的噪声整数点位集合；

s6针对不在噪声整数点位集合中的每个整数点位，根据其所在的整数点位分组确定面积差阈值和面积比阈值；根据该整数点位各个碱基检出信号强度包络面积，计算该整数点位的面积最大的两种碱基的面积差和面积比，若其面积差小于面积差阈值且面积比小于面积比阈值，则该点位判定为疑似杂合突变。

优选的，所述步骤s1包括以下步骤：

s1.1读取各类碱基对应的检出信号强度值所构成的正有理数序列，以下简称为检出序列；

s1.2读取当前整数点位所在的检出序列下标；所述点位的单位为bp；所述下标为正整数；

s1.3求当前整数点位与前一整数点位中的检出序列的下标差，取其倒数作为前一整数点位到当前整数点位的数据点位，简称为步进单位值；该步进单位值为正有理数；单位为bp；

s1.4横坐标从前一点位开始，按下标顺序取出检出序列数据，该数据为正有理数；

s1.5将检出序列数据作为纵坐标，以累加s1.3求出的步进单位值作为其横坐标，直至下个整数点位前；将这些横坐标和纵坐标对应的点所构成的序列推入坐标序列末尾；

s1.6若当前整数点位尚未超过被测序列长度，则将当前整数点位加一后返回s1.3步骤继续。若当前整数点位达到被测序列长度，则此时就获得了完整的坐标序列。

优选的，所述步骤s2包括以下步骤：

s2.1在当前整数点位上对于各个碱基的两侧各0.5bp区间截取序列，简称为截取序列；

s2.2在s2.1中当前整数点位两侧各0.5bp的两端点无数据时，可通过插值算法补正，提升包络面积计算精度；

s2.3计算各整数点位截取序列中每两组坐标构成的梯形面积，将其累加即得到各类碱基于各点位的包络面积，并将其推入面积序列末尾；

s2.4若当前整数点位尚未超过被测序列长度，则将当前整数点位加一后返回s2.1步骤继续；若当前整数点位达到被测序列长度，则此时就获得了完整的面积序列。

优选的，所述步骤s4包括以下步骤：

s4.1将每个点位处的各类碱基面积的最大值减去第二大值，作为最大包络面积差序列，简称为面积差序列；

s4.2将每个点位处的各类碱基面积的最大值除以第二大值，作为最大包络面积比序列，简称为面积比序列；

s4.3计算s4.1面积差序列的中位数，简称为面积差中位数；计算s4.2面积比序列的中位数，简称为面积比中位数；

s4.4面积差阈值等于面积差中位数乘以分组内波形峰值最大变化率；面积比阈值等于面积比中位数乘以分组内波形峰值最大变化率。

优选的，所述步骤s6包括以下步骤：

s6.1若当前点位在噪声点位中，则跳过此点位，继续分析下一个点位；

s6.2根据当前点位所在的按最大碱基包络面积，获取其所在整数点位分组，从而获知该分组的面积差阈值及面积比阈值；

s6.3将分析点位的面积差和面积比自序列取出，若分析点位的面积差小于等于面积差阈值，且分析点位的面积比小于等于面积比阈值，则认为该点位为疑似杂合突变点位。

本发明的有益效果：

本发明设计合理，使用方便，可通过利用计算机图形学方法来仿造人类观察杂合突变的过程，大大提升杂合突变的识别效率及准确率，降低人工及时间成本，具有广泛的应用前景。

附图说明

图1是不同碱基的连线图碱基荧光信号强度曲线示意图。

图2是图1的局部放大图。

图3是两组坐标点所围图形面积。

图4是每个整数点位范围内包络线与点位轴所围图形面积。

其中：

图1中顶端数字为整数点位，曲线上的字母为该数据点对应的碱基类型，下方的a、t、c、g字母为该整数点位识别出的碱基类型，最下方的数字为phred得分，用以评价此处测序质量。

图2中框出的313整数点位出现了c和t碱基的杂合突变特征信号曲线，因为该两种碱基对应的包络线与点位轴面积差和面积比接近。

图3是使用梯形面积公式计算两组坐标点所围图形面积。

图4为将某一类碱基对应碱基荧光曲线的每个整数点位范围内的多个梯形面积进行计算，并求和得到每个整数点位内的该类碱基对应碱基荧光曲线与点位轴的包络面积。

具体实施方式

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

如图1、2、3和4所示。

一种一代基因测序中自动识别杂合突变的方法，包括以下步骤：

s1将一代基因测序结果中各个碱基的检出信号数据转换为横坐标为正实数点位，其中正整数点位显示在图1、图2的上方，纵坐标为对应碱基检出信号强度值的坐标序列，因为在每个正实数点位的横坐标上会同时出现多种碱基检出信号，故在图1、图2所示的波形部分用对应碱基的字母进行标出；其中，所述正实数点位的单位为bp，例如图2框出的点位为313，其内部波形所对应的碱基检出信号主要由c和t构成；

s1.1本实施例以人类的lpl蛋白对应的核苷酸序列作为示例，首先读取各类碱基对应的检出信号强度值所构成的正有理数序列，例如在某次测序中a碱基对应的检出信号强度值序列为[50,99,203,389,679,3455,3816,4172…]，这个序列包含了5921个采样点上的检出信号强度值，以下简称为检出序列；

s1.2读取当前整数点位所在的检出序列下标；所述点位的单位为bp；所述下标为正整数，例如[3,39,44,58,76,94,101,109…]，这个序列包含了494个整数点位所对应的检出信号强度值数据下标，如点位1，它对应的数据点应该是下标为0，1，2的三条数据，在a碱基上也就是50，99和203；

s1.3求当前整数点位与前一整数点位中的检出序列的下标差，取其倒数作为前一整数点位到当前整数点位的数据点位，简称为步进单位值；该步进单位值为正有理数；单位为bp，接s1.1和s1.2，例如2点位到1点位的下标差为39-3=36，说明中间有36个均匀的样点，步进值为1/36，1点位比较特殊，其下标差与0相减为3-0=3，步进值为1/3；

s1.4横坐标从前一点位开始，按下标顺序取出检出序列数据，该数据为正有理数，接s1.3中例，其中整数点位1对应的检出信号强度值为50、99、203共计3个样点，整数点位2对应的检出信号强度值为389,679,3455,3816,4172…共计36个样点；

s1.5将检出序列数据作为纵坐标，以累加s1.3求出的步进单位值作为其横坐标，直至下个整数点位前；将这些横坐标和纵坐标对应的点所构成的序列推入坐标序列末尾，接s1.4，在整数点位1的三个样点中，其横坐标自0起，每个数据点累计1/3，从而构成坐标序列(1/3,50)、(2/3,99)，（3/3,203）,在整数点位2的36个样点中其横坐标自1起每个数据点累计1/36，从而构成坐标序列(1+1/36,389)、(1+2/36,679)、(1+3/36,3455)…，以此类推；

s1.6若当前整数点位尚未超过被测序列长度，则将当前整数点位加一后返回s1.3步骤继续。若当前整数点位达到被测序列长度，则此时就获得了完整的坐标序列；

s2基于所述各个碱基的坐标序列，计算各个碱基在各整数点位左右各0.5bp区间中的碱基检出信号强度点曲线与横坐标轴所围图形的面积，作为各整数点位上各个碱基检出信号强度包络面积；

s2.1在当前整数点位上对于各个碱基的两侧各0.5bp区间截取序列，简称为截取序列,接s1.4，使用整数点位1两侧各0.5bp区间即[0.5bp,1.5bp]，则整数点位1处截取序列为(2/3,99)，（3/3,203），(1+1/36,389)、(1+2/36,679)、(1+3/36,3455)…至横坐标1+18/36为止；

s2.2在s2.1中当前整数点位两侧各0.5bp的两端点无数据时，可通过插值算法补正，提升包络面积计算精度，以s2.1中数据为例，以1位当前整数点位，在右侧1.5bp为止有数据点，但是在左侧0.5bp为止恰好没有数据点，所以此时可以用斜率方程求得0.5bp处由(1/3,50)和(2/3,99)构成点的斜率，再由此斜率求得0.5bp处上的插值点为74，从而得到(0.5,74)；

s2.3如图3所示，利用梯形面积公式计算各整数点位截取序列中每两组坐标构成的梯形面积，将其累加即得到如图4所示的各类碱基于各点位的包络面积，并将其推入面积序列末尾，以s2.2中补正过的数据为例，(0.5,74)，(2/3,99)两点所围面积可以用(上底+下底)*高/2.0求得，也就是(74+99)x(2/3–0.5)/2.0=14；

s2.4若当前整数点位尚未超过被测序列长度，则将当前整数点位加一后返回s2.1步骤继续；若当前整数点位达到被测序列长度，则此时就获得了完整的面积序列；

s3以各整数点位上各碱基检出信号的最大包络面积作为分类依据，将各整数点位分组，使得每个分组内各碱基检出信号的最大包络面积数值相接近；常用的分组方法包括k-means，svm等；我们可以从图1和图2中看出，各碱基检出信号的最大包络峰值是起伏不定的，通过这种方法可以减少这种起伏对我们判断杂合突变带来的影响，从而提升正确率；以atcg碱基在1-4整数点位左右各0.5范围内的包络面积序列：(0,100,0,0),(98,0,91,0),(10,2,3,1),(2,50,49,1),(4,2,30,1),(0,8,0,8),…为例，其最大包络面积序列为(100,98,10,50,30,8)，如果分为高中低三组，则高组包含(0,100,0,0),(98,0,91,0),…中组包含(2,50,49,1),(4,2,30,1)，…低组包含(10,2,3,1),(0,8,0,8)…。

s4依据杂合突变基因中出现频率的统计特征及发生杂合突变的整数点位左右各0.5bp区间中的图形特征，在所述的每个整数点位分组内，各自确定用于识别杂合突变的面积差阈值和面积比阈值；

s4.1将每个点位处的各类碱基面积的最大值减去第二大值，作为最大包络面积差序列，简称为面积差序列；接s3，，其中的高组面积差序列为(100,7,…)，中组面积差序列为(1,28,…)，低组面积差序列为(7,0,…)；

s4.2将每个点位处的各类碱基面积的最大值除以第二大值，作为最大包络面积比序列，简称为面积比序列；接s3，其中的高组面积比序列为(正无穷,98/91,…)，中组面积差序列为(50/49,30/4,…)，低组面积差序列为(10/3,8/8,…)；

s4.3计算s4.1面积差序列的中位数，简称为面积差中位数，接s4.1，其中位数高组96、中组24，低组8；计算s4.2面积比序列的中位数，简称为面积比中位数，接s4.1，其中位数高组999、中组777，低组666；

s4.4面积差阈值等于面积差中位数乘以分组内波形峰值最大变化率；面积比阈值等于面积比中位数乘以分组内波形峰值最大变化率，此处可设各组内的波形峰值最大变化率为20%，面积差阈值分别为高组96x20%=19.2，中组24x20%=4.8和低组8x20%=1.6，面积比阈值分比为高组400x20%=80、中组280x20%=56，低组80x20%=16；

s5使用噪声识别算法标识出被测序列中的噪声整数点位集合；

可用基于测序质量分（phred分）的modifiedmotttrimming算法寻找到phred分数最大的低噪声序列，该序列两侧即为首尾噪声序列，其中元素均为噪声点位；

可用小波变换或傅里叶变换去除坐标序列中的噪音频谱；

可用基于神经网络的方法对于每个点位的图形质量进行评估，将评估中图形质量较差的识别为噪声；

也可用其他的噪声识别算法；

s6针对不在噪声整数点位集合中的每个整数点位，根据其所在的整数点位分组确定面积差阈值和面积比阈值；根据该整数点位各个碱基检出信号强度包络面积，计算该整数点位的面积最大的两种碱基的面积差和面积比，若其面积差小于面积差阈值且面积比小于面积比阈值，则该点位判定为疑似杂合突变，从图2中可以看出，在发生杂合突变的点位上面积差较小，面积比也较为接近于1，符合视觉规律；

s6.1若当前点位在噪声点位中，则跳过此点位，继续分析下一个点位，例如，噪声点位集合为(1,2,3,412,413)，当前点位为1，则直接跳过，因为2、3也是噪声点位，故从点位4开始分析；

s6.2根据当前点位所在的按最大碱基包络面积，获取其所在整数点位分组，从而获知该分组的面积差阈值及面积比阈值，如点位4、5和6的各类碱基荧光信号包络面积序列为(2,50,49,1),(4,2,30,1),(0,8,0,8)，其中(2,50,49,1),(4,2,30,1)都是属于中组，则它们的面积差阈值为4.8，面积比阈值为56，其中(0,8,0,8)属于低组，则它的面积差阈值为1.6，面积比阈值为16；

s6.3将分析点位的面积差和面积比自序列取出，若分析点位的面积差小于等于面积差阈值，且分析点位的面积比小于等于面积比阈值，则认为该点位为疑似杂合突变点位，接s6.2例点位4、5、6上的面积差序列为1、26、0，其中4、6符合小于面积差阈值条件，面积比序列为50/49，30/4和1，均小于对应阈值，故同时符合两个条件的整数点位4、6被识别为杂合突变点位。

本发明通过利用计算机图形学原理来仿造人类观察杂合突变的过程，大大提升了杂合突变的识别效率及准确率，降低人工及时间成本，因而，具有广泛的应用前景。

本发明未涉及部分均与现有技术相同或可采用现有技术加以实现。