上一篇文章生物信息中的Python 01 | 从零开始处理基因序列自己造轮子实现了序列的基础操作,但是在Python的世界里,一项工作只要重复的次数多了,那么一定就会有大神来开发相应的包来解决,这个包名就是 Biopython 。接下来我们试着使用它来实现简单的序列处理。

一、准备工作

1、 按照上一篇下载fasta文件的步骤,可以同理得到GeneBank的数据格式

2、现在我们的目录结构是这样的

搭建下面的目录结构参考:搭建 Python 高效开发环境: Pycharm + Anaconda

3、安装Biopython,这里有两种方案:

3.1 用pip安装Biopython,在cmd命令窗口输入
  • 下载Python的包管理工具:pip

https://pypi.org/project/pip/#files

  • 下载完,解压,进入解压目录

    • Linux 下输入
    sudo python setup.py install
    
    • windows 下,在下载目录,Shift+右键

    • 如下图所示,点击在此处打开命令窗口

    • 输入如下命令

    python setup.py install
    

  • 测试是否安装成功,出现下图所示的提示即表示安装成功

    pip -v
    

  • 进入 Pycharm 的Terminal 窗口,输入以下命令来安装 Biopython

pip install biopython

3.2 直接用安装包安装
  • window系统

    • 下载地址:http://biopython.org/DIST/biopython-1.72.tar.gz
    • 解压
    • 按住shift并点击右键
    • 在菜单栏点击在此处打开命令窗口,并输入如下命令:python setup.py install
  • Linux系统

    • 打开终端 (快捷键:Ctrl+Alt+T

    • 在终端输入以下命令

      $ wget http://biopython.org/DIST/biopython-1.72.tar.gz
      $ tar -zxvf biopython-1.72.tar.gz
      $ cd biopython-1.72/
      $ sudo python setup.py install
      
    • 测试是否安装成功

      $ python
      >>> from Bio.Seq import Seq
      >>> seq = Seq('ATCG')
      >>> seq
      Seq('ATCG')
      

二、Biopython 基础用法

1 读取常见的序列文件格式(fasta,gb)
from Bio import SeqIO# 读取包含单个序列 Fasta 格式文件
fa_seq = SeqIO.read("res/sequence1.fasta", "fasta")
# print fa_seq
# 读取包含多个序列的 fasta 格式文件
for fa in SeqIO.parse("res/multi.fasta", "fasta"):print (fa.seq)
# 一个多序列文件中的所有序列
seqs = [fa.seq for fa in SeqIO.parse("res/multi.fasta",  "fasta")]
print (seqs)
# 如果不想要seq对象中的字母表,可以用str()来强制类型转换
seqs = [str(fa.seq) for fa in SeqIO.parse("res/multi.fasta",  "fasta")]
print (seqs)# 读取包含单个序列的 gb 格式文件
gb_seq = SeqIO.read("res/sequence1.gb", "genbank")
print (gb_seq)
2 浏览 fasta 序列文件内容
from Bio import SeqIO# 读取包含单个序列 Fasta 格式文件
fa_seq = SeqIO.read("res/sequence1.fasta", "fasta")# =====获取详细的信息=====
# 提取基因ID,name
# Fasta 文件中序列名所在行的第一个词被作为 id 和 name
print ("id: ", fa_seq.id)
print ("name: ", fa_seq.name)
# 基因 Description 是fasta文件格式中的第一行
print ("description: ", fa_seq.description)
# 序列
print ("seq: ", fa_seq.seq)
# 序列来源库信息(NCBI的数据库信息会包括数据库交叉引用)
print ("dbxrefs: ", fa_seq.dbxrefs)
# 全部序列的注释信息
print ("annotations: ", fa_seq.annotations)
# 序列中每个字母的注释信息
print ("letter_annotations: ", fa_seq.letter_annotations)
# 部分序列的注释信息
print ("features: ", fa_seq.features)
3 浏览 genebank 序列文件内容
from Bio import SeqIO# 读取包含单个序列的 gb 格式文件
gb_seq = SeqIO.read("res/sequence1.gb", "genbank")
print (gb_seq)# =====获取详细的信息=====
# 提取基因ID,name
# gb文件中序列名包含比fasta更加详细的序列信息,下面分别是 id 和 name
print ("id: ", gb_seq.id)
print ("name: ", gb_seq.name)
# 基因 Description 是fasta文件格式中的第一行
print ("description: ",  gb_seq.description)
# 序列信息, 这里的序列信息是以 bioPython 中的seq对象存储
print ("seq: ", gb_seq.seq)
# 序列来源库信息(NCBI的数据库信息会包括数据库交叉引用)
print ("dbxrefs: ", gb_seq.dbxrefs)
# 全部序列的注释信息
print ("annotations: ", gb_seq.annotations)
# 序列中每个字母的注释信息
print ("letter_annotations: ", gb_seq.letter_annotations)
# 部分序列的注释信息,SeqFeature 对象的形式保存了features table中的所有entries(如genes和CDS等)
print ("features: ", gb_seq.features)
# 该基因的物种信息
print ("organism: ", gb_seq.annotations["organism"])
# 关于序列的注释信息,相关数据库的交叉引用号
print ("comment: ", gb_seq.annotations["comment"])
# 序列来源的物种名
print ("source: ", gb_seq.annotations["source"])
# 该基因的分类学信息
print ("taxonomy: ", gb_seq.annotations["taxonomy"])
# 该基因的整理后的注释信息
print ("structured_comment: ", gb_seq.annotations["structured_comment"])
# 该基因序列相关的关键词
print ("keywords: ", gb_seq.annotations["keywords"])
# 该基因的相关文献编号,或递交序列的注册信息
print ("references: ", gb_seq.annotations["references"])
# 该基因的入库时,给的基因编号,以及在染色体上的位点信息
print ("accessions: ", gb_seq.annotations["accessions"])
# 该基因的分子类型,一般为 DNA
print ("molecule_type: ", gb_seq.annotations["molecule_type"])
# 该基因的数据文件划分方式
print ("data_file_division: ", gb_seq.annotations["data_file_division"])
# 基因发布时间
print ("date: ", gb_seq.annotations["date"])
# 该基因的更新版本
print ("sequence_version: ", gb_seq.annotations["sequence_version"])
# 该基因的拓扑结构
print ("topology: ", gb_seq.annotations["topology"])

相信大家可以看到 GeneBank 比 fasta 格式更加详细和贴心,但是对于序列处理来说内存占用和运行时间比这些信息更加重要。这就使fasta成为我们一般在序列分析中常用的格式。

4 新建序列文件
from Bio.Seq import Seq# 新建一个DNA序列对象
dna_seq = Seq("GGATGGTTGTCTATTAACTTGTTCAAAAAAGTATCAGGAGTTGTCAAGGCAGAGAAGAGAGTGTTTGCA", IUPAC.unambiguous_dna)
# 新建一个RNA序列对象
rna_seq = Seq("GGATGGTTGTCTATTAACTTGTTCAAAAAAGTATCAGGAGTTGTCAAGGCAGAGAAGAGAGTGTTTGCA", IUPAC.unambiguous_rna)
# # 新建一个蛋白质序列对象
protein_seq = Seq("GGATGGTTGTCTATTAACTTGTTCAAAAAAGTATCAGGAGTTGTCAAGGCAGAGAAGAGAGTGTTTGCA", IUPAC.protein)

序列对象由一段字符串和其对应的编码表所定义。我们可以从上述的代码中看到,字符串内容一样,唯一不同的就是第二个参数IUPAC值不一样。IUPAC (International Union of Pure and Applied Chemistry ) 是一个制定化学相关标准的组织,Biopython 所使用的编码表就是由它制定的,想了解详细细节可以参考http://www.bioinformatics.org/sms2/iupac.html ,详细定义如下:

名称 编码表
ambiguous_dna_letters GATCRYWSMKHBVDN
unambiguous_dna_letters GATC
ambiguous_rna_letters GAUCRYWSMKHBVDN
unambiguous_rna_letters GAUC
protein ARNDCQEGHILKMFPSTWYV
5 修改序列文件

在生物学意义上,序列是不可以随便更改的,也就是不可变的。如果强行修改,那么就会报错TypeError: 'Seq' object does not support item assignment

dna_seq[0] = "G"

如果你执意修改也是可以的,但是不建议这么做

dna_seq_mutable = dna_seq.tomutable()
dna_seq_mutable[0] = "G"
6 操作序列文件
from Bio.Seq import Seq# 新建一个DNA序列对象
dna_seq = Seq("GGATGGTTGTCTATTAACTTGTTCAAAAAAGTATCAGGAGTTGTCAAGGCAGAGAAGAGAGTGTTTGCA", IUPAC.unambiguous_dna)
# 序列信息
print ("Sequence: ", dna_seq)
# 序列长度
print ("Length : ", len(dna_seq))
# 单个核苷酸计数
print ("G Counts: ", dna_seq.count("G"))
# 获取反向序列
print ("reverse: ", dna_seq[::-1])
# 获取反向互补序列
print ("Reverse complement: ", dna_seq.complement())
# 获取蛋白质的反向互补序列,这里显然是报错的,因为蛋白序列没有这一属性
print ("Protein reverse complement: ", protein_seq.complement())       
7 用 Biopython 将 DNA 翻译为 RNA
# =====转录=====
# 如果序列为编码链,那么直接转换
print ("rna: ", dna_seq.transcribe())
# 如果序列为模板链,就需要先转为编码链
transcribe_seq = dna_seq.reverse_complement().transcribe()
print ("rna: ", transcribe_seq)
8 用BioPython 将 RNA 翻译为 蛋白质
# =====翻译=====
print ("protein: ", transcribe_seq.translate())
# 如果翻译的是线粒体密码子,那么在参数中需要输入,其他参考 https://www.ncbi.nlm.nih.gov/Taxonomy/Utils/wprintgc.cgi?mode=c or ftp://ftp.ncbi.nlm.nih.gov/entrez/misc/data/gc.prt
print ("protein: ", transcribe_seq.translate(table="Vertebrate Mitochondrial"))
# 在现实生物世界中,一般在遇到终止密码子之后的序列不用翻译
print ("protein: ", transcribe_seq.translate(table="Vertebrate Mitochondrial", to_stop=True))
# 如果DNA序列为编码序列,可以直接翻译,DNA序列不是3的倍数时,报错
print ("protein: ", dna_seq.translate())
# 在细菌世界中,在细菌遗传密码中 GTG 是个有效的起始密码子,注意第一个密码子(正常情况下 GTG编码缬氨酸, 但是如果作为起始密码子,则翻译成甲硫氨酸)
bacterial_dna = Seq("GTGAAAAAGATGCAATCTATCGTACTCGCACTTTCCCTGGTTCTGGTCGCTCCCATGGCATAA", generic_dna)
print  ("protein: ", bacterial_dna.translate(table="Bacterial", to_stop=True))
print  ("protein: ", bacterial_dna.translate(table="Bacterial", cds=True))
7 做一些有意思的事
# =====寻找TATA框=====
# TATA框约在多数真核生物基因转录起始点上游约-30bp(-25~-32bp)处,基本上由A-T碱基对组成,是决定基因转录始的选择,为RNA聚合酶的结合处之一
print ("TA Counts: ", dna_seq.count("TA"))# =====GC含量=====
# (A+T)/(G+C)之比随DNA的种类不同而异。GC含量愈高,DNA的密度也愈高,同时热及碱不易使之变性,因此利用这一特性便可进行DNA的分离或测定。
print ("GC Contenten", 100 * float(dna_seq.count("G") + dna_seq.count("C")) / len(dna_seq))# =====得到promoter序列=====
# 在寻找基因的promoter时(一般promoter的位点不确定),但是可以通过将起始位点左右2kb基因视为promoter
# 这里训练切取,将切取设起始位点为前10bp
print ("Promoter seq: ",dna_seq[:10])

生物信息中的Python 03 | 自动化操作NCBI

生物信息中的Python 04 | 批量下载基因与文献

生物信息中的Python 02 | 用biopython解析序列相关推荐

  1. python中序列和列表区别细菌真菌病毒_生物信息中的Python 02 | 用biopython解析序列...

    上一篇文章生物信息中的Python 01 | 从零开始处理基因序列自己造轮子实现了序列的基础操作,但是在Python的世界里,一项工作只要重复的次数多了,那么一定就会有大神来开发相应的包来解决,这个包 ...

  2. 生物信息中的Python 01 | 从零开始处理基因序列

    一. 序列数据的下载 在开始了解序列的处理流程时,我们先要知道序列下载网址.其中一个知名的网站就是NCBI (National Center for Biotechnology Information ...

  3. 利用python处理dna序列_Python + 生物信息 02 :Biopython 分析序列

    Biopython 做序列分析 一.安装Biopython:如果环境已经有Biopython可以跳过这一步.这里有两种安装方案,一种通过pip快速安装,另一种通过安装包安装 1. 用pip安装Biop ...

  4. 生物信息中的Python 05 | 从 Genbank 文件中提取 CDS 等其他特征序列

    1 介绍 在基因结构分析或其他生物功能分析中会时常用到 CDS 序列,以及其他诸如 mRNA 序列,misc RNA序列等具有生物意义的序列片段.而NCBI 的基因库中已经包含有这些的信息,但是只有一 ...

  5. Python 中 xpath 语法 与 lxml 库解析 HTML/XML 和 CSS Selector

    The lxml.etree Tutorial :https://lxml.de/tutorial.html python3 解析 xml:https://www.cnblogs.com/deadwo ...

  6. c++中调用python脚本提示 error LNK2001: 无法解析的外部符号 __imp_Py_Initialize等错误的解决方法

    c++中调用python脚本提示 error LNK2001: 无法解析的外部符号 __imp_Py_Initialize等错误的解决方法 时间:2017-05-09 12:32:06阅读:234评论 ...

  7. python中set函数_python中set()函数简介及实例解析

    set函数也是python内置函数的其中一个,属于比较基础的函数.其具体介绍和使用方法,下面进行介绍. set() 函数创建一个无序不重复元素集,可进行关系测试,删除重复数据,还可以计算交集.差集.并 ...

  8. Python中sort和sorted函数代码解析

    Python中sort和sorted函数代码解析 本文研究的主要是Python中sort和sorted函数的相关内容,具体如下. 一.sort函数 sort函数是序列的内部函数 函数原型: L.sor ...

  9. python 可以根据元素值删除的是_python中删除某个元素的方法解析

    这篇文章主要介绍了python中删除某个元素的方法解析,文中通过示例代码介绍的非常详细,对大家的学习或者工作具有一定的参考学习价值,需要的朋友可以参考下 python中关于删除list中的某个元素,一 ...

最新文章

  1. elif在python中啥意思_python中elif 结构语句怎么判断?
  2. 霍夫变换概述和标准霍夫变换
  3. 如何创建 Code Snippet
  4. centos java发送邮件发不出去_传真机发不出传真怎么办 传真机发不出传真解决方法【详解】...
  5. 2020\Simulation_1\2.约数个数
  6. 对于大规模机器学习的理解和认识
  7. kubenetes中port、targetPort、nodePort、containerPort的区别与联系
  8. Zigbee费尽心思做mesh网究竟在智能家居中有什么用?
  9. mingw64+msys2下使用cmake问题
  10. 国防现代化的数据_Linux容器如何解决国防虚拟化问题
  11. java命令大全_Java命令行工具:javac、java、javap 的使用详解
  12. Mate 50系列首发?曝鸿蒙3.0用户版5月内测
  13. 对象的自身引用(Self-Reference) 动态绑定(Dynamic Binding)
  14. Spring教程动画文字版
  15. netty权威指南 微云_《Netty权威指南》(二)NIO 入门
  16. 开源音乐软件——落雪
  17. 6-32 表头插入法构造链表
  18. libxml2的参考手册
  19. Java中 VO、PO、DO、DTO、BO、QO、DAO、POJO的概念
  20. python unicode error_python-ValueError:操作参数必须为str或unicode

热门文章

  1. 百度语音合成 java 教程_【百度语音合成】JavaAPI方式语音合成示例
  2. 单点登录、域用户、常规登录、AD域
  3. java生成word文档 图片_java生成带有图片的word的文档-Go语言中文社区
  4. Linux的时间和时区设置
  5. Centos7安装Python
  6. su: warning: cannot change directory to /home/mysql: No such file or directory
  7. springboot整合全文搜索引擎Elasticsearch | Spring Boot 28
  8. Sa函数 与 sinc函数
  9. sql查找数据中某个字段是否有重复的值
  10. 多多自走棋改动_多多自走棋:20日更新,刺客、光羽修改,装备小幅调整