Fetch Many Files With Accession Number, Output File Format Is Coding Sequences In Fasta Per

Question

[Shihai Feng]

Hi all,

How do I fetch many many files according to accession numbers from Genbank. The output I would like is one file per accession numbers. The format of output file is Coding Sequences in Fasta. Any suggestion are highly appreciated.

Thanks, Shihai

Answer 1

[fransua]

perhaps you can use biomart at ensembl: http://www.ensembl.org/biomart/martview/ or directly from biomart: http://www.biomart.org/

Answer 2

[Fabian Bull]

If you are able to code have a look at: using Bioperl to retrieve multiple sequences

Just change the query variable accordingly.

Answer 3

[Woa]

You can download the nucleic acid sequences in the GENBANK format using the following piece of code:Run through a for loop for many acc. nos. Then one by one parse the files to get the CDS.

use Bio::DB::GenBank;
use Bio::SeqIO;
use strict;

my $acc='NM_001003933';#example entry
my $format='genbank';
my $file_name=$acc.'.gbk';

my $seqout = new Bio::SeqIO( -file => ">$file_name", -format=>$format);
    my $getseq = new Bio::DB::GenBank;
    my $seq = $getseq->get_Seq_by_acc($acc);
    $seqout->write_seq($seq);

This piece of code can parse the GENBANK files to get the CDS in fasta format.However where CDS is not availble like in case of pseudogenes or rRNA etc this can crash. The scalar(@cds_features) will return zero then. Use it if you have non coding sequences

    use strict;
    use warnings;
    use Bio::SeqIO;
    my $gb_file="NM_001003933.gbk";
    my @cds_features = grep { $_->primary_tag eq 'CDS' } Bio::SeqIO->new(-file => $gb_file)->next_seq->get_SeqFeatures;

    my @tags=qw/protein_id product translation/;
    my ($feat_object)=@cds_features;
    my ($id,$product,$translation)=map{$feat_object->get_tag_values($_);}@tags;

print ">",$id," ",$product,"\n",$translation,"\n";

[Shihai Feng]

Thanks, the only problem is this give me proteins and what I need is DNA. Since at Genbank web interface I can get a file looks like:

lcl|CP003084.1cdsidAER04537.1 [gene=TIIST4400005] [protein=tyrosine] [proteinid=AER04537.1] [location=complement(250..1506)] GTGTCGGATCTGTTGGACGACCTGGAGTGGAGGGGCCTCGTCGCTGACTCCACCGACCGTGAGGCGTT lcl|CP003084.1cdsidAER04538.1 [gene=TIIST4400010] [protein=arginine repressor] [proteinid=AER04538.1] [location=complement(1590..2138)] ATGAGTGAGGGGCAGTGTGCACCGACCTCACGTACGGCGCGGCATGCCCGGATTCGCGTCTTAATTGA . . . How can I it from efetch directly in Genbank?