序列模式识别课件

生物信息学：预测1.生物信息学最核心的问题：预测2.生物信息学工具的作用：预测3.生物信息学所有的分析：预测4.基本假设(贝叶斯的哲学理念)：我们能够通过对已知世界的观察，总结经验，并以此来预测未知世界已经存在或者即将发生的事物/事件5.在生物信息学中的应用：对现有的数据，使用合适的算法，进行训练，构建计算模型和计算工具，预测未知的现象本章内容提要1.统计学基础2.序列模式3.预测性能检验4.位点特异性打分矩阵(PSSM)5.模体发现：GibbsSampler等6.马尔科夫及隐马尔科夫模型7.模式识别的其他算法简介1.统计学基础排列组合从N个物品中取出k个物品的排列数（排序）：从N个物品中取出k个物品的组合数（不排序）：概率模型概率模型:一个能够通过不同的概率产生不同结果的模型。概率模型可以模拟或者仿真某一类型的所有事件，并且对每个事件赋予一个概率。色子模型：一个色子存在6个概率值：p1,p2,…,p6，其中，掷出i的概率为pi(i=1,2,…,6)。因此：pi≥0，且考虑三次连续的掷色子，结果为[1，6，3]，则总概率为：p1p6p3概率分布1.考虑连续变量x，例如：物体的重量。则当重量确切为1公斤时的概率，为0。2.变量的区间：P(x0≤x≤x1)3.当区间无限小->0时，上式：P(x-δx/2

≤x≤x+δx/2

)=f(x)δx4.f(x)称为概率密度函数5.因此：且二项分布1.事件只有两种可能出现的结果。例如掷硬币，正面记为“1”，反面记为“0”。2.则，掷硬币N次，有k次是1的概率为：二项分布的期望值期望值代表了随机变量的“平均”值。它是把每个可能取值乘以对应的概率，然后累加起来。期望值E(x)=μ二项分布的与标准方差标准方差描述了随机变量中具有正概率值的分散性。所有可能的值离期望值的距离的平方，再乘以对应的概率。方差VarX=σ2泊松分布1.稀有事件发生的概率：在一个连续的时间或空间中，稀有离散变量出现的概率2.N->∞,E(x)=μe=2.71828…泊松分布与二项分布的近似对于大的N及小的p值的二项分布，能够相当准确地用一个参数为μ=Np的泊松分布近似。当实验次数很多而概率很小时：二项分布~泊松分布例1：鸟枪法的覆盖率假设：需要测序的BAC长度200kbp;总共测序的序列数量：N;每次测序：500bp；每次测序的覆盖率p：500/200kbp=0.0025因此：总覆盖率μ=Np(每个点平均覆盖到的次数)k:测序能够覆盖到点X的次数。鸟枪法：覆盖率点X被覆盖k次的概率：（二项分布~泊松分布)当点X一次都不被覆盖时，k=0;此时的概率为：覆盖率vs.准确性泊松分布：例2Prof.Gene发现一条1mbp的序列上存在5个某种调控信号，该调控信号在人的基因组上平均每500kbp出现一个。那么，完全是随机产生该种情况的概率是多少？本例中，N=3.0*109bp->∞,E(x)=μ=2(1mbp)统计性显著：p-value<0.05非常显著：p-value<0.01超几何分布与二项分布的区别：不放回抽样。例：有N个球，其中红球M个，白球N-M个，每次拿出一个球再放回，总共n次，其中有m个球是红球的概率为(二项式分布)：p=M/N超几何分布(2)上例改为：有N个球，其中红球M个，白球N-M个，每次拿出一个球不放回，总共n次，其中有m个球是红球的概率为：并且，0≤m≤M<N超几何分布右尾概率上例再改为：有N个球，其中红球M个，白球N-M个，每次拿出一个球不放回，总共n次，其中有至少有m个球是红球的概率为：并且，0≤m≤M<N超几何分布左尾概率上例再改为：有N个球，其中红球M个，白球N-M个，每次拿出一个球不放回，总共n次，其中有最多有m个球是红球的概率为：并且，0≤m≤M<N超几何分布双尾概率方法一：所有出现概率<=观察表概率的概率之和方法二：双尾概率=2×min(左尾概率,右尾概率)超几何分布：例Prof.Gene从26873个人的蛋白质中预测了2264个具有某种特定功能的底物，并进行进一步的分析。其中，已知有421个人的蛋白质具有某种功能结构域D，而在预测的2264个底物中，有94个蛋白质具有结构域D。问：结构域D在2264个底物中是显著出现，显著不出现，还是随机出现？问题转化：在26873个人的蛋白质中有421个具有功能结构域D，任意取出2264个蛋白质，其中至少有94个具有功能结构域D的概率是多少？N=26873;n=2264;M=421;m=94;Fisher’sExactTest超几何分布的精确概率计算：2X2表B2：抽样B1：剩余A2：阳性A1：阴性超几何分布计算公式p-value==如上例a+b+c+d=26873c+d=2264b+d=421d=94/fisher.htmFisherExact.jarCMD下输入命令：java–jarFisherExact.jarNMnmFisher’sExactTest：再例假设，我们调查了100个学生，比较是否男生比女生更喜欢玩电子游戏。数据统计如下：玩游戏不玩游戏男生4515女生2713P-value=0.496854471943056>0.05统计性不显著！序列模式识别2.序列模式1.功能结构域，functionaldomain2.模块，BLOCK3.模体，motif4.模式，pattern/profile功能结构域1.具有完整的、独立的三级结构2.具有特定的生物学功能3.一般长度，几十到几百个氨基酸4.允许插入/缺失，即允许存在gap模块/BLOCK1.几个到几十个氨基酸2.无gap，从全局多序列比对的结果直接处理得到3.描述蛋白质家族或者一类蛋白质的序列保守性BLOCK模体/Motif1.不具有独立的三级结构2.具有特定的生物学功能：结合，修饰，细胞亚定位，维持结构，等3.长度一般几个到几十个氨基酸或者碱基；4.例如，SUMO化的序列模体：Ψ-K-X-E(Ψ:A,I,L,V,M,F,P;X:任意氨基酸)模式/Pattern/Profile1.在算法上用来描述一类功能结构域、模体或者模块的表示方式2.根据序列数据，构建的预测模型3.数据形式：正则表达式4.用来预测新的可能符合特定模式的序列5.例如，直接将Ψ-K-X-E视为SUMO化位点的，普适的“模式”，则可以预测所有包含该模式的蛋白质序列3.预测性能的计算和检验1.样本/检验数据：阳性数据(P)，阴性数据(N)a.阳性数据(P)：真实的，被实验所证实的数据b.阴性数据(N)：被实验所证明为无功能的数据2.对于预测结果的评测，定义：a.真阳性(TP):阳性数据中被预测为阳性的数据b.假阳性(FP):阴性数据中被预测为阳性的数据c.真阴性(TN):阴性数据中被预测为阴性的数据d.假阴性(FN):阳性数据中被预测为阴性的数据TPFPFNTNPositiveNegativeCutoff常用的检验指标1.灵敏度(Sensitivity,Sn)对于真实的数据，能够预测成“真”的比例是多少2.特异性(Specificity,Sp)对于阴性的数据，能够预测成“假”的比例是多少3.准确性(Accuracy,Ac)对于整个数据集(包括阳性和阴性数据)，预测总共的准确比例是多少4.马修相关系数(Mathewcorrelationcoefficient,MCC)

