序列模式識別課件

生物信息學：預測1.生物信息學最核心的問題：預測2.生物信息學工具的作用：預測3.生物信息學所有的分析：預測4.基本假設(貝葉斯的哲學理念)：我們能夠通過對已知世界的觀察，總結經驗，并以此來預測未知世界已經存在或者即將發(fā)生的事物/事件5.在生物信息學中的應用：對現有的數據，使用合適的算法，進行訓練，構建計算模型和計算工具，預測未知的現象本章內容提要1.統(tǒng)計學基礎2.序列模式3.預測性能檢驗4.位點特異性打分矩陣(PSSM)5.模體發(fā)現：GibbsSampler等6.馬爾科夫及隱馬爾科夫模型7.模式識別的其他算法簡介1.統(tǒng)計學基礎排列組合從N個物品中取出k個物品的排列數（排序）：從N個物品中取出k個物品的組合數（不排序）：概率模型概率模型:一個能夠通過不同的概率產生不同結果的模型。概率模型可以模擬或者仿真某一類型的所有事件，并且對每個事件賦予一個概率。色子模型：一個色子存在6個概率值：p1,p2,…,p6，其中，擲出i的概率為pi(i=1,2,…,6)。因此：pi≥0，且考慮三次連續(xù)的擲色子，結果為[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)稱為概率密度函數5.因此：且二項分布1.事件只有兩種可能出現的結果。例如擲硬幣，正面記為“1”，反面記為“0”。2.則，擲硬幣N次，有k次是1的概率為：二項分布的期望值期望值代表了隨機變量的“平均”值。它是把每個可能取值乘以對應的概率，然后累加起來。期望值E(x)=μ二項分布的與標準方差標準方差描述了隨機變量中具有正概率值的分散性。所有可能的值離期望值的距離的平方，再乘以對應的概率。方差VarX=σ2泊松分布1.稀有事件發(fā)生的概率：在一個連續(xù)的時間或空間中，稀有離散變量出現的概率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發(fā)現一條1mbp的序列上存在5個某種調控信號，該調控信號在人的基因組上平均每500kbp出現一個。那么，完全是隨機產生該種情況的概率是多少？本例中，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超幾何分布雙尾概率方法一：所有出現概率<=觀察表概率的概率之和方法二：雙尾概率=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個學生，比較是否男生比女生更喜歡玩電子游戲。數據統(tǒng)計如下：玩游戲不玩游戲男生4515女生2713P-value=0.496854471943056>0.05統(tǒng)計性不顯著！序列模式識別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.缺點：不能反映計算模型的穩(wěn)定性Leave-one-outvalidation每次從數據集中去掉一個重新進行訓練，構建預測模型，并對去除的數據進行預測。保證每個數據去掉一次n-foldcross-validation將數據集分成n組，并保證陽性數據與陰性數據的比例與原數據相同隨意將n-1組作為訓練數據，1組作為檢驗數據，計算性能重復若干次，例如，重復20次計算平均值缺點：每次計算結果有偏差預測性能及穩(wěn)定性1.Self-consistency:反映檢驗性能（對已知數據的預測能力）2.Leave-one-outvalidation&n-foldcross-validation:反映預測系統(tǒng)的穩(wěn)定性（對未知數據的預測能力）3.預測性能vs.檢驗性能a.差距較?。合到y(tǒng)穩(wěn)定b.差距過大：系統(tǒng)不穩(wěn)定，數據過訓練過訓練1.根據已知數據構建的模型只能很好的適用于訓練數據2.不適合用來預測未知數據3.對訓練數據的微小改變對于預測性能影響過大4.預測工具過訓練：只能很好的符合訓練數據，而對新數據則性能很差4.位點特異性打分矩陣(1)PositionSpecificScoringMatrix(PSSM)/WeightMatrixModel(WMM)(2)對蛋白質家族進行多序列比對分析，發(fā)現結果中保守的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.模體發(fā)現：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不再發(fā)生變化StrongMotifACGTAGCAGibbsSampler:總結1.模體發(fā)現的一種隨機算法(MonteCarlo)2.尋找次優(yōu)解的算法3.根據PSSM/WMM對隨機抽取的序列進行打分來調整采樣，直到結果收斂4.不能夠保證每次運算的結果一致：需要多運算幾次，并進行比較5.對蛋白質、DNA、RNA序列模體的發(fā)現有幫助期望最大化算法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.背景概率:根據剩下的序列計算四種堿基的概率分布似然性概率值的計算似然性概率值的計算(2)計算每條序列，在不同的起始位置，其似然性的概率值Estep:起始位置估計Z值：motif在不同位置起始的幾率值假設，motif在任意位置起始的概率相同，則Z值最大化，即為“最可能的起始位置”Mstep：P值最大化根據選擇的最大Z值，重新計算矩陣，并計算P值最大的motif；P值最大原先的motifEM算法：迭代Gibbs&EM:總結1.基本假設：所有序列都擁有，且僅擁有一個motif2.估算兩個關聯的函數：Gibbs(WMM&似然性)，EM(motif起始位置，Z值&似然性)3.利用兩個函數的其中之一修正另一個，采取迭代/反復計算的方法，使結果收斂4.不保證得到的結果為最優(yōu)，近似算法有待解決的問題1.給定的一組序列，可能的motif僅在部分序列中出現，怎么解決？2.給定一組序列，其中存在某種motif可能在序列上出現兩次以上，如何解決？6.馬爾科夫及隱馬爾科夫模型1870年，俄國有機化學家VladimirV.Markovnikov首次提出馬爾科夫模型馬爾科夫模型馬爾科夫鏈隱馬爾科夫模型VladimirV.Markovnikov馬爾科夫模型馬爾科夫模型：隨機過程的一種，主要特點為“無后效性”，即根據當前的狀態(tài)即可完全確定將來的狀態(tài)馬爾科夫性&馬爾科夫鏈1.定義：對于隨機變量X1,X2,X3…,這些變量的范圍，即他們所有可能取值的集合，被稱為“狀態(tài)空間”，而Xn的值則是在時間n的狀態(tài)。如果Xn+1對于過去狀態(tài)的條件概率分布僅是Xn的一個函數，則符合馬爾科夫性：2.具有馬爾科夫性的過程稱為馬爾科夫過程3.時間(先后順序)和狀態(tài)都離散的馬爾科夫過程稱為馬爾科夫鏈馬爾科夫模型：參數估計轉移概率：K-order馬爾科夫模型一階馬爾科夫模型：當前位置僅依賴前一位k階馬爾科夫模型：當前位置依賴前一位，而前一位依賴前兩位,…,前k-1位依賴前k位0階馬爾科夫模型：位點獨立Markov&PSSM1.對真實的數據進行訓練，PSSM=~0階馬爾科夫模型2.對新序列的掃描：從頭至尾，每次移動1~n位（窗口滑動的方法）3.分別計算窗口內的序列，是(+)和(-)的概率，計算log-oddsratio4.設定閾值，若高于閾值，則預測為陽性另外長度不確定！起始位置不知！Markovmodels&PSSM:Notwork!!!隱馬爾科夫模型(HMM)隱馬爾科夫模型：

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

