基因表达数据分析

第8章基因表达数据分析

基因芯片或DNA微阵列等高通量检测技术的发展，可以从全基因组水平定量或定性检测基因转录产物mRNA，获取基因表达的信息。由于生物体中的细胞种类繁多，同时基因表达具有时空特异性，因此，基因表达数据要比基因组数据更为复杂、数据量更大、数据的增长速度更快。基因表达数据中蕴含着基因调控的规律，可以反映细胞当前的生理状态，例如（？？）是否恶化、（？？）是否对药物有效等。对基因表达数据的分析是生物信息学的重大挑战之一，也是DNA微阵列能够推广应用的关键环节之一。

基因表达数据分析的对象是在不同条件下，全部或部分基因的表达数据所构成的数据矩阵。通过对数据矩阵的分析，回答一些生物学问题，例如，基因的功能是什么？在不同条件或不同细胞类型中，哪些基因的表达存在差异？在特定的条件下，哪些基因的表达发生了显著改变，这些基因受到哪些基因的调节，或者调控哪些其它的基因？哪些基因的表达是条件特异性的，根据它们的行为可以判断细胞的状态（正常或癌变）？？？？等等。对这些问题的回答，结合其他生物学知识和数据有助于阐明基因的调控路径和基因之间的调控网络。揭示基因调控路径和网络是生物学和生物信息学共同关注的目标，是系统生物学(Systems Biology，在附录中增加解释条目！)研究的核心内容。目前，对基因表达数据的分析主要是在三个逐渐复杂的层次上进行：1、分析单个基因的表达水平，根据在不同实验条件下，该基因表达水平的变化，来判断它的功能，例如可以确定肿瘤类型特异基因。采用的分析方法可以是统计学中的假设检验等。2、考虑基因组合，将基因分组，研究基因的共同功能、相互作用以及协同调控等。多采用聚类分析等方法。3、尝试推断潜在的基因调控网络，从机理上解释观察到的基因表达谱。多采用反工程的方法。

本章首先介绍基因表达数据的来源和预处理方法；然后介绍基因表达数据分析的主要方法，即表达差异分析和聚类分析；最后简单介绍从基因表达数据出发研究基因调控网络的一些经典模型。

8.1 基因表达数据的获取

基因表达数据反映的是直接或间接测量得到的基因转录产物mRNA在细胞中的拷贝数或者水平（转录？？），这些数据可以用于分析哪些基因的表达发生了改变，它们有何相关性，在不同条件下基因是如何受影响的。它们在医学临床诊断、药物疗效判断、揭示疾病发生机制等方面有重要的应用。目前检测mRNA水平的方法有DNA微阵列、基因芯片、基因表达串行化分析（Serial analysis of gene expression，SAGE）、RT-PCR、EST测序等。目前，最主要的表达数据来自于基因芯片或cDNA微阵列，它们的原理是相同的，利用4种核苷酸之间两两配对互补的特性，使两条在序列上互补的单链形成双链，这个过程被称为杂交。基本技术是：在一个约1cm2大小的玻璃片上，将称为探针的核苷酸片段固定在上面，这个过程称为芯片制备；从细胞或组织中提取mRNA，通过RT-PCR合成荧光标记的cDNA，与芯片杂交；用激光显微镜或荧光显微镜检测杂交后的芯片，获取荧光强度，分析细胞中的mRNA的相对水平。

8.1.1 cDNA微阵列

cDNA微阵列最早是在1995年，由斯坦福大学研制并应用于基因表达分析的。首先将细胞内的mRNA逆转录成cDNA并分离，然后将分离得到的所有或部分cDNA（通常大于200bp）作为探针，用机器手点到玻璃片上，玻璃片上的每一个点包含一种cDNA分子，这样就制成了cDNA微阵列。固定在玻片上的cDNA探针可以通过测序得到序列或者其来源是已知的。在使用cDNA微阵列时，首先是提取组织或细胞系的mRNA样本，逆转录成cDNA 并用荧光素标记；然后把标记混合物加到cDNA微阵列上，与探针杂交，杂交过程完成后，清洗微阵列；然后用激光扫描仪扫描并获取荧光图像，对图像进行分析，得到cDNA芯片上每一个点的荧光强度值。荧光强度值定量反映了样本中存在的与探针互补的mRNA量，也就是反映了探针对应基因的表达水平。

在制造cDNA微阵列时，点样点的大小是不能保证完全一样的，点的排列也是不规则的，这样要比较不同微阵列图像的荧光绝对强度是不合理的，因此通常使用双色荧光系统来纠正点之间的差异。在制备样本时，使用两个样本，一个称为控制样本或对照样本，其cDNA 用红色（Cy5）或绿色（Cy3）荧光素标记，另一个为测量样本，其cDNA用与对照样本不同的绿色或红色荧光素标记。这两个样本按1：1的比例混合，同时与微阵列杂交，杂交后用不同波长的激光扫描，分别获取荧光强度，并成像。来自两个样本的基因如果以相同水平表达则显示黄色，而如果表达水平有差异，则图像显示红色或绿色。因此，cDNA微阵列的实验数据反映了两个样本中基因的相对表达水平。通常，在cDNA微阵列实验中对获取的原始图像数据必须进行归一化，例如基于全局强度值调整、强度相关归一化、玻片之间的对比归一化等，通常这些工作由与微阵列扫描系统配套的软件自动完成。为什么要进行归一化？如果用不同荧光素标记的是相同的样本，那么比率Cy5/Cy3（ratio值）的期望值为1，但由于Cy3和Cy5的标记效率不相等，或存在系统噪声等原因，得到的Cy5/Cy3往往不等于1，所以通过归一化可以使之回到1，并调整其它的测量值。归一化方法包括总密度（假设两个样本中的总RNA是相等的）、线性回归、Ratio统计、迭代log(ratio)平均值中心化等。

cDNA微阵列实验得到的值反映了基因的相对表达水平，即测量样本与对照样本之间荧光信号强度的比率或者比率取对数，这是一个无量纲的值，可用于比较一组实验中的基因相对表达水平。如果对照样本的信号非常低，那么这个比率就可能很大，因为可能主要是噪声信号，因此它很可能是无意义的，对于这些数据往往看作是不确定的，在后续分析时要注意这些数据，根据需要确定是否保留以及如何赋值。（是否是自己的语言？？？，或用我们的文章，陆老师）

8.1.2 寡核苷酸芯片

又称为基因芯片、DNA芯片。它是在玻璃片上按阵列固定寡核苷酸探针，这些探针是在片原位合成的。现有产品中应用最广泛的是Affymetrix公司制造的GENECHIP?芯片，它使用一种光掩模技术和传统的DNA合成化学的组合以非常高的密度制造寡核苷酸阵列。例如，Affymetrix公司的Human Genome U133芯片包含了100万个不同的寡核苷酸探针，代表了33000个人类基因。寡核苷酸芯片主要用于DNA多态性检测和基因表达分析，还可以用于微生物基因组的再测序。

寡核苷酸探针的长度通常为20-25bp，在检测mRNA表达水平时可能存在寡核苷酸之间的非特异性交叉杂交的冗余信息，可能会掩盖杂交信号；此外，对于特定的寡核苷酸，信号

强度对于寡核苷酸的碱基组成是敏感的。对于第一个问题，通常是采用匹配/失配（PM/MM）探针对的方法，即在设计一个特异的寡核苷酸(匹配)时，同时设计一个非特异的寡核苷酸探针，仅仅在中间位置有一个碱基替换（失配），这样可以用PM与MM之间的差值作为信号强度。为了解决第二个问题，在设计探针时，对于每一个待检测的mRNA包含多个寡核苷酸探针，例如为每一个转录本设计11-20个探针对来检测。

