CN104298893A - 一种基因表达缺失数据的填补方法 - Google Patents

Abstract

本发明公开了一种基因表达缺失数据的填补方法，包括如下步骤：首先对于给定的目标基因，选取与目标基因最相似的k个近邻基因；其次利用选出的k个近邻基因构建相应的目标函数，并采用最小二乘准则和拉格朗日方法迭代地求解回归系数矩阵和对角权重矩阵；再根据不同的回归情况，分别采用不同的填补公式对缺失数据进行填补；最后利用一个不确定性度量对本发明提出的方法进行循环迭代，直到缺失值不再变化。本发明有效地解决了基因表达缺失数据的填补问题，使得对基因表达数据分析的结果更为可靠，为后续的基因表达数据分析提供更为有利的数据支持。

Description

一种基因表达缺失数据的填补方法

技术领域

本发明涉及生物信息技术领域，尤指一种基因表达缺失数据的填补方法。 

背景技术

随着生物技术的不断发展，特别是DNA微阵列技术的开发，已经允许人们在基因水平上揭示有机体的奥秘。由DNA微阵列技术产生的数据亦称为基因表达数据，目前已广泛地应用于致病基因的发现，疾病的辅助诊断，药物疗效的判断等领域。然而，由于各种各样的原因，基因表达数据存在大量的缺失值。根据调查，所有的基因表达数据都含有一定程度的缺失现象，有时候超过90％的基因都会受到缺失数据的影响。而且现有的大多数数据分析工具都是针对完整的数据设计的，无法直接应用于含有缺失的基因表达数据。因此，准确的填补基因表达缺失数据将具有重要的理论意义和实用价值。 

解决数据缺失的最简单方法就是直接删除或忽略缺失值。另外，还有一些比较简单的缺失数据处理方法，如使用“0”填充，或者使用数据的行平均值代替，或者列平均值代替。但由于这些方法没有考虑到数据之间的关联，也没考虑到基因表达数据本身具有的结构，因此这样填补的数据常常具有较差的性能和较高的估计误差。近年来，随着生物信息技术的深入研究与发展，人们尝试借鉴其他领域缺失数据填补的一些思想，提出了许多基因表达缺失数据的填补方法。按照算法思想分类，可将常用的填补算法分为五类：基于k近邻的方法(KNN)、基于奇异值分解的方法(SVD)、基于贝叶斯主成分分析的方法(BPCA)、基于最小二乘准则的方法以及基于聚类的方法等。 

基于最小二乘准则的方法，由于其思想简单，数学表达清楚，填补效果显著等特点而得到了广泛的关注，并取得了大量的研究成果，如局部最小二乘方法(LLS)，迭代的局部最小二乘方法(ILLS)、序列局部最小二乘方法(SLLS)、权重局部最小二乘方法(WLLS)、迭代的双聚类局部最小二乘方法(BI-ILS)等。最小二乘准则填补方法都采用两步实现：第一步，针对需要填补的基因，预选出最相似的k个近邻基因；第二步，运用预选的基因采用最小二乘准则进行缺失数据填补。多元回归分析是所有最小二乘准则填补方法的一个必须而重要的步骤，多元回归的结果进一步用于基因表达缺失数据的填补。然而，现存的所有最小二乘准则方法中，多元回归分析都没有考虑到不同近邻基因的不同重要性。为此，本发明考虑了目标基因的k个最近邻基因的不同重要性，提出了一种自动地加权不同近邻基因的缺失数据填补方法。 

发明内容

鉴于现有的基于最小二乘准则填补方法存在的问题，本发明的目的是提供一种考虑了不同近邻基因的不同重要性的缺失数据填补方法——迭代的局部自加权最小二乘填补方法(ILAW-LS)。 

具体实施步骤如下所示： 

一种基因表达缺失数据的填补方法，自动加权不同近邻基因并对基因表达缺失数据进行填补，包括如下步骤： 

步骤1：输入具有缺失数据的基因表达矩阵其中，m为基因的数量，n为样本的数量； 

步骤2：选取目前尚未被填补的缺失数据行作为目标基因，假设该目标基因具有p个缺失样本值；再选取与目标基因最相似的k个近邻基因，并构建相应的矩阵A和B，其中，矩阵A由k个近邻基因的p个缺失样本所在的列构成，矩阵B由k个近邻基因的q＝n-p个未缺失样本所在的列构成； 

