基于双正交非负矩阵三因式的肿瘤识别

谭青青，王 年，苏亮亮，方正文

（安徽大学 计算智能与信号处理教育部重点实验室，安徽 合肥 230039）

DNA微阵列是分子生物学在实验领域的一项重大突破，它已经成为分子诊断的一个重要手段.随着DNA微阵列技术的出现和不断发展，产生的基因表达谱数据规模日趋庞大，而这类数据的主要特点是基因维数大而样本维数小.因此，如何对肿瘤基因进行有效分析，继而从海量的数据中提取出有意义的生物学信息就成了生物信息学研究的重点问题之一［1－2］.通过检测基因表达的变化，可以为疾病的诊断和预测提供新的目标.到目前为止，已经有很多算法适合于处理高维数据集，如支持向量机［3］、随机森林［4］、贝叶斯法［5］等，这些算法都已经应用到了各种肿瘤样本的分类和聚类中.各种特征选择的方法也已经被提出并不断发展，以便于准确鉴定相关疾病基因.当前对于基因表达谱数据分析和挖掘的一般方法是对基因数据进行特征提取后，再用相应的分类或聚类方法对肿瘤样本进行识别［6－7］.1999年，Golub等［8］以“信噪比”为指标对白血病的两个亚型样本进行分类研究；Lee和Seung等首次提出了NMF［9］的概念，它是一种特征提取和降维的新方法.随着NMF在现实生活中的广泛应用，NMF算法的缺点也显现出来了.2008年，李乐等对各种非负矩阵分解算法进行了总结［10］；Cichocki等把稀疏惩罚项和SED（square of euclidian distance）结合作为目标函数，并对迭代规则进行了修改，得到了稀疏性增强的NMF［11］；Hoyer提出了可精确控制稀疏性的算法 NMFSC（NMF with sparseness constraints）［12］.Wang等把Fisher鉴别信息作为惩罚项加入到目标函数中，构造了FNMF（Fisher NMF）算法［13］，使得分解结果中蕴含较多的分类信息.

论文所介绍的正交非负矩阵三因式分解算法是一种改进的NMF方法［14］.一般的NMF算法采用随机初始化，而论文将K－均值聚类作为BONMTF算法的初始化步骤，这样做提高了算法的效率.实验采用了4组具有代表性的肿瘤基因表达谱数据进行研究，其结果证明了该方法的可行性和有效性.

1 NMF算法

1999年，Lee和Seung第一次提出了NMF的基本框架［9］.假设矩阵X＝（xij）m×n为一个m×n的矩阵，矩阵X中的所有元素值都是非负的.非负矩阵分解（NMF）算法就是将矩阵X分解成两个非负矩阵其中：矩阵U∈Rm×k＋，V∈Rk×n＋，通常k的值要远远小于m和n，则获得的子矩阵U和V的维数得到了降低，一个非负矩阵分解成了另外两个非负矩阵相乘的形式.

矩阵X的一行可以写成矩阵V的列的线性组合，即

其中：xj表示矩阵X的一行，vi表示矩阵V的一列，而uji表示矩阵U中的元素.此模型如图1所示，在这个模型中，可以将矩阵V的k行作为基向量，而将矩阵U的行作为编码系数.此模型可以用于肿瘤分类的特征提取中，特征为矩阵U的行.经过优化U和V使之与矩阵X线性近似.

矩阵分解实际上是数据最优化过程，因此需要定义一个目标函数来量化近似度，并通过更新迭代规则来减小此函数.此目标函数可定义为

当UV与X偏离越大时，R的值也随之增大；当UV越接近于X时，R的值随之减小.从而R的大小能够衡量UV与X的近似程度.使用迭代逼近的方法找到满足式（1）的两个矩阵，并使得R为非增的，则可使用如下的迭代规则

2 双正交非负矩阵三因式分解BONMTF

2.1 BONMTF算法

针对NMF算法的缺点，2006年，Chris Ding和Tao Li等第一次提出了双正交三因式非负矩阵分解和对称三因式非负矩阵分解的概念［14］，并给出了BONMTF迭代公式的计算规则.论文主要将正交非负矩阵三因式分解用于文件聚类，并给出了二因式NMF迭代规则的详细证明.

假设X＝（xij）m×n为一个m×n的矩阵，其每一行代表在同一样本中所有基因的表达值，X中的所有元素值都为非负的.而每一列代表了一个基因在不同样本中的表达值.BONMTF分解就是分解成3个矩阵相乘的形式

