以二階隱藏式馬可夫模型預測特定蛋白激酶磷酸化的位置
122
0
0
全文
(2) 以二階隱藏式馬可夫模型預測特定蛋白激酶磷酸化的位置 Predicting Protein Kinase-Specific Phosphorylation Sites Using Second-Order Hidden Markov Models. 研 究 生：葉智國. Student：Chih-Kuo Yeh. 指導教授：何信瑩. Advisor：Shinn-Ying Ho. 國 立 交 通 大 學 生 物 資 訊 研 究 所 碩 士 論 文 A Thesis Submitted to Institute of Bioinformatics Department of Biological Science National Chiao Tung University in partial Fulfillment of the Requirements for the Degree of Master in Bioinformatics July 2007 Hsinchu, Taiwan, Republic of China. 中華民國九十六年七月. II.
(3) 以 二 階 隱 藏 式 馬 可 夫 模 型 預 測 特 定 蛋 白 激 酶 磷 酸 化 的 位 置. 學生：葉智國. 指導教授：何信瑩. 國立交通大學生物資訊研究所碩士班. 摘要 蛋白質磷酸化是在蛋白質轉譯後修飾中很重要的機制，在調控基本的細胞進行過程 像是新陳代謝、訊號傳遞、細胞分化和細胞膜穿透性等扮演重要角色。在過去，要標記 已被磷酸化的蛋白質和被磷酸化的位置，即使透過如二維電泳分析和質譜儀分析等新技 術的幫助，仍然要耗費大量的人力與資源。因此發展使用蛋白質序列資訊的電腦輔助預 測軟體來預測磷酸化的位置與它們特定的激酶，可以提供一個關鍵的選擇步驟，用來減 少實驗中候選者的數目。HMMer 是一個用來做蛋白激酶特定磷酸化位置預測的很好軟 體，它使用一階隱藏式馬可夫模型(HMM-1)及 Plan7 的架構，在重要激酶 PKA、PKC、 CDK 等現有資料集上，可分別達到 82%、74% 和 82% 辨識率。 本論文經由對 HMMer 的分析，希望提出一套以二階隱藏式馬可夫模型(HMM-2)為 基底的改良演算法 iHMM (improved HMM)用來預測特定蛋白激酶磷酸化的位置，希望 從序列中取得更豐富的前後文資訊來提升預測正確率。藉由搭配使用貝氏資訊準則 (Bayesian Information Criterion)的模型參數選擇方法，使 HMM-2 在資料集不夠大時能 盡量避免眾所周知的過度適化(over fitting)問題。本論文將 Phospho.ELM 資料庫與 Swiss-Prot 資料庫結合，將有蛋白激酶註解的資料依照蛋白激酶屬性 PKA、PKC、CDK 等分類建立十八個資料集，然後分別用 iHMM 建立預測模型，以 5-fold 交互驗證做 30 次獨立測試的結果，並經由與 HMMer 和傳統 HMM-1 的效能比較來評估 iHMM。實驗 結果發現 iHMM 的辨識率與 HMMer 相比得到將近平均 4.3%的提升，而跟傳統 HMM-1 相比則得到將近平均 3.6%的提升。本論文並進一步探討比較 HMMer、傳統 HMM-1 和 iHMM 三種方法的特性與優缺點。. 關鍵字：磷酸化、激酶、隱藏式馬可夫模型、貝氏資訊準則. I.
(4) Predicting Protein Kinase-Specific Phosphorylation Sites Using Second-Order Hidden Markov Models. Student: Chih-Kuo Yeh. Advisor: Shinn-Ying Ho. Institute of Bioinformatics National Chiao Tung University. ABSTRACT Protein phosphorylation is an important mechanism of posttranslational modifications and it plays important roles in regulation of essential cellular processes such as metabolism, cell signaling, differentiation and membrane transportation. In the past, laboratory identification of phosphorylated proteins and phosphorylation sites is usually tedious and cumbersome. Recently, large-scale methods of two-dimensional gel analysis and mass spectrometry techniques were applied to efficiently detect phosphorylation sites. However, experimental identification of phosphorylation sites is still expensive. Therefore, computational prediction of phosphorylation sites with their specific kinases using protein's primary sequences can provide a crucial selection step to reduce the number of candidates. HMMer using first-order Hidden Markov Model (HMM-1) with the Plan7 architecture is a conventional tool for prediction of kinase-specific phosphorylation sites and the prediction accuracies of HMMer are 82%, 74% and 82% for existing data sets of the important kinases PKA, PKC and CDK, respectively. From the analysis of HMMer, this thesis aims to propose an improved algorithm iHMM using the second-order HMM (HMM-2) with more context information of sequences to advance prediction accuracy. With the use of Bayesian Information Criterion on the selection of model parameters, iHMM tries to avoid the known over-fitting problem when the sizes of data sets are not large. This thesis established 18 data sets of annotated kinases family such as PKA, PKC, CDK, etc from the Phospho.ELM database and Swiss-Prot database. The performance of iHMM is compared with those of HMMer and the conventional HMM-1 using 5-fold cross validation for 30 independent runs. Simulation results reveal that iHMM can improve the average accuracies of HMMer and HMM-1 near to 4.3% and 3.6%, respectively. Furthermore, this thesis investigated the advantages and disadvantages of iHMM, compared with those of HMMer and HMM-1.. Keywords：phosphorylation, kinase, Hidden Markov Model, Bayesian Information Criterion II.
(5) 致謝 首先，我要特別感謝我的指導教授－何信瑩老師，讓我有機會進入這個生物資訊的 領域，並且在研究陷入瓶頸時，不厭其煩很有耐心的指導我，在何老師身上學到了除了 研究論文、學習到許多生物資訊相關技術以外，還有更重要的是做人處事的態度。 此篇論文要感謝義雄，冠維與俊維學長在論文寫作上幫我檢查並提出修正與改進的 地方。還有特别感謝俊維，家達與凱迪在生物方面提供專業知識與寶貴意見。另外感謝 孝邦與鍾錡讓我有豐富的研究所生活，同時也要感謝廖芹的祝福，使我可以順利畢業。 感謝實驗室所有學長、學弟們，在每次實驗室開會時的批評與指教，以及精神上的鼓勵。 這短短的兩年研究生活，不管在學習、運動、玩樂大家總是互相照應，是我這一生難忘 的回憶。另外也要感謝黃憲達老師，以及李宗夷學長提供生物資料的協助，以讓我可以 完成這次研究。 最後我要感謝我的父母的栽培與無時無刻的關心，你們是我最好的依靠，有你們的 支持與鼓勵，讓我沒有後顧之憂，可以專心於學業上的研究！. III.
(6) 目錄 摘要 ............................................................................................................................................1 ABSTRACT .............................................................................................................................II 致謝 ......................................................................................................................................... III 目錄 ......................................................................................................................................... IV 圖目錄 ..................................................................................................................................... VI 表目錄 ....................................................................................................................................... X 符號說明 ................................................................................................................................. XI 第一章 緒論 ..............................................................................................................................1 1.1 研究動機 ..........................................................................................................................1 1.2 問題描述與研究方向 ......................................................................................................1 1.3 章節概要 ..........................................................................................................................3 第二章 蛋白激酶基本介紹 ......................................................................................................4 2.1 蛋白激酶的分類 ..............................................................................................................4 2.1 蛋白激酶的功能 ..............................................................................................................4 2.2 蛋白激酶的結構 ..............................................................................................................5 2.3 特定蛋白激酶磷酸化的位置 ..........................................................................................6 第三章 方法 ..............................................................................................................................8 3.1 First-order HMM ..............................................................................................................8 3.2 Profile HMM ....................................................................................................................8 3.3 Second-order HMM..........................................................................................................9 3.4 分數計算 ........................................................................................................................10 3.5 訓練演算法 .................................................................................................................... 11 第四章 提出iHMM方法論 ....................................................................................................16 4.1 想法起源 ........................................................................................................................16 4.2 iHMM 之設計 ...............................................................................................................16 第五章 實驗資料與結果討論 ................................................................................................22 5.1 實驗資料 ........................................................................................................................22 5.2 實驗結果 ........................................................................................................................22 5.3 相關系數分析 ................................................................................................................28 5.4 訓練兩類的結果 ............................................................................................................30 第六章 結論與展望 ................................................................................................................32 6.1 結論 ................................................................................................................................32 6.2 未來展望 ........................................................................................................................32 參考文獻 ..................................................................................................................................34 IV.
(7) 附錄: 詳細實驗數據...............................................................................................................37. V.
(8) 圖目錄 圖 1.1:系統架構示意圖 .............................................................................................................2 圖 2.1: PKA 三級結構 ..............................................................................................................5 圖 3.1: Plan7 結構 .....................................................................................................................9 圖 3.2:原始HMM-2..................................................................................................................10 圖 3.3: HMM-1 等價模型.......................................................................................................10 圖 4.1:狀態數目與辨識率之間的關係的示意圖 ...................................................................18 圖 4.2:狀態數目與BIC的示意圖 ............................................................................................20 圖 5.1: PKA sequence logos .....................................................................................................24 圖 5.2: PKA (S)資料的HMMer結構 .......................................................................................24 圖 5.3: PKA (S)資料的HMM-1 結構 ......................................................................................24 圖 5.4: PKA (S)資料的HMM-1 結構的符號觀測機率圖 ......................................................24 圖 5.5: PKA (S)資料的HMM-2 結構 ......................................................................................25 圖 5.6: PKA (S)資料的HMM-2 結構的符號觀測機率圖 ......................................................25 圖 5.7: HMM-1 等價展開圖....................................................................................................25 圖 5.8: PKC (S)資料的HMM-2 結構 ......................................................................................26 圖 5.9: PKC (S)資料的HMM-2 結構的符號觀測機率圖 ......................................................26 圖 5.10: CaM-KII (S)資料的HMM-2 結構 .............................................................................27 圖 5.11: CaM-KII (S)資料的HMM-2 結構的符號觀測機率圖 .............................................27 圖 5.12: CDK (S)資料的HMM-2 結構 ...................................................................................27 圖 5.13: CDK (S)資料的HMM-2 結構的符號觀測機率圖 ...................................................27 圖 5.14: CK1 (S)資料的相關係數分析圖...............................................................................28 圖 5.15: PKA (S) negative資料的HMM-2 結構......................................................................30 圖 5.16: PKA (S)兩類資料於訓練資料的ROC圖 ..................................................................31 圖 5.17: PKA (S)兩類資料於測試資料的ROC圖 ..................................................................31 圖A.1: PKA (S)資料的序列圖案 ............................................................................................37 圖A.2: PKA (S)資料於iHMM的門檻值與正確率對應圖 .....................................................38 圖A.3: PKA (S)資料於HMMer的門檻值與正確率對應圖 ...................................................38 圖A.4: PKA (S)於訓練資料的ROC圖 ....................................................................................39 圖A.5: PKA (S)於測試資料的ROC圖 ....................................................................................39 圖A.6: PKA (S)資料的相關係數分析圖 ................................................................................40 圖A.7: PKA (T)資料的序列圖案 ............................................................................................41 圖A.8: PKA (T)資料於iHMM的門檻值與正確率對應圖 .....................................................42 圖A.9: PKA (T)資料於HMMer的門檻值與正確率對應圖 ...................................................42 圖A.10: PKA (T)於訓練資料的ROC圖 ..................................................................................43 圖A.11: PKA (T)於測試資料的ROC圖 ..................................................................................43 圖A.12: PKA (T)資料的相關係數分析圖 ..............................................................................44 圖A.13: PKB (S)資料的序列圖案 ..........................................................................................45 圖A.14: PKB (S)資料於iHMM的門檻值與正確率對應圖 ...................................................46 圖A.15: PKB (S)資料於HMMer的門檻值與正確率對應圖 .................................................46 圖A.16: PKB (S)於訓練資料的ROC圖 ..................................................................................47 圖A.17: PKB (S)於測試資料的ROC圖 ..................................................................................47 圖A.18: PKB (S)資料的相關係數分析圖 ..............................................................................48 VI.
(9) 圖A.19: PKC (S)資料的序列圖案 ..........................................................................................49 圖A.20: PKC (S)資料於iHMM的門檻值與正確率對應圖 ...................................................50 圖A.21: PKC (S)資料於HMMer的門檻值與正確率對應圖 .................................................50 圖A.22: PKC (S)於訓練資料的ROC圖 ..................................................................................51 圖A.23: PKC (S)於測試資料的ROC圖 ..................................................................................51 圖A.24: PKC (S)資料的相關係數分析圖 ..............................................................................52 圖A.25: PKC (T)資料的序列圖案 ..........................................................................................53 圖A.26: PKC (T)資料於iHMM的門檻值與正確率對應圖 ...................................................54 圖A.27: PKC (T)資料於HMMer的門檻值與正確率對應圖 .................................................54 圖A.28: PKC (T)於訓練資料的ROC圖 ..................................................................................55 圖A.29: PKC (T)於測試資料的ROC圖 ..................................................................................55 圖A.30: PKC (T)資料的相關係數分析圖 ..............................................................................56 圖A.31: PKG (S)資料的序列圖案 ..........................................................................................57 圖A.32: PKG (S)資料於iHMM的門檻值與正確率對應圖 ...................................................58 圖A.33: PKG (S)資料於HMMer的門檻值與正確率對應圖 .................................................58 圖A.34: PKG (S)於訓練資料的ROC圖 ..................................................................................59 圖A.35: PKG (S)於測試資料的ROC圖 ..................................................................................59 圖A.36: PKG (S)資料的相關係數分析圖 ..............................................................................60 圖A.37: CDK (S)資料的序列圖案..........................................................................................61 圖A.38: CDK (S)資料於iHMM的門檻值與正確率對應圖...................................................62 圖A.39: CDK (S)資料於HMMer的門檻值與正確率對應圖.................................................62 圖A.40: CDK (S)於訓練資料的ROC圖..................................................................................63 圖A.41: CDK (S)於測試資料的ROC圖..................................................................................63 圖A.42: CDK (S)資料的相關係數分析圖..............................................................................64 圖A.43: CDK (T)資料的序列圖案..........................................................................................65 圖A.44: CDK (T)資料於iHMM的門檻值與正確率對應圖 ..................................................66 圖A.45: CDK (T)資料於HMMer的門檻值與正確率對應圖 ................................................66 圖A.46: CDK (T)於訓練資料的ROC圖 .................................................................................67 圖A.47: CDK (T)於測試資料的ROC圖 .................................................................................67 圖A.48: CDK (T)資料的相關係數分析圖..............................................................................68 圖A.49: CaM-KII (S)資料的序列圖案 ...................................................................................69 圖A.50: CaM-KII (S)於iHMM的門檻值與正確率對應圖 ....................................................70 圖A.51: CaM-KII (S)於HMMer的門檻值與正確率對應圖 ..................................................70 圖A.52: CaM-KII (S)於訓練資料的ROC圖 ...........................................................................71 圖A.53: CaM-KII (S)於測試資料的ROC圖 ...........................................................................71 圖A.54: CaM-KII (S)資料的相關係數分析圖 .......................................................................72 圖A.55: CK1 (S)資料的序列圖案...........................................................................................73 圖A.56: CK1 (S)資料於iHMM的門檻值與正確率對應圖....................................................74 圖A.57: CK1 (S)資料於HMMer的門檻值與正確率對應圖..................................................74 圖A.58: CK1 (S)於訓練資料的ROC圖 ..................................................................................75 圖A.59: CK1 (S)於測試資料的ROC圖 ..................................................................................75 圖A.60: CK1 (S)資料的相關係數分析圖...............................................................................76 圖A.61: CK2 (S)資料的序列圖案...........................................................................................77 VII.
(10) 圖A.62: CK2 (S)資料於iHMM的門檻值與正確率對應圖....................................................78 圖A.63: CK2 (S)資料於HMMer的門檻值與正確率對應圖..................................................78 圖A.64: CK2 (S)於訓練資料的ROC圖 ..................................................................................79 圖A.65: CK2 (S)於測試資料的ROC圖 ..................................................................................79 圖A.66: CK2 (S)資料的相關係數分析圖...............................................................................80 圖A.67: CK2 (S)資料的序列圖案...........................................................................................81 圖A.68: CK2 (T)資料於iHMM的門檻值與正確率對應圖 ...................................................82 圖A.69: CK2 (T)資料於HMMer的門檻值與正確率對應圖 .................................................82 圖A.70: CK2 (T)於訓練資料的ROC圖 ..................................................................................83 圖A.71: CK2 (T)於測試資料的ROC圖 ..................................................................................83 圖A.72: CK2 (T)資料的相關係數分析圖 ..............................................................................84 圖A.73: MAPK (S)資料的序列圖案.......................................................................................85 圖A.74: MAPK (S)資料於iHMM的門檻值與正確率對應圖................................................86 圖A.75: MAPK (S)資料於HMMer的門檻值與正確率對應圖..............................................86 圖A.76: MAPK (S)於訓練資料的ROC圖 ..............................................................................87 圖A.77: MAPK (S)於測試資料的ROC圖 ..............................................................................87 圖A.78: MAPK (S)資料的相關係數分析圖...........................................................................88 圖A.79: MAPK (T)資料的序列圖案 ......................................................................................89 圖A.80: MAPK (T)資料於iHMM的門檻值與正確率對應圖 ...............................................90 圖A.81: MAPK (T)資料於HMMer的門檻值與正確率對應圖 .............................................90 圖A.82: MAPK (T)於訓練資料的ROC圖 ..............................................................................91 圖A.83: MAPK (T)於測試資料的ROC圖 ..............................................................................91 圖A.84: MAPK (T)資料的相關係數分析圖 ..........................................................................92 圖A.85: ATM (S)資料的序列圖案 ..........................................................................................93 圖A.86: ATM (S)資料於iHMM的門檻值與正確率對應圖 ...................................................94 圖A.87: ATM (S)資料於HMMer的門檻值與正確率對應圖 .................................................94 圖A.88: ATM (S)於訓練資料的ROC圖..................................................................................95 圖A.89: ATM (S)於測試資料的ROC圖..................................................................................95 圖A.90: ATM (S)資料的相關係數分析圖 ..............................................................................96 圖A.91: EGFR (Y)資料的序列圖案 .......................................................................................97 圖A.92: EGFR (Y)資料於iHMM的門檻值與正確率對應圖 ................................................99 圖A.93: EGFR (Y)資料於HMMer的門檻值與正確率對應圖 ..............................................99 圖A.94: EGFR (Y)於訓練資料的ROC圖 .............................................................................100 圖A.95: EGFR (Y)於測試資料的ROC圖 .............................................................................100 圖A.96: EGFR (Y)資料的相關係數分析圖 .........................................................................101 圖A.97: INSR (Y)資料的序列圖案 ......................................................................................102 圖A.98: INSR (Y)資料於iHMM的門檻值與正確率對應圖 ...............................................103 圖A.99: INSR (Y)資料於HMMer的門檻值與正確率對應圖 .............................................103 圖A.100: INSR (Y)於訓練資料的ROC圖 ............................................................................104 圖A.101: INSR (Y)於測試資料的ROC圖 ............................................................................104 圖A.102: INSR (Y)資料的相關係數分析圖 ........................................................................105 圖A.103: SRC (Y)資料的序列圖案 ......................................................................................106 圖A.104: SRC (Y)資料於iHMM的門檻值與正確率對應圖 ...............................................107 VIII.
(11) 圖A.105: SRC (Y)資料於HMMer的門檻值與正確率對應圖 .............................................107 圖A.106: SRC (Y)於訓練資料的ROC圖..............................................................................108 圖A.107: SRC (Y)於測試資料的ROC圖..............................................................................108 圖A.108: SRC (Y)資料的相關係數分析圖 ..........................................................................109. IX.
(12) 表目錄 表 2.1:蛋白激酶磷酸化的位置保留性 .....................................................................................7 表 3.1:自然界胺基酸出現的頻率 ........................................................................................... 11 表 5.1:效能比較表 ...................................................................................................................23 表 5.2:相關係數比較表 ...........................................................................................................29 表 5.3: PKA (S)訓練兩類資料的正確率表 ............................................................................30 表A.1: PKA (S)的 30 次 5-CV 於測試資料的效能比較表...................................................37 表A.2: PKA (T)的 30 次 5-CV 於測試資料的效能比較表 ..................................................41 表A.3: PKB (S)的 30 次 5-CV 於測試資料的效能比較表...................................................45 表A.4: PKC (S)的 30 次 5-CV 於測試資料的效能比較表...................................................49 表A.5: PKC (T)的 30 次 5-CV 於測試資料的效能比較表 ..................................................53 表A.6: PKG (S)的 30 次 5-CV 於測試資料的效能比較表 ..................................................57 表A.7: CDK (S)的 30 次 5-CV 於測試資料的效能比較表 ..................................................61 表A.8: CDK (T)的 30 次 5-CV 於測試資料的效能比較表..................................................65 表A.9: CaM-KII (S)的 30 次 5-CV 於測試資料的效能比較表............................................69 表A.10: CK1 (S)的 30 次 5-CV 於測試資料的效能比較表.................................................73 表A.11: CK2 (S)的 30 次 5-CV 於測試資料的效能比較表 .................................................77 表A.12: CK2 (T) 的 30 次 5-CV 於測試資料的效能比較表...............................................81 表A.13: MAPK (S)的 30 次 5-CV 於測試資料的效能比較表.............................................85 表A.15: ATM (S)的 30 次 5-CV 於測試資料的效能比較表 ................................................93 表A.16: EGFR (Y)的 30 次 5-CV 於測試資料的效能比較表..............................................97 表A.17: INSR (Y)的 30 次 5-CV 於測試資料的效能比較表.............................................102 表A.18: SRC (Y)的 30 次 5-CV 於測試資料的效能比較表 ..............................................106. X.
(13) 符號說明 Σ. ：胺基酸字母的集合. |Σ|. ：胺基酸字母的總數. A. ：狀態轉移機率矩陣. B. ：狀態內部的符號觀測機率矩陣. π. ：初始狀態機率向量. λ. ：HMM 所有的參數的集合. aij. ：First-order HMM 的狀態轉移機率. aijk. ：Second-order HMM 的狀態轉移機率. bj(k) ：狀態 j 觀測到符號 k 的機率 S. ：HMM 所有狀態的集合. si. ：HMM 的第 i 個狀態. Q. ：觀測的狀態向量. qi. ：時間 i 時的狀態. Op. ：positive 序列樣本集合. On. ：negative 序列樣本集合. Oip (k ) ：positive 序列的第 k 條的第 i 個字母 (胺基酸) Oin (k ) ：negative 序列的第 k 條的第 i 個字母 (胺基酸). δ1. ：門檻值取訓練資料正確率最高的那一點來當測試資料的門檻值. δ2. ：門檻值取測試資料正確率最高的那一點來當測試資料的門檻值. XI.
(14) 第一章 緒論 1.1 研究動機 蛋白質磷酸化是在蛋白質轉移後修飾中很重要的機制，影響基本的細胞進程像是新 陳代謝、訊號傳遞、細胞骨架重組、細胞運動、細胞凋亡、分化等和細胞膜穿透性等等。 蛋白質常常被各式各樣的蛋白激酶(protein kinase) 磷酸化。在過去，要標記已被磷酸化 的蛋白質和被磷酸化的位置，時常是費時費力的。最近在技術上有很大的近步像是二維 電泳分析和質譜儀分析，然而，這些實驗依舊是很昂貴的。在這點上，發展僅使用蛋白 質一級序列的資訊來預測磷酸化的位置與它們特定的激酶的電腦輔助預測工具，可以提 供一個關鍵的第一步選擇，來減少實驗中候選的數目，並且可以取代領域專家用人工知 識來篩選。 運用機器學習(machine learning)的技術來預測磷酸化位置已成為生物資訊中熱門 的研究領域。而在這個問題中一些有名的預測系統，像是 Nikolaj Blom 等人發表的以類 神經網路方法為主的預測系統 NetPhos [1]，以及最近 Rune Linding 等人發表的 NetworKIN 系統[2]，然而這些系統在 sensitivity 與 specificity 相差都很大，都是其中一 個很高，另一個很低，而在眾多方法中以 2005 年黃博士等人發表的 KinasePhos 系統 [3]，所使用的 HMMer[4] 在 sensitivity 和 specificity 均有不錯的表現，HMMer 是一個 用來做蛋白激酶特定磷酸化位置的預測不錯的工具，但 HMMer 使用上有些限制跟缺 點，第一，它有結構限制，第二，它在訓練前需要先將資料做過序列對齊，第三，它基 本上無法訓練多類的資料，第四，它基本上是屬於一階隱藏式馬可夫模型(first-order Hidden Markov Model, HMM-1)，基於這些動機，本研究開發一個改良式的統計分析模 型 iHMM (improved HMM)來改良上例缺點。. 1.2 問題描述與研究方向 本論文題目主要以二階隱藏式馬可夫模型預測特定蛋白激酶磷酸化的位置，在這裡 我們的將此預測問題的定義是輸入一級蛋白質(primary structure)序列片段，然後系統要 能準確判斷出中間的位置是否會被特定的激酶磷酸化。 HMM 是一個富有彈性而且複雜的統計模型，HMM 發展自今以來已有相當多的變 形，從傳統的 HMM-1 到目前生物上大家常在用的 profile HMM，到更強大的二階隱藏 式馬可夫模型(second-order HMM, HMM-2)。本研究開發的 iHMM 主要是建構於 HMM-2 下的一個改良式 HMM。iHMM 在設計上使用一個簡單且快速的方法來幫助 Baum-Welch 演 算 法 加 速 收 斂 到 全 域 最 佳 解 。 以 及 運 用 貝 氏 資 訊 準 則 (Bayesian information criterion, BIC)的模型參數選擇方法，使 iHMM 在資料集不夠大時能盡量避 免眾所周知的過度適化(over fitting)問題。本論文會探討傳統 HMM-1 跟 HMMer 和本 篇論文主要使用的 iHMM 三種方法的特性與優缺點的比較，我們最後將用實驗結果來 呈現不受結構限制的完全圖的 iHMM 要比 first-order 的 Profile HMMs 的效能要好，不 過缺點就是在訓練階段時需要多一些計算時間。 圖 1.1 為我們的系統架構示意圖，本研究把 Phospho.ELM[5]資料庫與 Swiss-Prot 資料庫結合，將有蛋白激酶註解的資料依照蛋白激酶屬性 PKA、PKC、CDK 等分類十 八個資料集，然後分別用 iHMM 建立預測模型。本研究將完成蛋白激酶預測系統，並 且單純由一級序列來預測出特定蛋白激酶磷酸化的位置，當使用者輸入待測蛋白質序列 1.
(15) 時，預測程式將預測出哪些位置將會產生磷酸化。以期望未來生物學家藉此系統將能深 入了解蛋白質磷酸化過程之重要機制，亦或是以反應之重要特徵更深入研究。 S4 S0. S3 S1. S2. PKA S4 S0. S3 S1. S2. PKB S4 S0. S3 S1. S2. PKC S4 S0. S3 S1. S2. CDK Phospho.ELM database + Swiss-Prot database. S4 S0. S3 S1. S2. CK2 S4 S0. S3 S1. S2. CaM-KII S4 S0. S3 S1. S2. MAPK S4 S0. S3 S1. Src. 圖 1.1:系統架構示意圖. 2. S2.
(16) 1.3 章節概要 本論文章節概要如下： 第二章介紹什麼是蛋白激酶，以及目前蛋白激酶主要有哪些分類，還有它們的活 性，功能，結構做一個簡單的介紹。以及本研究所採用的激酶資料它的磷酸化的位置做 一個簡單的整理。 第三章主要探討 HMM 的原理與方法。並分別介紹 HMM-1 和 Profile HMMs 和 HMM-2，以及 HMM 常使用的符號與解說，而在此章最後會講到 HMM-1 與 HMM-2 的訓練演算法。 第四章要介紹本論文所提出的 iHMM 方法，以及一般在使用 HMM 上可能會遭遇 的困難與缺點，在這章節說明如何改良。以及 iHMM 如何將 HMM-2 與 BIC 結合，以 及它的特性與優點。 第五章為實驗結果與討論，包括生物資料的擷取來源，以及如何處理還有如何實 驗。最後我們將 HMM-1，HMMer，iHMM 三個方法做了正確率的分析與相關係數的比 較。並且以 PKA 資料為例，以視覺化呈現的方式展示三者方法的結構。 第六章為歸納了本論文的結論，並提出針對磷酸化問題未來還可以繼續研究的方向 與課題。. 3.
(17) 第二章 蛋白激酶基本介紹 2.1 蛋白激酶的分類 酵素上某些胺基酸的側鏈對於酵素反應的專一性及催化能力而言是很重要的，因為 在酵素固有的構形中，會迫使這些旁鏈集中在一起，形成活性區 (active site)。活性區包 含兩個重要的區域：受質結合區(binding domain)負責辨識並和受質結合，催化區(catalytic domain)則是在受質結合後負責催化反應。在某些酵素中，催化區就是受質結合區的一 部分，而在其它酵素，這兩個區域的結構就如同它們功能一樣有很大的差異。蛋白激酶 是激酶酵素，會在蛋白質的絲氨酸(serine, S)、蘇氨酸(threonine, T)或酪氨酸( tyrosine, Y) 殘基由化學修正加入磷酸根（磷酸化），直接或是間接來調控蛋白質的功能。蛋白激酶 作用絲氨酸(S)，蘇氨酸(T)和酪氨酸(Y)的磷酸化位置上的殘基具有保留性。它們有趣的 研究是在分子演化的部份。在序列和結構的基本上，這些酵素形成關係密切的超家族 (superfamily) 與組氨酸激酶和其它的磷酸轉移酵素截然不同 [6]。1988 年的時候 Steven K. Hanks 等人把真核細胞（eukaryotic cell）中的蛋白激酶分成五大類 AGC，CAMK， CMGC，PTK 和其它類[7, 8]。而在最近 G. Manning 等人的論文中[9]，則是把蛋白激酶 依照催化域結構的差別、序列相似度，和生物上的功能分成七大類，並建立了蛋白激酶 家族圖譜 ： I. AGC 激酶家族：底下主要代表的是 PKA、PKG 和 PKC 三種激酶。AGC 激酶是比 較寢向鹼性胺基酸的蛋白質，磷酸化主要做用於精氨酸(arginine, R)和賴氨酸(lysine, K)附近的絲氨酸(S)跟上蘇氨酸(T)上。 II. CAMK 激酶家族：底下包含了 CaMK、MARK 等許多激酶，CAMK 激酶是比較寢 向鹼性胺基酸的蛋白質，在演化上與 AGC 激酶家族相近。 III. CMGC 激酶家族：底下主要代表的是 CDK、MAPK、GSK3 和 CLK 等依賴於細胞 週期素的激酶，磷酸化主要做用於脯胺酸(proline, P)附近的絲氨酸(S)跟上蘇氨酸(T) 上。 IV. TK 激酶家族：底下包含了 EGFR、INSR 和 SRC 等許多作用於酪氨酸(Y) 的激酶。 V. TKL 激酶家族：底下包含了 MLK、LISK、IRAK、Raf、RIPK 和 STRK 等激酶。 VI. CK1 激酶家族：底下包含了 CK1、TTBK 和 VRK 等許多激酶。 VII. STE 激酶家族：底下包含了 MAP2K、MAP3K、MAP 4K 等以 MAPK 組成級聯反 應(cascade)的激酶。. 2.1 蛋白激酶的功能 蛋白激酶作用在許多重要的信息傳遞途徑上，催化作用的緊密調控對於真核細胞的 發展和維護至關重要，當它們活化時，不同蛋白激酶的催化區域採用顯著相似的結構， 相較之下，未活化的蛋白激酶的結晶結在面對與特定調控區域或蛋白質的相互作用允許 採用孑然不同的構型蛋白激酶的區域暴露出可易見的可塑性 [6]。 重要的蛋白質磷酸化在真核細胞信號中反應出實際蛋白激酶的區域 2%全部在的 真核基因中發現。特定絲氨酸(S)、蘇氨酸(T)或酪氨酸(Y)的磷酸化的空間與時間是控制 細胞成長和發展的關鍵，激酶活化在錯誤的地方或是在錯誤的時間會有災害的後果，常 常導致細胞轉型 (transformation) 或是癌症。其中以 TK 激酶家族與細胞信號通路有關， 對基本細胞功能比如分裂、分化和抗凋亡信號進行調節。這些激酶雖可促進細胞分裂， 但在非調控的啟動狀態下會導致各種腫瘤的形成。事實上已知的癌症基因中超過 70％就 4.
(18) 是與酪氨酸激酶有關[6, 10]。 當幾百個真核細胞蛋白激酶，全部含有相同的催化支架(catalytic scaffold)，一些非 常不一樣的調控機制已經演化成可以允許不同個體的家族輸入特殊訊號打開下游 (downstream) 的功能 [6]。 蛋白激酶是分子開關可以採用至少兩個極端的構型，分別是啟動狀態與停止狀態。 啟動狀態呈現最大活性，而停止狀態呈現最小活性。所有蛋白激酶的催化有同樣的反 應，就是將 ATP 上的磷酸基轉移到絲氨酸(S)、蘇氨酸(T)或酪氨酸(Y)的殘基上的氫氧 根群。這些活性他們全部只催化在結構非常相似的構型上。蛋白激酶在停止狀態時不需 要受到必須達到活化狀態的化學限制，而且不同類別的蛋白激酶演化成採用不同的方式 阻礙催化活性的構型來達成關閉狀態。. 2.2 蛋白激酶的結構 為了要說明活性區如何和與特定受質結合並促進這個已結合的受質發生化學變 化 ， 我 們 來 看 環 狀 單 磷 酸 腺 苷 依 賴 型 蛋 白 質 激 酶 (cyclic adenosine monophosphate dependent protein kinase)，現在通常被稱為 PKA，是真核細胞蛋白激酶領域中第一個用 X 射線結晶圖觀察到三級結構的例子，是第二訊息相關的酵素，主要做用在絲氨酸(S)，蘇 氨酸(T)上，並且在細胞內處理(cellular processes)，包括轉錄(transcription)，新陳代謝及 細胞的成長和凋亡扮演重要的角色。. 圖 2.1: PKA 三級結構 (來源是Protein Data Bank http://www.rcsb.org 編號 2CPK) 圖 2.1 來源是Protein Data Bank (http://www.rcsb.org) 編號 2CPK，是由Knighton 等 人 結 晶 出 來 的 PKA 三 級 結 構 [11] ， 此 圖 是 由 PyMol 工 具 產 生 ( http://www.pymol.org )。如圖 2.1，它的結構分成兩個小區域，小 N-端，或稱 N 葉(N lobe)，是一個五股的 β折板和一個明顯的 α 螺旋，稱為 helix αC。 C 葉(C lobe)是一 5.
(19) 個大且顯著地螺旋狀。ATP (adenosine triphosphate)是包在兩葉深入的裂縫中，位置位於 高保留的環 (loop) 連接者β1 和β2 兩股的下方。這個磷酸接合環，或 P 環，包含有保 留 的 富 含 glycine 基 序 (motif) (GxGxψG) 這 裡 的 ψ 通 常 是 酪 氨 酸 (Y) 或 是 苯 丙 氨 酸 (phenylalanine,F)。甘氨酸(glycine,G)殘基允許環靠 ATP 的磷酸非常接近且用骨架的交 互作用協調他們，表留的芳香族的側鏈覆蓋磷酸轉移位置。 在缺乏 ATP的情況下甘氨 酸(G)殘基使得 P 環變的非常易彎曲，事實上他促進小分子抑制物的接合。在環中，一 些抑制物由保留的芳香族殘基的相互作用引起大結構的扭曲。 多肽型受質接合在核苷酸前端的延伸的構型，與 ATP 的γ-磷酸基接近。中心區域 環是 “activation loop”，長通常約 20-30 個殘基，提供了一個平台給多肽型受質。 PKA 是最常見的蛋白激酶，當蛋白激酶活化時它的環是被磷酸化的。活化環的磷酸化使它穩 定在開放和延伸的構型中允許受質接合。 PKA 的活化區是位在催化次單元上一個長達 240 殘基的「激酶核心 (kinase core)」，激酶核心在所有蛋白激酶都具有高度保留性，主要功能是和受質(ATP 及目標 peptide 序列)結合，並且將 ATP 上的磷酸基轉移到絲氨酸(S)、蘇氨酸(T)或酪氨酸(Y) 的殘基上。激酶核心是由一個大的結構域及一個小結構域所組成，介於中間的是一個深 的裂縫，而兩個結構域上都有構成活性區的殘基。 蛋白激酶通常保持停止狀態，而且催化區域的活化通常是包埋在多層的控制中，範 圍從異位效應 (allosteric effectors) 的相互結合到亞細胞內位置轉移。當激酶在不活化與 活化兩個狀態之間切換時，活化環有能力進行較大的構型改變。例如，胰島素受器激酶 (isuline receptor kinase, IRK) 是一個被三個 Y 殘基磷酸化與它的活化環活化。在沒有磷 酸化狀態，活化環塌陷進入活化位置，凍結核苷酸和多肽型受質兩個接合。在磷酸化上 面，活化環延著催化中心移動並且採用一個允許接合與催化的構型。在許多激酶，包括 IRK，在活化環的 N 端因為曲軸上的移動造成這些構型的改變，造成在 Asp-Phe-Gly 序 列保留適當的方向有活化位置。. 2.3 特定蛋白激酶磷酸化的位置 蛋白激酶是激酶酵素，所以他與其它酵素一樣在活性區俱有專一性辨識與結合基質 的能力。我們將不同激酶特定磷酸化位置絲氨酸(S)、蘇氨酸(T)或酪氨酸(Y)的地方取出 左右兩邊範圍各 7 個殘基，總長 15 建出序列圖案(sequence logos)，放在附錄中。這些序 列圖案是由 WebLogo 工具 (http://weblogo.berkeley.edu [12] )產生出來。序列圖案是一 種基由 Shannon 訊息概念圖示化的技術，是用來呈現已經序列排列好常用的慣例，每一 個字母的高度是與相對胺基酸出現頻率成比例，表示在那個位置的保留程度，如果某個 位置全部都只出現那一個胺基酸，則 Shannon 訊息最高就是 log 2 20 = 4.32 bits，相反的 如果那個位置每個胺基酸出現的頻率都一樣，則就是 0 bits，這個觀念由 Schneider 在 1990 做了介紹[13]。我們將 13 個常見的激酶的位置保留性列在下表，紅色字體(S/T) 表示在那個位置磷酸化，X 表示這個位置可以是任一個胺基酸。舉例來說我們可以觀查 到 PKA 在位置 -2 和 -3 的地方是有很高頻率會出現精氨酸(arginine, R) 和賴氨酸 (lysine, K)，而 PKB 在位置 -3 地方是有很高頻率會出現精氨酸(R)，在位置 -5 的地方 有很高頻率會出現精氨酸(R) 和賴氨酸(K)，這些對於蛋白激酶特定磷酸化位置的預測是 一個很重要的特徵。. 6.
(20) 表 2.1:蛋白激酶磷酸化的位置保留性 激酶 PKA PKB PKC PKG CaM-KII CK1 CK2 CDK MAPK ATM EGFR INSR SRC. 位置保留性 (R/K)-(R/K)-X-(S/T) R-X-R-X-X-(S/T) R-R-X-(S/T)-X-(R/K) R-(R/K)-X-(S/T) R-X-X-(S/T) S-X-X-(S/T)-X-E-S (S/T)-(D/E)-(D/E)-(D/E) (S/T)-P-X-(R/K) P-X-(S/T)-P S-Q E-E-X-X-Y-V D-Y-M-X-M E-(E/D)-X-X-Y. 7.
(21) 第三章 方法 在眾多的預測模型中，像是決策樹(decision tree)、支持向量機器(support vector machine, SVM)以及類神經網路（neural network）等等，這些方法在解決不同類型的分 類問題時都有各自的優缺點，然而我們這裡主要採用 HMM，主要是因為我們今天要預 測的磷酸化問題與它周遭的胺基酸關係很密切，而 HMM 具備有良好處理時間序列的 能力，把序列特徵中相關的胺基酸以時間來區分，找出某段時間內的可預期行為特徵， 正好符合這些特性，因此我這裡主要探討以 HMM 為主的預測方法，分別是傳統的 HMM-1 和 Profile HMM 和 HMM-2。. 3.1 First-order HMM HMMs 是一個將有限狀態機與機率結合的一個統計模型，專門用來描述連續符號 序列的條件概率。HMM 技術在語音辨識領域中已經有相當多年的歷史，是馬爾可夫模 型的擴展。該模型由兩個隨機變量序列組成：一個是觀測不到的馬爾可夫鏈，另一個是 可以觀測到的隨機序列。在數學上完整描述一個 HMM 是以起始狀態機率，狀態轉移 機率，和狀態觀測機率三種機率矩陣所組成，分別以π、A 、B 表示，並且以 λ≡(A,B,π) 表是 HMM 所有的參數的集合[14]： (1) 起始狀態機率 起始狀態機率是一個隱藏的隨機程序，它決定以某個狀態開始的機率，在數學 N. 上一般是 π i ≡ P(q1 = S i ) 來表示，並且要遵循 ∑ π i = 1 這個機率限制。 i =1. (2) 狀態轉移機率 狀態轉移機率是一個隱藏的隨機程序，它決定停留在原狀態或是遷移到下個狀 態的可能性，狀態之間的連線表示它們的因果關係。在數學上一般是 A≡[aij]NxN 矩陣來 N. 表示，其中 aij ≡ P(qt = S j | qt −1 = S i ) ，所有狀態都要遵循 ∑ aij = 1 這個機率限制，如： j =1. 從某一狀態出發轉移到下個狀態的機率總和為 1。. (3) 狀態觀測機率 狀態觀測機率則描述在每個狀態觀察到某個現象的機率，在數學上一般是 N. B ≡ {b j (k )} 來表示，其中 b j (k ) ≡ P (v k at t | q t = S j ) ，並且要遵循 ∑ b j (k ) = 1 這個機率限 k =1. 制。. 3.2 Profile HMM Profile HMMs 是一種「機率模型」的多重序列對齊技術(statistical models of multiple sequence alignments)，它捕捉在多重序列比對後每一行的位置特異性(position - specific) 的資訊，包括有關保留性(conserve)，以及殘基的相似性。Profile HMMs 是由 Anders Krogh，David Haussler 等人在美國加州聖塔克魯斯大學首先引介給計算生物學 [15]， 然而在生物學上 Krogh/Haussler 等人並不是第一個使用 HMM，從歷史上來看，HMM 的 使 用 最 早 可 以 追 朔 到 1989 年 ， Churchill 為 第 一 位 將 HMM 用 在 異 質 基 因 型 8.
(22) (heterogeneous DNA) 序列的建模上[16]。生物學家一旦拿到未知功能的 DNA 或蛋白質 序列時，最常做的事情就是從資料庫中進行序列比對，而 Krogh 的論文提出了一個震 撼性的觀點，因為 HMM 技術是非常適合使用多重序列比對以廣泛性的 profile 方法進 行資料庫搜尋取代單一序列搜尋。如圖 3.1，一般而言，Profile HMMs 有三個主要的狀 態分別是 match (M)，insertion (I)和 deletion (D)，而機率的參數是從多重序列比對而 來。HMMer 是其中有名的 profile hidden Markov model 一個工具程式[4]。 J2. J3. J4 F. T. O. N2. N3. N4. N5. E2. E3. E4. E5. D. U. C. K. 圖 3.1: Plan7 結構. HMM 是一個對於內文相關的序列統計分析相當有力的方法。Profile HMMs 參數 的權重給定來自於多重序列比對，其演算法與 Plan7 架構的搭配可以說非常迅速而且乾 淨俐落，然而 profile HMMs 由於參數的權重給定演算法和架構的限制有一些缺點，其 中一個缺點是 HMMer 在訓練時無法直接處理長度不一的字串，簡單來說就是 HMMer 輸入訓練的資料必須是全部字串要都一樣長，若不一樣長則需要事先做多重序列比對， 插入 gap 將所有字串對齊後才可以丟進去訓練，而在本論文提出的 iHMM 是不須要有 這個限制，本論文提出的 HMM 訓練時可以接受長度不一的字串。另一個缺點是 Profile HMMs 是屬於 HMM-1 不能捕捉到微弱的相關訊號(high order dependency signals)。上 圖 3.1 來說，從 M3 走到 M4 在 HMM-1 裡是已經把資訊壓縮掉，也就是不考慮前狀 態，若我們把 M3→M4 的訊號放大來看，到底他之前的路徑 I2→M3→M4 還是 M2→M3→M4 還是 D2→M3→M4，然而在資料夠複雜的情況下，這些資訊的保留有時 是很有用的。. 3.3 Second-order HMM HMM-2 是 HMM-1 的擴展，它與 HMM-1 的差別在於狀態轉移機率，在 HMM-1 中 狀態轉移機率只有目前狀態與前一個狀態的關係。而在 HMM-2 中狀態轉移機率是目前 狀 態 與 前 二 個 狀 態 的 關 係 ， 它 的 狀 態 轉 移 機 率 定 義 是 ： A ≡ {a ijk } ， aijk ≡ P (q t = S k | qt −1 = S j , qt − 2 = S i ) 此狀態所有 1≦i≦N, 1≦j≦N 轉移機率都要符合 N. ∑a k =1. ijk. = 1 這個機率限制，這裡的 N 表示狀態數目，qt 表示在時間 t 的狀態。如果以維度. 來看，HMM-1 的狀態轉移機率是一個二維的矩陣 N × N ，而 HMM-2 狀態轉移機率是一 9.
(23) 個三維的矩陣 N × N × N 。每一個 HMM-2 其實等價於二維狀態 SxS 空間展開的 HMM-1，但是它不像 HMM-1 會增加狀態數目，如果有足夠的 order ，訓練資料所沒 有走過的路徑，在測試資料中也將不會出現，這可以有助於區別錯資料的能力。如下圖， 在這裡我引用 1997 年 Jean-Francois Mari 等人在論文中 [17] 所舉的例子，圖 3.2 是 HMM-2 的結構，圖 3.3 是與圖 3.2 結構等價的 HMM-1，這兩個結構的符號觀測機率有 其相同對印，也就是 S01 等於 S1，S12 與 S22 等於 S2，S23 與 S33 等於 S3，S34 等於 S4。待會我們會發現利用更高階的轉移機率來減少狀態數目，在訓練和預測時會有他的 優勢。. S0. S1. S2. S4. S3. 圖 3.2:原始 HMM-2 S01. a012. a123. S12. 1-a123. 1-a222. S23. 1-a234. S22. S33. a222. a333. a234. S34 1-a333. 圖 3.3: HMM-1 等價模型 在語音辨認上已有越來越多的論文指出連續型的 HMM-2 的效能要比 HMM-1 好 的多[18, 19]。生物資料這幾年正迅速的增加，higher-order 的統計模型也將越來越重要。 當資料足夠多時 higher-order HMM 可以解決 HMM-1 所不能辦到的事情，因為 higher-order HMM 的轉移機率不只會參考前一個走過的狀態，而且還會考量過去之前二 個以上所走過的狀態，higher-order HMM 可以比傳統 HMM-1 以更少的參數表達更複 雜的模型，我們希望透過這種更複雜的前後文資訊(contextual information)幫助我們找出 哪些蛋白激酶在哪些位置起作用。. 3.4 分數計算 分數的計算是用來做為未知蛋白質序列判定接受或拒絕的標準，高於門檻值即接 受，低於門檻值即拒絕。若 HMM 的參數 λ 全部都已知，給一條序列 O=<O1,O2,…,OT> 與相對應一條所走的路徑 Q=<q1,q2,…,qT>，在 HMM-1 中我們可以很快的計算這條路徑 下的機率：. P(Q | λ ) = π q1 aq1q2 aq2 q3 " aqT −1qT. (1). T. P (Q, O | λ ) = π qq bq1 (O1 )∏ aqt −1qt bqt (Ot ). (2). t =2. 10.
(24) 同樣的我們也可以很快的算出 HMM-2 一條路徑下的機率：. P(Q | λ ) = π q1 aq1q2 aq1q2 q3 " aqT −2 qT −1qT. (3) T. P (Q, O | λ ) = π qq bq1 (O1 )a q1q2 bq2 (O2 )∏ a qt − 2 qt −1qt bqt (Ot ). (4). t =3. 然而我們要把一條序列所有可能走過的路徑的機率累加起來才是我們要的 P(O | λ ) 機率，此機率的計算是 HMM 中第一個基本問題，要計算這個機率我們可以由動態規 劃(dynamic programming ) 方法來快速求解，alpha 函數的計算待會在 3.5 節會有詳細的 算法：. P (O | λ ) = ∑ P(O | Q, λ ) P(Q | λ ). (5). Q N. P (O | λ ) = ∑ α T (i ). (6). i =1. 這裡我們在計算分數時與傳統的 Profile HMMs 一樣 HMM-1 與 iHMM 均是採用 log-odds的分數算法，log-odds 分數是用來評估這條序列相對於虛無模型 (null model) 多少機率是由HMM產生。這個分數又稱做為log-likelihood ratio，這個公式描述在公式 (7)。在此公式中 P (O | λ P ) 是序列在positive HMM 中的機率。 P(O | Null ) 是序列在自然 界中的機率。在 HMMer 中，虛無模型是如同簡化一個狀態的 HMM，狀態裡的符號觀 測機率是由已經統計好的自然界胺基酸出現的頻率。在我們的 iHMM 中所使用的自然 機率，例在表 3.1，此自然機率表來自於 HMMer 原始碼。. score = log. P(O | λ ) P(O | Null ). (7). 表 3.1:自然界胺基酸出現的頻率 胺基酸 代號 A C D E F. 自然機率. 0.083534 0.014532 0.052473 0.061821 0.039814. 胺基酸 代號 G H I K L. 胺基酸 代號 M N P Q R. 自然機率. 0.070123 0.022714 0.0572 0.052905 0.098249. 自然機率. 0.023504 0.041622 0.05025 0.040249 0.055901. 胺基酸 代號 S T V W Y. 自然機率. 0.070339 0.055862 0.065857 0.013155 0.029897. score 輸出是一個機率分數，表示 HMM 的機率相對於自然界中的機率兩者的距 離，高的分數表示序列是傾向是會被磷酸化，反之底的分數表示序列是傾向是不會被磷 酸化。如果一個新的序列進來得到的分數比給定的門檻值還高，則我們說這條序列是 positive，反之這條序列是 negative 。. 3.5 訓練演算法 11.
(25) 我們令 Σ 是胺基酸字母的集合. |Σ| 是胺基酸字母的總數, 假設胺基酸的號碼是從 1 到 20，則標示的 |Σ| = 20. 我們令 positive 抽樣本 O P ≡ (O P (1) , O P ( 2 ) , O P ( 3) , " , O P ( k ) ) 表示 k 條序列，並且假設每條序列都相互獨立並且具有相同分配 (independent and identically distributed , i.i.d.)，其中 O P ( k ) ≡< O1P ( k ) , O2P ( k ) , O3P ( k ) ,", OTkP ( k ) > 表示第 k 條由字 母組成長度為 Tk 的觀查序列 (observation sequence)，其中 OTkP ( k ) ∈ Σ。而在己知的 HMM 參數 λ≡(A,B,π)之下，得到的聯合機率密度函數，稱之為 likelihood N. P (O | λ ) = ∏ P(Oi | λ ). (8). i =1. 最佳化機率模型的參數來精準描述所觀查到的序列是非常重要的。在這裡訓練的問 題是要挑一組 HMM 參數 λ≡(A,B,π) 使得最大化抽樣資料的概似率(likelihood) ，最大 概似率是求得一估計式使其抽樣資料的聯合機率密度函數的值為最大。正規的定義是 N. λ * = arg max ∏ P(O P (i ) | λ ) λ. (9). i =1. 其中期望值最大化演算法(expection maximization algorithm, EM algorithm)是一個廣 泛用來解決最大概似率估計的(maximum likelihood estimation, MLE) 技術。它的廣義步 驟如下： 步驟 1:隨機挑選一個初始參數 λ. [. 步驟 2:E-step，求出輔助函數 (auxiliary function) Q(λ , λ ) = ∑ P(Q | O, λ ) log P(O, Q | λ ). ]. Q. 步驟 3:M-step，根據訓練資料以步驟 2 的輔助函數計算 λ 期望值， λ = max Q(λ , λ ) λ. 步驟 4:以新的參數取代舊的參數 λ = λ ，然後重複步驟 2 跟 3，直到收斂為止. EM 演算法目前依然是解 MLE 最有效的方法，它保證每一次的疊代（iteration）後 一定是 P (O | λˆ ) ≥ P(O | λ ) 。在 HMM-1 裡，步驟 2 中輔助函數已被 Baum 和它的同僚 所解出。它的詳細演算法在 Rabiner [14]已有很好的整理，在這裡我直接引用 Rabiner 論 文裡所整理的 Baum-Welch 參數重估演算法，Baum-Welch 演算法主要是運用 forward probability (alpha 函數) 跟 backward probability (beta 函數) 來進行參數重估，forward probability 跟 backward probability 都可以用動態規劃( dynamic programming ) 方法來 求解，為了要方便計算我們另外還要 gamma 函數、跟 xi 函數。它們的完整公式如下 (1) alpha 函數 說明：alpha 函數，通稱為 forward probability，表示在模型λ底下，看到 O1O2…Ot 並且於時間 t 時停在狀態 i 之下的機率。而我們剛才提到的 P(O | λ ) 的機率可以由 N. 此 forward 演算法求得，即 P (O | λ ) = ∑ α T (i ) i =1. 定義： α t (i ) = P(O1O2 "Ot , qt = si | λ ) 初始值： α1 (i ) = π i bi (O1 ) ,1≦i≦N 12.
(26) N. 遞迴式： α t +1 ( j ) = ∑ α t (i )aij b j (Ot +1 ) i =1. (2) beta 函數 說明：beta 函數，通稱為 backward probability，表示在模型λ底下，已知時間 t 時停 在狀態 i 之下，看到 Ot+1Ot+2…OT 的機率。 定義： β t (i ) = P(Ot +1Ot + 2 "OT | qt = si , λ ) 初始值： β T (i ) = 1 ,1≦i≦N N. 遞迴式： β t (i ) = ∑ aij b j (Ot +1 ) β t +1 ( j ) ,t=T-1,T-2,…,1,1≦i≦N j =1. (3) gamma 函數 說明：gamma 函數，通稱為 forward-backward probability，表示在模型λ底下，時間 t 時停在狀態 i 之下的機率。 定義： γ t (i ) = P(qt = S i | O, λ ) α (i ) β t (i) 公式： γ t (i ) = t P (O | λ ) (4) xi 函數 說明：表示在模型λ底下，在時間 t 由狀態 i 轉移到狀態 j 的機率。 定義： ξ t (i, j ) = P(qt = s j , qt +1 = s k | O, λ ) 公式： ξ t (i, j ) =. α t (i, j )aij bk (Ot +1 ) β t +1 ( j , k ) P (O | λ ). 而由上面式子，我面又可以整理如下式 N. γ t (i ) = ∑ ξ t (i, j ) j =1. Baum-Welch 參數重估演算法，即藉由上例 alpha、beta、gamma 跟 xi 函數，進行 以下計算，並且以疊代方式，不斷重估 HMM 的參數，直到收斂為止。. πi =. γ 1 (i ) N. ∑γ i =1. 1. (10). (i ). T −1. aij =. ∑ ξ (i, j ) t =1 T −1. t. ∑γ. t. t =1. (i ). T. ∑γ. b j (k ) =. t =1 s .t .Ot = vk T. ∑γ t =1. (11). t. t. ( j) (12). ( j). 13.
(27) 在 HMM-2 這裡我們使用推廣的 Baum-Welch 演算法，來做為我們訓練演算法 [20]。為了要方便計算我們另外還要推廣的 gamma 函數、eta 函數跟 xi 函數。它們的 完整公式如下. (1) 推廣的 alpha 函數 說明：alpha 函數，通稱為 forward probability，表示在模型λ底下，看到 O1O2…Ot 並且於時間 t-1 跟時間 t 時分別停在狀態 j 跟狀態 k 之下的機率 定義： α t ( j , k ) = P(O1O2 " Ot , qt −1 = s j , qt = s k | λ ) 初始值： α 1 (i, j ) = π i bi (O1 )aij b j (O2 ) ,1≦i≦N N. 遞迴式： α t ( j , k ) = ∑ α t −1 ( j , k )aijk bk (Ot +1 ) i =1. (2) 推廣的 beta 函數 說明：beta 函數，通稱為 backward probability，表示在模型λ底下，已知時間 t-1 和 時間 t 分別停在狀態 i 跟狀態 j 之下，看到 Ot+1Ot+2…OT 的機率 定義： β t (i, j ) = P (Ot +1Ot + 2 " OT | qt =1 = si , qt = s j , λ ) 初始值： β T (i, j ) = 1 N. 遞迴式： β t (i, j ) = ∑ β t +1 ( j , k )aijk bk (Ot +1 ) i =1. (3) 推廣的 eta 函數 定義： η t (i, j , k ) = P(qt −1 = si , qt = s j , qt +1 = s k | O, λ ) 公式： η t (i, j , k ) =. α t (i, j )aijk bk (Ot +1 ) β t +1 ( j , k ) P(O | λ ). (4) 推廣的 xi 函數 說明：表示在模型λ底下，在時間 t 由狀態 i 轉移到狀態 j 的機率。 定義： ξ t (i, j ) = P(qt = s j , qt +1 = s k | O, λ ) N. 公式： ξ t (i, j ) = ∑η t (i, j , k ) k =1. (5) 推廣的 gamma 函數 說明：gamma 函數，通稱為 forward-backward probability，表示在模型 λ 底下，時間 t 時停在狀態 i 之下的機率。 N. 公式： γ t (i ) = ∑ ξ t (i, j ) j =1. 推廣的 Baum-Welch 參數重估演算法，即藉由上例推廣的 alpha、beta、gamma 、 xi 跟 eta 函數，進行以下計算，並且以疊代方式，不斷重估 HMM 的參數，直到收斂 為止。. π i = γ 1 (i ). (13) 14.
(28) aijk =. ∑ η (i, j, k ) ∑ η (i, j, k ) t. t. k ,t. T. ∑γ. b j (k ) =. (14). t. t =1 s .t .Ot = vk T. ∑γ t =1. t. t. ( j) (15). ( j). 15.
(29) 第四章 提出 iHMM 方法論 4.1 想法起源 HMMer 的演算法無論是在參數的權重給定演算法或是使用的 Plant7 結構，都具有 其生物意義，然而我們在實驗中發現，HMMer 在小量資料時，常常會選到太長的結構， 結果因為狀態太多，資料量太少，而導致過適化。用有限的資訊對未知現象作判斷與分 析乃是統計學所面臨的必然問題。而 HMM-2 的特性正好可以利用更多的轉移機率來 進行節省狀態的使用。我提出 iHMM 是將 HMM-2 加入二種機制，第一是並配合簡易 解參數最佳化問題，第二是加上 BIC 的機制來防止訓練中過適化的問題。另外我們 iHMM 在設計上可以很方便進行第二類資料的訓練與分數計算。 本論文所提出的 iHMM 在訓練階段時不需要多重序列比對而且沒有結構限制，缺 點是它需要花一些時間來訓練。然而，現在電腦速度這幾年越來越快，計算速度是越來 越不重要的。我們提出的 iHMM 要比 profile HMMs 有更多的優點。其中一個優點是 由完全圖可以捕捉到更多序列內文微弱的相關訊號。另一個優點是在少量的資料中， iHMM 可以比 profile HMMs 能更防止過適化的問題，基於這些原因，我們在測試資料 上可以提升辨識率。. 4.2 iHMM 之設計 Baum-Welch 演算法是一個眾所都知的用來解 HMM 問題的演算法，然而，在解 空間中可能存在許多的區域最佳解。 Baum-Welch 演算法不保證可以找得到全域最佳 解。Baum-Welch 演算法對於一開始給定的初始點是很重要的，因為 Baum-Welch 演算 法會與一開始參數的設定而掉入不同的區域最佳解 本論文提出一個族群搜尋式的 Baum-Welch 演算法。其主要的想法是使用多點式搜 尋來取代單點搜尋。它將有更多跳脫區域最佳解的能力，這個演算法的虛擬程式碼如下：. 16.
(30) Population-based Baum-Welch algorithm 1: Input: Training set, initial barrier ω, and population size k 2: Output: HMM model 3: Begin 4: 5: S={λ1, λ2, …, λk } 6: 7: for i=1 to |S| 8: randomize S(i) parameters 9: end i 10: while ( |S|>1) 11: { 12: δ← 0 13: for i=1 to |S| 14: do { 15: 16: S(i) do E-step and M-step 17: Δ ← (new likelihood of S(i)- old likelihood of S(i)) 18: } 19: while(Δ not converge to ω); 20: δ←δ + Δ; 21: end i 22: ω←δ 23: sort S by likelihood 24: remove the last half λi of S, i= |S|/2 to |S| 25: } 26: return S(1) 27: 28: End. 17.
(31) 這裡我們對演算法做一個描述，S 是一組 HMMs 候選解的集合，k 是一開始要灑 的粒子數，在 7 到 9 行中，我們是用亂數隨機在解空間進行灑點，k 越大將越可以找 到全域最佳解，但也相對要花越多時間搜尋，k 值是完全依據使用者自己選定。在 10 到 25 行 while loop 的迴圈每一次的重復，我每放入柵欄 ω 來同時控制所有點的行動。所 有點的都由 Baum-Welch 演算法移動，並且當 likelihood 收斂到柵欄 ω 值則停下來。在 演算法的第 16 行，E-step 和 M-step 是來自於 EM 演算法，它的補助函數 Q-step 己 經由 Baum 和它的同僚所解出，詳細算法可以參考第三章 3.5 節的訓練演算法的部份， 在演算法的第 22 行，我們每一代都會調控柵欄 ω 值，在演算法的第 24 行，為了加快 搜尋速度，我們以 likelihood 做排名，然後只保留前半部目前較好的點，拿掉後面一半 比較差的點，以減少計算時間，因為這些較差的點將來要再比最好的還要好的機率很 小。然後一直重復直到只剩下只有一個點。換句話說，每一個點必須為了生存而競爭， 這個如同達爾文的適者生存，不適者淘汰的觀念。這個演算法可以很方便改成平行處 理，我們相信這個方法是很快的，可以幫助加速收斂到全域最佳解，而且很穩定。 在圖形識別(pattern recognition)課題上，過適化(over fitting)的現象一直是最令人棘 手的問題，也是嚴重影響辨識率的主要原因之一。要挑選多少個狀態下去訓練是一個很 難的問題，因為 HMM 會隨著狀態數目越多而會越有過適化現像，如下圖 4.1，我們在 訓練時只能看的見訓練資料的正確率，就是圖 4.1 中實線的部份，虛線的變化是我們看 不到的部份，然而我們必須要選擇正確的模型才可以最小化錯誤預測能力。 Accuracy 1.0 0.8 0.6 Training set 0.4 Test set 0.2 1. 2. 3. 4. 5. 6. 7. 8. 9. State number 圖 4.1:狀態數目與辨識率之間的關係的示意圖. 過適化的現象原因出在於當我們的訓練資料抽樣不平均，或是不夠足以代表母體 樣本分佈，或著是資料量少於相對於模型的自由度，MLE的方法會讓我們很容易 fit 某 一類的訓練資料，這就是為什麼說HMM通常會特別偏愛某些條序列。原因是雖然都是 在同類 (positive) 的序列，但就是有些序列從HMM出來的機率特別的高，他們的分數 會把孤立的序列分數拉下來，會導致系統不對稱，使得訓練 出來的 HMM profile 將來 會特別偏愛某些條序列。有些論文指出用 max min P (O i | λ ) 的方法可能會比傳統 λ. 18. i.
(32) N. max ∏ P(O i | λ ) 的方法來的好，以用來避免過適化的問題[15, 21]。在 HMM 裡為了 λ. i =1. 避免MLE的方法因為少量的訓練資料遭到過適化的問題。另一個方式則是採用最大事後 機率(maximum a posteriori probability estimation, MAP)，取代MLE的方法，作為HMM 之 評估準則，假設每個線性組合係數都是隨機且其事前機率為一個高斯分佈，我們可以推 導出一組最大事後機率估測值。它的優點是可以避免因資訊的不足，而使得估測產生過 大的誤差。根據貝氏定理，我們的事後機率分佈如公式 (16):. P (λ | O ) =. P(O | λ ) P(λ ) P(O). (16). 不幸的是，在公式 (16) 這個 P(λ) 事前機率很難求得。而在 HMMer裡則使用混入 Dirichlet 事前機率[22]，並使用 Blocks9 表裡的芳香族、脂肪族，有極性，沒極性．．． 等九個自然界胺基酸事前機率單元，這個混入Dirichlet 事前機率的效用如同模擬了很多 與訓練資料可能的演化樣本，用來平滑化(smooth) HMM的參數，使得訓練出來的 HMM 不會遭到過適化。 廣義線性模型（generalized degrees of freedom）是傳統線性模型的延伸，使用 maximum likelihood 方法來估計參數，並假設反應變數的機率分佈屬於指數家族（exponential family），將廣義自由度的想法由常態分布推廣到更廣泛的指數族群分布，並研究廣義 線性模型（generalized linear model）的變數選取問題。不同於傳統的狀態數目選取方式， 我們不直接選取狀態數目，而是以廣義自由度方法選取具不同懲罰(penalty)參數的資訊 準則（information criterion），然後再依照所選出的準則執行 HMM 狀態數目選取。我 們用自由度的觀念修正變數選取過程對推論結果所產生的偏差。這個方法可隨資料背後 的真實模型而自動調適，因而兼顧。 基 本 上 有 兩 種 公 認 標 準 的 衡 量 方 法 叫 做 艾 凱 克 資 訊 準 則 (Akaike information criterion, AIC) [23] 和貝氏資訊準則(Bayesian information criterion, BIC) [24]。而在磷酸 化預測問題中，本論文首先提出離散型的 HMM-2 加上 BIC 模型選擇方法，來取代傳 統 HMM-1 以解決生物資訊的問題。AIC與 BIC在公式上很相像，差別在於採用不同的 懲罰方式，BIC在懲罰的地方則是比AIC多了一個訓練的資料個數的代價倍率，所以在 這裡我們採用 BIC 來做為我們挑選 HMM 狀態數目的 衡量準 則。BIC 定義如式(17)：. BIC = − log( L) + 0.5 × K × log( N ). (17). L 是統 計 模型的 likelihood 值 統 計 模型總共所使用的參數個數，或自由度 N 是所提供訓練的資料個數. 19.
(33) BIC. 1. 2. 3. 4. 5. 6. 7. 8. 9. State number 圖 4.2:狀態數目與 BIC 的示意圖. BIC 主要概念是運用少量參數，來建立一個統計模型，使它更具有強健性。如上圖 4.2，我們的狀態數目選取，則是採用 BIC 最低點的位置。用 BIC 來選擇模型，除了用 在 HMM-2 上，在許多文獻上 HMM-1 用 BIC 來選擇模型也有一也不錯的例子，像是 要最佳化以微陣列基因表現的 HMM 模型[25]，以及最佳化預測蛋白質二級結構的 HMM 模型[26]等等。另外 1999 年 Brand 也發表了一篇論文[27]，他在文章提到說，將結構的 一個參數趨近於零，可以減少模型一個自由度，並且減低模型的亂度(entropy)，模型的 亂度越低，越有助於抵抗過適化的問題。同時參數使用的越少，將可以使得模型越可解 讀。模型參數量的使用，其實跟信息熵有關，就好比將一份資料壓縮，我們希望在仍足 以表達某程度資訊內容條件下，進行減少該模型的參數使用。 在語音辨識中，常會因為聲音受到背景雜音干擾而降低辨識率，於是人們常需要從 不的錄音環境，如教室，車內或是室外中，將背景雜音與目標語者的聲音分離出來，所 以要提高辨識率的一個方法就是另外再從背景聲音建一個模型。生物資料中也有相同的 情況，positive 資料就好比是目標語者的聲音，而negative 資料就好比是背景雜音。理 論上來說，隨著訓練樣本類別數的增加，辨識結果可以更加近似於實際結果。假設從 positive 建立iHMM 統計模型以符號λP來表示，而從negative 資料建立iHMM 統計模 型，以符號λN來表示，如果 negative 資料足夠亂，則 P (Oi | λN ) 就成為常數，我們就只 要最大化 ∑ logP(Oi | λP ) 即可。相反的如果negative 資料不夠亂，則訓練negative的資料 i. 模型則會有用。兩個模型的分離成度可以用相對熵(relative entropy)來解釋，即有名的 Kullback-Leibler (KL)量度，KL是用以量測兩個機率分佈的分離成度 (divergence)，它的 定義在式(18)。. KL(λP , λN ) = H (λP || λN ) = ∑ P(Oi | λP ) log i. P(Oi | λP ) P(Oi | λN ). (18). 然而KL並不是真正的距離公式，因為他在數學上不滿足交換率，如式 (19)，而在 Matthew A. Siegler 的[28]論文中，則是把 KL 寫成對稱形式，如式(20)，即成為兩個模 型的距離。我們就可以用 KL2 來衡量positive 模型是否與negative 模型是否有存在區 別。 20.
(34) KL(λP || λN ) ≠ KL(λN || λP ). (19). KL 2(λP , λN ) = KL(λP , λN ) + KL(λN , λP ). (20). 21.
數據
+7
相關文件
Using a one-factor higher-order item response theory (HO-IRT) model formulation, it is pos- ited that an examinee’s performance in each domain is accounted for by a
volume suppressed mass: (TeV) 2 /M P ∼ 10 −4 eV → mm range can be experimentally tested for any number of extra dimensions - Light U(1) gauge bosons: no derivative couplings. =>
Using this formalism we derive an exact differential equation for the partition function of two-dimensional gravity as a function of the string coupling constant that governs the
• We have found a plausible model in AdS/CFT that captures some essential features of the quantum hall eect as well as a possible experimental prediction. • Even if this is the
For a deep NNet for written character recognition from raw pixels, which type of features are more likely extracted after the first hidden layer.
“Since our classification problem is essentially a multi-label task, during the prediction procedure, we assume that the number of labels for the unlabeled nodes is already known
The angle descriptor is proposed as the exterior feature of 3D model by using the angle information on the surface of the 3D model.. First, a 3D model is represented
D.Wilcox, “A hidden Markov model framework for video segmentation using audio and image features,” in Proceedings of the 1998 IEEE Internation Conference on Acoustics, Speech,
