一种基于二代测序技术的单基因或多基因拷贝数检测***及
方法
技术领域
本发明属于生物医药领域,具体是二代测序领域中的拷贝数变异检测相关的数据分析方法,尤其涉及一种基于二代测序技术的单基因或多基因拷贝数检测***;此外,本发明还涉及一种基于二代测序技术的单基因或多基因拷贝数检测方法。
背景技术
拷贝数变异是很多遗传病和罕见病发生发展的原因。比如α地中海贫血、脊髓性肌萎缩症、杜氏肌营养不良症等等,均是由一个或多个基因上的外显子发生了拷贝数异常导致的。因此,这些基因的拷贝数检测对遗传病的控制和出生缺陷的预防起着至关重要的作用。传统的检测方法比如qPCR(实时荧光定量核酸扩增检测***),MLPA(多重连接探针扩增技术),基因芯片等存在通量低,成本高,或者一次检测的基因数和变异类型有限的缺点。二代测序方法弥补了传统技术通量不足的特点,且可以一次性检测包括单核苷酸变异、短序列***缺失、拷贝数变异、结构变异在内的多种变异类型,极大提升了可检测变异范围,具有非常广阔的临床应用前景。
然而,二代测序方法在检测基因拷贝数变异方面存在着诸多挑战。比如,建库试剂盒品牌、批次之间,甚至不同时期的实验批次之间都存在着诸多差异和批次效应,这使得常用的基于测序深度来推断拷贝数变异的方法在实际使用过程中往往得不到很高的检测精度。另外由于基因之间的差异性也使得各计算方法之间难以通用,往往需要结合不同的方法才能达到比较好的精度。另外,对于罕见病以及发病率比较低的遗传病,由于阳性样本稀缺,也难以得到合适的数据做方法开发。为了解决或者避开这些问题,本发明方法基于LASSO的机器学习方法,利用大规模正常人群数据作为训练集得到不依赖于批次和试剂盒品牌的特征区间以及系数参数,结合层次转移模型(shifting level model)来推断阳性样本在单基因上的拷贝数状态,精度可达单一外显子,具有较高的灵敏度和特异性。
目前,LASSO机器学习方法尚未用于二代测序方法中单基因或多基因的拷贝数变异检测。LASSO回归是一种正则化线性模型,其特点是在拟合广义线性模型的同时进行变量筛选和复杂度调整。LASSO回归通过一系列参数控制模型,有选择性把变量放入模型中,不仅可以解决过拟合问题,而且在参数缩减过程中,将一些重复的没必要的参数直接缩减为零,达到提取有用特征的作用。对于线性模型来说,复杂度与模型的变量数有直接关系,变量数越多,模型复杂度越高。更多的变量在拟合时往往给出一个更好的模型,但是同时也面临过度拟合的危险。而LASSO回归复杂度调整的程度由参数λ来控制,λ越大对变量较多的线性模型的惩罚力度就越大,从而最终获得一个变量较少的模型。LASSO线性回归很容易通过随机梯度下降来更新数据模型,通过正则化来避免过度拟合,特别适用于参数数目缩减与参数的选择,因而用来估计稀疏参数的线性模型。目前,LASSO回归广泛应用于计算机、互联网、电信、临床医学、经济学、金融投资、生物学等各个领域。利用二代测序技术来检测拷贝数变异受批次效应等测序噪音影响很大,而LASSO线性回归的变量筛选和复杂度调整特征特别适用于这种变量特别多的数据类型。
中国发明专利申请CN2011100654064公开了检测拷贝数变异的方法和装置。该方法主要是基于全基因组测序数据,检测癌症和癌旁配对样本的拷贝数变异,局限于检测大片段的拷贝数变异。其缺点在于:需要配对样本,依赖于批次效应和试剂盒品牌效应,检测精度有限,适用范围局限于检测大片段的拷贝数变异。而本发明不需要配对样本,是通过对大量正常样本进行训练,利用机器学***;适用的范围更广,适用于所有二代测序技术的数据类型;灵敏度更高,不依赖于批次效应和试剂盒品牌效应。
发明内容
本发明所要解决的技术问题:本发明公开了一种基于二代测序技术的单基因或多基因拷贝数检测***及方法,解决了二代测序技术领域单基因或多基因拷贝数变异检测的批次效应强、检测精度低的问题。LASSO线性回归方法的优势在于通过对于大样本训练的计算资源和时间花费比较少,且适用于基因组学数据这种变量特别多的数据类型。通过LASSO训练可以筛选出一些不依赖于实验和批次效应的变量区域。结合层次转移模型,可以对目标区域做更好的矫正。此方法可广泛用于多种全外显子组或者基因面板的捕获数据,适用于单基因或者多基因的拷贝数状态检测。
本发明解决技术问题所采用的技术方案是:
在本发明的一方面,提供一种基于二代测序技术的单基因或多基因拷贝数检测***,包括:
序列比对模块,用于将测序读段使用比对软件比对到人参考基因组;
去除重复序列模块,与序列比对模块相连,用于去除PCR扩增过程产生的重复序列;
计算覆盖深度模块,与去除重复序列模块相连,用于将靶向捕获区间分成多个片段,计算i样本j片段上的平均覆盖深度dij;
标准化覆盖深度模块,与计算覆盖深度模块相连,用于将i样本j片段上的平均覆盖深度dij根据样本的平均测序深度meanDi进行标准化;
正则化线性回归训练模块,与标准化覆盖深度模块相连,利用正则化线性回归模型在目标基因拷贝数正常的样本上进行参数训练;
覆盖深度预测模块,与正则化线性回归训练模块相连,用于预测待测样本各目标片段上的覆盖深度期望值
断点检测和log2Ratio值矫正模块,与覆盖深度预测模块相连,用于根据待测样本在目标基因待测片段上的标准化覆盖深度实际值得到log2Ratio值;并利用层次转移模型,进一步消除测序过程中的误差以及对CNV拷贝数的断点进行识别及log2Ratio值矫正;
拷贝数状态推断模块,与断点检测和log2Ratio值矫正模块相连,用于设置log2Ratio经验阈值来推断拷贝数状态。
作为本发明优选的技术方案,所述计算覆盖深度模块中,将靶向捕获区间分成多个片段,每个片段的长度为100-150bp。
作为本发明优选的技术方案,所述标准化覆盖深度模块中,将i样本j片段上的平均覆盖深度dij根据样本的平均测序深度meanDi进行标准化,具体为:
定义xij为非待测基因的其他所有片段上的标准化覆盖深度,yij为目标基因片段上的标准化覆盖深度:
所述正则化线性回归模型包括LASSO回归、岭回归和弹性网络回归等,本发明中所述正则化线性回归训练模块优选为LASSO线性回归训练模块,所述正则化线性回归模型优选为LASSO线性回归模型。
作为本发明优选的技术方案,所述LASSO线性回归训练模块中,利用LASSO线性回归模型在目标基因拷贝数正常的样本上进行参数训练,具体为:
假定有m个训练样本,每个样本n个片段:
其中为各训练样本在非目标基因上的所有片段的标准化覆盖深度矩阵,为各训练样本在目标基因外显子或检测片段上的标准化覆盖深度矩阵;βtarget为系数矩阵;经10倍交叉验证训练得到损失函数最小时的超参数λ以及对应的系数矩阵βtarget;对于每一个待检基因target,均有一个对应的拟合系数矩阵βtarget。
作为本发明优选的技术方案,所述覆盖深度预测模块中,用于预测待测样本各目标片段上的覆盖深度期望值具体为:有待测样本p在非目标基因上的所有片段的标准化覆盖深度矩阵根据系数矩阵βtarget预测待测样本各目标片段上的覆盖深度期望值
作为本发明优选的技术方案,所述断点检测和log2Ratio值矫正模块中,所述根据待测样本在目标基因待测片段上的标准化覆盖深度实际值得到log2Ratio值,具体为:并利用层次转移模型,进一步消除测序过程中的误差以及对CNV拷贝数的断点进行识别及log2Ratio值矫正。
作为本发明优选的技术方案,所述拷贝数状态推断模块中,设置log2Ratio经验阈值来推断拷贝数状态,具体为:
在本发明的另一方面,还提供一种所述的基于二代测序技术的单基因或多基因拷贝数检测***的实现方法,该方法利用基于正则化线性回归模型的机器学习算法(优选为LASSO线性回归模型)以及层次转移模型推断单基因或者多基因外显子的拷贝数变异,具体步骤如下:
步骤1:序列比对:将测序得到的读段比对使用比对软件比对到人参考基因组;
步骤2:去除重复序列:将对比对得到的结果去除PCR扩增过程产生的重复序列;
步骤3:将靶向捕获区间分成特定长度的片段,计算i样本j片段上的平均覆盖深度dij;
步骤4:覆盖深度标准化:将i样本j片段上的平均覆盖深度dij根据样本的平均测序深度meanDi进行标准化;
步骤5:训练模型:利用正则化线性回归模型(优选为LASSO线性回归模型),在目标基因拷贝数正常的样本上进行参数训练;
步骤6:覆盖深度预测:有待测样本p在非目标基因上的所有片段的标准化覆盖深度矩阵根据系数矩阵βtarget预测待测样本各目标片段上的覆盖深度期望值
步骤7:断点检测和log2Ratio值矫正:对于待测样本在目标基因待测片段上的标准化覆盖深度实际值得到log2Ratio值:
利用层次转移模型,进一步消除测序过程中的误差以及对CNV拷贝数的断点进行识别及log2Ratio值矫正;
步骤8:拷贝数状态推断:设置log2Ratio经验阈值来推断拷贝数状态。
作为本发明优选的技术方案,步骤3中,所述特定长度为100-150bp。
与现有技术相比,本发明的有益效果在于:本发明解决了二代测序技术领域单基因或多基因拷贝数变异检测的批次效应强、检测精度低的问题。LASSO线性回归方法的优势在于通过对于大样本训练的计算资源和时间花费比较少,且适用于基因组学数据这种变量特别多的数据类型。通过LASSO训练可以筛选出一些不依赖于实验和批次效应的变量区域。结合层次转移模型,可以对目标区域做更好的矫正。此方法可广泛用于多种全外显子组或者基因面板的捕获数据,适用于单基因或者多基因的拷贝数状态检测。
本发明的关键点是利用机器学习方法对大规模二代测序数据进行训练,结合层次转移模型,目的是降低由于批次效应造成的技术和生物学误差,从而达到更好的拷贝数检测的准确性和精度。目前,LASSO机器学习方法尚未用于二代测序方法中拷贝数变异检测,LASSO机器学习模型也可用其他同类型模型替代,包括但不限于岭回归(RidgeRegression),弹性网络回归(ElasticNet Regression)等。LASSO回归、岭回归和弹性网络回归都是属于正则化线性回归模型,在数据集中变量之间的高纬度和多重线性的情况下都能很好地工作。岭回归在拟合数据的同时使模型权重尽可能小,它会缩小系数的值,但不会达到零,即其没有特征选择功能。相较于岭回归,LASSO回归倾向于完全消除不重要的权重,其在拟合模型的同时进行变量筛选和复杂度调整。弹性网络回归是在LASSO回归和岭回归中进行了折中,用于只有很少的权重非零的系数模型,可以在旋转下继承岭回归的一些稳定性。当多个特征和另一个特征相关的时候,LASSO回归效果并不理想,它只能随机的将变量组中的某一变量选入模型,其他的变量归零,而弹性网络回归所选变量的数量没有限制。如果只有少部分特征是有用的情况下,倾向于选择LASSO回归或弹性网络回归。本发明***及方法不但可以用于单基因拷贝数变异检测的结果,还可推广到多基因拷贝数检测。
附图说明
图1是本发明基于二代测序技术的单基因或多基因拷贝数检测***及方法的流程图。
图2是采用本发明方法为单基因拷贝数变异检测的实例结果图。其中,图a)表示3个女性样本DMD基因部分外显子发生拷贝数变异的检测结果;图b)表示3个男性样本DMD基因部分外显子发生拷贝数变异的检测结果;图c)表示2个女性样本SLC26A4基因部分外显子发生拷贝数变异的检测结果;图d)表示1个男性样本PTPRQ基因部分外显子发生拷贝数变异的检测结果。
具体实施方式
以下通过实施例对本发明进行更加具体的说明。应当理解,此处所描述的实施例是用于解释本发明,而非用于限定本发明。
如图1所示,本发明基于二代测序技术的单基因或多基因拷贝数检测***,包括:
序列比对模块,用于将测序读段使用比对软件比对到人参考基因组;
去除重复序列模块,与序列比对模块相连,用于去除PCR扩增过程产生的重复序列;
计算覆盖深度模块,与去除重复序列模块相连,用于将靶向捕获区间分成多个片段(每个片段的长度为100-150bp),计算i样本j片段上的平均覆盖深度dij;
标准化覆盖深度模块,与计算覆盖深度模块相连,用于将i样本j片段上的平均覆盖深度dij根据样本的平均测序深度meanDi进行标准化;
LASSO线性回归训练模块,与标准化覆盖深度模块相连,利用LASSO线性回归模型在目标基因拷贝数正常的样本上进行参数训练;
覆盖深度预测模块,与LASSO线性回归训练模块相连,用于预测待测样本各目标片段上的覆盖深度期望值
断点检测和log2Ratio值矫正模块,与覆盖深度预测模块相连,用于根据待测样本在目标基因待测片段上的标准化覆盖深度实际值得到log2Ratio值;并利用层次转移模型,进一步消除测序过程中的误差以及对CNV拷贝数的断点进行识别及log2Ratio值矫正;
拷贝数状态推断模块,与断点检测和log2Ratio值矫正模块相连,用于设置log2Ratio经验阈值来推断拷贝数状态。
如图1所示,本发明基于二代测序技术的单基因或多基因拷贝数检测方法利用基于LASSO的机器学习算法以及层次转移模型推断单基因或者多基因外显子的拷贝数变异,具体步骤如下:
步骤1:序列比对。将测序得到的读段比对使用比对软件BWA比对到人参考基因组;
步骤2:去除重复序列。与序列比对模块相连,用Picard软件将对比对得到的结果去除PCR扩增过程产生的重复序列;
步骤3:将靶向捕获区间分成特定长度100-150bp的片段,计算i样本j片段上的平均覆盖深度dij;
步骤4:覆盖深度标准化。将i样本j片段上的平均覆盖深度dij根据样本的平均测序深度meanDi进行标准化,定义xij为非待测基因的其他所有片段上的标准化覆盖深度,yij为目标基因片段上的标准化覆盖深度:
步骤5:训练模型。利用LASSO线性回归模型,在目标基因拷贝数正常的样本上进行参数训练。假定有m个训练样本,每个样本n个片段:
其中为各训练样本在非目标基因上的所有片段的标准化覆盖深度矩阵,为各训练样本在目标基因外显子或检测片段上的标准化覆盖深度矩阵。βtarget为系数矩阵。训练得到损失函数最小时的超参数λ以及对应的系数矩阵βtarget。对于每一个待检基因target,均有一个对应的拟合系数矩阵βtarget。
步骤6:覆盖深度预测。有待测样本p在非目标基因上的所有片段的标准化覆盖深度矩阵根据系数矩阵βtarget预测待测样本各目标片段上的覆盖深度期望值
步骤7:断点检测和log2Ratio值矫正。对于待测样本在目标基因待测片段上的标准化覆盖深度实际值可以得到log2Ratio值:
利用层次转移模型,进一步消除测序过程中的误差以及对CNV拷贝数的断点进行识别及log2Ratio值矫正。层次转移模型(Shifting level models,SLM)是一种用于评估未知断点和均值的迭代分割(segmentation)算法,可以模拟噪声序列过程,用于显示均值的骤变。这一模型最早用于分析微阵列比较基因组杂交(arrayCGH)的log2Ratio基因组图谱,也适用于对二代测序得到的基因组图谱进行分割,对基因组上的拷贝数增加和缺失的断点进行识别。
步骤8:拷贝数状态推断。设置log2Ratio经验阈值来推断拷贝数状态:
实施例1采用本发明方法为单基因拷贝数变异检测
分别使用500个正常男性和500个正常女性的数据作用训练集,来训练LASSO线性回归模型。待测阳性样本包括6个DMD基因部分缺失或增加样本,2个SLC26A4基因部分缺失样本,和1个PTRRA基因部分缺失样本。样本拷贝数变异状态详细信息如下表1,9个测试样本的一致性达到100%。
表1
检测结果如图2所示,图2中横坐标表示划分的捕获区间划分的区间,纵坐标表示log2Ratio,每个圆点表示一个区间的log2Ratio值,圆点上的折线表示矫正之后的log2Ratio值,上下两条虚线分别表示拷贝数增加和缺失的log2Ratio阈值,如果圆点上的折线高于上虚线,则表明该区间所对应的基因片段发生了拷贝数增加,而过圆点上的折线低于下虚线,则表明该区间所对应的基因片段发生了拷贝数缺失,虚线间的折线则表明该区间所对应的基因片段上拷贝数正常。图a)表示3个女性样本DMD基因部分外显子发生拷贝数变异的检测结果;图b)表示3个男性样本DMD基因部分外显子发生拷贝数变异的检测结果;图c)表示2个女性样本SLC26A4基因部分外显子发生拷贝数变异的检测结果;图d)表示1个男性样本PTPRQ基因部分外显子发生拷贝数变异的检测结果。
检测可视化结果表明,本发明使用的基于LASSO回归模型检测单基因拷贝数变异方法能将测试的所有拷贝数变异阳性样本的拷贝数状态检测出来,且灵敏度高,特异性强。
采用本发明方法同样可以为多基因拷贝数变异检测,将单基因检测组合即为多基因检测。
以上仅是本发明的具体应用范例,对本发明的保护范围不构成任何限制;对于所属领域的普通技术人员来说,在上述说明的基础上还可以做出其它不同形式的变化或变动。这里无需也无法对所有的实施方式予以列举说明。凡采用等同变换或者等效替换而形成的类似此种的技术方案,均落在本发明权利保护范围之内。