Bio/Seq/LargeLocatableSeq.pm

   1
   2 # BioPerl module for Bio::Seq::LargeLocatableSeq
   3 #
   4 # Cared for by Albert Vilella
   5 #
   6 #       based on the Bio::LargePrimarySeq module
   7 #       by Ewan Birney <birney@sanger.ac.uk>
   8 #
   9 #       and the Bio::LocatableSeq module
  10 #       by Ewan Birney <birney@sanger.ac.uk>
  11 #
  12 # Copyright Albert Vilella
  13 #
  14 # You may distribute this module under the same terms as perl itself
  15
  16 # POD documentation - main docs before the code
  17
  18 =head1 NAME
  19
  20 Bio::Seq::LargeLocatableSeq - LocatableSeq object that stores sequence as
  21 files in the tempdir
  22
  23 =head1 SYNOPSIS
  24
  25   # normal primary seq usage
  26     use Bio::Seq::LargeLocatableSeq;
  27     my $seq = Bio::Seq::LargeLocatableSeq->new(-seq => "CAGT-GGT",
  28                                               -id  => "seq1",
  29                                               -start => 1,
  30                                               -end   => 7);
  31
  32 =head1 DESCRIPTION
  33
  34 Bio::Seq::LargeLocatableSeq - object with start/end points on it that
  35 can be projected into a MSA or have coordinates relative to another
  36 seq.
  37
  38 This object, unlike Bio::LocatableSeq, stores a sequence as a series
  39 of files in a temporary directory. The aim is to allow someone the
  40 ability to store very large sequences (eg, E<gt> 100MBases) in a file
  41 system without running out of memory (eg, on a 64 MB real memory
  42 machine!).
  43
  44 Of course, to actually make use of this functionality, the programs
  45 which use this object B<must> not call $primary_seq-E<gt>seq otherwise
  46 the entire sequence will come out into memory and probably crash your
  47 machine. However, calls like $primary_seq-E<gt>subseq(10,100) will cause
  48 only 90 characters to be brought into real memory.
  49
  50 =head1 FEEDBACK
  51
  52 =head2 Mailing Lists
  53
  54 User feedback is an integral part of the evolution of this and other
  55 Bioperl modules. Send your comments and suggestions preferably to
  56 the Bioperl mailing list.  Your participation is much appreciated.
  57
  58   bioperl-l@bioperl.org                  - General discussion
  59   http://bioperl.org/wiki/Mailing_lists  - About the mailing lists
  60
  61 =head2 Reporting Bugs
  62
  63 Report bugs to the Bioperl bug tracking system to help us keep track
  64 of the bugs and their resolution. Bug reports can be submitted via
  65 the web:
  66
  67   http://bugzilla.open-bio.org/
  68
  69 =head1 AUTHOR - Albert Vilella
  70
  71 Email avilella-AT-gmail-DOT-com
  72
  73 =head1 APPENDIX
  74
  75 The rest of the documentation details each of the object methods.
  76 Internal methods are usually preceded with a _
  77
  78 =cut
  79
  80
  81 # Let the code begin...
  82
  83 package Bio::Seq::LargeLocatableSeq;
  84 use vars qw($AUTOLOAD);
  85 use strict;
  86
  87
  88 use base qw(Bio::Seq::LargePrimarySeq Bio::LocatableSeq Bio::Root::IO);
  89
  90
  91 =head2 new
  92
  93  Title   : new
  94  Usage   : my $obj = Bio::Seq::LargeLocatableSeq->new();
  95  Function: Builds a new Bio::Seq::LargeLocatableSeq object
  96  Returns : an instance of Bio::Seq::LargeLocatableSeq
  97  Args    :
  98
  99
 100 =cut
 101
 102 sub new {
 103     my ($class, %params) = @_;
 104
 105     # don't let PrimarySeq set seq until we have
 106     # opened filehandle
 107
 108     my $seq = $params{'-seq'} || $params{'-SEQ'};
 109     if($seq ) {
 110         delete $params{'-seq'};
 111         delete $params{'-SEQ'};
 112     }
 113     my $self = $class->SUPER::new(%params);
 114     $self->_initialize_io(%params);
 115     my $tempdir = $self->tempdir( CLEANUP => 1);
 116     my ($tfh,$file) = $self->tempfile( DIR => $tempdir );
 117
 118     $tfh     && $self->_fh($tfh);
 119     $file    && $self->_filename($file);
 120     $self->length(0);
 121     $seq && $self->seq($seq);
 122
 123     return $self;
 124 }
 125
 126
 127 =head2 length
 128
 129  Title   : length
 130  Usage   :
 131  Function:
 132  Example :
 133  Returns :
 134  Args    :
 135
 136
 137 =cut
 138
 139 sub length {
 140    my ($obj,$value) = @_;
 141    if( defined $value) {
 142       $obj->{'length'} = $value;
 143     }
 144
 145     return (defined $obj->{'length'}) ? $obj->{'length'} : 0;
 146 }
 147
 148 =head2 seq
 149
 150  Title   : seq
 151  Usage   :
 152  Function:
 153  Example :
 154  Returns :
 155  Args    :
 156
 157
 158 =cut
 159
 160 sub seq {
 161    my ($self, $data) = @_;
 162    if( defined $data ) {
 163        if( $self->length() == 0) {
 164            $self->add_sequence_as_string($data);
 165        } else {
 166            $self->warn("Trying to reset the seq string, cannot do this with a LargeLocatableSeq - must allocate a new object");
 167        }
 168    }
 169    return $self->subseq(1,$self->length);
 170 }
 171
 172
 173 =head2 subseq
 174
 175  Title   : subseq
 176  Usage   :
 177  Function:
 178  Example :
 179  Returns :
 180  Args    :
 181
 182
 183 =cut
 184
 185 sub subseq{
 186    my ($self,$start,$end) = @_;
 187    my $string;
 188    my $fh = $self->_fh();
 189
 190    if( ref($start) && $start->isa('Bio::LocationI') ) {
 191        my $loc = $start;
 192        if( $loc->length == 0 ) {
 193            $self->warn("Expect location lengths to be > 0");
 194            return '';
 195        } elsif( $loc->end < $loc->start ) {
 196            # what about circular seqs
 197            $self->warn("Expect location start to come before location end");
 198        }
 199        my $seq = '';
 200        if( $loc->isa('Bio::Location::SplitLocationI') ) {
 201            foreach my $subloc ( $loc->sub_Location ) {
 202                if(! seek($fh,$subloc->start() - 1,0)) {
 203                    $self->throw("Unable to seek on file $start:$end $!");
 204                }
 205                my $ret = read($fh, $string, $subloc->length());
 206                if( !defined $ret ) {
 207                    $self->throw("Unable to read $start:$end $!");
 208                }
 209                if( $subloc->strand < 0 ) {
 210 #                  $string = Bio::PrimarySeq->new(-seq => $string)->revcom()->seq();
 211                    $string = Bio::Seq::LargePrimarySeq->new(-seq => $string)->revcom()->seq();
 212                }
 213                $seq .= $string;
 214            }
 215        } else {
 216            if(! seek($fh,$loc->start()-1,0)) {
 217                $self->throw("Unable to seek on file ".$loc->start.":".
 218                             $loc->end ." $!");
 219            }
 220            my $ret = read($fh, $string, $loc->length());
 221            if( !defined $ret ) {
 222                $self->throw("Unable to read ".$loc->start.":".
 223                             $loc->end ." $!");
 224            }
 225            $seq = $string;
 226        }
 227        if( defined $loc->strand &&
 228            $loc->strand < 0 ) {
 229 #          $seq = Bio::PrimarySeq->new(-seq => $string)->revcom()->seq();
 230            $seq = Bio::Seq::LargePrimarySeq->new(-seq => $seq)->revcom()->seq();
 231        }
 232        return $seq;
 233    }
 234    if( $start <= 0 || $end > $self->length ) {
 235        $self->throw("Attempting to get a subseq out of range $start:$end vs ".
 236                     $self->length);
 237    }
 238    if( $end < $start ) {
 239        $self->throw("Attempting to subseq with end ($end) less than start ($start). To revcom use the revcom function with trunc");
 240    }
 241
 242    if(! seek($fh,$start-1,0)) {
 243        $self->throw("Unable to seek on file $start:$end $!");
 244    }
 245    my $ret = read($fh, $string, $end-$start+1);
 246    if( !defined $ret ) {
 247        $self->throw("Unable to read $start:$end $!");
 248    }
 249    return $string;
 250 }
 251
 252
 253 =head2 add_sequence_as_string
 254
 255  Title   : add_sequence_as_string
 256  Usage   : $seq->add_sequence_as_string("CATGAT");
 257  Function: Appends additional residues to an existing LargeLocatableSeq object.
 258            This allows one to build up a large sequence without storing
 259            entire object in memory.
 260  Returns : Current length of sequence
 261  Args    : string to append
 262
 263 =cut
 264
 265 sub add_sequence_as_string{
 266    my ($self,$str) = @_;
 267    my $len = $self->length + CORE::length($str);
 268    my $fh = $self->_fh();
 269    if(! seek($fh,0,2)) {
 270        $self->throw("Unable to seek end of file: $!");
 271    }
 272    $self->_print($str);
 273    $self->length($len);
 274 }
 275
 276
 277 =head2 _filename
 278
 279  Title   : _filename
 280  Usage   : $obj->_filename($newval)
 281  Function:
 282  Example :
 283  Returns : value of _filename
 284  Args    : newvalue (optional)
 285
 286
 287 =cut
 288
 289 sub _filename{
 290    my ($obj,$value) = @_;
 291    if( defined $value) {
 292       $obj->{'_filename'} = $value;
 293     }
 294     return $obj->{'_filename'};
 295
 296 }
 297
 298
 299 =head2 alphabet
 300
 301  Title   : alphabet
 302  Usage   : $obj->alphabet($newval)
 303  Function:
 304  Example :
 305  Returns : value of alphabet
 306  Args    : newvalue (optional)
 307
 308
 309 =cut
 310
 311 sub alphabet{
 312    my ($self,$value) = @_;
 313    if( defined $value) {
 314       $self->SUPER::alphabet($value);
 315     }
 316     return $self->SUPER::alphabet() || 'dna';
 317
 318 }
 319
 320
 321 =head2 end
 322
 323  Title   : end
 324  Usage   : $obj->end($newval)
 325  Function:
 326  Returns : value of end
 327  Args    : newvalue (optional)
 328
 329 =cut
 330
 331 sub end {
 332    my $self = shift;
 333    if( @_ ) {
 334       my $value = shift;
 335       my $string = $self->seq;
 336       if ($self->seq) {
 337           my $len = $self->_ungapped_len;
 338           my $id = $self->id;
 339           $self->warn("In sequence $id residue count gives end value $len.
 340 Overriding value [$value] with value $len for Bio::LargeLocatableSeq::end().")
 341               and $value = $len if $len != $value and $self->verbose > 0;
 342       }
 343
 344       $self->{'end'} = $value;
 345     }
 346
 347    return $self->{'end'} || $self->_ungapped_len;
 348 }
 349
 350
 351 sub DESTROY {
 352     my $self = shift;
 353     my $fh = $self->_fh();
 354     close($fh) if( defined $fh );
 355     # this should be handled by Tempfile removal, but we'll unlink anyways.
 356     unlink $self->_filename() if defined $self->_filename() && -e $self->_filename;
 357     $self->SUPER::DESTROY();
 358 }
 359
 360 1;