Bio/AlignIO/psi.pm

   1 # $Id$
   2 #
   3 # BioPerl module for Bio::AlignIO::psi
   4 #
   5 # Cared for by Jason Stajich <jason@bioperl.org>
   6 #
   7 # Copyright Jason Stajich
   8 #
   9 # You may distribute this module under the same terms as perl itself
  10
  11 # POD documentation - main docs before the code
  12
  13 =head1 NAME
  14
  15 Bio::AlignIO::psi - Read/Write PSI-BLAST profile alignment files
  16
  17 =head1 SYNOPSIS
  18
  19 This module will parse PSI-BLAST output of the format seqid XXXX
  20
  21 =head1 DESCRIPTION
  22
  23 This is a parser for psi-blast blocks.
  24
  25 =head1 FEEDBACK
  26
  27 =head2 Mailing Lists
  28
  29 User feedback is an integral part of the evolution of this and other
  30 Bioperl modules. Send your comments and suggestions preferably to
  31 the Bioperl mailing list.  Your participation is much appreciated.
  32
  33   bioperl-l@bioperl.org                  - General discussion
  34   http://bioperl.org/wiki/Mailing_lists  - About the mailing lists
  35
  36 =head2 Reporting Bugs
  37
  38 Report bugs to the Bioperl bug tracking system to help us keep track
  39 of the bugs and their resolution. Bug reports can be submitted via
  40 the web:
  41
  42   http://bugzilla.open-bio.org/
  43
  44 =head1 AUTHOR - Jason Stajich
  45
  46 Email jason@bioperl.org
  47
  48 =head1 APPENDIX
  49
  50 The rest of the documentation details each of the object methods.
  51 Internal methods are usually preceded with a _
  52
  53 =cut
  54
  55
  56 # Let the code begin...
  57
  58
  59 package Bio::AlignIO::psi;
  60 use vars qw($BlockLen $IdLength);
  61 use strict;
  62
  63 $BlockLen = 100;
  64 $IdLength = 13;
  65
  66 # Object preamble - inherits from Bio::Root::Root
  67
  68 use Bio::SimpleAlign;
  69 use Bio::LocatableSeq;
  70
  71 use base qw(Bio::AlignIO);
  72
  73 =head2 new
  74
  75  Title   : new
  76  Usage   : my $obj = Bio::AlignIO::psi->new();
  77  Function: Builds a new Bio::AlignIO::psi object
  78  Returns : Bio::AlignIO::psi
  79  Args    :
  80
  81 =cut
  82
  83 =head2 next_aln
  84
  85  Title   : next_aln
  86  Usage   : $aln = $stream->next_aln()
  87  Function: returns the next alignment in the stream
  88  Returns : Bio::Align::AlignI object
  89  Args    : NONE
  90
  91 See L<Bio::Align::AlignI>
  92
  93 =cut
  94
  95 sub next_aln {
  96     my ($self) = @_;
  97     my $aln;
  98     my %seqs;
  99     my @order;
 100     while( defined ($_ = $self->_readline ) ) {
 101         next if( /^\s+$/);
 102         if( !defined $aln ) {
 103             $aln = Bio::SimpleAlign->new();
 104         }
 105         my ($id,$s) = split;
 106         push @order, $id if( ! defined $seqs{$id});
 107         $seqs{$id} .= $s;
 108     }
 109     foreach my $id ( @order) {
 110     my $gaps = $seqs{$id} =~ tr/-/-/;
 111         my $seq = Bio::LocatableSeq->new(-seq => $seqs{$id},
 112                                         -id  => $id,
 113                                         -start => 1,
 114                                         -end   => length($seqs{$id}) - $gaps
 115                     );
 116         $aln->add_seq($seq);
 117     }
 118     return $aln;
 119 }
 120
 121 =head2 write_aln
 122
 123  Title   : write_aln
 124  Usage   : $stream->write_aln(@aln)
 125  Function: writes the NCBI psi-format object (.aln) into the stream
 126  Returns : 1 for success and 0 for error
 127  Args    : Bio::Align::AlignI object
 128
 129 L<Bio::Align::AlignI>
 130
 131 =cut
 132
 133 sub write_aln {
 134         my ($self,$aln) = @_;
 135         unless( defined $aln && ref($aln) &&
 136                           $aln->isa('Bio::Align::AlignI') ) {
 137                 $self->warn("Must provide a valid Bio::Align::AlignI to write_aln");
 138                 return 0;
 139         }
 140         my $ct = 0;
 141         my @seqs = $aln->each_seq;
 142         my $len = 1;
 143         my $alnlen = $aln->length;
 144         my $idlen = $IdLength;
 145         my @ids = map { substr($_->display_id,0,$idlen) } @seqs;
 146         while( $len < ($alnlen + 1) ) {
 147                 my $start = $len;
 148                 my $end   = $len + $BlockLen;
 149                 $end = $alnlen if ( $end > $alnlen );
 150                 my $c = 0;
 151                 foreach my $seq ( @seqs ) {
 152                         $self->_print(sprintf("%-".$idlen."s %s\n",
 153                                                                                  $ids[$c++],
 154                                                                                  $seq->subseq($start,$end)));
 155                 }
 156                 $self->_print("\n");
 157                 $len += $BlockLen+1;
 158         }
 159         $self->flush if $self->_flush_on_write && defined $self->_fh;
 160         return 1;
 161 }
 162
 163 1;