← Back to branch summary

~ubuntu-branches/ubuntu/saucy/bioperl/saucy-proposed

« back to all changes in this revision

Viewing changes to Bio/DB/Fasta.pm

  • Committer: Bazaar Package Importer
  • Author(s): Charles Plessy
  • Date: 2009-03-10 07:19:11 UTC
  • mfrom: (1.2.3 upstream)
  • Revision ID: james.westby@ubuntu.com-20090310071911-fukqzw54pyb1f0bd
Tags: 1.6.0-2
https://launchpad.net/bugs/324001
* Removed patch system (not used):
  - removed instuctions in debian/rules;
  - removed quilt from Build-Depends in debian/control.
* Re-enabled tests:
  - uncommented test command in debian/rules;
  - uncommented previously missing build-dependencies in debian/control.
  - Re-enabled tests and uncommented build-dependencies accordingly.
* Removed libmodule-build-perl and libtest-harness-perl from
  Build-Depends-Indep (provided by perl-modules).
* Better cleaning of empty directories using find -type d -empty -delete
  instead of rmdir in debian/rules (LP: #324001).

Show diffs side-by-side

added added

removed removed

Lines of Context:
Bio/DB/Fasta.pm
1
 
# $Id: Fasta.pm,v 1.44.4.3 2006/10/02 23:10:14 sendu Exp $
 
1
# $Id: Fasta.pm 15364 2009-01-13 17:18:10Z cjfields $
2
2
#
3
3
# BioPerl module for Bio::DB::Fasta
4
4
#
30
30
  my $db      = Bio::DB::Fasta->new('/path/to/fasta/files');
31
31
 
32
32
  my $obj     = $db->get_Seq_by_id('CHROMOSOME_I');
33
 
  my $seq     = $obj->seq;
34
 
  my $subseq  = $obj->subseq(4_000_000 => 4_100_000);
 
33
  my $seq     = $obj->seq; # sequence string
 
34
  my $subseq  = $obj->subseq(4_000_000 => 4_100_000); # string
 
35
  my $trunc   = $obj->trunc(4_000_000 => 4_100_000); # seq object
35
36
  my $length  = $obj->length;
36
37
  # (etc)
37
38
 
75
76
a sequence entry, all lines must be the same length except for the
76
77
last.
77
78
 
 
79
An error will be thrown if this is not the case.
 
80
 
78
81
The module uses /^E<gt>(\S+)/ to extract the primary ID of each sequence 
79
82
from the Fasta header.  During indexing, you may pass a callback routine to
80
83
modify this primary ID.  For example, you may wish to extract a
420
423
use constant DNA     => 1;
421
424
use constant RNA     => 2;
422
425
use constant PROTEIN => 3;
 
426
use constant DIE_ON_MISSMATCHED_LINES => 1; # if you want 
423
427
 
424
428
# Bio::DB-like object
425
429
# providing fast random access to a directory of FASTA files
427
431
=head2 new
428
432
 
429
433
 Title   : new
430
 
 Usage   : my $db = new Bio::DB::Fasta( $path, @options);
 
434
 Usage   : my $db = Bio::DB::Fasta->new( $path, @options);
431
435
 Function: initialize a new Bio::DB::Fasta object
432
436
 Returns : new Bio::DB::Fasta object
433
437
 Args    : path to dir of fasta files or a single filename
481
485
    # that contain whitespace.
482
486
    $path = Win32::GetShortPathName($path)
483
487
      if $^O =~ /^MSWin/i && eval 'use Win32; 1';
484
 
    $offsets = $self->index_dir($path,$opts{-reindex});
 
488
    $offsets = $self->index_dir($path,$opts{-reindex}) or return;
485
489
    $dirname = $path;
486
490
  } elsif (-f _) {
487
491
    $offsets = $self->index_file($path,$opts{-reindex});
519
523
  my %offsets;
520
524
  my $flags = $write ? O_CREAT|O_RDWR : O_RDONLY;
521
525
  my @dbmargs = $self->dbmargs;
522
 
  tie %offsets,'AnyDBM_File',$index,$flags,0644,@dbmargs 
523
 
         or $self->throw( "Can't open cache file $index: $!");
 
526
  eval {
 
527
      tie %offsets,'AnyDBM_File',$index,$flags,0644,@dbmargs 
 
528
          or die "Can't open sequence index file $index: $!";
 
529
  };
 
530
  warn $@ if $@;
524
531
  return \%offsets;
525
532
}
526
533
 
547
554
 
548
555
  # find all fasta files
549
556
  my @files = glob("$dir/$self->{glob}");
550
 
  $self->throw( "no fasta files in $dir") unless @files;
 
