目录
题目描述
0.05X低覆盖度测序的情况下每条染色体只能测到一些片段，所以对于杂合个体，测到的可能是亲本之一，所以无法确定其基因型；第二也存在测序错误。
参数设置
通过隐马尔可夫模型可以纠正测序错误，HMM模型包括五大参数：初始状态集；初始概率矩阵；观测集；转移概率矩阵；输出概率矩阵。其中转移概率矩阵需要自己设置一下。Parent-independent genotyping for constructing an ultrahigh-density linkage map based on population sequencing(Xie Weibo,2010) 这篇文献中将转移为除自身之外的状态均设为R，使得某一个状态的总转移概率为(1-R)+R+R=1+R>1，这是有一点问题的。实际一个状态转移为其它状态的概率应该看初始分布的比例，这个题目中，MZ:ZZ=1：2，则MM转移为MZ的概率应为R/3，转移为ZZ的概率为2R/3。（重组率为R，则非重组即仍转移为自身的概率为1-R）
输出概率矩阵比较容易得到，例如MM输出为m和z的概率分别为1-E和E（E为错误率）。其它三个参数题目中交代的比较清楚。
下面是具体代码实现。
Python代码
import numpy as np
#定义viterbi函数
def viterbi(obs_seq, state_names, initial_dist, trans_prob_matrix, emit_prob_matrix, ):
num_states = len(initial_dist)
# 初始化概率矩阵和回溯矩阵
probabilities = np.zeros((len(obs_seq), num_states))
backtrack = np.zeros((len(obs_seq), num_states), dtype=int)
# 初始化第一行的概率，基于初始分布和发射概率
probabilities[0, :] = initial_dist * emit_prob_matrix[:, int(obs_seq[0])] #第一个观测值为1，即要乘以输出概率的第二列（m）
# 前向传播
for i in range(1, len(obs_seq)):
for current_state in range(num_states):
# 计算当前观察值下每个状态的概率，基于上一个状态的情况
transition_probabilities = probabilities[i - 1, :] * trans_prob_matrix[:, current_state]
# 更新当前观察值下每个状态的概率，选择转移概率最大的，再计算当前情况下输出为某一观测值的概率
probabilities[i, current_state] = np.max(transition_probabilities) * emit_prob_matrix[current_state, int(obs_seq[i])]
# 记录回溯路径：记录最大概率的上一个状态的索引
backtrack[i, current_state] = np.argmax(transition_probabilities)
# 回溯以找到最可能的路径
path = [np.argmax(probabilities[-1])] # 从最后一个观察值中具有最大概率的状态开始
#print("Probabilities matrix is :\n",probabilities) #打印观测序列概率
#print(backtrack) #打印回溯路径
#将当前观察值下最可能的状态的索引（即前一个时间步的状态）加入到路径中，用于后续回溯找到整个序列的最可能路径。
for i in range(len(obs_seq) - 2, -1, -1):
path.append(backtrack[i + 1, path[-1]])
# 将索引转换为状态名称
path = [state_names[state_idx] for state_idx in path[::-1]]
return ''.join(path) #将列表 path 中的字符串元素连接成一个单一的字符串。
#定义五大参数
#参数一：观测值
obs_sequence = '1111111011111101111111000000001000000001000001000000000000001011001010101010111111101111111111101111111011111'
#参数二：初始状态分布,1为MM,2为MZ,0为ZZ
state_names = ['1', '2', '0']
#参数三：初始概率矩阵，分别对应MM,MZ,ZZ
initial_distribution = np.array([0.4, 0.2, 0.4])
#参数四：转移概率矩阵，行列分别为MM,MZ,ZZ
transition_matrix = np.array([[0.99, 0.01/3, 0.02/3],
[0.01/2, 0.99, 0.01/2],
[0.02/3, 0.01/3, 0.99]])
#参数五：输出概率矩阵，行为MM,MZ,ZZ;列为z,m
emission_matrix = np.array([[0.03, 0.97],
[0.5, 0.5],
[0.97, 0.03]])
result = viterbi(obs_sequence, state_names, initial_distribution, transition_matrix, emission_matrix, )
print('函数的输出结果:')
print(result)
# 参考结果
reference_answer = '1111111111111111111111000000000000000000000000000000000000002222222222222222111111111111111111111111111111111'
print('参考结果:')
print(reference_answer)运行结果如下：
R代码实现
# 定义viterbi函数
viterbi <- function(obs_num, state_names, initial_dist, trans_prob_matrix, emit_prob_matrix, symbols) {
# 获取状态数和观测数
num_states <- length(initial_dist)
n_obs <- length(obs_num)
# 初始化概率矩阵和回溯矩阵
probabilities <- matrix(0, nrow = num_states, ncol = n_obs)
backtrack <- matrix(0, nrow = num_states, ncol = n_obs)
# 初始化第一行的概率，基于初始分布和输出概率
probabilities[, 1] <- initial_dist * emit_prob_matrix[, which(symbols == obs_num[1])]
backtrack[, 1] <- 1:num_states
# Viterbi:向前传播
for (t in 2:n_obs) {
for (s in 1:num_states) {
# 计算最大概率和对应的状态
max_prob <- 0
max_state <- 0
for (s_prev in 1:num_states) {
prob <- probabilities[s_prev, t - 1] * trans_prob_matrix[s_prev, s] * emit_prob_matrix[s, which(symbols == obs_num[t])]
if (prob > max_prob) {
max_prob <- prob
max_state <- s_prev
}
}
# 更新动态规划表和回溯路径
probabilities[s, t] <- max_prob
backtrack[s, t] <- max_state
}
}
# 回溯最优路径
path <- numeric(n_obs)
# 从最后一列开始
path[n_obs] <- which.max(probabilities[, n_obs])
# 逆序回溯
for (t in (n_obs - 1):1) {
path[t] <- backtrack[path[t + 1], t + 1]
}
# 转换为数字格式
path <- as.character(path)
path[path == "1"] <- "1" # 路径位置为1的是MM，还是1不变
path[path == "3"] <- "0" # 路径位置为3的是ZZ，转换为0
path[path == "2"] <- "2" # 路径位置为2的是MZ，还是2不变
# 返回最优路径
return(path)
}
#定义五大参数
#参数一：观测值
obs_seq <- "1111111011111101111111000000001000000001000001000000000000001011001010101010111111101111111111101111111011111"
obs_num <- as.numeric(strsplit(obs_seq, "")[[1]])# 转换为数值向量
obs_num <- obs_num
#参数二：初始状态分布,1为MM,2为MZ,0为ZZ
state_names <- c("MM", "MZ", "ZZ") # 隐状态，即真实的基因型
symbols <- c("1", "0") # 观测符号，即测序获得的基因型
#参数三：初始概率矩阵，分别对应MM,MZ,ZZ
initial_dist <- c(0.4, 0.2, 0.4) # 初始概率，即后代中基因型为MM、MZ和ZZ的比例
#参数四：转移概率矩阵，行列分别为MM,MZ,ZZ
trans_prob_matrix <- matrix(c(0.99, 0.01/3, 0.02/3, # 转移概率矩阵，即相邻位点的基因型发生重组的概率
0.01/2, 0.99, 0.01/2,
0.02/3, 0.01/3, 0.99), nrow = 3, byrow = TRUE)
#参数五：输出概率矩阵，行为MM,MZ,ZZ;列为m,z
emit_prob_matrix <- matrix(c(0.97, 0.03, # 输出概率矩阵，即真实的基因型输出为测序结果的概率
0.5, 0.5,
0.03, 0.97), nrow = 3, byrow = TRUE)
# 调用viterbi函数
result <- viterbi(obs_num, state_names, initial_dist, trans_prob_matrix, emit_prob_matrix,symbols)
result_string <- paste(result, collapse = "") #使用paste连接成单一字符串
print("函数的输出结果:")
print(result_string)
reference_answer <- c("1111111111111111111111000000000000000000000000000000000000002222222222222222111111111111111111111111111111111")
print("参考结果:")
print(reference_answer)运行结果如下：
总结
隐马尔可夫模型（Hidden Markov Model, HMM）可用于解决多种问题，包括评估问题，给定模型，求某个观察值序列的概率，可用向前算法实现；解码问题，对给定模型和观察值序列，求可能性最大的状态序列，可用Viterbi算法实现；学习问题，对给定的一个观察值序列，调整参数，使得观察值出现概率最大，可用向前向后算法（EM）实现。
此外，HMM有很多应用：语音识别：声音信号可以被视为一个时间序列，而HMM可以用于建模声音信号中的语音单元，如音素。通过学习发射概率矩阵，HMM可以用于识别说话人所说的单词或短语。手写体识别：HMM被用于建模手写字母或数字的轨迹。每个字母或数字可以被视为一个隐藏状态序列，而观察序列是由手写轨迹产生的数据点。生物信息学：HMM被广泛应用于生物信息学中的序列分析，如基因识别、蛋白质结构预测等。DNA或蛋白质序列可以被视为隐藏的状态序列，而实验数据（例如DNA片段或蛋白质序列）则是观察序列。自然语言处理：HMM被用于词性标注、分词、命名实体识别等任务。文本数据中的单词或词组可以被视为隐藏的状态序列，而观察序列是文本中的单词。金融领域：HMM被用于建模金融时间序列，如股票价格变化。在这个应用中，隐藏状态可以表示市场的不同状态，而观察序列是金融市场的价格变动。天气预测：其中隐藏状态表示天气的不同状态（晴天、雨天、多云等），而观察序列是气象观测数据。
扫一扫
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
