基因表达数据分析

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

2、对数据矩阵的分析，回答一些生物学问题，例如，基因的功能是什么？在不同条件 或不同细胞类型中， 哪些基因的表达存在差异？在特定的条件下， 哪些基因的表达发生了显 著改变， 这些基因受到哪些基因的调节， 或者调控哪些其它的基因？哪些基因的表达是条件 特异性的， 根据它们的行为可以判断细胞的状态（正常或癌变）？等等。对这些问题 的回答， 结合其他生物学知识和数据有助于阐明基因的调控路径和基因之间的调控网络。 揭 示基因调控路径和网络是生物学和生物信息学共同关注的目标，是系统生物学 （Systems Biology ，在附录中增加解释条目！ ）研究的核心内容。目前，对基因表达数据的分析主要是 在三个逐

3、渐复杂的层次上进行： 1、分析单个基因的表达水平，根据在不同实验条件下，该 基因表达水平的变化， 来判断它的功能， 例如可以确定肿瘤类型特异基因。 采用的分析方法 可以是统计学中的假设检验等。 2、考虑基因组合，将基因分组，研究基因的共同功能、相 互作用以及协同调控等。多采用聚类分析等方法。3、尝试推断潜在的基因调控网络，从机理上解释观察到的基因表达谱。多采用反工程的方法。本章首先介绍基因表达数据的来源和预处理方法； 然后介绍基因表达数据分析的主要方 法，即表达差异分析和聚类分析； 最后简单介绍从基因表达数据出发研究基因调控网络的一 些经典模型。8.1 基因表达数据的获取基因表达数据反映的是直

4、接或间接测量得到的基因转录产物 mRNA 在细胞中的拷贝数 或者水平（转录？） ，这些数据可以用于分析哪些基因的表达发生了改变，它们有何相关 性，在不同条件下基因是如何受影响的。 它们在医学临床诊断、药物疗效判断、揭示疾病发 生机制等方面有重要的应用。目前检测 mRNA 水平的方法有 DNA 微阵列、基因芯片、基 因表达串行化分析 （ Serial analysis of gene expression ，SAGE ）、RT-PCR、EST 测序等。目前， 最主要的表达数据来自于基因芯片或 cDNA 微阵列，它们的原理是相同的，利用 4 种核苷 酸之间两两配对互补的特性，使两条在序列上互补的单

5、链形成双链，这个过程被称为杂交。 基本技术是：在一个约 1cm2 大小的玻璃片上，将称为探针的核苷酸片段固定在上面，这个 过程称为芯片制备；从细胞或组织中提取mRNA，通过RT-PCR合成荧光标记的 cDNA，与芯片杂交；用激光显微镜或荧光显微镜检测杂交后的芯片，获取荧光强度，分析细胞中的 mRNA 的相对水平。8.1.1 cDNA 微阵列cDNA 微阵列最早是在 1995 年，由斯坦福大学研制并应用于基因表达分析的。首先将 细胞内的 mRNA 逆转录成 cDNA 并分离，然后将分离得到的所有或部分 cDNA （通常大于 200bp）作为探针，用机器手点到玻璃片上，玻璃片上的每一个点包含一种c

6、DNA分子，这样就制成了 cDNA 微阵列。固定在玻片上的 cDNA 探针可以通过测序得到序列或者其来源 是已知的。在使用 cDNA 微阵列时， 首先是提取组织或细胞系的 mRNA 样本，逆转录成 cDNA 并用荧光素标记；然后把标记混合物加到 cDNA 微阵列上，与探针杂交，杂交过程完成后， 清洗微阵列；然后用激光扫描仪扫描并获取荧光图像，对图像进行分析，得到cDNA 芯片上每一个点的荧光强度值。荧光强度值定量反映了样本中存在的与探针互补的mRNA 量，也就是反映了探针对应基因的表达水平。在制造 cDNA 微阵列时， 点样点的大小是不能保证完全一样的， 点的排列也是不规则的， 这样要比较不同

7、微阵列图像的荧光绝对强度是不合理的， 因此通常使用双色荧光系统来纠正 点之间的差异。在制备样本时，使用两个样本，一个称为控制样本或对照样本，其cDNA用红色（Cy5）或绿色（Cy3）荧光素标记，另一个为测量样本，其cDNA用与对照样本不同的绿色或红色荧光素标记。这两个样本按1：1 的比例混合，同时与微阵列杂交，杂交后用不同波长的激光扫描， 分别获取荧光强度， 并成像。 来自两个样本的基因如果以相同水平 表达则显示黄色，而如果表达水平有差异，则图像显示红色或绿色。因此，cDNA 微阵列的实验数据反映了两个样本中基因的相对表达水平。通常，在cDNA 微阵列实验中对获取的原始图像数据必须进行归一化，

