0

基本上，GenBank 文件包含基因条目（由“基因”宣布，然后是其相应的“CDS”条目（每个基因只有一个），就像我在下面展示的两个一样。我想在制表符分隔中获取 locus_tag 与产品两列文件。“基因”和“CDS”总是前后有空格。

上一个问题建议了一个脚本。

问题在于，似乎因为“产品”有时在其名称中包含“/”字符，因此它与该脚本有冲突，据我所知，它使用“/”作为字段分隔符将信息存储在大批？

我想解决这个问题，要么修改这个脚本，要么构建另一个。

perl -nE'
  BEGIN{ ($/, $") = ("CDS", "\t") }
  say "@r[0,1]" if @r= m!/(?:locus_tag|product)="(.+?)"!g and @r>1
' file


 gene            complement(8972..9094)
                 /locus_tag="HAPS_0004"
                 /db_xref="GeneID:7278619"
 CDS             complement(8972..9094)
                 /locus_tag="HAPS_0004"
                 /codon_start=1
                 /transl_table=11
                 /product="hypothetical protein"
                 /protein_id="YP_002474657.1"
                 /db_xref="GI:219870282"
                 /db_xref="GeneID:7278619"
                 /translation="MYYKALAHFLPTLSTMQNILSKSPLSLDFRLLFLAFIDKR"
 gene            68..637
                 /locus_tag="HPNK_00040"
 CDS             68..637
                 /locus_tag="HPNK_00040"
                 /codon_start=1
                 /transl_table=11
                 /product="NinG recombination protein/bacteriophage lambda
                 NinG family protein"
                 /protein_id="CRESA:HPNK_00040"
                 /translation="MIKPKVKKRKCKCCGGEFKSADSFRKWCSAECGVKLAKIAQEKA
                 RQKAIEKRNREERAKIKATRERLKSRSEWLKDAQAIFNEYIRLRDKDEPCISCRRFHQ
                 GQYHAGHYRTVKAMPELRFNEDNVHKQCSACNNHLSGNITEYRINLVRKIGAERVEAL
                 ESYHPPVKWSVEDCKEIIKTYRAKIKELK"
4

2 回答 2

2

由于您的示例 GenBank 文件不完整，我上网查找了可以用于示例的示例文件，我找到了这个文件。

使用此代码和Bio::GenBankParser模块，它被解析为猜测您所追求的结构的哪些部分。在这种情况下，“特征”同时包含一个locus_tag字段和一个product字段。

use strict;
use warnings;
use feature 'say';
use Bio::GenBankParser;

my $file = shift;
my $parser = Bio::GenBankParser->new( file => $file );
while ( my $seq = $parser->next_seq ) {
    my $feat = $seq->{'FEATURES'};
    for my $f (@$feat) {
        my $tag = $f->{'feature'}{'locus_tag'};
        my $prod = $f->{'feature'}{'product'};
        if (defined $tag and defined $prod) {
            say join "\t", $tag, $prod;
        }
    }
}

用法：

perl script.pl input.txt > output.txt

输出：

MG_001  DNA polymerase III, beta subunit
MG_470  CobQ/CobB/MinD/ParA nucleotide binding domain-containing protein

同一输入的单线输出将是：

MG_001  DNA polymerase III, beta subunit
MG_470  CobQ/CobB/MinD/ParA nucleotide binding
                     domain-containing protein

当然，假设您将/s修饰符添加到正则表达式以考虑多行条目（leeduhem在评论中指出）：

m!/(?:locus_tag|product)="(.+?)"!sg
#                                ^---- this
于 2014-02-27T13:06:15.127 回答
1

阅读了您重复的问题http://www.biostars.org/p/94164/（请不要像这样重复发布），这是一个最小的 Biopython 答案：

import sys
from Bio import SeqIO
filename = sys.argv[1] # Takes first command line argument input filename
for record in SeqIO.parse(filename, "genbank"):
    for feature in record.features:
        if feature.type == "CDS":
            locus_tag = feature.qualifiers.get("locus_tag", ["???"])[0]
            product = feature.qualifiers.get("product", ["???"])[0]
            print("%s\t%s" % (locus_tag, product))

通过较小的更改，您可以将其写入文件。

于 2014-02-28T16:16:31.730 回答