使用biopython操作gff文件

Question

我有一个GFF文件，它是一个选项卡限制的9列文件。我的Gff文件如下所示：

chr1    GenBank region  1   2821361 .   +   1   ID=CP000253.1
chr1    S-MART  utr5    313 516     .   +   .   ID=CP000253.1|+313..516
chr1    GenBank gene    517 1878    .   +   1   ID=SAOUHSC_00001

。。。。。。。。。等等。在

问题陈述：

现在，我想合并满足条件的行。条件是第i行的第5列值应等于i+1行的第4列减去1。在

所以最后的结果应该是

^{pr2}$

为此，我编写了以下程序：

from BCBio import GFF
from Bio.SeqFeature import SeqFeature, FeatureLocation

in_file = "infile.gff"
out_file = "outfile.gff"

limit_info = dict(
        gff_type = ['CDS','exon','gene','mRNA','operon','rRNA','tRNA','utr3','utr5'])
new_qualifiers = {"source": "prediction","ID": "CP000253.1"}
new_sub_qualifiers = {"source": "prediction"}
new_top_feature = SeqFeature(FeatureLocation(0, 2821361), type="genomeRegion", strand=1,
                         qualifiers=new_qualifiers)
i=0

in_handle = open(in_file)
for rec in GFF.parse(in_handle, limit_info=limit_info):
    for i in range(10):
        if rec.features[i].location.end == rec.features[i+1].location.start :
            # print rec.features[i]
            new_top_feature.sub_features[i] =     
[SeqFeature(FeatureLocation(rec.features[i].location.start ,  
rec.features[i+1].location.end ,strand=rec.features[i].strand),  
type="Transcription_unit",  qualifiers=new_sub_qualifiers)]             

in_handle.close()

rec.features = [new_top_feature]

with open(out_file, "w") as out_handle:
    GFF.write([rec], out_handle)

我得到以下错误：

/usr/lib/python2.7/dist-packages/Bio/SeqFeature.py:171: BiopythonDeprecationWarning: Rather using f.sub_features, f.location should be a CompoundFeatureLocation
  BiopythonDeprecationWarning)
Traceback (most recent call last):
  File "/home/nkumar/workplacekepler/random/src/limit.py", line 26, in <module>
    new_top_feature.sub_features[i] = [SeqFeature(FeatureLocation(rec.features[i].location.start , rec.features[i+1].location.end ,strand=rec.features[i].strand), type="Transcription_unit",  qualifiers=new_sub_qualifiers)]
IndexError: list assignment index out of range

即使我不能计算出什么是错误的？在

in_handle = open(in_file)
for rec in GFF.parse(in_handle, limit_info=limit_info):
    for i in range(10):        
        if rec.features[i].location.end == rec.features[i+1].location.start :
            print 1          
        else:
            print rec.features[i]            
in_handle.close()

这一个完美的工作和打印所有的功能。在

Answer 1

您将新的“顶部”功能定义为：

type: genomeRegion
location: [0:2821361](+)
qualifiers: 
    Key: ID, Value: CP000253.1
    Key: source, Value: prediction

但它没有副功能

^{pr2}$

因此，new_top_feature.sub_features是一个空列表。不能直接分配给空列表：

>>> a = []
>>> a[0] = 3
Traceback (most recent call last):
  File "<input>", line 1, in <module>
IndexError: list assignment index out of range

这就是你要做的

new_top_feature.sub_features[i] =  .....

要将数据添加到此列表中，应使用append而不是索引。如果您需要按任意顺序在给定位置填充列表，您可以创建一个用零填充的足够大小的列表，然后在这些位置出现时将值分配给这些位置。在