从本文开始,到了序列比对的第二部分,主要是介绍HMM以及其在序列比对中的基础应用。
前面八篇文章介绍了动态规划在序列比对中的基础应用。从本文开始,开始介绍HMM(隐马尔可夫模型)。为什么要介绍它呢?因为该模型在发现Motif、预测CpG岛等多方面有广泛的应用。而这些又和序列比对息息相关。所以我们要了解该模型。不过,为了方便,这一部分的开头几篇文章都会以掷骰子为例来对HMM展开讨论。
(说明:由于公众号对Latex的公式支持不好,所以涉及到Latex公式的部分暂时先以图片显示)
首先我们先介绍马尔可夫链(Markov链),这个大家可能都熟悉,高中数学中就学过。简单来说,Markov链就是一系列状态构成,每一个状态出现的概率只与它前一个状态有关。
为了更形象说明,以投骰子为例来说,假设我们有下面这张图:
图片引自《生物序列分析》
这张图是说有两种骰子,一种是“公平的”(Fair,简写F),一种是“作弊的”(Loaded,简写L)。F骰子掷出1-6的概率是一样的,而L骰子掷出6的概率为0.5,其余1-5的概率都是0.1。 此外,如果这次使用F骰子,那么下次仍然使用F骰子的概率是0.95,换用L骰子的概率是0.05。相反,如果这次使用L骰子,那么下次仍然使用L骰子的概率是0.9,换用F骰子的概率是0.1。
为了后续深入学习隐马模型,我们首先得写一个程序,能根据转移矩阵以及发射矩阵生成一个随机的状态序列以及相应的符号序列。
效果如下:一共投掷了300次,Rolls代表掷出来的数字,Die代表投掷时使用的是公平骰子(F)还是作弊骰子(L)。
代码如下:
#include <stdio.h>
#include <stdlib.h>
#include <time.h>
typedef char State;
typedef char Result;
State state[] = {'F', 'L'}; // 所有的可能状态
Result result[] = {'1', '2', '3', '4', '5', '6'}; // 所有的可能符号
double init[] = {0.9, 0.1}; // 初始状态的概率向量
double emission[][6] = { // 发射矩阵:行对应着状态,列对应着符号
1.0/6, 1.0/6, 1.0/6, 1.0/6, 1.0/6, 1.0/6,
0.1, 0.1, 0.1, 0.1, 0.1, 0.5
};
double trans[][2] = { // 转移矩阵:行和列都是状态
0.95, 0.05,
0.1, 0.9
};
const int nstate = 2;
const int nresult = 6;
State* rst; // 一串随机状态序列
Result* rres; // 一串随机符号序列
int random(double* prob, const int n);
void randSeq(State* st, Result* res, const int n);
void printState(State* st, const int n);
void printResult(Result* res, const int n);
int main(void) {
int i;
int n = 300;
int ll = 60;
int nl, nd;
if ((rst = (State*) malloc(sizeof(State) * (n))) == NULL || \
(rres = (Result*) malloc(sizeof(Result) * (n))) == NULL) {
fputs("Error: out of space!\n", stderr);
exit(1);
}
randSeq(rst, rres, n);
nl = n / ll;
nd = n % ll;
for (i = 0; i < nl; i++) {
printf("Rolls\t");
printResult(rres + ll * i, ll);
printf("Die\t");
printState(rst + ll * i, ll);
printf("\n");
}
if (nd > 0) {
printf("Rolls\t");
printResult(rres + ll * i, nd);
printf("Die\t");
printState(rst + ll * i, nd);
printf("\n");
}
free(rst);
free(rres);
}
// 根据一个概率向量从0到n-1随机抽取一个数
int random(double* prob, const int n) {
int i;
double p = rand() / 1.0 / (RAND_MAX + 1);
for (i = 0; i < n - 1; i++) {
if (p <= prob[i])
break;
p -= prob[i];
}
return i;
}
// 根据转移矩阵和发射矩阵生成一串随机状态和符号
void randSeq(State* st, Result* res, const int n) {
int i, ls, lr;
srand((unsigned int) time(NULL));
ls = random(init, nstate);
lr = random(emission[ls], nresult);
st[0] = state[ls];
res[0] = result[lr];
for (i = 1; i < n; i++) {
ls = random(trans[ls], nstate);
lr = random(emission[ls], nresult);
st[i] = state[ls];
res[i] = result[lr];
}
}
void printState(State* st, const int n) {
int i;
for (i = 0; i < n; i++)
printf("%c", st[i]);
printf("\n");
}
void printResult(Result* res, const int n) {
int i;
for (i = 0; i < n; i++)
printf("%c", res[i]);
printf("\n");
}