序列模式識別課件

生物信息學：預測1.生物信息學最核心的問題：預測2.生物信息學工具的作用：預測3.生物信息學所有的分析：預測4.基本假設(貝葉斯的哲學理念)：我們能夠通過對已知世界的觀察，總結(jié)經(jīng)驗，并以此來預測未知世界已經(jīng)存在或者即將發(fā)生的事物/事件5.在生物信息學中的應用：對現(xiàn)有的數(shù)據(jù)，使用合適的算法，進行訓練，構(gòu)建計算模型和計算工具，預測未知的現(xiàn)象本章內(nèi)容提要1.統(tǒng)計學基礎2.序列模式3.預測性能檢驗4.位點特異性打分矩陣(PSSM)5.模體發(fā)現(xiàn)：GibbsSampler等6.馬爾科夫及隱馬爾科夫模型7.模式識別的其他算法簡介1.統(tǒng)計學基礎排列組合從N個物品中取出k個物品的排列數(shù)（排序）：從N個物品中取出k個物品的組合數(shù)（不排序）：概率模型概率模型:一個能夠通過不同的概率產(chǎn)生不同結(jié)果的模型。概率模型可以模擬或者仿真某一類型的所有事件，并且對每個事件賦予一個概率。色子模型：一個色子存在6個概率值：p1,p2,…,p6，其中，擲出i的概率為pi(i=1,2,…,6)。因此：pi≥0，且考慮三次連續(xù)的擲色子，結(jié)果為[1，6，3]，則總概率為：p1p6p3概率分布1.考慮連續(xù)變量x，例如：物體的重量。則當重量確切為1公斤時的概率，為0。2.變量的區(qū)間：P(x0≤x≤x1)3.當區(qū)間無限小->0時，上式：P(x-δx/2

≤x≤x+δx/2

)=f(x)δx4.f(x)稱為概率密度函數(shù)5.因此：且二項分布1.事件只有兩種可能出現(xiàn)的結(jié)果。例如擲硬幣，正面記為“1”，反面記為“0”。2.則，擲硬幣N次，有k次是1的概率為：二項分布的期望值期望值代表了隨機變量的“平均”值。它是把每個可能取值乘以對應的概率，然后累加起來。期望值E(x)=μ二項分布的與標準方差標準方差描述了隨機變量中具有正概率值的分散性。所有可能的值離期望值的距離的平方，再乘以對應的概率。方差VarX=σ2泊松分布1.稀有事件發(fā)生的概率：在一個連續(xù)的時間或空間中，稀有離散變量出現(xiàn)的概率2.N->∞,E(x)=μe=2.71828…泊松分布與二項分布的近似對于大的N及小的p值的二項分布，能夠相當準確地用一個參數(shù)為μ=Np的泊松分布近似。當實驗次數(shù)很多而概率很小時：二項分布~泊松分布例1：鳥槍法的覆蓋率假設：需要測序的BAC長度200kbp;總共測序的序列數(shù)量：N;每次測序：500bp；每次測序的覆蓋率p：500/200kbp=0.0025因此：總覆蓋率μ=Np(每個點平均覆蓋到的次數(shù))k:測序能夠覆蓋到點X的次數(shù)。鳥槍法：覆蓋率點X被覆蓋k次的概率：（二項分布~泊松分布)當點X一次都不被覆蓋時，k=0;此時的概率為：覆蓋率vs.準確性泊松分布：例2Prof.Gene發(fā)現(xiàn)一條1mbp的序列上存在5個某種調(diào)控信號，該調(diào)控信號在人的基因組上平均每500kbp出現(xiàn)一個。那么，完全是隨機產(chǎn)生該種情況的概率是多少？本例中，N=3.0*109bp->∞,E(x)=μ=2(1mbp)統(tǒng)計性顯著：p-value<0.05非常顯著：p-value<0.01超幾何分布與二項分布的區(qū)別：不放回抽樣。例：有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超幾何分布雙尾概率方法一：所有出現(xiàn)概率<=觀察表概率的概率之和方法二：雙尾概率=2×min(左尾概率,右尾概率)超幾何分布：例Prof.Gene從26873個人的蛋白質(zhì)中預測了2264個具有某種特定功能的底物，并進行進一步的分析。其中，已知有421個人的蛋白質(zhì)具有某種功能結(jié)構(gòu)域D，而在預測的2264個底物中，有94個蛋白質(zhì)具有結(jié)構(gòu)域D。問：結(jié)構(gòu)域D在2264個底物中是顯著出現(xiàn)，顯著不出現(xiàn)，還是隨機出現(xiàn)？問題轉(zhuǎn)化：在26873個人的蛋白質(zhì)中有421個具有功能結(jié)構(gòu)域D，任意取出2264個蛋白質(zhì)，其中至少有94個具有功能結(jié)構(gòu)域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：再例假設，我們調(diào)查了100個學生，比較是否男生比女生更喜歡玩電子游戲。數(shù)據(jù)統(tǒng)計如下：玩游戲不玩游戲男生4515女生2713P-value=0.496854471943056>0.05統(tǒng)計性不顯著！序列模式識別2.序列模式1.功能結(jié)構(gòu)域，functionaldomain2.模塊，BLOCK3.模體，motif4.模式，pattern/profile功能結(jié)構(gòu)域1.具有完整的、獨立的三級結(jié)構(gòu)2.具有特定的生物學功能3.一般長度，幾十到幾百個氨基酸4.允許插入/缺失，即允許存在gap模塊/BLOCK1.幾個到幾十個氨基酸2.無gap，從全局多序列比對的結(jié)果直接處理得到3.描述蛋白質(zhì)家族或者一類蛋白質(zhì)的序列保守性BLOCK模體/Motif1.不具有獨立的三級結(jié)構(gòu)2.具有特定的生物學功能：結(jié)合，修飾，細胞亞定位，維持結(jié)構(gòu)，等3.長度一般幾個到幾十個氨基酸或者堿基；4.例如，SUMO化的序列模體：Ψ-K-X-E(Ψ:A,I,L,V,M,F,P;X:任意氨基酸)模式/Pattern/Profile1.在算法上用來描述一類功能結(jié)構(gòu)域、模體或者模塊的表示方式2.根據(jù)序列數(shù)據(jù)，構(gòu)建的預測模型3.數(shù)據(jù)形式：正則表達式4.用來預測新的可能符合特定模式的序列5.例如，直接將Ψ-K-X-E視為SUMO化位點的，普適的“模式”，則可以預測所有包含該模式的蛋白質(zhì)序列3.預測性能的計算和檢驗1.樣本/檢驗數(shù)據(jù)：陽性數(shù)據(jù)(P)，陰性數(shù)據(jù)(N)a.陽性數(shù)據(jù)(P)：真實的，被實驗所證實的數(shù)據(jù)b.陰性數(shù)據(jù)(N)：被實驗所證明為無功能的數(shù)據(jù)2.對于預測結(jié)果的評測，定義：a.真陽性(TP):陽性數(shù)據(jù)中被預測為陽性的數(shù)據(jù)b.假陽性(FP):陰性數(shù)據(jù)中被預測為陽性的數(shù)據(jù)c.真陰性(TN):陰性數(shù)據(jù)中被預測為陰性的數(shù)據(jù)d.假陰性(FN):陽性數(shù)據(jù)中被預測為陰性的數(shù)據(jù)TPFPFNTNPositiveNegativeCutoff常用的檢驗指標1.靈敏度(Sensitivity,Sn)對于真實的數(shù)據(jù)，能夠預測成“真”的比例是多少2.特異性(Specificity,Sp)對于陰性的數(shù)據(jù)，能夠預測成“假”的比例是多少3.準確性(Accuracy,Ac)對于整個數(shù)據(jù)集(包括陽性和陰性數(shù)據(jù))，預測總共的準確比例是多少4.馬修相關系數(shù)(Mathewcorrelationcoefficient,MCC)