步骤3：采用最小二乘准则和拉格朗日方法，迭代地求解回归系数矩阵和对角权重矩阵

步骤3.1：随机初始化对角权重矩阵W⁰(k)，采用最小二乘方法计算回归系数矩阵Y⁰使得||W⁰(k)(AY-B)||_F达到最小，并计算回归残差矩阵的F范数值ErrorF＝||AY⁰-B||_F；设置最小的残差F范数值MinErrorF＝ErrorF和相应的回归系数矩阵Y_min＝Y⁰，并设定初始的迭代次数t＝0以及最大的迭代次数T； 

步骤3.2：令采用最小二乘方法计算回归系数矩阵Y^t+1使得 达到最小，并计算回归残差矩阵的F范数值ErrorF＝||AY^t+1-B||_F；如果ErrorF<MinErrorF，则修正最小的残差F范数值MinErrorF＝ErrorF和相应的回归系数矩阵Y_min＝Y^t+1；如果或者t达到最大的迭代次数T，则跳转到步骤4；否则，跳转到步骤3.3。其中，P(W(k),Y)为本发明的目标函数，ε为预先给定的较小的常数值； 

步骤3.3：令采用拉格朗日方法计算权重矩阵W^t+1使得带约束的目标函数 达到最小；如果则跳转到步骤4；否则，重置迭代 次数t＝t+1，并跳转到步骤3.2； 

步骤4：如果迭代次数t达到最大的迭代次数T，跳转到步骤5；否则，采用表达式u^T＝v^TY^t+1来填补相应的缺失数据，并且跳转到步骤6；其中，u^T为目标基因的p个缺失样本列，v^T为目标基因的q个未缺失样本列； 

步骤5：如果则采用表达式u^T＝v^TY^t+1来填补缺失数据，并跳转到步骤6；否则，采用表达式u^T＝v^TY_min来填补缺失数据，并跳转到步骤6； 

步骤6：如果所有的缺失数据都被填补完全，则以现有的填补值为基础，多次迭代的填补缺失数据直到填补达到稳定状态或者达到最大迭代次数T；否则，跳转到步骤2； 

步骤7：输出填补完全的基因表达矩阵； 

采用这样的处理，本发明方法能根据近邻基因的不同重要性自动地更新近邻基因的权重，从而自动地强调重要的近邻基因并淡化包含噪声的近邻基因。该方法在局部最小二乘法(LLS)的基础上，添加了一个新的步骤来计算近邻基因的权重，并采用拉格朗日法给出了相应的最优权重计算公式。为了提高该方法的收敛速度，本发明设计了一个加速策略，该策略保证在最差的情况下也能具有局部最小二乘方法(LLS)的收敛速度。 

此外，本发明还设计了一个迭代的填补框架来进一步提高填补的准确性。所述的步骤6中多次迭代的迭代框架可以采用以下规则来实现：对于某个给定的缺失样本，只有当不确定性下降的时候，才以现有的迭代结果替换上一次的值。其中，不确定性度量θ的计算公式如下： 

$\mathrm{\&theta;}=t_{\mathrm{\&alpha;}/2,u-r}\sqrt{v^T[{\left(A^{T}A\right)}^{-1}v+1]{\widehat{\mathrm{\&sigma;}}}^2}$

其中，t_α/2,u-r为置信水平为α自由度为u-r的t分布，u和r分别为矩阵A行和列的数目， 为多元回归方差的无偏估计值，其计算公式如下所示： 

${\widehat{\mathrm{\&sigma;}}}^2={(B_j-A{\hat{Y}}_j)}^T(B_j-A{\hat{Y}}_j)/(u-r)$

其中，B_j是矩阵B的第j列，是矩阵的第j列。 

与现有技术相比，本发明的积极效果是： 

一、区别于经典的局部最小二乘填补方法(LLS)，本发明考虑了不同近邻基因的不同重要性，并通过加权近邻基因，实现了自动地强调重要近邻基因而淡化包含噪声的近邻基因，使得本发明的模型具有一般性和灵活性。 

目前所有的基于最小二乘准则的填补方法都没有考虑到目标基因的不同近邻基因的不同重要性，而本发明在经典的局部最小二乘方法(LLS)的基础上，通过在经典LLS方法的目 标函数的左边乘以对角权重矩阵来构建新的目标函数，并利用最小二乘准则和拉格朗日方法迭代地求解相应的回归系数矩阵和对角权重矩阵，从而实现了自动加权近邻基因的目的，使得新的模型具有一般性和灵活性。一般性表现为：经典的局部最小二乘方法(LLS)是本发明的一种特殊情况；灵活性表现为：本发明的模型可以表达更多的权重最小二乘方法可能出现的情况。 

