一种基于枚举策略重建多倍体个体单体型的方法与流程

本发明属于生物信息学领域，具体涉及一种基于枚举策略重建多倍体个体单体型的方法。

技术背景

基因测序技术的快速发展产生了大量的基因数据，这使得对多倍体生物个体开展功能基因组研究工作成为可能。多倍体生物常常具有生长旺盛、存活率高、抗病力强、产量高、适应性和抗逆性强等优点，多倍体生物基因组的研究对探索多倍体物种的遗传特征意义重大，对提高多倍体物种的培育水平具有重要的理论价值和实用价值。单核苷酸多态性(Single Nucleotide Polymorphism，SNP)是基因组中最常见的一种遗传多态性，它是由染色体基因组在单个核苷酸碱基尺度上的变化而引起的DNA序列多态性。SNP的研究对探索多倍体生物的遗传特性和表型差异等方面发挥着重要的作用。由于连锁不平衡现象以及缺乏遗传重组事件，染色体上某一区域的一组SNP位点趋于共同遗传，这组位于一条染色体上的相关SNP位点序列即为单体型(Haplotype)，染色体上共同表现出来的SNP复合序列即组成了基因型(Genotype)。研究表明，单体型数据作为携带大量遗传疾病信息的载体，比单个SNP位点的研究和应用更深入和广泛，因此如何快速经济地获取个体单体型成为生物信息学关注的热点问题。由于在当前的实验技术水平下，直接通过生物学实验手段来测定一个个体的单体型的时间和金钱成本过高，因此利用计算技术来确定个体的单体型，即个体单体型重建问题，有极其重要的现实意义。个体单体型重建问题是指，给定来自同源染色体的由DNA测序方法测定的DNA片段数据，根据片段上的SNP位点取值，组装出同源染色体上的单体型。

个体单体型重建问题于2001年被首次提出，针对二倍体生物已经进行了深入的研究，并取得了不错的研究成果，对于多倍体单体型重建问题的研究较少。Li等对最少片段删除(Minimum Fragment Removal,MFR)模型进行研究，针对无空隙的SNP矩阵，提出时间复杂度为O(m²n+m^k+1)的多项式算法，其中m表示DNA片段个数，n表示SNP位点数，针对片段中洞的最大个数不超过t的SNP矩阵，提出时间复杂度为O(2^2tm²n+2^(k+1)tm^(k+1))的动态规划算法，该算法只能运用于m值(片段数目)较小的情况，当m值很大时，运行效率很低，不适用于实际应用中；Wang和Qian等人对最少错误更正(Minimum Error Correction,MEC)及带基因型数据的最少错误更正(MEC with Genotype Information,MEC/GI)模型进行研究，提出将其求解二倍体问题的遗传算法和粒子群算法进行扩展可以用来求解K-个体单体型重建问题，但这种扩展方法的性能受片段数目严重制约，且重建率不高；吴璟莉等提出求解MEC模型的三倍体个体单体型重建算法GTIHR，该方法通过引用新颖的染色体编码方式构造出了较小的解空间，但随着错误率和单体型长度的增加，重建率也不高，且随着单体型长度的增加，其效率不高；Emily等于2014年提出贝叶斯算法HapTree，该算法计算复杂且扩展性不好，在实际运用中受到一定限制；2015年，Sepideh等提出HapColor算法，该算法随着错误率的增加，MEC值呈明显上升趋势，重建效果下降明显。

技术实现要素：

本发明的目的是克服现有技术的不足而提供一种基于枚举策略重建多倍体个体单体型的方法，该方法需要的时间和金钱成本低，扩展性较好，能适用于多倍体个体单体型重建问题，重建率高,且其性能受片段数目、测序错误率和单体型长度等参数的影响较小。

本发明的方法是一种求解MEC/GI模型重建多倍体个体单体型的方法，下面给出一些符号定义和说明。

