杨华兰，庞朝阳，董体智，胡本琼，4

(1．四川师范大学数学与软件科学学院，四川成都610066; 2．四川师范大学计算机科学学院，四川成都610066; 3．东华软件股份公司，陕西西安710068; 4．成都理工大学管理科学学院，四川成都610059)

基于PCA的模糊C均值聚类算法识别AD候选致病基因

杨华兰1，庞朝阳2*，董体智3，胡本琼3，4

(1．四川师范大学数学与软件科学学院，四川成都610066; 2．四川师范大学计算机科学学院，四川成都610066; 3．东华软件股份公司，陕西西安710068; 4．成都理工大学管理科学学院，四川成都610059)

研究表明阿尔茨海默病(AD)的致病机理可能与基因有关．利用计算方法对AD基因表达数据进行挖掘，以获得AD候选致病基因，寻找治愈AD方法．结合生物信息理论应用基于主成分分析(PCA)方法的模糊C均值算法处理基因表达数据:观察到AD基因表达数据具有线性相关性后，先用PCA对数据降维，再利用一维分类方法对降维后的数据聚类，然后将结果提供给模糊C均值算法作为其初始聚类数目和聚类中心．通过算法，最终识别出9个AD候选致病基因．

基因表达数据;AD候选致病基因;模糊C均值算法;主成分分析

AD是一种危害性极大的痴呆症．目前为止，科研工作者已成功识别出3种不同的早发性AD致病基因:APP基因、PSEN1基因、PSEN2基因及一种晚发性AD致病风险因子APOE基因［1－3］．对致病基因进行研究后研发的相关药物，对治疗AD有一定的效果．研究表明AD可能为多基因遗传病，为彻底根治AD，扩充和健全AD致病基因系统显得必要．

由于已有的AD基因表达数据量大、维数高、噪音多．很难直接从中获取有用的生物信息．为此，本文利用计算方法对数据聚类．将与AD致病基因具有相同或相似功能的基因聚为一类，同一类中的基因表达模式类似称作共表达基因［4－5］．若一些基因在AD病发历程中(AD病发历程共包括正常阶段、轻度、中度、重度4个阶段)一直与致病基因同属一类，那么这样的基因就被识别出来作为AD候选致病基因．

1 预备知识

1．1 AD基因表达数据的组织形式 本文所使用的AD基因表达数据下载于 NCBI的GEO数据库［6］．原始数据是通过对9个正常、7个轻度、8个中度、7个重度阶段 AD样本(患者)中相同的22 283个基因实验得到．原始的基因表达数据组成4个不同的矩阵．(1)式显示的是正常阶段下基因表达数据组成的矩阵(n=22 283，p=9)．

式中，矩阵的行表示同一基因在p个样本中的表达．列是同一样本中n个基因的表达情况．类似，轻度阶段的基因表达数据构成22 283×7的矩阵，中度阶段的数据构成22 283×8的矩阵，重度阶段的基因表达数据构成22 283×7的矩阵．

1．2 AD基因表达数据预处理 由于基因表达数据矩阵中的每列数据是在不同样本中获得．为避免量纲对实验结果的影响，采用标准化处理数据的方式［7］预处理数据．记预处理后的矩阵为Y．

这里

2AD基因表达数据的处理

2．1 AD基因表达数据的特征 基因表达数据矩阵中，不同列的数据来源于不同的样本．但由于样本所处的状态相同，且矩阵的行数据是同一个基因的表达情况．因此，认为基因表达数据矩阵的不同列间的数据也应该存在某种相关性．事实上，这种相关性确实存在．

以正常阶段下基因表达数据矩阵(1)式X为例．将X中的每列数据都看作是对应一个变量．这样，X中任意2列数据就对应2个不同的变量．将这两个变量分别记为x，y．那么同一基因的2个变量所对应的基因表达数据就构成空间中的二维点(x，y)．图1［8］显示了X中第一、二列数据构成的二维点的分布情况．由图1可以看到所有的点都落在一条直线的周围，这意味着同一阶段下相同基因在不同的样本中的表达数据是具有线性相关性的．

