Alex Dou
autoasm@yahoo.com
[摘要] 隐马尔科夫模型及其变体在生物信息学领域有着广泛的应用，本文简要地介绍了隐马尔科夫模型，并讲解如何用BioJava实现一个具体的模型。
[关键词] 隐马尔科夫模型 色子问题 BioJava
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，色子的切换是个马尔科夫过程。如下图所示。
图3.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
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
346313546514332164351242356166641344615135266642261112465663
ffffffffffffffffffffffffffffffffffffffffffffffffffffffffffff
ffffffffffffffffffffffffffffffffffffffffffffffffffffffffffff
最上面的一行是隐马模型生成的300个投掷结果(数字代表点数).第二行是投掷时使用的状态(f代表正常的骰子,l代表作弊的骰子,我们用它们来创建SimpleEmissionState对象,这个对象代表骰子)。
最后一行和上一行很相似,但是却是从Viterbi算法的结果??状态路径v生成的。本次的结果是相当不错的，但是有些情况下，Viterbi算法的结果可能与实际情况有较大出入。
5． 结束语
本文所介绍的色子问题仅仅是一个简化的模型，在实际应用中，要获得模型中的发散概率和状态转移概率往往并不容易（以色子作弊问题为例，实际上我们很难准确地知道作投掷作了手脚的色子时，6出现的概率），需要对模型进行训练，以提高预测精度。
本文仅仅介绍了BioJava实现HMM的一个简单例子，实际上，HMM还提供了对Profile HMM的支持，使用BioJava的这些特性，可以开发出产品级的生物信息学应用系统。
参考文献
1. 基因组
2. Introduction to Bioinformatics
3. 生物序列分析
4. 隐马尔科夫模型在语音识别中的应用
5. Simple HMM with BioJava
原文连接>> http://autoasm.blog-city.com/read/615665.htm