用BioJava开发基于隐马尔科夫模型的应用

核心提示：1． 隐马尔科夫模型在生物信息学中的应用作为一种强大的概率模型，隐马尔科夫模型被广泛地应用到语音识别、计算机视觉等领域，用BioJava开发基于隐马尔科夫模型的应用，而在计算分子生物学和生物信息学领域，隐马尔科夫模型同样有着广泛的应用，我们还需要实例化相应的动态规划对象， //create the Dynamic Pr

　　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