lib/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
  61 use strict;
  62
  63 use base qw(Bio::AlignIO);
  64 our $WIDTH = 60;
  65 use Bio::LocatableSeq;
  66
  67 =head2 next_aln
  68
  69  Title   : next_aln
  70  Usage   : $aln = $stream->next_aln
  71  Function: returns the next alignment in the stream.
  72  Returns : Bio::Align::AlignI object - returns 0 on end of file
  73             or on error
  74  Args    : -width => optional argument to specify the width sequence
  75            will be written (60 chars by default)
  76
  77 See L<Bio::Align::AlignI>
  78
  79 =cut
  80
  81 sub next_aln {
  82     my $self = shift;
  83     my ($width) = $self->_rearrange( [qw(WIDTH)], @_ );
  84     $self->width( $width || $WIDTH );
  85
  86     my ($start, $end,      $name,     $seqname, $seq,  $seqchar,
  87         $entry, $tempname, $tempdesc, %align,   $desc, $maxlen
  88     );
  89     my $aln = Bio::SimpleAlign->new();
  90
  91     while ( defined( $entry = $self->_readline ) ) {
  92         chomp $entry;
  93         if ( $entry =~ s/^>\s*(\S+)\s*// ) {
  94             $tempname = $1;
  95             chomp($entry);
  96             $tempdesc = $entry;
  97             if ( defined $name ) {
  98                 $seqchar =~ s/\s//g;
  99                 $seqname = $name;
 100                 $start   = 1;
 101                 $end     = $self->_get_len($seqchar);
 102                 $seq     = Bio::LocatableSeq->new(
 103                     -seq         => $seqchar,
 104                     -display_id  => $seqname,
 105                     -description => $desc,
 106                     -start       => $start,
 107                     -end         => $end,
 108                     -alphabet    => $self->alphabet,
 109                 );
 110                 $aln->add_seq($seq);
 111                 $self->debug("Reading $seqname\n");
 112             }
 113             $desc    = $tempdesc;
 114             $name    = $tempname;
 115             $desc    = $entry;
 116             $seqchar = "";
 117             next;
 118         }
 119
 120         # removed redundant symbol validation
 121         # this is already done in Bio::PrimarySeq
 122         $seqchar .= $entry;
 123     }
 124
 125     #  Next two lines are to silence warnings that
 126     #  otherwise occur at EOF when using <$fh>
 127     $name    = "" if ( !defined $name );
 128     $seqchar = "" if ( !defined $seqchar );
 129     $seqchar =~ s/\s//g;
 130
 131     #  Put away last name and sequence
 132     if ( $name =~ /(\S+\/(\d+)-(\d+))$/ ) {
 133         $seqname = $1;
 134         $start   = $2;
 135         $end     = $3;
 136     }
 137     else {
 138         $seqname = $name;
 139         $start   = 1;
 140         $end     = $self->_get_len($seqchar);
 141     }
 142
 143     # This logic now also reads empty lines at the
 144     # end of the file. Skip this is seqchar and seqname is null
 145     unless ( length($seqchar) == 0 && length($seqname) == 0 ) {
 146         $seq = Bio::LocatableSeq->new(
 147             -seq         => $seqchar,
 148             -display_id  => $seqname,
 149             -description => $desc,
 150             -start       => $start,
 151             -end         => $end,
 152             -alphabet    => $self->alphabet,
 153         );
 154         $aln->add_seq($seq);
 155         $self->debug("Reading $seqname\n");
 156     }
 157     my $alnlen = $aln->length;
 158     foreach my $seq ( $aln->each_seq ) {
 159         if ( $seq->length < $alnlen ) {
 160             my ($diff) = ( $alnlen - $seq->length );
 161             $seq->seq( $seq->seq() . "-" x $diff );
 162         }
 163     }
 164
 165     # no sequences means empty alignment (possible EOF)
 166     return $aln if $aln->num_sequences;
 167     return;
 168 }
 169
 170
 171 =head2 write_aln
 172
 173  Title   : write_aln
 174  Usage   : $stream->write_aln(@aln)
 175  Function: writes the $aln object into the stream in fasta format
 176  Returns : 1 for success and 0 for error
 177  Args    : L<Bio::Align::AlignI> object
 178
 179 See L<Bio::Align::AlignI>
 180
 181 =cut
 182
 183 sub write_aln {
 184     my ($self,@aln) = @_;
 185     my $width = $self->width;
 186     my ($seq,$desc,$rseq,$name,$count,$length,$seqsub);
 187
 188     foreach my $aln (@aln) {
 189         if( ! $aln || ! $aln->isa('Bio::Align::AlignI')  ) {
 190             $self->warn("Must provide a Bio::Align::AlignI object when calling write_aln");
 191             next;
 192         }
 193         if( $self->force_displayname_flat ) {
 194             $aln->set_displayname_flat(1);
 195         }
 196         foreach $rseq ( $aln->each_seq() ) {
 197             $name = $aln->displayname($rseq->get_nse());
 198             $seq  = $rseq->seq();
 199             $desc = $rseq->description || '';
 200             $desc = ' '.$desc if $desc;
 201             $self->_print (">$name$desc\n") or return;
 202             $count = 0;
 203             $length = length($seq);
 204             if(defined $seq && $length > 0) {
 205                 $seq =~ s/(.{1,$width})/$1\n/g;
 206             } else {
 207                 $seq = "\n";
 208             }
 209             $self->_print($seq);
 210         }
 211     }
 212     $self->flush if $self->_flush_on_write && defined $self->_fh;
 213     return 1;
 214 }
 215
 216 =head2 _get_len
 217
 218  Title   : _get_len
 219  Usage   :
 220  Function: determine number of alphabetic chars
 221  Returns : integer
 222  Args    : sequence string
 223
 224 =cut
 225
 226 sub _get_len {
 227     my ($self,$seq) = @_;
 228     my $chars = $Bio::LocatableSeq::GAP_SYMBOLS.$Bio::LocatableSeq::FRAMESHIFT_SYMBOLS;
 229     $seq =~ s{[$chars]+}{}gi;
 230     return CORE::length($seq);
 231 }
 232
 233 =head2 width
 234
 235  Title   : width
 236  Usage   : $obj->width($newwidth)
 237            $width = $obj->width;
 238  Function: Get/set width of alignment
 239  Returns : integer value of width
 240  Args    : on set, new value (a scalar or undef, optional)
 241
 242
 243 =cut
 244
 245 sub width{
 246     my $self = shift;
 247
 248     return $self->{'_width'} = shift if @_;
 249     return $self->{'_width'} || $WIDTH;
 250 }
 251
 252 1;