CN109063417A - 一种构造隐马尔科夫链的基因型填补方法 - Google Patents

Abstract

本发明涉及一种构造隐马尔科夫链的基因型填补方法，利用样本具有的稀有位点上的少数等位基因型对基因型数据库进行预处理;创建局部单倍型簇模型；对局部单倍型簇模型的节点进行合并后续相似状态;利用向前‑向后算法对局部单倍型簇模型的通路进行剪枝;利用由蒙特卡洛模拟求得的单倍型概率，通过hardy‑weinberg平衡得到样本双倍型先念概率并结合测序似然值求得双倍型后验概率；由样本微单倍型筛选双倍型后念概率作为迭代重新构筑局部单倍型簇模型；循环N次，得到双倍体后验概率后并获得相应位点的基因型概率。本发明构建局部单倍型簇，从隐马尔科夫模型中产生双倍型的后验概率，按照样本的微单倍型将明显矛盾的滤去，重复迭代多次就能得到高准确性的判型。

Description

一种构造隐马尔科夫链的基因型填补方法

技术领域

本发明涉及一种构造隐马尔科夫链的基因型填补方法。

背景技术

基因是遗传的基本单元，携带有遗传信息的DNA（脱氧核糖核酸）或RNA（核糖核酸）序列，通过复制，把遗传信息传递给下一代，指导蛋白质的合成来表达自己所携带的遗传信息，从而控制生物个体的性状表达。基因检测是通过血液、其他体液、或细胞对DNA进行检测的技术，是取被检测者外周静脉血或其他组织细胞，扩增其基因信息后，通过特定设备对被检测者细胞中的DNA分子信息作检测，分析它所含有的基因类型和基因缺陷及其表达功能是否正常的一种方法。

在被测序的基因序列中，仍有部分位点因测序深度不够，位点覆盖度不高而没有被检测到。为了推测这些数据缺失的位点，引进了基因型填补方法。基因型填补是根据已经分型的位点的基因型对数据缺失的位点或者未分型的位点进行基因型预测的方法。基因型填充可以降低直接分型的成本，对不同基因型分型平台合并导致丢失的大量缺失的基因型填补，可助于对这些数据的联合分析。

由于HLA（细胞抗原）不同基因座位的某些等位基因经常连锁在一起遗传，而连锁的基因并非完全随机地组成单倍型，有些基因总是较多地一起出现，致使某些单倍型在群体中呈现较高的频率，从而引起连锁不平衡（LD）。而LD会随位点距离衰减，如果基因型数据在扩展区域定相时未考虑到这种定位，则将通过采样变异引入误差，导致远距离标记之间观察到明显的相关性，从而降低单倍型推断的准确性。基于位点间的不连锁平衡，目前的基因型填充主要分成两大类，一类是基于一般人群的GWAS（全基因组关联分析）中，利用非亲属个体的信息进行填补，另一类主要用于基于家族的GWAS,这两种方法的原理基本类似，由于采集样本有限，多数填充方法以第一类为基础。

目前市面上的基因型填充方法主要是利用参考群体提供的完整的基因型信息，构建出彼此连锁的标记的单倍型信息，然后利用单倍型信息，将目标群体基因型缺失位点的信息填充完整。根据是否利用家系信息，目前提出的主要填充方法包括两大类，一类是利用群体的连锁不平衡信息构建单倍型，相应的软件有：FAMHAP、fastPHASE、MaCH、Beagle、PLINK、BIMBAM、IMPUTE2、PHASEBOOK等。另一类是利用家系信息和标记连锁信息构建单倍型，相应的软件有Find-hap、Fimpute、AlphaImpute、PEDIMPUTE等。根据算法不同，又可将这些方法分为3类，即基于马尔可夫链（MCMC）算法、简化约束法（parsimony）以及期望最大化算法（EM算法）。

