maint: restructure to use Dist::Zilla
[bioperl-live.git] / bin / bp_fetch
blob86a5bbc9a4f5ada0b89a46098a4582331fdfab67
1 #!/usr/bin/perl
3 =head1 NAME
5 bp_fetch.pl - fetches sequences from bioperl indexed databases
7 =head1 SYNOPSIS
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
19 =head1 DESCRIPTION
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:
28 db:name
30 with the potential of putting in a 'meta' database type, being
32 meta::db:name
34 The meta information can be one of three types
36 local - local indexed flat file database
37 net - networked http: based database
38 ace - ACeDB database
40 This information defaults to 'local' for database names with no meta
41 db information
43 =head1 OPTIONS
45 -fmt <format> - Output format
46 Fasta (default), EMBL, Raw, swiss or GCG
47 -acc - string is an accession number, not an
48 id.
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)
57 =head1 ENVIRONMENT
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
63 SDBM_File
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:
70 Author Code
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.
81 =head1 EXTENDING IT
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
93 obviously).
95 =head1 FEEDBACK
97 =head2 Mailing Lists
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
110 web:
112 https://github.com/bioperl/bioperl-live/issues
114 =head1 AUTHOR
116 Ewan Birney E<lt>birney@ebi.ac.ukE<gt>
118 =cut
120 use strict;
121 use warnings;
122 use Getopt::Long;
125 # Dofus catcher for people who are trying this script without
126 # installing bioperl. In your own script, you can just go
128 # use Bio::Index::Fasta etc, rather than this
130 use Bio::SeqIO;
131 BEGIN {
132 eval {
133 require Bio::Index::Fasta;
134 require Bio::Index::EMBL;
135 require Bio::Index::GenBank;
136 require Bio::Index::Swissprot;
137 require Bio::Index::SwissPfam;
139 if ( $@ ) {
140 # one up from here is Bio directory - we hope!
141 push(@INC,"..");
142 eval {
143 require Bio::Index::Fasta;
144 require Bio::Index::EMBL;
145 require Bio::Index::GenBank;
146 require Bio::Index::Swissprot;
147 require Bio::Index::SwissPfam;
149 if ( $@ ) {
150 print STDERR ("\nbp_index cannot find Bio::Index::Fasta and Bio::Index::EMBL\nbp_index needs to have bioperl installed for it to run.\nBioperl is very easy to install\nSee http://bio.perl.org for more information\n\n");
151 exit(1);
152 } else {
153 print STDERR ("\nYou are running bp_index.pl without installing bioperl.\nYou have done it from bioperl/scripts, and so we can find the necessary information\nbut it is much better to install bioperl\n\nPlease read the README in the bioperl distribution\n\n");
157 eval {
158 require Bio::DB::GenBank;
159 require Bio::DB::GenPept;
160 require Bio::DB::EMBL;
161 require Bio::DB::SwissProt;
164 if ( $@ ) {
165 if ( !exists $ENV{'BIOPERL_SAVVY'} ) {
166 print STDERR ("\nbp_fetch cannot find Bio::DB::GenBank and Bio::DB::EMBL modules\nThis is most likely because LWP has not been installed\nThis does not effect local indexing\nset environment variable BIOPERL_SAVVY to suppress this message\n\n");
172 # Start processing the command line
175 my $dir = $ENV{'BIOPERL_INDEX'};
176 my $type = $ENV{'BIOPERL_INDEX_TYPE'};
177 my $fmt = 'Fasta';
178 my $useacc = 0;
179 my $ret = GetOptions('d|dir=s' => \$dir,
180 'f|fmt=s' => \$fmt ,
181 't|type=s' => \$type ,
182 'acc!' => \$useacc);
185 # print pod documentation if we have no arguments
188 exec('perldoc',$0) unless @ARGV;
190 my ($isnet,$db,$dbobj,$id,$seq,$seqio,$out,$meta);
193 # Reset the type if needed
196 if( $type ) {
197 $Bio::Index::Abstract::USE_DBM_TYPE = $type;
201 # Build at run time the SeqIO output
203 if ( $fmt !~ /swisspfam|pfam/ ) {
204 $out = Bio::SeqIO->new(-fh => \*STDOUT , -format => $fmt);
208 # Main loop over remaining arguments
211 for my $arg ( @ARGV ) {
212 $_= $arg;
213 # strip out meta:: if there
214 if ( /^(\w+)::/ ) {
215 $meta = $1;
216 s/^(\w+):://;
217 } else {
218 $meta = 'local';
221 # parse to db:id
223 /^(\S+)\:(\S+)$/ || do { warn "$_ is not parsed as db:name\n"; next; };
224 ($db,$id) = split/:/,$_,2;
226 # the eval block catches exceptions if they occur
227 # in the code in the block. The exception goes in $@
230 eval {
231 SWITCH : {
232 $_ = $meta;
233 /^net$/ && do {
234 if ( $db =~ /genbank/i ) {
235 $dbobj = Bio::DB::GenBank->new(-format => $fmt);
236 } elsif ( $db =~ /genpept/i ) {
237 $dbobj = Bio::DB::GenPept->new();
238 } elsif ( $db =~ /embl/i ) {
239 $dbobj = Bio::DB::EMBL->new();
240 } else {
241 die "Net database $db not available";
243 last SWITCH;
245 /^ace$/ && do {
246 # yank in Bio::DB::Ace at runtime
247 eval {
248 require Bio::DB::Ace;
250 if ( $@ ) {
251 die "Unable to load Bio::DB::Ace for ace::$db\n\n$@\n";
254 # db is server,port
255 my ($server,$port);
257 $db =~ /(\S+)\,(\d+)/ || die "$db is not server.name,port for acedb database";
258 $server = $1;
259 $port = $2;
260 # print STDERR "Connecting to $server,$port\n";
262 $dbobj = Bio::DB::Ace->new(-host => $server, -port => $port);
263 last SWITCH;
265 /^local$/ && do {
266 if ( !$dir ) {
267 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";
271 # $db gets re-blessed to the correct index when
272 # it is made from the abstract class. Cute eh?
275 $dbobj = Bio::Index::Abstract->new("$dir/$db");
276 last SWITCH;
278 die "Meta database $meta is not valid";
280 }; # end of eval to get db
281 if ( $@ ) {
282 warn("Database $db in $arg is not loadable. Skipping\n\nError $@");
283 next;
287 # We expect the databases to adhere to the BioSeqI
288 # the sequence index databases and the GenBank/GenPept do already
290 if ( $dbobj->isa("Bio::Index::SwissPfam") ) {
291 my $seq = $dbobj->fetch($id);
292 if ( $seq ) {
293 my $started;
294 while ( <$seq> ) {
295 last if ( /^\s+$/ );
296 print;
298 } else {
299 warn("Cannot find $id\n");
301 next;
303 if ( ! $dbobj->isa('Bio::DB::RandomAccessI') ) {
304 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");
306 eval {
307 if ( $useacc == 0 ) {
308 $seq = $dbobj->get_Seq_by_id($id);
309 } else {
310 $seq = $dbobj->get_Seq_by_acc($id);
313 if ( $@ ) {
314 warn("Sequence $id in Database $db in $arg is not loadable. Skipping.\n\nError $@");
315 next;
316 } elsif ( !defined $seq ) {
317 warn("Sequence $id in Database $db is not present\n");
318 next;
320 $out->write_seq($seq);