python 生物信息学_生物信息学算法之Python实现
众所周知 基因组的核酸链不可能是随机形成的。有时候许多物种基因组之间 存在一些保守序列 motif 这意味着它们可能具有重要功能。但是 我们如何确定这些序列不是随机形成的 DNA 片段呢一个常识是 越短的序列越容易随机形成 越长的序列
众所周知,基因组的核酸链不可能是随机形成的。有时候许多物种基因组之间,存在一些保守序列(motif),这意味着它们可能具有重要功能。但是,我们如何确定这些序列不是随机形成的 DNA 片段呢?
一个常识是:越短的序列越容易随机形成,越长的序列越难随机形成。如何对随机形成序列的概率进行量化,以及如何确定容易和不容易随机形成的序列的长度的阈值呢?这篇文章将对这个问题进行探索。
给定: 一段 DNA 序列,以及一系列假定的 GC 出现的概率。
需得: 在特定 GC 出现的概率的情况下,得到一条与给定 DNA 序列 GC 含量相同的序列的概率,并且将概率值取对数输出。
示例数据
ACGATACAA
0.129 0.287 0.423 0.476 0.641 0.742 0.783
示例结果
-5.737 -5.217 -5.263 -5.360 -5.958 -6.628 -7.009
Python 实现
本题思路参考自下述博客:
Rosalind – Introduction to Random Strings[1]
因为 DNA 有 4 种碱基,每一个位置都有 4 种可能。如果每一种碱基出现的概率都是 25%,那么一个 9bp 的序列,共有 4·4·4·4·4·4·4·4·4 = 49=262144 种可能性。但现在我们假定 GC 出现的概率是 0.129 而不是 0.5,那么 1-0.129,即 0.871 就是 A 或 T 出现的概率。在此情况下,现在计算一个分子有多大概率得到与所给序列相同的 GC 含量?
ACGATACAA
0.129 0.287 0.423 0.476 0.641 0.742 0.783
下面我们在 Python 中演示如何计算。
dna = 'ACGATACAA'
A = '0.129 0.287 0.423 0.476 0.641 0.742 0.783'
prob_gc = map(float, A.split())
由于 G、C 单独出现的概率是 GC 同时出现的概率的一半,A、T 单独出现的概率是 AT 同时出现的概率的一半;又由于每一种碱基都是独立出现的,因此一条序列出现的概率是所有碱基分别出现的概率之乘积。
总的概率值 = 第1个碱基的概率 * 第2个碱基的概率 ……*最后一个碱基的概率
用代码表示就是:
gc = 0.129
at = 1 - 0.129
prob = 1.0
for base in dna:
if base in 'GC':
prob *= gc*0.5
else:
prob *= at*0.5
由于概率值都小于 1,多个概率值相乘,结果会越来越小,不如对结果取对数,以方便表示。
print(log10(prob))
又因为对数具有如下特性:
log(x·y) = log(x) + log(y)
因此,可以修改上面的代码:
gc = 0.129
at = 1 - gc
prob = 0.0
for base in dna:
if base in 'GC':
prob += log10(gc*0.5)
else:
prob += log10(at*0.5)
上面的代码也等价于:
gc = 0.129
at = 1 - gc
gc_num = dna.count('G') + dna.count('C')
at_num = len(dna) - gc_num
prob = gc_num * log10(gc*0.5) + at_num * log10(at*0.5)
完整的代码是:
Introduction_to_Random_Strings.py
import sys
from math import log10
#dna = 'ACGATACAA'
#A = '0.129 0.287 0.423 0.476 0.641 0.742 0.783'
with open('rosalind_prob.txt') as fh:
dna, A = [l.strip() for l in fh.readlines()]
gc_prob = map(float, A.split())
gc_num = dna.count('G') + dna.count('C')
at_num = len(dna) - gc_num
B = []
for gc in gc_prob:
at = 1 - gc
p = gc_num * log10(gc*0.5) + at_num * log10(at*0.5)
B.append(p)
print(' '.join([str(round(i, 3)) for i in B]))
注意:
我们计算概率的序列,只是 GC 含量与给定序列相同,并没有考虑碱基顺序;
由结果可知:当 GC 出现概率与给定 DNA 的 GC 含量相同时,出现与给定 DNA 的 GC 含量相同的序列的概率最大。
参考资料
[1]
Rosalind – Introduction to Random Strings: https://recologia.com.br/2016/06/08/rosalind-introduction-to-random-strings/
喜欢文章请点个“赞”吧!或者点击“在看”让更多朋友看到,点击“阅读原文”可以在知乎专栏上给我留言,博客地址:https://jianzuoyi.github.io
以上信息来源于网络,如有侵权,请联系站长删除。
python 生物信息学_生物信息学算法之Python实现相关推荐
- 分类算法python程序_分类算法——k最近邻算法(Python实现)(文末附工程源代码)...
kNN算法原理 k最近邻(k-Nearest Neighbor)算法是比较简单的机器学习算法.它采用测量不同特征值之间的距离方法进行分类,思想很简单:如果一个样本在特征空间中的k个最近邻(最相似)的样 ...
- svm算法python实现_(转载)python应用svm算法过程
除了在Matlab中使用PRTools工具箱中的svm算法,Python中一样可以使用支持向量机做分类.因为Python中的sklearn库也集成了SVM算法,本文的运行环境是Pycharm. 一.导 ...
- pythoncookbook和流畅的python对比_为什么你学Python效率比别人慢?因为你没有这套完整的学习资料...
以下资源免费获取方式! 关注!转发!私信"资料"即可免费领取! 入门书籍 1.<Python基础教程>(Beginning Python From Novice to ...
- 零基础学python 视频_全网最全Python视频教程真正零基础学习Python视频教程 490集...
Python Web开发-进阶提升 490集超强Python视频教程 真正零基础学习Python视频教程 [课程简介] 这是一门Python Web开发进阶课程,手把手教你用Python开发完整的商业 ...
- 3 x 10的python表达式_这道数学题用PYTHON编程语言怎么写? 编程语言python是用
我觉着,这个应该这样解决比较符合计算机解题思路. 下面的回答的,思考的东西太多. # -*- coding: utf-8 -*- __author__ = 'lpe234' __date__ = '2 ...
- python图像分类_用于实现用python和django编写的图像分类的Keras UI
KerasUI是一种可视化工具,可以在图像分类中轻松训练模型,并允许将模型作为服务使用,只需调用API. https://github.com/zeppaman/KerasUI 主要特点: 用oaut ...
- 小象python培训班_小象最新Python机器学习升级版视频学习教程 共24节精品课
小象最新Python机器学习升级版视频学习教程 共24节精品课 本课程特点是从数学层面推导最经典的机器学习算法,以及每种算法的示例和代码实现(Python).如何做算法的参数调试.以实际应用案例分析各 ...
- python初学者_初学者的热门Python资源
python初学者 Python新手? 还是您已经是一位经验丰富的开发人员,希望增加和提高您的Python知识? 我们为希望学习Python编程的任何人编制了一份推荐资源的书包. 我们对这些资源进行了 ...
- 为什么要学python语言_我们为什么要学习Python语言?
原标题:我们为什么要学习Python语言? 聊到我们为什么要学习Python语言?小编不禁又想起大佬潘石屹准备开启Python学习旅程时所发布的微博. 我们为什么要学习Python语言? 在农业社会时 ...
- 下载python步骤_下载及安装Python详细步骤
安装python分三个步骤: *下载python *安装python *检查是否安装成功 1.下载python (1)python下载地址 (2)选择下载的版本 (3)点开download后,找到下载 ...
最新文章
- pytorch 优化GPU显存占用,避免out of memory
- [NodeJs] npm提供了哪些钩子?各有什么作用?
- box-sizing属性
- Missing iOS Distribution signing identity问题解决
- outlook域用户名怎么填_家谱制作软件怎么做成电子版
- python 数据结构转换层_python – 具有Maxpooling1D和channel_first的Keras模型
- 牛津词典 2018 年度词汇 ——「有毒」!
- Powerpoint 插件制作日记-1
- h3c 链路聚合测试_H3CSE学习之链路聚合
- arping命令解析
- 机械制造技术类毕业论文文献都有哪些?
- 图解计算机基础网站上线了
- selenium元素模糊定位xpath contains、starts-with和ends-with
- 蓝牙BLE协议分析【附代码实例】
- PO*创建标准采购订单
- 使用.mdf和.ldf文件还原sqlserver数据库
- python画树林_在Python 3中使用深度森林(Deep Forest)进行分类
- 杰卡德相似系数(Jaccardsimilarity coefficient)
- Bitmap 处理图片修改为透明背景,改变主颜色
- 分布式架构实现概述(大型网站技术架构-读后感)