其中，Findhap软件是简化约束法的代表，它能够根据系谱的单倍型信息或者群体的单倍型信息进行缺失位点的填充。其具体做法为，首先将依据群体信息构建的单倍型根据单倍型频率进行排序，然后将每条染色体按每100个标记作为一个区段进行分割，对某个个体，将每一区段的每一个标记的基因型与排序好的单倍型上的基因型进行比对，如果比对结果每个标记皆为非相反纯合子，则判定这一位点比对上，而后将这一区段中的缺失位点根据比对的单倍型进行纯合子填充，然后再比对下一个区段。如果存在某一区段比对不上任何该区段的单倍型的情况，那么该区段形成一个新的单倍型，加到单倍型排序的末位，而该区段中的高密度位点则全部为未知基因型。进行完以上基于群体单倍型信息的填充过程后，再查系谱，将该个体的父亲和母亲的单倍型信息找到，来进一步校正该个体填充出的基因型。简化约束法运算速度较快，但在准确性相对其他算法而言不高。PLINK是EM算法的代表，这种算法

广泛运用于各个领域的各项研究中，并且在处理偏离哈代-温伯格平衡的群体时仍很稳健。其具体过程分为E步和M步，E步就是通过参考群体的高密度标记信息估计缺失位点的期望值，然后重复利用观测信息和估计信息进行下一轮的期望估计，直到M步，也就是估计的期望值达到稳定状态，这种算法最大的缺点就是准确性不高。

同时，随着高通量测序的应用逐渐普及，从测序数据中获取单倍型比利用SNP芯片中获取的基因型在某些程度上更适合应用于基因型填补算法。首先，二代测序分阶段多态性位点的密度比SNP芯片测试结果要高的多，并且包含了更多的低频位点。其次，当测序层数较低时，基因型只能被部分探测到，并且每种基因型都有一定程度的不确定性，这取决于覆盖该位点的层数。二代测序在每个位点用基因型可能性（GL）表示这种不确定性，而现有的基因型填补算法中，Thunder、Beagle、Impute2等都能将基因型可能性作为输入进行运算。

到目前为止，一些没有基于连锁不平衡的判型方法都是同利用二代测序中存在的单倍型数据。二代测序有着点位众多，而且很多都是稀有位点，在低倍情况下，每个位点的基因型不准确，只能用概率来描述。本发明充分利用二代低倍测序以上属性，充分挖掘二代测序数据的潜力，使得全基因组广度的填补成为可能。

发明内容

有鉴于此，本发明的目的在于提供一种基于稀有位点构造隐马尔科夫链的基因型填补方法，解决基因型填补问题，能大幅度提高填补的准确率，并在一定程度上提高运行时间。

为实现上述目的，本发明采用如下技术方案：

一种构造隐马尔科夫链的基因型填补方法，其特征在于：包括以下步骤：

步骤S1: 利用待测样本具有的在稀有位点上的少数等位基因型对基因型数据库进行预处理，得到预处理后的基因型数据库；

步骤S2:定义并创建局部单倍型簇模型；

步骤S3:对局部单倍型簇模型的节点进行合并后续相似状态;

步骤S4:利用单倍型簇模型的边统计算法对局部单倍型簇模型的通路进行剪枝;

步骤S5:利用由蒙特卡洛模拟求得局部单倍型簇模型的单倍型概率，结合Hardy-Weinberg平衡并与待测样本对应双倍型位点上的基因型概率求得双倍型后验概率；

步骤S6：由双倍型后验概率和待测样本微单倍型求得其组成单倍型的附加权重；

步骤S7:减少本轮单倍型权重并与对应附加权重相加求得下轮单倍型权重；

步骤S8:循环N次步骤S2-S7，其中N为预设值，且N≥10，N次循环以后，得到双倍型的后验概率并获得相应位点的基因型后验概率。

进一步的，所述步骤S1具体为：

步骤S11:利用待测样本位点信息对基因型数据库的单倍型进行局部判型;

步骤S12:利用待测样本中存在的稀有位点上的少数等位基因型对判型后的单倍型进行判断，将拥有该稀有位点的单倍型的权重*u，其中u为常数;

进一步的，所述步骤S2具体为：

步骤S21:假设单倍型数据库中有M个单倍型类型，每种单倍型类型在数据库中的数量为i_m，数据库中这些单倍型没有缺失的等位基因;

步骤S22:通过数据库中单倍型类型和数量，得到局部单倍型簇聚类模型。

进一步的，所述局部单倍型簇聚类模型是有向无环图，该图需要满足以下四个特征：

1）、该图有一个没有传入边的根节点和一个没有传出边的末节点；

2）、该图在M+1个级别中保持一致，每个节点都有一个级别m，

节点的等级为0，终端节点的等级为M；

