前往小程序,Get更优阅读体验!
立即前往
首页
学习
活动
专区
工具
TVP
发布
社区首页 >专栏 >用 Python 玩转常用生物序列

用 Python 玩转常用生物序列

作者头像
生信菜鸟团
发布2021-04-30 16:07:52
1.7K0
发布2021-04-30 16:07:52
举报
文章被收录于专栏:生信菜鸟团生信菜鸟团

一、准备工作

1、获取感兴趣的基因,蛋白质,转录本等生物序列 FASTA 或 GenBank

这里举例,进入 NCBI 获取的GeneBank / FASTA 的数据格式

比如查看 POU5F1 基因:https://www.ncbi.nlm.nih.gov/gene/5460

2、搭建 Python 环境与项目目录

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

搭建目录结构及Python环境参考:https://blog.csdn.net/u011262253/article/details/105902060

二、操作生物序列

1、读取常见的序列文件格式(fasta,gb)
代码语言:javascript
复制
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 序列文件内容
代码语言:javascript
复制
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 文件详细信息
代码语言:javascript
复制
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、新建序列文件
代码语言:javascript
复制
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

代码语言:javascript
复制
dna_seq[0] = "G"

如果你一定要修改,当然也可以,但不建议这么做

代码语言:javascript
复制
dna_seq_mutable = dna_seq.tomutable()
dna_seq_mutable[0] = "G"
6、操作序列文件
代码语言:javascript
复制
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、将 DNA 翻译为 RNA
代码语言:javascript
复制
# =====转录=====
# 如果序列为编码链,那么直接转换
print ("rna: ", dna_seq.transcribe())
# 如果序列为模板链,就需要先转为编码链
transcribe_seq = dna_seq.reverse_complement().transcribe()
print ("rna: ", transcribe_seq)
8、将 RNA 翻译为 蛋白质
代码语言:javascript
复制
# =====翻译=====
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、做一些有意思的事
代码语言:javascript
复制
# =====寻找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])
本文参与 腾讯云自媒体同步曝光计划,分享自微信公众号。
原始发表:2021-04-28,如有侵权请联系 cloudcommunity@tencent.com 删除

本文分享自 生信菜鸟团 微信公众号,前往查看

如有侵权,请联系 cloudcommunity@tencent.com 删除。

本文参与 腾讯云自媒体同步曝光计划  ,欢迎热爱写作的你一起参与!

评论
登录后参与评论
0 条评论
热度
最新
推荐阅读
目录
  • 1、读取常见的序列文件格式(fasta,gb)
  • 2、浏览 fasta 序列文件内容
  • 3、读取 genebank 文件详细信息
  • 4、新建序列文件
  • 5、修改序列文件
  • 6、操作序列文件
  • 7、将 DNA 翻译为 RNA
  • 8、将 RNA 翻译为 蛋白质
  • 7、做一些有意思的事
相关产品与服务
对象存储
对象存储(Cloud Object Storage,COS)是由腾讯云推出的无目录层次结构、无数据格式限制,可容纳海量数据且支持 HTTP/HTTPS 协议访问的分布式存储服务。腾讯云 COS 的存储桶空间无容量上限,无需分区管理,适用于 CDN 数据分发、数据万象处理或大数据计算与分析的数据湖等多种场景。
领券
问题归档专栏文章快讯文章归档关键词归档开发者手册归档开发者手册 Section 归档