Python一阶马尔科夫链生成随机DNA序列实现示例

目录
  • 1. 原理
  • 2. 代码实现
  • 3. 运行结果

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。

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

2. 代码实现

以下代码运行于Jupyter Notebook (Python 3.7);代码功能是随机生成一定数量的DNA序列,统计序列长度并绘制分布图。若希望显示随机生成的序列,将代码# print(''.join(Seq))前的#删除即可。

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条序列

3. 运行结果

从以下4个序列长度分布统计可以看到,随着随机生成的序列数量增多,序列长度分布愈发集中,且长度为1bp的序列占比最多且逐渐增加。

图2 10,000条DNA序列的序列长度分布统计

10,000条DNA序列的序列中

DNA序列的最大长度为: 65

DNA序列长度的众数为: 1

图3 100,000条DNA序列的序列长度分布统计

100,000条DNA序列的序列中

DNA序列的最大长度为: 71

DNA序列长度的众数为: 1

以上就是Python实现一阶马尔科夫链生成随机DNA序列的详细内容,更多关于Python一阶马尔科夫DNA序列的资料请关注我们其它相关文章!

(0)

相关推荐

  • Python办公自动化批量处理文件实现示例

    目录 引言 需求分析 Python实现 结束语 引言 要说在工作中最让人头疼的就是用同样的方式处理一堆文件夹中文件,这并不难,但就是繁.所以在遇到机械式的操作时一定要记得使用Python来合理偷懒!今天我将以处理微博热搜数据来示例如何使用Python批量处理文件夹中的文件,主要将涉及: Python批量读取不同文件夹() Pandas数据处理() Python操作Markdown文件() 需求分析 首先来说明一下需要完成的任务,下面是我们的文件夹结构 因为微博历史热搜是没有办法去爬的,所以只能写

  • python从gbff文件中直接提取cds序列

    目录 什么是GBFF文件 每个序列条目所代表的意义 最后直接上代,更改输入和输出文件即可使用 什么是GBFF文件 GenBank纯文本文件格式(GenBank flatfile, 简称GBFF) GBFF是GenBank数据库的基本信息单位 GBFF序列文件由单个的序列条目组成. 序列条目由字段组成,每个字段由关键字起始,后面为该 字段的具体说明. 字段分若干次子字段,以次关键字或特性表说明符开始. 每个序列条目以双斜杠“//*作结束标记 每个序列条目所代表的意义 1. LOCUS(代码)序列的

  • Python办公自动化Word转Excel文件批量处理

    目录 前言 首先使用Python将Word文件导入 row和cell解析所需内容 内层解析循环 前言 大家好,今天有一个公务员的小伙伴委托我给他帮个忙,大概是有这样一份Word(由于涉及文件私密所以文中的具体内容已做修改) 一共有近2600条类似格式的表格细栏,每个栏目包括的信息有: 日期 发文单位 文号 标题 签收栏 需要提取其中加粗的这三项内容到Excel表格中存储,表格样式如下: 也就是需要将收文时间.文件标题.文号填到指定位置,同时需要将时间修改为标准格式,如果是完全手动复制和修改时间,

  • python中序列的逆序方式

    目录 序列的逆序方式 1. range 函数 2. reversed 函数 3. 其他方法 一个字符串的逆序函数 序列的逆序方式 1. range 函数 一般 for 循环中总会用到 range 函数来进行顺序遍历,同样的,range 也能表示序列的逆序. 在 range(start, end, step) 中,start 表示序列的起始索引(默认为0),end 表示终止索引,step 表示移动步长(默认为1).由于 range 函数是“顾头不顾尾” 的形式,因此实际上其遍历的索引终止位置是 e

  • 自学python求已知DNA模板的互补DNA序列

    目录 DNA序列 简述其代码 原始序列上进行替换 利用upper()输出大写结果 结尾 DNA序列 ACTGATCGATTACGTATAGTATTTGCTATCATACATATATATCGATGCGTTCAT 求其互补DNA序列. 在生物上DNA互补序列简述表达可以表示为:A与T,C与G互补,可以理解为将上述序列中现有的A用T代替,C用G代替,T用A代替,G用C代替,则其互补序列为: TGACTAGCTAATGCATATCATAAACGATAGTATGTATATATAGCTACGCAAGTA 根

  • Python一阶马尔科夫链生成随机DNA序列实现示例

    目录 1. 原理 2. 代码实现 3. 运行结果 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. 图1 DNA序列的一阶马尔科夫链 2. 代码实现 以下代码运行于Jupyter Notebook (Pyt

  • python实现隐马尔科夫模型HMM

    一份完全按照李航<<统计学习方法>>介绍的HMM代码,供大家参考,具体内容如下 #coding=utf8 ''''' Created on 2017-8-5 里面的代码许多地方可以精简,但为了百分百还原公式,就没有精简了. @author: adzhua ''' import numpy as np class HMM(object): def __init__(self, A, B, pi): ''''' A: 状态转移概率矩阵 B: 输出观察概率矩阵 pi: 初始化状态向量 '

  • python基于隐马尔可夫模型实现中文拼音输入

    在网上看到一篇关于隐马尔科夫模型的介绍,觉得简直不能再神奇,又在网上找到大神的一篇关于如何用隐马尔可夫模型实现中文拼音输入的博客,无奈大神没给可以运行的代码,只能纯手动网上找到了结巴分词的词库,根据此训练得出隐马尔科夫模型,用维特比算法实现了一个简单的拼音输入法.githuh地址:https://github.com/LiuRoy/Pinyin_Demo 原理简介隐马尔科夫模型 抄一段网上的定义: 隐马尔可夫模型 (Hidden Markov Model) 是一种统计模型,用来描述一个含有隐含未

  • Python实现隐马尔可夫模型的前向后向算法的示例代码

    本篇文章对隐马尔可夫模型的前向和后向算法进行了Python实现,并且每种算法都给出了循环和递归两种方式的实现. 前向算法Python实现 循环方式 import numpy as np def hmm_forward(Q, V, A, B, pi, T, O, p): """ :param Q: 状态集合 :param V: 观测集合 :param A: 状态转移概率矩阵 :param B: 观测概率矩阵 :param pi: 初始概率分布 :param T: 观测序列和状态

  • python实现马耳可夫链算法实例分析

    本文实例讲述了python实现马耳可夫链算法的方法.分享给大家供大家参考.具体分析如下: 在<程序设计实践>(英文名<The Practice of Programming>)的书中,第三章分别用C语言,C++,AWK和Perl分别实现了马耳可夫链算法,来通过输入的文本,"随机"的生成一些有用的文本. 说明: 1. 程序使用了字典,字典和散列可不是一个东西,字典是键值对的集合,而散列是一种能够常数阶插入,删除,不过可以用散列来实现字典. 2. 字典的setdef

  • 在Python上基于Markov链生成伪随机文本的教程

    首先看一下来自Wolfram的定义 马尔可夫链是随机变量{X_t}的集合(t贯穿0,1,...),给定当前的状态,未来与过去条件独立. Wikipedia的定义更清楚一点儿 ...马尔可夫链是具有马尔可夫性质的随机过程...[这意味着]状态改变是概率性的,未来的状态仅仅依赖当前的状态. 马尔可夫链具有多种用途,现在让我看一下如何用它生产看起来像模像样的胡言乱语. 算法如下, 找一个作为语料库的文本,语料库用于选择接下来的转换. 从文本中两个连续的单词开始,最后的两个单词构成当前状态. 生成下一个

  • python 用opencv实现霍夫线变换

    霍夫变换是一种检测任何形状的流行技术,可以检测形状,即使它被破坏或扭曲一点点. 一条线可以表示成y = mx + c或参数形式,像ρ=xcosθ+ysinθ,其中ρ是从原点到直线的垂直距离,θ角是由这条垂线和水平轴以逆时针的方向形成的(这个方向取决于你如何表示坐标系统,这种表示法在OpenCV中使用) OpenCV中的Hough变换 cv.HoughLines() 第一个参数,输入图像应该是一个二值图像,因此在应用hough变换之前应用阈值或使用Canny边缘检测. 第二和第三个参数分别是ρ和θ

  • Python实现将sqlite数据库导出转成Excel(xls)表的方法

    本文实例讲述了Python实现将sqlite数据库导出转成Excel(xls)表的方法.分享给大家供大家参考,具体如下: 1. 假设已经安装带有sliqte 库的Python环境 我的是Python2.5 2. 下载 python xls 写操作包(xlwt)并安装 下载地址: http://pypi.python.org/pypi/xlwt 3. 下面就是代码(db2xls.py): import sqlite3 as sqlite from xlwt import * #MASTER_COL

  • Python实现将目录中TXT合并成一个大TXT文件的方法

    本文实例讲述了Python实现将目录中TXT合并成一个大TXT文件的方法.分享给大家供大家参考.具体如下: 在网上下了一个dota的英雄攻略,TXT格式,每个英雄一个文件,看得疼,就写了一个小东西,合并一下. #coding=gbk import os import sys import glob def dirTxtToLargeTxt(dir,outputFileName): '''从dir目录下读入所有的TXT文件,将它们写到outputFileName里去''' #如果dir不是目录返回

  • python使用PythonMagick将jpg图片转换成ico图片的方法

    本文实例讲述了python使用PythonMagick将jpg图片转换成ico图片的方法.分享给大家供大家参考.具体分析如下: 这里使用到了PythonMagick模块,关于PythonMagick模块和ImageMagick的详细信息请参考:http://www.imagemagick.org/. 下面这段代码可以讲jpg图片转换成ico图标格式. # -*- coding: utf-8 -*- import PythonMagick img = PythonMagick.Image("c:/

随机推荐