|
我如何创建一个隐马模型序列谱(Profile HMM)?
隐马模型序列谱(例如HMMER程序使用的)是一种敏感度非常高的搜索模体的方法.通常使用Baum-Welch算法对一组序列进行训练,这组序列含有我们感兴趣模体.
这个算法能将模型参数不断优化直到符合终止标准.一旦序列谱被构建,则能用Viterbi算法决定最可能生成观测序列的状态序列.如果在观测序列中有足够多的
状态匹配或者使用一个定义了阈值的打分矩阵(比如说对数几率)且分值超过阈值则该序列被认为含有模体.下面的程序展示了这点.
首先创建一个隐马模型序列谱.
// 创建一个十二列的DNA隐马模型序列谱,用默认的分布工厂类创建转移和生成分布
ProfileHMM hmm = new ProfileHMM(DNATools.getDNA(),
12,
DistributionFactory.DEFAULT,
DistributionFactory.DEFAULT,
"my profilehmm");
// 为模型创建一个动态规划矩阵
dp = DPFactory.DEFAULT.createDP(hmm);
这时你要读取训练序列集.
// 存储训练集的数据库
SequenceDB db = new HashSequenceDB();
// 在这里编码,装载训练集
现在,以一致的数值初始化所有的模型参数。另外，参数可以随机化或者猜测最符合模型的参数值。
使用Baum-Welch算法优化序列。
// 将所有参数以同一数值训练
ModelTrainer mt = new SimpleModelTrainer();
// 注册模型
mt.registerModel(hmm);
// 使用空权重训练参数
mt.setNullModelWeight(1.0);
mt.train();
// 使用BW训练器，以动态规划矩阵为参数
BaumWelchTrainer bwt = new BaumWelchTrainer(dp);
// 使用匿名内类制定终止标准：20遍
StoppingCriteria stopper = new StoppingCriteria() {
public boolean isTrainingComplete(TrainingAlgorithm ta) {
return (ta.getCycle()>20);
}
};
// 使用空模型（null model)权重1.0和上面定义的终止标准，优化动态规划矩阵
bwt.train(db,1.0,stopper);
下面的例子对一条序列打分，并且输出状态序列
SymbolList test = null;
// 在这里编码，并且初始化待测序列
// 将待测序列装进一个数组中，这里使用数组是因为使用HMM的双序列比对需要在数组中的两条标志链。
SymbolList[] sla = {test};
// 解码最可能的状态序列，产生几率分值
StatePath path = dp.viterbi(sla,ScoreType.ODDS);
System.out.println("Log Odds = "+path.getScore());
// 输出状态序列
for(int i = 1; i<= path.length();i++){
System.out.println(path.symbolAt(StatePath.STATES, i).getName());
}
--BACK TO TOP
|
