基因表达数据分析

第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（/）是一个国际性的成员联盟，参与人员包括生物学家、计算机科学家、数据分析学家。

它的目标是促进由功能基因组学和蛋白组学研究产生的微阵列数据的共享。

当前集中于建立微阵列数据注释和交换的标准，推动微阵列数据库建设和相关软件来实现这些标准，促进高质量的、经过注释的基因表达数据在生命科学领域的共享。

该组织开发的微阵列数据标准称为
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（/geo）是由NCBI在2000年开发的一个基因表达和杂交微阵列数据仓库，同时作为获取来自不同生物体的基因表达数据的在线资源。

到2004年3月，数据仓库中包含内容605个Platforms，14391个Sample，816个Serial。

Platform是关于物理反应物的信息，例如核酸、抗体和组织微阵列和SAGE数据等的基因表达数据被接受、增加和归档作为公共数据集。

Series是关于样本集的信息，反映样本间的相关性和组织。

ArrayExpress(/arrayexpress/)是基于基因表达数据的微阵列公共知识库，目的是存储被很好注释的数据，当前包含多个基因表达数据集和与实验相关的原始图像集。

ArrayExpress数据库接受MAGE-ML格式的数据递交或者通过MIAMExpress的基于Web的数据注释和递交工具。

ArrayExpress提供一个简单的基于Web的数据查询界面，并直接与Expession Profiler数据分析工具相连，可以进行表达数据聚类，和其它类型的Web 数据发掘，并将进一步开发多个实验和数据库间的交叉查询。

ArrayExpress数据库中的数据将与所有相关的由EBI维护的或在线的数据库相联接。

斯坦福微阵列数据库(SMD，/ )是一个使用Oracle作为数据库管理软件的关系数据库。

SMD存储微阵列实验的原始、归一化数据和对应的图像文件。

自从2002年1月1日起，到现在包括85篇学术论文，超过3500个双色点样DNA微阵列的实验数据，每年增加1000个微阵列实验的数据。

另外，SMD提供数据获取、分析和可视化的界面，目前包括层次聚类和自组织映射等方法，还将加入k-平均聚类、单值分解和丢失值归纳等方法。

除了以上3个综合性的基因表达数据仓库外，还有一些专门的基因表达数据库，例如YMD (Yale Microarray Database，/microarray/)、ArrayDB （/arraydb/）、BodyMap（http://bodymap.ims.u-tokyo.ac.jp/）、ExpressDB(/ExpressDB/)、HuGE Index(Human Gene Expression Index，/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 微阵列实验，两个。