Ncbi下载的文件很多是注释完的Genebank格式文件,如果只有一个菌株的全基因组注释还能人工通过查删改,如果是一个文件含有大量的菌株,就要自动化解析了。
如果只有一条记录就用Bio.SeqIO.read()解析 如果有多条记录就用Bio.SeqIO.parse()遍历
Bio.SeqIO.read()返回的对象是SeqRecord数据结构
Bio.SeqIO.parse()返回的对象是GenBankIterator数据结构,遍历返回的是SeqRecord数据结构,也就是Genbank文件的元结构,是Bio模块赋予的
.seq返回Seq对象,也就是碱基序列,用bytes()函数可以得到测序碱基
Bio.SeqIO.write(Bio.SeqIO.parse(gbk_file, 'genbank'), "out_fasta.fasta", "fasta")
复制代码
最简单的把genebank格式转化成fasta格式是使用Bio.SeqIO.write()方法
但是为了了解一下SeqRecord结构
0bb0836ae2f6583b27b79548177570f.png
#!/usr/bin/env python
#usages: python genebank2fasta.py input_gbk_file > output_fasta_file
import Bio
from Bio import GenBank
from Bio import SeqIO
import sys
def abstract_fastas_from_gbk(gbk_file): #gbk文件生成一个key为note信息,value为Seq结构的字典
abstract_fasta = {}
for record in SeqIO.parse(gbk_file, 'genbank'):
for feature in record.features:
if feature.type == 'source' and 'note' in feature.qualifiers:
locus = "|".join(feature.qualifiers['note'])
abstract_fasta[locus] = record.seq
return abstract_fasta
def ouput_dict_as_txt(input_dict):
for key,value in input_dict.items():
print(">"+key)
print(value)
if __name__ == '__main__':
gbk_file = sys.argv[1]
abstract_fasta_dict = abstract_fastas_from_gbk(gbk_file)
ouput_dict_as_txt(abstract_fasta_dict)
复制代码
以上就是一个从Gnebank转换成muti-fasta的代码,每个SeqRecord都有很多feature
c0f5d2b2845f71223692220fab58335.png
也能通过修改以下代码把注释文件中的gene甚至翻译的氨基酸序列提取出来
for record in SeqIO.parse(gbk_file, 'genbank'):
for feature in record.features:
if feature.type == 'source' and 'note' in feature.qualifiers:
names.append(feature.qualifiers['note'][0])
if feature.type == 'CDS' and 'gene' in feature.qualifiers:
names.append(feature.qualifiers['gene'][0])
复制代码
所以genebank文件在Bio模块中以GenBankIterator》》SeqRecord》》SeqFeature》》(FeatureLocation,OrderedDict……)等对象依次组成
作者:尖尖家的大宝贝
链接:https://juejin.cn/post/7052135116508364830
来源:稀土掘金
著作权归作者所有。商业转载请联系作者获得授权,非商业转载请注明出处。