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


猜你喜欢
- PyCharm 是我用过的python编辑器中,比较顺手的一个。而且可以跨平台,在macos和windows下面都可以用,这点比较好。是py
- Windows10本机环境:win10未安装node,安装了nvm工具,尝试使用nvm安装我开发环境的node版本(10.9.0 or la
- 一个是没有对输入的数据进行过滤(过滤输入),还有一个是没有对发送到数据库的数据进行转义(转义输出)。这两个重要的步骤缺一不可,需要同时加以特
- 做项目的时候,一位同事导数据的时候,不小心把一个表中的数据全都搞重了,也就是说,这个表里所有的记录都有一条重复的。这个表的数据是千万级的,而
- 用途:图片经常使用onload来改变宽度,但这样会出现图片的闪烁,这个简单的类就是用来解决这个问题的。函数loadImage,用来加载图片,
- 目录Java实现上传Excel文件并导出到数据库1、导入依赖2、domain3、utils4、Controller5、xmlJava实现上传
- python3标准库里自带线程池ThreadPoolExecutor和进程池ProcessPoolExecutor。如果你用的是python
- 操作系统上正确配置python环境之后,pyuic5也是一个可以识别的命令行指令到.ui文件的目录下,直接cmd进入,输入pyuic5 -o
- 思路:先随机排序然后再分组就好了。1、创建表:CREATE TABLE `xdx_test` ( `id` int(11) NOT NULL
- 最近我遇到了一个 bug,我试着通过 Rails 在以“utf8”编码的 MariaDB 中保存一个 UTF-8 字符串,然后出现了一个离奇
- 本文实例为大家分享了SpringBoot整合Mybatis使用Druid数据库连接池的方法,具体内容如下在SpringBoot项目中,增加如
- 说明:操作系统:Windows Server 2003MySQL版本:5.5.25MySQL程序安装目录:D:\Program Files\
- 背景和目的:利用python request 编写脚本测试公司系统的文件上传接口。前端读取文件的大小然后文件分片传给后端,后端将每一片数据重
- sympy版本:1.2假设求解矩阵方程AX=A+2X其中求解之前对矩阵方程化简为(A−2E)X=A令B=(A−2E)使用qtconsole输
- 分区视图联接来自一组成员的水平分区数据,使数据看起来象来自同一张表。SQL Server 2000 区分本地分区视图和分布式分区视图。在本地
- 标量标量由普通小写字母表示(例如,x、y和z)。我们用 R \mathbb{R} R表示所有(连续)实数标量的空间。标量由只有一个元素的张量
- Macromedia Dreamweaver MX 2004提供了更多功能强劲的可视化设计工具、应用开
- 在这篇文章里,我们会聊一聊为什么 Python 决定不支持 switch 语句。为什么想要聊这个话题呢?主要是因为 switch 在其它语言
- 最近要使用python调用C++编译生成的DLL动态链接库,因此学习了一下ctypes库的基本使用。ctypes是一个用于Python的外部
- <!DOCTYPE html PUBLIC "-//W3C//DTD XHTML 1.0 Transitional//EN&