CN101256602A - 基于优化解集合的个体单体型重建方法 - Google Patents

Abstract

本发明公开了一种基于优化解集合的个体单体型重建方法，包括以下步骤：预处理单核苷酸多态性(single nucleotide polymorphism，SNP)矩阵，去掉对重建工作没有帮助的冗余信息，得到只含杂合位点的SNP矩阵。通过粒子群优化策略得到一个小规模的优化解集合，即只带有杂合位点的单体型对集合。最后的扩展阶段将预处理阶段删掉的SNPs重新加上，得到最终的单体型对集合。本发明提出基于小规模优化解集合求解MEC模型的个体单体型重建方法，该方法可以获得比以往相关方法更高的单体型重建率，并且在求解大规模问题时仍具有较高的执行效率。

Description

基于优化解集合的个体单体型重建方法

技术领域

本发明涉及生物信息学，特别涉及个体单体型的重建。

背景技术

人类基因组测序工作完成之后，遗传差异性研究已成为基因组的热点研究之一。众所周知，人类几乎有99.9％的基因是相同的，因此我们所呈现出的外部差异性仅仅是由于0.1％的基因差异引起的。在各种遗传变异之中，单核苷酸多态性(single nucleotide polymorphisms，SNPs)是最显著的一种形式，它是人类染色体某个位点上的碱基变化。研究SNP在阐明疾病易感性机制、设计个体化治疗方案和药物研制等方面都具有重要意义和实际应用价值。

然而，检测人类染色体上所有一千万个常见SNPs的费用极其昂贵，所幸的是，由于连锁不平衡现象以及缺乏重组事件，一些相邻的多态位点趋于在一起共同遗传，这些变异连锁的区域即为单体型(haplotype)，它定义为一条染色单体上某一区域的一组相关联的SNP位点。最近的研究表明，在与疾病相关的研究中，单体型数据通常比单个SNP携带更多的信息，但在当前的实验技术下，直接通过生物学实验手段来测定单体型既费钱又费时间，因此利用计算机技术来确定个体的单体型有极其重要的现实意义。

个体单体型重建问题可描述如下：给定一组来自某对同源染色体的由DNA测序方法得到的DNA片断，若只关注SNP位点，这些DNA片断即为SNP片断。单体型重建问题是要根据片断上SNP位点的状态信息将这些片断分成两个集合，每个集合中的片断组装成一条单体型。由于在DNA测序过程中会产生测序错误，而且当片断中存在错误时，无法准确地对片断进行分组。因此，在2002年，Lippert等提出了最少错误更正(the minimum error correction，MEC)模型，它要求通过更正最少的片断错误来重建单体型，目前求解该模型的方法主要有：

(1)王瑞省等提出的基于分支定界思想的方法，但由于MEC模型是NP难的，该方法无法求解大规模问题。

(2)王瑞省等提出两种动态聚类方法(文中称为DC1和DC2)以及一种基于遗传算法的启发式方法(文中称为GA)。

这些方法均致力于得到一对重建率最高的单体型数据。但是由于MEC模型及方法本身的原因，最优结果会在问题求解过程中被遗失，从而使结果单体型重建率并不高。

发明内容

为了解决上述基于MEC模型的个体单体型重建方法存在的技术问题，本发明提供了一种基于优化解集合求解MEC模型的重建方法。该方法能够生成一个小规模的优化解集合，且基于该优化解集，能够获得较以往方法更高重建率的单体型。

本发明基于MEC模型解决个体单体型重建问题，包括以下步骤：预处理SNP矩阵，得到只含杂合位点的SNP矩阵。通过粒子群优化策略得到一个小规模的优化解集合，即只具有杂合位点的单体型对集合。最后的扩展阶段将预处理阶段删掉的SNPs重新加上，得到最终的单体型对集合。

