致老婆的那些模型与算法之韩梅梅（转）-HMM， Hidden Markov Models，隐马模型
老婆我想一个人静静
静静？静静是谁！！！
韩梅梅，HMM， Hidden Markov Models，隐马模型。至于定义维基百科甚至百度皆可不再赘述。 现在哪里都能看见其踪影，部分因为韩梅梅真地太亲和，就像画图一样直观上让我们一窥序列本质。韩梅梅在生物里的用处之一，举个例子，给你一段序列：
CTTCATGTGAAAGCAGACGTAAGTCA如何推测可变剪切位点的位置？
在撕开韩梅梅的外衣前，我们谈谈大概怎么扯开？可变剪切位点是分开exon和intron的界限，而在该界限左右两边，exon和intron的属性本质上是不同的，最明显的不同是四种碱基的分配比例。所以找可变剪切位点的问题，就转化为怎样确定一个边界，能使得左右两边各自都非常非常像exon和intron，在各自碱基比例上都比较吻合exon和intron的属性。正经一点儿说，对于韩梅梅，exon和intron是两种状态（state)，需要做的是基于看得见的序列，逆向地去推断看不见的状态，这正是梅梅同学韩姓的来源（后续再扯）。
Multinomial Model生成DNA序列
在逆向推断之前，先正向来思考，基于已知碱基比例，如何生成序列。 假设这段序列四种碱基的概率分别是这样的：
PA=0.2PC=0.3PG=0.3PT=0.2
在R语言里，只需要这么做即可：
> nucleotides <- c("A", "C", "G", "T") # 四种碱基> probabilities1 <- c(0.2, 0.3, 0.3, 0.2) # 假设前提里的四种概率> seqlength <- 30 # 序列长度> sample(nucleotides, seqlength, rep=TRUE, prob=probabilities1) # 一句话生成序列
这是不是特别像玩幸运大转盘？
如果比例是这样的，
PA=0.1PC=0.41PG=0.39PC=0.1
那么还是类似的：
> probabilities2 <- c(0.1, 0.41, 0.39, 0.1) # 注意相对应的碱基> sample(nucleotides, seqlength, rep=TRUE, prob=probabilities2) # sample 用来抽样# rep 表示可放回# prob 则是我们的比例
到这里，老婆你懂得如何玩两个碱基比例不同的幸运大转盘了：
Markov Model生成DNA序列
上面这种固定的设置，自然里哪里有这么简单。人3G基因组，上1KB是一个碱基比例，下1KB也一定是这样么？明显不是，所以上述这种配置，图样图森破。
另一种理论是，在自然选择时，当前的碱基种类，是受上一个碱基影响的，比如连续出现A之后大概慢慢就要出现Terminator区域了。而这一个，正是一个不错的Markov Model。
于是，现在玩儿的就是四个转盘了。当前是A，与当前是G，随后的碱基比例是不一样的。
> nucleotides <- c("A", "C", "G", "T") > afterAprobs <- c(0.2, 0.3, 0.3, 0.2) > afterCprobs <- c(0.1, 0.41, 0.39, 0.1) > afterGprobs <- c(0.25, 0.25, 0.25, 0.25) > afterTprobs <- c(0.5, 0.17, 0.17, 0.17) > myTransitionMatrix <- matrix(c(afterAprobs, afterCprobs, afterGprobs, afterTprobs), 4, 4, byrow = TRUE) ＃ 构建好了Transition Matrix> rownames(myTransitionMatrix) <- nucleotides> colnames(myTransitionMatrix) <- nucleotides> myTransitionMatrix A C G TA 0.20 0.30 0.30 0.20C 0.10 0.41 0.39 0.10G 0.25 0.25 0.25 0.25T 0.50 0.17 0.17 0.17
该矩阵中，比如当前碱基是A，那转向到碱基C的概率是0.2； 若当前是C，转向到G的概率是0.41.借着这个韩梅梅矩阵，我们类似地，也来生成一条序列。
我们先写一个名字叫做generateMarkovSeq的函数。 在这个函数呢，myTransitionMatrix是核心，它决定了怎么转。另外，既然梅梅需要知道上一个碱基的状态，所以第一个碱基究竟是什么呢？所以我们定义一个initialProbs 来表示鸡生蛋蛋生鸡太初混沌时四种碱基的比例。一般地，我们不妨设定四者是相同的都为0.25。
> myInitialProbs <- c(0.25, 0.25, 0.25, 0.25)
现在看看具体的：
> generateMarkovSeq <- function(transitionMatrix, initialProbs, seqLength){ nucleotides <- c("A", "C", "G", "T") mysequence <- character() # 确定第一个碱基 # 正如玩一个幸运转盘的情况 firstnucleotide <- sample(nucleotides, 1, rep=TRUE, prob=initialProbs) mysequence[1] <- firstnucleotide # 开始继续生成序列 for (i in 2:seqLength){ prevnucleotide <- mysequence[i-1] # 在4个转盘里找到应该需要的那一个，比如上一个碱基是G，那么找到转盘G probabilities <- transitionMatrix[prevnucleotide,] # 回到玩一个转盘的情况 nucleotide <- sample(nucleotides, 1, rep=TRUE, prob=probabilities) mysequence[i] <- nucleotide } return(mysequence)}
这样一段基于梅梅生成的DNA序列就好了。
韩梅梅为什么姓韩
到此已经一窥Markov Model了，但Hidden Markov Model又从何谈起？
先前Markov Model里下一个碱基决定于上一个碱基的种类，因而有四种A、T、G、C的转盘。而现实中（或者更加接近现实情况），下一个碱基很大程度受到上一个碱基所处于的状态。每一个状态里，有各自的四种碱基比例。这种状态之间的转换，正是一个很好的HMM模型。
比如DNA序列里有两种状态，AT-rich和GC-rich。DNA序列就一直在这两种状态里切换来切换去。类似的，还可以在exon、剪切位点、intron三种状态里切换来切换去。因此一段DNA序列，包含着有两种信息：
第一种是你直观就看得见的，即碱基序列。比如第三碱基是T，紧随之后是A。 第二种则是你肉眼看不见的，即状态序列。第三个碱基是T，但它处于什么状态呢？是出自于AT－rich么？紧随其后的是A，它又处于什么状态呢？有没有可能是出自GC－rich呢？
那看不见的状态序列（State Path），梅梅同学的韩姓就来源于此。
切换和旋转转盘
现在正向体验一下HMM。假设我们已知AT－rich和GC－rich两种状态，两个转盘如下
：
每一次选定碱基之前所需要做的，首先确定用哪一种转盘，然后才能玩这个转盘。而第一步如何选转盘，这个问题之前你已经体验过了，就是如何利用Transition Matrix来玩那四个转盘。类似的，在AT－rich转盘和GC－rich转盘之间，存在一个相互切换的概率。如上图，我们假定已知，玩过AT－rich转盘之后，切换玩GC－rich转盘的概率是0.3继续玩儿AT－rich的概率则是1-0.3=0.7； 反过来则是0.1. 接下来很自然地，我们可以构建AT－rich与GC－rich转盘之间的切换概率。
> states <- c("AT-rich", "GC-rich") > ATrichprobs <- c(0.7, 0.3) > GCrichprobs <- c(0.1, 0.9) > thetransitionmatrix <- matrix(c(ATrichprobs, GCrichprobs), 2, 2, byrow = TRUE) > rownames(thetransitionmatrix) <- states> colnames(thetransitionmatrix) <- states> thetransitionmatrix AT-rich GC-richAT-rich 0.7 0.3GC-rich 0.1 0.9
Transition Matrix构建好了，类似地，针对每一种状态下，有它属于自己的一套碱基发生的概率。把所有状态的各自的概率表都构建起来，我们就构建了名叫Emission Matrix。表示在某个转盘下，转出各个碱基的概率。如AT－rich的转盘图，转出A的概率是0.39.
> nucleotides <- c("A", "C", "G", "T") > ATrichstateprobs <- c(0.39, 0.1, 0.1, 0.41) > GCrichstateprobs <- c(0.1, 0.41, 0.39, 0.1) > theemissionmatrix <- matrix(c(ATrichstateprobs, GCrichstateprobs), 2, 4, byrow = TRUE)> rownames(theemissionmatrix) <- states> colnames(theemissionmatrix) <- nucleotides> theemissionmatrix A C G TAT-rich 0.39 0.10 0.10 0.41GC-rich 0.10 0.41 0.39 0.10
韩梅梅生成DNA序列
> generatehmmseq <- function(transitionmatrix, emissionmatrix, initialprobs, seqlength) { nucleotides <- c("A", "C", "G", "T") states <- c("AT-rich", "GC-rich") # 记录着你看得见的碱基序列 mysequence <- character() # 记录着你看不到的状态序列 mystates <- character() # 太初之时，确定第一个碱基的状态 # 类似地，两种状态假定一开始都是0.5的概率 firststate <- sample(states, 1, rep=TRUE, prob=initialprobs) # 确定状态后，找到该状态转盘，即它的Emission Matrix probabilities <- emissionmatrix[firststate,] # 类似一开始的，玩转一个转盘你已经会了 firstnucleotide <- sample(nucleotides, 1, rep=TRUE, prob=probabilities) mysequence[1] <- firstnucleotide mystates[1] <- firststate for (i in 2:seqlength) { # 翻出上一个碱基的状态 prevstate <- mystates[i-1] # 找到该状态转盘的Transition Matrix stateprobs <- transitionmatrix[prevstate,] # 根据Transition Matrix，看是不是要切换转盘 # 类似的，玩转一个转盘，你早就会了 state <- sample(states, 1, rep=TRUE, prob=stateprobs) # 找到该状态下的Emission Matrix，开始生成序列 probabilities <- emissionmatrix[state,] nucleotide <- sample(nucleotides, 1, rep=TRUE, prob=probabilities) mysequence[i] <- nucleotide mystates[i] <- state # 随后就这样周而复始下去直到尽头 } # 我们打印出来看看生成的序列是什么效果 for (i in 1:length(mysequence)) { nucleotide <- mysequence[i] state <- mystates[i] print(paste("Position", i, ", State", state, ", Nucleotide = ", nucleotide)) } }
看看效果吧：
> theinitialprobs <- c(0.5, 0.5)> generatehmmseq(thetransitionmatrix, theemissionmatrix, theinitialprobs, 30)[1] "Position 1 , State AT-rich , Nucleotide = A"[1] "Position 2 , State AT-rich , Nucleotide = A"[1] "Position 3 , State AT-rich , Nucleotide = G"[1] "Position 4 , State AT-rich , Nucleotide = C"[1] "Position 5 , State AT-rich , Nucleotide = G"[1] "Position 6 , State AT-rich , Nucleotide = T"[1] "Position 7 , State GC-rich , Nucleotide = G"[1] "Position 8 , State GC-rich , Nucleotide = G"[1] "Position 9 , State GC-rich , Nucleotide = G"[1] "Position 10 , State GC-rich , Nucleotide = G"[1] "Position 11 , State GC-rich , Nucleotide = C"[1] "Position 12 , State GC-rich , Nucleotide = C"[1] "Position 13 , State GC-rich , Nucleotide = C"[1] "Position 14 , State GC-rich , Nucleotide = C"[1] "Position 15 , State GC-rich , Nucleotide = G"[1] "Position 16 , State GC-rich , Nucleotide = G"[1] "Position 17 , State GC-rich , Nucleotide = C"[1] "Position 18 , State GC-rich , Nucleotide = G"[1] "Position 19 , State GC-rich , Nucleotide = A"[1] "Position 20 , State GC-rich , Nucleotide = C"[1] "Position 21 , State GC-rich , Nucleotide = A"[1] "Position 22 , State AT-rich , Nucleotide = T"[1] "Position 23 , State GC-rich , Nucleotide = G"[1] "Position 24 , State GC-rich , Nucleotide = G"[1] "Position 25 , State GC-rich , Nucleotide = G"[1] "Position 26 , State GC-rich , Nucleotide = G"[1] "Position 27 , State GC-rich , Nucleotide = T"[1] "Position 28 , State GC-rich , Nucleotide = G"[1] "Position 29 , State GC-rich , Nucleotide = T"[1] "Position 30 , State GC-rich , Nucleotide = C"
从生成的结果看，GC－rich里的序列不都只是G和C，还有A和T。这正是因为在GC－rich的转盘里（即Emission Matrix）里A和T发生的概率不是0.
小结
至此，你已经正向地体验到了，HMM是怎么运转的。选定转盘方能玩转盘。选定转盘根据Transition Matrix来切换、玩转盘根据Emission Matrix来转出结果。
至于开篇的问题，如何根据看得见的碱基序列，去推测看不见的状态序列？未完待续。
Image credits
a Little Book of R for Bioinformatics参考
Hidden Markov Models, Chapter10, a Little Book of R for Bioinformatics.
Sean R Eddy. What is a hidden Markov model?, Nature Biotechonology, 2004
