← Back to branch summary

~ubuntu-branches/ubuntu/gutsy/bioperl/gutsy

« back to all changes in this revision

Viewing changes to scripts/blast_fetch_local.pl

  • Committer: Bazaar Package Importer
  • Author(s): Matt Hope
  • Date: 2004-04-18 14:24:11 UTC
  • mfrom: (1.2.1 upstream) (2.1.1 warty)
  • Revision ID: james.westby@ubuntu.com-20040418142411-gr92uexquw4w8liq
Tags: 1.4-1
* New upstream release
* Examples and working code are installed by default to usr/bin,
  this has been moved to usr/share/doc/bioperl/bin

Show diffs side-by-side

added added

removed removed

Lines of Context:
scripts/blast_fetch_local.pl
1
 
#!/usr/local/bin/perl -w
2
 
 
3
 
# This is a working script to extract sequence ids from a blast 
4
 
# result file and extract those sequences from a local db.
5
 
 
6
 
# This may not work if the blastdb to index is est depending on 
7
 
# wha your machine supports
8
 
 
9
 
# $Id: blast_fetch_local.pl,v 1.1 2001/08/01 17:47:05 jason Exp $
10
 
 
11
 
use strict;
12
 
use Bio::Tools::BPlite;
13
 
use Bio::SeqIO;
14
 
use Bio::Index::Fasta;
15
 
 
16
 
use Getopt::Long;
17
 
use vars qw($USAGE);
18
 
$USAGE = 'usage: blast_fetch_local.pl -b report.bls -o outseqs.fa -d /path/to/db -i /path/to/index -p pvalue';
19
 
 
20
 
my ($db, $index, $output, $blastfile,$pvalue);
21
 
 
22
 
$pvalue = 1e-5;
23
 
 
24
 
&GetOptions
25
 
    (
26
 
     'd|db:s'  => \$db,
27
 
     'i|index:s' => \$index,
28
 
     'o|output:s' => \$output,
29
 
     'b|blast:s' => \$blastfile,
30
 
     'p|pvalue:s' => \$pvalue);
31
 
 
32
 
die $USAGE if( ! $db || ! -r $db || ! $index || ! $output 
33
 
               || ! $blastfile || ! -r $blastfile );
34
 
if( ! -r $index ) {
35
 
    # skip if it is already created
36
 
    my $indexdb = new Bio::Index::Fasta('-filename' => $index,
37
 
                                        '-write_flag' => 1);
38
 
    $indexdb->id_parser( \&parse_ncbi_id );
39
 
    $indexdb->make_index($db);
40
 
    undef $indexdb;
41
 
}
42
 
 
43
 
my $seqout = new Bio::SeqIO(-file => ">$output");
44
 
my $indexdb = new Bio::Index::Fasta('-filename' => $index);
45
 
my $bplite = new Bio::Tools::BPlite( -file => $blastfile );
46
 
 
47
 
SBJCT: while( my $sbjct = $bplite->nextSbjct ) {
48
 
    HSP: while( my $hsp = $sbjct->nextHSP ) { 
49
 
        if( $hsp->P > $pvalue) {        
50
 
            # skip Sbjcts that don't meet the minimum P value
51
 
            print 'skipping Sbjct with a HSP with pvalue=', $hsp->P, "\n";
52
 
            next SBJCT;    
53
 
        } 
54
 
# could also accept any Subject with an HSP that has appropriate e-value
55
 
#       else {
56
 
#           last HSP;
57
 
#       }
58
 
    }
59
 
    my  ($id) = split(/\s+/, $sbjct->name);
60
 
    # get the last value
61
 
    my @ids = split(/\|/, $id);
62
 
    $id = pop @ids;
63
 
    my $seq = $indexdb->fetch($id);
64
 
    if( ! $seq ) {
65
 
        print "could not find id '$id'\n";
66
 
    } else {
67
 
        $seqout->write_seq($seq);
68
 
    }
69
 
}
70
 
 
71
 
sub parse_ncbi_id {
72
 
    my @retvals;
73
 
    if( $_[0] =~ /^>(\S+)/ ) {
74
 
        my $val = $1;
75
 
        my (@elements) = split(/\|/, $val);
76
 
        while( @elements ) {
77
 
            my $id = shift @elements;
78
 
            if( $id eq 'gb' || $id eq 'gi' ) {
79
 
                $id = $id . '|' . shift @elements;
80
 
            }
81
 
            push @retvals, $id;
82
 
        }
83
 
    }
84
 
    return @retvals;
85
 
}