database - 如何使用 Perl 从 NCBI 获取 FASTA 核苷酸格式的基因特征？

Question

我可以手动下载 FASTA 文件，如下所示:

>lcl|CR543861.1_gene_1...
ATGCTTTGGACA...
>lcl|CR543861.1_gene_2...
GTGCGACTAAAA...

通过单击“发送到”并选择“基因特征”，FASTA 核苷酸是 this page 上的唯一选项(这很好，因为这就是我想要的) .

使用这样的脚本:

#!/usr/bin/env perl
use strict;
use warnings;
use Bio::DB::EUtilities;

my $factory = Bio::DB::EUtilities->new(-eutil   => 'efetch',
                                       -db      => 'nucleotide',
                                       -id      => 'CR543861',
                                       -rettype => 'fasta');
my $file = 'CR543861.fasta';
$factory->get_Response(-file => $file);

我得到一个如下所示的文件:

>gi|49529273|emb|CR543861.1| Acinetobacter sp. ADP1 complete genome
GATATTTTATCCACA...

将整个基因组序列集中在一起。 如何获取第一个(手动下载)文件中的信息？

我看了其他几个帖子:

• how to download complete genome sequence in biopython entrez.esearch ( this answer seemed relevant )
• How can I download the entire GenBank file with just an accession number?

以及this section from EUtilities Cookbook .

我尝试获取并保存一个 GenBank 文件(因为我得到的 .gb 文件中的每个基因似乎都有单独的序列)，但是当我使用 Bio::SeqIO 处理它时，我只会得到 1 个大的顺序。

Answer 1

使用该登录号和返回类型，您将获得完整的基因组序列。如果您想获取单个基因序列，请指定您想要完整的 genbank 文件，然后解析出基因。这是一个例子:

#!/usr/bin/env perl

use 5.010;
use strict;
use warnings;
use Bio::SeqIO;
use Bio::DB::EUtilities;


my $factory = Bio::DB::EUtilities->new(-eutil   => 'efetch',
                                       -email   => 'foo@bar.com',
                                       -db      => 'nucleotide',
                                       -id      => 'CR543861',
                                       -rettype => 'gb');
my $file = 'CR543861.gb';
$factory->get_Response(-file => $file);

my @gene_features = grep { $_->primary_tag eq 'gene' } 
                    Bio::SeqIO->new(-file => $file)->next_seq->get_SeqFeatures;

for my $feat_object (@gene_features) {
    for my $tag ($feat_object->get_all_tags) {
        # open a filehandle here for writing each to a separate file
        say ">",$feat_object->get_tag_values($tag);
        say $feat_object->spliced_seq->seq;
        # close it!
    } 
}

这会将每个基因写入同一个文件(如果您重定向它，现在它只写入 STDOUT)但我指出了您可以在哪里进行小的更改以将它们写入单独的文件。解析 genbank 有时会有点棘手，因此阅读文档总是有帮助的，尤其是优秀的 Feature Annotation HOWTO .