注册 登录  
 加关注
   显示下一条  |  关闭
温馨提示!由于新浪微博认证机制调整,您的新浪微博帐号绑定已过期,请重新绑定!立即重新绑定新浪微博》  |  关闭

Puriney's Notes

Puriney=purine+Y, my Wonderland

 
 
 

日志

 
 

致老婆的那些模型与算法之韩梅梅 (1)  

2013-11-03 10:28:20|  分类: 默认分类 |  标签: |举报 |字号 订阅

  下载LOFTER 我的照片书  |

老婆我想一个人静静

静静?静静是谁!!!

韩梅梅,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) # 一句话生成序列

这是不是特别像玩幸运大转盘?

致老婆的那些模型与算法之韩梅梅 (1) - Puriney - Purineys Notes

如果比例是这样的,

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 则是我们的比例

到这里,老婆你懂得如何玩两个碱基比例不同的幸运大转盘了:

 致老婆的那些模型与算法之韩梅梅 (1) - Puriney - Purineys Notes

Markov Model生成DNA序列

上面这种固定的设置,自然里哪里有这么简单。人3G基因组,上1KB是一个碱基比例,下1KB也一定是这样么?明显不是,所以上述这种配置,图样图森破。

另一种理论是,在自然选择时,当前的碱基种类,是受上一个碱基影响的,比如连续出现A之后大概慢慢就要出现Terminator区域了。而这一个,正是一个不错的Markov Model。

于是,现在玩儿的就是四个转盘了。当前是A,与当前是G,随后的碱基比例是不一样的。

致老婆的那些模型与算法之韩梅梅 (1) - Puriney - Purineys Notes

> 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    T
A 0.20 0.30 0.30 0.20
C 0.10 0.41 0.39 0.10
G 0.25 0.25 0.25 0.25
T 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两种状态,两个转盘如下

: 致老婆的那些模型与算法之韩梅梅 (1) - Puriney - Purineys Notes

每一次选定碱基之前所需要做的,首先确定用哪一种转盘,然后才能玩这个转盘。而第一步如何选转盘,这个问题之前你已经体验过了,就是如何利用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-rich
AT-rich     0.7     0.3
GC-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    T
AT-rich 0.39 0.10 0.10 0.41
GC-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

参考

  1. Hidden Markov Models, Chapter10, a Little Book of R for Bioinformatics.
  2. Sean R Eddy. What is a hidden Markov model?, Nature Biotechonology, 2004
  评论这张
 
阅读(837)| 评论(0)
推荐 转载

历史上的今天

在LOFTER的更多文章

评论

<#--最新日志,群博日志--> <#--推荐日志--> <#--引用记录--> <#--博主推荐--> <#--随机阅读--> <#--首页推荐--> <#--历史上的今天--> <#--被推荐日志--> <#--上一篇,下一篇--> <#-- 热度 --> <#-- 网易新闻广告 --> <#--右边模块结构--> <#--评论模块结构--> <#--引用模块结构--> <#--博主发起的投票-->
 
 
 
 
 
 
 
 
 
 
 
 
 
 

页脚

网易公司版权所有 ©1997-2017