为使这种相关性得以更清楚的表达，计算并得出了X的相关系数矩阵T(见表1)．

表1 正常阶段下基因表达数据的相关系数值［8］Table 1 The correlation coefficients of gene expression data in the normal stage of AD

将上述数值构成相关系数矩阵T=(tij)9×9，T是通过将X中的列数据分别记为向量x1，…，x9后，计算

所得到．(xi×xj表示向量内积，‖xj‖表示向量的范数)．

表1表明X中任意两列数据间相关系数大于0．85．这进一步验证了同一阶段中相同基因在不同的样本中的表达数据是具有线性相关性．这种相关性的存在使得对数据降维有意义．

2．2 应用主成分分析方法对AD基因表达数据降维 AD基因表达数据量大，维数高．而列数据间线性相关性的存在使得对基因表达数据降维有意义．主分量分析(PCA)方法是常用的降维方法，该方法用少数新变量来解释原变量，新变量是原变量的线性组合且携带了原变量的多于85%的信息．

PCA对AD基因表达数据降维处理的具体步骤如下(以正常阶段数据为例)．

步骤一:计算公式(2)中Y的相关系数矩阵R．

这里

cov(i，j)是xi与xj的协方差是第j列的平均值，Sj是第j列的标准差．

步骤二:计算R的特征值和特征向量．相关系数矩阵R是p行p列的实对称矩阵．设它的p个特征值为(λ1，λ2，…，λp)．满足λ1＞λ2＞…＞λp．特征值所对应正交化的特征向量为

步骤三:计算贡献率及累计贡献率，确定主成分的个数．

贡献率:

按照累计贡献率大于85%的原则，选择前k个主成分．

步骤四:计算主成分的载荷和得分．载荷矩阵为p行p列的矩阵Z．

元素

得分矩阵C=X×Z，这里Z是p行p列的矩阵．按照上述步骤，利用PCA对正常、轻度、中度、重度4种不同阶段的AD基因表达数据进行处理后，发现数据呈现出2个特征(以正常阶段为例)．

特征一:第一主成分的贡献率远大于其他主成分的贡献率，如图2所示．这表明第一主成分携带多于85%的原始生物信息．因而后续思考或者计算用第一主成分来代替原始数据显得合理和有意义．

特征二:第一主成分得分数据的分布具有明显的分类特征且类与类间有比较明显的边界，如图3

3 应用基于PCA的模糊C均值算法对基因表达数据聚类

识别已知的AD致病基因的共表达基因是有生物学意义的，这是因为共表达基因表达模式类似，且在功能上具有相似性或相关性［4－5］．对基因表达数据聚类可以将具有共调控关系和功能上彼此关联的基因聚在一起．截止目前，一些经典的聚类算法比如 K－Means聚类算法［8］、层次聚类算法［9］等均被应用到 AD候选致病基因的识别过程中．

虽然在AD候选致病基因识别过程中应用这些聚类算法是有意义的，但文献［4］指出生命活动中，一些蛋白在功能上具有多样性，而其功能的发挥是通过与不同蛋白协调工作达到的．不同条件下，编码这些蛋白的基因与多组不同的基因共表达．这一生物特点的存在会使得聚类过程中出现各个类别重叠的情况．识别与多组不同基因共表达(表达模式相似)的基因用一般非模糊聚类算法是很难达到的．尤其是当数据是由不同实验条件下的数据合并而成时，这种缺陷更加明显．而模糊聚类算法的出现却较好地解决了该问题．模糊聚类算法借助隶属度概念对基因进行分类，这种聚类方法相比其他聚类方法柔性更好，更适合对基因进行聚类．模糊聚类算法中最经典的为模糊C均值算法—FCM算法．

