Biopython
Biofeatures
Biofeatures是BioIO中的一部分。
整个程序使用之前需要输入from Bio import SeqIO
from Bio import SeqIO
#先建立一个SeqRecord对象SeqIO.read(序列文件名称,序列形式)
record = SeqIO.read("AB042240.3.gb","genbank")
#获取文件中全部features,即所有的gene,CDS,tRNA等
features = record.features
#test
for f in features:
print(f)
输出的结果如下图:
图片中左侧页面为输出,右侧页面为输入文件形式。
可以看到每一个feature对象都有这样几个部分:
- type
- location
- qualifiers
根据上图我将阐述个人对这些参数的理解
1.type
type应该就是genebank文件中在FEATURES这一大块中,写在最前面的一部分”小标题”,即下图红框中内容。这些小标题应该分别表示了后面一段记录的内容是什么。
2.location
可以看到这一部分表示了对应的序列位置信息。
其形式为:
- 如果有多个片段(有编码序列和非编码序列穿插存在,对于CDS,只需要编码序列,其位置信息就会是一段一段的)
形式为:
join{[67271:67385](-), [88797:89018](-), [88229:88257](-)
}
{[片段1起始位点:片段1结束位点](“-”为complement,反向互补序列,“+”为正常序列)} - 只有一个片段
形式为:
[134213:134495](+)
[起始位点:结束位点](“-”为complement,反向互补序列,“+”为正常序列)
具体探究如何对于这些信息进行提取,如下。
首先,书本给出的定义是:
SeqFeature 对象主要用于描述相对于父序列中的位置(region)信息。Region用location对象表示,通常是两个position间的范围。为了区分location和position,我们定义如下:
position
– 表示位于序列中的单一位置, 可以是精确的也可以是不确定的位置(如5, 20, <100 和 >200 )。
location
– 介于两个positions间的区域。比如5..20 (5到20)。
ref:https://biopython-cn.readthedocs.io/zh_CN/latest/cn/chr04.html#positionslocations
但是,对于如何提取多片段中每一小段的始末位点,在上述链接中博主没有找到很好的解决。因此博主使自己编码对于这一部分进行了完善。
#测试location
loc_of_2 = features[1].location
print(loc_of_2)
print(loc_of_2.start)
print(loc_of_2.end)
print(str(loc_of_2))
str_loc_of_2 = str(loc_of_2)
if "join"in str_loc_of_2:
splits = str_loc_of_2.split(",")
for specie in splits:
k=0 #记录为起始位点还是结束位点,0为起始位点,1为结束位点
start=""
end = ""
is_comp = False #定义是否为反向互补序列
if "-" in specie:
is_comp = True
for char in specie:
if char.isdigit() and k ==0 :
start = start+char
if char == ":" and k==0:
k=1
if char.isdigit() and k ==1:
end = end+char
print(start,end,is_comp)
结果:
join{[67271:67385](-), [88797:89018](-), [88229:88257](-)}
67271
89018
join{[67271:67385](-), [88797:89018](-), [88229:88257](-)}
67271 67385 True
88797 89018 True
88229 88257 True
可以看到直接使用loc_of_2.start和loc_of_2.end只能得到这一大片的始末位点。
因此需要将得到的信息转换为string,通过对字符串匹配进行各个片段信息的提取。
3.qualifiers
通过如下代码尝试:
#测试qualifiers
qual_of_2 = features[1].qualifiers
print(qual_of_2)
print(qual_of_2['gene'])
#print(qual_of_2.gene) 报错
print(qual_of_2['gene'][0])
print(type(qual_of_2['gene'][0]))
得到结果:
OrderedDict([('gene', ['rps12']), ('codon_start', ['1']), ('transl_table', ['11']), ('product', ['ribosomal protein S12']), ('protein_id', ['BAB47079.1']), ('translation', ['MPTVKQLIRNARQPIRNARKTAALKGCPQRRGTCARVYKKPNSALRKVARVRLTSGFEITAYIPGIGHNLQEHSVVLVRGGRVKDLPGVRYRIIRGTLDAVAVKNRQQGRSRYGVKKPKK'])])
['rps12']
rps12
<class 'str'>
可以看到,qualifiers方法可以得到一个ordereddict格式的结果。具体使用方法与普通字典类似,OrderdedDict([(键1,值1),(键2,值2)])
以上内容可直接运行的代码见:https://download.csdn.net/download/wwwddd777/19360082?spm=1001.2014.3001.5503