基于非负矩阵分解的大脑不同区域基因表达数据分析

孔 薇 陶伟杰 牟晓阳

1（上海海事大学信息工程学院，上海 201306）2（美国罗文大学生物化学系，新泽西 08028，USA）

引言

后基因组时代研究的重心已由对单个基因或蛋白质功能的局部性研究，转移到以细胞内全部基因、蛋白质及代谢产物等为研究对象的“组学”时代。阿尔茨海默症（Alzheimer's disease，AD）是老年痴呆症中最普遍的一种，但其早期发现难、确诊难、发病机理复杂，至今对它的生物学机制研究仍然没有本质上的突破。随着DNA微阵列技术在生物学中的广泛应用，可以得到大量暗含着基因调控关系的基因表达数据。越来越多的学者正利用微阵列基因表达数据，对AD产生和发展的机理、基因调控关系进行探寻。利用生物信息学方法进行特征基因选择、信号传导通路分析及基因调控网络重构，是目前帮助生物学家缩短实验周期、减小实验范围的最有效手段之一。

微阵列基因表达数据分析可以从单基因、多基因和基因调控网络3个层次进行。单基因水平主要依据单个基因在不同实验条件下的表达水平，利用统计学的方法推测其功能；多基因水平主要利用聚类算法，从协同调控、共同功能、相互作用等方面研究基因组的表达情况；最高层的分析是通过探寻潜在基因调控的网络，从机制上解释观察到的基因表达数据。本课题将采用改进的无监督双向聚类方法、非负矩阵分解法，探寻AD基因表达数据中的潜在信息，识别与AD密切相关的特征基因、分析信号传导通路以及基因调控关系。

基因表达数据的基因数目多、样本过少，具有高噪声、高变异且信息变量隐藏很难分析的特点，考虑到这种不平衡容易造成许多经典模式识别和机器学习方法泛化能力不强、精度不够等缺陷，本研究除了使用在AD病理改变中最明显的人脑海马区（hippocampus，HIP）组织样本，还对内嗅区皮质（entorhinal cortex，EC）、颞中回（media temporal gyrus，MTG）及视 觉皮层（primary visual cortex，VCX）等区域的组织样本进行非负矩阵的分解分析，因为这些区域在AD患者脑内也有不同程度的病理改变。利用不同脑区域数据所提取的显著基因，构建了AD的基因表达调控网络；生物学分析证实，采用不同脑区域显著基因的联合分析，对连通作用下共表达基因模块的相互作用分析更加有利。本课题的研究为进一步指导分子生物学实验，协助生物学家分析AD致病机理，提供了有益的依据和新方法。

1 NMF的基因表达数据模型

1.1 NMF的应用

非负矩阵分解（non-negative matrix factorization，NMF）技术是一种矩阵分解方法，其计算模型与主成分分析（principal component analysis，PCA）、独立成分分析（independent component analysis，ICA）、矢量量化（vector quantization，VQ）分析及因子分析（factor analysis，FA）等的计算模型类似，所不同的是对数据进行分解的限制条件为非负性。NMF分解方法相对于传统的一些算法而言，具有实现上的简便性、分解形式和结果上的可解释性以及占用存储空间少等诸多优点，目前已在图像分析［1］、人脸识别［2］、文本聚类［3］、语音信号处理［4］等领域有了广泛应用。在医学信息处理中，由于NMF算法的局部性具有双向聚类功能，特别是在基因表达数据中，它能将一个基因聚类到不同的类别中，符合同一个基因能参与不同的生物过程或信号传导通路这一实际生物学特性，近年来已经被用于处理核医学动态连续图像［5］、癌症基因表达数据样本分类［6］、脑磁共振化学位移成像、fMRI 成像分析［7－8］等方面。

1.2 基因表达数据的NMF模型

1999年，Lee等提出了 NMF及其简易算法，并证明了这种算法的收敛性［9］。NMF可表示为

式中：V为n×m的基因表达值，n为基因个数，m为样本数或不同的实验条件；矩阵W为n×k维非负基向量，Hk×m为非负系数矩阵；NMF把数据分解为k个聚类，k代表W×H特征子空间的维数，也是降维的维数或分类的类别数。基因表达数据V则为非负基矩阵W和系数矩阵H的线性组合。