3）、 对同一等级为m级的子节点的每条边用等m个等位基因标记，源于同一父节点的两条边不能用同一个等位基因标记；

4）、 对于数据库中的每个单倍型，存在从根节点到终端节点的路径，使得单倍型的第m个等位基因是路径的第m个边的标记；而单倍型簇的每条边e表示由所有从图的初始节点到终端节点的路径遍历e的单倍型组成的单倍型集合。

进一步的，所述步骤S3具体为；

步骤S31:对当前级别的每对节点计算得分，当得分少于m(n_x-1 +n_y-1)1/2+b时，两个节点才能合并；

其中，m和b为规模参数和转换参数；

步骤S32:将处理好的局部单倍型簇模型的每个簇定义为双列伪标记，然后对每个标记对应性状，进行卡方检验算出初始概率P值

步骤S33:利用置换检验来对多个卡方检验的结果进行调整。

进一步的，所述步骤S4具体为：

步骤S41:由后向前，按级别计算各节点的单倍型的数量；

步骤S42:对于同一节点出发的两条边，如果权重超过了500：1，权重小的边将被剪去;

步骤S43:移除剪去边所对应的单倍型，修正单倍型簇模型各边的权重；

进一步的，所述步骤S5具体为：

步骤S51：对剪枝后的局部单倍型簇模型进行蒙特卡洛模拟，具体为从起点开始随机挑选一条路径一直到终点，而从统计出在P次模拟后的不同通路次数的比率，即单倍型的概率；

步骤S52：将单倍型概率根据Hardy-Weinburg平衡做笛卡儿积求得双倍型概率；

步骤S53：将待测样本对应双倍型位点上的基因型概率作为发射函数，与双倍型概率相乘求得双倍型的比率，并由贝叶斯公式求得双倍型的后验概率；

进一步的，所述步骤S6具体为：

步骤S61：将求得的双倍型的后验概率作为其组成单倍型的部分权值，总和所有同一单倍型的部分权值作为该单倍型的附加权值；

步骤S62：将单倍型与样本二代测序所得的微单倍型对比，若与其相悖则将单倍型的附加权值除以5；

进一步的，所述步骤S7具体为：

步骤S71：将本轮全部单倍型的权重除以2；

步骤S72：将单倍型权重与对应的附加权重相加作为下一轮的单倍型权重；

进一步的，所述步骤S8具体为：

步骤S81：若当前迭代次数小于N次，返回S2。

步骤S82：根据各双倍型的后验概率可以得到各基因位点基因型的后验概率，给出基因型后出错概率作为输出。

本发明与现有技术相比具有以下有益效果：

1、本发明使用百万级大数据的样本池作为隐马尔科夫模型的观测状态，增加算法准确度；

2、本发明将隐马尔可夫模型计算的节点设定为在目标数据中进行基因分型的标记簇，减少了存储器需求和计算时间；

3、本发明在填补过程前，利用稀有位点（MAF<1%）作为seed使数据库由百万级降为数百，从而在不损伤精度的情况下有效成倍地减少算法运行时间；

4、本发明在每次迭代过程采样时使用二代测序的部分单倍型数据对单倍型进行筛选，从而提高由基因型数据判型的单倍型的准确性，也进一步减少运算时间。

附图说明

图1是本发明流程图

图2是本发明一实施例数据图

图3是本发明一实施例原理图

图4是本发明另一实施例原理图。

具体实施方式

下面结合附图及实施例对本发明做进一步说明。

请参照图1，本发明提供一种构造隐马尔科夫链的基因型填补方法，其特征在于：包括以下步骤：

步骤S1: 利用待测样本具有的在稀有位点上的少数等位基因型对基因型数据库进行预处理，得到预处理后的基因型数据库；

步骤S2:定义并创建局部单倍型簇模型；

步骤S3:对局部单倍型簇模型的节点进行合并后续相似状态;

步骤S4:利用单倍型簇模型的边统计算法对局部单倍型簇模型的通路进行剪枝;

步骤S5:利用由蒙特卡洛模拟求得局部单倍型簇模型的单倍型概率，结合Hardy-Weinberg平衡并与待测样本对应双倍型位点上的基因型概率求得双倍型后验概率；

步骤S6：由双倍型后验概率和待测样本微单倍型求得其组成单倍型的附加权重；