当阳性数据的数量与阴性数据的数量差别较大时，能够更为公平的反映预测能力，值域[-1,1]常用的检验指标(2)ROCcurveX轴：1-SpY轴：SnROC的面积越大，表明其预测能力越强预测性能的计算Self-consistencyLeave-one-outvalidationn-foldcross-validationSelf-consistency1.将训练数据当成测试数据训练数据中所有的阳性数据为测试数据中的阳性数据训练数据中所有的阴性数据为测试数据中的阴性数据2.反映当前预测工具对目前已知的数据的预测能力3.假设：根据目前已知的数据所构建的计算模型能够反映未知的数据的模式4.缺点：不能反映计算模型的稳定性Leave-one-outvalidation每次从数据集中去掉一个重新进行训练，构建预测模型，并对去除的数据进行预测。保证每个数据去掉一次n-foldcross-validation将数据集分成n组，并保证阳性数据与阴性数据的比例与原数据相同随意将n-1组作为训练数据，1组作为检验数据，计算性能重复若干次，例如，重复20次计算平均值缺点：每次计算结果有偏差预测性能及稳定性1.Self-consistency:反映检验性能（对已知数据的预测能力）2.Leave-one-outvalidation&n-foldcross-validation:反映预测系统的稳定性（对未知数据的预测能力）3.预测性能vs.检验性能a.差距较小：系统稳定b.差距过大：系统不稳定，数据过训练过训练1.根据已知数据构建的模型只能很好的适用于训练数据2.不适合用来预测未知数据3.对训练数据的微小改变对于预测性能影响过大4.预测工具过训练：只能很好的符合训练数据，而对新数据则性能很差4.位点特异性打分矩阵(1)PositionSpecificScoringMatrix(PSSM)/WeightMatrixModel(WMM)(2)对蛋白质家族进行多序列比对分析，发现结果中保守的BLOCK(3)根据BLOCK序列推导相应的PSSM(4)不考虑gap的影响(5)BLOCK长度一般在几个~几十个残基/碱基BLOCK->PSSM代表每一列二十种氨基酸矩阵中的数值：当前位置上，某种氨基酸出现的频率的log值第二种PSSM每一个位置上显示每种氨基酸或者碱基出现的频率碱基的位置四种碱基第三种PSSM每一个位置显示氨基酸/碱基出现的概率PSSM矩阵使用P(S|+)，根据阳性训练数据计算出来的概率；未知序列：ACGGTACGG背景概率选择,P(S|-)1.负样本/阴性数据的概率计算2.计算方法：A.DNA序列，四种碱基出现的频率B.蛋白质序列，20种氨基酸出现的频率OddsRatioLog-oddsRatio计算流程：滑动窗口窗口宽度9bp，依次打分设定阈值（Threshold），凡是高于阈值的预测为阳性，低于阈值的预测为阴性5.模体发现：GibbsSamplerGibbsSampler是一种Monte-Carlo类的方法