上述的基于优化解集合的个体单体型重建方法，粒子群优化策略采用二进制串X(x₁，x₂，...，x_n)(x_i∈{0，1})和V(v₁，v₂，...，v_n)(v_i∈{0，1})来分别表示一个粒子的位置和速度，粒子位置代表一条只含杂合位点的单体型。

上述的基于优化解集合的个体单体型重建方法，某个粒子位置X对应的错误更正数E(X)的计算方式如下：

$E\left(X\right)=\underset{i=1}{\overset{m1}{\mathrm{\Σ}}}\min (S(f_i,X),D(f_i,X))$

上式中m1表示预处理后SNP矩阵的行数，f_i表示SNP矩阵中的行(SNP片断)，S(f_i，X)表示f_i和X对应位取值相同(同为1或同为0)的位数，D(f_i，X)表示f_i和X对应位取值相异(一个为1则另一个为0)的位数。

本发明的技术效果在于：本发明提出一种生成小规模优化解集合以降低最优解遗失概率的新研究思路。基于这种研究思路，针对SNP位点杂合率较低的特点，设计了一种短的粒子编码，给出求解MEC模型的粒子群优化方法。这种短粒子编码一方面能够有效控制解空间的大小，使本发明更容易获得最优解；另一方面它使粒子群优化策略能够使用小的群体规模，继而可以生成一个小规模的优化解集合，解集合中的解均与粒子群中的P_g有相同的适应值，但它们的重建率不一定相同。优化解集合的研究思路与短粒子编码方式的结合使用，使得本发明能够获得比以往相关方法更高的单体型重建率。另外，由于群体规模得到控制，本发明即使在求解大规模问题时，仍具有较高的执行效率，因此具有很高的实用价值。

附图说明

图1：本发明的流程图。

具体实施方式

下面结合附图对本发明的具体实施作进一步说明。参见图1，图1为本发明的流程图，虚线框部分表示粒子群优化方法。本发明中预处理SNP矩阵M_m×n，去掉对重建工作没有帮助的冗余信息，即删除M中所有满足条件f₀≤t或f₁≤t的列(在矩阵M中，令n_x为某一列中值为x的元素个数，且f_x＝n_x/(n_x+n_1-x))，这里t设置为0.2，若被删除的列中大部分非空元素值为0，则称其为0-列，否则称为1-列。将所有满足上述条件的列删除之后，某些行将变成空行(元素值全为-)，它们对于重建工作没有任何帮助，因此也将其删除。预处理后得到只含杂合位点的SNP矩阵M1_m1×n1。实施粒子群优化策略，将与P_g适应值相同的k个解均保留下来，并将k个解转换成只含杂合位点的单体型对H_i’＝(h_i1’，h_i2’)(i＝1，...，k)，以得到一个规模为k的优化解集合H’＝{H₁’，...，H_k’}。最后的扩展阶段将预处理阶段删掉的SNPs重新加上，对于优化解集H’中只含杂合位点的单体型对H_i’＝(h_i1’，h_i2’)，如果某个已被删除的同合位点为0-列(1-列)，则将0(1)插回到单体型对(h_i1’，h_i2’)的相应位置，以此得到扩展后的单体型对H_i＝(h_i1，h_i2)，扩展结束后得到最终的单体型对集合H＝{H₁，...，H_k}。

基于这种生成小规模优化解集的研究思路，本发明设计了一种短的粒子编码，从而提出求解MEC模型的粒子群优化方法。这种短的粒子编码采用二进制串X(x₁，x₂，...，x_n1)(x_i∈{0，1})来表示一个粒子的位置，它代表一条只含杂合位点的单体型。如前所述，新矩阵M1中的片断只保留了杂合位点，则由它们构建的单体型必定也只具有杂合位点。由于一对单体型在其杂合位点上的值是相异的(一个值为0(1)，另一个值则为1(0))，因此对于这样一对只具有杂合位点的单体型，可以通过其中一条推导出另一条。所以，由一个粒子的位置可以推导出一对只含杂合位点的单体型。