557
#  $self->throw( "no fasta files in $dir") unless @files;
 
558
  return unless @files;
551
559
 
552
560
  # get name of index
553
561
  my $index = $self->index_name($dir,1);
732
740
  my $fh = IO::File->new($file) or $self->throw( "Can't open $file: $!");
733
741
  binmode $fh;
734
742
  warn "indexing $file\n" if $self->{debug};
735
 
  my ($offset,$id,$linelength,$type,$firstline,$count,$termination_length,$seq_lines,$last_line,%offsets);
 
743
  my ($offset,@id,$linelength,$type,$firstline,$count,
 
744
      $termination_length,$seq_lines,$last_line,%offsets);
 
745
  my ($l3_len,$l2_len,$l_len)=(0,0,0);
 
746
 
736
747
  while (<$fh>) {               # don't try this at home
737
 
    $termination_length ||= /\r\n$/ ? 2 : 1;  # account for crlf-terminated Windows files
 
748
    $termination_length ||= /\r\n$/ ? 2 : 1; # account for crlf-terminated Windows files
738
749
    if (/^>(\S+)/) {
739
750
      print STDERR "indexed $count sequences...\n" 
740
751
        if $self->{debug} && (++$count%1000) == 0;
741
752
      my $pos = tell($fh);
742
 
      if ($id) {
 
753
      if (@id) {
743
754
        my $seqlength    = $pos - $offset - length($_);
744
755
        $seqlength      -= $termination_length * $seq_lines;
745
 
        $offsets->{$id}  = &{$self->{packmeth}}($offset,$seqlength,
746
 
                                        $linelength,$firstline,
747
 
                                        $type,$base);
 
756
        my $ppos = &{$self->{packmeth}}($offset,$seqlength,
 
757
                                       $linelength,$firstline,
 
758
                                       $type,$base);
 
759
        for my $id (@id) { $offsets->{$id}  = $ppos }
748
760
      }
749
 
      $id = ref($self->{makeid}) eq 'CODE' ? $self->{makeid}->($_) : $1;
 
761
      @id = ref($self->{makeid}) eq 'CODE' ? $self->{makeid}->($_) : $1;
750
762
      ($offset,$firstline,$linelength) = ($pos,length($_),0);
751
763
      $self->_check_linelength($linelength);
 
764
      ($l3_len,$l2_len,$l_len)=(0,0,0);
752
765
      $seq_lines = 0;
753
766
    } else {
 
767
      $l3_len= $l2_len; $l2_len= $l_len; $l_len= length($_); # need to check every line :(
 
768
      if (DIE_ON_MISSMATCHED_LINES &&
 
769
          $l3_len>0 && $l2_len>0 && $l3_len!=$l2_len) {
 
770
        my $fap= substr($_,0,20)."..";
 
771
        $self->throw("Each line of the fasta entry must be the same length except the last.
 
772
    Line above #$. '$fap' is $l2_len != $l3_len chars.");
 
773
  }
754
774
      $linelength ||= length($_);
755
775
      $type       ||= $self->_type($_);
756
776
      $seq_lines++;
760
780
 
761
781
  $self->_check_linelength($linelength);
762
782
  # deal with last entry
763
 
  if ($id) {
 
783
  if (@id) {
764
784
    my $pos = tell($fh);
765
785
    my $seqlength   = $pos - $offset;
766
 
 
767
786
    if ($linelength == 0) { # yet another pesky empty chr_random.fa file
768
787
      $seqlength = 0;
769
788
    } else {
770
789
      if ($last_line !~ /\s$/) {
771
 
        $seq_lines--;
 
790
        $seq_lines--;
772
791
      }
773
792
      $seqlength -= $termination_length * $seq_lines;
774
793
    };
775
 
    $offsets->{$id} = &{$self->{packmeth}}($offset,$seqlength,
776
 
                                   $linelength,$firstline,
777
 
                                   $type,$base);
778
 
}
 
794
    my $ppos = &{$self->{packmeth}}($offset,$seqlength,
 
795
                                           $linelength,$firstline,
 
796
                                           $type,$base);
 
797
    for my $id (@id) { $offsets->{$id}  = $ppos }
 
798
  }
779
799
  $offsets->{__termination_length} = $termination_length;
780
800
  return \%offsets;
781
801
}
1115
1135
 
1116
1136
sub length {
1117
1137
  my $self = shift;
1118
 
  return $self->{db}->length($self->{id});
 
1138
  #return $self->{db}->length($self->{id}); # wrong because ignores sequence start and stop values
 
1139
  return length($self->seq);
 
1140
 
1119
1141
}
1120
1142
 
1121
1143
sub description  {