3．1 模糊C均值聚类算法简介 FCM算法作为传统的硬划分(HCM)的改进，是一种模拟人类思维的方法，该算法融合模糊集中的隶属度概念［10－11］．相对HCM算法来说FCM算法本身具有一定的柔性，理论上更适用于生物计算．FCM算法利用隶属度uij判别基因的所属类别，uij表示第j个基因属于第i个类的可能性．给定一个相对合理的阈值，如果隶属度uik大于阈值则将第k个基因划分到第i个类．FCM算法的最优聚类在目标函数J(U，V)取最小值时得到．

(4)式中，模糊指数m(m＞1)，U是c×n的隶属度矩阵，uij是第j个样本隶属于第i类的隶属度值(i= 1，…，c;j=1，…，n)，V是s×c的聚类中心矩阵，dij为样本点xj到中心点vi的距离(dij=|xj－xi|)，目标函数是关于自变量(U，V)的优化约束问题．利用KT极值的必要条件有(6)和(7)式所显示的迭代方程的存在:

记

利用FCM算法对数据聚类，需要先初始化聚类中心或隶属度矩阵，然后利用(6)和(7)式不断更新聚类中心和隶属度矩阵，直到满足设置的终止条件，算法具体步骤如下:

Step 1:设置模糊指数m，收敛值ε＞0，迭代次数k=0，初始化聚类中心V(0);

Step 2:利用隶属矩阵的迭代公式计算u(0);

Step 3:计算目标函数J(0)，k=k+1;

……

Step k:更新聚类中心得到V(m);

Step k+1:更新隶属度矩阵得到U(m);

Step k+2:计算目标函数J(m)．若

则算法停止．输出聚类中心和隶属度矩阵．否则，返回step k，算法继续．

FCM算法是一种具有柔性的模糊划分方式，在生物学基因识别方面应用合理．但是，该算法的缺点在于要求事先确定初始聚类中心．而一般情况下，初始聚类中心的选择是随机的，具有不确定性．这会给最后的实验结果带来误差．为了解决该问题，我们对FCM算法做了一定的改进:引用文献［8］提出的一维分类方法来确定FCM算法的初始聚类中心和聚类数目．

3．2 一维分类算法确定FCM初始聚类中心和聚类数目 主成分分析方法处理数据后，第一主成分被选为主要的研究对象．观察第一主成分的得分数据发现它们具有明显的分类特征(图3)．为确定类的个数．可对这些数据进行聚类，聚类遵循的原则:类间样本点间距小，类类样本点间距大．基于此，考虑到了第一主成分得分数据的间距且获得它们由大到小排列的分布直方图(图4)．观察图像，间距由大到小的变化过程可被看做是类类间距向类间间距的过度过程，而判据应该集中在曲线与直线的交点附近．为找到该交点，确定类的数目．我们引用文献［9］提出的一维分类算法．该方法具体实验步骤如下:

Step 1:设第一主成分得分集为X*．

Step 2:对X*升序排列，记排列后的得分集为Y，索引号为ID．

Step 3:对Y做向前差分，得到前后相邻2个基因间的距离矩阵Y1．

Step 4:计算所有距离的平均值，筛选出大于平均值的距离，将其个数记做l．将它们组合成新的矩阵Y2．

Step 5:定义并计算门限值Lepso．

Step 6:利用Lepso对X*分类．设n1=1，若di＜Lepso，则nj+1=nj．否则，

便得到初始分类矩阵N．

依据索引号ID={idx1，idx2，…，idxn}，将N中索引号对应到X*，取出相应的分类集合N*．表示第i个类别，k表示类别数目．

Step 7:计算N*中每个类的类别中心n(0)，

应用基于PCA的模糊C均值算法对AD基因表达数据聚类，由如下步骤得到聚类结果．

3．3 基于PCA的模糊C均值算法对AD基因表达数据聚类的具体步骤 以正常阶段的基因表达数据为例，其他3种情况类似．

Step 1:设正常阶段下的数据矩阵为X，预处理后得到矩阵Y．

Step 2:对Y作PCA降维处理．降维后的第一主成分得分集为

Step 3:用一维分类法对X*聚类(详细过程见3．2)，得到聚类结果