步骤S7:减少本轮单倍型权重并与对应附加权重相加求得下轮单倍型权重；

步骤S8:循环N次步骤S2-S7，其中N为预设值，且N≥10，N次循环以后，得到双倍型的后验概率并获得相应位点的基因型后验概率。

在本发明一实施例中，所述步骤S1具体为：

步骤S11:利用待测样本位点信息对基因型数据库的单倍型进行局部判型;

步骤S12:利用待测样本中存在的稀有位点上的少数等位基因型对判型后的单倍型进行判断，将拥有该稀有位点的单倍型的权重*u，其中u为常数;

在本发明一实施例中， 进一步的，所述步骤S2具体为：

步骤S21:假设单倍型数据库中有M个单倍型类型，每种单倍型类型在数据库中的数量为i_m，数据库中这些单倍型没有缺失的等位基因;

步骤S22:通过数据库中单倍型类型和数量，得到局部单倍型簇聚类模型。

在本发明一实施例中，进一步的，所述局部单倍型簇聚类模型是有向无环图，该图需要满足以下四个特征：

1）、该图有一个没有传入边的根节点和一个没有传出边的末节点；

2）、该图在M+1个级别中保持一致，每个节点都有一个级别m，

节点的等级为0，终端节点的等级为M；

3）、 对同一等级为m级的子节点的每条边用等m个等位基因标记，源于同一父节点的两条边不能用同一个等位基因标记；

4）、 对于数据库中的每个单倍型，存在从根节点到终端节点的路径，使得单倍型的第m个等位基因是路径的第m个边的标记；而单倍型簇的每条边e表示由所有从图的初始节点到终端节点的路径遍历e的单倍型组成的单倍型集合。

在本发明一实施例中，进一步的，所述步骤S3具体为；

步骤S31:对当前级别的每对节点计算得分，当得分少于m(n_x-1 +n_y-1)1/2+b时，两个节点才能合并；

其中，m和b为规模参数和转换参数；

步骤S32:将处理好的局部单倍型簇模型的每个簇定义为双列伪标记，然后对每个标记对应性状，进行卡方检验算出初始概率P值

步骤S33:利用置换检验来对多个卡方检验的结果进行调整。

在本发明一实施例中，进一步的，所述步骤S4具体为：

步骤S41:由后向前，按级别计算各节点的单倍型的数量；

步骤S42:对于同一节点出发的两条边，如果权重超过了500：1，权重小的边将被剪去;

步骤S43:移除剪去边所对应的单倍型，修正单倍型簇模型各边的权重；

在本发明一实施例中，进一步的，所述步骤S5具体为：

步骤S51：对剪枝后的局部单倍型簇模型进行蒙特卡洛模拟，具体为从起点开始随机挑选一条路径一直到终点，而从统计出在P次模拟后的不同通路次数的比率，即单倍型的概率；

步骤S52：将单倍型概率根据Hardy-Weinburg平衡做笛卡儿积求得双倍型概率；

步骤S53：将待测样本对应双倍型位点上的基因型概率作为发射函数，与双倍型概率相乘求得双倍型的比率，并由贝叶斯公式求得双倍型的后验概率；

进一步的，所述步骤S6具体为：

步骤S61：将求得的双倍型的后验概率作为其组成单倍型的部分权值，总和所有同一单倍型的部分权值作为该单倍型的附加权值；

步骤S62：将单倍型与样本二代测序所得的微单倍型对比，若与其相悖则将单倍型的附加权值除以5；

在本发明一实施例中，进一步的，所述步骤S7具体为：

步骤S71：将本轮全部单倍型的权重除以2；

步骤S72：将单倍型权重与对应的附加权重相加作为下一轮的单倍型权重；

进一步的，所述步骤S8具体为：

步骤S81：若当前迭代次数小于N次，返回S2。

步骤S82：根据各双倍型的后验概率可以得到各基因位点基因型的后验概率，给出基因型后出错概率作为输出。

为了让一般技术人员更好的理解本发明的技术方案，以下结合附图对本发明进行详细介绍。

