jpred/lib/PSIBLAST.pm

   1 package PSIBLAST;
   2
   3 use strict;
   4 use warnings;
   5 use Carp;
   6 use IO::String;
   7
   8 use base qw(Root Read);
   9
  10 =head1 NAME
  11
  12 PSIBLAST - Yet Another PSIBLAST Parser
  13
  14 =head1 SYNOPSIS
  15
  16   use PSIBLAST;
  17
  18   $psi = PSIBLAST->new;
  19   $psi->read_file("psiblast.output");
  20   $psi->get_all;
  21
  22 =head1 DESCRIPTION
  23
  24 A module to parse PSIBLAST output. It only grabs the last itteration, and there are simple methods to retrieve the ids of the sequences, the aligments and the start positions of the alignments.
  25
  26 Now gets the scores for the last iteration as well.
  27
  28 =head1 METHODS
  29
  30 =head2 $psi->read(\*FH)
  31
  32 Reads a PSIBLAST file from a filehandle.
  33
  34 =cut
  35
  36 sub read {
  37         my ($self, $fh) = @_;
  38
  39         croak "No filehandle passed to method read" unless ref $fh eq 'GLOB';
  40
  41         # Make sure that $/ is a defined value
  42         local $/ = "\n";
  43
  44         $self->_blastpgp_parse($fh);
  45 }
  46
  47 sub _blastpgp_parse {
  48         my ($self, $fh) = @_;
  49
  50         local $_;
  51         my $version = <$fh>;
  52         my $pos;
  53
  54         #die "Wrong version of BLASTP" if $version !~ /BLASTP 2.2.6 \[Apr-09-2003\]/;
  55
  56         # Find the position of the last itteration
  57         while (<$fh>) {
  58                 <$fh>, $pos = tell $fh if /^Sequences used in model and found again:/
  59         }
  60
  61         # and make sure we found it...
  62         unless ($pos) { $self->{psi} = undef, return }
  63         seek $fh, $pos, 0;
  64
  65         # Now read in the last itteration
  66         # @alignments holds the lines containing the aligment
  67         # @order is a list of the IDs
  68         # $flag indicates that the ID's have been collected
  69         #
  70         # All of this has to be list based as there may be multiple local
  71         # alignments against the same sequence
  72         my (@alignments, @order, $flag);
  73
  74         while (<$fh>) {
  75                 # Have we reached the end of the alignment
  76                 last if /^\s*Database: /;
  77
  78                 # Skip babble
  79                 next if /^Sequences not found previously or not previously below threshold:/;
  80                 #next if /^\s*$/;
  81
  82                 chomp;
  83                 # Capture E scores and check to see if we've reached the start of the alignment
  84                 unless (/^QUERY\s*/ or @alignments) {
  85                         next if /^$/;
  86                         my ($id, $bit, $e) = unpack "A70A3A*", $_;
  87                         ($id) = (split /\s+/, $id)[0];
  88                         s/\s*//g for $bit, $e;
  89
  90                         # Store the scores
  91                         $self->{$id}{bit} = $bit;
  92                         $self->{$id}{evalue} = $e;
  93                         next;
  94                 }
  95
  96                 # Have we reached the end of the first block of the
  97                 # alignment
  98                 $flag = 1, next if /^$/;
  99
 100                 # Add IDs if we're still in the first block
 101                 push @order, (split)[0] unless $flag;
 102                 # Add alignment lines
 103                 push @alignments, $_;
 104         }
 105
 106         $self->size(scalar @order);
 107
 108         # Add sequences and start positions to self
 109         for (0..$#alignments) {
 110                 my $i = $_ % @order;
 111                 my @fields = split ' ', $alignments[$_], 4;
 112
 113                 # blastpgp output allways has a ID at the begining of the line, but can
 114                 # have a variable number of fields after that. e.g.:
 115 # $seqid $start $seq                                                         $end
 116 # O24727        ------------------------------------------------------------
 117 # O24727   0    ------------------------------------------------------------
 118 # O24727   26   -Y--A---A-----L-----R----S-------L---------G-------P----V--- 35
 119
 120                 if (@fields == 4 or @fields == 3) {
 121                         $self->{psi}{$i}{seq} .= $fields[2];
 122                         $self->{psi}{$i}{start} = $fields[1]
 123                                 unless defined $self->{psi}{$i}{start};
 124                 }
 125                 elsif (@fields == 2) { $self->{psi}{$i}{seq} .= $fields[1] }
 126                 else {
 127                         croak "Fatal error! Wrong number of fields in BLAST alignment:\n$alignments[$_]\n"
 128                 }
 129         }
 130
 131         # Add IDs
 132         for (0..$#order) {
 133                 $self->{psi}{$_}{id} = $order[$_];
 134         }
 135
 136         1;
 137 }
 138
 139 =head2 @ids = $psi->get_ids
 140
 141 The IDs that were in the last iteration in the order they appeared.
 142
 143 =cut
 144
 145 sub get_ids {
 146         my ($self) = @_;
 147         #return map { $self->{psi}{$_}{id} } sort { $a <=> $b } keys %{$self->{psi}}
 148         return map { $self->{psi}{$_}{id} } 0..$self->size - 1;
 149 }
 150
 151 =head2 @alignments = $psi->get_alignments
 152
 153 The alignment from the last iteration in the same order as they appeared.
 154
 155 =cut
 156
 157 sub get_alignments {
 158         my ($self) = @_;
 159         #return map { $self->{psi}{$_}{seq} } sort { $a <=> $b } keys %{$self->{psi}};
 160         return map { $self->{psi}{$_}{seq} } 0..$self->size - 1;
 161 }
 162
 163 =head2 @starts = $psi->get_starts
 164
 165 The start position of the first residue from the alignments, returned in the same order as they appeared.
 166
 167 =cut
 168
 169 sub get_starts {
 170         my ($self) = @_;
 171         return map { $self->{psi}{$_}{start} } sort { $a <=> $b } keys %{$self->{psi}};
 172 }
 173
 174 =head2 @all = $psi->get_all
 175
 176 Returns an array of arrays containing the ID, start position and alignment from the PSIBLAST output in the order they appeared.
 177
 178 =cut
 179
 180 sub get_all {
 181         my ($self) = @_;
 182
 183         my @ids = $self->get_ids;
 184         my @starts = $self->get_starts;
 185         my @alignments = $self->get_alignments;
 186
 187         return map { [ $ids[$_], $starts[$_], $alignments[$_] ] } 0..$#ids
 188 }
 189
 190 =head2 $psi->evalue($label)
 191
 192 =cut
 193
 194 sub evalue {
 195         my ($self, $label) = @_;
 196         defined $label or confess "No label passed\n";
 197         exists $self->{evalue}{$label} or warn "No such label '$label'\n" and return undef;
 198         return $self->{evalue}{$label};
 199 }
 200
 201 =head2 $psi->bit($label)
 202
 203 =cut
 204
 205 sub bit {
 206         my ($self, $label) = @_;
 207         defined $label or confess "No label passed\n";
 208         exists $self->{bit}{$label} or warn "No such label '$label'\n" and return undef;
 209         return $self->{bit}{$label};
 210 }
 211
 212 =head2 $psi->size
 213
 214 Returns the number of sequences in the alignments.
 215
 216 =cut
 217
 218 sub size {
 219         my ($self, $size) = @_;
 220
 221         unless (defined $size) {
 222                 exists $self->{size} and return $self->{size};
 223         }
 224         else {
 225                 $self->{size} = $size;
 226         }
 227 }
 228
 229 =head2 BUGS
 230
 231 Doesn't aquire any other information.
 232
 233 =head2 AUTHOR
 234
 235 Jonathan Barber <jon@compbio.dundee.ac.uk>
 236
 237 =cut
 238
 239 1;