NMF是利用非负性约束来获取数据表示的一种方法。非负性是对矩阵分解非常有效的限制条件，它使得原始数据可以基于“部分”的表示形式，即样本数据只允许加性的和非负的组合，而不存在减运算。算法所得到的非负基向量组W具有一定的线性无关性和稀疏性，从而使其对原始数据的特征及结构具有相当的表达能力，反映了由部分构成整体最直观的解释。

但是，传统NMF方法的性质造成了它分解不唯一的缺陷。近年来有许多改进算法，大致分为两类：第1类主要使算法能收敛到全局最优点，第2类是提高算法的稀疏度。近几年通过提高分解矩阵的稀疏性，学者们提出了许多改进方法，如局部NMF 方法（LNMF）［10］、稀疏 NMF 方法（SNMF）［11］、稀疏约束的NMF方法（NMFSC）［12－13］和非平滑NMF（non-smooth NMF，nsNMF）［14］等。

由于混合模型的乘法本质，如果要求 W和H同时具有一定稀疏度，那么必然会恶化模型对数据的拟合效果。非平滑NMF（nsNMF）算法则是引入了一个“平滑”矩阵S到分解模型中，使W和H矩阵都具有稀疏性，达到了全局的稀疏性，可表示为

式中，S∈Rk×k，是一个正平滑矩阵，定义为

式中，I 表示单位矩阵，I表示元素为1的列向量，K是一个k×k矩阵（所有元素为1），参数?θ满足0≤θ≤1，控制模型的稀疏度。

平滑矩阵S的调节作用可解释如下：假设X为一个非负向量，乘以平滑矩阵S得到Y=SX，如果θ=0，则Y=X，矩阵X为非平滑的；而随着 θ→1，Y趋向于常数向量。矩阵Y中的每一个元素趋近于X元素的平均值，这是最平滑的向量，因为向量Y中所有的数值都等于相同的非零值，而不是有一些大的数值和接近于0的离散值。进一步观察nsNMF模型，可以等效为

这样的结合方式在算法迭代过程中十分有利。当迭代H时，可用WS代替W，同理，迭代W时，用SH替换H。这样在算法收敛时，就可以同时保证矩阵W和H均具有较高的稀疏性。采用以下的散度（熵）为目标函数式中，i表示4个不同数据集的基因个数，j表示4个不同数据集的样本量。

分别对W、H和S求偏导，得到W、H和S的迭代式为

式中，I表示4个不同数据集的基因个数，u表示4个不同数据集的样本量，a表示基向量的列。

式中，k表示基向量的列。

式中，v表示基向量的列。

nsNMF算法的优势在于同时增强了基向量矩阵W和系数矩阵H的稀疏性。在分解结果中，利用矩阵H中对应样本的系数值高低，即可将m个样本进行分类。在NMF分解后，第i个原样本可表示为基向量（k个W的列）与对应系数矩阵（H的第i列）的线性组合。显然，在基向量及系数矩阵中表达值大的基因及系数对于原数据值起着重要作用。

2 仿真实验

2.1 数据准备

本课题使用的数据集来自NCBI（National Center for Biotechnology Information）GEO DataSet中的GSE5281系列。此数据集采集了161个独立个体的大脑样本，共54 675个基因探针的表达谱数据，把大脑与AD及衰老相关的脑区域分为6个数据块。由于人脑海马区担负着将短期记忆转化为长期记忆以及空间定位的重要功能，同时临床研究表明AD患者最初有Aβ淀粉样蛋白沉积及脑组织病变的区域就是海马区，所以本研究首先选取了人脑海马区（HIP）组织样本的基因表达数据。另外，考虑到内嗅皮质区（EC）在人脑定位、学习和记忆中起着重要作用，颞中回（MTG）的结构对记忆功能起着重要作用，以及视觉皮层区（VCX）大脑皮层与长期记忆有密切关系，还选取了人脑上述3个区域的基因表达数据来进行特征基因提取实验。

海马区（HIP）数据块共23个样本，包括13个正常对照样本和10个AD患病样本；内嗅区皮质区（EC）数据块与海马区相同，具有13个正常对照样本和10个AD患病样本；颞中回（MTG）包括 12个正常对照样本和16个 AD患病样本；而视觉皮层（VCX）包括 12个正常对照样本和19个AD患病样本。经过对各组数据进行t检验及微阵列显著性分析（significance analysis of microarray，SAM）方法的去噪和预处理，将 HIP、EC、MTG和 VCX数据集的基因个数分别筛减为12 000、6 200、10 000和7 600。

