用BioJava開發基于隱馬爾科夫模型的應用

2019-11-18 12:05:07

字體：大中小

來源：轉載

供稿：網友

　　1．隱馬爾科夫模型在生物信息學中的應用
　　作為一種強大的概率模型，隱馬爾科夫模型被廣泛地應用到語音識別、計算機視覺等領域。而在計算分子生物學和生物信息學領域，隱馬爾科夫模型同樣有著廣泛的應用。例如，基于隱馬爾科夫模型的基因發現和序列比對算法已經被開發出來。
　　
　　在生物信息學的教科書中，使用隱馬爾科夫模型的一個典型的例子是CpG島的發現。CpG島是一段短的DNA序列，該序列中的CG含量比整個基因組中的平均含量要高。許多基因的上游有有CpG島，例如，56%的人類基因上游含有CpG島[1]。由于CpG島與基因的非凡關系，假如能準確地定位一段長序列中的CpG島，則對于基因發現是有重要意義的。
　　
　　本文的下面部分將介紹隱馬爾科夫模型，并介紹一個具體問題的Biojava實現。
　　
　　2．馬爾科夫鏈與隱馬爾科夫模型
　　馬爾科夫鏈
　　馬爾科夫鏈是一個狀態序列，對于每個由狀態i和狀態j組成的狀態對，有轉換概率aij，∑jaij=1[2]。馬爾科夫鏈的要害特性是狀態i出現的概率依靠且僅依靠其前一個狀態[3]。
　　
　　隱馬爾科夫模型
　　隱馬爾科夫模型（以下簡稱HMM）是在馬爾科夫鏈基礎上發展而來。
　　實際問題往往比馬爾科夫鏈所描述的更復雜，觀察到的事件并不與狀態一一對應，而是通過一組概率分布對應，這樣的模型即為HMM[4]。
　　
　　HMM的強大之處在于在觀察到的事件與內在的狀態間建立了一種概率模型，使用Vertbi算法，能夠根據一個給定的觀察序列和一個模型，在最佳的意義上確定內部狀態序列[4]。也就是說，根據可觀察的事件序列，來推測不可觀察的內部狀態序列。
　　
　　文獻[3，4]中有關于馬爾科夫鏈和HMM的更具體描述。
　　
　　3．作弊的色子問題
　　作弊色子問題是學習HMM的典型例子。
　　
　　假設有一個賭場，在絕大多數情況下，他們時候正常的色子，但是有些情況下，他們也使用做過弊的（Loaded）色子。當使用正常的色子時，出現1至6的概率均為1/6；而使用作弊的色子時，6出現的概率為1/2，其他數字出現的概率均為1/10。同時，我們假定，在每次投擲中，從正常模式切換到作弊模式的概率為0.05，而由作弊模式切換到正常模式的概率為0.1，色子的切換是個馬爾科夫過程。如下圖所示。
　　
　　事實上，盡管我們非常關心賭場是否使用了作弊的色子，但是恐怕除了賭場的工作人員，沒有人知道真相，也就是說，內部狀態是隱藏的。我們可以做的，就是使用HMM建模，并運用相應算法計算出最可能的內部狀態序列。
　　
　　4．作弊賭場問題的BioJava實現
　　作為一個通用的基礎庫，BioJava提供了對HMM的支持。BioJava目錄下的demos/dp目錄下的Dice.java是色子問題的實現。
　　
　　該例子程序的核心是createCasino()方法。該方法中創建了實現描述問題的模型的馬爾科夫模型類的實例。
　　
　　　public static MarkovModel createCasino() {
　　 Symbol[] rolls=new Symbol[6];
　　
　　 //set up the dice alphabet
　　 SimpleAlphabet diceAlphabet=new SimpleAlphabet();
　　　　diceAlphabet.setName("DiceAlphabet");
　　
　　　　for(int i=1;i<7;i++) {
　　　　 try {
　　　　　　 rolls[i-1]= AlphabetManager.createSymbol((char)('0'+i),""+i,Annotation.EMPTY_ANNOTATION);
　　　　　　 diceAlphabet.addSymbol(rolls[i-1]);
　　　　　 } catch (Exception e) {
　　　　　　throw new NestedError(
　　　　　　　e, "Can't create symbols to rePResent dice rolls"
　　　　　　);
　　　　　}
　　　　}
　　
　　首先創建一個用于保存AlphabetManager產生的符號的符號數組，這些符號表示了投擲色子的結果。
　　
　　接下來，創建表示正常和作弊狀態下發散概率(emission probability)的分布以及狀態對象。
　　
　　　　int [] advance = { 1 };
　　　　Distribution fairD;
　　　　Distribution loadedD;
　　　　try {
　　　　　fairD = DistributionFactory.DEFAULT.createDistribution(diceAlphabet);
　　　　　loadedD = DistributionFactory.DEFAULT.createDistribution(diceAlphabet);
　　　　} catch (Exception e) {
　　　　　throw new NestedError(e, "Can't create distributions");
　　　　}
　　　　EmissionState fairS = new SimpleEmissionState("fair", Annotation.EMPTY_ANNOTATION, advance, fairD);
　　　　EmissionState loadedS = new SimpleEmissionState("loaded", Annotation.EMPTY_ANNOTATION, advance, loadedD);
　　
　　下面代碼創建HMM對象。
　　
　　　　SimpleMarkovModel casino = new SimpleMarkovModel(1, diceAlphabet, "Casino");
　　　　try {
　　　　　casino.addState(fairS);
　　　　　casino.addState(loadedS);
　　　　} catch (Exception e) {
　　　　　throw new NestedError(e, "Can't add states to model");
　　　　}
　　
　　下一步，我們還要設置狀態間的轉移規則
　　
　　　　//set up transitions between states.
　　　　try {
　　　　　casino.createTransition(casino.magicalState(),fairS);
　　　　　casino.createTransition(casino.magicalState(),loadedS);
　　　　　casino.createTransition(fairS,casino.magicalState());
　　　　　casino.createTransition(loadedS,casino.magicalState());
　　　　　casino.createTransition(fairS,loadedS);
　　　　　casino.createTransition(loadedS,fairS);
　　　　　casino.createTransition(fairS,fairS);
　　　　　casino.createTransition(loadedS,loadedS);
　　　　} catch (Exception e) {
　　　　　throw new NestedError(e, "Can't create transitions");
　　　　}
　　
　　當然，也還需要設置不同狀態下的發散概率.
　　
　　　　try {
　　　　　for(int i=0;i<rolls.length;i++) {
　　　　　　fairD.setWeight(rolls[i],1.0/6.0);
　　　　　　loadedD.setWeight(rolls[i], 0.1);
　　　　　}
　　　　　loadedD.setWeight(rolls[5],0.5);
　　　　} catch (Exception e) {
　　　　　throw new NestedError(e, "Can't set emission probabilities");
　　　　}
　　
　　最后，我們還需要初始化狀態轉移概率:
　　
　　　　//set up transition scores.
　　　　try {
　　　　　Distribution dist;
　　
　　　　　dist = casino.getWeights(casino.magicalState());
　　　　　dist.setWeight(fairS, 0.8);
　　　　　dist.setWeight(loadedS, 0.2);
　　
　　　　　dist = casino.getWeights(fairS);
　　　　　dist.setWeight(loadedS, 0.04);
　　　　　dist.setWeight(fairS, 0.95);
　　　　　dist.setWeight(casino.magicalState(), 0.01);
　　
　　　　　dist = casino.getWeights(loadedS);
　　　　　dist.setWeight(fairS, 0.09);
　　　　　dist.setWeight(loadedS, 0.90);
　　　　　dist.setWeight(casino.magicalState(), 0.01);
　　　　} catch (Exception e) {
　　　　　throw new NestedError(e, "Can't set transition probabilities");
　　　　}
　　
　　現在，描述本問題的HMM對象已經構造完成。
　　
　　在創建了HMM對象后，我們還需要實例化相應的動態規劃對象。
　　
　　 //create the Dynamic Programming matrix for the model.
　　 DP dp= DPFactory.DEFAULT.createDP(casino);
　　
　　現在, 至少我們有些東西可用了.用這個模型生成一條投擲骰子的序列:
　　
　　 StatePath obs_rolls = dp.generate(300);
　　　
　　 SymbolList roll_sequence = obs_rolls.symbolListForLabel(
　　　　 StatePath.SEQUENCE);
　　
　　generate()方法通過這個模型生成了300個標志的一個路徑,在第二行我們將這個路徑轉化成標志鏈.
　　
　　下面,我們用DP對象檢測一種DP算法.我們想處理我們剛生成的投擲序列標志鏈,使用Viterbi方法猜測每次投擲使用的是哪一個骰子。
　　
　　為了做到這點,我們創建包含投擲序列(roll_sequence)的多條標志鏈數組.我們使用單向隱馬模型.應用Viterbi算法,使用模型概率(ScoreType.PROBABILITY)來計算(你也可以使用空模型或者基于對數幾率的概率).這將生成模型支持的狀態路徑,這個路徑就是通過這個模型猜測的每次投擲使用哪個骰子的結果。
　　
　　 SymbolList[] res_array = {roll_sequence};
　　 StatePath v = dp.viterbi(res_array, ScoreType.PROBABILITY);
　　
　　輸出結果形如：
　　544552213525245666363632432522253566166546666666533666543261
　　fffffffffffflllllllllllfffffffffffflllllllllllllllllllffffff
　　ffffffffffffffffffffffffffffffffffllllllllllllllllllllffffff
　　
　　363546253252546524422555242223224344432423341365415551632161
　　ffffffffffffffffffffffffffffffffffffffffffffffffffffffffffff
　　ffffffffffffffffffffffffffffffffffffffffffffffffffffffffffff
　　
　　144212242323456563652263346116214136666156616666566421456123
　　fffffflllfffffffffffffffffffffffflfllllllllllllllllfffffffff
　　fffffffffffffffffffffffffffffffffffllllllllllllllllfffffffff
　　
　　34631354651433

上一篇：使用emacs-jde進行Java開發

下一篇：事件監聽器: 將JavaBeans接通起來的方法