CN106611106B

CN106611106B - 基因变异检测方法及装置

Info

Publication number: CN106611106B
Application number: CN201611110748.2A
Authority: CN
Inventors: 何光铸; 王东辉; 蔡文君; 颜芹
Original assignee: UNITED ELECTRONICS CO Ltd
Current assignee: Ronglian Technology Group Co ltd
Priority date: 2016-12-06
Filing date: 2016-12-06
Publication date: 2019-05-03
Anticipated expiration: 2036-12-06
Also published as: CN106611106A

Abstract

本发明公开了一种基因变异检测方法及装置，包括：从基因比对结果中统计每个位点的比对信息；考虑碱基变异和插入缺失变异，创建16基因型模型；使用所述16基因型模型搜索候选变异位点；使用随机森林对候选变异位点进行分类与筛选，并输出筛选后的候选变异结果。本发明提供的基因变异检测方法及装置，能够同时检测单碱基变异和插入缺失变异，效率较高。

Description

基因变异检测方法及装置

技术领域

本发明涉及数据处理技术领域，特别是指一种基因变异检测方法及装置。

背景技术

基因组变异检测，这里指的是从二代测序数据的比对结果中，找出与参考基因组不同的碱基或序列片段，即单碱基变异(SNV)和插入缺失变异(INDEL)。

目前被广泛应用的10基因型模型只考虑了单碱基变异类型，插入缺失变异一般要单独检测，这使得现有模型的基因变异检测不够简便。

发明内容

有鉴于此，本发明的目的在于提出一种能够同时检测单碱基变异和插入缺失变异的基因变异检测方法及装置。

基于上述目的本发明提供的基因变异检测方法，包括：

从基因比对结果中统计每个位点的比对信息；

考虑碱基变异和插入缺失变异，创建16基因型模型；

使用所述16基因型模型搜索候选变异位点；

使用随机森林对候选变异位点进行分类与筛选，并输出筛选后的候选变异结果。

在一些可选实施方式中，所述从基因比对结果中统计每个位点的比对信息，具体包括以下比对信息：

碱基类型和对应的每个碱基类型的比对质量值、等位基因型及其Reads支持数量、正负链数量、插入缺失数量及插入序列信息，和/或，软剪切位点数量。

在一些可选实施方式中，所述考虑碱基变异和插入缺失变异，创建16基因型模型，具体包括：

假设样品是一个二倍体生物样品，碱基类型有ATCG四种，则二倍体基因型的统计类型有{AA,AC,AG,AT,CC,CG,CT,GG,GT,TT,AX,CX,GX,TX,XX,XY},其中X和Y分别代表有最多比对reads支持和第二多reads支持的插入或缺失。

在一些可选实施方式中，所述使用所述16基因型模型搜索候选变异位点，具体包括：

通过贝叶斯模型计算出每个位点最大可能的基因型；

将所述最大可能的基因型与参考基因组的对应位点的参考信息进行比较，得到所述候选变异位点。

在一些可选实施方式中，所述使用随机森林对候选变异位点进行分类与筛选，并输出筛选后的候选变异结果，具体包括：

定义真实变异位点和伪变异位点；

建立随机森林模型；

经过随机森林模型从所述候选变异位点中筛选得到更加可信的候选变异位点；

将所述更加可信的候选变异位点以VCF格式输出，并且直接应用于下游的分析工具。

本发明的另一方面，提供了一种基因变异检测装置，包括：

统计模块，用于从基因比对结果中统计每个位点的比对信息；

模型创建模块，用于考虑碱基变异和插入缺失变异，创建16基因型模型；

搜索模块，用于使用所述16基因型模型搜索候选变异位点；

分类与筛选模块，用于使用随机森林对候选变异位点进行分类与筛选，并输出筛选后的候选变异结果。

在一些可选实施方式中，所述模型创建模块，具体用于：

在一些可选实施方式中，所述搜索模块，具体用于：

通过贝叶斯模型计算出每个位点最大可能的基因型；

在一些可选实施方式中，所述分类与筛选模块，具体用于：

定义真实变异位点和伪变异位点；

建立随机森林模型；

从上面所述可以看出，本发明提供的基因变异检测方法及装置，通过考虑碱基变异和插入缺失变异，创建了16基因型模型，使得整体计算更加方便而且大幅提高了准确性和灵敏度；同时，利用随机森林对检测结果进行修正，使得检测结果更加精确。

附图说明

图1为本发明提供的基因变异检测方法的一个实施例的流程示意图；

图2为本发明提供的基因变异检测装置的一个实施例的模块结构示意图；