Beagle 4.1版本可以将数百万单倍型作为参考组，利用参考单倍型的数量弥补算法过程中的近似。但是根据基因概率型算法的模块（应用于二代低倍测序）还是只能使用双倍型状态，所以该模块的计算时间正比于参考组单倍型数量的平方，计算时间上无法满足全基因组广度的填补的任务。本算法和Beagle相比，主要有以下几方面的改进，从而使全基因组广度的填补成为可能：

1）. 隐马尔科夫链的状态一直使用单倍型，并用蒙特卡洛方法确定单倍型的人群概率，在最后计算发射函数时，才确定样本双倍型的概率是单倍型概率的笛卡尔积，从而大大减少了状态和运行时间。

2). 在初始化人群单倍型簇模型时，通过比较目标样本在稀有位点上的基因型增加和目标样本相符合的单倍型的权重，从而减少了无效的探索，快速迭代优化。

3). 利用二代测序的读数覆盖多个位点形成了微单倍型，有效排除了错误的单倍型赋值。

在本发明一实施例中，将用本算法与其他现有的Beagle算法和基于单窗口的GeneImp算法进行比较，其中参考单倍型为国际千人基因组第三版，一共26个民族2504个样本，只填补常染色体部分：

将R^2的平均值作为MAF的函数值，每个函数值都是由在MAF库内的每个最小基因型概率的R^2的值求平均而得， 我们的方法所得值与BEAGLE相差无几，而基于单窗口的GeneImp的R^2的平均值则与另外两种相差较大。

将三种方法的R2平均值、运行时间以及占用内存，具体数据见表1。由表1可知， 基于单窗口的GeneImp的R2与另外两种方法相差较大。而我们的算法在确保R2平均值相差不大的情况下，大幅度的缩短了运行时间，同时也将运行时占用的CPU缩小为BEAGLE运行时CPU的一半。

表1 三种算法的参数比较

算法	R<sup>2</sup>平均值	运行时间	占用CPU
				本算法	0.933	8.9	16个并行进程
BEAGLE	0.939	136.7	32个并行进程
				单窗口的GeneImp	0.901	1.6	16个并行进程

图2-4给出如何优化单倍簇群。假设一共4个基因位点，1表示和参考基因组相同，2表示和参考基因组不同，人群单倍型频率统计如图2所示，直接构造的单倍簇群如图3所示，这儿用实线表示1，虚线表示2，结合表格计算每条边的单倍型总数，可以得到n(e_A’)=311, n(e_B’)=289, n(e_c’)=195，n(e_D’)=116, n(e_E’)=289，n(e_F’)=100, n(e_G’)=95, n(e_H’)=116, n(e_I’)=137, n(e_J;)=152, n(e_K’)=21, n(e_L)=79，n(e_m’)=95, n(e_n’)=116, n(e_o’)=25, n(e_p’)=112, n(e_q’)=152。

然后从终结点分层向后，合并节点，最后一层可以合并（从上往下数）1节点和4节点，2节点和5节点。最后第二层可以合并1节点和3节点。合并的节点大致分叉权重比相当，合并以后如图4所示，这样就减少了节点数量。

比较图3和图4不难得出，第三层第1第3节点合并，对应于该人群的染色体重组。

以上所述仅为本发明的较佳实施例，凡依本发明申请专利范围所做的均等变化与修饰，皆应属本发明的涵盖范围。

Claims (10)

1.一种构造隐马尔科夫链的基因型填补方法，其特征在于：包括以下步骤：

步骤S1: 利用待测样本具有的在稀有位点上的少数等位基因型对基因型数据库进行预处理，得到预处理后的基因型数据库；

步骤S2:定义并创建局部单倍型簇模型；

步骤S3:对局部单倍型簇模型的节点进行合并后续相似状态;

步骤S4:利用单倍型簇模型的边统计算法对局部单倍型簇模型的通路进行剪枝;

步骤S5:利用由蒙特卡洛模拟求得局部单倍型簇模型的单倍型概率，结合Hardy-Weinberg平衡并与待测样本对应双倍型位点上的基因型概率求得双倍型后验概率；

步骤S6：由双倍型后验概率和待测样本微单倍型求得其组成单倍型的附加权重；

步骤S7:减少本轮单倍型权重并与对应附加权重相加求得下轮单倍型权重；

步骤S8:循环N次步骤S2-S7，其中N为预设值，且N≥10，N次循环以后，得到双倍型的后验概率并获得相应位点的基因型后验概率。

