维特比算法的实现与应用 HMM 隐马尔可夫链作业分析

最新推荐文章于 2021-09-30 17:49:50 发布

原创

最新推荐文章于 2021-09-30 17:49:50 发布 · 645 阅读

0 ·

CC 4.0 BY-SA版权

文章标签：

#python #机器学习

本文介绍了如何利用隐马尔可夫模型（HMM）和维特比算法来分析基因组序列，特别是针对GC组成的不均一性。问题设定为一个有两个状态（A和B）的HMM，分别对应不同的碱基组成。通过HMM的参数，如发射概率和转移概率，以及维特比算法，对给定的基因序列进行分割，确定属于状态A和状态B的区域。

问题提出：

大多数基因组序列的碱基组成不均一。特别是，GC组成可能会因地区而异。例如，在人类基因组中，存在“ CpG岛”，其显示出非常强的统计信号（在二核苷酸组成水平上甚至比单核苷酸组成水平更强），并且倾向于标记许多基因的5’末端。因此，自然出现将基因组序列客观地划分为不同组成区域的问题。分割基因组的一种方法是使用HMM算法。

假定可以将基因组建模为具有两个状态A和B的两个状态的HMM。假定此HMM的参数如下。状态A对于{A，C，G，T}的概率具有AT偏移的发射分布{0.35，0.15，0.15，0.35}。状态B具有GC偏向的发射分布{0.15，0.35，0.35，0.15}。状态A以0.001的概率切换到状态B。状态B以0.01的概率切换回状态A。假设HMM的初始分布是均匀的； 0.5，两个状态为0.5。

输入：

首先读模型结构 输入的格式是这样的
2 4 ACGT
.5 .5
.999 .001 .35 .15 .15 .35
.01 .99 .15 .35 .35 .15
给出基因序列判断哪些是状态a , 哪些是状态b.

思路


"_@AUTHOR_ Lavinia ZHANGMENGFAN"
"STUDENT ID 55955735"

import numpy as np

"""
get the parameters we need about markov model from the .hmm file we have 
"""
def get_hmm_model(model_path="example.hmm"):
    f=open(model_path)
    hmm_model=f.readlines()         #reading files
    count=1
    matrix=[]
    for line in hmm_model:      
        parameters=line.split()    
        
        if count == 1 :
            state_size=int(parameters[0])
            symbols_size=int(parameters[1])
            symbols=[]
            for i in parameters[2]:
                symbols.append(i)

        elif count == 2 :
            primary_p=np.array(list(map(float,parameters)))

        else:
            matrix.append(np.array(list(map(float,line.split()))))
        count=count+1
    matrix=np.array(matrix)
    transition_p=matrix[:,:2]      
    emission_p=matrix[:,2:]
    f=open("hmm_optimal_path.txt",'a')
    f.writelines(['The markov model have ',str(state_size), ' states ：state A and state B respectively.'])
    f

最低0.47元/天解锁文章