當陽性數(shù)據(jù)的數(shù)量與陰性數(shù)據(jù)的數(shù)量差別較大時，能夠更為公平的反映預測能力，值域[-1,1]常用的檢驗指標(2)ROCcurveX軸：1-SpY軸：SnROC的面積越大，表明其預測能力越強預測性能的計算Self-consistencyLeave-one-outvalidationn-foldcross-validationSelf-consistency1.將訓練數(shù)據(jù)當成測試數(shù)據(jù)訓練數(shù)據(jù)中所有的陽性數(shù)據(jù)為測試數(shù)據(jù)中的陽性數(shù)據(jù)訓練數(shù)據(jù)中所有的陰性數(shù)據(jù)為測試數(shù)據(jù)中的陰性數(shù)據(jù)2.反映當前預測工具對目前已知的數(shù)據(jù)的預測能力3.假設：根據(jù)目前已知的數(shù)據(jù)所構(gòu)建的計算模型能夠反映未知的數(shù)據(jù)的模式4.缺點：不能反映計算模型的穩(wěn)定性Leave-one-outvalidation每次從數(shù)據(jù)集中去掉一個重新進行訓練，構(gòu)建預測模型，并對去除的數(shù)據(jù)進行預測。保證每個數(shù)據(jù)去掉一次n-foldcross-validation將數(shù)據(jù)集分成n組，并保證陽性數(shù)據(jù)與陰性數(shù)據(jù)的比例與原數(shù)據(jù)相同隨意將n-1組作為訓練數(shù)據(jù)，1組作為檢驗數(shù)據(jù)，計算性能重復若干次，例如，重復20次計算平均值缺點：每次計算結(jié)果有偏差預測性能及穩(wěn)定性1.Self-consistency:反映檢驗性能（對已知數(shù)據(jù)的預測能力）2.Leave-one-outvalidation&n-foldcross-validation:反映預測系統(tǒng)的穩(wěn)定性（對未知數(shù)據(jù)的預測能力）3.預測性能vs.檢驗性能a.差距較小：系統(tǒng)穩(wěn)定b.差距過大：系統(tǒng)不穩(wěn)定，數(shù)據(jù)過訓練過訓練1.根據(jù)已知數(shù)據(jù)構(gòu)建的模型只能很好的適用于訓練數(shù)據(jù)2.不適合用來預測未知數(shù)據(jù)3.對訓練數(shù)據(jù)的微小改變對于預測性能影響過大4.預測工具過訓練：只能很好的符合訓練數(shù)據(jù)，而對新數(shù)據(jù)則性能很差4.位點特異性打分矩陣(1)PositionSpecificScoringMatrix(PSSM)/WeightMatrixModel(WMM)(2)對蛋白質(zhì)家族進行多序列比對分析，發(fā)現(xiàn)結(jié)果中保守的BLOCK(3)根據(jù)BLOCK序列推導相應的PSSM(4)不考慮gap的影響(5)BLOCK長度一般在幾個~幾十個殘基/堿基BLOCK->PSSM代表每一列二十種氨基酸矩陣中的數(shù)值：當前位置上，某種氨基酸出現(xiàn)的頻率的log值第二種PSSM每一個位置上顯示每種氨基酸或者堿基出現(xiàn)的頻率堿基的位置四種堿基第三種PSSM每一個位置顯示氨基酸/堿基出現(xiàn)的概率PSSM矩陣使用P(S|+)，根據(jù)陽性訓練數(shù)據(jù)計算出來的概率；未知序列：ACGGTACGG背景概率選擇,P(S|-)1.負樣本/陰性數(shù)據(jù)的概率計算2.計算方法：A.DNA序列，四種堿基出現(xiàn)的頻率B.蛋白質(zhì)序列，20種氨基酸出現(xiàn)的頻率OddsRatioLog-oddsRatio計算流程：滑動窗口窗口寬度9bp，依次打分設定閾值（Threshold），凡是高于閾值的預測為陽性，低于閾值的預測為陰性5.模體發(fā)現(xiàn)：GibbsSamplerGibbsSampler是一種Monte-Carlo類的方法

