/usr/bin/bp_fetch is in bioperl 1.6.924-3.
This file is owned by root:root, with mode 0o755.
The actual contents of the file can be viewed below.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295 296 297 298 299 300 301 302 303 304 305 306 307 308 309 310 311 312 313 314 315 316 317 318 319 320 321 322 323 324 | #!/usr/bin/perl
eval 'exec /usr/bin/perl -S $0 ${1+"$@"}'
if 0; # not running under some shell
=head1 NAME
bp_fetch.pl - fetches sequences from bioperl indexed databases
=head1 SYNOPSIS
bp_fetch.pl swiss:ROA1_HUMAN
bp_fetch.pl net::genbank:JX295726
bp_fetch.pl net::genpept:ROA1_HUMAN
bp_fetch.pl ace::myserver.somewhere.edu,21000:X56676
bp_fetch.pl -fmt GCG swiss:ROA1_HUMAN
=head1 DESCRIPTION
Fetches sequences using the DB access systems in Bioperl. The most
common use of this is to bp_fetch sequences from bioperl indices built
using bpindex.pl, or to fetch sequences from the NCBI website
The format for retrieving sequences is delibrately like the
GCG/EMBOSS format like the following:
db:name
with the potential of putting in a 'meta' database type, being
meta::db:name
The meta information can be one of three types
local - local indexed flat file database
net - networked http: based database
ace - ACeDB database
This information defaults to 'local' for database names with no meta
db information
=head1 OPTIONS
-fmt <format> - Output format
Fasta (default), EMBL, Raw, swiss or GCG
-acc - string is an accession number, not an
id.
options only for expert use
-dir <dir> - directory to find the index files
(overrides BIOPERL_INDEX environment varaible)
-type <type> - type of DBM file to open
(overrides BIOPERL_INDEX_TYPE environment variable)
=head1 ENVIRONMENT
bp_index and bp_fetch coordinate where the databases lie using the
environment variable BIOPERL_INDEX. This can be overridden using the
-dir option. The index type (SDBM or DB_File or another index file)
is controlled by the BIOPERL_INDEX_TYPE variable. This defaults to
SDBM_File
=head1 USING IT YOURSELF
bp_fetch is a wrapper around the bioperl modules which support
the Bio::DB::BioSeqI abstract interface. These include:
Author Code
James Gilbert - Fasta indexer, Abstract indexer
Aaron Mackay - GenBank and GenPept DB access
Ewan Birney - EMBL .dat indexer
Many people - SeqIO code
These modules can be used directly, which is far better than using
this script as a system call or a pipe to read from. Read the
source code for bp_fetch to see how it is used.
=head1 EXTENDING IT
bp_fetch uses a number of different modules to provide access to
databases. Any module which subscribes to the Bio::DB::BioSeqI
interface can be used here. For flat file indexers, this is
best done by extending Bio::Index::Abstract, as is done in
Bio::Index::EMBL and Bio::Index::Fasta. For access to other
databases you will need to roll your own interface.
For new output formats, you need to add a new SeqIO module. The
easiest thing is to look at Bio::SeqIO::Fasta and figure out
how to hack it for your own format (call it something different
obviously).
=head1 FEEDBACK
=head2 Mailing Lists
User feedback is an integral part of the evolution of this and other
Bioperl modules. Send your comments and suggestions preferably to
the Bioperl mailing list. Your participation is much appreciated.
bioperl-l@bioperl.org - General discussion
http://bioperl.org/wiki/Mailing_lists - About the mailing lists
=head2 Reporting Bugs
Report bugs to the Bioperl bug tracking system to help us keep track
of the bugs and their resolution. Bug reports can be submitted via the
web:
https://github.com/bioperl/bioperl-live/issues
=head1 AUTHOR
Ewan Birney E<lt>birney@ebi.ac.ukE<gt>
=cut
use strict;
use warnings;
use Getopt::Long;
#
# Dofus catcher for people who are trying this script without
# installing bioperl. In your own script, you can just go
#
# use Bio::Index::Fasta etc, rather than this
#
use Bio::SeqIO;
BEGIN {
eval {
require Bio::Index::Fasta;
require Bio::Index::EMBL;
require Bio::Index::GenBank;
require Bio::Index::Swissprot;
require Bio::Index::SwissPfam;
};
if ( $@ ) {
# one up from here is Bio directory - we hope!
push(@INC,"..");
eval {
require Bio::Index::Fasta;
require Bio::Index::EMBL;
require Bio::Index::GenBank;
require Bio::Index::Swissprot;
require Bio::Index::SwissPfam;
};
if ( $@ ) {
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");
exit(1);
} else {
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");
}
}
eval {
require Bio::DB::GenBank;
require Bio::DB::GenPept;
require Bio::DB::EMBL;
require Bio::DB::SwissProt;
};
if ( $@ ) {
if ( !exists $ENV{'BIOPERL_SAVVY'} ) {
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 supress this message\n\n");
}
}
}
#
# Start processing the command line
#
my $dir = $ENV{'BIOPERL_INDEX'};
my $type = $ENV{'BIOPERL_INDEX_TYPE'};
my $fmt = 'Fasta';
my $useacc = 0;
my $ret = GetOptions('d|dir=s' => \$dir,
'f|fmt=s' => \$fmt ,
't|type=s' => \$type ,
'acc!' => \$useacc);
#
# print pod documentation if we have no arguments
#
exec('perldoc',$0) unless @ARGV;
my ($isnet,$db,$dbobj,$id,$seq,$seqio,$out,$meta);
#
# Reset the type if needed
#
if( $type ) {
$Bio::Index::Abstract::USE_DBM_TYPE = $type;
}
#
# Build at run time the SeqIO output
#
if ( $fmt !~ /swisspfam|pfam/ ) {
$out = Bio::SeqIO->new(-fh => \*STDOUT , -format => $fmt);
}
#
# Main loop over remaining arguments
#
for my $arg ( @ARGV ) {
$_= $arg;
# strip out meta:: if there
if ( /^(\w+)::/ ) {
$meta = $1;
s/^(\w+):://;
} else {
$meta = 'local';
}
# parse to db:id
/^(\S+)\:(\S+)$/ || do { warn "$_ is not parsed as db:name\n"; next; };
($db,$id) = split/:/,$_,2;
#
# the eval block catches exceptions if they occur
# in the code in the block. The exception goes in $@
#
eval {
SWITCH : {
$_ = $meta;
/^net$/ && do {
if ( $db =~ /genbank/i ) {
$dbobj = Bio::DB::GenBank->new(-format => $fmt);
} elsif ( $db =~ /genpept/i ) {
$dbobj = Bio::DB::GenPept->new();
} elsif ( $db =~ /embl/i ) {
$dbobj = Bio::DB::EMBL->new();
} else {
die "Net database $db not available";
}
last SWITCH;
};
/^ace$/ && do {
# yank in Bio::DB::Ace at runtime
eval {
require Bio::DB::Ace;
};
if ( $@ ) {
die "Unable to load Bio::DB::Ace for ace::$db\n\n$@\n";
}
# db is server,port
my ($server,$port);
$db =~ /(\S+)\,(\d+)/ || die "$db is not server.name,port for acedb database";
$server = $1;
$port = $2;
# print STDERR "Connecting to $server,$port\n";
$dbobj = Bio::DB::Ace->new(-host => $server, -port => $port);
last SWITCH;
};
/^local$/ && do {
if ( !$dir ) {
die "\nNo directory specified for index\nDirectory must be specified by the environment varaible BIOPERL_INDEX or --dir option\ngo bp_index with no arguments for more help\n\n";
}
#
# $db gets re-blessed to the correct index when
# it is made from the abstract class. Cute eh?
#
$dbobj = Bio::Index::Abstract->new("$dir/$db");
last SWITCH;
};
die "Meta database $meta is not valid";
}
}; # end of eval to get db
if ( $@ ) {
warn("Database $db in $arg is not loadable. Skipping\n\nError $@");
next;
}
#
# We expect the databases to adhere to the BioSeqI
# the sequence index databases and the GenBank/GenPept do already
#
if ( $dbobj->isa("Bio::Index::SwissPfam") ) {
my $seq = $dbobj->fetch($id);
if ( $seq ) {
my $started;
while ( <$seq> ) {
last if ( /^\s+$/ );
print;
}
} else {
warn("Cannot find $id\n");
}
next;
}
if ( ! $dbobj->isa('Bio::DB::RandomAccessI') ) {
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");
}
eval {
if ( $useacc == 0 ) {
$seq = $dbobj->get_Seq_by_id($id);
} else {
$seq = $dbobj->get_Seq_by_acc($id);
}
};
if ( $@ ) {
warn("Sequence $id in Database $db in $arg is not loadable. Skipping.\n\nError $@");
next;
} elsif ( !defined $seq ) {
warn("Sequence $id in Database $db is not present\n");
next;
}
$out->write_seq($seq);
}
|