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 Support
  37
  38 Please direct usage questions or support issues to the mailing list:
  39
  40 L<bioperl-l@bioperl.org>
  41
  42 rather than to the module maintainer directly. Many experienced and
  43 reponsive experts will be able look at the problem and quickly
  44 address it. Please include a thorough description of the problem
  45 with code and data examples if at all possible.
  46
  47 =head2 Reporting Bugs
  48
  49 Report bugs to the Bioperl bug tracking system to help us keep track
  50 the bugs and their resolution.  Bug reports can be submitted via the
  51 web:
  52
  53   http://bugzilla.open-bio.org/
  54
  55 =head1 AUTHORS
  56
  57 Chris Fields
  58
  59 =head1 APPENDIX
  60
  61 The rest of the documentation details each of the object
  62 methods. Internal methods are usually preceded with a _
  63
  64 =cut
  65
  66 # Let the code begin...
  67
  68 package Bio::AlignIO::xmfa;
  69 use strict;
  70
  71 use base qw(Bio::AlignIO);
  72 our $WIDTH = 60;
  73
  74 =head2 next_aln
  75
  76  Title   : next_aln
  77  Usage   : $aln = $stream->next_aln
  78  Function: returns the next alignment in the stream.
  79  Returns : Bio::Align::AlignI object - returns 0 on end of file
  80             or on error
  81  Args    : -width => optional argument to specify the width sequence
  82            will be written (60 chars by default)
  83
  84 See L<Bio::Align::AlignI>
  85
  86 =cut
  87
  88 sub next_aln {
  89     my $self = shift;
  90     my ($width) = $self->_rearrange([qw(WIDTH)],@_);
  91     $self->width($width || $WIDTH);
  92
  93     my ($name, $tempname, $seqchar);
  94     my $aln = Bio::SimpleAlign->new();
  95     my $seqs = 0;
  96     # alignments
  97     while (defined (my $entry = $self->_readline) ) {
  98         chomp $entry;
  99         if ( index($entry, '=') == 0 ) {
 100             if (defined $name && $seqchar) {
 101                 my $seq = $self->_process_seq($name, $seqchar);
 102                 $aln->add_seq($seq);
 103             }
 104             if ($aln && $entry =~ m{score\s*=\s*(\d+)}) {
 105                 $aln->score($1);
 106             }
 107             $seqchar = '';
 108             undef $name;
 109             last;
 110         } elsif ( $entry =~ m{^>.+$}xms) {
 111             if ( defined $name ) {
 112                 my $seq = $self->_process_seq($name, $seqchar);
 113                 $aln->add_seq($seq);
 114             }
 115             $seqchar = '';
 116             $name = $entry;
 117         } else {
 118             $seqchar .= $entry;
 119         }
 120     }
 121
 122     # this catches last sequence if '=' is not present (Mauve)
 123     if ( defined $name ) {
 124         my $seq = $self->_process_seq($name, $seqchar);
 125         $aln->add_seq($seq);
 126     }
 127     $aln->no_sequences ? return $aln : return;
 128 }
 129
 130 =head2 write_aln
 131
 132  Title   : write_aln
 133  Usage   : $stream->write_aln(@aln)
 134  Function: writes the $aln object into the stream in xmfa format
 135  Returns : 1 for success and 0 for error
 136  Args    : L<Bio::Align::AlignI> object
 137
 138 See L<Bio::Align::AlignI>
 139
 140 =cut
 141
 142 sub write_aln {
 143     my ($self,@aln) = @_;
 144     my $width = $self->width;
 145     my ($seq,$desc,$rseq,$name,$count,$length,$seqsub,$start,$end,$strand,$id);
 146
 147     foreach my $aln (@aln) {
 148         if( ! $aln || ! $aln->isa('Bio::Align::AlignI')  ) {
 149             $self->warn("Must provide a Bio::Align::AlignI object when calling write_aln");
 150             next;
 151         }
 152         #if( $self->force_displayname_flat ) {
 153         #    $aln->set_displayname_flat(1);
 154         #}
 155         my $seqct = 1;
 156         foreach $rseq ( $aln->each_seq() ) {
 157             ($start, $end, $strand, $id) = ($rseq->start, $rseq->end, $rseq->strand || 0,
 158                                             $rseq->display_id);
 159             $strand = ($strand == 1)  ? '+' :
 160                       ($strand == -1) ? '-' :
 161                       '';
 162             $name = sprintf("%d:%d-%d %s %s",$seqct,$start,$end,$strand,$id);
 163             $seq  = $rseq->seq();
 164             $desc = $rseq->description || '';
 165             $self->_print (">$name $desc\n") or return ;
 166             $count = 0;
 167             $length = length($seq);
 168             if(defined $seq && $length > 0) {
 169             $seq =~ s/(.{1,$width})/$1\n/g;
 170             } else {
 171             $seq = "\n";
 172             }
 173             $self->_print($seq) || return 0;
 174             $seqct++;
 175         }
 176         my $alndesc = '';
 177         $alndesc = "score = ".$aln->score if ($aln->score);
 178         $self->_print("= $alndesc\n") || return 0;
 179
 180     }
 181     $self->flush if $self->_flush_on_write && defined $self->_fh;
 182     return 1;
 183 }
 184
 185 =head2 _get_len
 186
 187  Title   : _get_len
 188  Usage   :
 189  Function: determine number of alphabetic chars
 190  Returns : integer
 191  Args    : sequence string
 192
 193 =cut
 194
 195 sub _get_len {
 196     my ($self,$seq) = @_;
 197     $seq =~ s/[^A-Z]//gi;
 198     return CORE::length($seq);
 199 }
 200
 201 =head2 width
 202
 203  Title   : width
 204  Usage   : $obj->width($newwidth)
 205            $width = $obj->width;
 206  Function: Get/set width of alignment
 207  Returns : integer value of width
 208  Args    : on set, new value (a scalar or undef, optional)
 209
 210
 211 =cut
 212
 213 sub width{
 214     my $self = shift;
 215
 216     return $self->{'_width'} = shift if @_;
 217     return $self->{'_width'} || $WIDTH;
 218 }
 219
 220 ####### PRIVATE #######
 221
 222 sub _process_seq {
 223     my ($self, $entry, $seq) = @_;
 224     my ($start, $end, $strand, $seqname, $desc, $all);
 225     # put away last name and sequence
 226     if ( $entry =~ m{^>\s*\d+:(\d+)-(\d+)\s([+-]{1})(?:\s+(\S+)\s*(\S\.*)?)?} ) {
 227         ($start, $end, $seqname, $desc) = ($1, $2, $4, $5);
 228         $strand = ($4 eq '+')  ?  1  : -1;
 229     } else {
 230         $self->throw("Line does not comform to XMFA format:\n$entry");
 231     }
 232     my $seqobj = Bio::LocatableSeq->new(
 233              -nowarnonempty => 1,
 234              -strand      => $strand,
 235              -seq         => $seq,
 236              -display_id  => $seqname,
 237              -description => $desc || $all,
 238              -start       => $start,
 239              -end         => $end,
 240              );
 241     $self->debug("Reading $seqname\n");
 242     return $seqobj;
 243 }
 244
 245 1;