8、 例如基于全局强度值调整、 强度相关归一化、 玻片之间的对 比归一化等，通常这些工作由与微阵列扫描系统配套的软件自动完成。为什么要进行归一 化？如果用不同荧光素标记的是相同的样本，那么比率 Cy5/Cy3 （ ratio 值）的期望值为 1， 但由于 Cy3 和 Cy5 的标记效率不相等，或存在系统噪声等原因，得到的 Cy5/Cy3 往往不等 于 1，所以通过归一化可以使之回到 1，并调整其它的测量值。归一化方法包括总密度（假 设两个样本中的总 RNA是相等的）、线性回归、Ratio统计、迭代log（ratio）平均值中心化等。cDNA 微阵列实验得到的值反映了基因的相对表达水平，即测量样本与

9、对照样本之间荧 光信号强度的比率或者比率取对数， 这是一个无量纲的值， 可用于比较一组实验中的基因相 对表达水平。 如果对照样本的信号非常低， 那么这个比率就可能很大， 因为可能主要是噪声 信号， 因此它很可能是无意义的， 对于这些数据往往看作是不确定的， 在后续分析时要注意 这些数据，根据需要确定是否保留以及如何赋值。（是否是自己的语言？，或用我们的文章，陆老师）8.1.2 寡核苷酸芯片又称为基因芯片、 DNA 芯片。它是在玻璃片上按阵列固定寡核苷酸探针，这些探针是 在片原位合成的。现有产品中应用最广泛的是 Affymetrix 公司制造的 GENECHIP ?芯片，它 使用一种光掩模技术和

10、传统的 DNA 合成化学的组合以非常高的密度制造寡核苷酸阵列。 例 如， Affymetrix 公司的 Human Genome U133 芯片包含了 100 万个不同的寡核苷酸探针，代 表了 33000 个人类基因。寡核苷酸芯片主要用于 DNA 多态性检测和基因表达分析，还可以 用于微生物基因组的再测序。寡核苷酸探针的长度通常为20-25bp，在检测mRNA表达水平时可能存在寡核苷酸之间的非特异性交叉杂交的冗余信息，可能会掩盖杂交信号； 此外，对于特定的寡核苷酸， 信号强度对于寡核苷酸的碱基组成是敏感的。对于第一个问题，通常是采用匹配/失配（PM/MM ）探针对的方法，即在设计一个特异的寡核

11、苷酸（匹配）时，同时设计一个非特异的寡核苷酸探针，仅仅在中间位置有一个碱基替换（失配），这样可以用PM与MM之间的差值作为信号强度。为了解决第二个问题，在设计探针时，对于每一个待检测的mRNA包含多个寡核苷酸探针，例如为每一个转录本设计11-20个探针对来检测。与cDNA微阵列不同的是，与寡核苷酸芯片杂交的是测量样本，而不是cDNA微阵列实验中的测量样本与对照样本的混合物。对于基因芯片的检测结果有两种，一种是P/A/M，表示有/无/不确定，另一种是信号强度。前者的结果主要是用来判断样本中有无特定基因的 表达，这个结果对于部分实验，特别是一些定性实验是有意义的，例如判断肿瘤与正常情况下的细胞基因

12、表达差异。 当需要对几个不同条件下的基因表达情况进行分析时，对基因表达的相对变化更感兴趣，所以多采用第二种方式。有时基因表达数据的信号强度是负值，这是由于测量的信号小于背景信号或者背景/阴性控制样本的定义不正确造成的，对于前者，一般把负值做为0考虑，现在的Asymetrix的芯片分析系统已不产生负值。（？）在考虑基因表达谱时，所采用的数据与cDNA微阵列数据一样，也是一系列测量样本与对照样本之间的信号强度比率或比率的对数值。实验得到的信号强度也是经过规格化的数值，规格化的方法很多，但归一化过程一般都包含在芯片扫描系统的图像处理软件中。cDNA微阵列或基因芯片（以下统称微阵列）在用于基因表达分析

13、时的一个最大优点是 高通量性，在一次芯片实验中可以对成千上万个基因的表达进行并行测量。由于实验环节较多，虽然在设计芯片时可以通过添加阴性和阳性探针等手段来保证数据的可靠，但是需要提醒的是，数据的可靠性仍然是对数据进行后续分析时必须考虑的一个问题。8.1.3基因表达数据的网络资源大量基于微阵列实验的基因表达数据是公开在In ternet网上的，尤其是学术机构在发表论文时所用的实验数据都能免费提供给全世界的研究人员下载使用。作为学术论文的补充资料在网上发布的数据主要是文本文件或Excel格式的文件，这些数据往往都是经过归一化处理后的Ratio值或log2（Ratio），对于寡核苷酸芯片数据有的是P