隨機抽樣對于輸入序列，找到一個最大的似然函數(shù)GibbsSampler算法(1)1.從每條序列上隨機的抽取一段序列，序列長度固定所有序列motifGibbsSampler算法(2)2.構(gòu)建PSSM/權(quán)重矩陣GibbsSampler算法(3)3.隨機挑選一條序列GibbsSampler算法(4)4.用構(gòu)建好的PSSM對該序列上所有可能的motif進行打分(窗口滑動，每次1個氨基酸或者堿基)GibbsSampler算法(5)5.根據(jù)似然性的計算，得到似然值最大的模體，即新的motifGibbsSampler算法(6)6.更新PSSM矩陣GibbsSampler算法(7)7.反復迭代計算，直到似然性結(jié)果與PSSM不再發(fā)生變化StrongMotifACGTAGCAGibbsSampler:總結(jié)1.模體發(fā)現(xiàn)的一種隨機算法(MonteCarlo)2.尋找次優(yōu)解的算法3.根據(jù)PSSM/WMM對隨機抽取的序列進行打分來調(diào)整采樣，直到結(jié)果收斂4.不能夠保證每次運算的結(jié)果一致：需要多運算幾次，并進行比較5.對蛋白質(zhì)、DNA、RNA序列模體的發(fā)現(xiàn)有幫助期望最大化算法1.ExpectationMaximizationAlgorithm2.已開發(fā)工具：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.背景概率:根據(jù)剩下的序列計算四種堿基的概率分布似然性概率值的計算似然性概率值的計算(2)計算每條序列，在不同的起始位置，其似然性的概率值Estep:起始位置估計Z值：motif在不同位置起始的幾率值假設，motif在任意位置起始的概率相同，則Z值最大化，即為“最可能的起始位置”Mstep：P值最大化根據(jù)選擇的最大Z值，重新計算矩陣，并計算P值最大的motif；P值最大原先的motifEM算法：迭代Gibbs&EM:總結(jié)1.基本假設：所有序列都擁有，且僅擁有一個motif2.估算兩個關聯(lián)的函數(shù)：Gibbs(WMM&似然性)，EM(motif起始位置，Z值&似然性)3.利用兩個函數(shù)的其中之一修正另一個，采取迭代/反復計算的方法，使結(jié)果收斂4.不保證得到的結(jié)果為最優(yōu)，近似算法有待解決的問題1.給定的一組序列，可能的motif僅在部分序列中出現(xiàn)，怎么解決？2.給定一組序列，其中存在某種motif可能在序列上出現(xiàn)兩次以上，如何解決？6.馬爾科夫及隱馬爾科夫模型1870年，俄國有機化學家VladimirV.Markovnikov首次提出馬爾科夫模型馬爾科夫模型馬爾科夫鏈隱馬爾科夫模型VladimirV.Markovnikov馬爾科夫模型馬爾科夫模型：隨機過程的一種，主要特點為“無后效性”，即根據(jù)當前的狀態(tài)即可完全確定將來的狀態(tài)馬爾科夫性&馬爾科夫鏈1.定義：對于隨機變量X1,X2,X3…,這些變量的范圍，即他們所有可能取值的集合，被稱為“狀態(tài)空間”，而Xn的值則是在時間n的狀態(tài)。如果Xn+1對于過去狀態(tài)的條件概率分布僅是Xn的一個函數(shù)，則符合馬爾科夫性：2.具有馬爾科夫性的過程稱為馬爾科夫過程3.時間(先后順序)和狀態(tài)都離散的馬爾科夫過程稱為馬爾科夫鏈馬爾科夫模型：參數(shù)估計轉(zhuǎn)移概率：K-order馬爾科夫模型一階馬爾科夫模型：當前位置僅依賴前一位k階馬爾科夫模型：當前位置依賴前一位，而前一位依賴前兩位,…,前k-1位依賴前k位0階馬爾科夫模型：位點獨立Markov&PSSM1.對真實的數(shù)據(jù)進行訓練，PSSM=~0階馬爾科夫模型2.對新序列的掃描：從頭至尾，每次移動1~n位（窗口滑動的方法）3.分別計算窗口內(nèi)的序列，是(+)和(-)的概率，計算log-oddsratio4.設定閾值，若高于閾值，則預測為陽性另外長度不確定！起始位置不知！Markovmodels&PSSM:Notwork!!!隱馬爾科夫模型(HMM)隱馬爾科夫模型：