随机抽样对于输入序列，找到一个最大的似然函数GibbsSampler算法(1)1.从每条序列上随机的抽取一段序列，序列长度固定所有序列motifGibbsSampler算法(2)2.构建PSSM/权重矩阵GibbsSampler算法(3)3.随机挑选一条序列GibbsSampler算法(4)4.用构建好的PSSM对该序列上所有可能的motif进行打分(窗口滑动，每次1个氨基酸或者碱基)GibbsSampler算法(5)5.根据似然性的计算，得到似然值最大的模体，即新的motifGibbsSampler算法(6)6.更新PSSM矩阵GibbsSampler算法(7)7.反复迭代计算，直到似然性结果与PSSM不再发生变化StrongMotifACGTAGCAGibbsSampler:总结1.模体发现的一种随机算法(MonteCarlo)2.寻找次优解的算法3.根据PSSM/WMM对随机抽取的序列进行打分来调整采样，直到结果收敛4.不能够保证每次运算的结果一致：需要多运算几次，并进行比较5.对蛋白质、DNA、RNA序列模体的发现有帮助期望最大化算法1.ExpectationMaximizationAlgorithm2.已开发工具：MultipleEMforMotifElicitation(MEME)3.motif大致的位置与长度是确定的4.重点：确定motif在每条序列上的起始位置5.分为两步：Estep:估计motif起始位置的期望最大化Mstep:motif似然性的期望最大化期望最大化算法(2)1.例，假设10条序列，长度20个碱基2.进行多序列比对，大致确定motif的位置3.待找motif长度为4个碱基Motif的概率vs.背景概率1.计算motif中每个位置的碱基的概率分布2.背景概率:根据剩下的序列计算四种碱基的概率分布似然性概率值的计算似然性概率值的计算(2)计算每条序列，在不同的起始位置，其似然性的概率值Estep:起始位置估计Z值：motif在不同位置起始的几率值假设，motif在任意位置起始的概率相同，则Z值最大化，即为“最可能的起始位置”Mstep：P值最大化根据选择的最大Z值，重新计算矩阵，并计算P值最大的motif；P值最大原先的motifEM算法：迭代Gibbs&EM:总结1.基本假设：所有序列都拥有，且仅拥有一个motif2.估算两个关联的函数：Gibbs(WMM&似然性)，EM(motif起始位置，Z值&似然性)3.利用两个函数的其中之一修正另一个，采取迭代/反复计算的方法，使结果收敛4.不保证得到的结果为最优，近似算法有待解决的问题1.给定的一组序列，可能的motif仅在部分序列中出现，怎么解决？2.给定一组序列，其中存在某种motif可能在序列上出现两次以上，如何解决？6.马尔科夫及隐马尔科夫模型1870年，俄国有机化学家VladimirV.Markovnikov首次提出马尔科夫模型马尔科夫模型马尔科夫链隐马尔科夫模型VladimirV.Markovnikov马尔科夫模型马尔科夫模型：随机过程的一种，主要特点为“无后效性”，即根据当前的状态即可完全确定将来的状态马尔科夫性&马尔科夫链1.定义：对于随机变量X1,X2,X3…,这些变量的范围，即他们所有可能取值的集合，被称为“状态空间”，而Xn的值则是在时间n的状态。如果Xn+1对于过去状态的条件概率分布仅是Xn的一个函数，则符合马尔科夫性：2.具有马尔科夫性的过程称为马尔科夫过程3.时间(先后顺序)和状态都离散的马尔科夫过程称为马尔科夫链马尔科夫模型：参数估计转移概率：K-order马尔科夫模型一阶马尔科夫模型：当前位置仅依赖前一位k阶马尔科夫模型：当前位置依赖前一位，而前一位依赖前两位,…,前k-1位依赖前k位0阶马尔科夫模型：位点独立Markov&PSSM1.对真实的数据进行训练，PSSM=~0阶马尔科夫模型2.对新序列的扫描：从头至尾，每次移动1~n位（窗口滑动的方法）3.分别计算窗口内的序列，是(+)和(-)的概率，计算log-oddsratio4.设定阈值，若高于阈值，则预测为阳性另外长度不确定！起始位置不知！Markovmodels&PSSM:Notwork!!!隐马尔科夫模型(HMM)隐马尔科夫模型：

