bin/bp_fetch

   1 #!/usr/bin/perl
   2
   3 =head1 NAME
   4
   5 bp_fetch.pl - fetches sequences from bioperl indexed databases
   6
   7 =head1 SYNOPSIS
   8
   9   bp_fetch.pl swiss:ROA1_HUMAN
  10
  11   bp_fetch.pl net::genbank:JX295726
  12
  13   bp_fetch.pl net::genpept:ROA1_HUMAN
  14
  15   bp_fetch.pl ace::myserver.somewhere.edu,21000:X56676
  16
  17   bp_fetch.pl -fmt GCG swiss:ROA1_HUMAN
  18
  19 =head1 DESCRIPTION
  20
  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
  24
  25 The format for retrieving sequences is delibrately like the
  26 GCG/EMBOSS format like the following:
  27
  28   db:name
  29
  30 with the potential of putting in a 'meta' database type, being
  31
  32   meta::db:name
  33
  34 The meta information can be one of three types
  35
  36   local - local indexed flat file database
  37   net   - networked http: based database
  38   ace   - ACeDB database
  39
  40 This information defaults to 'local' for database names with no meta
  41 db information
  42
  43 =head1 OPTIONS
  44
  45   -fmt  <format> - Output format
  46                    Fasta (default), EMBL, Raw, swiss or GCG
  47   -acc           - string is an accession number, not an
  48                    id.
  49
  50 options only for expert use
  51
  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)
  56
  57 =head1 ENVIRONMENT
  58
  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
  64
  65 =head1 USING IT YOURSELF
  66
  67 bp_fetch is a wrapper around the bioperl modules which support
  68 the Bio::DB::BioSeqI abstract interface. These include:
  69
  70   Author          Code
  71
  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
  76
  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.
  80
  81 =head1 EXTENDING IT
  82
  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.
  89
  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).
  94
  95 =head1 FEEDBACK
  96
  97 =head2 Mailing Lists
  98
  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.
 102
 103   bioperl-l@bioperl.org                  - General discussion
 104   http://bioperl.org/wiki/Mailing_lists  - About the mailing lists
 105
 106 =head2 Reporting Bugs
 107
 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:
 111
 112   https://github.com/bioperl/bioperl-live/issues
 113
 114 =head1 AUTHOR
 115
 116 Ewan Birney E<lt>birney@ebi.ac.ukE<gt>
 117
 118 =cut
 119
 120 use strict;
 121 use warnings;
 122 use Getopt::Long;
 123
 124 #
 125 # Dofus catcher for people who are trying this script without
 126 # installing bioperl. In your own script, you can just go
 127 #
 128 # use Bio::Index::Fasta etc, rather than this
 129 #
 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;
 138   };
 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;
 148     };
 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");
 154     }
 155   }
 156
 157   eval {
 158     require Bio::DB::GenBank;
 159     require Bio::DB::GenPept;
 160     require Bio::DB::EMBL;
 161     require Bio::DB::SwissProt;
 162   };
 163
 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");
 167     }
 168   }
 169 }
 170
 171 #
 172 # Start processing the command line
 173 #
 174
 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);
 183
 184 #
 185 # print pod documentation if we have no arguments
 186 #
 187
 188 exec('perldoc',$0) unless @ARGV;
 189
 190 my ($isnet,$db,$dbobj,$id,$seq,$seqio,$out,$meta);
 191
 192 #
 193 # Reset the type if needed
 194 #
 195
 196 if( $type ) {
 197    $Bio::Index::Abstract::USE_DBM_TYPE = $type;
 198 }
 199
 200 #
 201 # Build at run time the SeqIO output
 202 #
 203 if ( $fmt !~ /swisspfam|pfam/ ) {
 204   $out = Bio::SeqIO->new(-fh => \*STDOUT , -format => $fmt);
 205 }
 206
 207 #
 208 # Main loop over remaining arguments
 209 #
 210
 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';
 219   }
 220
 221   # parse to db:id
 222
 223   /^(\S+)\:(\S+)$/ || do { warn "$_ is not parsed as db:name\n"; next; };
 224   ($db,$id) = split/:/,$_,2;
 225   #
 226   # the eval block catches exceptions if they occur
 227   # in the code in the block. The exception goes in $@
 228   #
 229
 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";
 242         }
 243         last SWITCH;
 244       };
 245       /^ace$/ && do {
 246         # yank in Bio::DB::Ace at runtime
 247         eval {
 248           require Bio::DB::Ace;
 249         };
 250         if ( $@ ) {
 251           die "Unable to load Bio::DB::Ace for ace::$db\n\n$@\n";
 252         }
 253
 254         # db is server,port
 255         my ($server,$port);
 256
 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";
 261
 262         $dbobj = Bio::DB::Ace->new(-host => $server, -port => $port);
 263         last SWITCH;
 264       };
 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";
 268         }
 269
 270         #
 271         # $db gets re-blessed to the correct index when
 272         # it is made from the abstract class. Cute eh?
 273         #
 274
 275         $dbobj = Bio::Index::Abstract->new("$dir/$db");
 276         last SWITCH;
 277       };
 278       die "Meta database $meta is not valid";
 279     }
 280     };                          # end of eval to get db
 281   if ( $@ ) {
 282     warn("Database $db in $arg is not loadable. Skipping\n\nError $@");
 283     next;
 284   }
 285
 286   #
 287   # We expect the databases to adhere to the BioSeqI
 288   # the sequence index databases and the GenBank/GenPept do already
 289   #
 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;
 297       }
 298     } else {
 299       warn("Cannot find $id\n");
 300     }
 301     next;
 302   }
 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");
 305   }
 306   eval {
 307     if ( $useacc == 0 ) {
 308       $seq = $dbobj->get_Seq_by_id($id);
 309     } else {
 310       $seq = $dbobj->get_Seq_by_acc($id);
 311     }
 312   };
 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;
 319   }
 320   $out->write_seq($seq);
 321 }