5 bp_fetch.pl - fetches sequences from bioperl indexed databases
9 bp_fetch.pl swiss:ROA1_HUMAN
11 bp_fetch.pl net::genbank:JX295726
13 bp_fetch.pl net::genpept:ROA1_HUMAN
15 bp_fetch.pl ace::myserver.somewhere.edu,21000:X56676
17 bp_fetch.pl -fmt GCG swiss:ROA1_HUMAN
21 Fetches sequences using the DB access systems in Bioperl. The most
22 common use of this is to bp_fetch sequences from bioperl indices built
23 using bpindex.pl, or to fetch sequences from the NCBI website
25 The format for retrieving sequences is delibrately like the
26 GCG/EMBOSS format like the following:
30 with the potential of putting in a 'meta' database type, being
34 The meta information can be one of three types
36 local - local indexed flat file database
37 net - networked http: based database
40 This information defaults to 'local' for database names with no meta
45 -fmt <format> - Output format
46 Fasta (default), EMBL, Raw, swiss or GCG
47 -acc - string is an accession number, not an
50 options only for expert use
52 -dir <dir> - directory to find the index files
53 (overrides BIOPERL_INDEX environment variable)
54 -type <type> - type of DBM file to open
55 (overrides BIOPERL_INDEX_TYPE environment variable)
59 bp_index and bp_fetch coordinate where the databases lie using the
60 environment variable BIOPERL_INDEX. This can be overridden using the
61 -dir option. The index type (SDBM or DB_File or another index file)
62 is controlled by the BIOPERL_INDEX_TYPE variable. This defaults to
65 =head1 USING IT YOURSELF
67 bp_fetch is a wrapper around the bioperl modules which support
68 the Bio::DB::BioSeqI abstract interface. These include:
72 James Gilbert - Fasta indexer, Abstract indexer
73 Aaron Mackay - GenBank and GenPept DB access
74 Ewan Birney - EMBL .dat indexer
75 Many people - SeqIO code
77 These modules can be used directly, which is far better than using
78 this script as a system call or a pipe to read from. Read the
79 source code for bp_fetch to see how it is used.
83 bp_fetch uses a number of different modules to provide access to
84 databases. Any module which subscribes to the Bio::DB::BioSeqI
85 interface can be used here. For flat file indexers, this is
86 best done by extending Bio::Index::Abstract, as is done in
87 Bio::Index::EMBL and Bio::Index::Fasta. For access to other
88 databases you will need to roll your own interface.
90 For new output formats, you need to add a new SeqIO module. The
91 easiest thing is to look at Bio::SeqIO::Fasta and figure out
92 how to hack it for your own format (call it something different
99 User feedback is an integral part of the evolution of this and other
100 Bioperl modules. Send your comments and suggestions preferably to
101 the Bioperl mailing list. Your participation is much appreciated.
103 bioperl-l@bioperl.org - General discussion
104 http://bioperl.org/wiki/Mailing_lists - About the mailing lists
106 =head2 Reporting Bugs
108 Report bugs to the Bioperl bug tracking system to help us keep track
109 of the bugs and their resolution. Bug reports can be submitted via the
112 https://github.com/bioperl/bioperl-live/issues
116 Ewan Birney E<lt>birney@ebi.ac.ukE<gt>
124 use Bio
::Index
::EMBL
;
125 use Bio
::Index
::Fasta
;
126 use Bio
::Index
::GenBank
;
127 use Bio
::Index
::SwissPfam
;
128 use Bio
::Index
::Swissprot
;
132 require Bio
::DB
::GenBank
;
133 require Bio
::DB
::GenPept
;
134 require Bio
::DB
::EMBL
;
135 require Bio
::DB
::SwissProt
;
137 if ($@
&& ! exists $ENV{'BIOPERL_SAVVY'}) {
138 print STDERR
("\nbp_fetch cannot find Bio::DB::GenBank and Bio::DB::EMBL modules\nThis does not effect local indexing\nset environment variable BIOPERL_SAVVY to suppress this message\n\n");
142 # Start processing the command line
145 my $dir = $ENV{'BIOPERL_INDEX'};
146 my $type = $ENV{'BIOPERL_INDEX_TYPE'};
149 my $ret = GetOptions
('d|dir=s' => \
$dir,
151 't|type=s' => \
$type ,
155 # print pod documentation if we have no arguments
158 exec('perldoc',$0) unless @ARGV;
160 my ($isnet,$db,$dbobj,$id,$seq,$seqio,$out,$meta);
163 # Reset the type if needed
167 $Bio::Index
::Abstract
::USE_DBM_TYPE
= $type;
171 # Build at run time the SeqIO output
173 if ( $fmt !~ /swisspfam|pfam/ ) {
174 $out = Bio
::SeqIO
->new(-fh
=> \
*STDOUT
, -format
=> $fmt);
178 # Main loop over remaining arguments
181 for my $arg ( @ARGV ) {
183 # strip out meta:: if there
193 /^(\S+)\:(\S+)$/ || do { warn "$_ is not parsed as db:name\n"; next; };
194 ($db,$id) = split/:/,$_,2;
196 # the eval block catches exceptions if they occur
197 # in the code in the block. The exception goes in $@
204 if ( $db =~ /genbank/i ) {
205 $dbobj = Bio
::DB
::GenBank
->new(-format
=> $fmt);
206 } elsif ( $db =~ /genpept/i ) {
207 $dbobj = Bio
::DB
::GenPept
->new();
208 } elsif ( $db =~ /embl/i ) {
209 $dbobj = Bio
::DB
::EMBL
->new();
211 die "Net database $db not available";
216 # yank in Bio::DB::Ace at runtime
218 require Bio
::DB
::Ace
;
221 die "Unable to load Bio::DB::Ace for ace::$db\n\n$@\n";
227 $db =~ /(\S+)\,(\d+)/ || die "$db is not server.name,port for acedb database";
230 # print STDERR "Connecting to $server,$port\n";
232 $dbobj = Bio
::DB
::Ace
->new(-host
=> $server, -port
=> $port);
237 die "\nNo directory specified for index\nDirectory must be specified by the environment variable BIOPERL_INDEX or --dir option\ngo bp_index with no arguments for more help\n\n";
241 # $db gets re-blessed to the correct index when
242 # it is made from the abstract class. Cute eh?
245 $dbobj = Bio
::Index
::Abstract
->new("$dir/$db");
248 die "Meta database $meta is not valid";
250 }; # end of eval to get db
252 warn("Database $db in $arg is not loadable. Skipping\n\nError $@");
257 # We expect the databases to adhere to the BioSeqI
258 # the sequence index databases and the GenBank/GenPept do already
260 if ( $dbobj->isa("Bio::Index::SwissPfam") ) {
261 my $seq = $dbobj->fetch($id);
269 warn("Cannot find $id\n");
273 if ( ! $dbobj->isa('Bio::DB::RandomAccessI') ) {
274 warn("$db in $arg does not inherit from Bio::DB::RandomAccessI, so is not expected to work under the DB guidlines. Going to try it anyway");
277 if ( $useacc == 0 ) {
278 $seq = $dbobj->get_Seq_by_id($id);
280 $seq = $dbobj->get_Seq_by_acc($id);
284 warn("Sequence $id in Database $db in $arg is not loadable. Skipping.\n\nError $@");
286 } elsif ( !defined $seq ) {
287 warn("Sequence $id in Database $db is not present\n");
290 $out->write_seq($seq);