lib/Bio/AlignIO/msf.pm

   1 #
   2 # BioPerl module for Bio::AlignIO::msf
   3 #   based on the Bio::SeqIO::msf module
   4 #       by Ewan Birney <birney@ebi.ac.uk>
   5 #       and Lincoln Stein  <lstein@cshl.org>
   6 #
   7 #       and the SimpleAlign.pm module of Ewan Birney
   8 #
   9 # Copyright Peter Schattner
  10 #
  11 # You may distribute this module under the same terms as perl itself
  12 # _history
  13 # September 5, 2000
  14 # POD documentation - main docs before the code
  15
  16 =head1 NAME
  17
  18 Bio::AlignIO::msf - msf sequence input/output stream
  19
  20 =head1 SYNOPSIS
  21
  22 Do not use this module directly.  Use it via the L<Bio::AlignIO> class.
  23
  24 =head1 DESCRIPTION
  25
  26 This object can transform L<Bio::Align::AlignI> objects to and from msf
  27 flat file databases.
  28
  29 =head1 FEEDBACK
  30
  31 =head2 Support
  32
  33 Please direct usage questions or support issues to the mailing list:
  34
  35 I<bioperl-l@bioperl.org>
  36
  37 rather than to the module maintainer directly. Many experienced and
  38 reponsive experts will be able look at the problem and quickly
  39 address it. Please include a thorough description of the problem
  40 with code and data examples if at all possible.
  41
  42 =head2 Reporting Bugs
  43
  44 Report bugs to the Bioperl bug tracking system to help us keep track
  45 the bugs and their resolution.  Bug reports can be submitted via the
  46 web:
  47
  48   https://github.com/bioperl/bioperl-live/issues
  49
  50 =head1 AUTHORS - Peter Schattner
  51
  52 Email: schattner@alum.mit.edu
  53
  54
  55 =head1 APPENDIX
  56
  57 The rest of the documentation details each of the object
  58 methods. Internal methods are usually preceded with a _
  59
  60 =cut
  61
  62 # Let the code begin...
  63
  64 package Bio::AlignIO::msf;
  65
  66 use vars qw(%valid_type);
  67 use strict;
  68
  69 use Bio::SeqIO::gcg; # for GCG_checksum()
  70 use Bio::SimpleAlign;
  71
  72 use base qw(Bio::AlignIO);
  73
  74 BEGIN {
  75         %valid_type = qw( dna N rna N protein P );
  76 }
  77
  78 =head2 next_aln
  79
  80  Title   : next_aln
  81  Usage   : $aln = $stream->next_aln()
  82  Function: returns the next alignment in the stream. Tries to read *all* MSF
  83            It reads all non whitespace characters in the alignment
  84            area. For MSFs with weird gaps (eg ~~~) map them by using
  85            $aln->map_chars('~','-')
  86  Returns : Bio::Align::AlignI object
  87  Args    : NONE
  88
  89 =cut
  90
  91 sub next_aln {
  92     my $self = shift;
  93     my $entry;
  94     my (%hash,$name,$str,@names,$seqname,$start,$end,$count,$seq);
  95
  96     my $aln =  Bio::SimpleAlign->new(-source => 'gcg' );
  97
  98     while( $entry = $self->_readline) {
  99         $entry =~ m{//} && last; # move to alignment section
 100         $entry =~ /Name:\s+(\S+)/ && do { $name = $1;
 101                                           $hash{$name} = ""; # blank line
 102                                           push(@names,$name); # we need it ordered!
 103                                       };
 104         # otherwise - skip
 105     }
 106
 107     # alignment section
 108
 109     while( $entry = $self->_readline) {
 110         next if ( $entry =~ /^\s+(\d+)/ ) ;
 111         $entry =~ /^\s*(\S+)\s+(.*)$/ && do {
 112             $name = $1;
 113             $str = $2;
 114             if( ! exists $hash{$name} ) {
 115                 $self->throw("$name exists as an alignment line but not in the header. Not confident of what is going on!");
 116             }
 117             $str =~ s/\s//g;
 118             $str =~ s/~/-/g;
 119             $hash{$name} .= $str;
 120         };
 121     }
 122
 123     return if @names < 1;
 124
 125     # now got this as a name - sequence hash. Let's make some sequences!
 126
 127     for $name ( @names ) {
 128         if( $name =~ m{(\S+)/(\d+)-(\d+)} ) {
 129             $seqname = $1;
 130             $start = $2;
 131             $end = $3;
 132         } else {
 133             $seqname = $name;
 134             $start = 1;
 135             $str = $hash{$name};
 136             $str =~ s/[^0-9A-Za-z$Bio::LocatableSeq::OTHER_SYMBOLS]//g;
 137
 138             $end = length($str);
 139         }
 140
 141         $seq = Bio::LocatableSeq->new('-seq'        => $hash{$name},
 142                                       '-display_id' => $seqname,
 143                                       '-start'      => $start,
 144                                       '-end'        => $end,
 145                                       '-alphabet'   => $self->alphabet,
 146                                       );
 147         $aln->add_seq($seq);
 148
 149 #  If $end <= 0, we have either reached the end of
 150 #  file in <> or we have encountered some other error
 151
 152     }
 153     return $aln if $aln->num_sequences;
 154     return;
 155 }
 156
 157
 158 =head2 write_aln
 159
 160  Title   : write_aln
 161  Usage   : $stream->write_aln(@aln)
 162  Function: writes the $aln object into the stream in MSF format
 163            Sequence type of the alignment is determined by the first sequence.
 164  Returns : 1 for success and 0 for error
 165  Args    : Bio::Align::AlignI object
 166
 167
 168 =cut
 169
 170 sub write_aln {
 171         my ($self,@aln) = @_;
 172         my $msftag;
 173         my $type;
 174         my $count = 0;
 175         my $maxname;
 176         my ($length,$date,$name,$seq,$miss,$pad,%hash,@arr,$tempcount,$index);
 177         foreach my $aln (@aln) {
 178                 if( ! $aln || ! $aln->isa('Bio::Align::AlignI')  ) {
 179                         $self->warn("Must provide a Bio::Align::AlignI object when calling write_aln");
 180                         next;
 181                 }
 182                 $date = localtime(time);
 183                 $msftag = "MSF";
 184                 $type = $valid_type{$aln->get_seq_by_pos(1)->alphabet};
 185                 $maxname = $aln->maxdisplayname_length();
 186                 $length  = $aln->length();
 187                 $name = $aln->id();
 188                 if( !defined $name ) {
 189                         $name = "Align";
 190                 }
 191
 192                 $self->_print (sprintf("\n%s   MSF: %d  Type: %s  %s  Check: 00 ..\n\n",
 193                                $name,  $aln->num_sequences, $type, $date));
 194
 195     my $seqCountFormat = "%".($maxname > 20 ? $maxname + 2: 22)."s%-27d%27d\n";
 196     my $seqNameFormat  = "%-".($maxname > 20 ? $maxname : 20)."s  ";
 197
 198                 foreach $seq ( $aln->each_seq() ) {
 199                         $name = $aln->displayname($seq->get_nse());
 200                         $miss = $maxname - length ($name);
 201                         $miss += 2;
 202                         $pad  = " " x $miss;
 203
 204                         $self->_print (sprintf(" Name: %s%sLen:    %d  Check:  %d  Weight:  1.00\n",$name,$pad,length $seq->seq(), Bio::SeqIO::gcg->GCG_checksum($seq)));
 205
 206                         $hash{$name} = $seq->seq();
 207                         push(@arr,$name);
 208                 }
 209         # ok - heavy handed, but there you go.
 210         #
 211         $self->_print ("\n//\n\n\n");
 212
 213         while( $count < $length ) {
 214                         # there is another block to go!
 215                         $self->_print (sprintf($seqCountFormat,' ',$count+1,$count+50));
 216                         foreach $name  ( @arr ) {
 217                                 $self->_print (sprintf($seqNameFormat,$name));
 218
 219                                 $tempcount = $count;
 220                                 $index = 0;
 221                                 while( ($tempcount + 10 < $length) && ($index < 5)  ) {
 222
 223                                         $self->_print (sprintf("%s ",substr($hash{$name},
 224                                                                                                                                         $tempcount,10)));
 225
 226                                         $tempcount += 10;
 227                                         $index++;
 228                                 }               #
 229                                 # ok, could be the very last guy ;)
 230                                 #
 231                                 if( $index < 5) {
 232                                         # space to print!
 233                                         #
 234                                         $self->_print (sprintf("%s ",substr($hash{$name},$tempcount)));
 235                                         $tempcount += 10;
 236                                 }
 237                                 $self->_print ("\n");
 238                         }
 239                         $self->_print ("\n\n");
 240                         $count = $tempcount;
 241         }
 242         }
 243         $self->flush if $self->_flush_on_write && defined $self->_fh;
 244         return 1;
 245 }
 246
 247 1;