针对不同特征基因挖掘方法的特征基因功能一致性分析

邹 晶 高 磊 李 晋 戴静珠 李 霞，*

1(首都医科大学生物医学工程学院生物信息学系，北京 100069)

2(哈尔滨医科大学生物信息科学与技术学院，哈尔滨 150086)

引言

随着基因芯片技术的发展，生物学家可以检测高通量的基因表达谱数据。在对于基因表达谱数据的分析中发现，基因数远远多于样本数，而只有少部分基因差异表达并可以此分类不同的样本类，其中多数的被检测基因与疾病状态并无关联。提取特征基因并以此推断该基因是否与某种组织类别、疾病状态或临床表现相关联，是基因表达谱分析的基础工作。用少量的特征基因代替高通量的数据分析不仅可以提高准确率，而且可以降低计算复杂度和成本。

为此，研究者提出了各种特征基因挖掘方法，大体上可将其分为三类［1］:过滤法 (filter)［2－3］、缠绕法(wrapper)［4］以及过滤法与缠绕法的混合法(hybrid)［5］。

过滤法主要针对单个基因，用某个指标对单个基因在各类样本间表达水平的差异程度进行打分，进而滤除分数较低的基因。过滤法运行速度快，对含有大量基因少量样本的表达谱尤为适宜，但是过滤法所提取的基因的样本分类准确率较低，而且对各基因独立评测损失了基因间的信息。t检验是经典的过滤法，其他还有信息增益［6］、χ2检验［7］、Relief－F［8］、马尔可夫毯(Markov Blanket)［9］和边际过滤法［10］等。

缠绕法是在特征基因的选取过程中结合样本分类，挑选部分样本作为训练集进行分类，在测试集中评价所训练的分类器的优劣。缠绕法相较于过滤法对样本分类的准确率较高，且考虑了基因之间的相互关系，但存在滤除某些特征基因的风险，且计算量也较大。常用的缠绕法有决策树、支持向量机、PMBGA(probabilistic model building genetic algorithm)等［11］。

过滤法与缠绕法的混合方法是一种结合过滤法和缠绕法优劣势的折中方法，它运用样本分类提高了精度，但损失了一部分运算速度，如基于决策树的集成决策方法［1］。

与此同时，新方法的诞生必然产生一个与其他方法所不同的特征基因集合［12－14］，这也增大了研究者对特征挖掘方法的选择难度。不论基于何种理论的特征基因挖掘方法所得的特征基因都应倾向于样本间差异表达，而功能相似的基因倾向于共表达，则各方法所挖掘的特征基因集合可能在功能上也具有一致性。生物学上，基因功能是指基因及其产物在生物过程中所发挥的作用。在基因表达谱的分析中，基因功能可以狭隘的定义为该基因表达水平差异对样本的分类准确性，即为基因的样本分类能力。基于这种情况，本研究试图分别从生物功能和基因样本分类能力的角度，评价特征基因挖掘方法所挖掘的特征基因的功能一致性，并将此方案应用到 RankGene［14］中的八种常用的特征排秩准则，如:传统的t检验、一维支持向量机等。

1 材料与方法

1.1 基因表达谱数据及其预处理

四套基因表达谱数据均于NCBI(National Center for Biotechnology Information)的GEO(Gene expression omnibus)数据库下载，其基本信息如表1所示。由于RankGene中的一维支持向量机和t检验只能处理两类样本的数据，故研究中的基因表达谱数据均只有两类样本。

表1 基因表达谱数据的基本信息(GDS2250中原有47个样本分类于4类中，此处只选取了2类样本量大的样本类)Tab.1 Information ofgeneexpression profiles(In GDS2250，there are four classed with 47 samples，only two classes with larger samples adopted)