使得

然而只有当三因式NMF不能转化为二因式NMF时，此研究才有意义.也就是说需要给三因式NMF施加一定的约束.显然式（7）没有与之相对应的二因式，称之为双正交三因式分解，且其为论文的研究重点.可以通过使用以下的更新规则来计算D的值

2.2 论文算法

论文将BONMTF算法应用于基因表达谱数据的分析.论文算法的基本步骤为：

（1）利用Fisher判别准则对基因表达谱数据进行预处理，得到去除无关基因后的基因子集，并对所得的基因子集数据进行归一化处理.

对基因表达谱数据进行深入分析和挖掘，是基因表达谱数据分析的根本内容［1］.但是，由于癌症基因表达谱数据的复杂性，所得的微阵列数据中只有较少的基因包含和类别相关的信息，因此要对基因表达谱数据进行预处理.简单的处理方法即按照某种记分准则来对每一个基因进行记分，所得分值越大则说明该基因越重要；再按照分值的大小对基因进行降序排列，选取排在前面的一定数量的基因作为筛选的基因子集，从而初步去除了无关基因和噪声，即实现了样本维数的约减.论文选取Fisher判别准则对基因表达谱数据进行预处理，其基因的记分准则表示为

（2）利用BONMTF算法提取特征向量，以此来进一步剔除冗余基因，得到能够表征样本属性的矩阵.

3 用SVM分类器实现肿瘤样本类型的识别

在二因式NMF的基础上，给出了BONMTF算法的正确性和收敛性的证明，即旨在解决以下优化问题

使得

其中：X为非负输入矩阵.

根据约束优化问题的标准理论，引入了拉格朗日乘数λ，并最小化拉格朗日函数L1（F），如下式

L1（F）的梯度为

由非负Fik的KTT互补条件，得

为了求解式（15）的耦合方程以及求解满足约束条件FTF＝I的F和λ，可以利用非线性的方法如牛顿法求解.然而非线性方程一般是很难解的，在这里可以利用一个更简单的算法即利用式（9）的更新规则来求解.但是在这里出现了两个问题，即该算法是否是正确且收敛的.因此，下面给出了该算法正确性和收敛性的证明.

正确性：首先给出F的一个初始化预测，然后通过利用以下公式的连续更新而收敛到一个局部最小值.

正确性是由收敛性来保障的，其解将满足式（15），拉格朗日乘数λ将在接下来的式（27）中给出.

收敛性：单调性定理保证了收敛性.论文利用命题1和定理1［14］来证明收敛性.

证明 令Sip＝S′ipuip，则差值Δ＝LHS－RHS可以写成

由于矩阵A和B都为对称矩阵，则上式等价于

当B＝I，S为一个列向量.命题1为证明定理1做好了准备.

定理1 根据式（16）的迭代规则，假定SGTGSΤ＋λ≥0时，拉格朗日函数L1（F）为单调递减（非增）函数.

证明 在这里使用辅助函数的方法［15］，一个函数Z（H，H，）若满足以下条件，则被称为L（H）的辅助函数

对于任意的H和，定义

Z（F，F′）为L1（F）的一个辅助函数.由命题1可得，其满足式（20）的条件.现在，根据式（21）可知，需要找到固定F′时，f（F）＝Z（F，F′）的全局最小值，而局部极小值可由下式给出

求解Fik，极小值为

由于海赛矩阵∂Z（F，F′）／Fik∂Fjl是正定的.因此，这是一个凸函数，且极小值即为全局最小值.在式（16）中，有（－FΤXGSΤ＋SGΤGSΤ＝λ）kk＝0.因此，可以获得拉格朗日的对角元素为

拉格朗日因子的非对角元素是通过设定∂L／∂Fik＝0近似获得的，由式（14），可以得出下式

结合式（25）和（26），可以得到拉格朗日乘数

将式（27）代入式（16）中，可以获得式（9）的迭代规则.

这里值得注意的是，由于拉格朗日因子Λ＝λkl的非对角元素是近似获得的，因此最后的解F并不恰好满足FΤF＝I.这实际上却是一个优势，如果FΤF正好等于I，根据非负性，F的每一行只能有一个非零元素.而对FΤF＝I的轻微偏离却能够允许一行中的所有元素非零（尽管大多数值都很小）.

类似地，对于给定的F和S，可以使用式（8）的迭代规则来更新G，证明的方法与上述方法类似.综上所述，论文可以通过式（8）、（9）来证明式（7）的最小化过程；对于给定的F和G，可使用式（10）来更新矩阵S的迭代规则，下文给出了证明［14］.