2.根据权利要求1所述的一种构造隐马尔科夫链的基因型填补方法，其特征在于：所述步骤S1具体为：

步骤S11:利用待测样本位点信息对基因型数据库的单倍型进行局部判型;

步骤S12:利用待测样本中存在的稀有位点上的少数等位基因型对判型后的单倍型进行判断，将拥有该稀有位点的单倍型的权重*u，其中u为常数。

3.根据权利要求1所述的一种构造隐马尔科夫链的基因型填补方法，其特征在于：所述步骤S2具体为：

步骤S21:假设单倍型数据库中有M个单倍型类型，每种单倍型类型在数据库中的数量为i_m，数据库中这些单倍型没有缺失的等位基因;

步骤S22:通过数据库中单倍型类型和数量，得到局部单倍型簇聚类模型。

4.根据权利3所述的一种构造隐马尔科夫链的基因型填补方法，其特征在于：所述局部单倍型簇聚类模型是有向无环图，该图需要满足以下四个特征：

1）、该图有一个没有传入边的根节点和一个没有传出边的末节点；

2）、该图在M+1个级别中保持一致，每个节点都有一个级别m，

节点的等级为0，终端节点的等级为M；

3）、 对同一等级为m级的子节点的每条边用等m个等位基因标记，源于同一父节点的两条边不能用同一个等位基因标记；

4）、 对于数据库中的每个单倍型，存在从根节点到终端节点的路径，使得单倍型的第m个等位基因是路径的第m个边的标记；而单倍型簇的每条边e表示由所有从图的初始节点到终端节点的路径遍历e的单倍型组成的单倍型集合。

5.根据权利要求1所述的一种基于稀有位点构造隐马尔科夫链的基因型填补方法，其特征在于：所述步骤S3具体为；

步骤S31:对当前级别的每对节点计算得分，当得分少于m(n_x-1 +n_y-1)1/2+b时，两个节点才能合并；

其中，m和b为规模参数和转换参数；

步骤S32:将处理好的局部单倍型簇模型的每个簇定义为双列伪标记，然后对每个标记对应性状，进行卡方检验算出初始概率P值

步骤S33:利用置换检验来对多个卡方检验的结果进行调整。

6.根据权利要求7所述的一种构造隐马尔科夫链的基因型填补方法，其特征在于：所述步骤S4具体为：

步骤S41:由后向前，按级别计算各节点的单倍型的数量；

步骤S42:对于同一节点出发的两条边，如果权重超过了500：1，权重小的边将被剪去;

步骤S43:移除剪去边所对应的单倍型，修正单倍型簇模型各边的权重。

7.根据权力要求1所述的一种构造隐马尔科夫链的基因型填补方法，其特征在于：所述步骤S5具体为：

步骤S51：对剪枝后的局部单倍型簇模型进行蒙特卡洛模拟，具体为从起点开始随机挑选一条路径一直到终点，而从统计出在P次模拟后的不同通路次数的比率，即单倍型的概率；

步骤S52：将单倍型概率根据Hardy-Weinburg平衡做笛卡儿积求得双倍型概率；

步骤S53：将待测样本对应双倍型位点上的基因型概率作为发射函数，与双倍型概率相乘求得双倍型的比率，并由贝叶斯公式求得双倍型的后验概率。

8.根据权力要求1所述的一种构造隐马尔科夫链的基因型填补方法，其特征在于：所述步骤S6具体为：

步骤S61：将求得的双倍型的后验概率作为其组成单倍型的部分权值，总和所有同一单倍型的部分权值作为该单倍型的附加权值；

步骤S62：将单倍型与样本二代测序所得的微单倍型对比，若与其相悖则将单倍型的附加权值除以5。

9.根据权力要求所述的一种构造隐马尔科夫链的基因型填补方法，其特征在于：所述步骤S7具体为：

步骤S71：将本轮全部单倍型的权重除以2；

步骤S72：将单倍型权重与对应的附加权重相加作为下一轮的单倍型权重。

10. 根据权力要求 所述的一种构造隐马尔科夫链的基因型填补方法，其特征在于：所述步骤S8具体为：

步骤S81：若当前迭代次数小于N次，返回S2；

步骤S82：根据各双倍型的后验概率可以得到各基因位点基因型的后验概率，给出基因型后出错概率作为输出。