14、/A/M （ Present/Absent/Don'tKnow）的表示或基因绝对表达值。因为这些数据文件没有包含原始的实验方案、实验材料、原始扫描图像、图像处理方法和数据归一化方法等信息，对于要比较、集成和整合分析来自不同研究小组的基因表达数据是非常困难的。主要原因是微阵列并不是在任何客观的个体上测量基因表达水平，大多数测量值仅仅是基因表达的相对变化，而且使用的并不是一个标准化的对照样本。同时，基因表达数据比基因组序列数据要复杂的多，这些数据仅仅在有具体的关于实验条件的描述时才是有意义的，对于不同的细胞类型， 在不同的条件下都有一套转录本。因此，基于微阵列的基因表达数据存储量是非常大的

15、，对于具有20000个探针的微阵列实验，以10um的分辨率扫描，产生 3千万个离散的数据点，如果以tiff文件贮存，将占用60Mb的硬盘空间。一方面是基因表达数据量非常庞大，数据中蕴含着丰富的生物学知识，另一方面是这些数据没有注释，迫切需要一种标准来描述和存贮微阵列基因表达数据，同时建立公共的微阵列数据仓库。欧洲生物信息学研究所（EBI）与德国肿瘤研究中心（DKFZ）在1999年成立了MGED 讨论组（The Microarray Gene Expression Data） 。MGED （ / ）是一 个国际性的成员联盟，参与人员包括生物学家、计算机科学家、

16、数据分析学家。它的目标是 促进由功能基因组学和蛋白组学研究产生的微阵列数据的共享。当前集中于建立微阵列数据注释和交换的标准，推动微阵列数据库建设和相关软件来实现这些标准，促进高质量的、经过注释的基因表达数据在生命科学领域的共享。该组织开发的微阵列数据标准称为MIAME(the minimum information about a microarray experiment),是对于解释和验证结果所必需的微阵列实验的最小信息描述。MIAME不是微阵列实验必须遵循的教条，而是一组指导方针，它将帮助微阵列数据库和数据分析工具的开发。同时，MGED组织开发了微阵列基因表达标记语言 (MAGE-ML

17、，Microarray Gene Expression - Markup Language )，它是一种语 言，用来描述和基于实验的微阵列信息的通讯，它基于 XML，可以描述微阵列设计、微阵 列制造信息，微阵列实验组织和实施信息，基因表达数据和数据表达结果。MIMAE标准和MAGE-ML 语言受到了广泛关注。美国NCBI的Gene Expression Omnibus (GEO)、英国的EBI的ArrayExpress数据库都采用了该标准，斯坦福微阵列数据库(Stanford MicroarrayDatabase, SMD)也正在兼容该标准。目前收集、存贮微阵列基因表达数据的最有影响的数据库和

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

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

20、联接。斯坦福微阵列数据库 (SMD , / )是一个使用 Oracle作为数 据库管理软件的关系数据库。SMD存储微阵列实验的原始、归一化数据和对应的图像文件。 自从2002年1月1日起，到现在包括 85篇学术论文，超过3500个双色点样DNA微阵列的实验数据，每年增加1000个微阵列实验的数据。另外，SMD提供数据获取、分析和可视化的界面，目前包括层次聚类和自组织映射等方法，还将加入k-平均聚类、单值分解和丢失值归纳等方法。除了以上3个综合性的基因表达数据仓库外，还有一些专门的基因表达数据库，例如YMD (Yale Microarr

21、ay Database ，/microarray/)、 ArrayDB(/arraydb/ )、 BodyMap ( http:/bodymap.ims.u-tokyo.ac.jp/ )、 ExpressDB( /ExpressDB/ )、 HuGE Index(Human Gene Expression Index， /welcome/index.html) 等，这些数据库收集的数据往往具有物 种特