定理2 令F和G为固定的矩阵，则

是单调递减的.

证明 首先需要证明正确性.利用Sik的非负性和KTT互补条件，可以得出

由于式（10）的迭代规则满足上式，则证明了式（10）的正确性.接下来需要考虑式（10）的收敛性.

定理3 目标函数J2（S）在式（10）的迭代规则下是为非增的.

证明 同样可以使用辅助函数的方法来证明.其中

为J2（S）的一个辅助函数，根据式（21），当固定S′＝S（t）时，通过最小化J（S，S′）可以得到S（t＋1），最小值为

根据式（21），可以更新式（10）.在这个更新规则下，J2（S）为单调递减的.

由于一般的NMF算法采用随机初始化，这导致算法的效率低，且算法的复杂度较高.为了提高算法的效率及降低算法复杂度，论文算法是首先使用K－均值聚类作为BONMTF的初始化步骤，即利用K－均值聚类产生矩阵F和G的初始矩阵，这样做的好处是减小了最优解的寻找区间，即可以节省时间，且使得结果更加稳定.对矩阵S，采用下式来进行初始化

由于在初始化过程中，F、G以及矩阵X都是已知的，从而根据上式可以求得初始矩阵S.

4 实验结果与分析

4.1 数据实验

选用4组公开的且具有代表性的数据集：① 白血病数据（Leukemia），共52个样本，其中24个样本被确定为急性淋巴白血病（acute lymphoblastic leukemia，简称ALL），28个样本被确定为急性粒细胞白血病（acute myelold leukemia，简称AML），每个样本含有12 564个基因表达数据；② 前列腺癌数据集（prostate cancer），共102个样本，其中50个样本为正常样本，52个样本被确定为肿瘤样本，每个样本含有12 600个基因表达数据；③ 弥漫大B细胞淋巴瘤（diffuse large B－cell lymphoma，简称DLBCL）数据集，共77个样本，其中58个样本为正常样本，19个样本被确定为肿瘤样本，每个样本含有5 469个基因；④ 结肠癌数据（colon cancer），共62个样本，其中22个样本为正常样本，40个被确定为肿瘤样本，每个样本含有2 000个基因表达数据.

利用Fisher准则分别对白血病、前列腺癌、结肠癌以及弥漫大B细胞瘤的基因表达谱数据计算每个基因的得分，如图2所示.由图2可以看出有些基因的得分高而有些基因的得分低，而分值越高则表明含有越多的分类信息，所以需要保留得分高的基因，去除得分低的基因，以此来对基因集合进行初选，去除无关基因.

4.2 实验分析

论文采用留一法进行测试，即在整个数据集中，始终保留一个样本作为测试样本，剩下的样本作为训练样本，直至所有的样本都作为一次测试样本为止.

癌症数据集的实验按照如下步骤进行：先利用Fisher准则计算所有基因的得分，并通过对所计算的得分进行降序排列，再选取排在前面的一定数量的基因作为筛选的基因子集.通过大量实验发现选取前90个基因作为基因子集时，可以获得较好的实验结果；将所得的基因子集数据进行标准化处理，再利用BONMTF提取特征向量，以此来进一步剔除冗余基因，得到能够表征样本属性的矩阵F，然而F的维数k的不同将会影响最终的分类结果.因此要选择合适的k值，这样不仅使得类别信息得到凸显，还使得冗余信息得到消除.图3给出了在不同k值的情况下4个癌症数据集的分类情况.最后利用SVM分类器实现样本肿瘤类型的识别.

由图3可知，对于白血病、结肠癌和弥漫大B细胞淋巴瘤而言，当k值为2时，在实验中取得的样本分类的正确识别率相对较高，随着k值的增加，样本分类的正确识别率慢慢降低，最终趋于稳定.而对前列腺癌来说，当k值为3时分类正确率相对较高.这表明，当k值为2和3时，样本属性矩阵F最能反映样本的分类特性.如在白血病的实验中，正确识别率接近100%，这说明在样本属性矩阵F中包含了完整的分类信息，可以对样本集合中的所有样本进行正确分类.因此此时利用BONMTF算法获得的矩阵F为最佳的分类特征矩阵.

