序列比对(九)——从掷骰子谈HMM

从本文开始，到了序列比对的第二部分，主要是介绍HMM以及其在序列比对中的基础应用。

前面八篇文章介绍了动态规划在序列比对中的基础应用。从本文开始，开始介绍HMM（隐马尔可夫模型）。

Markov链

首先我们先介绍马尔可夫链（Markov链），这个大家可能都熟悉，高中数学中就学过。简单来说，Markov链就是一系列状态构成，每一个状态出现的概率只与它前一个状态有关。
具体地说，假设有一系列状态$\pi$,依次由${\pi }_1,{\pi }_2,...,{\pi }_L$构成，那么在${\pi }_1,{\pi }_2,...,{\pi }_{i-1}$已经发生的情况下，${\pi }_i$发生的概率只与${\pi }_{i-1}$相关。即：
$P({\pi }_i|{\pi }_1,{\pi }_2,...,{\pi }_{i-1})=P({\pi }_i|{\pi }_{i-1})$
由此，我们可以推导出：
$P(\pi )=P({\pi }_1,{\pi }_2,{\pi }_3,...{\pi }_L)$
$=P({\pi }_1)P({\pi }_2|{\pi }_1)P({\pi }_3|{\pi }_1,{\pi }_2)...P({\pi }_L|{\pi }_1,{\pi }_2,...,{\pi }_{L-1})$
$=P({\pi }_1)P({\pi }_2|{\pi }_1)P({\pi }_3|{\pi }_2)...P({\pi }_L|{\pi }_{L-1})$

为了建模的方便，一般还会假设一个初始状态${\pi }_0$，则上面的公式变成：
$P(\pi )=P({\pi }_0)P({\pi }_1|{\pi }_0)P({\pi }_2|{\pi }_1)P({\pi }_3|{\pi }_2)...P({\pi }_L|{\pi }_{L-1})$

为了更形象说明，以投骰子为例来说，假设我们有下面这张图：

图片引自《生物序列分析》

这张图是说有两种骰子，一种是“公平的”（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。

状态向量

上面的介绍中，使用骰子的过程就是一个简单的Markov链。令$\pi$为一系列掷骰子的状态，则 π i ∈ { F , L } \pi_i\in\{F, L\} πi​∈{F,L}。F为公平骰子，L为作弊骰子。Markov链中所有可能的状态构成的集合（向量），即此处的 { F , L } \{F, L\} {F,L}。称为状态向量。

符号向量

上图中，不同骰子状态下，掷出1-6的概率不一样，假设在第$i$次投掷中，使用的骰子状态是${\pi }_i$，而掷出的结果（符号）为$x_i$，那么 x i ∈ { 1 , 2 , 3 , 4 , 5 , 6 } x_i\in\{1,2,3,4,5,6\} xi​∈{1,2,3,4,5,6}。Markov链中所有可能的结果（符号）构成的集合（向量），即此处的 { 1 , 2 , 3 , 4 , 5 , 6 } \{1,2,3,4,5,6\} {1,2,3,4,5,6}，称为符号向量。

转移概率

为了方便，我们记$a_{kl}=P({\pi }_i=l|{\pi }_{i-1}=k)$，称为从状态$k$到状态$l$的转移概率。比如：
$a_{FF}=P({\pi }_i=F|{\pi }_{i-1}=F)=0.95;a_{FL}=P({\pi }_i=L|{\pi }_{i-1}=F)=0.05$
$a_{LL}=P({\pi }_i=L|{\pi }_{i-1}=L)=0.9;a_{LF}=P({\pi }_i=F|{\pi }_{i-1}=L)=0.1$

转移矩阵

不同状态之间互相转移的概率可以构成一个矩阵，成为转移矩阵。比如：

	F	L
F	0.95	0.05
L	0.1	0.9

初始向量

有了转移概率的记号，那么$P(\pi )$的计算公式可以写成：
$P(\pi )=P({\pi }_0)P({\pi }_1|{\pi }_0)P({\pi }_2|{\pi }_1)P({\pi }_3|{\pi }_2)...P({\pi }_L|{\pi }_{L-1})$
$=P({\pi }_0)a_{{\pi }_0{\pi }_1}{a}_{{\pi }_1{\pi }_2}{a}_{{\pi }_2{\pi }_3}...a_{{\pi }_{L-1}{\pi }_L}$
$=P({\pi }_0){\prod }_{i=0}^{L-1}{a}_{{\pi }_i{\pi }_{i+1}}$
其中，${\pi }_0$就是状态0，实际上它并不存在，只是为了建模的方便假设它存在。从状态0到其他状态的转移概率可以称为初始概率，所有可能的初始概率构成初始向量。比如：

	F	L
0	0.9	0.1

发射概率

不同骰子状态下某一符号出现的概率是不一样的，以骰子为例，公平骰子和作弊骰子掷出6的概率是不一样的，不光是6，掷出1-5的概率也是不一样的：
$P(x_i=k|{\pi }_i=F)=1/6,k=1,2,...,6$
$P(x_i=k|{\pi }_i=L)=0.1,k=1,2,...,5$
$P(x_i=6|{\pi }_i=L)=0.5$
我们将某一状态$l$下符号$k$出现的概率称为发射概率，记为$e_l(k)$，即
$e_l(k)=P(x_i=k|{\pi }_i=l)$

发射矩阵

不同状态下不同符号的发射概率可以构成一个矩阵，成为发射矩阵，如：

	1	2	3	4	5	6
F	1/6	1/6	1/6	1/6	1/6	1/6
L	0.1	0.1	0.1	0.1	0.1	0.5

隐马尔可夫模型

什么叫隐马尔可夫模型呢？还以上图为例。假设别人一共投掷了300次，那么我们可以看到一个长度为300的符号序列$x$，其中 x i ∈ { 1 , 2 , 3 , 4 , 5 , 6 } x_i\in\{1,2,3,4,5,6\} xi​∈{1,2,3,4,5,6}。但是我们并不知道每一次投掷用的是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");
}

（公众号：生信了）