二、区别于经典的权重回归方法，本发明能给出最优的权重计算公式和回归系数矩阵。 

经典的权重回归方法大多数采用经验公式或者回归残差的方差来确定权重值，采用最优化的方法来确定回归权重的研究比较少。相对于大多数权重回归方法，本发明通过最小二乘准则和拉格朗日方法迭代地求解带约束的目标函数，进而给出相应的权重矩阵和回归系数矩阵，具有更科学的权重选择方式和回归系数矩阵计算方法。 

附图说明

图1是不同缺失比例下的对比结果图 

图2是加权策略鲁棒性实验结果图 

图3是收敛性实验结果图 

具体实施方式

具体实施步骤如下所示： 

步骤1：输入具有缺失数据的基因表达矩阵其中，m为基因的数量，n为样本的数量。 

步骤2：选取目前尚未被填补的缺失数据行作为目标基因，假设该目标基因具有p个缺失样本值；然后选取与目标基因最相似的k个近邻基因，并构建相应的矩阵A和B。 

所述的最相似的k个近邻基因的计算方法如下：首先，选取p个缺失样本所在的列全都完整的行作为候选基因，并对候选基因其他列的缺失部分采用行平均进行填补；然后，忽略p个缺失样本所在的列，计算目标基因与候选基因之间的相似度，其中相似度可以采用欧式距离或者皮尔逊相关系数等来度量；最后，选取相似度绝对值最大的k个候选基因作为与目标基因最相似的k个近邻基因。 

所述的矩阵A由k个近邻基因的p个缺失样本所在的列构成，矩阵B由k个近邻基因的q＝n-p个未缺失样本所在的列构成。 

步骤3：采用最小二乘准则和拉格朗日方法，迭代地求解回归系数矩阵和对角权重矩阵

步骤3.1：随机的初始化对角权重矩阵W⁰(k)，采用最小二乘方法计算回归系数矩阵Y⁰使得||W⁰(k)(AY-B)||_F达到最小，并计算回归残差矩阵的F范数值ErrorF＝||AY⁰-B||_F。设置最小的残差F范数值MinErrorF＝ErrorF和相应的回归系数矩阵Y_min＝Y⁰，并设定初始的迭代次数t＝0和最大的迭代次数T。 

步骤3.2：令采用最小二乘方法计算回归系数矩阵Y^t+1使得 达到最小，并计算回归残差矩阵的F范数值ErrorF＝||AY^t+1-B||_F。如果ErrorF<MinErrorF，则修正最小的残差F范数值MinErrorF＝ErrorF和相应的回归系数矩阵Y_min＝Y^t+1。如果或者t达到最大的迭代次数T，则跳转到步骤4；否则，跳转到步骤3.3。其中，P(W(k),Y)为本方法的目标函数，ε为预先给定的较小的常数值。 

所述的本方法的目标函数是为了计算最优的权重矩阵而设计，其设计可以采用以下方法： 

令 $A=\left(\begin{array}{cccc}{a}_{11}& a_{12}& ...& a_{1q}\\ a_{21}& a_{22}& ...& a_{2q}\\ .& .& .& .\\ .& .& .& .\\ .& .& .& .\\ a_{k1}& a_{k2}& ...& a_{\mathrm{kq}}\end{array}\right),$ $B=\left(\begin{array}{cccc}{a}_{11}& a_{12}& ...& a_{1p}\\ a_{21}& a_{22}& ...& a_{2p}\\ .& .& .& .\\ .& .& .& .\\ .& .& .& .\\ a_{k1}& a_{k2}& ...& a_{\mathrm{kp}}\end{array}\right),$ 则本方法的目标函数为： 

$\underset{W\left(k\right),Y}{\min }P(W\left(k\right),Y)=\left|\right|W\left(k\right)(\mathrm{AY}-B){\left|\right|}_F^2=\underset{i=1}{\overset{k}{\mathrm{\&Sigma;}}}{w}_i^{2\mathrm{\&beta;}}\left(k\right)\underset{j=1}{\overset{p}{\mathrm{\&Sigma;}}}{(\underset{i=1}{\overset{q}{\mathrm{\&Sigma;}}}{a}_{\mathrm{il}}{y}_{\mathrm{lj}}-b_{\mathrm{ij}})}^2$

