Bio/SeqIO/largefasta.pm

   1 # $Id$
   2 # BioPerl module for Bio::SeqIO::largefasta
   3 #
   4 # Cared for by Jason Stajich
   5 #
   6 # Copyright Jason Stajich
   7 #
   8 # You may distribute this module under the same terms as perl itself
   9 # _history
  10 #
  11 # POD documentation - main docs before the code
  12
  13 =head1 NAME
  14
  15 Bio::SeqIO::largefasta - method i/o on very large fasta sequence files
  16
  17 =head1 SYNOPSIS
  18
  19 Do not use this module directly.  Use it via the Bio::SeqIO class.
  20
  21 =head1 DESCRIPTION
  22
  23 This object can transform Bio::Seq objects to and from fasta flat
  24 file databases.
  25
  26 This module handles very large sequence files by using the
  27 Bio::Seq::LargePrimarySeq module to store all the sequence data in
  28 a file.  This can be a problem if you have limited disk space on your
  29 computer because this will effectively cause 2 copies of the sequence
  30 file to reside on disk for the life of the
  31 Bio::Seq::LargePrimarySeq object.  The default location for this is
  32 specified by the L<File::Spec>-E<gt>tmpdir routine which is usually /tmp
  33 on UNIX.  If a sequence file is larger than the swap space (capacity
  34 of the /tmp dir) this could cause problems for the machine.  It is
  35 possible to set the directory where the temporary file is located by
  36 adding the following line to your code BEFORE calling next_seq. See
  37 L<Bio::Seq::LargePrimarySeq> for more information.
  38
  39     $Bio::Seq::LargePrimarySeq::DEFAULT_TEMP_DIR = 'newdir';
  40
  41 =head1 FEEDBACK
  42
  43 =head2 Mailing Lists
  44
  45 User feedback is an integral part of the evolution of this and other
  46 Bioperl modules. Send your comments and suggestions preferably to one
  47 of the Bioperl mailing lists.  Your participation is much appreciated.
  48
  49   bioperl-l@bioperl.org                  - General discussion
  50   http://bioperl.org/wiki/Mailing_lists  - About the mailing lists
  51
  52 =head2 Reporting Bugs
  53
  54 Report bugs to the Bioperl bug tracking system to help us keep track
  55 the bugs and their resolution.  Bug reports can be submitted via the web:
  56
  57   http://bugzilla.open-bio.org/
  58
  59 =head1 AUTHORS - Jason Stajich
  60
  61 Email: jason@bioperl.org
  62
  63 =head1 APPENDIX
  64
  65 The rest of the documentation details each of the object
  66 methods. Internal methods are usually preceded with a _
  67
  68 =cut
  69
  70 # Let the code begin...
  71
  72 package Bio::SeqIO::largefasta;
  73 use vars qw($FASTALINELEN);
  74 use strict;
  75
  76 use Bio::Seq::SeqFactory;
  77
  78 $FASTALINELEN = 60;
  79 use base qw(Bio::SeqIO);
  80
  81 sub _initialize {
  82   my($self,@args) = @_;
  83   $self->SUPER::_initialize(@args);
  84   if( ! defined $self->sequence_factory ) {
  85       $self->sequence_factory(Bio::Seq::SeqFactory->new
  86                               (-verbose => $self->verbose(),
  87                                -type => 'Bio::Seq::LargePrimarySeq'));
  88   }
  89 }
  90
  91 =head2 next_seq
  92
  93  Title   : next_seq
  94  Usage   : $seq = $stream->next_seq()
  95  Function: returns the next sequence in the stream
  96  Returns : Bio::Seq object
  97  Args    : NONE
  98
  99 =cut
 100
 101 sub next_seq {
 102     my ($self) = @_;
 103 #  local $/ = "\n";
 104     my $largeseq = $self->sequence_factory->create();
 105     my ($id,$fulldesc,$entry);
 106     my $count = 0;
 107     my $seen = 0;
 108     while( defined ($entry = $self->_readline) ) {
 109         if( $seen == 1 && $entry =~ /^\s*>/ ) {
 110             $self->_pushback($entry);
 111             return $largeseq;
 112         }
 113 #       if ( ($entry eq '>') || eof($self->_fh) ) { $seen = 1; next; }
 114         if ( ($entry eq '>')  ) { $seen = 1; next; }
 115         elsif( $entry =~ /\s*>(.+?)$/ ) {
 116             $seen = 1;
 117             ($id,$fulldesc) = ($1 =~ /^\s*(\S+)\s*(.*)$/)
 118                 or $self->warn("Can't parse fasta header");
 119             $largeseq->display_id($id);
 120             $largeseq->primary_id($id);
 121             $largeseq->desc($fulldesc);
 122         } else {
 123             $entry =~ s/\s+//g;
 124             $largeseq->add_sequence_as_string($entry);
 125         }
 126         (++$count % 1000 == 0 && $self->verbose() > 0) && print "line $count\n";
 127     }
 128     return unless $seen;
 129     return $largeseq;
 130 }
 131
 132 =head2 write_seq
 133
 134  Title   : write_seq
 135  Usage   : $stream->write_seq(@seq)
 136  Function: writes the $seq object into the stream
 137  Returns : 1 for success and 0 for error
 138  Args    : Bio::Seq object
 139
 140
 141 =cut
 142
 143 sub write_seq {
 144    my ($self,@seq) = @_;
 145    foreach my $seq (@seq) {
 146      my $top = $seq->id();
 147      if ($seq->can('desc') and my $desc = $seq->desc()) {
 148          $desc =~ s/\n//g;
 149          $top .= " $desc";
 150      }
 151      $self->_print (">",$top,"\n");
 152      my $end = $seq->length();
 153      my $start = 1;
 154      while( $start < $end ) {
 155          my $stop = $start + $FASTALINELEN - 1;
 156          $stop = $end if( $stop > $end );
 157          $self->_print($seq->subseq($start,$stop), "\n");
 158          $start += $FASTALINELEN;
 159      }
 160    }
 161
 162    $self->flush if $self->_flush_on_write && defined $self->_fh;
 163    return 1;
 164 }
 165
 166 1;