图像损坏、芯片上的小划痕或灰尘等细微因素都可能导致表达谱数据的部分表达值缺失，在表达谱数据的聚类算法中，缺失值可能对聚类结果造成严重的影响，且很多基因表达谱数据的分析均要求完整的数据，而由于时间、花费等问题使得重做表达谱的可行性较低，缺失值的填充成了一种比较折中的方法。RankGene中各种特征基因挖掘方法其本质上为分类器，需要在预处理中进行缺失值的填充。现阶段最常用的缺失值填充方法有K近邻法、取均值、补零等方法。王栋等［19］研究了各种填充缺失值的方法对分类的影响，分析发现，在缺失率＜40%情况下，补零的方法与其他方法的效果差别不大，在其研究的数据中基于差异表达基因的样本分类准确率基本达到80%以上，故研究中对缺失值进行了补零处理。另外，基因表达谱数据均进行了标准化，以使表达水平的均值为0，标准差为1。

1.2 特征基因挖掘方法

RankGene整合了八种特征基因挖掘方法:信息增益、二分规则、少数类总和、少数类极大值、Gini指数以、方差总和、t检验以及一维支持向量机。前六种方法试图以不同的准则评测基因的样本预测力以对基因排序(单纯的根据表达水平将样本分为上调组和下调组)［20］，具体准则如表 2 所示［21］。此六种方法的评价准则中，少数类极大值只考虑左右分支中较大的纯化度，其他准则均综合左右分支的纯化度。少数类总和的准则本质上为单纯的错误分类的样本计数，其有个比较明显的缺陷(假设隶属两类的100个样本，排列方式为50个第一类，20个第二类，30个第二类，则无论如何划分其纯化度均为20)。信息增益中计算的为信息熵的含量，即为左右分支信息量与总体信息量间的差异度。二分规则、Gini指数和方差总和则采用不同方式计算左右分支的纯化度总和。t检验以基因的t检验值对基因进行降序排秩。一维支持向量机(One Dimensional Support Vector Machine，1D－SVM)的排序准则为最优支持向量的值［22］。

表2 前6种特征基因挖掘方法的评价准则，其中k为总类数，n为总的表达值数目，nl、nr分别为被分类到左右分支的数量，li和 ri为分类到左右分支的特征值，ci为第 i个样本的类标签Table 2 Measure of predictability of the first six gene mining methods(k is the total number of classes;n is the total number of expression values;nl(resp.，nr)is the number of values in the left(resp.，right)partition;li(resp.，ri)is the number of values that belong to class i in the left(resp.，right)partition;and ciis the class of the ith sample.)

对于每个特征基因，八种方法中各有一个秩(对于有其他基因同秩的情况，取原秩次的期望值作为该基因的秩)，取基于八个秩的期望值作为集成方法中该基因的秩。用新的秩次对各特征基因进行排序，从而得到了新的特征基因挖掘结果，即为集成法。

1.3 生物功能一致性分析

GO(Gene Ontology)集成了3个结构化功能类库(生物过程、细胞成分和分子功能)以分析基因及其产物在各种生物机能中的作用［23］。基因可以注释到GO功能类以分析其生物功能。研究中采用Ochiai/Otsuka系数作为功能相似性评测指标［24］，对于已知的基因集合Gi，其所注释的GO功能类集合为Fi，则两个基因集合 Gi和 Gj的功能相似性系数FSC(Gi，Gj)(function similarity coefficient)定义为:

其中 ni和 nj分别为 Fi和 Fj中功能类数目，ni，j为Fi与Fj中共有的功能类数目。