$s.t\underset{i=1}{\overset{k}{\mathrm{\&Sigma;}}}{w}_i\left(k\right)=\mathrm{1,0}\&le;w_i\left(k\right)\&le;1$

其中，为第i个近邻的权重，β为权重的光滑系数。 

步骤3.3：令采用拉格朗日方法计算权重矩阵W^t+1使得带约束的目标函数 达到最小。如果则跳转到步骤4；否则，重置迭代次数t＝t+1，并跳转到步骤3.2。 

步骤4：如果t达到最大的迭代次数T，跳转到步骤5；否则，采用表达式u^T＝v^TY^t+1来填补相应的缺失数据，并且跳转到步骤6。 

其中，u^T为目标基因的p个缺失样本列，v^T为目标基因的q个未缺失样本列。 

步骤5：如果则采用表达式u^T＝v^TY^t+1来填补缺失数据， 并跳转到步骤6；否则，采用表达式u^T＝v^TY_min来填补缺失数据，并跳转到步骤6。 

步骤6：如果所有的缺失数据都被填补完全，则以现有的填补值为基础，多次迭代的填补缺失数据直到填补达到稳定状态或者达到最大迭代次数T；否则，跳转到步骤2。 

所述的迭代的填补缺失数据可以采用以下规则来实现：对于某个给定的缺失样本，只有当不确定性下降的时候，才以现有的迭代替换上一次的值。其中，不确定性度量θ的计算公式如下： 

$\mathrm{\&theta;}=t_{\mathrm{\&alpha;}/2,u-r}\sqrt{v^T[{\left(A^{T}A\right)}^{-1}v+1]{\widehat{\mathrm{\&sigma;}}}^2}$

其中，t_α/2,u-r为置信水平为α自由度为u-r的t分布，u和r分别为矩阵A行和列的数目， 为多元回归方差的无偏估计值，其计算公式如下所示： 

${\widehat{\mathrm{\&sigma;}}}^2={(B_j-A{\hat{Y}}_j)}^T(B_j-A{\hat{Y}}_j)/(u-r)$

其中，B_j是矩阵B的第j列，是矩阵的第j列。 

步骤7：输出填补完全的基因表达矩阵。 

本发明中设计了一个加速策略来有效地提高本发明的收敛速度，即使在最差的情况下本发明也具有局部最小二乘方法(LLS)的收敛速度，该加速策略可以通过如下规则实现的：当目标函数P(W(k),Y)达到最优的时候，采用u^T＝v^TY^t+1来填补缺失数据；如果目标函数P(W(k),Y)未能达到最优的时候，则在先前的权重迭代解中选择满足残差F范数值MinErrorF最小的解作为回归系数矩阵，并采用u^T＝v^TY_min来填补缺失数据。 

有益结果 

为了验证本发明的有效性，本发明从1998年Spellman等人公布的酿酒酵母细胞周期的基因表达数据中选取alpha因子数据作为本发明的验证数据(简记为SP.Aphla)，并使用NRMSE(标准均方根误差)作为本发明的评估度量，其计算公式如下所示： 

$\mathrm{NRMSE}=\frac{\sqrt{\frac{1}{\left|S\right|}{\mathrm{\&Sigma;}}_{(i,j)\&Element;S}{g}_{\mathrm{ij}}-{{\hat{g}}_{\mathrm{ij}}}^2}}{\sqrt{\frac{1}{\left|S\right|-1}{\mathrm{\&Sigma;}}_{(i,j)\&Element;S}{g}_{\mathrm{ij}}-{\stackrel{\&OverBar;}{g}}^2}}$

其中，g_ij是真实的基因表达值，是估计值，S是缺失样本构成的集合，|S|表示缺失样本的数量，是所有缺失样本的平均值。 

试验一 

在SP.Aphla数据集上，分别以1％，5％，10％，15％，20％的缺失比例对本发明的方法进行试验，并与四种流行的基因表达缺失数据填补方法进行了对比，相关结果如图1所示。这四种流行的方法分别是：基于k近邻的方法(KNN)、局部最小二乘方法(LLS)、迭代局部最小二乘方法(ILLS)和双聚类基于的局部最小二乘方法(Bi-ILS)。由NRMSE的定义可知，对比结果中的NRMSE值越低，表示填补方法的准确度越高。从图1中可以看出，对试验中给出的所有缺失比例的数据，本发明给出的方法都具有最好的效果，并且在较低的缺失比例下，本发明还具有更好的填补效果。如以NRMSE为评估度量，在1％的缺失比例下，本发明提出的方法比其他四种方法至少提高33.2％的填补精度。 