图3为采用本发明提供的基因变异检测方法及装置实施例后得到的预测正确率与真实比例的对比示意图。

具体实施方式

为使本发明的目的、技术方案和优点更加清楚明白，以下结合具体实施例，并参照附图，对本发明进一步详细说明。

需要说明的是，本发明实施例中所有使用“第一”和“第二”的表述均是为了区分两个相同名称非相同的实体或者非相同的参量，可见“第一”“第二”仅为了表述的方便，不应理解为对本发明实施例的限定，后续实施例对此不再一一说明。

基于上述目的，本发明实施例的第一个方面，提供了一种能够同时检测单碱基变异和插入缺失变异的基因变异检测方法的实施例。如图1所示，为本发明提供的基因变异检测方法的一个实施例的流程示意图。

所述基因变异检测方法，包括：

步骤101：从基因比对结果中统计每个位点的比对信息。这里，基因比对结果可以通过任意的基因比对软件的比对处理而得到，具体比对过程不再赘述。

在一些可选实施方式中，所述步骤101——从基因比对结果中统计每个位点的比对信息，具体包括以下比对信息：

具体地，测序数据经过质量值校正(score recalibration)、序列比对(alignment)、去重复(de-duplication)和重比对(realignment)等一系列处理后，需要收集每个位点的一套详细的统计信息以用于变异检测分析。每个位点的统计信息如下表1：

表1位点统计信息

对于等位基因型及其Reads支持数量(加权和)的统计信息：

在一个成功比对read中(Reads，读长，是高通量测序中获得的测序序列，每一个read是一段碱基序列)，每个碱基都会包含一个重校准的质量值，且质量值范围为0到40之间。为了储存碱基的质量值，我们为不同的质量值范围分配相应的权重，如下表2所示：

表2

碱基质量值	参数	权重
			0–10	[0–Weight0]	0
11–13	(Weight0–Weight1]	1
			14–17	(Weight1–Weight2]	2
18–20	(Weight2–Weight3]	3
			21–40	(Weight3–40]	4

表1中，为了将碱基质量值转化为权重值所设定的范围值参数，所述参数一列跟前面的碱基质量值一列的范围是相互对应的，这里的weight0\1\2\3,分别对为10、13、17、20。

每个成功配对的碱基使相应的等位基因型增加一个权重计数，如一个质量值为25的碱基A，其相应的等位基因型计数加4，若其质量值为5则计数加0。

对于正负链数量的统计信息：

依据比对结果，每个成功比对的碱基使对应的等位基因的正链或者负链计数加一。与权重计数不同，无论碱基的重校准质量值是多少，这里都增加一个计数。例如一个碱基的被多条碱基质量值小于10的reads覆盖，它的权重计数为零而正负链计数则确切地反映了成功比对的reads的条数。

对于插入缺失数量及插入序列信息的统计信息：

若比对结果中存在插入缺失，其信息将会被记录下来，格式为‘mI’或‘nD’，其中m和n分别表示插入和缺失的片段长度。除了不同类型的插入缺失的数量，插入的片段信息也会存储到一个动态分配的数据结构中且高质量与低质量片段信息分别记录在两个计数器中。

对于软剪切位点数量的统计信息：

如果比对结果中出现软剪切位点，其数量将会被同时记录下来。软剪切的方向也会被记录下来以区分头端剪切和末端剪切。

步骤102：考虑碱基变异和插入缺失变异，创建16基因型模型。

对于每个位点而言，我们需要依据S1中收集的比对信息来推测该位点的真实基因型并与参考基因组作比较，从而找出那些发生变异的位点，即候选变异位点。为了实现对一个位点真实基因型的推测，首先我们需要构建相应的基因型模型。

因此，在一些可选实施方式中，所述步骤102——考虑碱基变异和插入缺失变异，创建16基因型模型，具体包括以下步骤：

假设样品是一个二倍体生物样品，碱基类型有ATCG四种，则二倍体基因型的统计类型有{AA,AC,AG,AT,CC,CG,CT,GG,GT,TT,AX,CX,GX,TX,XX,XY},其中X和Y分别代表有最多比对reads支持和第二多reads支持的插入或缺失(reads支持越多，可信度越高)。

与被广泛应用的10基因模型不同，这里提出的16基因型模型在二倍体的背景中同时考虑了碱基变异和插入缺失变异，该16基因型模型统一了A,C,G,T和INDEL(插入缺失),这个统一的模型不仅使得计算方便而且大幅提高了准确性和灵敏度。

步骤103：使用所述16基因型模型搜索候选变异位点。

