使用biopython处理序列数据

欢迎关注”生信修炼手册”!
序列是基因组学数据的基本单位,对于序列先关信息的存储,有以下两种常用的文件格式

1. fasta

2. genebank

通过biopython, 我们可以方便的读取这些格式的文件,并提取其中的信息。具体地,通过以下3个子模块来处理序列数据

1. Bio.Seq

2. Bio.SeqRecore

3. Bio.SeqIO

其中Bio.Seq表示最原始的序列对象,是最核心的模块,提供了序列的格式化,反向互补,碱基计数等基本功能;Bio.SeqRecord表示序列记录,在序列对象的基础上,进一步添加了序列的id, 名称,属性等各种注释信息;Bio.SeqIO模块则用于读取特定的文件格式,返回 SeqRecord对象。

1. Bio.Seq

Bio.Seq提供了最核心的序列对象,即由基本字符构成的序列,比如核酸序列和蛋白质序列,初始化方式如下

>>> from Bio.Seq import Seq
>>> my_seq = Seq('ATCGTACGATCT')
>>> my_seq
Seq('ATCGTACGATCT')
在该模块中,为序列对象提供了python字符的基础操作,比如比较,大小写转换,切片,切分,连接, 格式化等操作,具体用法如
# 切片
>>> my_seq[1]
'T'
>>> my_seq[1:3]
Seq('TC')
>>> my_seq[::-1]
Seq('TCTAGCATGCTA')
# 小写转换
>>> my_seq.lower()
Seq('atcgtacgatct')
# 大写转换
>>> my_seq.upper()
Seq('ATCGTACGATCT')
# split, 序列分隔
>>> my_seq.split('A')
[Seq(''), Seq('TCGT'), Seq('CG'), Seq('TCT')]
# join, 序列连接
>>> my_seq2 = Seq('ACGACTGACTAGCT')
>>> Seq('NNN').join([my_seq, my_seq2])
Seq('ATCGTACGATCTNNNACGACTGACTAGCT')
# 格式化
>>> 'id:1,seq:{}'.format(my_seq)
'id:1,seq:ATCGTACGATCT'

除了基本用法之外,还提供了生物学序列特有的功能,比如互补,反向互补,转录,翻译等,具体用法如下

# 互补
>>> my_seq.complement()
Seq('TAGCATGCTAGA')
# 反向互补
>>> my_seq.reverse_complement()
Seq('AGATCGTACGAT')
# 转录
>>> my_seq.transcribe()
Seq('AUCGUACGAUCU')
# 翻译
>>> my_seq.translate()
Seq('IVRS')

2. Bio.SeqRecord

Bio.SeqRecord在序列的基础上,进一步存储了相关的注释信息,初始化的方式如下

>>> from Bio.SeqRecord import SeqRecord
>>> my_seq = Seq('AGCTACGT')
>>> my_seqrecord = SeqRecord(my_seq)

SeqRecord具有id, name等多个属性,示例如下

>>> my_seqrecord
SeqRecord(seq=Seq('AGCTACGT'), id='<unknown id>', name='<unknown name>', description='<unknown description>', dbxrefs=[])
>>> my_seqrecord.seq
Seq('AGCTACGT')
>>> my_seqrecord.id
'<unknown id>'
>>> my_seqrecord.name
'<unknown name>'
>>> my_seqrecord.description
'<unknown description>'

除了以上基本属性外,还具有annotations和lette_annotations两个属性,进一步丰富了注释信息,annotations属性是一个字典结构,通过key=value的形式可以存储不同类别的注释信息,letter_annotations属性也是一个字典结构,但是其中的value值是长度等于序列长度的列表,主要用于存储每个碱基对应的信息,示例如下

>>> my_seqrecord.annotations['organ'] = 'human'
>>> my_seqrecord.annotations
{'organ': 'human'}
>>> my_seqrecord.letter_annotations['quality'] = [20, 20, 20, 20, 20, 20, 20, 20]

3. Bio.SeqIO

Bio.SeqIO用于文件的读写,支持多种文件格式,对于序列的存储格式fasta和genebank而言,读取的方式如下

>>> from Bio import SeqIO
>>> for seq in SeqIO.parse('input.fasta', 'fasta'):
...     print(seq.id, seq.seq)
>>> for seq in SeqIO.parse('input.gb', 'genebank'):
...     print(seq.id, seq.seq)

在每个for循环中,返回的是SeqRecord对象,可以通过SeqRecord对象的方法来访问各种信息。除了for循环的遍历,也可以直接返回列表,示例如下

>>> records = list(SeqIO.parse('input.fasta', 'fasta'))
>>> records[0]
SeqRecord(seq=Seq('CGATCGATCGACT'), id='1', name='1', description='1', dbxrefs=[])

该模块也支持序列对象的写入操作,最典型的应用就是序列格式的转换,genebank转换为fasta格式,代码如下

>>> records = SeqIO.parse("input.gb", "genbank")
>>> SeqIO.write(records, "out.fasta", "fasta")

write方法提供了输出功能,将序列对象输出到指定格式的文件中,针对格式转换这一常见场景,用法如下

>>> count = SeqIO.convert("input.gb", "genbank", "out.fasta", "fasta")

以上3个子模块层层渐进,构建了biopython处理序列数据的完整生态,对于使用者而言,通过简单的几句代码,就可以完成基本的序列操作,对于开发者而言,其class的抽象设计,方法编写都值得参考借鉴。

·end·
—如果喜欢,快分享给你的朋友们吧—

原创不易,欢迎收藏,点赞,转发!生信知识浩瀚如海,在生信学习的道路上,让我们一起并肩作战!

本公众号深耕耘生信领域多年,具有丰富的数据分析经验,致力于提供真正有价值的数据分析服务,擅长个性化分析,欢迎有需要的老师和同学前来咨询。
  更多精彩
  • KEGG数据库,除了pathway你还知道哪些

  • 全网最完整的circos中文教程

  • DNA甲基化数据分析专题

  • 突变检测数据分析专题

  • mRNA数据分析专题

  • lncRNA数据分析专题

  • circRNA数据分析专题

  • miRNA数据分析专题

  • 单细胞转录组数据分析专题

  • chip_seq数据分析专题

  • Hi-C数据分析专题

  • HLA数据分析专题

  • TCGA肿瘤数据分析专题

  • 基因组组装数据分析专题

  • CNV数据分析专题

  • GWAS数据分析专题

  • 2018年推文合集

  • 2019年推文合集

  写在最后
转发本文至朋友圈,后台私信截图即可加入生信交流群,和小伙伴一起学习交流。

扫描下方二维码,关注我们,解锁更多精彩内容!

使用biopython处理序列数据

一个只分享干货的

生信公众号


生物医学科研方法

正式成立丨国家卫生健康委医疗机构感染防控专家委员会

2020-12-30 1:47:24

生物医学科研方法

国家科技进步奖“最神秘获奖人”,任科技部副部长!

2020-12-30 2:05:29