1.表示状态的可观察符号出现概率已知

2.状态之间的转移概率未知与马尔可夫模型的本质区别： 隐马模型观察到的符号并不是与状态一一对应，而是通过一组概率分布相联系ProfileHMM1.多序列比对的结果中，氨基酸之间存在的关系有匹配(M),插入(I)和缺失(D):三种状态2.HMM：三种状态之间的转换关系未知->hidden->转移概率3.每个位置上的氨基酸/碱基以及插入、缺失的频率/概率可以通过观测求得->nothidden4.模型训练：通过训练，估算转移概率例：CpG岛的HMM1.CpG岛：在人的基因组中，如果双碱基对CG出现，则C通常被甲基化。并且，甲基化的C很快会突变成T。因此基因组中CpG岛非常少。然而，在基因的起始位置，例如promotor区域，因为功能的保守性，其序列很少突变，CpG的含量能够保持在40~60%2.Howtopredict?PSSM&Markovarenotworkatall!CpG岛：HMM存在两种状态：是CpG岛(CpGIsland,I)，不是CpG岛(Genome,G)CpG岛：HMM1.Hidden:对当前未知的碱基，跳转到下一个位置，究竟是I还是G的概率，未知2.Observable:I和G中的四种碱基分布的概率能够通过实际数据的观测进行计算转移概率发散概率预测CpGIsland:Viterbi算法1.给定序列：ATCGCA,预测CpG的位置?初始概率：0.5CpGIsland:Viterbi算法(1)vATCGCAβ1C+G+A+0.1T+C-G-A-0.15T-0.5*0.20.5*0.3CpGIsland:Viterbi算法(2)vATCGCAβ1C+G+A+0.1T+0.015C-G-A-0.15T-0.02250.15*0.5*0.20.15*0.5*0.30.1*0.5*0.30.1*0.5*0.2CpGIsland:Viterbi算法(3)vATCGCAβ1C+0.0034G+A+0.1T+0.015C-0.00225G-A-0.15T-0.02250.0225*0.5*0.20.015*0.5*0.30.0225*0.5*0.30.015*0.5*0.2CpGIsland:Viterbi算法(4)vATCGCAβ1C+0.0034G+0.0005A+0.1T+0.015C-0.00225G-0.00034A-0.15T-0.02250.00225*0.5*0.20.00225*0.5*0.30.0034*0.5*0.20.0034*0.5*0.3CpGIsland:Viterbi算法(4)vATCGCAβ1C+0.00340.000075G+0.0005A+0.1T+0.015C-0.002250.00005G-0.00034A-0.15T-0.02250.0005*0.5*0.30.0005*0.5*0.20.0034*0.5*0.30.0034*0.5*0.2CpGIsland:Viterbi算法(5)vATCGCAβ1C+0.00340.000075G+0.0005A+0.10.0000075T+0.015C-0.002250.00005G-0.00034A-0.150.0000112T-0.02250.000075*0.5*0.20.000075*0.5*0.30.00005*0.5*0.30.00005*0.5*0.3CpGIsland:Viterbi算法(6)vATCGCAβ1C+0.00340.000075G+0.0005A+0.10.0000075T+0.015C-0.002250.00005G-0.00034A-0.150.0000112T-0.0225CpGIsland:预测结果1.ATCGCA:其中，CGC被预测为CpGIslandATCGCA2.Viterbi算法：求出在当前结果最大的概率值，以及保存相应的路线3.递归算法：动态规划的算法4.该例中，我们假设状态转移概率矩阵已知5.如何推算状态的概率矩阵？参数估计：

Baum-Welch(EM)算法目的：给定观察值序列O，通过计算确定一个模型H，使得P(O|H)最大算法步骤： 1.初始模型（待训练模型）H0, 2.基于H0

以及观察值序列O，训练新模型

H； 3.如果log

P(O