任意两个基因集合的功能一致性可以用FSC表示，FSC的取值在0和1之间，FSC值越高的两个基因集合的功能相似性越高。对于n个基因集合，两两集合间均可有一个FSC值，即可得到一个n×n的功能相似性系数矩阵(FSCM)，其中矩阵的每一行为该集合与其他各集合的生物功能相似性系数向量，且任意两个向量间的相关系数 CORR(correlation coefficient)亦可反映此两基因集合的生物功能一致程度。对任意两个基因集合Gi、Gj，其相关系数用 CORR(Gi，Gj)表示，CORR(Gi，Gj)不仅考虑到 FSC(Gi，Gj)的大小，同时考虑了 Gi、Gj与其他基因集合的FSC相关程度。相较于FSC的直接相似性而言，CORR是一种间接相似性，它全局地考虑了基因集合与其他集合的相似性。如基因SETDB1和BMP2间无相同的GO生物过程功能节点，即FSC=0，但其分别与另一个癌症相关的基因 RUNX3［24］有相同的GO节点，即 CORR≠0。有研究表明，基因 SETDB1 和 BMP2 也与癌症相关［25－26］，CORR 即可反映此间接功能相似性。这种间接相似性可以清晰反映基因集合间的生物功能一致程度，可采用CORR作为聚类距离，以生物功能相似性矩阵FSCM为特征，对基因集合进行聚类分析，进一步分析基因集合间的生物功能一致性程度。本研究中将各方法挖掘的特征基因注释到GO的生物过程体系，进而分析各方法挖掘的特征基因集合的生物功能一致性。

1.4 样本分类能力分析

基于所提取的n个特征基因，用K均值聚类的方法，对样本进行聚类分析，将聚类结果与原始类标签进行比较，以分析各方法的样本分类能力。考虑到排序方法与聚类距离可能存在特异性，故分别采用了欧式距离、相关系数和余弦3种距离。

研究中采用了约当指数(Rand index，RI)作为样本分类能力的评价指标，RI的定义如下

RI的取值为［0，1］，当 RI值为 1时，样本分类准确度为100%，即聚类结果与原数据类完全一致;反之，当RI值为0时，样本分类准确度为0%，即聚类结果与原数据类无相同类信息。Rand index值越大，该方法所挖掘的特征基因的分类能力越好。

2 结果分析及讨论

2.1 生物功能一致性分析

将各方法挖掘的特征基因注释到GO的生物过程功能类中，计算任意两种方法所得的特征基因集合的功能相似系数FSC。特征基因数目可能直接影响研究结果，实际研究中往往根据不同研究重点选择合适的特征基因数。研究中的特征基因数n取10到600，步长为10的梯度值，不同数目的特征基因数时集合间的FSC并不一致，但其FSC值变化较小(变异系数均小于0.05)。对特征基因数n为100时的FSCM进行层次聚类分析(聚类距离为相关系数)，其结果如图1所示。

图1 不同数据中，特征基因数n为100时的FSCM层次聚类结果(IG:信息增益，TR:二份规则，MS:少数类总和，MM:少数类极大值，GI:Gini指数，SV:方差总和，TT:t检验，1D:一维支持向量机，IN:集成法;深灰色区域具有高的FSC值，而浅灰色区域的FSC值较低)(a)GDS90;(b)GDS807;(c)GDS1949;(d)GDS2250Fig.1 Hierarchical cluster result with FSCM in different datasets while feature number equals to 100(IG:information gain，TR:twoing rule，MS:sum minority，MM:max minority，GI:Gini index，SV:sum of variances，TT:t-statistics，1D:1D-SVM，IN:integrated method;Darker areas are with higher FSC while lighter ones with lower FSC).(a)GDS90;(b)GDS807;(c)GDS1949;(d)GDS2250

分析图1可知，4套表达谱数据中，各方法所挖掘的特征基因集合间SFC值均在0.5以上(多数不低于0.7)，其中Gini指数、方差总和和二分规则所挖掘的特征基因生物功能一致性最好(FSC均在0.85以上)，且该三种方法与信息增益的一致性较好;t检验与一维支持向量机的生物功能一致性程度相对较高，而少数类极大值和集成法与各方法间的一致性结果在各数据间的变化较大，不可明确比较其与何种方法的一致性更高。