试验二 

随着近邻基因数量的增加，加权近邻基因的策略是否依然有效呢？为了验证这个问题，本发明给出了近邻基因加权策略的鲁棒性实验。本发明以5％缺失比例的SP.Aphla数据为实验对象，逐渐增加近邻基因的数量。从图2中可以看出，随着近邻数量的增加，本发明提出的方法是逐渐趋于稳定的，这种现象也在某种程度上说明了加权策略的鲁棒性。 

试验三 

为了验证本发明中给出的所有迭代算法的收敛性。本发明以20％缺失比例的SP.Aphla数据为实验数据，对比了两个流行的迭代缺失数据填补方法，相关结果如图3所示。从图3中可以看出，本发明给出的算法与Bi-ILSimpute都是收敛的，而ILLSILSimpute是不稳定的。 

Claims (2)

1.一种基因表达缺失数据的填补方法，自动加权不同近邻基因并对基因表达缺失数据进行填补，包括如下步骤： 

步骤1：输入具有缺失数据的基因表达矩阵其中，m为基因的数量，n为样本的数量； 

步骤2：选取目前尚未被填补的缺失数据行作为目标基因，假设该目标基因具有p个缺失样本值；再选取与目标基因最相似的k个近邻基因，并构建相应的矩阵A和B，其中，矩阵A由k个近邻基因的p个缺失样本所在的列构成，矩阵B由k个近邻基因的q＝n-p个未缺失样本所在的列构成； 

步骤3：采用最小二乘准则和拉格朗日方法，迭代地求解回归系数矩阵和对角权重矩阵

步骤3.1：随机初始化对角权重矩阵W⁰(k)，采用最小二乘方法计算回归系数矩阵Y⁰使得||W⁰(k)(AYB)||_F达到最小，并计算回归残差矩阵的F范数值ErrorF＝||AY⁰-B||_F；设置最小的残差F范数值MinErrorF＝ErrorF和相应的回归系数矩阵Y_min＝Y⁰，并设定初始的迭代次数t＝0以及最大的迭代次数T； 

步骤3.2：令采用最小二乘方法计算回归系数矩阵Y^t+1使得 达到最小，并计算回归残差矩阵的F范数值ErrorF＝||AY^t+1-B||_F；如果ErrorF<MinErrorF，则修正最小的残差F范数值MinErrorF＝ErrorF和相应的回归系数矩阵Y_min＝Y^t+1；如果或者t达到最大的迭代次数T，则跳转到步骤4；否则，跳转到步骤3.3。其中，P(W(k),Y)为本发明的目标函数，ε为预先给定的较小的常数值； 

步骤3.3：令采用拉格朗日方法计算权重矩阵W^t+1使得带约束的目标函数 达到最小；如果则跳转到步骤4；否则，重置迭代次数t＝t+1，并跳转到步骤3.2； 

步骤4：如果迭代次数t达到最大的迭代次数T，跳转到步骤5；否则，采用表达式u^T＝v^TY^t+1来填补相应的缺失数据，并且跳转到步骤6；其中，u^T为目标基因的p个缺失样本列，v^T为目标基因的q个未缺失样本列； 

步骤5：如果则采用表达式u^T＝v^TY^t+1来填补缺失数据，并跳转到步骤6；否则，采用表达式u^T＝v^TY_min来填补缺失数据，并跳转到步骤6； 

步骤6：如果所有的缺失数据都被填补完全，则以现有的填补值为基础，多次迭代的填补缺失数据直到填补达到稳定状态或者达到最大迭代次数T；否则，跳转到步骤2； 

步骤7：输出填补完全的基因表达矩阵。 

2.根据权利要求1所述的一种基因表达缺失数据的填补方法，其特征在于：所述的步骤6中多次迭代的迭代框架可以采用以下规则来实现：对于某个给定的缺失样本，只有当不确定性下降的时候，才以现有的迭代结果替换上一次的值。其中，不确定性度量θ的计算公式如下： 

其中，t_α/2,u-r为置信水平为α自由度为u-r的t分布，u和r分别为矩阵A行和列的数目， 为多元回归方差的无偏估计值，其计算公式如下所示： 

其中，B_j是矩阵B的第j列，是矩阵的第j列。