CN107016261B - 基于联合约束非负矩阵分解的差异表达基因辨识方法 - Google Patents

Abstract

本发明公开一种基于联合约束非负矩阵分解的差异表达基因辨识方法。1.用非负矩阵X表示癌症基因表达数据集；2.构造对角矩阵Q和元素全1矩阵E；3.在经典的非负矩阵分解方法中引入流形学习，对系数矩阵G施加正交约束稀疏和约束，得到联合约束非负矩阵分解目标函数；4.求解此目标函数，得到基矩阵F和系数矩阵G的迭代公式；5.对非负数据集X进行半监督非负矩阵分解，迭代收敛后得到基矩阵F和系数矩阵G；6.根据基矩阵F得到评估向量对评估向量中的元素从大到小进行排序，得到差异表达基因；7.通过GO工具对辨识的差异表达基因进行检测分析。本发明能够有效提取癌症数据集中的差异表达基因，能用于发现人类疾病基因数据库中的差异特征，对癌症的早期诊断和针对性治疗有重要的临床意义。

Description

基于联合约束非负矩阵分解的差异表达基因辨识方法

技术领域

本发明公开一种基于联合约束非负矩阵分解的差异表达基因辨识方法，属于模式识别技术领域，可用于辨识癌症差异表达基因，为癌症的早期诊断治疗提供依据。

背景技术

近年来癌症发病率逐年增加，早期诊断率低，死亡率高，发病机制非常复杂。挖掘发现癌症基因表达数据中所包含的相关信息有助于人们深入了解疾病相关表达基因及其调控网络。

随着深度测序技术和基因芯片技术的快速发展，产生了大规模的基因表达谱数据,其中只有少数基因跟细胞癌变有关。很多情况下原始的基因表达谱数据为高维小样本的数值矩阵,因此如何从大量的基因数据中剔除冗余基因，辨识出人类癌症基因表达数据中代表样本属性的差异表达基因成为疾病诊断治疗的重要依据。这就需要解释能力强、高效的矩阵分解方法，目前有很多数据分析方法用于基因表达数据的处理，如主成分分析、奇异值分解等，但效果并不十分理想。例如：尽管原始的基因表达数据是非负的，但是传统的降维方法不能保证因子矩阵元素也是非负的，这就使低秩因子矩阵的生物学解释复杂化。著名的科学杂志《Nature》于1999年刊登了两位科学家D.D.Lee和H.S.Seung提出的非负矩阵分解(Non-negative Matrix Factorization，NMF)算法。NMF算法为人类处理大规模数据提供了一种新的途径，通过对因子矩阵施加非负约束,用一个包含所有基因信息的低秩矩阵来近似表示原始的基因表达数据矩阵，能够实现矩阵维数削减和大规模数据压缩。NMF算法是一种优化机制，能够利用局部特征表达整体特性，为数据的处理提供了一条新的思路，能够得到基于部分的表达,分解结果直观，容易从物理上解释。在生物信息学领域，可以利用分解结果的局部特征来选择DNA阵列中的差异表达基因。

但是非负矩阵分解方法仍然存在改进的空间，例如：由于非负矩阵分解在处理基因表达数据过程中缺乏稀疏性，使其在差异表达特征选择方面没有优势；人类癌症基因表达数据中通常包含一些异常值和噪声，传统的非负矩阵分解方法不能有效的处理异常值和噪声的影响。

发明内容

发明目的：针对已有非负矩阵分解方法的不足，提出了一种联合约束非负矩阵分解的差异表达基因辨识方法，来保证数据降维之后仍然能保持原始数据内在的几何结构，使分解得到的基矩阵更加稀疏，减小了原始基因表达数据中异常值和噪声的影响，提高了算法的鲁棒性，从而更加准确的从癌症基因表达数据集中辨识差异表达基因。

技术思路：受非负矩阵分解及其改进方法启发，首先在目标函数中引入图正则约束，稀疏约束，正交约束得到更加稀疏合理的数据表达；其次利用L_2,1范数约束误差函数，提高算法的鲁棒性，最终形成联合约束非负矩阵分解的目标函数。对该目标函数进行优化求解，得到基矩阵和系数矩阵的迭代准则。然后利用这种新的非负矩阵分解方法来辨识癌症数据集中的差异表达基因。实现方案如下：

一种基于联合约束非负矩阵分解的差异表达基因辨识方法，其特征在于包括如下步骤：

(1)将癌症数据集中的数据进行预处理，表示为非负矩阵其中x_j表示第j个样本的表达基因构成的列向量，R₊表示正实数集合，p表示矩阵X中基因的个数，n表示矩阵X中样本的个数，将X中的所有元素都归一化到(0,1)上；

