Bio/AlignIO/fasta.pm

   1 #
   2 # BioPerl module for Bio::AlignIO::fasta
   3 #
   4 # Copyright Peter Schattner
   5 #
   6 # You may distribute this module under the same terms as perl itself
   7 # POD documentation - main docs before the code
   8
   9 =head1 NAME
  10
  11 Bio::AlignIO::fasta - fasta MSA Sequence input/output stream
  12
  13 =head1 SYNOPSIS
  14
  15 Do not use this module directly.  Use it via the L<Bio::AlignIO>
  16 class.
  17
  18 =head1 DESCRIPTION
  19
  20 This object can transform L<Bio::SimpleAlign> objects to and from
  21 fasta flat files.  This is for the fasta alignment format, not
  22 for the FastA sequence analysis program.  To process the alignments from
  23 FastA (FastX, FastN, FastP, tFastA, etc) use the Bio::SearchIO module.
  24
  25 =head1 FEEDBACK
  26
  27 =head2 Support
  28
  29 Please direct usage questions or support issues to the mailing list:
  30
  31 I<bioperl-l@bioperl.org>
  32
  33 rather than to the module maintainer directly. Many experienced and
  34 reponsive experts will be able look at the problem and quickly
  35 address it. Please include a thorough description of the problem
  36 with code and data examples if at all possible.
  37
  38 =head2 Reporting Bugs
  39
  40 Report bugs to the Bioperl bug tracking system to help us keep track
  41 the bugs and their resolution.  Bug reports can be submitted via the
  42 web:
  43
  44   https://github.com/bioperl/bioperl-live/issues
  45
  46 =head1 AUTHORS
  47
  48 Peter Schattner
  49
  50 =head1 APPENDIX
  51
  52 The rest of the documentation details each of the object
  53 methods. Internal methods are usually preceded with a _
  54
  55 =cut
  56
  57 # Let the code begin...
  58
  59 package Bio::AlignIO::fasta;
  60 use strict;
  61
  62 use base qw(Bio::AlignIO);
  63 our $WIDTH = 60;
  64 use Bio::LocatableSeq;
  65
  66 =head2 next_aln
  67
  68  Title   : next_aln
  69  Usage   : $aln = $stream->next_aln
  70  Function: returns the next alignment in the stream.
  71  Returns : Bio::Align::AlignI object - returns 0 on end of file
  72             or on error
  73  Args    : -width => optional argument to specify the width sequence
  74            will be written (60 chars by default)
  75
  76 See L<Bio::Align::AlignI>
  77
  78 =cut
  79
  80 sub next_aln {
  81     my $self = shift;
  82     my ($width) = $self->_rearrange( [qw(WIDTH)], @_ );
  83     $self->width( $width || $WIDTH );
  84
  85     my ($start, $end,      $name,     $seqname, $seq,  $seqchar,
  86         $entry, $tempname, $tempdesc, %align,   $desc, $maxlen
  87     );
  88     my $aln = Bio::SimpleAlign->new();
  89
  90     while ( defined( $entry = $self->_readline ) ) {
  91         chomp $entry;
  92         if ( $entry =~ s/^>\s*(\S+)\s*// ) {
  93             $tempname = $1;
  94             chomp($entry);
  95             $tempdesc = $entry;
  96             if ( defined $name ) {
  97                 $seqchar =~ s/\s//g;
  98                 $seqname = $name;
  99                 $start   = 1;
 100                 $end     = $self->_get_len($seqchar);
 101                 $seq     = Bio::LocatableSeq->new(
 102                     -seq         => $seqchar,
 103                     -display_id  => $seqname,
 104                     -description => $desc,
 105                     -start       => $start,
 106                     -end         => $end,
 107                     -alphabet    => $self->alphabet,
 108                 );
 109                 $aln->add_seq($seq);
 110                 $self->debug("Reading $seqname\n");
 111             }
 112             $desc    = $tempdesc;
 113             $name    = $tempname;
 114             $desc    = $entry;
 115             $seqchar = "";
 116             next;
 117         }
 118
 119         # removed redundant symbol validation
 120         # this is already done in Bio::PrimarySeq
 121         $seqchar .= $entry;
 122     }
 123
 124     #  Next two lines are to silence warnings that
 125     #  otherwise occur at EOF when using <$fh>
 126     $name    = "" if ( !defined $name );
 127     $seqchar = "" if ( !defined $seqchar );
 128     $seqchar =~ s/\s//g;
 129
 130     #  Put away last name and sequence
 131     if ( $name =~ /(\S+)\/(\d+)-(\d+)$/ ) {
 132         $seqname = $1;
 133         $start   = $2;
 134         $end     = $3;
 135     }
 136     else {
 137         $seqname = $name;
 138         $start   = 1;
 139         $end     = $self->_get_len($seqchar);
 140     }
 141
 142     # This logic now also reads empty lines at the
 143     # end of the file. Skip this is seqchar and seqname is null
 144     unless ( length($seqchar) == 0 && length($seqname) == 0 ) {
 145         $seq = Bio::LocatableSeq->new(
 146             -seq         => $seqchar,
 147             -display_id  => $seqname,
 148             -description => $desc,
 149             -start       => $start,
 150             -end         => $end,
 151             -alphabet    => $self->alphabet,
 152         );
 153         $aln->add_seq($seq);
 154         $self->debug("Reading $seqname\n");
 155     }
 156     my $alnlen = $aln->length;
 157     foreach my $seq ( $aln->each_seq ) {
 158         if ( $seq->length < $alnlen ) {
 159             my ($diff) = ( $alnlen - $seq->length );
 160             $seq->seq( $seq->seq() . "-" x $diff );
 161         }
 162     }
 163
 164     # no sequences means empty alignment (possible EOF)
 165     return $aln if $aln->num_sequences;
 166     return;
 167 }
 168
 169
 170 =head2 write_aln
 171
 172  Title   : write_aln
 173  Usage   : $stream->write_aln(@aln)
 174  Function: writes the $aln object into the stream in fasta format
 175  Returns : 1 for success and 0 for error
 176  Args    : L<Bio::Align::AlignI> object
 177
 178 See L<Bio::Align::AlignI>
 179
 180 =cut
 181
 182 sub write_aln {
 183     my ($self,@aln) = @_;
 184     my $width = $self->width;
 185     my ($seq,$desc,$rseq,$name,$count,$length,$seqsub);
 186
 187     foreach my $aln (@aln) {
 188         if( ! $aln || ! $aln->isa('Bio::Align::AlignI')  ) {
 189             $self->warn("Must provide a Bio::Align::AlignI object when calling write_aln");
 190             next;
 191         }
 192         if( $self->force_displayname_flat ) {
 193             $aln->set_displayname_flat(1);
 194         }
 195         foreach $rseq ( $aln->each_seq() ) {
 196             $name = $aln->displayname($rseq->get_nse());
 197             $seq  = $rseq->seq();
 198             $desc = $rseq->description || '';
 199             $desc = ' '.$desc if $desc;
 200             $self->_print (">$name$desc\n") or return;
 201             $count = 0;
 202             $length = length($seq);
 203             if(defined $seq && $length > 0) {
 204                 $seq =~ s/(.{1,$width})/$1\n/g;
 205             } else {
 206                 $seq = "\n";
 207             }
 208             $self->_print($seq);
 209         }
 210     }
 211     $self->flush if $self->_flush_on_write && defined $self->_fh;
 212     return 1;
 213 }
 214
 215 =head2 _get_len
 216
 217  Title   : _get_len
 218  Usage   :
 219  Function: determine number of alphabetic chars
 220  Returns : integer
 221  Args    : sequence string
 222
 223 =cut
 224
 225 sub _get_len {
 226     my ($self,$seq) = @_;
 227     my $chars = $Bio::LocatableSeq::GAP_SYMBOLS.$Bio::LocatableSeq::FRAMESHIFT_SYMBOLS;
 228     $seq =~ s{[$chars]+}{}gi;
 229     return CORE::length($seq);
 230 }
 231
 232 =head2 width
 233
 234  Title   : width
 235  Usage   : $obj->width($newwidth)
 236            $width = $obj->width;
 237  Function: Get/set width of alignment
 238  Returns : integer value of width
 239  Args    : on set, new value (a scalar or undef, optional)
 240
 241
 242 =cut
 243
 244 sub width{
 245     my $self = shift;
 246
 247     return $self->{'_width'} = shift if @_;
 248     return $self->{'_width'} || $WIDTH;
 249 }
 250
 251 1;