将粒子群优化方法运用于离散问题时，需要对粒子的速度表示及粒子间的运算操作进行定义：

(a)粒子的速度V定义为其两次位置X₁和X₂之间的距离。

V＝X₁-X₂＝(v₁，...，v_n)，

$v_i=\left\{\begin{array}{c}0,x_{1i}=x_{2i}\\ 1,x_{1i}\≠x_{2i}\end{array}\right.i=\mathrm{1,2},...,n$

(b)速度V₁和V₂间的加法操作定义为其相应位的逻辑加，结果为速度V。

V＝V₁+V₂＝(v₁，...，v_n)，v_i＝v_1iORv_2i，i＝1，2，...，n}.

(c)粒子速度V₁与概率C的乘积，结果为速度V。

V＝CV₁＝(v₁，...，v_n)，

(d)速度V和位置X₁间的加法操作定义为其相应位的逻辑异或，结果为位置X。

X＝X₁+V＝(x₁，...，x_n)，x_i＝x_1iXORv_i，i＝1，2，...，n}.

适应度函数用于评价粒子的搜索性能，指导粒子群的搜索过程。给定某个粒子位置X及矩阵M1中的所有片断f_i(i＝1，...，m1)，X的适应度函数Fitness(X)定义为：

$\mathrm{Fitness}\left(X\right)=1-\frac{E\left(X\right)}{m1\×n1},$

$E\left(X\right)=\underset{i=1}{\overset{m1}{\mathrm{\Σ}}}\min (S(f_i,X),D(f_i,X))$

其中，粒子位置X表示一对仅含杂合位点的单体型(h₁’，h₂’)中的一条，例如h₁’。于是S(f_i，X)表示片断f_i与单体型h₁’间等位基因相同的位点个数，即片断f_i与单体型h₂’间等位基因相异的位点个数；D(f_i，X)表示片断f_i与单体型h₁’间等位基因相异的位点个数；E(X)表示对应于单体型对(h₁’，h₂’)的最少错误更正数。

综上所述，优化解集合的研究思路与短粒子编码方式的结合使用可以有效降低最优解的遗失概率，从而获得具有更高重建率的单体型。

利用计算机模拟真实生物数据的特征生成测试数据集进行实验测试。实验在一台安装了Windows XP Professional操作系统的IBM工作站(Intel PentiumIV 2.0GHz，内存为512MB)上进行，程序编译器为Microsoft Visual C++6.0。

本发明中用“重建率”和“运行时间”来测试本发明方法的性能。在本发明的优化解集合H中选取重建率最大的单体型对作为结果，即该方法结果单体型对的重建率为max{RR(H₁)，...，RR(H_k)}(H₁，...，H_k∈H)。表1到表5中的每个计算结果均为100次重复测试的平均值。在下面的实验中，本发明的参数设置如下：w＝0.8，C₁＝C₂＝0.7，群体规模N为20，迭代次数M-ITER为100。

表1至表3的实验结果显示本发明能够获得较现有方法更高重建率的单体型对，这说明优化解集的研究思路能有效避免最优解的遗失。表中k表示优化解集合H中平均单体型的对数，其平均值均不超过4对，满足了解集合规模不宜太大的要求。本发明引入的粒子编码较短，使得本发明能够采用小规模种群，这为生成小规模的优化解集合奠定了基础。

表1重建率的比较(c＝5，n＝100)

表2重建率的比较(n＝100，P_s＝0.05)

表3重建率的比较(c＝5，P_s＝0.05)

表4运行时间比较(n＝100，P_s＝0.05)

表5运行时间比较(c＝5，P_s＝0.05)

表4和表5的结果显示，GA方法的运行时间最长，方法DC1和DC2运行速度较快，最长时间不超过0.1秒，虽然本发明比这两种动态聚类方法运行得慢，但最多也只需花费几秒钟，因此具有很高的实用价值。