上述方法确定出来的n(0)和类别数k，作为FCM算法的初始聚类中心和聚类数目．由于FCM算法中的模糊指数m控制算法的柔性．目前对m的最佳取值，并没有统一的规定．研究者在对m进行选择时，通常都是结合各自的研究目的和数据特点给出m的取值或取值范围:文献［12］给出［1，1．5］，聚类效果最佳;而文献［13］又从物理学角度出发解释m=2，聚类结果最优;而文献［14］从汉字识别的应用背景出发得到［1．25，1．75］聚类效果最好;随后，文献［15］再一次从收敛角度出发提出m的取值与

Step 4:用FCM算法对基因进行聚类(详细过程见3．1)，得到聚类结果

4AD候选致病基因的判定依据

由于共表达基因分享同一表达模式，属于同一类．这种类别的属性是不随外界条件改变而改变的．因此，本文中那些与致病基因同属一类并在AD发展进程中一直与致病基因属于同一类的基因便为AD候选致病基因．换言之，若基因H0与致病基因ZAD同属一类，H0被判定为AD候选致病基因的依据为

H0∈F(ZAD)=

Fc(ZAD)∩Fi(ZAD)∩Fm(ZAD)∩Fs(ZAD)，F(ZAD)表示 AD候选致病基因集合．Fc(ZAD)、Fi(ZAD)、Fm(ZAD)、Fs(ZAD)分别表示正常、轻度、中度、重度4种不同阶段下与致病基因ZAD同属一类的基因的集合．

5 实验结果

通过基于PCA的模糊C均值算法对AD基因表达数据聚类及AD候选致病基因的判定依据．本文最终筛选出9个符合条件的AD候选致病基因(表2)．

表2 识别出的候选致病基因列表Table 2 The list of identified candidate genes of AD

识别出的9个候选致病基因中，ap2b1、Rere、PUF60、TBCB作为APOE的共表达基因在AD进程中一直与APOE同属一类．而eif4g2、rpl24、Sepw1、trim2、loc12729作为APP的共表达基因在AD病发四个阶段一直与APP属于同一类．由于共表达基因功能上的相关性，因此认为识别出的候选致病基因在AD整个过程中与致病基因一样或许发挥了相同或者相似的功能．因而对候选基因进行生物学验证具有意义．

6 总结和讨论

本文从已知的AD候选致病基因出发，对已有的基因表达数据进行分析．企图获取在AD病发整个历程中均与已知的AD致病基因同属一类的基因．另外，由于基因功能上的多面性，文章应用基于PCA的模糊C均值聚类算法对基因聚类:观察到AD基因表达数据具有线性相关性后，先用主成分分析(PCA)方法对数据降维，选择蕴含大部分原始数据信息的第一主成分作为主要的研究数据，将第一主成分投影到一维子空间发现第一主成分的得分数据具有明显的分类特征．为此，利用文献［8］提出的一维分类算法对数据进行聚类，并将结果提供给FCM算法作为其初始聚类数目和聚类中心．

应用算法，文章最终识别出9个AD候选致病基因分别为 ap2b1、Rere、PUF60、TBCB、eif4g2、rpl24、Sepw1、trim2、loc12729．

由于候选基因与已知致病基因功能上的相关性，后续工作可以对识别出的AD候选致病基因进一步分析或验证．

［1］LARS B，ROUDOLPH E T．Thirty years of Alzheimer’s disease genetics:the implications of systematic meta－analyses［J］．Nature Reviews Neuroscience，2008，9(10):768－778．

［2］HAROLD D，ABRAHAM R．Genome－wide association study identifies variants at CLU and PICALM associated withAlzheimer’s disease［J］．Nature Genetics，2009，41(10):1088－1093．

［3］LAMBERT J C，HEATH S，EVEN G，et al．Genome－wide association study identifies variants at CLU and CR1 associated with Alzheimer’s disease［J］．Nature Genetics，2009，41(10):1094－1099．

［4］岳峰，孙亮．基因表达数据的聚类分析研究进展［J］．自动化学报，2008，34(2):113－120