(2)对基矩阵F和系数矩阵G进行非负约束，对系数矩阵G矩阵进行图正则约束，尽可能多的保留原始基因表达数据的内在几何信息，对系数矩阵G矩阵进行正交约束，可以得到更加有效稀疏的分解结果，利用L_2,1范数约束误差函数来减小异常值和噪声的影响，构建联合约束非负矩阵分解的目标函数：

s.t.F≥0,G≥0，G^TG＝I，将正交约束引入该目标函数，可得：其中，F为基矩阵，G为系数矩阵，L为图拉普拉斯矩阵，I为单位阵，||||_2,1表示矩阵的L_2,1范数，||||₁表示矩阵的L₁范数，Tr(·)表示矩阵的迹，α表示图正则约束项Tr(G^TLG)的权重，γ表示稀疏约束项的权重，β为正交约束的权重，α,β和γ均为正数；

(3)对联合约束非负矩阵分解的目标函数进行优化求解，实现对非负矩阵进行分解，得到相应的基矩阵和系数矩阵其中表示原始基因表达数据的基空间，其中f_i表示第i个样本包含所有的基因信息，表示低维实数空间R^p×K中全体非负矩阵构成的子集，p为基因的个数，K为降维的维数；为系数矩阵，表示实数空间R^K×n中全体非负矩阵构成的子集，n为样本的个数，K为分解后的矩阵维数；

(4)根据分解得到的基矩阵F，对其中每一行求绝对值的和，得到评估向量

(5)对评估向量中的元素从大到小进行排序，中的元素数值越大差异表达程度越高，从而辨识得到差异表达基因；

(6)利用基因本体分析工具(Gene Ontology，GO)对辨识的差异表达基因进行功能解释。

特别地，步骤(3)中，目标函数的优化求解包括如下步骤：

3a)将基矩阵F和系数矩阵G中的元素全部初始化为(0,1)之间的随机数；

3b)设定降维之后的维数K的值，图正则约束权重α，稀疏约束权重γ，正交约束权重β；

3c)设定最大迭代次数,利用迭代公式更新基矩阵F中的元素F_ik，其中Q∈R^n×n为对角阵，其对角元素为ε为无穷小的正数，利用更新系数矩阵G中的元素G_kj，其中i＝1,2,...,p，j＝1,2,...,n，k＝1,2,...,K，Ε∈R^n×k为全1矩阵，W为图正则约束中构造的权重矩阵，为对称阵，B为对角阵，且其对角线元素采用交替更新基矩阵F系数矩阵G的方法，得到局部最优解；

3d)算法收敛后得到基矩阵F和系数矩阵G。

本发明具有以下优点：

本发明将图正则约束施加到目标函数中，能够保留原始基因表达数据的几何结构；并且增加了正交约束和稀疏约束，从而能够得到更加稀疏更能体现局部特性的基矩阵，能够有效的辨识差异表达基因；利用L_2,1范数约束误差函数减小了癌症基因表达数据集中异常值和噪声的影响，提高了算法的鲁棒性。

附图说明

图1为本发明的流程图；

图2为本发明RGSON和经典NMF方法辨识误差收敛性；图2说明本发明收敛速度比经典NMF方法快，收敛后辨识误差比NMF方法小；

图3为降维的维数对辨识误差的影响。图3说明随着降维维数k的增加辨识误差快速减小，k＝10后收敛速度变慢。

具体实施方式

随着深度测序技术和基因芯片技术的快速发展，涌现出大量的基因表达谱数据，因此，寻找合适的数据分析方法来处理庞大的基因表达谱数据成为了生物信息学的研究热点。由于实验条件限制，通常的实验样本通常只有几十个，而基因测序技术可以同时监控上万个基因，因此，分析基因表达谱数据是统计学上典型的奇异值问题——高维小样本问题。通常可以采用降维的方法降低数据的复杂性，提高分析结果的准确性。很多维数处理技术，如主成分分析PCA，奇异值分解SVD等算法已经得到广泛应用。但是它们还存在一些不足之处，PCA的主成分和SVD的特征向量都不是稀疏的，因此这些方法用于特征选择时，就不能给出合理的生物学解释。

NMF作为一种有效的降维方法，因其能产生非负和低秩的结果，已经广泛用于处理海量数据。但是NMF分解结果缺乏稀疏性导致其在差异特征选择方面并不具备优势。在高通量数据应用中，一般认为高维数据镶嵌在低维流形上，因此考虑原始数据的几何结构是十分必要的，因此本方法结合稀疏性、图正则约束和NMF算法进行分析；基因表达谱数据通常是高维的并且含有噪声和冗余信息，利用L_2,1范数再产生稀疏结果的同时减小了噪声和异常值的影响，提高了算法的鲁棒性。

