Bio/AlignIO/xmfa.pm

   1 # $Id$
   2 #
   3 # BioPerl module for Bio::AlignIO::xmfa
   4 #
   5 # Copyright Chris Fields
   6 #
   7 # You may distribute this module under the same terms as perl itself
   8 # POD documentation - main docs before the code
   9
  10 =head1 NAME
  11
  12 Bio::AlignIO::xmfa - XMFA MSA Sequence input/output stream
  13
  14 =head1 SYNOPSIS
  15
  16 Do not use this module directly.  Use it via the L<Bio::AlignIO>
  17 class.
  18
  19 =head1 DESCRIPTION
  20
  21 This object can transform L<Bio::SimpleAlign> objects from
  22 XMFA flat file databases.  For more information, see:
  23
  24   http://gel.ahabs.wisc.edu/docserver/mauve/files.stx
  25
  26 This module is based on the AlignIO::fasta parser written by
  27 Peter Schattner
  28
  29 =head1 TODO
  30
  31 Finish write_aln(), clean up code, allow LargeLocatableSeq (ie for
  32 very large sequences a'la Mauve)
  33
  34 =head1 FEEDBACK
  35
  36 =head2 Reporting Bugs
  37
  38 Report bugs to the Bioperl bug tracking system to help us keep track
  39 the bugs and their resolution.  Bug reports can be submitted via the
  40 web:
  41
  42   http://bugzilla.open-bio.org/
  43
  44 =head1 AUTHORS
  45
  46 Chris Fields
  47
  48 =head1 APPENDIX
  49
  50 The rest of the documentation details each of the object
  51 methods. Internal methods are usually preceded with a _
  52
  53 =cut
  54
  55 # Let the code begin...
  56
  57 package Bio::AlignIO::xmfa;
  58 use strict;
  59
  60 use base qw(Bio::AlignIO);
  61 our $WIDTH = 60;
  62
  63 =head2 next_aln
  64
  65  Title   : next_aln
  66  Usage   : $aln = $stream->next_aln
  67  Function: returns the next alignment in the stream.
  68  Returns : Bio::Align::AlignI object - returns 0 on end of file
  69             or on error
  70  Args    : -width => optional argument to specify the width sequence
  71            will be written (60 chars by default)
  72
  73 See L<Bio::Align::AlignI>
  74
  75 =cut
  76
  77 sub next_aln {
  78     my $self = shift;
  79     my ($width) = $self->_rearrange([qw(WIDTH)],@_);
  80     $self->width($width || $WIDTH);
  81
  82     my ($name, $tempname, $seqchar);
  83     my $aln = Bio::SimpleAlign->new();
  84     my $seqs = 0;
  85     # alignments
  86     while (defined (my $entry = $self->_readline) ) {
  87         chomp $entry;
  88         if ( index($entry, '=') == 0 ) {
  89             if (defined $name && $seqchar) {
  90                 my $seq = $self->_process_seq($name, $seqchar);
  91                 $aln->add_seq($seq);
  92             }
  93             if ($aln && $entry =~ m{score\s*=\s*(\d+)}) {
  94                 $aln->score($1);
  95             }
  96             $seqchar = '';
  97             $name = $entry;
  98             last;
  99         } elsif ( $entry =~ m{^>.+$}xms) {
 100             if ( defined $name ) {
 101                 my $seq = $self->_process_seq($name, $seqchar);
 102                 $aln->add_seq($seq);
 103             }
 104             $seqchar = '';
 105             $name = $entry;
 106         } else {
 107             $seqchar .= $entry;
 108         }
 109     }
 110
 111     # this catches last sequence if '=' is not present (Mauve)
 112     if ( defined $name ) {
 113         my $seq = $self->_process_seq($name, $seqchar);
 114         $aln->add_seq($seq);
 115     }
 116     $aln->no_sequences ? return $aln : return;
 117 }
 118
 119 =head2 write_aln
 120
 121  Title   : write_aln
 122  Usage   : $stream->write_aln(@aln)
 123  Function: writes the $aln object into the stream in xmfa format
 124  Returns : 1 for success and 0 for error
 125  Args    : L<Bio::Align::AlignI> object
 126
 127 See L<Bio::Align::AlignI>
 128
 129 =cut
 130
 131 sub write_aln {
 132     my ($self,@aln) = @_;
 133     my $width = $self->width;
 134     my ($seq,$desc,$rseq,$name,$count,$length,$seqsub,$start,$end,$strand,$id);
 135
 136     foreach my $aln (@aln) {
 137         if( ! $aln || ! $aln->isa('Bio::Align::AlignI')  ) {
 138             $self->warn("Must provide a Bio::Align::AlignI object when calling write_aln");
 139             next;
 140         }
 141         #if( $self->force_displayname_flat ) {
 142         #    $aln->set_displayname_flat(1);
 143         #}
 144         my $seqct = 1;
 145         foreach $rseq ( $aln->each_seq() ) {
 146             ($start, $end, $strand, $id) = ($rseq->start, $rseq->end, $rseq->strand,
 147                                             $rseq->display_id);
 148             $strand = ($strand == 1)  ? '+' :
 149                       ($strand == -1) ? '-' :
 150                       '';
 151             $name = sprintf("%d:%d-%d %s %s",$seqct,$start,$end,$strand,$id);
 152             $seq  = $rseq->seq();
 153             $desc = $rseq->description || '';
 154             $self->_print (">$name $desc\n") or return ;
 155             $count = 0;
 156             $length = length($seq);
 157             if(defined $seq && $length > 0) {
 158             $seq =~ s/(.{1,$width})/$1\n/g;
 159             } else {
 160             $seq = "\n";
 161             }
 162             $self->_print($seq) || return 0;
 163             $seqct++;
 164         }
 165         my $alndesc = '';
 166         $alndesc = "score = ".$aln->score if ($aln->score);
 167         $self->_print("= $alndesc\n") || return 0;
 168
 169     }
 170     $self->flush if $self->_flush_on_write && defined $self->_fh;
 171     return 1;
 172 }
 173
 174 =head2 _get_len
 175
 176  Title   : _get_len
 177  Usage   :
 178  Function: determine number of alphabetic chars
 179  Returns : integer
 180  Args    : sequence string
 181
 182 =cut
 183
 184 sub _get_len {
 185     my ($self,$seq) = @_;
 186     $seq =~ s/[^A-Z]//gi;
 187     return CORE::length($seq);
 188 }
 189
 190 =head2 width
 191
 192  Title   : width
 193  Usage   : $obj->width($newwidth)
 194            $width = $obj->width;
 195  Function: Get/set width of alignment
 196  Returns : integer value of width
 197  Args    : on set, new value (a scalar or undef, optional)
 198
 199
 200 =cut
 201
 202 sub width{
 203     my $self = shift;
 204
 205     return $self->{'_width'} = shift if @_;
 206     return $self->{'_width'} || $WIDTH;
 207 }
 208
 209 ####### PRIVATE #######
 210
 211 sub _process_seq {
 212     my ($self, $entry, $seq) = @_;
 213     my ($start, $end, $strand, $seqname, $desc, $all);
 214     # put away last name and sequence
 215     if ( $entry =~ m{^>\s*\d+:(\d+)-(\d+)\s([+-]{1})(?:\s+(\S+)\s*(\S\.*)?)?} ) {
 216         ($start, $end, $seqname, $desc) = ($1, $2, $4, $5);
 217         $strand = ($4 eq '+')  ?  1  : -1;
 218     } else {
 219         $self->throw("Line does not comform to XMFA format:\n$entry");
 220     }
 221     my $seqobj = Bio::LocatableSeq->new(
 222              -nowarnonempty => 1,
 223              -strand      => $strand,
 224              -seq         => $seq,
 225              -display_id  => $seqname,
 226              -description => $desc || $all,
 227              -start       => $start,
 228              -end         => $end,
 229              );
 230     $self->debug("Reading $seqname\n");
 231     return $seqobj;
 232 }
 233
 234 1;