从 genbank 文件中按特征提取 DNA 序列

Question

我有由多个带注释的重叠群组成的 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]

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

Answer 1

好的，没关系 - 更多的谷歌搜索和实验让我找到了更好的解决方案。不确定这是否是最好的解决方案，如果有人有更好的想法，我会留下这个问题并且不予回答。我为获得 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()

不是最漂亮的解决方案，但似乎可行...

Answer 2

我使用位置对象的 "extract" 方法。

dna_seq = feature.location.extract(record)

这很好用。

从 genbank 文件中按特征提取 DNA 序列

Pull dna sequence by feature from genbank file

biopython

python-2.7