lib/Bio/SeqIO/fasta.pm

   1 # BioPerl module for Bio::SeqIO::fasta
   2 #
   3 # Please direct questions and support issues to <bioperl-l@bioperl.org>
   4 #
   5 # Cared for by Ewan Birney <birney@ebi.ac.uk>
   6 #          and Lincoln Stein <lstein@cshl.org>
   7 #
   8 # Copyright Ewan Birney & Lincoln Stein
   9 #
  10 # You may distribute this module under the same terms as perl itself
  11 # _history
  12 # October 18, 1999  Largely rewritten by Lincoln Stein
  13
  14 # POD documentation - main docs before the code
  15
  16 =head1 NAME
  17
  18 Bio::SeqIO::fasta - fasta sequence input/output stream
  19
  20 =head1 SYNOPSIS
  21
  22 Do not use this module directly.  Use it via the Bio::SeqIO class.
  23
  24 =head1 DESCRIPTION
  25
  26 This object can transform Bio::Seq objects to and from fasta flat
  27 file databases.
  28
  29 =head1 FEEDBACK
  30
  31 =head2 Mailing Lists
  32
  33 User feedback is an integral part of the evolution of this and other
  34 Bioperl modules. Send your comments and suggestions preferably to one
  35 of the Bioperl mailing lists.  Your participation is much appreciated.
  36
  37   bioperl-l@bioperl.org                  - General discussion
  38   http://bioperl.org/wiki/Mailing_lists  - About the mailing lists
  39
  40 =head2 Support
  41
  42 Please direct usage questions or support issues to the mailing list:
  43
  44 I<bioperl-l@bioperl.org>
  45
  46 rather than to the module maintainer directly. Many experienced and
  47 reponsive experts will be able look at the problem and quickly
  48 address it. Please include a thorough description of the problem
  49 with code and data examples if at all possible.
  50
  51 =head2 Reporting Bugs
  52
  53 Report bugs to the Bioperl bug tracking system to help us keep track
  54 the bugs and their resolution.  Bug reports can be submitted via the
  55 web:
  56
  57   https://github.com/bioperl/bioperl-live/issues
  58
  59 =head1 AUTHORS - Ewan Birney & Lincoln Stein
  60
  61 Email: birney@ebi.ac.uk
  62        lstein@cshl.org
  63
  64 =head1 CONTRIBUTORS
  65
  66 Jason Stajich, jason-at-bioperl.org
  67
  68 =head1 APPENDIX
  69
  70 The rest of the documentation details each of the object
  71 methods. Internal methods are usually preceded with a _
  72
  73 =cut
  74
  75 # Let the code begin...
  76
  77 package Bio::SeqIO::fasta;
  78
  79 use strict;
  80 use warnings;
  81
  82 use Bio::Seq::SeqFastaSpeedFactory;
  83
  84 use parent qw(Bio::SeqIO);
  85
  86 sub _initialize {
  87     my ($self, @args) = @_;
  88     $self->SUPER::_initialize(@args);
  89
  90     ## Initialize fasta specific parameters
  91     ## There are some problems with _rearrange. If there's no value for one of
  92     ## the parameters, it will return an empty value (not undef). This means we
  93     ## can't just merge two hashes since the empty values would override the
  94     ## defaults anyway.
  95     my (%defs) = (
  96                   "width"             => 60,
  97                   "block"             => "", # default is same as width
  98                   "preferred_id_type" => "display",
  99                   );
 100     foreach my $param (keys %defs) {
 101         $self->$param( $self->_rearrange([$param], @args) ||
 102                        $defs{$param});
 103     }
 104
 105     unless ( defined $self->sequence_factory ) {
 106         $self->sequence_factory(Bio::Seq::SeqFastaSpeedFactory->new());
 107     }
 108 }
 109
 110 =head2 next_seq
 111
 112  Title   : next_seq
 113  Usage   : $seq = $stream->next_seq()
 114  Function: returns the next sequence in the stream
 115  Returns : Bio::Seq object, or nothing if no more available
 116  Args    : NONE
 117
 118 =cut
 119
 120 sub next_seq {
 121     my( $self ) = @_;
 122     my $seq;
 123     my $alphabet;
 124     local $/ = "\n>";
 125     return unless my $entry = $self->_readline;
 126
 127     # Replacing chomp for s///, since chomp is not working in some cases
 128     $entry =~ s/\n$//;
 129     $entry =~ s/\r$//;
 130     if ($entry =~ m/\A\s*\Z/s)  { # very first one
 131         return unless $entry = $self->_readline;
 132         chomp($entry);
 133     }
 134
 135     # this just checks the initial input; beyond that, due to setting $/ above,
 136     # the > is part of the record separator and is removed
 137     $self->throw("The sequence does not appear to be FASTA format ".
 138         "(lacks a descriptor line '>')") if $. == 1 && $entry !~ /^>/;
 139
 140     $entry =~ s/^>//;
 141
 142     my ($top,$sequence) = split(/\n/,$entry,2);
 143     defined $sequence && $sequence =~ s/>//g;
 144     #my ($top,$sequence) = $entry =~ /^>?(.+?)\n+([^>]*)/s
 145     #    or $self->throw("Can't parse fasta entry");
 146
 147     my ($id,$fulldesc);
 148     if( $top =~ /^\s*(\S+)\s*(.*)/ ) {
 149         ($id,$fulldesc) = ($1,$2);
 150     }
 151
 152     if (defined $id && $id eq '') {$id=$fulldesc;} # FIX incase no space
 153                                                    # between > and name \AE
 154     defined $sequence && $sequence =~ tr/ \t\n\r//d; # Remove whitespace
 155
 156     # for empty sequences we need to know the mol.type
 157     $alphabet = $self->alphabet();
 158     if(defined $sequence && length($sequence) == 0) {
 159         if(! defined($alphabet)) {
 160             # let's default to dna
 161             $alphabet = "dna";
 162         }
 163     }# else {
 164         # we don't need it really, so disable
 165         # we want to keep this if SeqIO alphabet was set by user
 166         # not sure if this could break something
 167         #$alphabet = undef;
 168     #}
 169
 170     $seq = $self->sequence_factory->create(
 171         -seq         => $sequence,
 172         -id          => $id,
 173         # Ewan's note - I don't think this healthy
 174         # but obviously to taste.
 175         #-primary_id  => $id,
 176         -desc        => $fulldesc,
 177         -alphabet    => $alphabet,
 178         -direct      => 1,
 179     );
 180
 181     # if there wasn't one before, set the guessed type
 182     #unless ( defined $alphabet ) {
 183         # don't assume that all our seqs are the same as the first one found
 184         #$self->alphabet($seq->alphabet());
 185     #}
 186     return $seq;
 187
 188 }
 189 =head2 next_seq_fast
 190
 191  Title   : next_seq_fast
 192  Usage   : $seq = $stream->next_seq_fast()
 193  Function: returns the next sequence in the stream
 194            Favors speed over perfection.
 195  Returns : Bio::Seq object, or nothing if no more available
 196  Args    : NONE
 197
 198 =cut
 199
 200 sub next_seq_fast {
 201     my( $self ) = @_;
 202     my $seq;
 203     my $alphabet;
 204     local $/ = "\n>";
 205     return unless my $entry = $self->_readline;
 206
 207     # Replacing chomp for s///, since chomp is not working in some cases
 208     $entry =~ s/\n(?:\r)$//;
 209     if ($entry =~ m/\A\s*\Z/s)  { # very first one
 210         return unless $entry = $self->_readline;
 211         chomp($entry);
 212     }
 213
 214     # this just checks the initial input; beyond that, due to setting $/ above,
 215     # the > is part of the record separator and is removed
 216     $self->throw("The sequence does not appear to be FASTA format ".
 217         "(lacks a descriptor line '>')") if $. == 1 && $entry !~ /^>/;
 218
 219     $entry =~ s/\s*>\s*//g;
 220     #$entry =~ tr/>//;
 221
 222     my ($top,$sequence) = split(/\n/,$entry,2);
 223     $sequence ||= "";
 224     $sequence =~ s/\s+//g;
 225
 226     my ($id,$fulldesc) = split(/\S+/, $top, 2);
 227     $id       ||= $top;
 228     $fulldesc ||= "";
 229
 230     # for empty sequences we need to know the mol.type
 231     $alphabet = $self->alphabet();
 232     if(defined $sequence && length($sequence) == 0) {
 233         if(! defined($alphabet)) {
 234             # let's default to dna
 235             $alphabet = "dna";
 236         }
 237     }
 238
 239     $seq = $self->sequence_factory->create(
 240         -seq         => $sequence,
 241         -id          => $id,
 242         # Ewan's note - I don't think this healthy
 243         # but obviously to taste.
 244         #-primary_id  => $id,
 245         -desc        => $fulldesc,
 246         -alphabet    => $alphabet,
 247         -direct      => 1,
 248     );
 249
 250     return $seq;
 251
 252 }
 253
 254 =head2 write_seq
 255
 256  Title   : write_seq
 257  Usage   : $stream->write_seq(@seq)
 258  Function: Writes the $seq object into the stream
 259  Returns : 1 for success and 0 for error
 260  Args    : Array of 1 or more Bio::PrimarySeqI objects
 261
 262 =cut
 263
 264 sub write_seq {
 265     my ($self,@seq) = @_;
 266     my $width = $self->width;
 267     my $block = $self->block;
 268
 269     ## take a reference for single string (the sequence) and add the whitespace
 270     local *format_str = sub {
 271         my $str = $_[0];
 272         my @lines = unpack ("(A$width)*", $$str);
 273         if ($block >= $width) {
 274             $$str = join ("\n", @lines)."\n";
 275         } else {
 276             $$str = "";
 277             $$str .= join (" ", unpack ("(A$block)*", $_)) . "\n" foreach (@lines);
 278         }
 279     };
 280
 281     foreach my $seq (@seq) {
 282         $self->throw("Did not provide a valid Bio::PrimarySeqI object")
 283           unless defined $seq && ref($seq) && $seq->isa('Bio::PrimarySeqI');
 284
 285         # Allow for different ids
 286         my $top;
 287         my $id_type = $self->preferred_id_type;
 288         if( $id_type =~ /^acc/i ) {
 289             $top = $seq->accession_number();
 290             if( $id_type =~ /vers/i ) {
 291                 $top .= "." . $seq->version();
 292             }
 293         } elsif($id_type =~ /^displ/i ) {
 294             $self->warn("No whitespace allowed in FASTA ID [". $seq->display_id. "]")
 295               if defined $seq->display_id && $seq->display_id =~ /\s/;
 296             $top = $seq->display_id();
 297             $top = '' unless defined $top;
 298             $self->warn("No whitespace allowed in FASTA ID [". $top. "]")
 299               if defined $top && $top =~ /\s/;
 300         } elsif($id_type =~ /^pri/i ) {
 301             $top = $seq->primary_id();
 302         }
 303
 304         if ($seq->can('desc') and my $desc = $seq->desc()) {
 305             $desc =~ s/\n//g;
 306             $top .= " $desc";
 307         }
 308
 309         if( $seq->isa('Bio::Seq::LargeSeqI') ) {
 310             $self->_print(">$top\n");
 311             # for large seqs, don't call seq(), it defeats the
 312             # purpose of the largeseq functionality.  instead get
 313             # chunks of the seq, $width at a time
 314             my $buff_max = 200000;
 315             my $buff_size = int($buff_max/$width)*$width; #< buffer is even multiple of widths
 316             my $seq_length = $seq->length;
 317             my $num_chunks = int(($seq_length-1)/$buff_size+1);
 318             for( my $c = 0; $c < $num_chunks; $c++ ) {
 319                 my $buff_end = $buff_size*($c+1);
 320                 $buff_end = $seq_length if $buff_end > $seq_length;
 321                 my $buff = $seq->subseq($buff_size*$c+1,$buff_end);
 322                 if($buff) {
 323                     format_str (\$buff);
 324                     $self->_print($buff);
 325                 } else {
 326                     $self->_print("\n");
 327                 }
 328             }
 329         } else {
 330             my $str = $seq->seq;
 331             if(defined $str && length($str) > 0) {
 332                 format_str (\$str);
 333             } else {
 334                 $str = "\n";
 335             }
 336             $self->_print (">",$top,"\n",$str) or return;
 337         }
 338     }
 339
 340     $self->flush if $self->_flush_on_write && defined $self->_fh;
 341     return 1;
 342 }
 343
 344 =head2 width
 345
 346  Title   : width
 347  Usage   : $obj->width($newval)
 348  Function: Get/Set the line width for FASTA output (not counting whitespace).
 349  Returns : value of width
 350  Args    : newvalue (optional)
 351
 352 =cut
 353
 354 sub width {
 355     my ($self,$value) = @_;
 356     if (defined $value) {
 357         $self->{'width'} = $value;
 358     }
 359     return $self->{'width'};
 360 }
 361
 362 =head2 block
 363
 364  Title   : block
 365  Usage   : $obj->block($newval)
 366  Function: Get/Set the length of each block for FASTA output. Sequence blocks
 367            will be split with a space. Configuring block, to a value of 10 for
 368            example, allows one to easily identify a position in a sequence by eye.
 369  Default : same value used for width.
 370  Returns : value of block
 371  Args    : newvalue (optional)
 372
 373 =cut
 374
 375 sub block {
 376     my ($self,$value) = @_;
 377     if (defined $value) {
 378         $self->{'block'} = $value;
 379     }
 380     return $self->{'block'} || $self->width;
 381 }
 382
 383 =head2 preferred_id_type
 384
 385  Title   : preferred_id_type
 386  Usage   : $obj->preferred_id_type('accession')
 387  Function: Get/Set the preferred type of identifier to use in the ">ID" position
 388            for FASTA output.
 389  Returns : string, one of values defined in @Bio::SeqIO::fasta::SEQ_ID_TYPES.
 390  Default : display
 391  Args    : string when setting. This must be one of values defined in
 392            @Bio::SeqIO::fasta::SEQ_ID_TYPES. Allowable values:
 393            accession, accession.version, display, primary
 394  Throws  : fatal exception if the supplied id type is not in @SEQ_ID_TYPES.
 395
 396 =cut
 397
 398 our @SEQ_ID_TYPES = qw(accession accession.version display primary);
 399
 400 sub preferred_id_type {
 401     my ($self,$type) = @_;
 402     if (defined $type) {
 403         if( ! grep lc($type) eq $_, @SEQ_ID_TYPES) {
 404             $self->throw(-class=>'Bio::Root::BadParameter',
 405                          -text=>"Invalid ID type \"$type\". Must be one of: @SEQ_ID_TYPES");
 406         }
 407         $self->{'_seq_id_type'} = lc($type);
 408     }
 409     $self->{'_seq_id_type'};
 410 }
 411
 412 1;