计算生物学讲座.ppt

计算生物学讲座,李伍举军事医学科学院基础医学研究所2019/12/15,基本概念主要类容数据库序列比较RNA二级结构预测外源基因高效表达数学模型蛋白质结构预测进化树构建基因预测：codingregion,noncodingRNAgenes基因表达谱分析计算工具：数据挖掘、模式识别等BioSun软件介绍,Whatiscomputationalbiology,Computationalbiologyisanewfieldofresearchwhichdevelopsmodelsandsoftwareimplementationsforcomputationalproblemsinmolecularbiology,biotechnology,andgenetics.Suchmeansareneededbothinthebasicresearchandintheindustrialapplicationsofbiotechnology.Computer-aidedDNAsequencing,sequencecomparison,predictionofproteinstructures,dockingofmolecules,andtheinterpretationofelectronmicroscopyorNMRdataareexamplesoftypicalcomputationallyintensivetasks.Bioinformaticsreferstothedevelopmentanduseofthe(molecular)biologicaldatabases,计算生物学定义,计算生物学是一门概念性学科，以生物信息为基础，以计算为工具，解决生物学问题。与生物信息学的定义类似，只是侧重点有所不同。计算生物学侧重于计算与问题，通过计算解决问题；生物信息学侧重于数据的管理与数据库的构建。,计算生物学/生物信息学定义解析,数据库,目前，各式各样的生物数据库很多，比较著名的有美国NCBI提供的Genbank，欧洲生物信息学研究所的EBI，日本国家遗传学研究所的DDBJ等。每年的NucleicAcidRes.数据库专刊了解目前的生物信息资源为计算生物学研究提供基础,序列比较,计算分子生物学中基本技术、非常重要多方面应用：蛋白质结构预测、RNA二级结构预测、进化树构建、进化谱构建、序列模式构建、基因功能预测和基因预测等。比较形式多样：两个序列比较、多个序列比较；全局比较、局部比较；最优比较、次优比较；扫描数据库等。流行的比较：blast，参见NCBI,EBI网页目前的难点：多序列的最优比较,序列比较中的基本参数,空格罚分：单独空格、延伸空格Wx=g+rx基本的相似性矩阵：残基与残基之间的相似性，碱基与碱基之间相似性。最简单的就是单位矩阵。,Needleman-Wunsch算法,Seq1=ATTGCTTSeq2=GTTCTT,ATTGCTT:GTTCTT,RNA二级结构预测,为什么要研究RNA二级结构预测,tRNA的转运氨基酸功能核酶的催化功能5和3端的非翻译区(UTR)的结构对基因表达的影响,RNA二级结构中的基本结构,螺旋区(helicalregions):稳定发卡环(hairpinloop):不稳定内部环(interiorloop):不稳定膨胀圈(bulgeloop):不稳定多分支环(multibranchloop):不稳定评价指标:自由能,RNA二级结构自由能计算,RNA二级结构自由能=各个基本结构的自由能之和基本结构的自由能由自由能计分系统确定具体计算方法见：军事医学科学院院刊，1995，19：293,自由能计分系统,Salser自由能系统Turner自由能系统/zukerm/rna/碱基对之间的堆积能+环区的不稳定能,自由能计算,5-AC-3-1.85-CG-3-3.43-UG-53-GC-55-GU-3-1.8H(4,10)=4.43-CA-5TotalG=-1.8-3.4-1.8+4.4=-2.6(Kcal/Mol),RNA二级结构预测方法分类,DynamicprogrammingalgorithmsKineticfoldingalgorithmsGeneticalgorithmsComparativemethodsHiggs.QuarterlyReviewsofBiophysics33:199-253,2000,Dynamicprogrammingalgorithms,基于自由能计分系统，以最低自由能为目标函数，首先采用迭代方法求出所有可能子片断的二级结构自由能，再用回归方法求出RNA序列的最低自由能结构。Maximummatchingmodel(Nussinovetal.1980)Minimumfreeenergy(Zuker,NAR1981),ZukersMinimumfreeenergy,V(i,j)表示i与j配对时的子片断最低自由能;W(i,j)表示不论i与j配对与否时的子片断最低自由能；基于W，寻找多分枝环，基于V，寻找分枝环的结构。,Kineticfoldingalgorithms,鉴于RNA的二级结构最终是由若干个螺旋区组合而成，所以RNA二级结构的预测目标就是想办法找出有关的螺旋区。动力学折叠算法的主要思想就是沿着自由能降低的折叠方向，以最低自由能为目标函数，去模拟RNA二级结构。,如果不考虑所有可能的折叠路径，尽管动力学折叠算法以最低自由能为目标去模拟RNA二级结构，并不能保证最终获得的结构为最低自由能结构。,Geneticalgorithms,遗传算法，本质上是一种模仿生物进化的优化算法，对一个复杂问题，如果没有数学上严格的最优解获得办法，可用遗传算法来获得比较优化的解，但并没有严格的数学证明，保证获得的解是最优解，从这个意义上讲，是Monte-Carlo模拟方法的特殊情况。,对RNA二级结构预测来说，由于已有动态编程算法来求RNA的最低自由能结构，所以，可视遗传算法为动力学折叠算法的一种。,Comparativemethods,基本假定：来自不同物种的RNA序列如Phe-tRNA序列，如果它们拥有相同的功能，那么这些序列就应该拥有相同或类似的结构。因此，必须有多个相关的序列，方可使用此方法来研究RNA二级结构预测。,收集来自不同物种的具有相同功能的RNA序列,多序列比较,识别所有可能的螺旋区，并按一定的计分方法排序,RNA二级结构的组装,Juan(JMB,1999,289:935只要有5个序列即可。,排序的指标：自由能和保守性；Likelihood-ratiotest.,将排序后的螺旋区依次加入到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；1996PredictionofRNAsecondarytructurebasedonhelicalregionsdistributionBioinformatics14(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螺旋区最佳长度为3bpRS的预测精度为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，自构建以来，得到了广泛应用，利用该载体已成功实现了多种细胞因子等外源基因的高效表达，但也有表达水平比较低的情况，给我们提供了丰富的数据资源。所以，我们以之为基础，来定量研究外源基因表达水平与其他因素之间的关系。,数据收集,Table1Relateddataof22foreigngenescarriedbypBV220vectorNo.NameLevelClassDRef.No.NameLevelClassDRef.1PCG121.9619312PCG1423.4121132PCG115.66115313HIFN242513HIL4_CD5-1016414MIL425-302594HIV1815515PCG1827.202835RB10-1516616HIL6_T2825106PCG1611.32111317PCG1328.932737PCG1711.9216318HBV30.4128118HCV_NS31418719PCG1530.7821039NAPIL818.515820HIL430-40261210HIL22025121GMCSF_W40261311GMCSF_Z2025222HIL6712514,方法,RNA二级结构预测:采用基于螺旋区随机堆积的RNA二级结构预测方法。密码子偏性:采用CAI指标。判别分析:采用Bayes判别分析方法。,RNA二级结构与高(低)表达关系,原核启动子含有两个保守区域：-10与-35区,当然也包含了Gold所证实的富含核糖体结合位点信息的-2521这个区域,那究竟以哪个区域为标准进行计算呢?另外,3端二级结构又用哪个区域进行计算呢?,为了确定与高(低)表达具有显著性统计学意义的区域，我们围绕起始密码子AUG和终止密码子TAA附近，随机选取了多个片段，并运用基于螺旋区随机堆积的RNA二级结构预测方法，分析了每个片段的二级结构，然后运用判别分析方法考察每个片段二级结构自由能与高(低)表达之关系，结果发现三个区间21,89,23,95和28,93的自由能与高(低)表达具有显著的统计学意义,结果见下表。,表:区间选取与高(低)表达之间关系区间F值概率判别符合率21,899.71950.001290.5%23,9510.41150.000986.4%28,937.29980.004473.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端条件,要实现外源基因的高效表达,根据判别准则,下式必须成立;LES5HES5LES3HES3由上面的两个不等式可得到G5,G3的取值范围为:G5-4.0(Kcal/mol)(7)-17.21G3-11.38(kcal/mol)(8),SD序列、ATG与TAA在二级结构中位置与高(低)表达关系,表：SD序列、ATG与TAA在二级结构中参与配对的碱基数目ClassAUGSDTAAClassAUGSDTAA1300211010302000113321001030230010032000103023001200200013002100100021302023210020422330,通过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值概率判别符合率X1Y35.45720.036273.3%X3Y33.04620.085173.3%X4Y33.06860.083973.3%X9Y33.20010.076973.3%X10Y23.16870.078580.0%X10Y33.45540.065366.7%X11Y33.07820.083466.7%X12Y32.98870.088566.7%,由上表可以看出:P值均在0.05左右摆动,以Y3的P值最小(因X1=1),并得到下列判别方程:LEC3=-0.44192+30.73696*Y3(9)HEC3=-2.35737+70.99093*Y3(10)P0.0362,由判别准则可得到外源基因高效表达条件是:Y30.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.LESCHESCDCOCNo.LESCHESCDCOC118.413419.8864211219.135819.288522219.533118.9596111328.651636.16722236.61242.3394111414.436415.891822412.191210.8165111518.955219.43802254.48140.8667111613.856014.774722619.677618.8400111718.557819.766822718.150115.0622111822.222126.185522818.404614.9546111918.955219.43802295.29701.7473112019.449820.2457221018.337022.4037222118.616822.9754221126.544727.6698222218.086919.421622,一般判别函数验证,表：4个外源基因判别结果No.NameG5G3X6LevelOCDC1HCV_L-4.8-11.40.2291810112HCV_W-4.8-12.80.2291811113JSNA1.0-14.00.5896260224TNFA-6.5-14.90.264193022,有关论文,李伍举，吴加金：pBV220载体中外源基因表达水平定量分析病毒学报，13:126-133；1997。李伍举，吴加金：pBV220载体中外源基因二级结构与表达水平关系生物技术通讯,7:149-151,1996。LiWuju,WuJiajin:GeneDn:forhigh-levelexpressiondesignofheterologousgenesinaprokaryoticsystemBioinformatics14: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%以上，推测其原因，有一些因素我们可能没有考虑到：如表达序列的长短，被表达蛋白质的稳定性，编码区内部的稀有密码子串联作用等。其它多例实验数据的验证：表达水平均较低，通过我们分析，均找出了原因。,思考：表达模型,对原核中的其它载体？对真核(Yeastsystem)中的载体？我们的目标：对某一个特定的外源基因，能够在实验之前,就能从理论上知道该基因的表达水平。分子生物学中的其它实验能否达到计算机模拟?,蛋白质结构预测,二级结构预测：神经网络方法、Chou-Fasman、Garnier等，目前最好的方法是Garnier方法、预测精度在78%左右，不到80%。三级结构预测：从头预测，同源模建等应用：抗原表位分析、小分子药物设计等方面/,进化树构建,收集一组相关的序列多序列比较，计算序列之间相似性基于相似性矩阵，可分别采用下列方法构建Maximumparsimonymethod(Minimumevolutionmethod)Distancemethods(Neighbor-joiningmethods)Maximumlikelihoodapproach/phylip.html,基于全基因组的物种进化树构建,从统计学上讲：对物种进化树的真实情况是无知的，因此，只能通过样本情况推测总体情况：基于不同基因家簇给出的进化树进行叠加基于全基因含量物种进化树构建基于若干保守的蛋白家簇联合,基因预测,编码区预测:GeneScan，Geneie，ncRNA基因预测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,为什么要研究基于基因表达谱的样本分型与分类问题,人类基因组测序,模式生物基因组测序,。,SequenceDatabase这些序列的功能是什么？,基因芯片技术,蛋白芯片技术,其它高通量技术？,基因芯片技术,Gene1,Gene2,Gene3,Genem,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。我们称以此为基础进行聚类分析所获得的样本分型方案为样本的一致分型方案。为了获得最佳的样本分型方案和对应的基因集合，考虑了样本的一致分型方案与各个推定的样本分型方案之间的距离，并称距离最小者为最佳的样本分型方案。,程序设计,运用Matlab程序设计语言，构建了基于基因表达谱的样本类型发现系统Samcluster，,结果,基于原始基因表达数据的样本分型没有进行基因变量选择,基于整合系统的样本分型研究进行基因变量选择,表1：基于Samcluster系统的结肠癌基因表达谱COLON的样本分型情况,图1：基于Samcluster系统的结肠癌基因表达谱COLON的样本分型情况,作者方法变量选择方法基于临床分型误判数,Alon双向聚类t-检验是8Xiong谱系聚类Fisher方法是5Samcluster谱系聚类CV计算与t-检验否6Samcluster谱系聚类CV计算与t-检验是5,表2：基于Samcluster系统的白血病基因表达谱LEUKEMIA72的样本分型情况,图2：基于Samcluster系统的白血病基因表达谱LEUKEMIA72的样本分型情况,系统方法变量选择方法基于临床分型误判数,CLIFF标准化分割Bayer误差是3聚类算法InformationgainrankingMarkovblanketfilteringCLIFF标准化分割Bayer误差否？聚类算法InformationgainrankingMarkovblanketfilteringSamcluster谱系聚类CV计算与t-检验否2Samcluster谱系聚类CV计算与t-检验是1,表3：基于Samcluster系统的白血病基因表达谱LEUKEMIA38的样本分型情况,图3：基于Samcluster系统的白血病基因表达谱LEUKEMIA38的样本分型情况,对LEUKEMIA38数据集来说，最好的样本分型方案给出的分型精度94.7%(36/38)。如果结合样本的临床表型来选择样本的分型方案，则所有样本都被正确分型。在Golub的结果中，他们使用自组图方法SOM对LEUKEMIA38进行分型研究，在38个样本分为两个类型时，有4个样本被误判（1AML，3ALLs）。因此，对LEUKEMIA38数据集来说，Samcluster系统要优于SOM。,从上表的D值分布来看，当CVi位于区间0.5，1.0时，样本分型往往达到最佳效果，因此，我们可以认为，当CVi值太高时，有许多有助于样本分型的基因变量将被舍弃；当CVi值太低时，有许多不利于样本分型的基因变量将被引进。而且，三个P值（0.01，0.05和0.1）均可以用来作为基因表达水平是否具有显著性差异的检验指标，其中以0.05和0.1的效果较好。,图4：基于Samcluster系统的卵巢癌基因表达谱OVARIAN的样本分型情况,讨论,小结,通过这一部分工作，我们构建了基于基因表达谱的样本类型自动识别系统Samcluster，其中心思想是将基于变异系数和t-检验的变量选择方法整合到常用的谱系聚类算法中，并提出了一致样本类型的概念，在此基础上获得了基于基因表达谱的样本的最佳分型方案。通过对4个基因表达数据的测试，结果表明，获得的基因水平上的样本分型与样本的临床分型具有较高的一致性。因此，这一工作对以基因芯片技术为基础的相关研究（如中药现代化等）具有较大意义。,与CLIFF的基于基因表达谱的样本分型系统比较,比较项目SamclusterCLIFF结果,聚类方法：谱系聚类标准化的分割聚类非监督学习CV计算Bayererror不能确定最佳值监督学习t-检验Informationgainranking不能确定最佳值Markovblanketfiltering不能确定最佳值基因表达分布无混合Gaussian分布一致分型有无参数敏感性不敏感敏感测试数据集4个1个LEUKEMIA721个误判3个误判,第二部分,基于基因表达谱的样本分类系统构建,数据与方法,1、基因表达数据集：COLON2、Fisher线性判别分析方法,类内离差矩阵W和类间离差矩阵B,M类样本的基因表达矩阵,计算第l类样本和总体样本的平均值,基因变量选择方法,1、主成分分析2、t-检验3、邻近关系分析中P值4、MonteCarlo模拟方法5、逐步优化方法6、全局优化方法,分类精度的稳定性分析,在基于基因表达谱的样本分类研究中，用不同的基因变量选择方法所选出的基因集合常常不一致，那究竟选择哪一个基因集合？即使用同一种方法，常常选出具有相同或相近分类精度的基因集合，那如何选择较好的基因集合？另外，在分类研究中，常常将整个样本按一定比例分为训练组与试验组，为此，我们不得不考虑样本分配情况对分类精度的影响。基于这些考虑，我们提出了分类稳定性概念。,程序设计,结果,、基因变量个数与分类精度之间的关系,主成分个数与分类精度关系,基于t-检验的基因个数与分类精度之间关系,基于P值的基因集合的分类精度,基于Monte-Carlo模拟的基因集合的分类精度,基于逐步优化的基因集合的分类精度,多种基因变量选择方法的比较分析,基于主成分的样本分类，尽管利用了基因表达谱的全部信息，但实际分类精度并不高；另外，为了对新样本进行判别，必须要知道2000个基因的表达水平。因此，无论是实际应用，还是理论分析，此种方法价值不大。另外,通过上述结果可以看出：无论是采用T-检验，还是P值方法，随着基因数目的增加，分类精度呈下降趋势，而用较少的基因，其分类精度却不是很高。因此，在基于基因表达谱的样本(肿瘤)分类研究中，这些方法不是很好。,Monte-Carlo模拟方法虽然运用较少的基因而达到较高的分类精度，但是，仍存在下列缺点：计算机运行时间较长，不能保证所得的基因集合是最佳的，特别严重的是，每次计算所得的基因集合并不固定，因此，此方法也不宜使用。然而，逐步优化方法与全组合方法却能克服上述缺点，既保证运用较少的基因获得较高的分类精度，又能保证每次运行获得相同的基因集合，因此，值得推广使用。在实际应用中，可以采用下列策略：如果基因表达谱中基因数目较少，采用全组合方法，寻求最佳基因集合；如基因数目很多，应用逐步优化方法寻找最优基因集合。,小结,通过这一部分工作，我们采用Fisher线性判别分析方法和多种变量选择方法构建了基于基因表达谱的样本分类系统Tclass,并以结肠癌基因表达谱数据为例说明了Tclass系统的功能,结果表明：运用较少的基因就可以达到较高的分类精度。另外，应用所设计的分类系统，还分析了由文献和合作实验室提供的乳腺癌，白血病，皮肤癌和肾癌的基因表达谱数据，均获得了满意结果。,方法比较,目前，有多种判别分析方法，如距离判别分析、Bayes判别分析、支持向量机等均可用于基于基因表达谱的样本分类研究，事实上，在Tclass系统中，已包含了距离判别分析和Bayes判别分析，通过多个数据的运行，结果表明：Fisher线性判别分析、距离判别分析和Bayes判别分析的效果类似。这里，我们选择Fisher线性判别分析，而没有选择支持向量机，主要是基于下列考虑：1、由于要进行基因变量选择和分类精度的稳定性分析，要求速度比较快，而Fisher线性判别分析就比较好地满足这一要求；2、Fisher线性判别分析可同时用于两类或多类样本的情形，逻辑比较简单，只需考虑新的样本个体与其它各个样本空间的