Python一阶马尔科夫链生成随机DNA序列实现示例
作者:BioChem 发布时间:2021-06-23 07:42:22
标签:Python,一阶马尔科夫,DNA
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
来源:https://www.jianshu.com/p/2abbba19bf4e
0
投稿
猜你喜欢
- 后台数据库: [Microsoft Access] 与 [Microsoft Sql Server] 更换之后,ASP代码应注意要修改的一些
- 前言:本篇文章基于卷积神经网络CNN,使用PyTorch实现MNIST数据集手写数字识别。一、PyTorch是什么?PyTorch 是一个
- 前言人类都是视觉动物,不管是男生还是女生看到漂亮的小姐姐、小哥哥就想截图保存下来。可是截图会对画质会产生损耗,截取的画面不规整,像素不高等问
- monfs :我想知道javascript是否可以实现这样的功能来改变本地的IP地址,例如我本地设置的IP地址是192.168.0.1,我想
- 环境:RHEL 5.4 x86 , oracle 11.2 1.设定环境变量 在/home/oracle编辑 # .bash_profile
- 现在很多以内容为核心的网站上都在文章底部添加了社会化分享按钮,能让浏览用户在发现一篇有价值的文章时,可以通过社会化网络快速分享给自己的好友,
- 问题重现Installing golang.org/x/tools/cmd/guru FAILED Installing golang.or
- 什么是1433端口 1433端口,是SQL Server默认的端口,SQL Server服务使用两个端口:TCP-1433、UDP-1434
- 代码如下:DECLARE @c INT DECLARE @c2 INT SELECT @c = COUNT(1) FROM dbo.Spli
- 一、用HTTP头信息 也就是用PHP的HEADER函数。PHP里的HEADER函数的作用就是向浏览器发出由HTTP协议规定的本来应该通过WE
- 登录、注销和登录限制:登录在使用authenticate进行验证后,如果验证通过了。那么会返回一个user对象,拿到user对象后,可以使用
- python的多重继承的理解Python和C++一样,支持多继承。概念虽然容易,但是困难的工作是如果子类调用一个自身没有定义的属性,它是按照
- 如何使用Iframe实现本页提交?例:chunfeng.html< html>< head>&n
- 我发布了一个通过FTP自动优化新图像的教程。这次我们将抓取整个网站,并在本地优化我们遇到的图像,按URL组织。请注意,这个简短但中级的脚本不
- 这篇文章主要介绍了python isinstance函数用法详解,文中通过示例代码介绍的非常详细,对大家的学习或者工作具有一定的参考学习价值
- 最近老师布置了个作业,爬取豆瓣top250的电影信息。按照套路,自然是先去看看源代码了,一看,基本的信息竟然都有,心想这可省事多了。简单分析
- 在循环对象和函数对象中,我们了解了循环器(iterator)的功能。循环器是对象的容器,包含有多个对象。通过调用循环器的next()方法 (
- 本文实例讲述了PHP数据库表操作的封装类及用法。分享给大家供大家参考,具体如下:数据库表结构:CREATE TABLE `test_user
- CSS重设就是由于各种浏览器解释CSS样式的初始值有所不同,导致设计师在没有定义某个CSS属性时,不同的浏览器会按照自己的默认值来为没有定义
- 目录1、方式一:自动创建2、方式二:纯手动创建3、方式三:半自动创建1、方式一:自动创建# django orm 自动帮我们创建第三张表,我