22、异性，使用比较方便。8.2基因表达数据预处理一次微阵列实验能获得细胞在某一条件下的全基因组表达数据，包含成千上万个基因在细胞中的相对或绝对丰度，不同条件(细胞周期的不同阶段、药物作用时间、肿瘤类型、不 同病人等)下的全基因组表达数据就构成了一个G N的数据矩阵 M，通常情况下G N，其中每一个元素Xj表示第i个基因在第j个条件下的表达水平值 (在多数应用情况下，表示的是 Ratio值或log（Ratio）值），行向量Xj. =（Xixi2,，xiN ）代表基因i在N个条件下的表达水平，称为基因i的表达谱，列向量xj=（x1i,x2i,xGi）T代表某一条件下的各基因的表达水平。XiiX21xG

23、1X12X22aXG2XlNX2N(8-1)注意排版，统一。公式背景为白底。对基因表达数据进行分类、 聚类等数据分析之前， 往往需要进行预处理， 包括对丢失数 据进行填补、清除不完整的数据或合并重复数据等数据清洗，根据分析的目的进行数据过滤，以及针对分析方法选择适当的数据转换等预处理方法。数据清洗是数据分析前必须进行的一项工作，对于基因表达数据，目的是去除表达水平是负值或很小的数据、或者明显的噪声数据（单个异常大或小的尖峰信号 ），同时处理缺失数据。微阵列实验得到的数据一般是经过归一化处理的，每个点的信号强度是前景信号减去背景信号，因此有时会出现负值或很小的值，显然负值是没有生物学意义的。对于

24、这些数据点， 通过数据过滤步骤可以置为缺失或赋予统一的数值，例如对于寡核苷酸芯片数据，将低于100的数据全部设置为100。微阵列表达数据由于实验条件和芯片的因素，检测得到的信号 强度往往与细胞中实际的mRNA丰度之间没有对应关系，因此，通常是采用两个条件下的信号强度的比值，例如在 cDNA微阵列双色实验中，最后得到的往往是Ratio值。而寡核苷酸单色实验的结果是信号强度，然而在处理一组数据时，也往往选择一个样本作为对照样本，将实验数据转换成 Ratio值。在计算Ratio值时，如果参考样本的信号强度很小，就可能得 到很大的 Ratio。如果一个基因谱中仅仅存在单个特别大的Ratio值，称之为异

25、常数据点（outlier），这往往是由于噪声造成的。对于这个异常数据点，必须进行去除。数据的缺失对 于某些后续数据分析方法（例如层次式聚类和PCA ）来说有着非常大的影响，甚至是致命性的，这时必须采取相应的方法。一种方法是直接过滤掉这些存在缺失数据项的行向量或列 向量。另一种方法是设定阈值， 计算一个基因表达谱中的缺失项数目，如果达到该阈值，则将该基因表达谱从数据矩阵M中删除；如果没有达到阈值但存在缺失项，对这些缺失项可以进行插值。以0代替或用基因表达谱的平均值或中值进行代替，这些方法比较简单，但是否与真实值接近，很难进行评估。较为复杂和可靠的方法是，分析基因表达谱的模式， 从中得到相邻数据点

26、之间的关系，根据这种关系，利用相邻数据点估算得到缺失值。这种方法类似于k近邻方法，需要有足够的完整的模式来发现有缺失值的相邻模式，需要有足够的值来确定它们的邻居。在细胞中，基因表达有时空特异性，在某一条件下，发生表达的基因占基因总数的少部分，而大多数基因仅维持基础转录或不转录，转录本丰度很小，因此微阵列实验得到的数据矩阵中存在大量的基因表达谱曲线是平坦的，即基因表达水平变化很小。对于这些基因，往往不是生物学家所关心的，而它们的存在，却会大大增加数据分析的复杂性，而且会对一些分析方法的结果有干扰。 对这些数据进行过滤是非常有必要的，可以给出一定的比例， 使存在的基因占总数的多少， 这是与分析目的

27、相密切相关的，例如是分析细胞周期， 可以多保留一些基因，而对于肿瘤特异基因表达谱分析，可以少保留一点基因。 过滤这些基因所采用的标准有：基因表达谱中最大值与最小值的差；标准差；均方根；绝对值大于阈值的数据个数等。根据分析的对象和目的，可以选择以上一个或多个标准，确定阈值，来选择基 因表达谱。基因表达谱数据经过过滤，在进行聚类分析等操作前，往往还需要进行数据转换。数据 变换是将数据转换为适合数据挖掘的形式，可以根据需要构造出新的数据属性以帮助理解分析数据的特点，或者将数据规范化， 使之落在一个特定的数据区间中。 因此，数据转换包括 对数转换和标准化两个过程。许多DNA微阵列实验的结果是测量样本与

28、对照样本间信号强度的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为底)，这三