从以上实验数据看来，应用本发明方法的得到的重建率和花费的运行时间均比较理想，这是因为本发明成功地将优化解集合的研究思路与粒子群优化方法相结合，通过引入短的粒子编码方式，使粒子群优化策略能够使用小的群体规模，这使得通过生成优化解集合来降低最优解的地遗失概率成为可能，且由于群体规模得到控制，本发明即使在求解大规模问题时，仍具有较高的执行效率，因此具有很高的实用价值。另外短粒子编码还能够有效控制解空间的大小，使本发明更容易获得最优解。

Claims (1)

1.一种基于优化解集合的个体单体型重建方法，包括以下步骤：

(1)预处理SNP矩阵M_m×n，去掉对重建工作没有帮助的冗余信息，即删除M中所有满足条件.f₀≤t或.f₁≤t的列，在矩阵M中，令n_x为某一列中值为x的元素个数，且.f_x＝n_x/(n_x+n_1-x)，t设置为0.2，若被删除的列中大部分非空元素值为0，则称其为0-列，否则称为1-列，将所有满足上述条件的列删除之后，得到只含杂合位点的SNP矩阵M1_m1×n1；

(2)通过粒子群优化策略得到一个小规模的优化解集合，即只具有杂合位点的单体型对集合，将与P_g适应值相同的k个解均保留下来，并将k个解转换成只含杂合位点的单体型对H_i’＝(h_i1’，h_i2)(i＝1，...，k)，以得到一个规模为k的优化解集合H’＝{H₁’，...，H_k’}；

采用二进制串X(x₁，x₂，...，x_n)(x_i∈{0，1})和V(v₁，v₂，...，v_n)(v_i∈{0，1})来分别表示一个粒子的位置和速度，粒子的速度表示及粒子间的运算操作定义如下：

(a)粒子的速度V定义为其两次位置X₁和X₂之间的距离；

V＝X₁-X₂＝(v₁，...，v_n)，

$v_i=\left\{\begin{array}{c}0,x_{1i}=x_{2i}\\ 1,x_{1i}\≠x_{2i}\end{array}\right.i=\mathrm{1,2},...,n$

(b)速度V₁和V₂间的加法操作定义为其相应位的逻辑加，结果为速度V；

V＝V₁+V₂＝(v₁，...，v_n)，v_i＝v_1iORv_2i，i＝1，2，...，n}.

(c)粒子速度V₁与概率C的乘积，结果为速度V。

V＝CV₁＝(v₁，...，v_n)，

(d)速度V和位置X₁间的加法操作定义为其相应位的逻辑异或，结果为位置X。

X＝X₁+V＝(x₁，...，x_n)，x_i＝x_1iXORv_i，i＝1，2，...，n}.

适应度函数用于评价粒子的搜索性能，指导粒子群的搜索过程；某个粒子位置X及矩阵M1中的所有片断f_i(i＝1，...，m1)，X的适应度函数Fitness(X)定义为：

$\mathrm{Fitness}\left(X\right)=1-\frac{E\left(X\right)}{m1\×n1},$

$E\left(X\right)=\underset{i=1}{\overset{m1}{\mathrm{\Σ}}}\min (S(f_i,X),D(f_i,X))$

其中，粒子位置X表示一对仅含杂合位点的单体型(h₁’，h₂’)中的一条，S(f_i，X)表示f_i和X对应位取值相同的位点个数，D(f_i，X)表示表示f_i和X对应位取值相同的位点个数；E(X)表示对应于单体型对(h₁’，h₂’)的最少错误更正数；

(3)最后的扩展阶段将预处理阶段删掉的SNPs重新加上，对于优化解集H’中只含杂合位点的单体型对H_i’＝(h_i1’，h_i2’)，如果某个已被删除的同合位点为0-列或1-列，则将0或1插回到单体型对(h_il’，h_i2’)的相应位置，以此得到扩展后的单体型对H_i＝(h_i1，h_i2)，扩展结束后得到最终的单体型对集合H＝{H₁，...，H_k}。

Images