Python一阶矩马尔可夫过程形成任意DNA序列完成实例

Python实现一阶矩马尔可夫过程生成DNA序列

小编小本本今天要给大家介绍的是Python完成一阶矩马尔可夫过程形成任意DNA序列实例详细说明。一阶矩马尔可夫过程是指当前状态仅与上一状态有关，因此对DNA序列而言，一阶矩马尔可夫过程可以看作现阶段碱基对的种类仅在于上一位碱基对种类。

基本原理

以一个编码序列开始，它的第一位有可能是A、T、G、C四种碱基对（且概率同等，均为0.25），如果此处是A，那么下一个碱基对可选A、T、G、C的概率为分别为0.25、0.20、0.20、0.20，如果此处没有碱基对了（即编码序列完毕），那么下一个状态（终点）的概率为0.15。

下面的图1显示了DNA序列的一阶马尔科夫链。

图1  DNA序列的一阶马尔科夫链

代码实现

以下代码在Jupyter Notebook (Python 3.7)中运行，功能是随机生成一定数量的DNA序列，记录序列长度并绘制长度分布图。

import numpy
import random
import seaborn as sns
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.append(START)#将起始碱基添加至Seq中
        while START!="E":
            if START=="A":
                change=numpy.random.choice(transitionName[0],p=transitionMatrix[0])
                #以transitionMatrix矩阵第一行的概率分布随机抽取transitionName第一行包含的事件
            if change=="AA":
                START="A"
                #如果转移状态是AA(即A碱基接下来的碱基是A，则将起始碱基设为A)
            elif change=="AG":
                START="G"
            elif change=="AC":
                START="C"
            elif change=="AT":
                START="T"
            elif change=="AE":
                START="E"
            elif START=="G":
                change=numpy.random.choice(transitionName[1],p=transitionMatrix[1])
            if change=="GA":
                START="A"
            elif change=="GG":
                START="G"
            elif change=="GC":
                START="C"
            elif change=="GT":
                START="T"
            elif change=="GE":
                START="E"
            elif START=="C":
                change=numpy.random.choice(transitionName[2],p=transitionMatrix[2])
            if change=="CA":
                START="A"
            elif change=="CG":
                START="G"
            elif change=="CC":
                START="C"
            elif change=="CT":
                START="T"
            elif change=="CE":
                START="E"
            elif START=="T":
                change=numpy.random.choice(transitionName[3],p=transitionMatrix[3])
            if change=="TA":
                START="A"
            elif change=="TG":
                START="G"
            elif change=="TC":
                START="C"
            elif change=="TT":
                START="T"
            elif change=="TE":
                START="E"
            if START!="E":
                Seq.append(START)#如果状态转移后不为End(E)，则将转移后的碱基加到Seq序列中
        i+=1
        Len.append(len(Seq))
        if len(Seq)>max_len:
            max_len=len(Seq)
        #print(''.join(Seq))
        Seq.clear()
    plt.hist(numpy.array(Len),bins=max_len,edgecolor="white")
    #显示横轴标签
    plt.xlabel("DNA Sequence Length")
    #显示纵轴标签
    plt.ylabel("Frequency")
    #显示图标题
    plt.title("Histogram of frequency distribution of DNA sequence length")
    plt.show()
    print("DNA序列的最大长度为:",max_len)
    print("DNA序列长度的众数为:",max(Len,key=Len.count))
%matplotlib noteBook#若未使用Jupyter Notebook，此句不需要
RandomDNAs(1000)#1000表示随机生成1000条序列

运行结果

以下4个序列长度分布统计展示了随着随机生成的序列数量增多，DNA序列长度分布逐渐集中，且长度为1bp的序