［5］EISEN M B，SPELLMAN P T，BROWN P O，et al．Cluster analysis and display of genome－wide expression patterns［J］．Proceed National AcADemy Sci United States Am，1998，95(25):14863－14868．

［6］美国国家生物信息技术中心．阿尔茨海默病基因数据［EB/OL］．［2013－9］Available:http://www．ncbi．nlm．nih．gov/ gds/?term=alzheimer(2013－10－10)

［7］PANG C Y，HU W，HU B Q，et al．A special local clustering algorithm for identifying the genes associated with Alzheimer’s disease［J］．IEEE Trans Nanobioscience，2010，9(1):44－50．

［8］PANG C Y，YANG L，ZHANG D X，et al．The strong correlation of gene expression data on Alzheimer’s disease and co－regulation of gene［C］//Taiwang:2011 IEEE International Conference on Granular Computing，2011:855－858．

［9］PANG C Y，LIU S Q，LI Y，et al．The nonlinear correlation character of gene expression data on Alzheimer’s disease and hierarchy clustering of co－regulated gene［C］//Taiwang:2011 IEEE International Conference on Granular Computing，2011:859－862．

［10］莫智文，舒兰，许彪．模糊数学理论及其应用评述［J］．四川师范大学学报(自然科学版)，1998，21(3):330－335

［11］庄刘，曾艳．基于模糊C－均值聚类的最优量化器设计［J］．四川师范大学学报(自然科学版)，2010，33(4):559－562．

［12］BEZDEK J C．Pattern Recognition with Fuzzy Objective Function Algorithms［M］．New York:Plenum Press，1981．

［13］BEZDEK J C．A physical interpretation of fuzzy ISODATA［J］．IEEE Trans Syst Man Cybern，1976，6(3):387－390．

［14］CHENG Y S，CHAN K P．Modified fuzzy ISODATA for the classification of handwriting chinese characters［C］//Singapore:Proc Int Conf Chinese Comput，1986:361－364．

［15］BEZDEK J C，HATHAWAY R．Convergence theory for fuzzy c－means:counterexamples and repairs［J］．IEEE Trans Syst Man Cybern，1987，17(5):873－877．

［16］FUKUNAGA K，HOSTETLER L．The Estimation of the G radient of a Density Function with Application in Pattern Recognition［J］．IEEE Transactions on Information Theory，1975，21(1):32－40．

Applying FCM Algorithm Based on PCA to Identify the Candidate Genes of Alzheimer’s Disease

YANG Hualan1，PANG Chaoyang2，DONG Tizhi3，HU Benqiong3，4

(1．College of Mathematics and Software Science，Sichuan Normal University，Chengdu 610066，Sichuan; 2．College of Computer Science，Sichuan Normal University，Chengdu 610066，Sichuan; 3．Donghua Software company，Xi’an 710068，Shaanxi; 4．College of ManagementScience，Chengdu University of Technology，Chengdu 610059，Sichuan)

Researches show that the pathogenesis of Alzheimer’s disease(AD)may be associated with genes．It’s significant to identify the candidate genes of AD and find the way of curing AD at last by mining the AD gene expression data．Combing the theory of biological information，the Fuzzy C-Means(FCM)clustering algorithm based on Principal Component Analysis(PCA)is chosen to process the gene expression data．After observing the fact that there is a linear correlation between AD gene expression data，we use the way of PCA to reduce the dimensions of data．Then a dimensional classification algorithm is put forward and is applied to cluster the AD gene expression data．Simultaneously，the clustering results are provided to FCM as its initial clustering centers and numbers．Finally，9 candidate genes of AD are identified．

gene expression;the candidate genes of AD;Fuzzy C-Means clustering;principal component analysis

O24

A

1001－8395(2016)04－0496－07

10．3969/j．issn．1001－8395．2016．04．006

(编辑 陶志宁)

2014－03－29

中国航空科学基金(2012ZD11)

*通信作者简介:庞朝阳(1973—)，男，教授，主要从事基因计算与量子力学的研究，E－mail:cypang402@126．com

2010 MSC:62H25;62H30