2.2 数据仿真实验

在NMF分解模型中，最佳样本聚类数k是给定范围内最佳因式分解的秩，也就是经NMF算法分解后所得到的特征子空间的维数，分解结果是对数据集的最佳聚类数和每个实验条件下聚类效果的最佳判断。在进行NMF分解之前，需先确定 k的值，最佳因式分解的秩对应其中最稳定的聚类，而稳定性采用共表型相关性系数来判断。定义一个矩阵C，存储每次运行程序产生的样本分类结果，如果两个样本在同一个类则记为1，否则为0；计算矩阵C的平均值得到平均值矩阵 averC，如果类别归属确定，则平均值矩阵averC的各项将接近于0或者1；用平均值矩阵averC所诱导的样本间的距离矩阵与I-averC间的Pearson相关系数去度量平均值矩阵averC的分散程度，则Pearson相关系数即为共表型相关性系数（cophenetic correlation coefficient，CCC）。用共表型相关性系数来度量k取不同值时分类的稳定性，该系数的范围从1到0，越接近1，说明分类越稳定。观察该系数随k的变化情况，则可选择到合适的数值如图1所示（每个小图的横坐标为k的取值，纵坐标为CCC值），其中，（a）～（d）分别为大脑4个不同区域传统NMF方法下不同k对应的 CCC取值，（e）～（h）则为大脑4个不同区域 nsNMF方法下不同k对应的 CCC取值。可以看出，当k取2时，CCC值均为最高。如表1所示为CCC的具体取值，就是图1所示的大脑4个不同区域基因表达数据当k分别取2～5时传统 NMF及 nsNMF方法的共表型相关性系数曲线。

表1 HIP、EC、MTG及 VCX区 NMF分解共表型相关性系数CCC值Tab.1 The CCC of the NMF result of HIP，EC，MTG and VCX

图1 大脑不同区域基因表达数据不同k值下不同方法的CCC值。（a）HIP区的传统NMF方法；（b）EC区的传统NMF方法；（c）MTG区的传统NMF方法；（d）VCX区的传统NMF方法；（e）HIP区的nsNMF方法；（f）EC区的nsNMF方法；（g）MTG区的nsNMF方法；（h）VCX区的nsNMF方法Fig.1 The CCC of gene expression data of 4 brain areas under different k by NMF and nsNMF.（a）CCC of HIP by NMF；（b）CCC of EC by NMF；（c）CCC of MTG by NMF；（d）CCC of VCX by NMF；（e）CCC of HIP by nsNMF；（f）CCC of EC by nsNMF；（g）CCC of MTG by nsNMF；（h）CCC of VCX by nsNMF

在图1中，（a）～（d）为传统 NMF在4个不同大脑区域的值，第2排（e）～（h）为 nsNMF在4个不同大脑区域的值。可以看出，在每个子图中k取2时，共表型相关性系数 CCC均为最高，同时nsNMF方法的CCC值略高于传统NMF法，因此取k=2，将nsNMF方法分别应用于4个不同大脑区域的基因表达数据，分别迭代2 000次达到收敛。通过nsNMF获得的基向量矩阵W具有较好的稀疏性，即只有小部分项具有非零值，使基因表达局部化。之后，分别根据W的列从大到小排序，将原数据矩阵V对应W不同的列进行基因重新排序，就可以得到基因及样本的局部表达特性。系数矩阵H的列与样本相对应，当某一类的样本与某些局部生物过程密切相关时，系数值将较高，反之较低，因此可以根据矩阵H的值来进行样本分类。图2为海马区（HIP）基因表达数据 nsNMF分解结果 H（2×23，2行分解系数，23个样本）的基因表达值曲线，其中横坐标为样本序号，纵坐标则为对应第1行和第2行H基因表达值的大小。可以看出，前13个样本和后10个样本的数值可以明显区分开。例如，图2中前13个样本对应的H第1行的值较大，而 H第2行的值较小；后10个样本对应的H第1行的值较小，而H第2行的值较大；分类效果较好。按照排序后V基因表达值的大小，本实验提取了海马区（HIP）与AD密切相关的900个显著基因。

图2 HIP区的nsNMF分解结果中两行H的基因表达值Fig.2 The gene expression of H by nsNMF on HIP

同理，分别得到了内嗅皮质区（EC）、颞中回（MTG）及视觉皮层（VCX）基因表达数据的nsNMF矩阵分解结果，如图3～图5所示。