在一些可选实施方式中，所述步骤103——使用所述16基因型模型搜索候选变异位点，具体包括以下步骤：

通过贝叶斯模型计算出每个位点最大可能的基因型；

具体地，16基因型的后验概率的计算，使用了贝叶斯模型：

P(G|F)∝P(F|G)P(G)

其中，F表示观察到的{A,C,T,G,X,Y}各自的加权计数(weighted count)，P(G)表示某种基因型G的先验概率，P(F|G)表示的是基因型为G时观察到的F的概率，P(G|F)表示的是观察到F的基因型G的概率。

一般有如下几个原因导致我们观察到某个位置的碱基跟参考基因组上的不一样：

测序错误(bad base call or primary analysis)，比对错误(bad alignment)，基因变异(variant allele)。

一般质量值校正，可以一定程度修正第1类错误(即测序错误)。这里，我们设置两种错误概率：PS表示单碱基等位基因概率，PID表示插入缺失等位基因概率。一般经验，PS设置会大于PID。

如果一个错误(测序错误或比对错误)发生，假设：

1){A,C,G,T}每种碱基被观察到的概率相同，为PS；

2){X，Y}每个被观察到的概率相同，为PID。

定义错误率为：

P_err＝mP_s+nP_ID

其中，m为基因型G中的单个碱基{A,T,C,G}的数量，n为基因型G中{X,Y}的数量。

默认的设置：

P_S＝0.01

P_ID＝0.005

当我们观察到纯合的基因型时，我们会期望观察到接近100％的纯合位点。当观察到杂合的位点时，我们期望观察到50％的两个等位基因。为了检测观察到的reads覆盖深度分布与预期匹配的好坏，我们使用双尾费舍尔精确检验(Two-tailed Fisher’s ExactTest(FET))来检测，计算公式如下：

计算的p-value会当作某种基因型G的概率。[p-value越小表示可能性越大]。

具体计算P(F|G)的过程如下：

当观察到加权计数F＝{FA,FC,FG,FT,FX,FY}，

一个纯合基因型G＝AA的概率的计算，表示如下：

P(F|AA,P_err)＝P_hom(F_A)·P_e(F_C,F_G,F_T,F_X,F_Y)

一个杂合基因型G＝CG的概率计算，表示如下：

P(F|CG,P_err)＝P_het(F_C,F_G)·P_e(F_A,F_T,F_X,F_Y)

其中，P_hom为观察到纯合基因型的概率：

P_het为观察到杂合基因型的概率：

P_e为观察到基因型G以外的等位基因：

定义：

θ表示两个不相关单倍体单个碱基不同的频率，ω表示两个不相关单倍体单个插入却是不同的频率，ε表示转换颠换比(T_i/T_v)。

先验概率可表示如下表3：

表3

默认值：

θ＝0.001

ω＝0.0001

ε＝2.1

最终输出的基因型G_max，为有最大后验概率的基因型：

G_max＝argmax{P(G|F,P_err)}。

至此，我们通过贝叶斯模型计算出每个位点最大可能的基因型G_max，将这个基因型与参考基因组该位点的参考信息作比较，就能初步地得到我们想要的候选变异位点。而这些搜索到的候选变异位点还需要进一步的筛选，去除一些假阳性的变异位点，我们将在下一步使用随机森林的模型来实现。

步骤104：使用随机森林对候选变异位点进行分类与筛选，并输出筛选后的候选变异结果。

在一些可选实施方式中，所述步骤104——使用随机森林对候选变异位点进行分类与筛选，并输出筛选后的候选变异结果，具体包括以下步骤：

定义真实变异位点和伪变异位点；

建立随机森林模型；

具体地，变异分类的目的是为了给每一个检测出来的候选变异一个更加精确的预测正确率(Probability of a“true site”)，并基于这一正确率的估计值筛选出一个高准确率的变异位点的集合；这里的预测正确率，可参考表6中的预测正确率，是模型经过计算后给出预测正确的概率，模型的使用者依据这个概率来判断一个候选变异位点是否是真实的。随机森林是一种常用的机器学习的分类方法，我们的变异位点分类即是利用随机森林模型来对候选变异是真实的遗传变异(genetic variant)而非测序及分析导致的人为误差(artifact)的概率和变异指标之间的关系做一个连续的共变的估计，模型基于的分类依据如下：

1)真实变异位点(true sites)，一般来说这些位点在SNP(Single NucleotidePolymorphisms，单核苷酸的多态性)数据库(如dbSNP v129,HapMap 3,Omni2.5M SNP chiparray and Mills,1000G gold standard indels)中呈现多态性。

