《计算生物学讲座》PPT课件.ppt

计算生物学讲座 李伍举 军事医学科学院基础医学研究所 * 基本概念 主要类容 数据库 序列比较 RNA二级结构预测 外源基因高效表达数学模型 蛋白质结构预测 进化树构建 基因预测：coding region, noncoding RNA genes 基因表达谱分析 计算工具：数据挖掘、模式识别等 BioSun软件介绍 What is computational biology Computational biology is a new field of research which develops models and software implementations for computational problems in molecular biology, biotechnology, and genetics. Such means are needed both in the basic research and in the industrial applications of biotechnology. Computer-aided DNA sequencing, sequence comparison, prediction of protein structures, docking of molecules, and the interpretation of electron microscopy or NMR data are examples of typical computationally intensive tasks. Bioinformatics refers to the development and use of the (molecular) biological databases 计算生物学定义 计算生物学是一门概念性学科，以生物信 息 为基础，以计算为工具，解决生物学问题。 与生物信息学的定义类似，只是侧重点有 所 不同。计算生物学侧重于计算与问题，通过计 算 解决问题；生物信息学侧重于数据的管理与数 据 库的构建。 组成部分表现形式创新点相关背景 数据数据库数据 数据库管 理系统 算法程序高效算法 数学、物理 、计算机 解决问题理论问题 提出新的 问题 生物学 计算生物学/生物信息学定义解析 数据库 目前，各式各样的生物数据库很多， 比较著名的有美国NCBI提供的Genbank， 欧洲生物信息学研究所的EBI，日本国家 遗 传学研究所的DDBJ等。 每年的Nucleic Acid Res.数据库专刊 了解目前的生物信息资源 为计算生物学研究提供基础 序列比较 计算分子生物学中基本技术、非常重要 多方面应用：蛋白质结构预测、RNA二级结构 预测、进化树构建、进化谱构建、序列模式构 建、基因功能预测和基因预测等。 比较形式多样：两个序列比较、多个序列比较 ；全局比较、局部比较；最优比较、次优比 较；扫描数据库等。 流行的比较：blast，参见NCBI,EBI网页 目前的难点：多序列的最优比较 序列比较中的基本参数 空格罚分：单独空格、延伸空格 Wx=g+rx 基本的相似性矩阵：残基与残基之间的 相似性，碱基与碱基之间相似性。最简 单的就是单位矩阵。 Needleman-Wunsch算法 A T T G C T T G 0001000 T 0110011 T 0110011 C 0000100 T 0110011 T 0110011 Seq1=ATTGCTT Seq2=GTTCTTSeq1=ATTGCTT Seq2=GTTCTT A TT G C TT G 543 4210 T454 3221 T344 3221 C 202 2 2310 T122 1121 T011 0011 ATTGCTTATTGCTT : : : : GTT CTTGTT CTT RNA二级结构预测 为什么要研究RNA二级结构预测 tRNA的转运氨基酸功能 核酶的催化功能 5和3端的非翻译区(UTR) 的结构 对基因表达的影响 RNA二级结构中的基本结 构 螺旋区(helical regions):稳定 发卡环(hairpin loop):不稳定 内部环(interior loop):不稳定 膨胀圈(bulge loop):不稳定 多分支环(multibranch loop):不稳 定 评价指标:自由能 RNA二级结构自由能计算 RNA二级结构自由能=各个基本结构的 自由能之和 基本结构的自由能由自由能计分系统确 定 具体计算方法见： 军事医学科学院院刊，1995，19：293 自由能计分系统 Salser自由能系统 Turner自由能系统 /zukerm/rna/ 碱基对之间的堆积能+环区的不稳定能 自由能计算 5-AC-3 -1.8 5-CG-3 -3.4 3-UG-5 3-GC-5 5-GU-3 -1.8 H(4,10)=4.4 3-CA-5 TotalG=-1.8-3.4-1.8+4.4=-2.6 (Kcal/Mol) RNA二级结构预测方法分类 Dynamic programming algorithms Kinetic folding algorithms Genetic algorithms Comparative methods Higgs. Quarterly Reviews of Biophysics 33:199-253,2000 Dynamic programming algorithms 基于自由能计分系统，以最低自由能 为目标函数，首先采用迭代方法求出所有 可能子片断的二级结构自由能，再用回归 方法求出RNA序列的最低自由能结构。 Maximum matching model (Nussinov et al. 1980) Minimum free energy (Zuker, NAR 1981) Zukers Minimum free energy V(i,j)表示i 与j配对时的子片断最低自由能;W(i,j)表示不 论i 与j配对与否时的子片断最低自由能；基于W，寻找 多分枝环，基于V，寻找分枝环的结构。 Kinetic folding algorithms 鉴于RNA的二级结构最终是由若干个螺旋 区组合而成，所以RNA二级结构的预测目标就 是想办法找出有关的螺旋区。动力学折叠算 法的主要思想就是沿着自由能降低的折叠方 向，以最低自由能为目标函数，去模拟RNA二 级结构。 如果不考虑所有可能的折叠路径，尽 管 动力学折叠算法以最低自由能为目标去模 拟RNA二级结构，并不能保证最终获得的 结 构为最低自由能结构。 Genetic algorithms 遗传算法，本质上是一种模仿生物进化的 优 化算法，对一个复杂问题，如果没有数学上严 格 的最优解获得办法，可用遗传算法来获得比较 优 化的解，但并没有严格的数学证明，保证获得 的 解是最优解，从这个意义上讲，是Monte-Carlo 模拟方法的特殊情况。 对RNA二级结构预测来说，由于已有 动 态编程算法来求RNA的最低自由能结构， 所 以，可视遗传算法为动力学折叠算法的一 种。 Comparative methods 基本假定：来自不同物种的RNA序列如Phe- tRNA序列，如果它们拥有相同的功能，那么 这些序列就应该拥有相同或类似的结构。因 此，必须有多个相关的序列，方可使用此方 法来研究RNA二级结构预测。 收集来自不同物种的具有相同功能的RNA序列 多序列比较 识别所有可能的螺旋区，并按一定的计分方法排序 RNA二级结构的组装 Juan(JMB,1999,289:935 只要有5个序列即可。 排序的指标：自由能和保守 性；Likelihood-ratio test. 将排序后的螺旋 区依次加入到 RNA二级结构中 去，直至形成一 个稳定结构 比较方法预测RNA二级结构流程 两种最低自由能预测方法比 较 Zuker的动态规划算法与Pipas的螺旋区组合算法均 可以用来求最低自由能结构，在多数情况下，结果是一 致的，导致不一致的原因有： 在Pipas算法中，通常只考虑全长的螺旋区，而对其 子螺旋区不予考虑，在有些情况下，可能是其子螺 旋区在最终的最低自由能结构中。 在Zuker算法中，是以单点配对为基础来求最低自由 能结构，而在最终的最低自由能结构中，要将之过 滤掉。 哪一个结构最好 无论是改进的Zuker算法（MFOLD），还是 螺 旋区组合类算法，对一个特定的RNA序列来说， 均可以求出多个RNA二级结构，如果不考虑实验 数据，究竟哪一个结构比较合理？目前，一般 认 为由比较方法求出的结构比较合理。但如果没 有 相关的序列数据，如何寻找合理的结构？ 我们的工作 螺旋区 每个螺旋区H(S,E,L)用三个参数表示： S:螺旋区起点；E:螺旋区终点; L:螺旋区长度 一级螺旋区概念 设有一个RNA序列，长度为N,Hk（Sk,Ek,Lk） (k=1,2,n)为区间i,j （1ijN）上的n个螺旋区，并满足下列关系 ： iS1E1S2E2SkEkSnEnj 则称这n个螺旋区为区间i,j上的一级螺旋区。 螺旋区之间的关系 RNA二级结构预测 根据上述概念，提出了下列算法： 基于螺旋区随机堆积的RNA二级结构预测 生物物理学报，12:213-218；1996 Prediction of RNA secondary tructure based on helical regions distribution Bioinformatics 14(8):700-706, 1998） RNA二级结构预测系统构建 生物化学与生物物 理进展，23:449-453；1996。 中心思想 RNA分子在溶液中可以有多种结构与之 对应，但每种结构出现的频率不同，有的结 构出现频率较大，起主导作用。为了求出主 导结构，以最低自由能为目标函数，采用非 决定性Monte-Carlo模拟方法获得一定数目的 二级结构，然后进行统计分析。 RNA二级结构预测的主要步骤 给定RNA序列 求出所有可能的螺旋区列表 进行随机堆积，获得一定数目的RNA二级结构 统计每个结构的出现频率 出现频率0.4，获 得主导结构 出现频率0.4，求每 个螺旋区的出现频率 ，然后反复迭代，最 终获得主导结构。 Phe-tRNA的三叶草结构 Phe-tRNA的最低自由能结构 自由能：-20.8Kcal/Mol 基于1161tRNA序列的RS, HD和MFold比较 主要结论 RS和HD:不容许螺旋区端点的GU配对 RS方法优于HD 如果螺旋区端点的GU配对容许,T3.0;如 果螺旋区端点的GU配对不容许,T25 螺旋区最佳长度为3 bp RS的预测精度为54.65% (523/957)，HD 为52.14 (499/957)，Zuker方法为32.92 (315/957) 。 BJRNAFold程序 鉴于RS与HD方法较慢，为此，基于 上 述比较获得的优化参数，对Zuker方法进 行 了改进，新的程序命名为BJRNAFold 对1139tRNA序列来说，BJRNAFold优 于RS方法。 对其它家簇来说，BJRNAFold与MFold 预测结果相当。 外源基因高效表达 数学模型构建 实现外源基因在原核系统或真核系统 的高效表达具有重要的理论和实际意义。 如一些重要的细胞因子，在体内含量甚 微，单靠提取方法，难以获得足量的细胞 因子，从而阻碍了结构与功能关系研究， 更谈不上临床的大规模应用。 目前，有一些指导性原则可用于外源 基因的高效表达设计，如利用表达系统的 优势密码子替换外源基因中的稀有密码子 或改变RNA二级结构以提高表达水平，但 是，这些原则都是定性的，其次，这些因 素是分别考虑的，没有综合考虑这些因素 的共同作用。 当时，国内的情况是：国内科学家自 己组建的原核高效表达载体pBV220，自构 建以来，得到了广泛应用，利用该载体已 成功实现了多种细胞因子等外源基因的高 效表达，但也有表达水平比较低的情况， 给我们提供了丰富的数据资源。所以，我 们以之为基础，来定量研究外源基因表达 水平与其他因素之间的关系。 数据收集 Table 1 Related data of 22 foreign genes carried by pBV220 vector No. Name Level Class D Ref. No. Name Level Class D Ref. 1 PCG12 1.96 1 9 3 12 PCG14 23.41 2 11 3 2 PCG11 5.66 1 15 3 13 HIFN 24 2 5 1 3 HIL4_CD 5-10 1 6 4 14 MIL4 25-30 2 5 9 4 HIV1 8 1 5 5 15 PCG18 27.20 2 8 3 5 RB 10-15 1 6 6 16 HIL6_T 28 2 5 10 6 PCG16 11.32 1 11 3 17 PCG13 28.93 2 7 3 7 PCG17 11.92 1 6 3 18 HBV 30.41 2 8 11 8 HCV_NS3 14 1 8 7 19 PCG15 30.78 2 10 3 9 NAPIL8 18.5 1 5 8 20 HIL4 30-40 2 6 12 10 HIL2 20 2 5 1 21 GMCSF_W 40 2 6 13 11 GMCSF_Z 20 2 5 2 22 HIL6 71 2 5 14 方法 RNA二级结构预测:采用基于螺旋区随机 堆积的RNA二级结构预测方法。 密码子偏性:采用CAI指标。 判别分析:采用Bayes判别分析方法。 RNA二级结构与高(低)表达 关系 原核启动子含有两个保守区域：-10与-35区,当然也包 含 了Gold所证实的富含核糖体结合位点信息的-2521这个区 域, 那究竟以哪个区域为标准进行计算呢 ?另外, 3端二级结构 又 用哪个区域进行计算呢? 为了确定与高(低)表达具有显著性统计学意 义的区域，我们围绕起始密码子AUG和终止密 码 子TAA附近，随机选取了多个片段，并运用基 于 螺旋区随机堆积的RNA二级结构预测方法，分 析 了每个片段的二级结构，然后运用判别分析方 法 考察每个片段二级结构自由能与高 (低)表达之 关系，结果发现三个区间 21,89,23,95和 28,93的自由能与高(低)表达具有显著的统计 学意义,结果见下表。 表: 区间选取与高(低)表达之间关系 区间 F值 概率 判别符合率 21,89 9.7195 0.0012 90.5% 23,95 10.4115 0.0009 86.4% 28,93 7.2998 0.0044 73.5% 5和3端联合判别函数 其中以区间21,89判别效果最好,由此得到5和 3端的联合判别函数为: LES=-10.8036-0.4732*G5-1.8649*G3 (1) HES=-17.1970+0.1559*G5-2.6214*G3 (2) E=HES-LES 上式中G5表示5端-3039区域(即区间21,89) 的二级结构自由能,G3表示3端30 -39区域 的 二级结构自由能,判别准则是:当E0时,判断 样品为低表达;当E0时,判断样品为高表达。 5端判别函数 5端判别函数: LES5=-2.4489-0.9339*G5 (3) HES5=-0.6790-0.4918*G5 (4) 统计量F=5.8383,P0.02540.05,判别符 合率为:59.1% 3端判别函数: LES3=-0.0535*G32-1.0971*G3- 6.7433 (5) HES3=-0.2130*G32-5.6572*G3-37.9814 (6) 统计量F=9.4318, P0.006,判别符合率为 68.2%。 3端判别函数 5和3端条件 要实现外源基因的高效表达,根据判别准 则, 下式必须成立; LES5HES5 LES3HES3 由上面的两个不等式可得到G5,G3的取值 范围为: G5-4.0 (Kcal/mol) (7) -17.21G3-11.38 (kcal/mol) (8) SD序列、ATG与TAA在二级结构 中位置与高(低)表达关系 表：SD序列、ATG与TAA在二级结构中参与配对的碱基数目 Class AUG SD TAA Class AUG SD TAA 1 3 0 0 2 1 1 0 1 0 3 0 2 0 0 0 1 1 3 3 2 1 0 0 1 0 3 0 2 3 0 0 1 0 0 3 2 0 0 0 1 0 3 0 2 3 0 0 1 2 0 0 2 0 0 0 1 3 0 0 2 1 0 0 1 0 0 0 2 1 3 0 2 0 2 3 2 1 0 0 2 0 4 2 2 3 3 0 通过SAS软件分析,AUG,SD与TAA 在 二级结构中参与配对的碱基数目与高(低) 表 达无显著关系。不过，由于所选数据均是 表达数据，通过分析上表可以看出:AUG, SD与TAA最好为0。 局部密码子偏性与高(低)表达关 系 为了探讨5与3端局部密码子偏性 与高(低)表达关系,利用Goldkey软件 计算下列值:X1,X2,.，X30；Y1， Y2，.，Y30。并且运用判别分析方法 考察Xi与Yj的所有可能900种组合与高(低) 表达的关系,结果发现:有8种组合,P值较小, 详见下表。 表：5和3端局部密码子偏性联合作用与高(低)表达关系 5偏性 3偏性 F值 概率 判别符合率 X1 Y3 5.4572 0.0362 73.3% X3 Y3 3.0462 0.0851 73.3% X4 Y3 3.0686 0.0839 73.3% X9 Y3 3.2001 0.0769 73.3% X10 Y2 3.1687 0.0785 80.0% X10 Y3 3.4554 0.0653 66.7% X11 Y3 3.0782 0.0834 66.7% X12 Y3 2.9887 0.0885 66.7% 由上表可以看出: P值均在0.05左右摆动,以Y3的 P值最小(因X1=1), 并得到下列判别方程: LEC3=-0.44192+30.73696*Y3 (9) HEC3=-2.35737+70.99093*Y3 (10) P0.0362,由判别准则可得到外源基因高效表达 条件是: Y 3 0.04758,即要求外源基因3端的 3个密码子(包括TAA)必须是大肠杆菌的优势密码 。 SD序列与起始密码子ATG之间 碱基数(D)与高(低)表达关系 由于外源基因与载体连接后,距离D就是定值,因 此,只能根据表 1中的数值作判别分析,得到的判 别函数为: LED=-4.4464+1.1273*D (11) HED=-3.1267+0.9453*D (12) 统计量F=1.2325,P0.2801,判别符合率为59.1%, 因此,外源基因的高(低)表达与D值无显著关系。 由表1可以看出,D值在83范围内较为合适。 一般判别函数构建与回顾性分 析 通过综合分析，最终得出X6,G5与G3联合作用 的 判别函数: LESC=-13.6401+12.8546*X6-0.3612*G5-2.0377*G3 (13) HESC=-21.8205+16.4293*X6+0.2990*G5-2.8423*G3 (14) 统计量F=6.8386，P0.0029，下表为22个外源 基因利用方程(13)和(14)的判别情况。从表中可 以看出:22个外源基因中只有第1个被误判。 通 过分析有关数据可知：将第1个样品误判的原因 是5端自由能太高了。 表：22个外源基因回顾性分析 No. LESC HESC DC OC No. LESC HESC DC OC 1 18.4134 19.8864 2 1 12 19.1358 19.2885 2 2 2 19.5331 18.9596 1 1 13 28.6516 36.1672 2 2 3 6.6124 2.3394 1 1 14 14.4364 15.8918 2 2 4 12.1912 10.8165 1 1 15 18.9552 19.4380 2 2 5 4.4814 0.8667 1 1 16 13.8560 14.7747 2 2 6 19.6776 18.8400 1 1 17 18.5578 19.7668 2 2 7 18.1501 15.0622 1 1 18 22.2221 26.1855 2 2 8 18.4046 14.9546 1 1 19 18.9552 19.4380 2 2 9 5.2970 1.7473 1 1 20 19.4498 20.2457 2 2 10 18.3370 22.4037 2 2 21 18.6168 22.9754 2 2 11 26.5447 27.6698 2 2 22 18.0869 19.4216 2 2 一般判别函数验证 表：4个外源基因判别结果 No. Name G5 G3 X6 Level OC DC 1 HCV_L -4.8 -11.4 0.22918 10 1 1 2 HCV_W -4.8 -12.8 0.22918 11 1 1 3 JSNA 1.0 -14.0 0.58962 60 2 2 4 TNFA -6.5 -14.9 0.26419 30 2 2 有关论文 李伍举，吴加金：pBV220载体中外源基因表达水平 定量分析 病毒学报，13:126-133；1997。 李伍举，吴加金：pBV220载体中外源基因二级结构 与表达水平关系 生物技术通讯,7:149-151,1996。 Li Wuju, Wu Jiajin: GeneDn: for high-level expression design of heterologous genes in a prokaryotic system Bioinformatics 14:884-885,1998。 模型的实验验证 进行Ricin-A链在E.coli 中的高效表达设 计，结果获得了高效表达（表达水平20%） 裴吴红 沈倍奋 李伍举等：细胞与分子免 疫学杂志，1998，14（1）：33 人FKBP12在E.coli中的高效表达设计，结 果获得了高效表达（表达水平20%） 裴武红 胡美茹 李伍举等：中国生物化学 与分子生物学报， 2000，16（3）：322 人SCF在E.coli 中的高效表达设计，结果获得了 高效表达（表达水平30%） 洪海燕等，待发表。 宋晓国等构建的高效原核融合表达载体pBVIL1， 在保持5和3端结构的情况下，改变编码区内 部 结构以实现不同基因的融合，按照我们的数学 模 型，这些基因的表达水平均应在20%以上，结果 得到了实验验证。细胞与分子免疫学杂志， 2001,17(3):231 与大连医科大学合作，进行人NMDA受体靶片断在 E.coli中的高效表达设计，结果获得了高效表达， 表达水平从未改造前的6%上升到29%。军事医学科 学院院刊，2002，第三期。 其它情况：有两例设计与预期不符，一例是关于人 Insulin，一例是关于NATO基因，它们的表达水平 均在15%左右，而没有达到预期的20%以上，推测其 原因，有一些因素我们可能没有考虑到：如表达序 列的长短，被表达蛋白质的稳定性，编码区内部的 稀有密码子串联作用等。 其它多例实验数据的验证：表达水平均较低，通过 我们分析，均找出了原因。 思考：表达模型 对原核中的其它载体？ 对真核(Yeast system)中的载体？ 我们的目标：对某一个特定的外源基因， 能够在实验之前,就能从理论上知道该 基因的表达水平。 分子生物学中的其它实验能否达到计算 机模拟? 蛋白质结构预测 二级结构预测： 神经网络方法、Chou-Fasman、 Garnier等，目前最好的方法是Garnier方法 、预测精 度在78%左右，不到80%。 三级结构预测： 从头预测，同源模建等 应用：抗原表位分析、小分子药物设计等方 面 / 进化树构建 收集一组相关的序列 多序列比较，计算序列之间相似性 基于相似性矩阵，可分别采用下列方法构建 Maximum parsimony method (Minimum evolution method) Distance methods (Neighbor-joining methods) Maximum likelihood approach /phylip.html 基于全基因组的物种进化树构 建 从统计学上讲：对物种进化树的真实情况 是无知的，因此，只能通过样本情况推测 总体情况： 基于不同基因家簇给出的进化树进行叠加 基于全基因含量物种进化树构建 基于若干保守的蛋白家簇联合 基因预测 编码区预测: GeneScan，Geneie， nc RNA 基因预测 ncRNA基因预测的专用方法 ncRNA基因预测的通用方法 ncRNA基因预测的专用方 法 主要是为识别某个特定的ncRNA基因家簇的 新成员而设计的一类方法，一般采用启发式算 法，即根据特定的ncRNA基因家簇的一级结构和 二 级结构特征，发现一些规则，然后根据这些规则扫 描基因组序列，并从基因组序列中发现符合这些规 则的基因片段，这些基因片段即为该ncRNA基因 家 簇的可能新成员，如能通过比较基因组学在相近物 种的基因组中发现类似的保守片段，便可进一步确 证这些新成员的身份，当然，最终要通过实验来验 证。 一级结构和二级结构特征分 析 为了寻找某一特定的ncRNA基因家簇的一级 结 构和二级结构特征，通常对该家簇的所有已知成 员 （序列）或部分已知成员（序列）进行多序列比 较，然后，基于多序列比较形式提取该ncRNA基 因 家簇的一级结构和二级结构信息，通常有三种方 法：加权矩阵即通常的频数矩阵、模式分析和 Eddy提出的用于RNA序列分析的共变化模型。 ncRNA基因预测的通用方 法 基本的出发点是希望能够找出一种通用 的预测方法，不依赖于某个特定的ncRNA 基 因家簇信息，将ncRNA基因从基因组序列 中 识别出来。目前已发展了下列四种方法：碱 基组成方法、神经网络方法、比较基因组学 方法和转录起始位点与终止位点预测法。 碱基组成方法 基本设想是将一个基因组中ncRNA基因序 列当作信号，基因组序列当作背景，然后利用 ncRNA基因的碱基组成与基因组中的碱基组成 的 差别将ncRNA基因识别出来，这种差别越大， 越 有利于ncRNA基因的识别。 神经网络方法 基本思想是将整个基因组序列分为三个部分：编 码蛋白质的基因部分、ncRNA的基因部分（第一类） 和 其它的非编码基因间区（第二类），并进一步假设第二 类中只有小部分含有ncRNA基因，然后以大肠杆菌为 例，运用神经网络方法探讨了第一类与第二类的区分问 题。鉴于第一类的长度要远远小于第二类，于是从第二 类中随机选出一个与第一类大小相当的部分来进行训 练，然后，运用获得的神经网络预测第二类中的其余部 分，从而识别出第二类中的含有ncRNA基因的序列片 段，为了预测结果可靠，将上述过程多次重复进行。 比较基因组学方法 基本假设是ncRNA基因在相近物种的基因组中， 不 仅一级结构有一定的保守性，更重要的是其二级结构也 非常保守，根据这个设想，Rivas和Eddy提出了一个 识 别ncRNA基因的自动方法，其相应程序为QRNA。该 方 法的核心部分是提出了三个概率模型，它们分别是编码 区模型COD、ncRNA基因模型RNA和零假设模型 OTH，然后基于序列的比较形式，采用这三种模型分 别 计算 Bayes后验概率，并根据数值的大小判定被比较 序 列片段是编码区或ncRNA基因或其它的序列。 转录起始位点与终止位点预测 法 在已经注释的基因组序列基础上，通过预测转录起 始 位点与终止位点来识别ncRNA基因的一种方法，主要包 含 四个步骤。首先，从待分析的基因组中找出“空白”的 基 因间区，在这空白的基因间区中，不包含任何已注释的 基 因（ORF，tRNA，rRNA等）；其次，进行转录起始位 点 预测；再次，进行转录终止位点预测；最后，采用序列 比 较方法进行保守性分析,找出可能的ncRNA基因。但对 真 核生物来说,转录起始位点的预测并非易事。 小结 上述四种方法均可用于ncRNA基因的 识别问题,但是预测结果的可靠性均不及蛋 白质编码区预测算法。通过对大肠杆菌基 因组序列测试表明，从精度与被证实的 ncRNA基因的数目来说，以比较基因组学 方法最好。 基于基因表达谱的生物信息 学 计算机辅助寡核苷酸微阵列探针设计 基于基因表达谱的分类系统构建基于基因表达谱的分类系统构建 基于基因表达谱的分型系统构建基于基因表达谱的分型系统构建 基于基因表达谱的调控网络构建 综合性数据库的构建等 李伍举：基因表达谱的生物信息学 军事医学科学院院刊,26:73, 2002 为什么要研究基于基因表达谱 的样本分型与分类问题 人类基因组测序模式生物基因组测序。 Sequence Database 这些序列的功能是什么？ 基因芯 片技术 蛋白芯 片技术 其它高通 量技术？ 基因芯片技术 Gene 1 Gene 2 Gene 3 Gene m Tissues 基因表达谱 基于基因表达谱，可以同时考察在特 定生理或病理过程中细胞内基因群的动态表 达水平，从而将基因的活动状态比较完整地 展现出来，使研究人员能够在基因组水平上 以系统的、全局的观念去研究生命现象及其 本质。但是，从数据处理的角度来看，通过 基因芯片实验直接获得的是一个基因表达 谱，相关的数据分析和挖掘已经无法仅仅通 过简单的计算来进行，辅助数据分析系统的 建立已经成为基因芯片技术的必要工具。 基于基因表达谱的样本分型研 究 聚类分析：谱系聚类、 K-平均值 变量选择：标准差、 比值法 样本分型研究 样本分型方案的评价 聚类分析中的变量选择 Lukashin（2001）采用下列规则来对酵母细胞的周 期数据进行聚类：要求所有17个时间点上的表达数据的 绝对值大于或等于100，并且表达水平至少有2.5倍的变化 Welsh（2001）采用的基因变量选择方法是要求每个 基因表达谱的标准差大于或等于250，最终选出1243个基 因对上皮细胞卵巢癌表达数据进行聚类分析 Perou(1999)则采用R/G比值方法，要求在26个上皮或 乳癌样本中，至少有三个样本的表达水平在3以上，最终 选出1247个基因来对样本进行聚类。然而，上述基因变 量的选择方法并不具有普遍适用的特点 为了较好地解决聚类分析中的变量选择问题， Xing（2001）提出了一个通用的用于聚类分析的基因 变 量选择方法CLIFF，其主要思想为：在假定样本表型未 知情况下，采用某种聚类方法获得初步的样本分型方 案，然后以之为基础，采用监督学习方法，选择一定数 目的基因来对样本进行进一步的分型，该过程反复迭 代，直到最后样本的分型方案没有变化为止。在Xing 方 案的监督学习中，基因变量的个数是靠经验来确定的， 最终的样本分型方案与选中的基因变量个数有关，其 次，并没有从理论上给出最佳的样本分型个数。因此， 到目前为止，基于基因表达谱的样本分型问题仍没有得 到很好解决。 我们在详细研究他人算法的基础上，尝试 了多种变量选择方法后发现，以变异系数和t- 检 验为基础进行基因变量选择是一个行之有效的 方 法，通过这种基因变量选择方法与谱系聚类、 K- 平均值方法和自组图方法的整合，构建了基于 基 因表达谱的样本分型的整合系统SamCluster， 较 好地解决了聚类分析中的变量选择问题。 基于基因表达谱的样本分类研 究 分类方法：Fisher、 距离判别 变量选择：逐步优化、 t-检验 样本分类研究 样本分类方案的评价 主要应用：基于基因表 达谱的肿瘤分类研究。 目前，样本类型预测研究主要应用于肿瘤分类，利 用基因表达谱在肿瘤与正常组织中的差异对肿瘤进行分 类与诊断已形成共识。当前的肿瘤分类技术高度依赖于 病理学工作者对肿瘤组织的主观判断，而基于基因芯片 技术，即使一些组织没有显著变化，利用基因表达谱也 可以对之做出早期诊断；另外，特别重要的一点是可以 根据基因表达谱的变化来区分形态学上相似的肿瘤，这 样对肿瘤类型的精确识别有助于制定配套的最佳方案， 从而达到增加疗效、降低毒性的目的；另外，基于基因 表达谱的肿瘤分类对肿瘤发生机制的理解以及征服这些 肿瘤提供了重要思路。 究竟采用多少个基因以及如何选择这些基 因，特别是对多类肿瘤情况如何选择这些基因， 到目前为止，还没有定论；不过，采用最少的基 因达到最高的预测精度将是追求的目标。为此， 以我们提出的分类稳定性概念为基础，探讨了 Fisher线性判别分析方法与Monte-carlo模拟和逐 步优化等基因变量选择方法的整合，以便对某个 特定的基因表达谱，寻找最好的基因集合，这将 对肿瘤发生机制与药物作用机制的理解具有一定 意义，并且，所找出的最佳基因或基因集合可用 于肿瘤诊断与药物靶基因识别等方面，为进一步 的分子生物学实验提供思路。 第一部分 基于基因表达谱的 样本分型系统构建 数据与方法 四个数据集说明 为了说明样本分型系统Samcluster的 性 能，并与其它分型方法进行比较 ，主要 运 用四个数据集：COLON、 OVARIAN 、 LEUKEMIA72、LEUKEMIA38 。 1、COLON：该基因表达谱包含2000个 基 因和62例样本，其中有22例正常组织 ， 其余40例样本为结肠癌组织。 2、LEUKEMIA72: 该基因表达谱包含 6817 个基因和72例样本，其中有47例急性 淋 巴细胞白血病（ALL），其余25例为 急 性粒细胞白血病（AML）。 3、LEUKEMIA38: 该基因表达谱包含6817个 基因和38例样本，在38例样本中，有 27例急 性淋巴细胞白血病（ALL）和11 例急性粒细 胞白血病（AML）, 此数据由LEUKEMIA72 提取而 来。 4、OVARIAN: 该基因表达谱包含7129个基因 和 36例样本，在36例样本中，有27例卵巢癌组 织、5例正常组织和4例恶性的卵巢癌细胞系 。 变异系数计算（CV） 设G=gij(mn) 表示基 因表达矩阵，其中m和n 分别表示基因和样本的个 数，gij表示第i个基因在 第j个样本中的表达水平， 对第i个基因来说，用 gi=gi1, gi2, ,gin表示 该基因在n个样本中的表 达向量。 两点特别考虑 1、如果在一个基因表达矩阵中，由于标 准 化或对数化等原因导致某个基因的最 小 表达水平(gmin)为负值，则要对该基因 的表达向量进行调整，调整的方式为 该 基因的表达向量加上gmin的绝对值，调 整后的基因向量的各个分量均为非负 值，从而便于均值与标准误的计算。 2、由于实验误差等原因，常常导致某个基因在 某个样本中的表达水平呈现异常（偏大或偏 小），为了克服这种异常对变异系数的影 响，在计算某个基因的变异系数时，不考虑 该基因向量的最大与最小分量。 CV=cv1, cv2, , cvm 用MCV和SCV表示变异系数的均值与标准误 。 t-检验 谱系聚类 聚类，就是物以类聚之意，为了刻画 样本之间的相似性并进行聚类，必须对样 本之间的相似性或距离进行定量，然后， 按特定的方式进行聚类。 1、标准化变换 标准化变换，首先对基因变量进行标 准化，然后对样本变量进行标准化，标准 化之后其均值为0，方差为1。 2、构建相关系数矩阵 构建相关系数矩阵，就是采用标准的 相关系数计算方法，计算n个样本之间的 两 两相关系数，从而获得相关系数矩阵。 3、谱系聚类 首先将n个样本视作n个类，从相关 系 数矩阵中寻找最大值元素即最相似的两个 类，并将它们合并成一类，由此原来的n 个 类便转化为n-1类，对新的n-1类继续构建 相 关系数矩阵，此过程反复进行，直到最 后，所有的n个样本在1类为止。 类与类之间的距离（或相关系数）定 义有多种方式，比较常用的有如下6种方 式：最短距离法、最长距离法、中间距离 法、重心法、类平均法和变差平方和法， 这在一般的统计学课本上均有介绍。在构 建Samcluster过程中，我们采用了类平均 法。 基于基因表达谱的样本类 型发现的整合方案 两个假设 1、要求基因在不同样本中的表达水平必须有波动。 变异系数是一个很好的指标，利用它可以刻画 不 同数据集中的各个基因的表达水平的波动情况。为了 选出用于样本类型发现的基因变量，必须要求有关基 因在各个样本中的表达水平有波动，且其变异系数大 于某个指定的阈值,如阈值太低，选出的基因变量集 合 中将包含一定数目的噪声基因变量（即不利于样本类 型发现的基因变量），如阈值太高，选出的基因变量 集合中将缺少一定数目的信号基因变量（即有利于样 本类型发现的基因变量）。为此，我们采用下列公式 来选择变异系数的阈值。 CVth = MCV+CiSCV 2、要求基因的表达水平在不同的样本分型中（指推定 的 样本分型），其表达水平的差异具有一定的显著性 。 由于在我们的研究中，总是假定两种可能的样 本 类型，因此，可用t-检验来刻画。通过第一个假设 ， 我们可以获得某个特定的基因变量集合，以之为基 础 进行聚类分析，可得到两个推定的样本类型,对此进 行 t-检验分析，可找出表达水平呈现一定显著性差异 的 一些基因（可假定P=0.01、0.05或0.1等），从而 将 表达水平不显著的一些噪声基因变量去除，在保留 的 基因变量集合基础上，进行进一步的聚类分析，此 过 程反复迭代，直到最后没有可剔除的基因变量为止 ， 最终获得一个稳定的样本分型方案。 一致样本类型构建 在固定P值的情况下，多次改变标准误的系 数因子Ci的值，对每一次改动，我们将获得一 个 推定的样本分型方案，并由此构建样本之间的 关 系矩阵Snn。我们称以此为基础进行聚类分 析 所获得的样本分型方案为样本的一致分型方案 。 为了获得最佳的样本分型方案和对应的基因集 合，考虑了样本的一致分型方案与各个推定的 样 本分型方案之间的距离，并称距离最小者为最 佳 的样本分型方案。 变异系数计算 基因表达矩阵，给定P值 选定CVth 聚类分析 两个推定的样本类型 t-检验 获得推定的样本分型方案 改变CVth，获得多个 样本分型方案 构