从 genbank 文件中按特征提取 dna 序列答案

【问题标题】：Pull dna sequence by feature from genbank file从 genbank 文件中按特征提取 dna 序列
【发布时间】：2015-10-26 12:24:41
【问题描述】：

我有由多个带注释的重叠群组成的 genbank 文件。我想做的是将其分离到一个数据库中，该数据库具有包含每个“CDS”特征的单个基因记录，以及它的 dna 和氨基酸序列。到目前为止，这工作得很好：

for record in SeqIO.parse(open_file, 'gb'):
   for feature in record.features:
        if feature.type == 'CDS':
            gene_record = {
                'locus_tag':feature.qualifiers['locus_tag'][0],
                'translation':feature.qualifiers['translation'][0],
                }

我遇到的麻烦是获取 dna 序列。 genbank 文件的格式如下：

FEATURES             Location/Qualifiers
     source          1..29869
                     /organism="Arthrobacter"
                     /mol_type="genomic DNA"
                     /strain="strain_name"
     gene            complement(4..462)
                     /locus_tag="ArthroDRAFT_00001"
     CDS             complement(4..462)
                     /locus_tag="ArthroDRAFT_00001"
                     /product="hypothetical protein"
                     /translation="LSTGKELLNYQSALNDIHDEFSRAQQSDAGVSHLSVAKITEKLS
                     YLKATALQMDDLFSVLRKQGVSLRSTGLADWASVPTIQDDKEEGKTEPSLAKKEISSR
                     TTSKPNKIDFPKFEYPDHGQPTNKIRVGTILDTFSESAFSYEWINVALQD"
     gene            complement(1126..1842)
                     /locus_tag="ArthroDRAFT_00002"
     CDS             complement(1126..1842)
                     /locus_tag="ArthroDRAFT_00002"
                     /product="hypothetical protein"
                     /translation="VPRAFIYGSCVGGDTANVFPSDWDRPTYVARQSIISAAFGPTSV
                     EGDIELTSAFQRSMLEGDIEATAFPRLRQELPTHDVLILDIVDERLGVYELAPGKYLT
                     RSMELISSKLIGKQPVTPRLIEFGSDEHYGLWTRSVDMLVDVVKHGGIPVFALLPPWS
                     EKSIQGEDLTWHSVSVDLMNNKYARYNEYLVQSEFTVVTVPDEEALGDAEHKWGLAPF
                     HYTESVYESLRDQILVGVSS"
...etc

ORIGIN
        1 ccctcaatcc tgaagagcca cattgatcca ttcgtatgag aatgcagatt ccgagaatgt
       61 atcaagaatt gttcctacac gtattttgtt tgtcggctgg ccgtggtccg gatactcaaa
      121 ttttggaaag tcaatcttgt ...

所以，我有每个特征的位置，但我需要以某种方式解析 dna 序列以提取正确的 dna 序列，而我有点不知所措。 feature.location 给了我看起来有用的输出：

[3:462](-)
[1125:1842](-)
[2159:3755](-)
[5190:5532](-)
[6226:6493](+)

但我知道使用它的唯一方法是使用正则表达式解析它，然后执行record.seq[start_num:finish_num]。然后，如果它在 - 链上，则通过 reverse_complement 运行任何内容。

这似乎太令人费解了，我认为 biopython 必须有一个更有效的方法，但我似乎找不到它。我注意到feature.location 似乎已经从其中减去了 1 以说明 python 的 0 索引，所以这个数字必须是可用的......对吧？

编辑：SeqFeature.location 似乎是我想要使用的，但不知道如何使用：http://biopython.org/DIST/docs/api/Bio.SeqFeature.FeatureLocation-class.html

EDIT2：看起来我可以做到

record.seq[feature.location.start:feature.location.end]

虽然这还不适用于负链的东西......

【问题讨论】：

标签： python-2.7 biopython

【解决方案1】：

我使用位置对象的“提取”方法。

dna_seq = feature.location.extract(record)

这很好用。

【讨论】：

几个月后我发现了这个...忘了回来回答这个问题 - 谢谢！

【解决方案2】：

好吧，没关系 - 更多的谷歌搜索和实验让我找到了更好的解决方案。不确定这是否是最好的解决方案，如果有人有更好的想法，我会留下问题并没有回答。我为获取 DNA 序列所做的工作如下：

if feature.location.strand == 1:
    dna_seq = record.seq[
        feature.location.start:feature.location.end
    ]
elif feature.location.strand == -1:
    dna_seq = record.seq[
        feature.location.start:feature.location.end
    ].reverse_complement()

不是最漂亮的解决方案，但它似乎工作......

【讨论】：