集成法为各方法间的均衡方法，当某些方法间的生物功能有差异时，其很难和某种方法保持稳定的一致性。少数类极大值考虑的为左右分支的最大的纯化度，与其他考虑左右分支纯化度的方法相比，其稳定性较差，这可能也是其在各数据间变化较大的原因。

图2 GDS90、GDS807及GDS1949中不同特征基因挖掘方法分别在欧氏距离、相关系数、余弦距离聚类距离下的RI曲线。(a)～(c)GDS90;(d)～(f)GDS807;(g)～(i)GDS1949Fig.2 RI curves of different feature selecting methods in GDS90，GDS807 and GDS1949 with Euclidean，correlation and cosine distance seperately.(a)～(c)GDS90;(d)～(f)GDS 807;(g)～(i)GDS1949

在功能相似系数FSC的定义中，简单的将功能类节点的关系分为0和1两种，即完全相同的功能类间相似系数为1，不同的为0，如果考虑不同功能类间的关系(如节点间的距离、节点间最近父节点的深度等)，则功能类节点间相似系数范围为［0，1］，两个基因集合的功能类相似系数较现有测度将有所增加。在相对简单的功能相似性系数FSC的评测下，各方法所得的特征基因集合间的功能相似性系数FSC均较高，如考虑更深层次的功能类间的关系，可以更精确地反映基因集合功能相似性，也能得到更高的生物功能相似性评测结果。

2.2 样本分类能力的一致性分析

2.2.1 四组数据的聚类结果

对于每个基因表达谱数据，选取不同的特征基因n，分别用三种距离进行K均值聚类，进而得到相应的RI值。对每种聚类距离，分别以RI为纵坐标，特征基因数n为横坐标作图，如图2所示。其中，中空圈线为随机500次选取n个基因作为特征基因进行RI分析。对于数据 GDS2250，不论采用何种距离，n取10到600间何值，聚类结果和原数据类完全一致，即RI值为1。

由图2可知RI值随着n值增大而趋于稳定，但不能明确分辨各方法RI值的高低;无论采用何种聚类距离，二分规则、Gini指数和方差总和的 RI曲线基本重合，且除GDS807数据的欧式距离外，各数据结果中各方法的RI值趋于稳定，其稳定值如表3所示。另外，随机情况下，RI值随着 n值增大而缓慢增大，且各方法的RI值均明显高于随机情况。

GDS807数据中的RI值相对不稳定(特别是欧氏距离)，但大体上1D－SVM的RI值最小，其他各方法间的差异较小。而对于其他三组表达谱数据，各方法的RI值基本一致。

表3 各数据在不同特征提取方法及不同聚类距离下的RI稳定值Tab.3 Stable value of Rand index with different feature selecting methods and cluster distances

2.2.2 样本分类能力分析

基因表达谱中，每个基因均包含一定的样本分类信息，且信息含量高的基因倾向于差异表达，特征基因挖掘的作用就是挖掘尽可能小的特征基因集合以包含尽可能多信息量，故信息量较高的基因倾向于被挑选为特征基因。基因数越多，信息量越多，该基因集合的分类效果越好，即RI值越大。但表达谱中的总信息量是有限的，基因间内也存在冗余信息，当特征基因数增大到一定值时，该基因集合基本包含了所有该挖掘方法所能挖掘到的所有信息，故RI曲线在特征基因数n大于一定值后趋于稳定，而随机情况下的RI值持续缓慢增大。基于此理论，即可分析比较各方法所能挖掘的总信息量以分析各方法所挖掘特征基因的样本分类能力，即分析比较稳定的RI值。

分析不同距离下的RI曲线发现，相关系数和余弦距离的一致度较高，且与欧式距离的结果有一定差异，欧式距离下的 RI值低于其他两种距离，且在GDS807中不稳定。与其他两种形式的距离相比，欧氏距离较易受到特异值(过大或过小)的影响，表达谱中的噪点对欧氏距离的影响较大，可能直接影响结果的稳定性，故研究中重点分析了相关系数和余弦距离的结果。