多倍体生物体细胞中含有k个染色体组(k为倍体数)，其染色单体上一段连续区域内相关的SNP序列就构成了一条单体型，即构成了k条单体型。SNP位点一般具有二态性，即每个SNP位点只有两种碱基组成，所以采用字符集{0,1}上的序列来表示单体型和基因型数据，而无需使用真正的碱基字符。如果染色体在某SNP位点的碱基取值相同，则这个位点上的基因型称为纯合子，否则称为杂合子。假设m条SNP片段是由k倍体中k条长度为n的单体型测序产生，由此定义SNP片段矩阵M_m×n，矩阵的行M_i-代表SNP片段，矩阵的列M_-j代表SNP位点，元素m_ij∈{0,1,-}(i＝1,2,…,m,j＝1,2,…,n)(-表示片段i未覆盖第j个SNP位点，或者片段i在第j个位点的取值未知)。令g＝g₁,g₂,...,g_n表示基因型矩阵，g_j＝(g_j1,g_j2,...,g_jk)^T表示第j个SNP位点的基因型，g_ji∈{0,1}(i＝1,...,k，j＝1,2,...,n)，k表示倍体数。

给定x，y∈{0,1,-}，定义

给定两个SNP片段M_u-＝(m_u1,m_u2,…,m_un)和M_v-＝(m_v1,m_v2,…,m_vn)，定义HD为片段M_u-和片段M_v-对应位取值相异的位数总和，即海明距离，如公式(2)所示：

HD可用于求解两条SNP片段之间的距离、SNP片段和单体型之间的距离以及两条单体型之间的距离。当HD表示两条SNP片段之间的距离时，若HD大于0，表示片段M_u-和M_v-冲突，否则表示它们相容。片段冲突意味着它们来自不同单体型或是片段数据中存在着测序错误。如果所有片段均没有测序错误或者错误已被更正，矩阵M中行可被划分成k(k与倍体数相对应)个彼此不相交且均包含相容片段的子集，每个子集通过片段重叠的部分推测出一条单体型，这时称片段矩阵M可行或无误差的。给定k条单体型h＝h₁,h₂,...,h_k(h_i＝h_i1,h_i2,...,h_in，i＝1,...,k)和基因型矩阵g＝g₁,g₂,...,g_n，g_j＝(g_j1,g_j2,...,g_jk)^T(g_ji∈{0,1}，i＝1,...,k，j＝1,2,...,n)，当时，称单体型h与基因型g相容。

根据上述定义，带基因型数据的最少错误更正(MEC with Genotype Information,MEC/GI)表述如下：给定SNP片段矩阵M和基因型矩阵g，修改矩阵M中最少数目的元素(0改为1或将1改为0),以使得修改过的SNP片段矩阵可行，且k个重建的单体型与基因型g相容。

实现本发明目的的技术方案是：

一种基于枚举策略重建多倍体个体单体型的方法，具体包括以下步骤：

1)预处理：对片段矩阵M和基因型矩阵g进行预处理，删掉对单体型重建工作无用的纯合SNP位点，预处理后只含杂合位点的片段矩阵和基因型矩阵分别用M'和g'表示；

2)进行枚举策略阶段：针对步骤1)得到的只含杂合位点的片段矩阵M'和基因型矩阵g'，按照SNP位点次序，依次根据每个杂合SNP位点的基因型枚举单体型在该位点上的所有可能取值，利用指标值计算函数计算SNP片段针对所有枚举值的指标值，并选取对应最大指标值的枚举值作为该位点的取值，以此重建出仅含杂合位点的单体型h'；

3)扩展阶段：将步骤1)中删除的纯合SNP位点重新加回，对于重建的仅含杂合位点的单体型h'，若某个已删除位点的取值记为0或1，则将单体型h'的相应SNP位点处插入0或1，从而扩展成最终的单体型h；

经过上述步骤，得到最终的单体型。

