1. 原理
对于DNA序列,一阶马尔科夫链可以理解为当前碱基的类型仅取决于上一位碱基类型。如图1所示,一条序列的开端(由B开始)可能是A、T、G、C四种碱基(且可能性相同,均为0.25),若序列的某一位是A,则下一位碱基是A、T、G、C的概率分别为0.25、0.20、0.20、0.20,下一位无碱基(即序列结束,状态为E)的概率为0.15。
2. 代码实现
以下代码运行于Jupyter Notebook (Python 3.7);代码功能是随机生成一定数量的DNA序列,统计序列长度并绘制分布图。若希望显示随机生成的序列,将代码# print(''.join(Seq))前的#删除即可。
import numpy
import random
import matplotlib.pyplot as plt
# 状态空间
states = ["A","G","C","T","E"]
# 可能的事件序列
transitionName = [["AA","AG","AC","AT","AE"],
["GA","GG","GC","GT","GE"],
["CA","CG","CC","CT","CE"],
["TA","TG","TC","TT","TE"],]
# 概率矩阵(转移矩阵)
transitionMatrix = [[0.25,0.20,0.20,0.20,0.15],
[0.20,0.25,0.20,0.20,0.15],
[0.20,0.20,0.25,0.20,0.15],
[0.20,0.20,0.20,0.25,0.15]]
def RandomDNAs(Num):
max_len = 0
i = 0
Seq = [] #创建列表(Seq)用于添加碱基,以组成DNA序列
Len = [] #创建列表(Len)用于记录每条生成序列的长度
while i != Num:
Base = ["A","G","C","T"]
START = random.choice(Base) #随机从碱基中选择一个作为序列的起始碱基
Seq.a