从 GenBank 文件中输出基因位置

Posted 2023-02-23

【中文标题】从 GenBank 文件中输出基因位置

【英文标题】：Output gene positions from GenBank file

【发布时间】：2015-09-23 15:34:14

【问题描述】：

是否可以输出 CDS 特征的基因位置，还是我需要自己解析“位置”或“补充”字段？

例如，

seq = Sequence.read(genbank_fp, format='genbank')
for feature in seq.metadata['FEATURES']:
    if feature['type_'] == 'CDS':
        if 'location' in feature:
            print 'location = ', feature['location']
        elif 'complement' in feature:
            print 'location = ', feature['complement']
        else:
            raise ValueError('positions for gene %s not found' % feature['protein_id'])

会输出：

位置 =

位置 = 687..3158

对于 this 样本 GenBank 文件。

此功能在 BioPython 中是可能的（请参阅 this thread），我可以在其中输出已解析的位置（例如 start = 687，end = 3158）。

谢谢！

【参考方案1】：

对于示例，您可以使用以下代码仅获取该功能的 Sequence 对象：

# column index in positional metadata
col = feature['index_']
loc = seq.positional_metadata[col]
feature_seq = seq[loc]
# if the feature is on reverse strand
if feature['rc_']:
    feature_seq = feature_seq.reverse_complement()

注意：GenBank 解析器是在开发分支中新增的。

【讨论】：

感谢您的回复！但是我正在寻找基因的实际开始和结束位置（而不是基因序列本身）。为了扩展您的答案，我需要找到 loc 系列中第一个和最后一个 True 元素的索引：start_pos = loc[loc == True].index[0]; end_pos = loc[loc == True].index[-1]