相较其他三套数据，GDS807数据有较大的样本量(60个样本)，GDS90、GDS949和GDS2250的样本数分别为20，26，28，这可能直接导致 GDS807数据结果相对不稳定。但大体上GDS807中的一维支持向量机的分类能力相对较低，其他方法间的差异不大。对于其他三套数据，各方法所挖掘的特征基因均完全正确的对样本进行分类(除GDS90中的一维支持向量机)。当各方法挖掘的特征基因的样本分类能力均较高时，集成各方法优劣的集成法很难提高其分类准确度，故其RI值并不高于其他方法。综上分析，此九种方法所挖掘特征基因的分类能力明显优于随机情况，且各方法间的样本分类能力的一致性较高，但不能明确区分其样本分类能力的优劣。

t检验与一维支持向量机仅仅针对两类样本的特征基因提取，如用其处理多样本数据，则需要将多类样本两两类分化，此时可采用前六种方法进行多样本的特征基因挖掘。GDS807数据结果中，信息增益和t检验的RI值相对稍高，而在其他数据结果中基本无差异，在表达谱分析中可优先考虑采用信息增益和t检验进行特征基因挖掘。

3 结论

生物功能分析中，各方法所挖掘特征基因集合间的功能相似系数FSC随特征基因数的变化不大，且各集合间的生物功能一致性较高。样本分类能力分析中，各特征基因集合的RI值随着特征基因数的增加趋于稳定，且各集合的样本分类能力的一致性较高。无论是生物功能分析还是样本分类能力分析，Gini指数、方差总和和二分规则的一致性最好。个别方法所挖掘的特征基因集合的生物功能一致性较差，但是不能明确区分九种特征基因挖掘方法的样本分类能力的优劣。

研究中分析了RankGene中的8种方法及基于该8种方法的集成法，也可将此分析比较方法推广到更多的特征基因挖掘方法，如χ2检验、SAM等，同时也可以结合多种聚类和分类方法(支持向量机、贝叶斯、K近邻等等)综合评价特征基因挖掘方法所挖掘的基因集合的功能一致性。

［1］Li Xia，Rao Shaoqi，Wang Yadong，et al.Gene mining:a novel and powerful ensemble decision approach to hunting for disease genes using microarray expression profiling［J］.Nucleic Acids Research，2004，32(9):2685－2694.

［2］Mills JC，Gordon JI.A new approach for filtering noise from high－density oligonucleotide microarray datasets ［J］.Nucleic Acids Res，2001，29(15):E72－2.

［3］Xu Wenlong，Wang Minghui，Zhang Xianghua，et al.SDED:A novel filter method for cancer－related gene selection ［J］.Bioinformation，2008，2(7):301–303.

［4］Kohavi R，John GH.Wrappers for feature subset selection［J］.Artificial Intelligence，1997，97:273－324.

［5］Xing EP，Jordan MI，Karp RM.Feature selection for highdimensional genomic microarray data［A］.In:Proceedings of the Eighteenth Internatlonal Conference on Machine Learning［C］.Massachusetts:Morgan Kaufmann，2001.601－608.

［6］Golub TR，Slonim DK，Tamayo P，etal. Molecular classification of cancer:class discovery and class prediction by gene expression monitoring［J］.Science，1999，286(5439):531–537.

［7］Liu Huan，Setiono R.Chi2:feature selection and discretization of numeric attributes［A］.In:Proceedings of the Seventh International Conference on Tools with Artificial Intelligence［C］.Herndon:IEEE Computer Society，1995.388– 391.

［8］Kononenko I.Estimating attributes:analysis and extensions of relief［A］.In:Proceedings of the 7th European Conference on Machin Learning［C］.Berlin:Springer Verlag，1994.171－182.

［9］Hall M.Correlation－based feature selection for machine learning［D］.Hamilton:University of Waikato，1998.