所述步骤1)中，删除对单体型重建工作无用的纯合SNP位点即删除基因型矩阵g中满足g_j1＝g_j2＝...＝g_jk(j＝1,2,...,n，k表示倍体数)的列，并将片段矩阵M中的对应列删除，记该列取值为g_j1，经过列删除后，矩阵M中会产生某些空行，即元素值全为-的行，将这些行也删除。预处理后得到只含杂合SNP位点的片段矩阵M'和基因型矩阵g'，类似的，M'中的行M'_i-表示SNP片段，列M'_-j表示SNP位点，元素m'_ij∈{0,1,-}，-表示片段i未覆盖第j个SNP位点，或者片段i在第j个位点的取值未知(i＝1,2,...,m'，j＝1,2,...,n')。基因型矩阵g'＝g'₁,g'₂,...,g'_n'，g'_j＝(g'_j1,g'_j2,...,g'_jk)^T,g'_ji∈{0,1}，k表示倍体数(i＝1,...,k，j＝1,2,...,n')。

所述步骤2)中，h'表示重建出的仅含杂合位点的单体型，即h'＝h'₁,h'₂,...,h'_k，k为倍体数，其中h'_i＝h'_i1,h'_i2,...,h'_in',h'_ij∈{0,1}(i＝1,...,k，j＝1,...,n')，指标值计算函数可以采用不同的距离函数，在实施枚举策略时，按照SNP位点次序重建仅含杂合位点的单体型。假定前j-1个位点已重建好，则根据第j个待重建SNP位点的基因型g'_j，枚举出单体型在第j个位点上所有可能的r个取值，并根据指标值计算函数得到r个对应的指标值E₁,…,E_r，令Em＝max{E₁,…,E_r}，则选取对应指标值为Em的枚举值作为第j个位点的单体型取值。

本发明的技术效果在于：本发明提出一种基于枚举策略重建多倍体个体单体型的方法，该方法需要的时间和金钱成本低，其性能受片段数目、测序错误率和单体型长度等参数的影响较小，能快速地重建出具有较高重建率的多倍体个体单体型。同时，在该枚举方法的框架下，指标值可以使用不同的距离函数来计算，算法灵活性强，具有较好的扩展性。

附图说明

图1：本实施例为三倍体单体型的示意图；

图2：本实施例为预处理后三倍体SNP片段矩阵M'的示意图；

图3：本实施例为预处理后三倍体基因型矩阵g'的示意图；

图4：一种基于枚举策略重建多倍体个体单体型的方法的流程图。

具体实施方式

下面结合附图和实施例对本发明内容作进一步的阐述，但不是对本发明的限定。

实施例：

一种基于枚举策略重建多倍体个体单体型的方法，具体包括以下步骤：

1)预处理：对片段矩阵M和基因型矩阵g进行预处理，删掉对单体型重建工作无用的纯合SNP位点。预处理后只含杂合位点的的片段矩阵和基因型矩阵分别用M'和g'表示；

2)进行枚举策略阶段：针对步骤1)得到的只含杂合位点的片段矩阵M'和基因型矩阵g'，按照SNP位点次序，依次根据每个杂合SNP位点的基因型枚举单体型在该位点上的所有可能取值，利用指标值计算函数计算SNP片段针对所有枚举值的指标值，并选取最大指标值对应的枚举值作为该位点的取值，以此重建出仅含杂合位点的单体型h'；

3)扩展阶段：将步骤1)中删除的纯合SNP位点重新加回，对于重建的仅含杂合位点的单体型h'，若某个已删除位点的取值记为0或1，则将单体型h'的相应SNP位点处插入0或1，从而扩展成最终的单体型h；

经过上述步骤，得到最终的单体型。

所述步骤1)中，删除对单体型重建工作无用的纯合SNP位点即删除基因型矩阵g中满足g_j1＝g_j2＝...＝g_jk(j＝1,2,...,n，k表示倍体数)的列，并将片段矩阵M中的对应列删除，记该列取值为g_j1。经过列删除后，矩阵M中会产生某些空行，即元素值全为-的行，将这些行也删除。预处理后得到只含杂合SNP位点的片段矩阵M'和基因型矩阵g'。类似的，M'中的行M'_i-表示SNP片段，列M'_-j表示SNP位点，元素m'_ij∈{0,1,-}，-表示片段i未覆盖第j个SNP位点，或者片段i在第j个位点的取值未知(i＝1,2,...,m'，j＝1,2,...,n')。基因型矩阵g'＝g'₁,g'₂,...,g'_n'，g'_j＝(g'_j1,g'_j2,...,g'_jk)^T,g'_ji∈{0,1}，k表示倍体数(i＝1,...,k，j＝1,2,...,n')。

所述步骤2)中，h'表示重建出的仅含杂合位点的单体型，即h'＝h₁',h'₂,...,h'_k，k为倍体数，其中h'_i＝h'_i1,h'_i2,...,h'_in',h'_ij∈{0,1}(i＝1,...,k，j＝1,...,n')。指标值计算函数可以采用不同的距离函数，这里介绍两种距离函数：

