binaries/src/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   # Find the position of the last itteration
  55   while (<$fh>) {
  56     <$fh>, $pos = tell $fh if /^Sequences used in model and found again:/;
  57   }
  58
  59   # and make sure we found it...
  60   unless ($pos) { $self->{psi} = undef, return }
  61   seek $fh, $pos, 0;
  62
  63   # Now read in the last itteration
  64   # @alignments holds the lines containing the aligment
  65   # @order is a list of the IDs
  66   # $flag indicates that the ID's have been collected
  67   #
  68   # All of this has to be list based as there may be multiple local
  69   # alignments against the same sequence
  70   my ( @alignments, @order, $flag );
  71
  72   while (<$fh>) {
  73
  74     # Have we reached the end of the alignment
  75     last if /^\s*Database: /;
  76
  77     # Skip babble
  78     next if /^Sequences not found previously or not previously below threshold:/;
  79
  80     chomp;
  81
  82     # Capture E scores and check to see if we've reached the start of the alignment
  83     unless ( /^QUERY\s*/ or @alignments ) {
  84       next if /^$/;
  85       my ( $id, $bit, $e ) = unpack "A70A3A*", $_;
  86       ($id) = ( split /\s+/, $id )[0];
  87       s/\s*//g for $bit, $e;
  88
  89       # Store the scores
  90       $self->{$id}{bit}    = $bit;
  91       $self->{$id}{evalue} = $e;
  92       next;
  93     }
  94
  95     # Have we reached the end of the first block of the alignment
  96     $flag = 1, next if /^$/;
  97
  98     # Add IDs if we're still in the first block
  99     push @order, (split)[0] unless $flag;
 100
 101     # Add alignment lines
 102     push @alignments, $_;
 103   }
 104
 105   $self->size( scalar @order );
 106
 107   # Add sequences and start positions to self
 108   for ( 0 .. $#alignments ) {
 109     my $i = $_ % @order;
 110     my @fields = split ' ', $alignments[$_], 4;
 111
 112     # blastpgp output allways has a ID at the begining of the line, but can
 113     # have a variable number of fields after that. e.g.:
 114     # $seqid $start $seq                                                         $end
 115     # O24727        ------------------------------------------------------------
 116     # O24727   0    ------------------------------------------------------------
 117     # O24727   26   -Y--A---A-----L-----R----S-------L---------G-------P----V--- 35
 118
 119     if ( @fields == 4 or @fields == 3 ) {
 120       $self->{psi}{$i}{seq} .= $fields[2];
 121       $self->{psi}{$i}{start} = $fields[1]
 122         unless defined $self->{psi}{$i}{start};
 123     } elsif ( @fields == 2 ) {
 124       $self->{psi}{$i}{seq} .= $fields[1];
 125     } else {
 126       croak "Fatal error! Wrong number of fields in BLAST alignment:\n$alignments[$_]\n";
 127     }
 128   }
 129
 130   # Add IDs
 131   for ( 0 .. $#order ) {
 132     $self->{psi}{$_}{id} = $order[$_];
 133   }
 134
 135   1;
 136 }
 137
 138 =head2 @ids = $psi->get_ids
 139
 140 The IDs that were in the last iteration in the order they appeared.
 141
 142 =cut
 143
 144 sub get_ids {
 145   my ($self) = @_;
 146
 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
 160   #return map { $self->{psi}{$_}{seq} } sort { $a <=> $b } keys %{$self->{psi}};
 161   return map { $self->{psi}{$_}{seq} } 0 .. $self->size - 1;
 162 }
 163
 164 =head2 @starts = $psi->get_starts
 165
 166 The start position of the first residue from the alignments, returned in the same order as they appeared.
 167
 168 =cut
 169
 170 sub get_starts {
 171   my ($self) = @_;
 172   return map { $self->{psi}{$_}{start} } sort { $a <=> $b } keys %{ $self->{psi} };
 173 }
 174
 175 =head2 @all = $psi->get_all
 176
 177 Returns an array of arrays containing the ID, start position and alignment from the PSIBLAST output in the order they appeared.
 178
 179 =cut
 180
 181 sub get_all {
 182   my ($self) = @_;
 183
 184   my @ids        = $self->get_ids;
 185   my @starts     = $self->get_starts;
 186   my @alignments = $self->get_alignments;
 187
 188   return map { [ $ids[$_], $starts[$_], $alignments[$_] ] } 0 .. $#ids;
 189 }
 190
 191 =head2 $psi->evalue($label)
 192
 193 =cut
 194
 195 sub evalue {
 196   my ( $self, $label ) = @_;
 197   defined $label or confess "No label passed\n";
 198   exists $self->{evalue}{$label} or warn "No such label '$label'\n" and return undef;
 199   return $self->{evalue}{$label};
 200 }
 201
 202 =head2 $psi->bit($label)
 203
 204 =cut
 205
 206 sub bit {
 207   my ( $self, $label ) = @_;
 208   defined $label or confess "No label passed\n";
 209   exists $self->{bit}{$label} or warn "No such label '$label'\n" and return undef;
 210   return $self->{bit}{$label};
 211 }
 212
 213 =head2 $psi->size
 214
 215 Returns the number of sequences in the alignments.
 216
 217 =cut
 218
 219 sub size {
 220   my ( $self, $size ) = @_;
 221
 222   unless ( defined $size ) {
 223     exists $self->{size} and return $self->{size};
 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;