头颈鳞癌也称为鳞状细胞癌(简称鳞癌)，起源于表皮或附件如皮脂腺导管、毛囊、汗腺管等的角朊细胞，可以扩散到身体的其他部位，如淋巴结或肺。如果转移，癌症预后较差，可能是致命的。头颈部鳞癌对化疗药物不敏感和早期转移是预后不良的重要原因，因此肿瘤转移及化疗敏感性相关的关键基因，一直是人们关注的焦点和热点。本发明以头颈鳞癌为例，在公开的癌症数据基因表达集——头颈鳞癌数据基因表达集上进行实验。

结合图1，本发明的具体实施步骤如下：

步骤1,将人类头颈鳞癌基因表达数据进行预处理，表示为非负矩阵p表示样本，n表示基因。人类头颈鳞癌基因表达数据共包含418个样本，20502个基因：

1a)将矩阵X中的值进行去负处理；

1b)将矩阵X中所有元素归一化为(0,1)之间的数。

步骤2,构建近邻图权值矩阵W∈R^n×n：

2a)初始化基矩阵F和系数矩阵G为(0,1)之间的随机数；

2b)初始化E矩阵为元素全为1的矩阵；

2c)结合图3设定低维维数K的值，以及图正则约束权重α，稀疏约束权重γ，正交约束权重β的值，设定迭代次数；

2d)构建近邻图权值矩阵W∈R^n×n：

即数据x_i，x_j看作近邻图中的两个节点，如果这两个数据在数据集中相邻，则在这两个节点之间连边，边的权值为1，否则为0。B是一个对角阵，且其对角线元素等于W的所有行元素或列元素之和，即对角线元素从而得到图拉普拉斯矩阵L＝B-W，L∈R^n×n。

步骤3，构建联合约束的非负矩阵分解的目标函数O：

3a)为了减小异常值和噪声对误差的影响，采用L_2,1范数约束误差项。由非负矩阵X、基矩阵F和系数矩阵G构造误差项||X-FG^T||_2,1；

3b)为了考虑包含在原始数据中的几何结构，在误差函数中引入图正则约束项。由拉普拉斯矩阵L和系数矩阵G构建图正则约束项Tr(G^TLG)，其中Tr(·)表示矩阵的迹；

3c)为了得到更加稀疏的分解结果，增加稀疏性的调节能力，由系数矩阵G构建稀疏约束项||G||₁，其中||·||₁表示矩阵的L₁范数；

3d)将误差项||X-FG^T||_2,1、图正则约束项Tr(G^TLG)、稀疏约束项||G||₁加权构成联合约束的非负矩阵分解的目标函数O＝||X-FG^T||_2,1+αTr(G^TLG)+r||G||₁；

3e)将对系数矩阵G的正交约束写入目标函数：

其中，表示非负数据矩阵，和分别表示对非负矩阵X进行非负分解得到的基矩阵和系数矩阵，L为图拉普拉斯矩阵，I为单位阵，||||_2,1表示矩阵的L_2,1范数，||||₁表示矩阵的L₁范数，Tr(·)表示矩阵的迹，α为控制平滑行的正则化参数，γ为控制稀疏度的调节稀疏，可以控制稀疏约束的强度，β为正交约束参数，可以调节平滑性和稀疏性，α,β和γ均为正数。

步骤4，对步骤1得到的非负矩阵利用本发明方法进行分解，得到对应的基矩阵F和系数矩阵G，其中基矩阵F中包含所有的基因表达信息。本发明采用乘性迭代更新的办法对目标函数进行求解。步骤如下：

4a)引入拉格朗日乘子，构造拉格朗日函数，并对其求一阶偏导数，利用KKT条件可得到基矩阵F和系数矩阵G的迭代准则

其中F_ik为基矩阵F中第i行第j列元素，G_kj为系数矩阵G中第k行第j列元素，i＝1,2,...,p，j＝1,2,...,n，k＝1,2,...,K，Ε∈R^n×k为全1矩阵，W为图正则约束中构造的权重矩阵，为对称阵，B为对角阵，且其对角线元素

4b)利用对基矩阵F中的元素进行迭代更新，

4c)利用更新系数矩阵G中的元素，其中Q∈R^n×n为对角阵，其对角元素为ε为无穷小的正数；