1)支持度函数

其中，

公式(3)可计算得出某个枚举值的支持度。在公式(3)中，m'表示片段矩阵M'的行，j(j＝2,...,n')表示待重建的位点，k表示倍体数，m'_ij表示片段矩阵M'的元素，h'_lj表示仅含杂合位点的单体型在第j个位点的取值，c(m'_it,m'_ij,h'_lt,h'_lj)表示m'_it和h'_lt及m'_ij和h'_lj两两比较都相等的数目。

2)差异度函数

其中，

公式(5)可计算得出某个枚举值的差异度。在公式(5)中，j(j＝2,...,n')表示待重建的位点，r(j)表示覆盖第j个位点的片段集，l(i)表示矩阵M'中第i行元素m'_ij≠-的最小列号，HD(h'_l,M'_i-,l(i),j)表示从l(i)位点到第j个位点，h'_l和第i行M'_i-对应位取值相异的位点数。

在实施枚举策略时，按照SNP位点次序重建仅含杂合位点的单体型。假定前j-1个位点已重建好，则根据第j个待重建SNP位点的基因型g'_j，枚举出单体型在第j个位点上所有可能的r个取值，并根据公式(3)(或公式(5))得到r个对应的支持度S₁,…,S_r(或差异度D₁,…,D_r)。令符号E₁,…,E_r分别表示r个枚举值对应的指标值，且E_i＝S_i(或E_i＝-D_i)(i＝1,…,r)。令Em＝max{E₁,…,E_r}，则选取对应指标值为Em的枚举值作为第j个位点的单体型取值。

如图1所示，是三倍体单体型为例的单体型示意图，图1中三倍体对应的三条单体型分别为CTC、CAT和ATC。

如图2所示，是预处理后三倍体SNP片段矩阵M'的示意图。

如图3所示，是预处理后三倍体基因型矩阵g'的示意图。

如图4所述，一种基于枚举策略重建多倍体个体单体型的方法，经过如下步骤得到最终的单体型：

执行S101：开始；

进一步，执行S102：输入片段矩阵M和基因型矩阵g；

进一步，执行S103：对S102输入的片段矩阵M和基因型矩阵g进行预处理，删掉对单体型重建工作无用的纯合SNP位点，得到只含杂合位点的片段矩阵M'和基因型矩阵g'；

进一步，执行S104：将g'₁的取值作为j＝1位点单体型的取值，初始化j＝2；

进一步，执行S105：根据g'_j产生r个枚举值；

进一步，执行S106：计算指标值E_i(i＝1,…,r)；

进一步，执行S107：选定E_m对应的枚举值，确定第j个位点单体型取值；

进一步，执行S108：j＝j+1；

进一步，执行S109：j≤n'；

若满足S109的条件，跳转执行S105，然后继续执行后续的步骤；

若不满足S109的条件，进一步，执行S110：进行扩展；

进一步，执行S111，输出最终的单体型；

进一步，执行S112，结束。

下面以三倍体个体为例，介绍其实施枚举策略的过程：

按照SNP位点序列，重建只含杂合SNP位点的三倍体个体单体型h'，假设h'的前3个SNP位点，即已经构建好，现在重建第4个SNP位点若第4个SNP位点的基因型g'₄取值为(001)^T,则枚举该杂合SNP位点的三种取值可能为若第4个SNP位点的基因型g'₄取值为(011)^T,则枚举该杂合SNP位点的三种取值可能为以附图3的基因型矩阵举例，第4个位点的基因型取值为(101)^T,依据附图2的SNP片段矩阵，利用支持度函数(或差异度函数)计算SNP片段对三种枚举值的支持度(或差异度)为S₁＝6，S₂＝7，S₃＝8(或D₁＝3，D₂＝2，D₃＝1)，则指标值为E₁＝6，E₂＝7，E₃＝8(或E₁＝-3，E₂＝-2，E₃＝-1)，选择Em＝E₃对应的第3个枚举值作为第4个位点的取值，则重建出的含前4个杂合SNP位点的单体型为依次类推，直至重建出全部杂合位点的取值，得到仅含杂合位点的单体型h'。