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