与cDNA微阵列不同的是，与寡核苷酸芯片杂交的是测量样本，而不是cDNA微阵列实验中的测量样本与对照样本的混合物。对于基因芯片的检测结果有两种，一种是P/A/M，表示有/无/不确定，另一种是信号强度。前者的结果主要是用来判断样本中有无特定基因的表达，这个结果对于部分实验，特别是一些定性实验是有意义的，例如判断肿瘤与正常情况下的细胞基因表达差异。当需要对几个不同条件下的基因表达情况进行分析时，对基因表达的相对变化更感兴趣，所以多采用第二种方式。有时基因表达数据的信号强度是负值，这是由于测量的信号小于背景信号或者背景/阴性控制样本的定义不正确造成的，对于前者，一般把负值做为0考虑，现在的Affymetrix的芯片分析系统已不产生负值。（？？）在考虑基因表达谱时，所采用的数据与cDNA微阵列数据一样，也是一系列测量样本与对照样本之间的信号强度比率或比率的对数值。实验得到的信号强度也是经过规格化的数值，规格化的方法很多，但归一化过程一般都包含在芯片扫描系统的图像处理软件中。

cDNA微阵列或基因芯片（以下统称微阵列）在用于基因表达分析时的一个最大优点是高通量性，在一次芯片实验中可以对成千上万个基因的表达进行并行测量。由于实验环节较多，虽然在设计芯片时可以通过添加阴性和阳性探针等手段来保证数据的可靠，但是需要提醒的是，数据的可靠性仍然是对数据进行后续分析时必须考虑的一个问题。

8.1.3 基因表达数据的网络资源

大量基于微阵列实验的基因表达数据是公开在Internet网上的，尤其是学术机构在发表论文时所用的实验数据都能免费提供给全世界的研究人员下载使用。作为学术论文的补充资料在网上发布的数据主要是文本文件或Excel格式的文件，这些数据往往都是经过归一化处理后的Ratio值或log2(Ratio)，对于寡核苷酸芯片数据有的是P/A/M（Present/Absent/Don’t Know）的表示或基因绝对表达值。因为这些数据文件没有包含原始的实验方案、实验材料、原始扫描图像、图像处理方法和数据归一化方法等信息，对于要比较、集成和整合分析来自不同研究小组的基因表达数据是非常困难的。主要原因是微阵列并不是在任何客观的个体上测量基因表达水平，大多数测量值仅仅是基因表达的相对变化，而且使用的并不是一个标准化的对照样本。同时，基因表达数据比基因组序列数据要复杂的多，这些数据仅仅在有具体的关于实验条件的描述时才是有意义的，对于不同的细胞类型，在不同的条件下都有一套转录本。因此，基于微阵列的基因表达数据存储量是非常大的，对于具有20000个探针的微阵列实验，以10um的分辨率扫描，产生3千万个离散的数据点，如果以tiff文件贮存，将占用~60Mb的硬盘空间。

一方面是基因表达数据量非常庞大，数据中蕴含着丰富的生物学知识，另一方面是这些数据没有注释，迫切需要一种标准来描述和存贮微阵列基因表达数据，同时建立公共的微阵列数据仓库。欧洲生物信息学研究所（EBI）与德国肿瘤研究中心(DKFZ)在1999年成立了MGED讨论组(The Microarray Gene Expression Data)。MGED（https://www.360docs.net/doc/a24138694.html,/）是一个国际性的成员联盟，参与人员包括生物学家、计算机科学家、数据分析学家。它的目标是促进由功能基因组学和蛋白组学研究产生的微阵列数据的共享。当前集中于建立微阵列数据注释和交换的标准，推动微阵列数据库建设和相关软件来实现这些标准，促进高质量的、经过注释的基因表达数据在生命科学领域的共享。该组织开发的微阵列数据标准称为

MIAME(the minimum information about a microarray experiment)，是对于解释和验证结果所必需的微阵列实验的最小信息描述。MIAME不是微阵列实验必须遵循的教条，而是一组指导方针，它将帮助微阵列数据库和数据分析工具的开发。同时，MGED组织开发了微阵列基因表达标记语言（MAGE-ML，Microarray Gene Expression - Markup Language），它是一种语言，用来描述和基于实验的微阵列信息的通讯，它基于XML，可以描述微阵列设计、微阵列制造信息，微阵列实验组织和实施信息，基因表达数据和数据表达结果。MIMAE标准和MAGE-ML语言受到了广泛关注。美国NCBI的Gene Expression Omnibus (GEO)、英国的EBI的ArrayExpress数据库都采用了该标准，斯坦福微阵列数据库（Stanford Microarray Database，SMD）也正在兼容该标准。

目前收集、存贮微阵列基因表达数据的最有影响的数据库和网站是GEO、ArrayExpress 和SMD。

GEO（https://www.360docs.net/doc/a24138694.html,/geo）是由NCBI在2000年开发的一个基因表达和杂交微阵列数据仓库，同时作为获取来自不同生物体的基因表达数据的在线资源。到2004年3月，数据仓库中包含内容605个Platforms，14391个Sample，816个Serial。Platform是关于物理反应物的信息，例如核酸、抗体和组织微阵列和SAGE数据等的基因表达数据被接受、增加和归档作为公共数据集。Series是关于样本集的信息，反映样本间的相关性和组织。