图3所示为内嗅皮质区（EC）基因表达数据nsNMF分解结果 H（2×23，2行分解系数，23个样本）的表达值曲线，横坐标为样本序号，纵坐标则为对应第1和第2行H基因表达值的大小。可以看出，前13个样本和后10个样本分别对应于H第1和第2行值的分类效果很好。

同样，图4～图5分别为颞中回（MTG）及视觉皮层（VCX）基因表达数据nsNMF分解结果H的基因表达曲线。在图4中，颞中回（MTG）的前12个正常对照样本和后16个AD患病样本，以及视觉皮层（VCX）的前 12个正常对照样本和后19个AD患病样本，均能够根据H值进行较好的分类。这表明，W中所示的表达差异的基因与类别具有明显关联，这些基因与AD密切相关。分别从内嗅区皮质区（EC）、颞中回（MTG）和视觉皮层（VCX），提取了700、800和1 400个表达显著的基因。

图3 EC区的nsNMF分解结果中两行H的基因表达值Fig.3 The gene expression of H by nsNMF on EC

图4 MTG区的nsNMF分解结果中两行H的基因表达值Fig.4 The gene expression of H by nsNMF on MTG

图5 VCX区的nsNMF分解结果中两行H的基因表达值Fig.5 The gene expression of H by nsNMF on VCX

3 生物学分析及讨论

本研究利用nsNMF方法，根据算法结果的局部特性和双向聚类的效果，共提取了3 800个显著基因。对比最新的国际阿尔茨海默氏症研究的官方网站（http：//www.alzgene.org/TopResults.asp）上已给出的确定与AD致病机理有关联的42个基因，利用基因表达数据，在HIP区识别出了 CR1，在 EC区提取出了 PGBD1，在 MTG区提取出了 CLU，在EXOC3L2、MTHFR、THRA及 VCX区识别出了TNK1、SORCS1、MTHFR、THRA共10个基因。

为了建立所提显著基因与已知致病基因之间的信号传导及调控关系，将所提取的3 800个基因，联合已知的42个基因，利用 Cytoscape软件构建了AD的基因表达调控网络。由于Cytoscape数据库中不一定涵盖所有导入的显著基因数据，因此会对输入的显著基因进行筛选和剔除，基因调控网络如图6所示。

图6 大脑4个不同区域所提显著基因构件的AD基因调控网络Fig.6 The gene regulatory network reconstructed by the significant genes extracted by 4 brain areas

从nsNMF的分解结果可以看出，所提取的显著基因在细胞凋亡、生物合成、代谢过程、神经系统等相关生物过程中起着重要的作用。在显著基因中还发现，大量与炎症相关的基因在AD样本中表现为上调。许多显著基因的过表达伴随着细胞的炎症反应，导致神经元损害，引起记忆减退和认知障碍，产生痴呆症状。许多抑制细胞生长和促进细胞凋亡的基因表达均上调，如Programmed cell death 6、deleted in liver cancer 1、cyclin N-terminal domain containing 2、killer cell immunoglobulin-like receptor、P21、P53等，而那些促进细胞周期、细胞分裂和细胞修复的基因均表现为下调。同时，还提取出了与线粒体功能障碍、转录翻译、氧化磷酸化过程相关的基因簇。此外，NMF算法提取出了 ATF4、AXIN1、CALCOCO1、 CD81、 CHP、 CREBBP、 CSNK1A1、CTBP1、CTBP2、EEF1D、FGF1、FGF2、FGFR3、FKBP1A、G3BP2、GNG12、HRAS、JUND、KIT、KRAS、LITAF、MAP2K1、MAP2K4、MAP3K5、MAP4K4、MAPK9、MAPKAPK2、MAPT、MAX、NFATC4 等与 AD密切相关的信号传导通路的基因，以及MAPK信号通路、NF-κB信号传导通路、Wnt信号传导通路等。另外，还提取出了直接参与AD信号传导通路的重要基因，包括 APP、APOE、THRA、TF、BIN1、CLU、SORL1、DAPK1、NFKBIA、SNCA、BACE1、NAE1 等。