2.狀態(tài)之間的轉移概率未知與馬爾可夫模型的本質區(qū)別： 隱馬模型觀察到的符號并不是與狀態(tài)一一對應，而是通過一組概率分布相聯系ProfileHMM1.多序列比對的結果中，氨基酸之間存在的關系有匹配(M),插入(I)和缺失(D):三種狀態(tài)2.HMM：三種狀態(tài)之間的轉換關系未知->hidden->轉移概率3.每個位置上的氨基酸/堿基以及插入、缺失的頻率/概率可以通過觀測求得->nothidden4.模型訓練：通過訓練，估算轉移概率例：CpG島的HMM1.CpG島：在人的基因組中，如果雙堿基對CG出現，則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:對當前未知的堿基，跳轉到下一個位置，究竟是I還是G的概率，未知2.Observable:I和G中的四種堿基分布的概率能夠通過實際數據的觀測進行計算轉移概率發(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:預測結果1.ATCGCA:其中，CGC被預測為CpGIslandATCGCA2.Viterbi算法：求出在當前結果最大的概率值，以及保存相應的路線3.遞歸算法：動態(tài)規(guī)劃的算法4.該例中，我們假設狀態(tài)轉移概率矩陣已知5.如何推算狀態(tài)的概率矩陣？參數估計：

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

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

H； 3.如果log

P(O