2)伪变异位点(false sites)，每个候选变异位点，若5个用于伪变异位点筛选的参数指标(Strand bias；Read position bias；Total depth；Left average basequality；Right average base quality)中有3个以上落在最差的5％内，则这个位点被归为伪变异位点。5％指的是这个候选变异在所有的由上一步贝叶斯模型检测出的候选变异位点中落在最差的5％。参照表5，对于链偏差，链偏差值取值范围(0,1],最差的5％则指的是链偏差值最小的5％的候选变异位点；对于Read位置偏差，位置偏差取值范围[-1,1],最差的5％为绝对值最大的5％；对于各等位基因深度总和，总测序深度越深越好，越少的read支持数，可信度越差，最差5％指的是深度最少的5％；对于位点左、右侧碱基平均质量值，碱基质量值取值范围为[0,40]，值越大越好，越小越差，也就越不可信，最差5％指的是质量值最小的5％。

之后，这一自适应性误差模型就可用于变异检测出的候选变异位点真实性的概率估算。

模型训练使用的特性如表4所示。

表4模型训练所用到的特性

用于挑选伪变异位点的特性如表5所示。

表5用于挑选伪变异位点的特性

模型训练详情及结果：

应用上文步骤102中介绍的16基因型模型从NA82178样本50×150bp的双端测序数据中搜索单碱基变异和插入缺失变异位点，再利用SNP数据库(dbSNP v137,IndelDB,1000Gand Mills)从这些候选变异位点中挑选真实变异位点。

这样，我们总共得到1,813,021个“true sites”和31,588个“false sites”。我们使用31,588个“false sites”和26,501个随机选取的“true sites”组成58,089个位点的训练集合。用这个训练集合建立了一个有96棵决策树的随机森林模型。

模型的可靠性分析如下表6所示：

表6

其中，Probability of a“true site”为随机森林模型给出的变异候选位点的预测正确率，预测正确率也就是模型经过计算后给出的预测得到的正确概率，即候选变异位点为真实变异位点(true site)的概率。模型的使用者依据这个概率来判断一个候选变异位点是否是真实可靠的。“比例”为训练集中“true sites”所占的真实比例，预测正确率与真实比例的对比如图3所示。从表6和图3可以看出我们的随机森林模型预测的正确率与“true sites”所占的真实的比例非常接近，可以说明我们的模型可以有效的区分候选变异位点是否为真实变异位点。

经过第三步的候选变异分类，我们进一步筛选出了更加可信的候选变异位点。最终的候选变异位点将以VCF(Variant Calling File)的格式输出，并且可以直接应用于下游的分析工具(如snpEff,VEP,GATK)和在线数据库(如Ingenuity,GenomeTrax)。

其中，所述输出结构中还可以包括每一个变异的质量值，每一个变异的质量值计算公式如下：

其中P_opt(G|F)是最大的后验概率，P_subOpt(G|F)是第二大后验概率。一般来说，质量值q越大，这一位点的最大概率基因型的不确定性越小，G_max也就越可信。

从上述实施例可以看出，本发明实施例提供的基因变异检测方法，通过考虑碱基变异和插入缺失变异，创建了16基因型模型，使得整体计算更加方便而且大幅提高了准确性和灵敏度；同时，利用随机森林对检测结果进行修正，使得检测结果更加精确。

本发明实施例的第二个方面，提供了一种基因变异检测装置的实施例。如图2所示，为本发明提供的基因变异检测装置的一个实施例的模块结构示意图。

所述基因变异检测装置，包括：

统计模块201，用于从基因比对结果中统计每个位点的比对信息；

模型创建模块202，用于考虑碱基变异和插入缺失变异，创建16基因型模型；

搜索模块203，用于使用所述16基因型模型搜索候选变异位点；

分类与筛选模块204，用于使用随机森林对候选变异位点进行分类与筛选，并输出筛选后的候选变异结果。

从上述实施例可以看出，本发明实施例提供的基因变异检测装置，通过考虑碱基变异和插入缺失变异，创建了16基因型模型，使得整体计算更加方便而且大幅提高了准确性和灵敏度；同时，利用随机森林对检测结果进行修正，使得检测结果更加精确。

在一些可选实施方式中，所述模型创建模块202，具体用于：

在一些可选实施方式中，所述搜索模块203，具体用于：

通过贝叶斯模型计算出每个位点最大可能的基因型；

在一些可选实施方式中，所述分类与筛选模块204，具体用于：

定义真实变异位点和伪变异位点；