4d)当目标函数值小于规定值，或者达到预定义的最大迭代次数时迭代停止，得到分解后的基矩阵F和系数矩阵G，否则返回4b)。

通过上述迭代得到的基矩阵代表非负矩阵X的基空间，其中表示第i个样本包含所有的基因信息，表示p维向量空间R^p中全体非负向量构成的子集合，p为基因的个数，K为降维的维数；系数矩阵为所有训练样本低维表示的集合， 表示K维向量空间R^K中全体非负向量构成的子集合，n为样本的个数。

步骤5,利用基矩阵F辨识差异表达基因：

5a)对基矩阵F中的每行元素求和得到评估向量

5b)对评估向量中元素素进行降序排列，数值越大，对应基因的差异表达程度越高；

5c)选择评估向量中前100个最大的元素跟疾病基因特征名数据集相比对，得到差异表达基因；

步骤6，本发明的效果通过实验进行说明，在头颈鳞癌基因表达数据集上进行了差异表达基因提取，并对辨识得到的差异表达基因进行分析：

6a)实验使用matlab仿真软件，根据实验结果参数α,β和γ的变化对实验结果影响不大，本实验中参数设置如下α＝0.5，β＝0.9，γ＝0.5，最大迭代次数100次；

6b)利用基因本体分析工具(Gene Ontology，GO)对辨识的差异表达基因进行功能解释；

6c)将联合约束的非负矩阵分解方法GRSON和其他三种非负矩阵分解方法(非负矩阵分解NMF,鲁棒性图正则非负矩阵分解RGNMF,判决图正则非负矩阵分解GDNMF)辨识结果相比较。为了公平，每种方法都选择差异表达程度最高的前100个基因，利用在线工具ToppFun(https://toppgene.cchmc.org/input_enrichment.jsp)对辨识得到的差异表达基因进行GO分析，ToppFun的阈值参数设置如下：p值小于等于0.01,最小的基因数设为2，具体结果见下表。从表中不难看出对于联合约束的非负矩阵分解方法辨识得到前10个基因相应的p值明显比其他三种方法小，因此对于头颈鳞癌基因表达数据集，利用联合约束的非负矩阵分解方法辨识差异表达基因更有效。

图2为本发明RGSON和经典NMF方法辨识误差收敛性,图2说明本发明收敛速度比经典NMF方法快，收敛后辨识误差比NMF方法小。图3为降维的维数对辨识误差的影响。图3说明随着降维维数k的增加辨识误差快速减小，k＝10后收敛速度变慢。

Claims (1)

1.一种基于联合约束非负矩阵分解的差异表达基因辨识方法，其特征在于包括如下步骤：

(1)将癌症数据集中的数据进行预处理，表示为非负矩阵j∈[1,n]，其中x_j表示第j个样本的表达基因构成的列向量，R₊表示正实数集合，p表示矩阵X中基因的个数，n表示矩阵X中样本的个数，将X中的所有元素都归一化到(0,1)上；

(2)对基矩阵F和系数矩阵G进行非负约束，对系数矩阵G矩阵进行图正则约束，能够保留原始基因表达数据的内在几何信息，对系数矩阵G矩阵进行正交约束，可以得到更加有效稀疏的分解结果，利用L_2,1范数约束误差函数来减小异常值和噪声的影响，构建联合约束非负矩阵分解的目标函数：s.t.F≥0,G≥0,G^TG＝I，将正交约束引入该目标函数，可得：其中，F为基矩阵，G为系数矩阵，L为图拉普拉斯矩阵，I为单位阵，|| ||_2,1表示矩阵的L_2,1范数，|| ||₁表示矩阵的L₁范数，Tr(·)表示矩阵的迹，α表示图正则约束项Tr(G^TLG)的权重，γ表示稀疏约束项的权重，β为正交约束的权重，α,β和γ均为正数；

(3)对联合约束非负矩阵分解的目标函数进行优化求解，实现对非负矩阵进行分解，得到相应的基矩阵和系数矩阵其中表示原始基因表达数据的基空间，其中f_i表示第i个样本包含所有的基因信息，表示低维实数空间R^p×K中全体非负矩阵构成的子集，p为基因的个数，K为降维的维数；为系数矩阵，表示实数空间R^K×n中全体非负矩阵构成的子集，n为样本的个数，K为分解后的矩阵维数；

(4)根据分解得到的基矩阵F，对其中每一行求绝对值的和，得到评估向量

(5)对评估向量中的元素从大到小进行排序，中的元素数值越大差异表达程度越高，从而辨识得到差异表达基因；

(6)利用基因本体分析工具GO对辨识的差异表达基因进行功能解释。