/usr/bin/bp_seqret is in bioperl 1.7.1-2.
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 | #!/usr/bin/perl
# -*-Perl-*- mode (for emacs)
=head1 NAME
bp_seqret - bioperl implementation of sequence fetch from local db (like EMBOSS seqret)
=head1 USAGE
bp_seqret [-f/--format outputformat] [-o/--out/--outfile outfile] [-d/--db dbname] [-i/--id/-s/--seqname seqname1]
Example usage:
bp_seqret -f fasta -db db.fa -i seq1 -i seq2 > output.fas
bp_seqret db.fa:seq1 output.fas
bp_seqret db.fa:seq1 -o output.fas
bp_seqret -db db.fa -o output.fas seq1 seq2 seq3
bp_seqret -db db.fa seq1 seq2 seq3 output.fas
bp_seqret -db db.fa seq1 seq2 seq3 - > output.fas
The DB is expected to be a Fasta formatted sequence file with multiple
sequences.
Output format is Fasta by default.
If no output filename is provided then output is written to STDOUT.
Providing '-' as the output filename will accomplish the same thing.
=head1 AUTHOR
Jason Stajich jason_AT_bioperl-dot-org
=cut
use strict;
use warnings;
use Bio::DB::Fasta;
use Bio::SeqIO;
use Getopt::Long;
my $dbname;
my @names;
my $format = 'fasta';
my $outfile;
my ($start,$end);
GetOptions(
'f|format:s' => \$format,
'o|out|outfile:s' => \$outfile,
's|sbegin|begin|start:s' => \$start,
'e|send|end|stop:s' => \$end,
'd|db|dbname:s' => \$dbname,
'i|id|seqname:s' => \@names);
if( ! $dbname ) {
die "need a dbname\n" unless @ARGV;
$dbname = shift @ARGV;
if( $dbname =~ s/^([^:]+):// ) {
push @names, $dbname;
$dbname = $1;
}
}
my $db = Bio::DB::Fasta->new($dbname, -glob => "*.{fa,fas,fsa,fasta,pep,aa,seq,cds,peps}");
if( ! $outfile ) {
$outfile = pop @ARGV;
}
my $out;
if( $outfile ) {
$out = Bio::SeqIO->new(-format => $format,
-file => ">$outfile");
} else {
$out = Bio::SeqIO->new(-format => $format);
}
for my $nm ( @names ) {
my $seq;
if( $start || $end ) {
$seq = $db->seq($nm, $start => $end);
} else {
$seq = $db->seq($nm);
}
if( $seq ) {
my ($id,$desc) = split(/\s+/,$db->header($nm),2);
if( $start && $end ) {
$id = sprintf("%s_%d-%d",$id,$start || 0,$end || 0);
}
$out->write_seq(Bio::PrimarySeq->new(-display_id => $id,
-description => $desc,
-seq => $seq));
} else {
warn("$nm not found\n");
}
}
|