然而在基因表达谱数据分类的方法中，很多方法并没有普遍适用性，即可能针对某种数据的分类效果很好，但是对于其他数据，其分类效果并不理想.为了进一步验证论文的可行性、有效性以及普遍适用性，作者将论文方法与传统的S2N－SVM［16］法、Cluster－S2N［17］方法以及传统的 NMF法进行比较.S2N－SVM以“信噪比”为指标来选取特征基因，再利用SVM分类器对样本进行分类.Cluster－S2NSVM方法中的特征选取和分类器有关，先从初始特征集合中选择一个特征子集，并在此子集所包含的信息的基础上构建新的特征集合，最后结合SVM实现分类.如表1所示.

表1 论文方法与另外3种方法的实验结果比较Tab.1 Comparisons of experiment results with other three methods

由表1可知，对于4组不同的肿瘤数据集，论文算法的识别率都达到了90%以上，明显优于其他3种方法.这是因为该算法能够提取到1组数量更少却更具样本分类能力的特征基因.虽然传统NMF方法对于白血病的识别率也接近100%，但是对于其他3组数据集，其分类效果不如论文的方法.这说明论文算法具有相当的可行性以及普遍适用性.

5 结束语

NMF已经被广泛地用于图像分析、论本聚类以及癌症的发现与分类.论文结合双正交和三因式的思想，提出了一种改进的NMF算法，称之为BONMTF算法，并将其应用于基因表达数据分析.该方法包括高维基因表达数据的基因选择和特征提取两个步骤，最后利用SVM进行分类.论文计算了BONMTF算法的更新规则，并证明了此算法的收敛性和正确性.通过4组癌症数据的实验，证明了该方法在预测人体组织中正常样本和肿瘤样本的有效性和高效性，即有效地提高了预测的精度并降低了算法的时间复杂度，具有广泛的应用前景.

［1］黄德双.基因表达谱数据挖掘方法研究［M］.北京：科学出版社，2009：20－66.

［2］庄振华，王年，李学俊，等.癌症基因表达数据的熵度量分类方法［J］.安徽大学学报：自然科学版，2010，34（2）：73－76.

［3］Buturovic L J.PCP：aprogram for supervised classification of gene expression profiles［J］.Bioinformatics，2006，22（2）：245－247.

［4］Nannapaneni P，Hertwig F，Depke M，et al.Defining the structer of the general stress regulon of Bacillus subtilis using targeted microarray analysis and random forest classification［J］.Microbiology，2012，158（3）：696－707.

［5］Nanni L，Brahnam S，Lumini A.Combining multiple approaches for gene microarray classification［J］.Bioinformatics，2012，28（8）：1151－1157.

［6］阮小刚，李晓明，王金莲.边介数聚类算法在肿瘤基因表达谱中的应用［J］.北京工业大学学报：自然科学版，2008，34（7）：696－700.

［7］Ancherla K，Mukkamala S.Feature selection for lung cancer detection using SVM based recursive feature elimination method［M］.Berlin Heidelberg：Springer，2012，168－176.

［8］Golub T R，Slonim D K，Tamayo，et al.Class discover and class prediction by gene expression monitoring［J］.Molecular classification of Cancer，1999，256（5439）：531－527.

［9］Lee D D，Seung H S.Learning the parts of objects by nonnegative matrix factorization［J］.Nature，1999，401（6755）：788－791.

［10］李乐，章毓晋.非负矩阵分解算法综述［J］.电子学报，2008，36（4）：737－743.

［11］Cichocki A，Amory S，Zdunek R，et al.Extended SMART algorithms for non－negative matrix factorization［M］.Berlin Heidelberg：Springer，2006：548－562.

［12］Hoyer P O.Non－negative factorization with sparseness constraints［J］.J of Mach Learning Res，2004，5（9）：1457－1469.

［13］Wang Y，Jia Y，Hu C，et al.Fisher non－negative matrix factorization for learning local features［C］／／Proc of Asian Conf on Comp Vision，2004：27－30.

［14］Ding C，Li T，Peng W，et al.Orthogonal nonnegative matrix tri－factorization for clustering［C］／／Proceedings of the 12th ACM SIGKDD International Conference on Knowledge Discovery and Data Mining，2006：126－135.

［15］Lee D D，Seung H S.Algorithms for non－negative matrix factorization［J］.Advances in Neural Information Processing Systems，2001，13（1）：556－562.

［16］Singh D，Febbo P G，Ross K，et al.Gene expression correlates of clinical prostate cancer behavior［J］.Cancer Cell，2002，1（2）：203－209.

［17］阮小刚，晁浩.肿瘤识别过程中特征基因的选取［J］.控制工程，2007，14（4）：373－375.