29、个数据点分别为 0、1.0、-1.0,上调2倍与下调2倍是关于0对 称的。因此，对数转换可以使小于1的值变大，大于1的值变小，从而使它们关于0对称化， 这种变换是否反映了一定的生物学意义，能更直观的了解基因的上调或下调的幅度？尚没有定论，但是对于大多数基因表达数据分析过程，都是在log空间中进行的。数据点图8.1表达数据的 Ratio和log2(Ratio)表示数据的标准化是将所有的数据转换到同一个范围内， 这样做的好处是方便比较和计算相 关系数，缺点是在标准差接近 0的时候，会产生大的噪声，这也是首先要进行数据过滤的一 个重要理由。数据标准化按如下公式进行，XijXij -XiN_(Xj -

30、Xi)2j壬8.2_NXi =N a Xj8.3j 1通过标准化，使得每个基因表达谱的平均值为0,标准差为1。如果要求所有的数据在0,1之间，还需要进行如下转换X =(X - Xmin ) /(xmax - xmin )xmin = min为公2,XnXmax = maxxi,X2,，Xn 8.4而要求数据满足a,b，则变换如下:-(b-a)(x8.5xaXmax - xmin还有一种数据标准化方法是数据的中心化。对于来自细胞系的大量肿瘤样本与一个共同的对照样本比较，对于每一个基因，都有一系列的Ratio值，相对于对照样本中那个基因的表达水平。因为对照样本通常对实验没有什么帮助，对照样本中的基

31、因表达量是独立于分析的。这样可以通过调整每一个基因的数值来反映系列观察值的变化，例如平均值或者中值。 这就是平均值/中值中心化，中心化可以减少参考样本的影响。中心化数据同样可以用于去 除某些类型的偏差。许多双色荧光杂交实验的结果没有校正Ratio值的系统偏差，它们是由于RNA数量差异，标记效率和图像获取参数所造成的。这样的偏差对于所有的基因与一个 固定数值的Ratio有放大的效应。在log空间的平均值和中值中心化有校正这种偏差的效果。 数据中心化是基于这样的一种假设，在特定的实验中，基因的平均值期望比率是1.0（在log空间中为0）。通常，更多的是使用中值中心化。目前对数据预处理这种策略的作用

32、还不是很清楚，还没有人进行系统的研究， 提供有说服力的证据来帮助研究人员针对特定的任务选择特定的数据预处理的策略和方法。在具体应用时，往往是根据分析目的和个人经验选择不同的方法。8.3基因表达差异的显著性分析在检测基因表达的微阵列实验中，有很大一部分是比较实验，目的是比较两个条件下的基因表达差异，从中识别出与条件相关的特异性基因，例如识别肿瘤特异性基因、 药物特异响应基因等。为了提高实验的可靠性，对于两个条件，往往有两个以上的重复实验，但是由 于微阵列实验仍然很昂贵，不可能重复足够的次数来满足实验数据分析的要求，因此需要采用一些比较复杂的方法来分析这些数据。对这些表达数据的分析目的就是要识别在

33、两个条件下有显著表达差异的基因。 何谓显著表达差异？通常是指一个基因在两个条件中表达水平的 检测值在排除实验、 检测等因素外，达到一定的差异，具有统计学意义，同时也具有生物学 意义。分析方法有三类，一类称之为倍数分析， 简单估计在两个条件中每个基因的表达水平 的比值，设定阈值得到表达差异显著的基因；第二类方法是估计表达差异的置信度，采用的方法是t检验和方差分析；第三类是建模的方法，通过确定两个条件下的模型参数是否相同 来判断表达差异的显著性，例如贝叶斯方法。8.3.1倍数分析早期基于CDNA微阵列技术的比较实验，用倍数来分析基因表达水平差异，即计算两 个条件下的表达水平的Ratio值。用Xgi

34、表示基因g在条件i中的表达水平测量值，因此,rg =Xg1/Xg2表示基因g在条件1和2中的表达水平比率。对于cDNA微阵列实验，两个条件的样本同时与同一微阵列进行杂交实验，得到的是成对数据，对于每次实验得到的数据计算rg，最后计算多次重复实验的平均rg。而对于寡核苷酸芯片，首先分别计算两个样本的重复微阵列实验的归一化表达水平的平均值，然后计算其比率。当rg=l时，基因g的表达水平没有改变，而 rg<1或rg>1意味着基因g在两个条件下存在表达差异，特别是山<1表示基因在条件1是下调的，而rg>l，表示在条件1是上调的。在具体应用中，如果一个基因的平均表达水平在两个条件

35、下的变化超过一个常数，典型的常数是2，即rg>2或<1/2，则认为该基因的表达差异是显著的。然而，对表达数据仔细考察后提出，这样简单的2倍法并不能产生最优的结果，因为因子2在不同的表达水平上有相当不同的显著性。对于低表达水平的基因，其信噪比太低，用2倍法作为判断条件太宽松，而对于高表达基因， 条件又太苛刻，往往小于2就有生物学意义上的表达差异显著性。在具体应用中，并没有明确的阈值，往往根据分析的具体要求由数据分析者自行确定，目前也没有关于这类简单的倍数分析方法的假阳性率和假阴性率的深入研究报道。832 t 检验对于两个条件下的多次重复实验， 为了判断基因的表达差异是否具有显著性，

36、在应用中 较多的采用假设检验，包括两个条件下的 t检验和多个条件下的方差分析( ANOVA)，这里 仅仅介绍t检验。零假设为H0：Ug1 =Ug2，与之对应的备选假设是H! ： Ug! = Ug2。 t统计量的计算公式如下:8.6t _xg1 -xg2gJsgj/m +Sg22 / n2_nj_2其中 Xgi =為 Xgij /口，Sgij 1n.1 二(Xgij ni 一1 j 总-Xgi)2,ni为某一条件下的重复实验次数。如果t超过了某个由给定的置信水平确定的阈值，就拒绝零假设，即认为基因g在两个条件下的表达差异是显著的。因为在t检验中，两个总体平均值之间的距离被经验得到的标准差归一化，

37、可以克服固定倍数阈值方法的一些缺点。然而，对于微阵列数据的t检验的基本问题是，即使用当前的高通量检测技术，实验仍然是花费很大或者实验过程很冗长，重复次数ni经常较小，n产2、3的小样本仍然非常普通。由于样本量小，导致总体方差被严重低估，得到的t值就较大，因此会导致较高的假发现率(FDR, False Discovery Rate)，即通过t检验得到的结果中表达差异不显著的基因数目较多。这样，需要更好的分析方法来克服这 些缺点。在t假设检验中，经常使用的显著性水平是p=0.01，其意思是在零假设正确的情况，进行100次抽样，会有1次错误地拒绝了零假设。对于微阵列实验，检测的基因数目巨大，如果微阵

38、列上有10000个基因，采用p=0.01，将会有100个基因是由于偶然性而被错误认为是有表达差异显著的。 这个数目已经对后续的生物学分析可能会产生很大的干扰，从而导致t检验分析结果的不可靠或失去意义。为了解决这个问题，可以对 t检验进行了改进，降低由于分母上方差小而带来的错误， 因此对t检验的计算公式修改如下：, Xg1 Xg2Sg 二 ap n1Xgii Xgi2 亠二 n2Xg2i -Xg228.78.8a = (1/ n11/n2)/(n n 2 -2)8.9变化可以较高。为了保证 dg独立于基因表达，在分母上增加S0， S0的选择疋最小化 d g的变化系数。通过对设计的一组对照样本的分

39、析，可以确定阈值，dg大于阈值的基因被认为是假设dg的分布是独立于基因表达水平的。因为较低的表达水平会使 sg的值较小，导致dg的表达差异显著的。8.3.3贝叶斯分析由于微阵列数据噪声大、波动大，而且在大量数据的背后还有很多相关变量不能被观察 到，因此贝叶斯方法可以用来分析微阵列表达数据。贝叶斯分析可以简单描述如下：P(M | D) =P(D | M )P(M )/P(D)8.10其中，P(M|D)表示由观测数据集 D得到参数化模型 M = M (w)为真的概率，称为后验概率；P(M)称为先验概率，表示在没有得到任何数据之前所估计的模型M为真的概率。P(D|M)是指似然度，表示从模型 M得到一

40、个观测数据集 D的概率。贝叶斯推断是通过参数估计和 模型选择来实现任务的，最常用的方法是最大后验概率(MAP)估计和最大似然(ML)估计。在用贝叶斯方法分析表达数据时， 首先假设在给定条件下， 一个基因的表达水平测量值 是独立的，并满足正态分布。根据经验， 这一假设是合理的， 特别是表达水平的对数大致服 从对数正态分布。对于重复实验，也可以引入伽玛分布、高斯/伽玛混合分布等。一个基因在一种条件下的表达测量值可以用一个正态分布N(x； ",二2)来建模。对每个基因和每个条件，都对应有一个双参数模型w = C12)，似然函数可以由下式给出：P(D | 点2) pi N%; ,2) 乂仟2

41、)"2©" z白心 8.11ii取遍所有的重复测量，C表示归一化常数。似然度取决于充分统计量n、x和s2。先验概率分布 P(,;2)的选择有几种，一般采用共扼先验分布。先验分布的四个超参数构成向量二=(%, o'.0,；0 )P(*2) =C(；2)32W由叱8.122对于微阵列数据，采用一个和二 相互不独立的先验分布很有意义。经过一些代数运 算，可以推导出后验分布具有与先验分布相同的函数形式:2 2 2 28.13P(),二 | D,：)二 N()； J,二 / n)】(匚；'n,；n)其中Jn8.14222-on. 2八 0匚0 (n -1)s

42、- (x - Jo)九0 + n后验分布后验分布的参数以一种合理的方式将先验分布的信息和数据信息结合了起来。是贝叶斯分析的基本对象，它包含了"和二2所有可能取值的相关信息，可以通过多种方法 进行估计。对于两个条件下的每个基因的表达测量值都与两个模型相关，根据模型参数的不同可以判断基因的表达差异是否显著。贝叶斯方法部分地克服了由于实验重复次数少造成的t检验的缺陷，如果实验次数很少(2或3次)，贝叶斯方法比t检验法效果好，如果有 5次重复实 验，两者结果相似。在重复次数为2时，贝叶斯方法的假阳性率会提高。8.4基因表达谱聚类分析对于基因表达谱数据的分析是目前生物信息学的研究热点和难点。转

43、化为数学问题，分析任务是从数据矩阵M中找出显著性结构，结构类型包括全局模型(model)和局部结构(pattern)，这些结构可以对基因或者条件形成亚类，结合生物学解释，这些亚类是细胞周期 特异性的或者是肿瘤类型特异性的。因此，对基因表达谱数据的分析是数据挖掘问题，所采用的方法包括通过可视化进行探索性数据分析(Exploratory Data Analysis )、描述建模(descriptive modeling)、分类、聚类和回归等。目前，基因表达谱分析所采用的方法主要是聚类，其目的就是将基因或条件分组，对于条件之间存在时间依赖关系的基因表达谱，更多的是对基因进行分组；而对于不同来源的肿瘤

44、样本，首先是确定特征基因，然后获取这些特征基因的表达谱数据子集，对条件进行分类，从而可以获取与肿瘤类型相关的特征，这些特征也称为基因组指纹(genomic fingerprint)或签名(sig nature)，它们可以作为肿瘤的临床诊断标准或药效的评价标准。因为从数学的角 度，对行或对列进行聚类所用的方法是一样的，不同的是对列进行分析时维度较大，样本数较小，结合生物学背景知识，通常需要降维，这个过程也称特征基因的选择过程，可以通过行向量的一些统计特征，例如方差、信息增益等，也可以用遗传算法等优化搜索算法来获取 特征基因。因此，以下部分主要是针对基因分组来介绍聚类算法。从数学的角度，聚类得到的

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

46、中，通过连坐(guilt-by-association)证据，根据一个简单的有相似表达模式的基因聚类，可以指派未知基因的功能。聚类分析是模式识别和数据挖掘中普遍使用的一种方法，是基于数据的知识发现的有效方法，特别适用于模式分类数不知道的情况。从机器学习的角度来看，聚类分析可以分为两类，即有监督聚类和无监督聚类。在有监督聚类中，对于基因表达谱数据，首先对基因进行功能分类，然后选择这些基因的表达谱作为训练集，通过有监督学习获取每一功能类的特征模式，利用这些特征模式，可以构建分类器，对未知功能的基因根据表达谱进行分类，从而能够对未知功能的基因预测功能，常用的方法有支持向量机、人工神经网络、朴素贝叶斯

47、分类、多元logistic回归等。无监督聚类不需要任何先验领域知识，根据数学特征提取分类标准，对数据进行分类，这种数学特征的例子有统计平均值、相关系数、协方差矩阵的本征值及本征向量等。无监督聚类在基因表达数据分析中应用的很多，主要有层次式聚类、k平均、自组织特征映射网络、谱聚类等。本节将对在基因表达数据分析中常用的聚类方法进行 具体介绍。841相似性度量函数对基因表达谱进行聚类分析之前，必须首先确定反映不同基因表达谱相似程度的度量函 数，根据该函数可以将具有相似程度高的基因聚为一类。在实际计算中，还可以用距离代替相似的概念，相似性度量被转化为两个基因表达谱之间的距离。距离越小，表达模式越相近，

48、反之，则表达模式差异大。常见的相似性度量有距离、点积、相关系数( correlation coefficient )、互信息(mutual in formation )等。假设两个基因表达谱分别为 X= (xi,x 2,x m)和丫= (yi,y 2,y m),距离 函数d(X , Y)必须满足如下条件：d(X , Y)仝 0d(X , Y) = d(Y , X) d(X , Y) = 0 if X = Y d(X , Y) W d(X, Z) + d(Z , Y)欧氏距离(Euclidean distanee)是一个通常采用的距离定义，与测度的选择无关，它是 在m维空间中两个点之间的真正距离

49、，两个基因表达谱之间的欧氏距离计算公式如下：1 |l_mD(X,Y)二以-yj28-15m Y y相关系数也是常用的相似性度量函数，计算公式如下：m8.16S(X,Y)二計(宁)(牛)i=d(Gi 'Goffseti 4,m8.17其中，GOffset是G的各分量的均值，G是标准方差。用上述两种相似性度量，可以找出表达模型相同或者变化趋势相同的基因，如图8.2（a）、（b）所示。欧氏距离、相关系数可以反映基因之间的共表达关系，两个基因表达谱间的距离小于或相关系数大于某个给定的阈值，就可以认为它们之间是共表达的。距离和相关系数之间存在相关性，在具体应用时，可以根据需要进行转换，例如 d=

50、1-r。距离和相关系数反映的都是基因表达谱之间的相似性，这种相似性反映了基因的共表达行为，而基因的行为是复杂的，它们之间存在调控和被调控的关系，或者存在调控链，例如基因A调控B, B调控C,调控还有正性调控和负性调控之分。对于这些调控关系，它们的 表达谱往往是不相似的，或者存在时延、或者存在反相，而基因表达的幅度也可能不相等。如何从数据中发现这些复杂的基因关系呢？互信息可能是一种有用的度量指标，其定义如下mH(X) p(Xi)log2P(xJi38.18MI（X,Y）是向量X和Y的互信息，H （X是X的熵。（c）两个基因的调控制输入一样，但是调控结果不一样，甚至相反。MI（ X, Y）=H（X

51、）+H（ Y -H（ X，Y）NMI （ X，Y）=MI （X，Y）/maxH （ X），H（ Y）8.19NMI称为归一化互信息，它独立于单个信息熵，抓住了模式上的相似性。互信息聚类分 析，没有规则上的约束。不象欧氏距离，MI不仅能确定负相关和非线性相关，同样可以反应正线性相关。因此，基于MI的聚类可能共享输入，但是对这些输入有不同的动力学过程。 例如A,B得到C输入，A是上调，B是下调。目前，还没有理论来指导如何选择最好的相似性度量，也许一个“正确”的距离在表达模式空间是不存在的，选择依赖于我们要问的问题。已知的共调控的基因在不同生物和基因 调控网络中标准集可能有利于发现理论上的相似性度量

52、。842无监督聚类方法（？）对于基因表达谱的聚类问题，由于目前对基因表达的系统行为了解得不全面，没有聚类的先验知识，所以通常采用无监督聚类方法。在基因表达数据分析方面，层次式聚类、k均值、自组织映射神经网络是应用中的常用方法。下面主要介绍这几种常用的聚类方法，并简单介绍一些其它方法。& 4。 2。1 .简单聚类假设有G个基因，它们的表达谱分别用向量表示为XX2.，,XG.。令任意一个基因的表达向量为第一个聚类的中心z1 ,假设选择zX1。然后计算X2.到z1的距离D1，如果D21大于给定的阈值 T,则说明X2.不属于第一类，应该分到另外的类。在这种情况下，建立一个新的聚类中心Z2 =

53、X2.。如果8小于阈值T,则将X2.分到第一类中。接着处理其它基因，在处理第i个基因时，首先计算该基因的表达谱与现有各类中心的距离，假设与第j类的距离D最小，并且 D<T,则将基因i分配到第j类；否则生成一个新类，该类的中心为第 i 个基因的表达向量。简单聚类算法的结果与第一个聚类中心的选择、基因的顺序、阈值T以及基因表达谱在其空间的分布有关。该方法对于给定的一组表达数据模式进行初步分类提供了一种快速的算 法。& 4。2。2 .层次聚类法层次聚类法，在统计分析中也称为系统聚类法，原理与算法与第六章所介绍的系统发生树连锁构造方法类似，所不同的只是将所分析的数据由生物分子序列换成了这

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

55、以被认为是一种局部最优解的方法。& 4。2。3. K均值聚类K均值聚类在数据划分上不考虑类的分层结构问题，该算法使待聚类的所有样本到聚类中心的距离平方和最小，这是在误差平方和准则的基础上得到的。K均值聚类算法的基本过程如下：(1) 任意选取K个基因表达向量作为初始聚类中心Zi, Z2,Zk，在没有先验知识的情况下，一般选择前 K个基因；(2) 反复迭代计算。在第I次迭代过程中，如果| X-Zj(1)|：| X-Zj(1)|(i=1,2,K,i j)，则将X所代表的基因归于第j类。按照上述办法处理所 有的基因；(3) 经过上述处理，聚类可能发生变化，因此需要重新计算K个新聚类中心：1Zj

56、(l 1)Xj=1,2，,K8.20N j Xfj(l)其中fj(l)为第I次迭代中第j个聚类的基因集合，N为该集合中基因的个数。(4)对于所有的聚类中心，如果Zj(l+1)= Z(l)(j=1,2,K)，则迭代结束，得到最后的聚类结果；否则转第 2步，继续进行迭代计算。聚类中心的个数 K、初始聚类中心的选择、基因排列的顺序以及基因表达数据的分布影 响聚类的结果，当基因表达模式类别之间分离较远时，该算法可以取得令人满意的聚类分析结果。& 4。2。4 .自组织映射神经网络人工神经网络技术在模式识别方面有着独特的优势，在生物信息学中的应用也非常广泛，如基因识别、蛋白质结构预测等。神经网络能