［10］Blum A，Langley P.Selection of relevant features and examples in machie learning［J］.Artificial Intelligence，1997，97:245 －271.

［11］Paul T，Iba H.Gene selection for classification of cancers using probabilistic model building genetic algorithm ［J］.Biosystems，2005，82(3):208－225.

［12］Wei Guan，Alexander G，Sham N，et al.Discovering ovarian cancer biomarkers using Gene Ontology based microarray analysis［A］.In:Proceedings of the Seventh International Workshop on Data Mining in Bioinformatics［C］.San Jose:KDD 2007 Workshops，2007.78 －87.

［13］Draminski M，Rada－Iglesias A，Enroth S，et al.Monte Carlo feature selection for supervised classification ［J］.Bioinformatics，2008，24(1):110－117.

［14］Su Yang，MuraliTM，Pavlovic V，etal. RankGene:identification of diagnostic genes based on expression data［J］.Bioinformatics，2003，19(12):1578－15799.

［15］Perou C，Jeffrey S，Van De Rijn M，et al.Distinctive gene expression patterns in human mammary epithelial cells and breast cancers［J］.Proc Natl Acad Sci USA，1999，96(16):9212 －9217.

［16］Ma Xiaojun，Wang Zuncai，Ryan PD，et al.A two－gene expression ratio predictsclinicaloutcome in breastcancer patients treated with tamoxifen ［J］.Cancer Cell，2004，5(6):607－616.

［17］Gilks CB，Vanderhyden BC，Zhu S，et al.Distinction between serous tumors of low malignant potential and serous carcinomas based on global mRNA expression profiling［J］.Gynecologic Oncology，2005，96(3):684－694.

［18］Richardson AL，Wang ZC，De Nicolo A，et al.X chromosomal abnormalities in basal－like human breast cancer ［J］.Cancer cell，2006，9(2):121－132.

［19］Wang Dong，Lv Yingli，Li Xia，et al.Effects of replacing the unreliablecDNA microarray measurements on the disease classification based on gene expression profiles and functional modules［J］，Bioinformatics，2006，22(23):2883－2889.

［20］Murthy SK，Kasif S，Salzberg S.A system for induction of oblique decision trees［J］.J Artif Intell Res，1994，2:1 –32.

［21］Su Yang，Murali TM，Pavlovic V，et al.Rankgene:a program to rank genesfrom expression data ［EB/OL］. http://genomics10.bu.edu/yangsu/rankgene/，2002－11－18/2009－9－17.

［22］Su Yang，Murali TM，Pavlovic V，et al.Training support vector machines in1D ［EB/OL］，http://genomics10.bu.edu/yangsu/rankgene/oned－svm.pdf，2002－9 －8/2009－9 －17.

［23］Ashburner M，Ball CA，Blake JA，et al.Gene ontology:tool forthe unification of biology gene ontology ［J］.Nature Genetics，2000，25(1):25－29.

［24］Rieck K，Laskov P，Sonnenburg S.Computation of similarity measures for sequential data using generalized suffix trees［J］.The Journal of Machine Learning Research，2008，9:23－48.

［25］Kim TY，Lee HJ，Hwang KS，et al.Methylation of RUNX3 in various types of human cancers and premalignant stages of gastric carcinoma ［J］. LaboratoryInvestigation，2004，84:479–484.

［26］Li Hongwei，Rauch T，Chen Zhaoxia，et al.The histone methyltransferase SETDB1 and the DNA methyltransferase DNMT3A interact directly and localize to promoters silenced in cancer cells［J］，J Biol Chem，2006，281:19489－19500.

［27］Reinholz MM，Iturria SJ，Ingle JN，et al.Differential gene expression of TGF－beta family members and osteopontin in breast tumor tissue:analysis by real－time quantitative PCR［J］.Breast Cancer Res Treat，2002，74(3):255 －269.