|
基因芯片筛选差异表达基因方法比较
单文娟, 童春发, 施季森
摘要: 使用计算机模拟数据和真实的芯片数据, 对8 种筛选差异表达基因的方法进行了比较分析, 旨在比较不同方法对基因芯片数据的筛选效果。模拟数据分析表明, 所使用的8 种方法对均匀分布的差异表达基因有很好的识别、检出作用。算法方面,SAM 和Wilcoxon 秩和检验方法较好; 数据分布方面, 正态分布的识别效果较好,卡方分布和指数分布的识别效果较差。杨树cDNA 芯片分析表明, SAM、Samroc 和回归模型方法相近, 而Wilcoxon 秩和检验方法与它们有较大差异。 关键词: 基因芯片; 杨树; 差异表达 基因芯片技术是随着“人类基因组计划”(Humangenome project, HGP)发展起来的一项新技术, 可广泛应用于基因序列分析、基因突变检测和多态性分析以及疾病的基因诊断等领域。基因芯片技术可以同时检测生物样本成千上万个基因的表达水平, 从而得到基因组水平的基因表达谱数据。通过对这些数据的分析, 能够了解基因的功能甚至基因之间的相互作用。为了达到这一目的, 最关键的步骤是从表达谱中筛选出潜在的差异表达基因[1], 常用的是两种不同实验条件下差异表达基因的筛选。不同的实验条件包括: 不同的生物组织、不同的试剂或药物处理、不同的发育时期、不同的实验处理等。 基因表达数据可以在3 个层次上加以分析: (1)单基因水平, 即分析基因在两种实验条件之间的差异; (2)分析功能相似的基因及基因之间的相互作用、共调控等; (3)基于基因网络及蛋白质网络的研究[2]。 本文只进行第一层次上单基因水平差异表达基因的研究。差异表达通常指一个基因在两种试验条件下表达水平值在排除实验、检测等因素后达到一定的差异, 具有统计学意义, 同时也具有生物学意义[3]。常用的统计方法包括参数方法和非参数方法两类。Newton 等[4]提出了用于分析差异表达的贝叶斯模型;Baldi 和Long 等[2]使用了贝叶斯模型方案分析正态分布的数据; Lonnstedt 和Speed [5]使用了参数贝叶斯方法分析芯片数据, Smyth 等[6]在此基础上提出了更加实用的系统模型。 以上的参数方法都假设基因芯片数据服从正态分布, 这种假设与实际情况存在差异。非参数方法直接估计检验统计量的分布情况, 不对数据的分布做任何假设, 在这一点上优于参数方法。非参数方法包括Tusher 等[7]提出的SAM,Dudoit 等[8]提出的调整p 值的T-检验及Pan等[9]提出的混合模型方法(MMM)等。Nykter 等[10]指出由于目前缺乏生物学背景的准确信息, 无法对各种算法进行有效的评价, 因此使用已知样本数据特征的计算机模拟数据(Simulateddatasets)可以解决以上难题。计算机模拟数据根据实验的需要进行设计, 是检测基因芯片数据分析方法的常用手段。模拟数据包括两部分: 差异表达基因的模拟表达值和非差异表达基因的模拟表达值, 具有两种或两种以上的实验条件(组织样本)及n 个(n≥3)生物学重复。差异表达基因一般占基因总数的5%。根据Kim 等[11]、Perelman等[12]和Shaik 等[13]等的设计思想, 以正态分布为例, 非差异表达基因在不同实验条件下其基因表达值的均值和方差都没有差异; 而差异表达基因的均值都存在差异, 方差可以相同也可以不同。 目前芯片模拟数据一般只考虑正态分布的情况,如kim 等[11]比较了T-检验、B-统计量(B-statistic)、贝叶斯T-检验3 种参数方法和SAM、Samroc、Zhao-Pan 3 种非参数方法对正态分布的筛选效果。Hunter 等[14]、Thomas等[15]、Pan 等[16]、Craig 等[17]、Giles 等[18]、Liu 等[19]、Zhao 和Pan 等[20]指出, 即使基因芯片表达数据经过预处理后仍然不满足正态分布。Neuhauser 等[21]使用了均匀分布、卡方分布、柯西分布和指数分布的数据对T-检验、Wilcoxon 秩和检验、FP 检验(Fisher-Pitman permutation test)、B 检验(Baumgartner-Wei β-Schindler test)4 种统计方法进行了分析。 本文使用8 种筛选差异表达基因方法(详见1.3)分别对4 种分布的计算机模拟数据以及真实的芯片数据进行分析。参照Neuhauser 所采用的数据分布, 本文的芯片模拟数据包括正态分布、均匀分布、卡方分布和指数分布, 以便全面分析各种方法对不同分布数据差异表达基因的识别情况。  1 材料和方法 1.1 模拟数据 使用R 软件[22]生成模拟数据, 包括4 种分布: 正态分布、均匀分布、卡方分布和指数分布。以正态分布为例, 非差异表达基因的表达值设定3 种方案, 差异表达基因的表达值在对照样本和实验样本中也分别设定3 种方案, 详细方案见表1。我们产生10 000 个基因的表达数据, 其中9 500 个为非差异表达基因数据, 500 个为差异表达基因数据。每个非差异表达基因的表达数据从已设定的3 种方案中任取一种产生, 共产生12 个数据, 前6 个数据为实验(样本)数据, 后6个为对照(样本)数据(说明: 使用R 软件模拟数据时为一次性同时产生6 个符合条件的数据集, 即一次性产生实验数据集或对照数据集)。同样, 每个差异表达基因的表达数据从已设定的3对方案中任取一对, 用表1实验栏中的分布产生6 个数据作为实验数据, 用表1对照栏中的分布产生6 个数据作为对照数据。如果一个基因的12 个数据占据一行, 那么这120 000 个数据构成一个10 000行12 列的数据表格或矩阵, 其中每列数据表示一个芯片上10 000 个基因的表达数据。该实验方案的生物学重复数为6。采用不同的统计方法对此分析, 并记录分析结果。 以上所作的模拟相当于一次实际实验, 所得结论不具有统计学意义, 我们将以上模拟重复了500次, 对每种分析方法统计500 次并分析结果, 以鉴别不同分析方法的优劣性。1.2 真实芯片数据 本文采用的真实芯片数据为Andersson等[23]发表文章中的杨树cDNA 芯片数据。该实验研究了杨树秋天叶片衰老的机制, 芯片的13 490 个克隆是由7 个cDNA 文库中的36 354 个杨树EST 序列拼接而成, 其中有12 376(92%)个功能已知。实验材料为瑞典北部生长期超过30 年的雄性白杨树Aspen(Populus tremula)叶片, 对照样本和实验样本的生物学重复数分别为3 和4。 1.3 差异表达基因筛选方法 使用以下8 种筛选差异表达基因的方法: (1)T-检验(T-test)[11],是较常用的统计方法, 用于判断某一基因在两个样本中其表达是否有显著性差异, 不要求等方差; (2)贝叶斯T-检验(Bayes T-test)[2, 24], 这种方法的最大优点是结合了一个基因样本间的差异和部分基因之间的差异信息, 即对于某一基因方差的估计, 既考虑到了基因本身样本间差异又利用了具有相似表达水平的部分基因(定义为邻近基因, Neighborhood gene)的合并方差, 且两者之间使用加权值, R 语言开发的Cyber-T 软件实现了这一统计量的应用, 网址http://www.genomics.uci.edu/;(3)SAM(Significance analysis of microarrays), 由Tusher、Tibshirani和Chu 提出[7, 11], SAM 是非参数的统计方法, 与T-检验相似, 为了避免T-检验中由于方差太小而误判某基因为差异表达基因的不足, SAM在分母上增加了一个较小的正值; (4)Zhao和Pan 方法[20], 非参数的统计方法, 其关键在于估计检验统计量的零分布(Null distribution) g Z ,并直接构建了零统计量(Null statistic) g z 。非参数方法的早期研究包括Efron 的经典贝叶斯方法[25]、Tusher的SAM[7]、Pan 的MMM 方法[9], Zhao 和Pan 指出它们存在的一个共同问题是g Z 及g z 的分子与分母之间都是非独立的, 因此提出了新的g Z 、g z 计算方法, 详见参考文献[20]; (5)Samroc方法, 由Broberg[26]提出,使用这一统计量的R 语言软件是SAG, 网址http://home.swipnet.se/pibroberg; (6) 回归模型方法(Regression modeling approach), Thomas 等[15]于2001年提出了这种鉴定两种不同实验条件差异表达基因的方法, 其给出的回归方程不仅包括筛选差异表达基因部分, 还包括数据的预处理部分, 本文只讨论筛选差异表达基因部分;(7)PaGE[27,28], 是使用perl 语言编写的用于分析基因芯片表达谱数据的生物学软件, 最初由Manduchi 等提出[29], 其界面简单, 使用方便, 该算 法经过了不断改进, 使用FDR 估计置信水平, 并对FDR的估计方法也给予了修改, 使其估计值更加准确;(8) Wilcoxon 秩和检验(Wilcoxon rank sum test), 又称Mann-Whitney 检验[30], 也是一种非参数的检验方法,该方法要比T-检验更加稳健, 更适合非正态分布的数据, 要求两样本数据分布相同。 2 结果与分析 2.1 模拟数据结果分析 由于模拟数据产生时差异表达基因占基因总数的5%, 即10 000 个基因中有500 个基因为差异表达基因, 因此判断某一方法检测差异基因的有效性主要取决于: 按统计量值排序后前500 个基因中检测出的差异表达基因的个数和平均秩。  表2 中每项包含两个数据, 如对应于正态分布和T-检验的数据为136.90、202.43, 第一个数据表示按统计量(从大到小)排序, 位于前500 个基因中的差异表达基因总数, 第二个数据表示这些差异表达基因的平均秩。例如T-检验在分析正态分布的某个模拟样本数据时, 按统计量值排序后前500 个基因中差异表达基因检测到137 个, 这样分析500 个模拟样本数据后得到的平均值为136.90;“202.43” 表示500 个模拟样本数据检测到的差异表达基因的平均秩的平均值为202.43,平均秩是指排在前500 个基因中差异表达基因的秩和除以差异表达基因个数。筛选出的差异表达基因个数越多越好; 而平均秩越低, 则说明差异表达基因的排位越靠前, 使用相应的检测方法筛选出差异表达基因的可能性越大。 对于正态分布, SAM和Bayes T-test两种方法最佳,在前500 个基因中都检测到了157 个以上的差异表达基因, 平均秩也较低。其次为Wilcoxon 秩和检验、回归模型方法、T-检验和Samroc。PaGE 的平均秩最小,但其检测到的差异表达基因个数较少。 对于均匀分布, 所有方法的平均秩都在200 左右, 差异不大, 因此可以不考虑此值对检测效果带来的影响。SAM 和PaGE 两种方法最佳, 在前500个基因中都检测到了较多的差异表达基因, 约占差异表达基因总数的80%。其次为Bayes T-test 和Wilcoxon 秩和检验, 再次为回归模型方法、T-检验和Samroc。 对于卡方分布, SAM 最佳, 检测到的差异表达基因数量较多且平均秩较低。其次为Wilcoxon 秩和检验和Bayes T-test。Wilcoxon 秩和检验检测到的差异表达基因数量最多但平均秩最高; Bayes T-test 检测到的差异表达基因数量较少但其平均秩较低。再次为Samroc、回归模型方法和T-检验。 对于指数分布, SAM 最佳, 检测到的差异表达基因数量最多且平均秩较低。Wilcoxon 秩和检验检测到的差异表达基因数量也较多, 但不足的是平均秩最高。PaGE 的平均秩非常低, 检测到的差异表达基因数目与T-test、Samroc 和回归模型方法相近, 说明PaGE 要优于它们。与SAM 相比较, PaGE 的不足是检测到的差异表达基因偏少, 这或许与其检测到的差异表达基因数目少有一定关系。 从整体上看, 均匀分布的检测结果最好,Zhao-Pan 检测到的差异表达基因最少, 却也检测到了286 个, 占全部差异表达基因的57.2%。其次为正态分布, 其余两种分布的检测结果都不佳, 检测到的差异表达基因都不超过100 个。为了充分考察筛选差异表达基因的各种方法, 本文在设计模拟数据时差异表达基因与非差异表达基因的区别较小, 但以上结果一定程度的说明八种差异表达基因筛选方法对均匀分布有很好的识别、检出作用, 对正态分布识别效果较好, 而对卡方分布和指数分布识别效果较差。从检测方法来看, SAM 最好, 其次为Wilcoxon秩和检验, 不足的是其平均秩有时较高。回归模型和T-检验对正态分布和均匀分布的检测结果非常相近。回归模型和Samroc 对卡方分布和指数分布的检测结果非常相近, 且这两种方法比较稳定, 具有较好的检测结果。Bayes T-test 比较适合检测正态分布和均匀分布的数据。PaGE 不适合检测卡方分布的数据。Zhao-Pan的检测结果平均秩一般较高, 检测到的差异表达基因数目也较少。 2.2 杨树cDNA 芯片数据结果分析 Andersson 等使用SAM 对杨树cDNA 芯片数据进行了差异表达基因的筛选, 共得到874 个差异表达基因。本文根据模拟数据分析的结果, 采用了Wilcoxon 秩和检验、回归模型和Samroc 3 种比较稳定的方法, 对其3 792 个有显著表达的基因再次做了筛选分析。我们发现3 种分析得到的差异表达基因个数都很多, 均在2 000 个以上, 于是我们根据统计量分别选取了排在前面的874 个基因。结果表明,回归模型与SAM 都检测为差异表达的基因共有649个; Samroc 与SAM 都检测为差异表达的基因共有687 个; Wilcoxon 秩和检验与SAM 都检测为差异表达的基因共有381 个; SAM、Samroc 和回归模型3种方法都检测为差异表达的基因共有622 个, 占总数的71%; SAM、Samroc、Wilcoxon 秩和检验及回归模型4 种方法都检测为差异表达的基因共有305个。以上结果说明SAM、Samroc 和回归模型方法有较多的共同点, 而Wilcoxon 秩和检验方法与它们有较大的差异。  表3 为SAM 筛选到的排在前20 位的基因在Samroc、Wilcoxon 秩和检验和回归模型3 种方法检测结果中的统计量值或P 值及其秩。从表3 可以看出, Wilcoxon 秩和检验方法有一个缺点, 由于它本身依靠排秩的方法来筛选基因, 使得某些基因的统计量值相同, 无法进一步区分它们。另外, Samroc和回归模型方法具有更大的相似性, 例如它们都将A024P46(Clone ID)排在了第一位, 将I027P15 排在了第6 位, 还将F065P41、I004P59、I059P90 都排在了百位之后。 3 讨 论 Kim 等[11]对模拟数据做分析时指出, 在其研究的所有情况下SAM 和Samroc 表现最好, 而且当样本数较小时Samroc 比SAM 更好。在本研究中SAM表现好再次得到证实, 但其要明显好于Samroc, 这或许与使用的生物学重复数不同有关, 因为Kim 所使用的小样本重复数为4, 本文使用的重复数为6。Neuhauser 等[21]指出当生物学重复数为6 及以上时,其所研究的4 种方法才能够较好的筛选差异表达基因, 这也正是本文生物学重复数为6 的原因。另外,Kim 只对正态分布的情况做了分析。他还指出,Bayes T-test 在检测小样本数据时表现最好,T-检验筛选差异表达基因的效果也不错。本文得出正态分布数据的分析结果与其相似, 只是SAM的筛选结果比Bayes T-test更好。 我们发现Bayes T-test不稳定,对卡方和指数分布的数据不敏感, 对这两种分布数据的检测结果不太理想。SAM、Samroc 和回归模型方法在分析真实芯片数据时得到了比较理想的结果, 但在分析卡方分布及指数分布的模拟数据时表现不佳, 其原因应该与真实芯片数据的分布特征有关, 因此对芯片数据的分布特征做一个简单的检测, 再有针对性地选择适合相应分布的分析方法, 或许有利于提高筛选差异表达基因的准确率。另外, SAM 在4 种分布中除了卡方分布的筛选结果没有达到最大外, 其余3 种分布均筛选到最多的差异表达基因, 这可能与SAM 通过调整△(SAM的一个统计量)的值来控制FDR 以得到更准确的结果有关, 并且可能与SAM算法的不断改进有关。根据Zimmerman 等[31]以及Blair 等 [32]的报道, Wilcoxon 秩和检验方法要好于T-检验, 本文的结果与其相同, 而且Wilcoxon 秩和检验方法还是一种比较稳定的筛选方法。 本文研究发现, 基因芯片数据经过分析后得到的差异表达基因的数量有时较多, 这使得生物学家无法从分析结果中挑选自己感兴趣的基因, 最简单的方法是根据统计量筛选适量的基因进行后续分析。另外可以使用多种方法分析, 结果取其交集, 这样可以一定程度上避免假阳性基因带给后期验证分析的麻烦。此外, 研究数据分析方法的特点有助于做好数据挖掘工作, 改进模拟数据方案有利于充分检验算法的特征,并推进基因芯片数据分析的发展。 |