建立随机森林模型；

需要特别指出的是，上述装置的实施例仅采用了所述方法的实施例来具体说明各模块的工作过程，本领域技术人员能够很容易想到，将这些模块应用到所述方法的其他实施例中。当然，由于所述方法实施例中的各个步骤均可以适当地进行相互交叉、替换、增加、删减，因此，这些合理的排列组合变换之于所述装置也应当属于本发明的保护范围，并且不应将本发明的保护范围局限在所述实施例之上。

所属领域的普通技术人员应当理解：以上任何实施例的讨论仅为示例性的，并非旨在暗示本公开的范围(包括权利要求)被限于这些例子；在本发明的思路下，以上实施例或者不同实施例中的技术特征之间也可以进行组合，步骤可以以任意顺序实现，并存在如上所述的本发明的不同方面的许多其它变化，为了简明它们没有在细节中提供。

另外，为简化说明和讨论，并且为了不会使本发明难以理解，在所提供的附图中可以示出或可以不示出与集成电路(IC)芯片和其它部件的公知的电源/接地连接。此外，可以以框图的形式示出装置，以便避免使本发明难以理解，并且这也考虑了以下事实，即关于这些框图装置的实施方式的细节是高度取决于将要实施本发明的平台的(即，这些细节应当完全处于本领域技术人员的理解范围内)。在阐述了具体细节(例如，电路)以描述本发明的示例性实施例的情况下，对本领域技术人员来说显而易见的是，可以在没有这些具体细节的情况下或者这些具体细节有变化的情况下实施本发明。因此，这些描述应被认为是说明性的而不是限制性的。

尽管已经结合了本发明的具体实施例对本发明进行了描述，但是根据前面的描述，这些实施例的很多替换、修改和变型对本领域普通技术人员来说将是显而易见的。例如，其它存储器架构(例如，动态RAM(DRAM))可以使用所讨论的实施例。

本发明的实施例旨在涵盖落入所附权利要求的宽泛范围之内的所有这样的替换、修改和变型。因此，凡在本发明的精神和原则之内，所做的任何省略、修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims

1.一种基因变异检测方法，其特征在于，包括：

从基因比对结果中统计每个位点的比对信息；

考虑碱基变异和插入缺失变异，创建16基因型模型，具体包括：样品是一个二倍体生物样品，碱基类型有ATCG四种，则二倍体基因型的统计类型有{AA,AC,AG,AT,CC,CG,CT,GG,GT,TT,AX,CX,GX,TX,XX,XY},其中X和Y分别代表有最多比对reads支持和第二多reads支持的插入或缺失；

使用所述16基因型模型搜索候选变异位点；

使用随机森林对候选变异位点进行分类与筛选，并输出筛选后的候选变异结果，具体包括：定义真实变异位点和伪变异位点；建立随机森林模型；经过随机森林模型从所述候选变异位点中筛选得到更加可信的候选变异位点；将所述更加可信的候选变异位点以VCF格式输出，并且直接应用于下游的分析工具。

2.根据权利要求1所述的方法，其特征在于，所述从基因比对结果中统计每个位点的比对信息，具体包括以下比对信息：

3.根据权利要求1所述的方法，其特征在于，所述使用所述16基因型模型搜索候选变异位点，具体包括：

通过贝叶斯模型计算出每个位点最大可能的基因型；

4.一种基因变异检测装置，其特征在于，包括：

模型创建模块，用于考虑碱基变异和插入缺失变异，创建16基因型模型；具体用于：样品是一个二倍体生物样品，碱基类型有ATCG四种，则二倍体基因型的统计类型有{AA,AC,AG,AT,CC,CG,CT,GG,GT,TT,AX,CX,GX,TX,XX,XY},其中X和Y分别代表有最多比对reads支持和第二多reads支持的插入或缺失；

搜索模块，用于使用所述16基因型模型搜索候选变异位点；

分类与筛选模块，用于使用随机森林对候选变异位点进行分类与筛选，并输出筛选后的候选变异结果；具体用于：定义真实变异位点和伪变异位点；建立随机森林模型；经过随机森林模型从所述候选变异位点中筛选得到更加可信的候选变异位点；将所述更加可信的候选变异位点以VCF格式输出，并且直接应用于下游的分析工具。

5.根据权利要求4所述的装置，其特征在于，所述从基因比对结果中统计每个位点的比对信息，具体包括以下比对信息：

6.根据权利要求4所述的装置，其特征在于，所述搜索模块，具体用于：

通过贝叶斯模型计算出每个位点最大可能的基因型；