在此基础上，利用GO数据库，对每个基因的生物学过程、分子功能及细胞组分的相关信息进行采集并分类。分析结果表明，在显著基因信号传导通路，与APP有直接相关通路关系的基因有 SHC3、F12、DAB1、SNCB、CLU、IDE、LDLRAP1、APPBP2、CHRNA7、 TGFB2、 CALR、 MAPK8IP1、 NUMBL、KLC1、HSPG2、NF1、TP53BP2、TM2D1、MME、NUMB、COL1A2等，与炎症反应密切相关的有 UBE2D3、BTRC、UCP2、MRCL3、BCL3、FBXW11、TNIP2、TIMM50、TRAF6、RUVBL2、SUCLG1、IL2RA、IL8、TPR、HMGA2、TP53BP1、NCOR2、MAP3K7IP1 等。这些显著基因信号传导及调控关系的识别，为生物学实验和AD致病机理的分析提供了依据和基础。

4 结论

由于基因表达数据的高噪声、高维性、高冗余以及数据分布不均匀等特点，使得在提取其潜在生物信息过程中仍然有很多挑战性问题。基于此，笔者将一种无监督学习方法——非负矩阵分解方法，应用到基因表达数据中；考虑到算法的稀疏性和解不唯一等问题，引入平滑矩阵来改进传统NMF算法提取与AD相关的特异性基因，探寻AD的致病机理。这一方法解决了传统聚类算法只能将一个基因聚类到某一类的缺陷，而且nsNMF算法中平滑矩阵的使用同时增强了两个分解矩阵的稀疏性，有效增强了生物过程、信号传导通路的局部特性。同时，将nsNMF算法应用于 AD患者大脑中与记忆、认知及定位密切相关的海马区（HIP），内嗅皮质区（EC），颞中回（MTG）和视觉皮层区（VCX）4个不同区域，进行显著基因提取及生物过程分析，通过与42个已知AD致病基因的比对以及基因调控网络的生物学分析，可知大脑多区域数据集的联合分析能更全面地探寻疾病的信号传导关系及基因调控方式，不仅能够帮助生物学家分析AD的致病机理，有效指导分子生物学实验，也为早期诊断、治疗、预防AD以及新药物的研发提供有价值的参考。

［3］Zhao Weizhong，Ma Huifang，Fu Yanxiang.An effective nonnegative Matrix Factorization algorithm and its application in document clustering［J］.Advanced Materials Research，2011，268－270：1508－1513.

［4］Shahnaz FB，Michael W，Pauca VP，et al.Document clustering using nonnegative matrix factorization ［J］. Information Processing and Management March，2006，42（2）：373－386.

［5］Woolfe F，Gerdes M，Bello M，et al.Autofluorescence removal by non-negative matrix factorization［J］.IEEE Transactions on Image Processing，2011，20（4）：1085－1093.

［6］Devarajan K，Ebrahimi N.Class discovery via nonnegative matrix factorization［J］. American JournalofMathematicaland Management Sciences，2008，28（3－4）：457－467.

［7］Gao Yuan，George C.Improving molecular cancer class discovery through sparse non-negative matrix factorization ［J］.Bioinformatics，2005，21（21）：3970 －3975.

［8］Lee JH，Seongbuk G.Investigation of spectrally coherent restingstate networks using non-negative matrix factorization for functional MRI data［J］.International Journal of Imaging Systems and Technology，2011，21（2）：211 －222.

［9］Lee DD，Seung HS.Learning the parts of objects by nonnegative matrix factorization［J］.Nature，1999，401：788 －793.

［10］SajdaP， Du S， Brown TR， etal. Nonnegative matrix factorization for rapid recovery of constituent spectra in magnetic resonance chemical shift imaging of the brain［J］.IEEE Trans Med Imaging.2004，23（12）：1453－1465.

［11］Montcuquet AS，Hervé L，Navarro F，et al.In vivo fluorescence spectra unmixing and autofluorescence removal by sparse nonnegative matrix factorization ［J］.IEEE Trans Biomed Eng，2011，58（9）：2554－2565.

［12］Kopriva I. Single-frame multichannel blind deconvolution by nonnegative matrix factorization with sparseness constraints［J］.Opt Letter，2005，30（23）：3135 －3137.

［13］Türkan M，Guillemot C.Image prediction based on neighborembedding methods［J］.IEEE Trans Image Process，2012，21（4）：1885－1898.

［14］Pascual-Montano A，Carazo JM，Kochi K，et al.Non-smooth non-negative matrix factorization（nsNMF） ［J］. IEEE Transactions on Pattern Analysisand Machine Intelligence，2006，28：403－415.