ArrayExpress(https://www.360docs.net/doc/a24138694.html,/arrayexpress/)是基于基因表达数据的微阵列公共知识库，目的是存储被很好注释的数据，当前包含多个基因表达数据集和与实验相关的原始图像集。ArrayExpress数据库接受MAGE-ML格式的数据递交或者通过MIAMExpress的基于Web的数据注释和递交工具。ArrayExpress提供一个简单的基于Web的数据查询界面，并直接与Expession Profiler数据分析工具相连，可以进行表达数据聚类，和其它类型的Web 数据发掘，并将进一步开发多个实验和数据库间的交叉查询。ArrayExpress数据库中的数据将与所有相关的由EBI维护的或在线的数据库相联接。

斯坦福微阵列数据库(SMD，https://www.360docs.net/doc/a24138694.html,/ )是一个使用Oracle作为数据库管理软件的关系数据库。SMD存储微阵列实验的原始、归一化数据和对应的图像文件。自从2002年1月1日起，到现在包括85篇学术论文，超过3500个双色点样DNA微阵列的实验数据，每年增加1000个微阵列实验的数据。另外，SMD提供数据获取、分析和可视化的界面，目前包括层次聚类和自组织映射等方法，还将加入k-平均聚类、单值分解和丢失值归纳等方法。

除了以上3个综合性的基因表达数据仓库外，还有一些专门的基因表达数据库，例如YMD (Yale Microarray Database，https://www.360docs.net/doc/a24138694.html,/microarray/)、ArrayDB （https://www.360docs.net/doc/a24138694.html,/arraydb/）、BodyMap（http://bodymap.ims.u-tokyo.ac.jp/）、ExpressDB(https://www.360docs.net/doc/a24138694.html,/ExpressDB/)、HuGE Index(Human Gene Expression Index，https://www.360docs.net/doc/a24138694.html,/welcome/index.html)等，这些数据库收集的数据往往具有物种特异性，使用比较方便。

8.2 基因表达数据预处理

一次微阵列实验能获得细胞在某一条件下的全基因组表达数据，包含成千上万个基因在细胞中的相对或绝对丰度，不同条件（细胞周期的不同阶段、药物作用时间、肿瘤类型、不

G?的数据矩阵M，通常情况下同病人等）下的全基因组表达数据就构成了一个N

x表示第i个基因在第j个条件下的表达水平值（在多数应用情G>>，其中每一个元素

N

ij

况下，表示的是Ratio 值或log(Ratio)值），行向量),,,(21.iN i i i x x x =x 代表基因i 在N 个条件下的表达水平，称为基因i 的表达谱，列向量T G i i i i x x x ),,,(21. =x 代表某一条件下的各基因的表达水平。

????????????=GN G G N N x x x x x x x x x 2122221

11211M （8-1）

注意排版，统一。公式背景为白底。

对基因表达数据进行分类、聚类等数据分析之前，往往需要进行预处理，包括对丢失数据进行填补、清除不完整的数据或合并重复数据等数据清洗，根据分析的目的进行数据过滤，以及针对分析方法选择适当的数据转换等预处理方法。

数据清洗是数据分析前必须进行的一项工作，对于基因表达数据，目的是去除表达水平是负值或很小的数据、或者明显的噪声数据(单个异常大或小的尖峰信号)，同时处理缺失数据。微阵列实验得到的数据一般是经过归一化处理的，每个点的信号强度是前景信号减去背景信号，因此有时会出现负值或很小的值，显然负值是没有生物学意义的。对于这些数据点，通过数据过滤步骤可以置为缺失或赋予统一的数值，例如对于寡核苷酸芯片数据，将低于100的数据全部设置为100。微阵列表达数据由于实验条件和芯片的因素，检测得到的信号强度往往与细胞中实际的mRNA 丰度之间没有对应关系，因此，通常是采用两个条件下的信号强度的比值，例如在cDNA 微阵列双色实验中，最后得到的往往是Ratio 值。而寡核苷酸单色实验的结果是信号强度，然而在处理一组数据时，也往往选择一个样本作为对照样本，将实验数据转换成Ratio 值。在计算Ratio 值时，如果参考样本的信号强度很小，就可能得到很大的Ratio 。如果一个基因谱中仅仅存在单个特别大的Ratio 值，称之为异常数据点(outlier)，这往往是由于噪声造成的。对于这个异常数据点，必须进行去除。数据的缺失对于某些后续数据分析方法（例如层次式聚类和PCA ）来说有着非常大的影响，甚至是致命性的，这时必须采取相应的方法。一种方法是直接过滤掉这些存在缺失数据项的行向量或列向量。另一种方法是设定阈值，计算一个基因表达谱中的缺失项数目，如果达到该阈值，则将该基因表达谱从数据矩阵M 中删除；如果没有达到阈值但存在缺失项，对这些缺失项可以进行插值。以0代替或用基因表达谱的平均值或中值进行代替，这些方法比较简单，但是否与真实值接近，很难进行评估。较为复杂和可靠的方法是，分析基因表达谱的模式，从中得到相邻数据点之间的关系，根据这种关系，利用相邻数据点估算得到缺失值。这种方法类似于k 近邻方法，需要有足够的完整的模式来发现有缺失值的相邻模式，需要有足够的值来确定它们的邻居。

在细胞中，基因表达有时空特异性，在某一条件下，发生表达的基因占基因总数的少部分，而大多数基因仅维持基础转录或不转录，转录本丰度很小，因此微阵列实验得到的数据矩阵中存在大量的基因表达谱曲线是平坦的，即基因表达水平变化很小。对于这些基因，往往不是生物学家所关心的，而它们的存在，却会大大增加数据分析的复杂性，而且会对一些分析方法的结果有干扰。对这些数据进行过滤是非常有必要的，可以给出一定的比例，使存在的基因占总数的多少，这是与分析目的相密切相关的，例如是分析细胞周期，可以多保留一些基因，而对于肿瘤特异基因表达谱分析，可以少保留一点基因。过滤这些基因所采用的标准有：①基因表达谱中最大值与最小值的差；②标准差；③均方根；④绝对值大于阈值的数据个数等。根据分析的对象和目的，可以选择以上一个或多个标准，确定阈值，来选择基

因表达谱。

基因表达谱数据经过过滤，在进行聚类分析等操作前，往往还需要进行数据转换。数据变换是将数据转换为适合数据挖掘的形式，可以根据需要构造出新的数据属性以帮助理解分析数据的特点，或者将数据规范化，使之落在一个特定的数据区间中。因此，数据转换包括对数转换和标准化两个过程。

许多DNA 微阵列实验的结果是测量样本与对照样本间信号强度的Ratio 值，对于Ratio 值，在大多数情况下是转换到对数(log)空间中进行处理，常用的对数底为2, e, 10。考虑时间序列上的基因表达数据，实验结果是相对于0时刻的表达水平。如图8.1所示，假设在时间点1，基因的表达水平没有改变，在时间点2，上调2倍，而时间点3，下调2倍，原始的比率值分别为1.0、2.0、0.5。在大多数应用中，需要把上调2倍和下调2倍看作是变化的相同幅度，只是方向不同。在Ratio 空间中，时间点1和2之间的差异是+1.0，而时间点1和3之间是-0.5，从数学角度看，上调2倍的数值是下调2倍的2倍。而在log 空间中，（为了简化，用2为底），这三个数据点分别为0、1.0、-1.0，上调2倍与下调2倍是关于0对称的。因此，对数转换可以使小于1的值变大，大于1的值变小，从而使它们关于0对称化，这种变换是否反映了一定的生物学意义，能更直观的了解基因的上调或下调的幅度？尚没有定论，但是对于大多数基因表达数据分析过程，都是在log 空间中进行的。

图8.1 表达数据的Ratio 和log2(Ratio)表示

数据的标准化是将所有的数据转换到同一个范围内，这样做的好处是方便比较和计算相关系数，缺点是在标准差接近0的时候，会产生大的噪声，这也是首先要进行数据过滤的一个重要理由。数据标准化按如下公式进行， ∑=---=N j i ij i

ij ij x x N x x x 1

2)(11 8.2 ∑==N j ij i x

x 11

8.3

通过标准化，使得每个基因表达谱的平均值为0，标准差为1。如果要求所有的数据在[0，1]之间，还需要进行如下转换

)/()(min max min x x x x x --= },,,min{21min N x x x x = },,,max{21max N x x x x =8.4 而要求数据满足[a,b]，则变换如下：

a x x x x a

b x +---=min max min ))(( 8.5

还有一种数据标准化方法是数据的中心化。对于来自细胞系的大量肿瘤样本与一个共同的对照样本比较，对于每一个基因，都有一系列的Ratio 值，相对于对照样本中那个基因的表达水平。因为对照样本通常对实验没有什么帮助，对照样本中的基因表达量是独立于分析的。这样可以通过调整每一个基因的数值来反映系列观察值的变化，例如平均值或者中值。这就是平均值/中值中心化，中心化可以减少参考样本的影响。中心化数据同样可以用于去除某些类型的偏差。许多双色荧光杂交实验的结果没有校正Ratio 值的系统偏差，它们是由于RNA 数量差异，标记效率和图像获取参数所造成的。这样的偏差对于所有的基因与一个固定数值的Ratio 有放大的效应。在log 空间的平均值和中值中心化有校正这种偏差的效果。数据中心化是基于这样的一种假设，在特定的实验中，基因的平均值期望比率是1.0(在log 空间中为0)。通常，更多的是使用中值中心化。

目前对数据预处理这种策略的作用还不是很清楚，还没有人进行系统的研究，提供有说服力的证据来帮助研究人员针对特定的任务选择特定的数据预处理的策略和方法。在具体应用时，往往是根据分析目的和个人经验选择不同的方法。

8.3 基因表达差异的显著性分析

在检测基因表达的微阵列实验中，有很大一部分是比较实验，目的是比较两个条件下的基因表达差异，从中识别出与条件相关的特异性基因，例如识别肿瘤特异性基因、药物特异响应基因等。为了提高实验的可靠性，对于两个条件，往往有两个以上的重复实验，但是由于微阵列实验仍然很昂贵，不可能重复足够的次数来满足实验数据分析的要求，因此需要采用一些比较复杂的方法来分析这些数据。对这些表达数据的分析目的就是要识别在两个条件下有显著表达差异的基因。何谓显著表达差异？通常是指一个基因在两个条件中表达水平的检测值在排除实验、检测等因素外，达到一定的差异，具有统计学意义，同时也具有生物学意义。分析方法有三类，一类称之为倍数分析，简单估计在两个条件中每个基因的表达水平的比值，设定阈值得到表达差异显著的基因；第二类方法是估计表达差异的置信度，采用的方法是t 检验和方差分析；第三类是建模的方法，通过确定两个条件下的模型参数是否相同来判断表达差异的显著性，例如贝叶斯方法。

8.3.1 倍数分析

早期基于cDNA 微阵列技术的比较实验，用倍数来分析基因表达水平差异，即计算两个条件下的表达水平的Ratio 值。用gi x 表示基因g 在条件i 中的表达水平测量值，因此，21/g g g x x r =表示基因g 在条件1和2中的表达水平比率。对于cDNA 微阵列实验，两个

条件的样本同时与同一微阵列进行杂交实验，得到的是成对数据，对于每次实验得到的数据计算g r ，最后计算多次重复实验的平均g r 。而对于寡核苷酸芯片，首先分别计算两个样本的重复微阵列实验的归一化表达水平的平均值，然后计算其比率。当g r =1时，基因g 的表达水平没有改变，而g r <1或g r >1意味着基因g 在两个条件下存在表达差异，特别是g r <1表示基因在条件1是下调的，而g r >1，表示在条件1是上调的。在具体应用中，如果一个基因的平均表达水平在两个条件下的变化超过一个常数，典型的常数是2，即g r >2或<1/2，则认为该基因的表达差异是显著的。然而，对表达数据仔细考察后提出，这样简单的2倍法并不能产生最优的结果，因为因子2在不同的表达水平上有相当不同的显著性。对于低表达水平的基因，其信噪比太低，用2倍法作为判断条件太宽松，而对于高表达基因，条件又太苛刻，往往小于2就有生物学意义上的表达差异显著性。在具体应用中，并没有明确的阈值，往往根据分析的具体要求由数据分析者自行确定，目前也没有关于这类简单的倍数分析方法的假阳性率和假阴性率的深入研究报道。

8.3.2 t 检验

对于两个条件下的多次重复实验，为了判断基因的表达差异是否具有显著性，在应用中较多的采用假设检验，包括两个条件下的t 检验和多个条件下的方差分析（ANOV A ），这里仅仅介绍t 检验。

零假设为210:g g u u H =，与之对应的备选假设是211:g g u u H ≠。t 统计量的计算公式如下：

2221212

1//n s n s x x t g g g g g +-= 8.6

其中 ∑==

j n j i gij gi n x x 1/，∑=--=i n j gi gij i gi x x n s 1

22)(11，i n 为某一条件下的重复实验次数。 如果t 超过了某个由给定的置信水平确定的阈值，就拒绝零假设，即认为基因g 在两个条件下的表达差异是显著的。因为在t 检验中，两个总体平均值之间的距离被经验得到的标准差归一化，可以克服固定倍数阈值方法的一些缺点。然而，对于微阵列数据的t 检验的基本问题是，即使用当前的高通量检测技术，实验仍然是花费很大或者实验过程很冗长，重复次数i n 经常较小，i n =2、3的小样本仍然非常普通。由于样本量小，导致总体方差被严重低估，得到的t 值就较大，因此会导致较高的假发现率(FDR, False Discovery Rate)，即通过t 检验得到的结果中表达差异不显著的基因数目较多。这样，需要更好的分析方法来克服这些缺点。

在t 假设检验中，经常使用的显著性水平是p=0.01，其意思是在零假设正确的情况，进行100次抽样，会有1次错误地拒绝了零假设。对于微阵列实验，检测的基因数目巨大，如

果微阵列上有10000个基因，采用p=0.01，将会有100个基因是由于偶然性而被错误认为是有表达差异显著的。这个数目已经对后续的生物学分析可能会产生很大的干扰，从而导致t 检验分析结果的不可靠或失去意义。

为了解决这个问题，可以对t 检验进行了改进，降低由于分母上方差小而带来的错误，因此对t 检验的计算公式修改如下：

021s s x x d g g g g +-= 8.7

}][][{22221211∑∑-+-=n g i g n g i g g x x x x a s

8.8 )2/()/1/1(2121-++=n n n n a

8.9 假设g d 的分布是独立于基因表达水平的。因为较低的表达水平会使g s 的值较小，导致g d 的变化可以较高。为了保证g d 独立于基因表达，在分母上增加s 0，s 0的选择是最小化g d 的变化系数。通过对设计的一组对照样本的分析，可以确定阈值，g d 大于阈值的基因被认为是表达差异显著的。

8.3.3 贝叶斯分析

由于微阵列数据噪声大、波动大，而且在大量数据的背后还有很多相关变量不能被观察到，因此贝叶斯方法可以用来分析微阵列表达数据。贝叶斯分析可以简单描述如下：

)(/)()|()|(D P M P M D P D M P = 8.10

其中，P(M|D)表示由观测数据集D 得到参数化模型)(w M M =为真的概率，称为后验概率；P(M)称为先验概率，表示在没有得到任何数据之前所估计的模型M 为真的概率。P(D|M)是指似然度，表示从模型M 得到一个观测数据集D 的概率。贝叶斯推断是通过参数估计和模型选择来实现任务的，最常用的方法是最大后验概率(MAP)估计和最大似然(ML)估计。

在用贝叶斯方法分析表达数据时，首先假设在给定条件下，一个基因的表达水平测量值是独立的，并满足正态分布。根据经验，这一假设是合理的，特别是表达水平的对数大致服从对数正态分布。对于重复实验，也可以引入伽玛分布、高斯/伽玛混合分布等。一个基因在一种条件下的表达测量值可以用一个正态分布),;(2

σμx N 来建模。对每个基因和每个条件，都对应有一个双参数模型),(2σμ=w ，似然函数可以由下式给出： 2222/])1()

([2/222)(),;(),|(σσσμσμs n u x n n i i e C x N D P -+---==∏ 8.11

i 取遍所有的重复测量，C 表示归一化常数。似然度取决于充分统计量n 、x 和s 2。

先验概率分布),(2σμP 的选择有几种，一般采用共扼先验分布。先验分布的四个超

参数构成向量),,,(20000σνλμα=

220020002/])([)12/(212)(),(σσνμλνσσσμ+--+--=u e C P 8.12

对于微阵列数据，采用一个μ和2σ相互不独立的先验分布很有意义。经过一些代数运算，可以推导出后验分布具有与先验分布相同的函数形式：

),;()/,;(),|,(2222n n n n N D P σνσλσμμασμΓ= 8.13

其中

2

00022002000000)()1(μλλσνσνννλλλμλλμ-++-+=+=+=+++=x n

n s n n n

m n n n

n n n n n 8.14

后验分布的参数以一种合理的方式将先验分布的信息和数据信息结合了起来。后验分布是贝叶斯分析的基本对象，它包含了μ和2σ所有可能取值的相关信息，可以通过多种方法

进行估计。

对于两个条件下的每个基因的表达测量值都与两个模型相关，根据模型参数的不同可以判断基因的表达差异是否显著。贝叶斯方法部分地克服了由于实验重复次数少造成的t 检验的缺陷，如果实验次数很少（2或3次），贝叶斯方法比t 检验法效果好，如果有5次重复实验，两者结果相似。在重复次数为2时，贝叶斯方法的假阳性率会提高。

8.4 基因表达谱聚类分析

对于基因表达谱数据的分析是目前生物信息学的研究热点和难点。转化为数学问题，分析任务是从数据矩阵M 中找出显著性结构，结构类型包括全局模型(model)和局部结构(pattern)，这些结构可以对基因或者条件形成亚类，结合生物学解释，这些亚类是细胞周期特异性的或者是肿瘤类型特异性的。因此，对基因表达谱数据的分析是数据挖掘问题，所采用的方法包括通过可视化进行探索性数据分析（Exploratory Data Analysis ）、描述建模(descriptive modeling)、分类、聚类和回归等。

目前，基因表达谱分析所采用的方法主要是聚类，其目的就是将基因或条件分组，对于条件之间存在时间依赖关系的基因表达谱，更多的是对基因进行分组；而对于不同来源的肿瘤样本，首先是确定特征基因，然后获取这些特征基因的表达谱数据子集，对条件进行分类，从而可以获取与肿瘤类型相关的特征，这些特征也称为基因组指纹（genomic fingerprint ）或签名（signature ），它们可以作为肿瘤的临床诊断标准或药效的评价标准。因为从数学的角度，对行或对列进行聚类所用的方法是一样的，不同的是对列进行分析时维度较大，样本数较小，结合生物学背景知识，通常需要降维，这个过程也称特征基因的选择过程，可以通过行向量的一些统计特征，例如方差、信息增益等，也可以用遗传算法等优化搜索算法来获取

特征基因。因此，以下部分主要是针对基因分组来介绍聚类算法。

从数学的角度，聚类得到的基因分组，一般是组内各成员在某数学特征上彼此相似，但与其它组中的成员不同。从生物学的角度，聚类分析方法所隐含的生物学意义（或基本前提）是，组内基因的表达谱相似，它们可能有相似的功能。当然，功能这个词太普通以至不能精确和定量，太广以至不能特指和特定意义。产物有相同功能的编码基因（例如对其它蛋白质有磷酸化作用），不一定共享相似的转录模式。相反，有不同功能的基因可能因为巧合或随机扰动而有相似的表达谱。尽管有许多意外的情况存在，大量功能相关的基因的确在相关的一组条件下有非常相似的表达谱，特别是被共同的转录因子共调控的基因，或者它们的产物构成同一个蛋白复合体，或者参与相同的调控路径。这就是在具体应用中，通过连坐（guilt-by-association ）证据，根据一个简单的有相似表达模式的基因聚类，可以指派未知基因的功能。

聚类分析是模式识别和数据挖掘中普遍使用的一种方法，是基于数据的知识发现的有效方法，特别适用于模式分类数不知道的情况。从机器学习的角度来看，聚类分析可以分为两类，即有监督聚类和无监督聚类。在有监督聚类中，对于基因表达谱数据，首先对基因进行功能分类，然后选择这些基因的表达谱作为训练集，通过有监督学习获取每一功能类的特征模式，利用这些特征模式，可以构建分类器，对未知功能的基因根据表达谱进行分类，从而能够对未知功能的基因预测功能，常用的方法有支持向量机、人工神经网络、朴素贝叶斯分类、多元logistic 回归等。无监督聚类不需要任何先验领域知识，根据数学特征提取分类标准，对数据进行分类，这种数学特征的例子有统计平均值、相关系数、协方差矩阵的本征值及本征向量等。无监督聚类在基因表达数据分析中应用的很多，主要有层次式聚类、k 平均、自组织特征映射网络、谱聚类等。本节将对在基因表达数据分析中常用的聚类方法进行具体介绍。

8.4.1 相似性度量函数

对基因表达谱进行聚类分析之前，必须首先确定反映不同基因表达谱相似程度的度量函数，根据该函数可以将具有相似程度高的基因聚为一类。在实际计算中，还可以用距离代替相似的概念，相似性度量被转化为两个基因表达谱之间的距离。距离越小，表达模式越相近，反之，则表达模式差异大。

常见的相似性度量有距离、点积、相关系数（correlation coefficient ）、互信息（mutual information ）等。假设两个基因表达谱分别为X=（x 1,x 2,…,x m ）和Y=（y 1,y 2,…,y m ）,距离函数d(X ，Y)必须满足如下条件：

d(X ，Y) ≧ 0

d(X ，Y) = d(Y ，X)

d(X ，Y) = 0 if X = Y

d(X ，Y) ≦ d(X ，Z) + d(Z ，Y)

欧氏距离（Euclidean distance ）是一个通常采用的距离定义，与测度的选择无关，它是在m 维空间中两个点之间的真正距离，两个基因表达谱之间的欧氏距离计算公式如下：

8.15

相关系数也是常用的相似性度量函数，计算公式如下：

8.16 ∑=-=m

i i i y x m Y X D 12)(1),(∑=Φ-Φ-=m i Y y X x m Y offset i X offset i Y X S 1

1))((),(

8.17

其中，G offset是G的各分量的均值，ΦG是标准方差。用上述两种相似性度量，可以找出表达模型相同或者变化趋势相同的基因，如图8.2(a)、(b)所示。欧氏距离、相关系数可以反映基因之间的共表达关系，两个基因表达谱间的距离小于或相关系数大于某个给定的阈值，就可以认为它们之间是共表达的。距离和相关系数之间存在相关性，在具体应用时，可以根据需要进行转换，例如d=1-r。

距离和相关系数反映的都是基因表达谱之间的相似性，这种相似性反映了基因的共表达行为，而基因的行为是复杂的，它们之间存在调控和被调控的关系，或者存在调控链，例如基因A调控B，B调控C，调控还有正性调控和负性调控之分。对于这些调控关系，它们的表达谱往往是不相似的，或者存在时延、或者存在反相，而基因表达的幅度也可能不相等。如何从数据中发现这些复杂的基因关系呢？互信息可能是一种有用的度量指标，其定义如下

8.18

MI(X,Y)是向量X和Y的互信息，H（X）是X的熵。

MI（X，Y）=H（X）+H（Y）-H（X，Y）

NMI（X，

Y）=MI（X，Y）/max[H（X），H（Y）] 8.19

NMI称为归一化互信息，它独立于单个信息熵，抓住了模式上的相似性。互信息聚类分析，没有规则上的约束。不象欧氏距离，MI不仅能确定负相关和非线性相关，同样可以反应正线性相关。因此，基于MI的聚类可能共享输入，但是对这些输入有不同的动力学过程。例如A,B得到C输入，A是上调，B是下调。

目前，还没有理论来指导如何选择最好的相似性度量，也许一个“正确”的距离在表达模式空间是不存在的，选择依赖于我们要问的问题。已知的共调控的基因在不同生物和基因调控网络中标准集可能有利于发现理论上的相似性度量。

8.4.2 无监督聚类方法（？？）

对于基因表达谱的聚类问题，由于目前对基因表达的系统行为了解得不全面，没有聚类的先验知识，所以通常采用无监督聚类方法。在基因表达数据分析方面，层次式聚类、k均（a）（b）（c）

图8.2 两个表达模式之间的关系（a）相似；（b）变化趋势一致；

（c）两个基因的调控制输入一样，但是调控结果不一样，甚至相反。

∑

=

-

=

Φ

m

i

offset

i

G m

G

G

,1

2

)

(

)

(

log

)

(

)

(

1

2i

m

i

i

x

p

x

p

X

H∑

=

-

=

值、自组织映射神经网络是应用中的常用方法。下面主要介绍这几种常用的聚类方法，并简单介绍一些其它方法。

8。4。2。1．简单聚类

假设有G 个基因，它们的表达谱分别用向量表示为..2.1,,,G x x x 。令任意一个基因的表达向量为第一个聚类的中心1z ,假设选择.11x z =。然后计算.2x 到1z 的距离D 21，如果D 21大于给定的阈值T ，则说明.2x 不属于第一类，应该分到另外的类。在这种情况下，建立一个新的聚类中心.22x z =。如果D 21小于阈值T ，则将.2x 分到第一类中。接着处理其它基因，在处理第i 个基因时，首先计算该基因的表达谱与现有各类中心的距离，假设与第j 类的距离D ij 最小，并且D ij

简单聚类算法的结果与第一个聚类中心的选择、基因的顺序、阈值T 以及基因表达谱在其空间的分布有关。该方法对于给定的一组表达数据模式进行初步分类提供了一种快速的算法。

8。4。2。2．层次聚类法

层次聚类法，在统计分析中也称为系统聚类法，原理与算法与第六章所介绍的系统发生树连锁构造方法类似，所不同的只是将所分析的数据由生物分子序列换成了这里的基因表达谱。该方法在基因表达谱聚类分析中是常用方法，它的优点是容易理解和实现，所得到的结果以树状图的形式表示，可以直观地观察基因之间的相互关系，尤其是类与类之间的关系。但是，基因表达谱的数量很多，往往要多于系统发生树分析时的物种数量，而且基因之间的相互关系信息也没有物种之间的多，所以对聚类结果的后续分析要比系统发生树分析复杂的多。对于表达谱聚类的结果还需要进一步分析基因的功能或者基因的序列特征，要通过剪枝才能得到分类，而剪枝的过程带有更多的主观性，这会导致丢失一些重要的信息或包括一些无关的信息。此外，在构建系统树时，已被合并的向量不再参与以后的分类，这会导致聚类结果与向量的次序有关，所以被认为是一种局部最优解的方法。

8。4。2。3．K 均值聚类

K 均值聚类在数据划分上不考虑类的分层结构问题，该算法使待聚类的所有样本到聚类中心的距离平方和最小，这是在误差平方和准则的基础上得到的。

K 均值聚类算法的基本过程如下：

（1） 任意选取K 个基因表达向量作为初始聚类中心z 1, z 2,…, z k ，在没有先验知识

的情况下，一般选择前K 个基因；

（2） 反复迭代计算。在第l 次迭代过程中，如果||)1(||||)1(||i j z x z x -<-

(i=1,2,…,K,i ≠j)，则将X 所代表的基因归于第j 类。按照上述办法处理所有

的基因；

（3） 经过上述处理，聚类可能发生变化，因此需要重新计算K 个新聚类中心： j=1,2,…,K 8.20 ∑

∈=+)(1)1(l f j j j N l Z x x

其中f j (l)为第l 次迭代中第j 个聚类的基因集合，N j 为该集合中基因的个数。

（4） 对于所有的聚类中心，如果Z j (l+1)=Z j (l)(j=1,2,…,K)，则迭代结束，得到最

后的聚类结果；否则转第2步，继续进行迭代计算。

聚类中心的个数K 、初始聚类中心的选择、基因排列的顺序以及基因表达数据的分布影响聚类的结果，当基因表达模式类别之间分离较远时，该算法可以取得令人满意的聚类分析结果。

8。4。2。4．自组织映射神经网络

人工神经网络技术在模式识别方面有着独特的优势，在生物信息学中的应用也非常广泛，如基因识别、蛋白质结构预测等。神经网络能够进行非线性数据处理，发现复杂的数据关系，其中，自组织映射神经网络（Self-Organizing Map ，SOM ）可以对模式数据进行自动聚类。

自组织特征映射是Kohonen 在1990年提出的类似大脑思维的一种人工神经网络方法，它是一种竞争学习算法，可以被认为是一种从N 维模式空间各点到输出空间少数点的映射。这一映射由系统本身完成，没有外部的监督，即聚类是以自组织的方式实现的。SOM 采用无教师学习训练，训练完成后，分类信息存储在网络个节点连接权值向量中，具有与权值向量相似的输入向量将分为一类。SOM 包括一维和二维模型，二维SOM 也称为KFM(Kohonen Feature Mapping)。它们的区别在于KFM 考虑邻近神经元的相互作用，即获胜神经元对周围神经元由于距离的不同会产生不同的影响。

KFM 的结构如图8.3所示，输出单元呈二维排列。Kohonen 认为，神经网络中邻近的各个神经元通过侧向交互作用彼此相互竞争，自适应地发展成检测不同信号的特殊检测器。该思想来自于大脑神经元的特性，即大脑皮层的神经元成二维空间排列，不同区域分工不同，各自对输入模式的不同特征敏感，神经元之间存在信息交互，主要是侧向交互，即某主响应神经元对周围神经元之间存在有近到远的不同的影响作用，而不是简单的侧抑制作用。对邻近神经元的交互作用函数有巴拿马草帽型或矩形型。 ?????>≤-=R R d R d j F cj c d 0 1)( 8.21

???>≤=R R d j F c d 0 1)(

图8.3 SOM 的结构图。

其中d cj 是输出单元c 与邻近单元j 之间在神经元平面上的距离，R 是交互作用半径。

KFM 的学习过程是：

（1） 用随机数设定权值初始值，并进行权向量归一化计算，在以后的每次修正权向量之后

也要进行归一化，使其满足 12

2==∑i j ij w W ；

（2） 反复进行以下运算，直到达到预定学习次数或每次学习后权值改变量小于某一阈值。

（a ）输入一个样本向量x ，计算各输出单元强度 ∑=

=i

i ij T j j x w W net x （b ）找出主兴奋单元c ，使)(max j j c net net = （c ）确定各输出单元在主兴奋单元影响下的兴奋度 y j =F c (j)

（d ）计算各权值修正量i j j ij x net y w )(-=?η

（e ）根据学习次数更新学习步长η和邻域交互作用半径R

学习过程可以采用从全局到局部的策略，在这种策略下，学习初期可设定较大的交互作用半径R ，随着学习过程的不断推进，逐步减小R ，直至不考虑对邻近单元的影响。邻域的形状可以为正方形、正六边形或圆形。

KFM 的聚类结果与k 平均法相似，它的优点是自动提取样本数据中的信息，同时也是一种全局的决策方法，能避免陷入局部最小，缺点在于必须实现人为设定类的数目与学习参数，而且学习时间较长。KFM 方法克服了K-均值聚类的一些缺点：它应用类间的全局关系，能提供大数据集内相似性关系的综合看法，便于研究数据变量值的分布及发现类结构。而且，它具有更稳健更准确的特点，对噪声稳定，一般不依赖于数据分布的形状。但它是一种拓扑保留的神经网络，产生不均衡分类。若不相关数据（如不变的“平坦”谱）或某种特定形式的谱过多，从而感兴趣的数据只能处于少数类，分辨率就可能很低。因此，在应用KFM 前，一般要对数据进行筛选。

8。4。2。5. 其它聚类方法

聚类方法是数据挖掘中的基本方法，方法很多，在基因表达谱的分析中，除了以上常用方法外，还有一些方法也被应用。然而，由于对聚类结果尚没有一种有效的方法进行评价，尤其是对聚类结果的进一步生物学知识发现尚没有新的分析思路和成功应用，因此这些方法更多的是方法研究，而没有进入到实用阶段。但是这些方法有不同的原理，能够提取不同数据模式的特征，有可能对具体的数据得到有意义的结果，发现更多的生物学知识。这里，简单介绍这些方法的原理，更详细的介绍请参看相关文献。

不用公式，简述核心思想！

（1）模糊聚类分析方法

这是一种模拟人类的思维方法，通过隶属度函数来反映某一对象属于某一类的程度。基本过程包括：计算两两基因表达谱之间的相似系数构建模糊相似矩阵，利用模糊数学中的传递闭包计算方法得到模糊等价矩阵，选择不同的置信水平就可以从模糊等价矩阵中构建动态聚类图。对于特定的置信水平，可以对基因进行分类。该方法的优点是利用了模糊数学中的隶属度概念，能够更好的反映基因之间的相互关系，而且它是一种全局的优化方法，与向量的顺序无关。

（2）模糊C 均值算法

该方法同样将模糊数学中的隶属度概念引入到k 均值聚类方法中，C 均值方法也称为硬分类。判断一个基因表达谱是否属于某一类时，以隶属度来确定属于每一类的可能性，因此

在计算聚类中心时也必须考虑所有表达谱向量的隶属度。聚类中心向量m 和隶属度μ的计算公式如下：

[][],,,2,1,)()(11

c j n

i b i j n i i b i j j ==∑∑==x μ

x x m 8.22 c j n i c k b k i b j i i j ,,2,1,,2,1,)1()

1()(1)1(12

)1(12 ==--=∑=--m x m x x μ 8.23 {}n i x i ,,2,1, =是n 个样本组成的样本集合，c 为预定的类别数目，c i m i ,,2,1, =为每个聚类的中心，)(i j x μ是第i 个样本对于第j 类的隶属度函数。聚类的评价是最小化聚类损失函数，计算公式如下：

[]∑∑==-=c j n i j i b i j f m x x J 112)

(μ， 8.24

其中b>1是一个可以控制聚类结果的模糊程度的函数。当算法收敛时，就得到了各类的聚类中心和各个样本对于各类的隶属度值，从而完成了模糊聚类划分。如果需要，还可以将模糊聚类结果进行去模糊化，即用一定的规则把模糊聚类划分转化为确定性分类。

（3）超顺磁性聚类（superparamagnetic clustering ）

这是一种基于模拟非均匀铁磁物质的物理特性的聚类方法，将数据聚类问题视为检验不均匀Potts 模型的平衡特性。根据表达谱向量之间的距离矩阵构建图，顶点是数据点，如果两个点之间的距离满足K-mutual-neighbor 准则，则称为邻居。数据点的聚类等价于有权图的划分。给每一数据点分配一个Potts 自旋子，相邻数据点间引入强度随距离下降的相互作用函数。非均匀Potts 模型系统随温度变化表现出三相：在低温下，所有自旋子呈现完全有序的排列，系统为铁磁相；随着温度的升高，小区域自旋子形成磁化“颗粒”，附属于同一“颗粒”者相互间产生强耦合，而无关者间相互作用很弱，不同“颗粒”的排列呈无序状态，为超顺磁相；在高温下，系统不表现任何有序性，为顺磁相。在超顺磁相的转换温度下，磁化率表现出显著的峰值。原则上，超顺磁相可以有一系列的转换点。随着温度的升高，系统可以首先分裂为两类，其中每一类又可以分裂为更多的子类，这样，数据就分层组织为类。SPC 算法的优点是对噪声及初始化不敏感，因为类由系统的综合性质产生。由磁化率的峰值很容易鉴别主要的分界，从而能清楚显示类的构成和分界，并且在每一个分辨率上能自动确定类数。

（4）谱聚类（Spectral clustering ）

k 均值聚类对于聚类中心相距较远的数据样本具有很好的聚类效果，而对于具有同心圆特征的数据样本很难得到好的效果，而谱聚类能够很好的对这类样本进行聚类。谱聚类是一种基于矩阵特征向量的方法，也是一种能根据顶点之间的权值对图进行划分的方法。用图表示基因表达谱矩阵，基因表达谱可以看作是一组顶点，连接顶点的边的权值反映了两个表达

谱之间的相似性，这样就得到有权无向图G=(V,A)，聚类过程等价于将G划分为不连接的子集，这可以通过简单地删除连接边来实现。聚类过程包括两个步骤：第1步是将表达谱空间转化为相似度矩阵的特征向量空间；第2步是应用简单的k均值法对特征向量空间的特征向量进行聚类，得到的结果就对应了基因表达谱的聚类结果。该方法可以根据特征值自动确定分类数目。谱聚类在聚类过程中进行了特征空间的转换，可以将大的空间转化为较小的空间，从而可以更快速的处理大规模的数据。

（5）双向聚类法

以上所述方法主要是对基因表达谱进行聚类，将具有相似表达模式的基因分为一组，这是在应用中要解决的主要问题。数据矩阵中列表示了实验条件，可以是时间序列也可以是不同的肿瘤样本或病人样本。从数学的角度，数据矩阵中的每一列对基因表达谱的行聚类结果有着不同的影响，挑选部分列进行行聚类的结果与选择所有列进行聚类的结果是不同的。从生物学应用的角度，对列进行聚类可以发现各实验条件之间的相互关系，例如同一肿瘤类型的样本可以聚成一类，它们具有相似的基因表达模式，这样就可以基于表达谱对肿瘤进行分类。将两者结合起来，就称为双向聚类法，在行和列两个方向上进行聚类分析，通常采用贪婪迭代检索的方法来发现子矩阵或稳定的类，这些子矩阵中感兴趣的模式具有特定的生物学意义。在应用中，可以根据具体的目的，确定以降低基因维还是以降低样本维为主，通过迭代得到稳定的若干样本分类或基因分组。

8.4.3 有监督聚类方法（分类？？）

如上所述，对于基因表达谱数据矩阵M，在对行或列向量进行聚类时很少或根本不需要生物学知识，所得到的结果是否具有生物学意义，还需要进一步的分析。然而，对于特定的数据矩阵，都有特定的生物学背景，例如基因之间存在相关性，可以分成若干功能类，样本之间也存在某种关系，可能来自相同肿瘤的病人。将这些知识应用于聚类方法中，挑选部分样本作为训练样本，可以提取特征构建分类器，然而对其它的样本进行分类，实现聚类的目的。这种方法称为有监督聚类方法，可以分类两类，一类是基于模型的方法，另一类是不基于模型的方法，包括人工神经网络和属于机器学习的支持向量机。这里主要介绍支持向量机的原理和基本方法。

8。4。3。1 支持向量机（support vector machine，SVM）

统计学习理论是目前针对小样本统计估计和预测学习的最佳理论，它从理论上较系统地研究了风险最小化原则成立的条件，有限样本下经验风险与期望风险的关系，以及如何利用这些理论找到新的学习原则和方法等问题。支持向量机是实现统计学习理论思想的具体方法，可以用于解决小样本、非线性及高维学习问题，不需要足够的样本来构建特征空间，而仅仅需要少量的样本来构建分类界面，这些用来构建分类界面的样本就称为支持向量。

考虑图示的二维两类线性可分情况，图8.4中实心点和空心点分别表示两类的训练样本，H是把两类没有错误地分开的分类线，H1，H2分别为过各类样本中离分类线最近的点，且平行于H的直线，H1、H2之间的距离叫做两类的分类空隙或分类间隔，最优分类线就是要求分类线不但能将两类无错误地分开，而且要使两类的分类空隙最大，前者是保证经验风险最小。使分类空隙最大，实际上就是使推广性的界中的置信范围最小，从而使真实风险最小。对于高维空间，最优分类线就成为最优分类面。

图8.4 两维两类的分类线 假设线性可分样本集为),(i i y x ,i=1,2,…,n,}1,1{-+∈i y ,d

R ∈x ，则分类判别函数为b g +?=x w x )(,分类面方程为0=+?b x w 。将判别函数归一化，使两类所有样本都满足1|)(|≥x g ，即使离分类面最近的样本的判别值为1，这样分类间隔||||/2w =。要使间隔最大，必须使||||w 最小，要求分类线对所有样本正确分类，就要求它满足：

01][≥-+?b y i i x w 8.26

因此满足上述条件且使||||w 最小的分类面就是最优分类面。过两类样本中离分类面最近的点且平行于最优分类面的超平面H 1,H 2上的训练样本就是使等号成立的那些样本，叫支持向量（Support Vectors ），它们支撑最优分类面。

那么，如何求最优分类面？根据上述讨论，最优分类面问题可以表示成在条件8.26的约束下。求函数)(2/1||||2

1)(2w w w w ?==φ的最小值。为此，可以构建Lagrange 函数 ∑=-+?-=n i i i i b y b L 1

2}1][{||||21),,(x w w x w α， 8.27 其中0>i α，为Lagrange 系数，问题是对w 和b 求Lagrange 函数的极小值，分别对w 和b 求偏微分并令它们等于0，在约束条件

∑==n i i i

y 10α，n i i ,,1,0 =≥α 8.28

之下，对i α求解函数的最大值，

∑∑==?-=n

j i j i j i j i n i i y y Q 1,1

)(21)(x x αααα 8.29 若*i α为最优解，则∑=i i i y x w *

*α，即最优分类面的权系数向量是训练样本向量的线性

组合。

求解上述问题后得到的最优分类函数是

})(sgn{})sgn{()(****∑+?=+?=b y b x f i i i x x x w α 8.30

sgn()为符号函数，由于非支持向量对应的i α均为0，因此式中的求和实际上只对支持向量进行，而b *是分类的阈值，可以由任意一个支持向量求得，或通过两类中任意一对支持向量取中值求得。从最终的分类判别函数中可以看到，只含有待分类样本与训练样本中的支持向量的内积运算。

对于二类情况的分类，就是要构建最优分类面或超平面。对于大多数真实世界的问题，包含不可分的数据，也就是在数据空间中不存在一个超平面可以成功的区分二类样本，一个解决办法是将数据映射到更高维的空间，在那里定义一个可分隔的超平面。这个更高维的空间称为特征空间，而对应的训练样本所在的空间叫输入空间。问题就成为如何选择一个足够维度的特征空间，使得训练数据可分？这里有计算和学习理论的代价。还有一个风险是人为分隔数据可能使学习系统有过学习的风险。SVM 可以避开这两个问题，避免过学习的问题，通过在特征空间中从众多可分隔超平面中选择最大边缘分隔超平面。

SVM 的基本思想是首先通过非线性变换将输入空间变换到一个高维的特征空间，然后在这个新空间中求取最优线性分类面，而这种非线性变换是通过定义适当的内积函数实现的。用于超平面分类点的决策函数仅仅包含特征空间中点之间的点积，因此寻找特征空间中可分隔超平面的算法完全可以用输入空间中的向量和特征空间中的点积来陈述。SVM 能够定位超平面，不需要明确的表示特征空间，只需要简单的定义称为核心函数的函数，它在特征空间中起着点积的作用，这个技术避免了明确表示特征向量的计算负担。

SVM 的判别函数为：})(sgn{

)(1**∑=+?=n i i i i b K y x f x x α 8.31

一些简单的核心函数有： 多项式形式的内积函数q K )1(),(+?=y x , q =1,2,3,… 8.32

核函数型内积:)2/||||exp(),(22αy x --=K 8.33 α是高斯宽度，等于正样本与最近负样本间欧氏距离的中值，得到的SVM 是一种径向基函数分类器。

8。4。3。2 基于模型的聚类方法

在基因表达数据显著性分析时，我们提到贝叶斯方法的应用，它是根据重复实验数据分别构建模型，通过比较两个条件下的模型参数来确定表达差异的显著性。同样，基因表达谱分析也可以引入建模的方法，假定每一个基因表达谱是由一种内在的概率模型产生的，它满

足一定概率分布或者有限组合，例如多元正态分布，而所有的基因表达谱是由若干个概率模型产生的，这样就可以通过确定这些概率模型，来实现对基因表达谱聚类的目的。高斯混合模型在许多聚类应用中都是很有用的工具。与启发式聚类算法比较的最大的优势是，没有严格的方法来确定类的数目或最佳的聚类方法。

基于混合模型聚类假定数据中蕴含的每一组(成分)由一种内在的概率分布混合产生，例如，高斯混合模型中，独立多元观测y 1,y 2,…,y n 组成的数据Y 中每一成分k 由参数为u k （均值矢量）和Σk （协方差矩阵）的多元正态分布模拟，

)

2det()}

()(2/1exp{),|(1k k i k T k i k k i k f Σu y Σu y Σu y π---=- 8.34 算法的目标是由数据预测参数u k 和Σk ，并确定相应于这些参数预测的类。为了便于参数预测，对Σk 进行特征值分解：T k k k k k D A D Σλ=， D 为特征向量组成的正交矩阵，确定成

分的取向，A 为对角矩阵，元素与Σk 的特征值成正比，确定成分的形状，λ为标量，确定成分的体积。约束和改变部分参数，可以得到这一框架下的一组模型，如等体积球模型（EI ：I Σλ=k ）I 为单位矩阵、不等体积球模型(VI:I k k λ=∑)等，以适应数据特征的变化。每一种协方差矩阵与类数组合情况相应于不同的概率模型，由贝叶斯信息准则(BIC)估计每一种模型下数据被观测到的概率，计算结果的BIC 的得分，最后选择BIC 的得分做大的模型和参数。

8.4.4 聚类结果的可视化

对微阵列基因表达谱的聚类分析得到的一般是基因的分组信息和基因表达谱在某种关系上的相互关系，如何进一步了解同组基因的表达谱差异？不同组基因之间是否存在相互关系？哪些基因是上调或下调的，幅度变化有多大，在何种条件下发生较大的变化？等等生物学家普遍关注的问题，需要通过简单、直观的方法来得到这些答案。数据对于大多数人来说太抽象了，而文字的描述又不能反映数据的本质，所以采用可视化的方法可以大大方便对表达谱分析结果的理解，有利于对这些问题的回答。

在数据挖掘中，可视化方法可以用来发现数据中的结构，在数据探索中占有特别重要的地位，这是因为人类的眼睛和大脑具有强大的结构探测能力——这是长期进化的成果。可视化方法就是以各种可以发挥出人类在模式处理方面的特殊能力的方式来显示数据。可视化方法在数据挖掘中具有非常重要的地位，它是筛选数据寻找未知数据关系的理想方法。可视化方法是生物信息学中常用的一种表示和发现生物学知识的方法，典型的例子有蛋白质三维构象的表示、转录因子DNA 结合位点一致序列的logo 表示、基因调控网络的图形表示等。 基因表达谱的数据特点是基因数目多，样本少，而每一个数据点所表示的是基因在特定条件下的表达数据，经过数据预处理后，往往是与参考样本表达水平的比值的对数值，负值表示基因表达水平的下调，正值表示上调，值的绝对值反映了变化的幅度。因此，对于表达谱数据及其聚类结果通常用两种方式表示：

①彩色盒图和树状图(dendrogram)。所谓的彩色盒图是用不同的颜色小方盒来表示基因表达谱数据，一个小方盒表示一个数据，绿色表示数据小于0，红色表示数据大于0，相对应其生物学意义是对应的基因在相应的条件下是上调还是下调，而颜色的深浅反映了数据的绝对值大小。基因表达谱的很多聚类结果是以树状图的形式表示的，例如层次式聚类、模糊聚类和SPC 等，对于这些结果的可视化是通过树状图和彩色盒图相结合表达的（如图8.5a 所示）。