1.表示狀態(tài)的可觀察符號出現(xiàn)概率已知

2.狀態(tài)之間的轉(zhuǎn)移概率未知與馬爾可夫模型的本質(zhì)區(qū)別： 隱馬模型觀察到的符號并不是與狀態(tài)一一對應，而是通過一組概率分布相聯(lián)系ProfileHMM1.多序列比對的結(jié)果中，氨基酸之間存在的關系有匹配(M),插入(I)和缺失(D):三種狀態(tài)2.HMM：三種狀態(tài)之間的轉(zhuǎn)換關系未知->hidden->轉(zhuǎn)移概率3.每個位置上的氨基酸/堿基以及插入、缺失的頻率/概率可以通過觀測求得->nothidden4.模型訓練：通過訓練，估算轉(zhuǎn)移概率例：CpG島的HMM1.CpG島：在人的基因組中，如果雙堿基對CG出現(xiàn)，則C通常被甲基化。并且，甲基化的C很快會突變成T。因此基因組中CpG島非常少。然而，在基因的起始位置，例如promotor區(qū)域，因為功能的保守性，其序列很少突變，CpG的含量能夠保持在40~60%2.Howtopredict?PSSM&Markovarenotworkatall!CpG島：HMM存在兩種狀態(tài)：是CpG島(CpGIsland,I)，不是CpG島(Genome,G)CpG島：HMM1.Hidden:對當前未知的堿基，跳轉(zhuǎn)到下一個位置，究竟是I還是G的概率，未知2.Observable:I和G中的四種堿基分布的概率能夠通過實際數(shù)據(jù)的觀測進行計算轉(zhuǎn)移概率發(fā)散概率預測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:預測結(jié)果1.ATCGCA:其中，CGC被預測為CpGIslandATCGCA2.Viterbi算法：求出在當前結(jié)果最大的概率值，以及保存相應的路線3.遞歸算法：動態(tài)規(guī)劃的算法4.該例中，我們假設狀態(tài)轉(zhuǎn)移概率矩陣已知5.如何推算狀態(tài)的概率矩陣？參數(shù)估計：

Baum-Welch(EM)算法目的：給定觀察值序列O，通過計算確定一個模型H，使得P(O|H)最大算法步驟： 1.初始模型（待訓練模型）H0, 2.基于H0

以及觀察值序列O，訓練新模型

H； 3.如果log

P(O