← Back to branch summary

~ubuntu-branches/ubuntu/vivid/bioperl/vivid

« back to all changes in this revision

Viewing changes to Bio/DB/SeqFeature/Store/berkeleydb3.pm

  • Committer: Package Import Robot
  • Author(s): Charles Plessy
  • Date: 2013-09-22 13:39:48 UTC
  • mfrom: (3.1.11 sid)
  • Revision ID: package-import@ubuntu.com-20130922133948-c6z62zegjyp7ztou
Tags: 1.6.922-1
http://bugs.debian.org/722910
* New upstream release.
* Replaces and Breaks grinder (<< 0.5.3-3~) because of overlaping contents.
  Closes: #722910
* Stop Replacing and Breaking bioperl ( << 1.6.9 ): not needed anymore. 

Show diffs side-by-side

added added

removed removed

Lines of Context:
Bio/DB/SeqFeature/Store/berkeleydb3.pm
94
94
    $numeric_cmp->{compare} = sub { $_[0] <=> $_[1] };
95
95
 
96
96
    for my $idx ($self->_index_files) {
97
 
        my $path    = $self->_qualify("$idx.idx");
98
 
        my %db;
99
 
        my $dbtype  = $idx eq 'locations' ? $numeric_cmp
100
 
                     :$idx eq 'summary'   ? $numeric_cmp
 
97
        my $path    = $self->_qualify("$idx.idx");
 
98
        my %db;
 
99
        my $dbtype  = $idx eq 'locations' ? $numeric_cmp
 
100
                     :$idx eq 'summary'   ? $numeric_cmp
101
101
                     :$idx eq 'types'     ? $numeric_cmp
102
102
                     :$idx eq 'seqids'    ? $DB_HASH
103
103
                     :$idx eq 'typeids'   ? $DB_HASH
104
104
                     :$string_cmp;
105
105
 
106
 
        tie(%db,'DB_File',$path,$flags,0666,$dbtype)
107
 
            or $self->throw("Couldn't tie $path: $!");
108
 
        %db = () if $create;
109
 
        $self->index_db($idx=>\%db);
 
106
        tie(%db,'DB_File',$path,$flags,0666,$dbtype)
 
107
            or $self->throw("Couldn't tie $path: $!");
 
108
        %db = () if $create;
 
109
        $self->index_db($idx=>\%db);
110
110
    }
111
111
 
112
112
}
160
160
sub _seq_ids {
161
161
    my $self = shift;
162
162
    if (my $fa = $self->{fasta_db}) {
163
 
        if (my @s = eval {$fa->ids}) {
164
 
            return @s;
165
 
        }
 
163
        if (my @s = eval {$fa->ids}) {
 
164
            return @s;
 
165
        }
166
166
    } 
167
167
    my $l = $self->seqid_db or return;
168
168
    return grep {!/^\./} keys %$l;
188
188
}
189
189
 
190
190
sub _update_type_index {
191
 
  my $self = shift;
192
 
  my ($obj,$id,$delete) = @_;
193
 
  my $db = $self->index_db('types')
194
 
    or $self->throw("Couldn't find 'types' index file");
 
191
    my $self = shift;
 
192
    my ($obj,$id,$delete) = @_;
 
193
    my $db = $self->index_db('types')
 
194
        or $self->throw("Couldn't find 'types' index file");
195
195
 
196
 
  my $key         = $self->_obj_to_type($obj);
197
 
  my $typeid      = $self->add_typeid($key);
198
 
  $self->update_or_delete($delete,$db,$typeid,$id);
 
196
    my $key    = $self->_obj_to_type($obj);
 
197
    my $typeid = $self->add_typeid($key);
 
198
    $self->update_or_delete($delete,$db,$typeid,$id);
199
199
}
200
200
 
201
201
sub _obj_to_type {
212
212
sub types {
213
213
    my $self = shift;
214
214
    eval "require Bio::DB::GFF::Typename" 
215
 
        unless Bio::DB::GFF::Typename->can('new');
 
215
        unless Bio::DB::GFF::Typename->can('new');
216
216
    my $db   = $self->typeid_db;
217
217
    return grep {!/^\./} map {Bio::DB::GFF::Typename->new($_)} keys %$db;
218
218
}
223
223
 
224
224
    my $db = $self->typeid_db;
225
225
    while (my($key,$id) = each %$db) {
226
 
        next if $key =~ /^\./;
227
 
        return $key if $id == $wanted_id;
 
226
        next if $key =~ /^\./;
 
227
        return $key if $id == $wanted_id;
228
228
    }
229
229
    return;
230
230
}
240
240
    my @all_types;
241
241
 
242
242
    for my $type (@types) {
243
 
        my ($primary_tag,$source_tag);
244
 
        if (ref $type && $type->isa('Bio::DB::GFF::Typename')) {
245
 
            $primary_tag = $type->method;
246
 
            $source_tag  = $type->source;
247
 
        } else {
248
 
            ($primary_tag,$source_tag) = split ':',$type,2;
249
 
        }
250
 
        if (defined $source_tag) {
251
 
            my $id = $db->{lc "$primary_tag:$source_tag"};
252
 
            $result{$id}++ if defined $id;
253
 
        } else {
254
 
            @all_types  = $self->types unless @all_types;
255
 
            $result{$db->{$_}}++ foreach grep {/^$primary_tag:/} @all_types;
256
 
        }
 
243
        my ($primary_tag,$source_tag);
 
244
        if (ref $type && $type->isa('Bio::DB::GFF::Typename')) {
 
245
            $primary_tag = $type->method;
 
246
            $source_tag  = $type->source;
 
247
        } else {
 
248
            ($primary_tag,$source_tag) = split ':',$type,2;
 
249
        }
 
250
        if (defined $source_tag) {
 
251
            my $id = $db->{lc "$primary_tag:$source_tag"};
 
252
            $result{$id}++ if defined $id;
 
253
        } else {
 
254
            @all_types  = $self->types unless @all_types;
 
255
            $result{$db->{$_}}++ foreach grep {/^$primary_tag:/} @all_types;
 
256
        }
257
257
    }
258
258
    return \%result;
259
259
}
260
260
 
261
261
sub _update_location_index {
262
 
  my $self = shift;
263
 
  my ($obj,$id,$delete) = @_;
264
 
 
265
 
  my $db = $self->index_db('locations')
266
 
    or $self->throw("Couldn't find 'locations' index file");
267
 
 
268
 
  my $seq_id      = $obj->seq_id || '';
269
 
  my $start       = $obj->start  || '';
270
 
  my $end         = $obj->end    || '';
271
 
  my $strand      = $obj->strand;
272
 
  my $bin_min     = int $start/BINSIZE;
273
 
  my $bin_max     = int $end/BINSIZE;
274
 
 
275
 
  my $typeid      = $self->add_typeid($self->_obj_to_type($obj));
276
 
  my $seq_no      = $self->add_seqid($seq_id);
277
 
 
278
 
  for (my $bin = $bin_min; $bin <= $bin_max; $bin++ ) {
279
 
    my $key = $seq_no * MAX_SEQUENCES + $bin;
280
 
    $self->update_or_delete($delete,$db,$key,pack("i5",$id,$start,$end,$strand,$typeid));
281
 
  }
 
262
    my $self = shift;
 
263
    my ($obj,$id,$delete) = @_;
 
264
 
 
265
    my $db = $self->index_db('locations')
 
266
        or $self->throw("Couldn't find 'locations' index file");
 
267
 
 
268
    my $seq_id      = $obj->seq_id || '';
 
269
    my $start       = $obj->start  || '';
 
270
    my $end         = $obj->end    || '';
 
271
    my $strand      = $obj->strand;
 
272
    my $bin_min     = int $start/BINSIZE;
 
273
    my $bin_max     = int $end/BINSIZE;
 
274
 
 
275
    my $typeid      = $self->add_typeid($self->_obj_to_type($obj));
 
276
    my $seq_no      = $self->add_seqid($seq_id);
 
277
 
 
278
    for (my $bin = $bin_min; $bin <= $bin_max; $bin++ ) {
 
279
        my $key = $seq_no * MAX_SEQUENCES + $bin;
 
280
        $self->update_or_delete($delete,$db,$key,pack("i5",$id,$start,$end,$strand,$typeid));
 
281
    }
282
282
 
283
283
}
284
284
 
285
285
sub _features {
286
 
  my $self = shift;
287
 
  my ($seq_id,$start,$end,$strand,
288
 
      $name,$class,$allow_aliases,
289
 
      $types,
290
 
      $attributes,
291
 
      $range_type,
292
 
      $iterator
293
 
     ) = rearrange([['SEQID','SEQ_ID','REF'],'START',['STOP','END'],'STRAND',
294
 
                    'NAME','CLASS','ALIASES',
295
 
                    ['TYPES','TYPE','PRIMARY_TAG'],
296
 
                    ['ATTRIBUTES','ATTRIBUTE'],
297
 
                    'RANGE_TYPE',
298
 
                    'ITERATOR',
299
 
                   ],@_);
300
 
 
301
 
  my (@from,@where,@args,@group);
302
 
  $range_type ||= 'overlaps';
303
 
 
304
 
  my @result;
305
 
  unless (defined $name or defined $seq_id or defined $types or defined $attributes) {
306
 
      my $is_indexed = $self->index_db('is_indexed');
307
 
      @result = $is_indexed ? grep {$is_indexed->{$_}} keys %{$self->db}
308
 
                            : grep { !/^\./ }keys %{$self->db};
309
 
  }
310
 
 
311
 
  my %found = ();
312
 
  my $result = 1;
313
 
 
314
 
  if (defined($name)) {
315
 
    # hacky backward compatibility workaround
316
 
    undef $class if $class && $class eq 'Sequence';
317
 
    $name     = "$class:$name" if defined $class && length $class > 0;
318
 
    $result &&= $self->filter_by_name($name,$allow_aliases,\%found);
319
 
  }
320
 
 
321
 
  if (defined $seq_id) { # location with or without types
322
 
      my $typelist = defined $types ? $self->_matching_types($types) : undef;
323
 
      $result &&= $self->filter_by_type_and_location($seq_id,$start,$end,$strand,$range_type,
324
 
                                                     $typelist, \%found);
325
 
  }
326
 
 
327
 
  elsif (defined $types) { # types without location
328
 
      $result &&= $self->filter_by_type($types,\%found);
329
 
  }
330
 
 
331
 
  if (defined $attributes) {
332
 
    $result &&= $self->filter_by_attribute($attributes,\%found);
333
 
  }
334
 
 
335
 
  push @result,keys %found if $result;
336
 
  return $iterator ? Bio::DB::SeqFeature::Store::berkeleydb::Iterator->new($self,\@result)
337
 
                   : map {$self->fetch($_)} @result;
 
286
    my $self = shift;
 
287
    my ($seq_id,$start,$end,$strand,
 
288
        $name,$class,$allow_aliases,
 
289
        $types,
 
290
        $attributes,
 
291
        $range_type,
 
292
        $iterator
 
293
       ) = rearrange([['SEQID','SEQ_ID','REF'],'START',['STOP','END'],'STRAND',
 
294
                      'NAME','CLASS','ALIASES',
 
295
                      ['TYPES','TYPE','PRIMARY_TAG'],
 
296
                      ['ATTRIBUTES','ATTRIBUTE'],
 
297
                      'RANGE_TYPE',
 
298
                      'ITERATOR',
 
299
                     ],@_);
 
300
 
 
301
    my (@from,@where,@args,@group);
 
302
    $range_type ||= 'overlaps';
 
303
 
 
304
    my @result;
 
305
    unless (defined $name or defined $seq_id or defined $types or defined $attributes) {
 
306
        my $is_indexed = $self->index_db('is_indexed');
 
307
        @result = $is_indexed ? grep {$is_indexed->{$_}} keys %{$self->db}
 
308
                              : grep { !/^\./ }keys %{$self->db};
 
309
    }
 
310
 
 
311
    my %found = ();
 
312
    my $result = 1;
 
313
 
 
314
    if (defined($name)) {
 
315
        # hacky backward compatibility workaround
 
316
        undef $class if $class && $class eq 'Sequence';
 
317
        $name     = "$class:$name" if defined $class && length $class > 0;
 
318
        $result &&= $self->filter_by_name($name,$allow_aliases,\%found);
 
319
    }
 
320
 
 
321
    if (defined $seq_id) { # location with or without types
 
322
        my $typelist = defined $types ? $self->_matching_types($types) : undef;
 
323
        $result &&= $self->filter_by_type_and_location(
 
324
            $seq_id, $start, $end, $strand, $range_type, $typelist, \%found
 
325
        );
 
326
    }
 
327
 
 
328
    elsif (defined $types) { # types without location
 
329
        $result &&= $self->filter_by_type($types,\%found);
 
330
    }
 
331
 
 
332
    if (defined $attributes) {
 
333
      $result &&= $self->filter_by_attribute($attributes,\%found);
 
334
    }
 
335
 
 
336
    push @result,keys %found if $result;
 
337
    return $iterator ? Bio::DB::SeqFeature::Store::berkeleydb::Iterator->new($self,\@result)
 
338
                     : map {$self->fetch($_)} @result;
338
339
}
339
340
 
340
341
sub filter_by_type {
341
 
  my $self = shift;
342
 
  my ($types,$filter) = @_;
343
 
  my @types = ref $types eq 'ARRAY' ?  @$types : $types;
344
 
 
345
 
  my $index = $self->index_db('types');
346
 
  my $db    = tied(%$index);
347
 
 
348
 
  my @results;
349
 
 
350
 
  for my $type (@types) {
351
 
    my ($primary_tag,$source_tag);
352
 
    if (ref $type && $type->isa('Bio::DB::GFF::Typename')) {
353
 
      $primary_tag = $type->method;
354
 
      $source_tag  = $type->source;
355
 
    } else {
356
 
      ($primary_tag,$source_tag) = split ':',$type,2;
357
 
    }
358
 
    $source_tag ||= '';
359
 
    $primary_tag  = quotemeta($primary_tag);
360
 
    $source_tag    = quotemeta($source_tag);
361
 
    my $match = length $source_tag ? "^$primary_tag:$source_tag\$" : "^$primary_tag:";
362
 
    my $key      = lc "$primary_tag:$source_tag";
363
 
    my $value;
364
 
 
365
 
    # If filter is already provided, then it is usually faster to
366
 
    # fetch each object.
367
 
    if (%$filter) {  
368
 
        for my $id (keys %$filter) {
369
 
            my $obj = $self->_fetch($id) or next;
370
 
            push @results,$id if $obj->type =~ /$match/i;
371
 
        }
372
 
 
373
 
    }
374
 
 
375
 
    else {
376
 
        my $types   = $self->typeid_db;
377
 
        my @typeids = map {$types->{$_}} grep {/$match/} keys %$types;
378
 
        for my $t (@typeids) {
379
 
            my $k = $t;
380
 
            for (my $status = $db->seq($k,$value,R_CURSOR);
381
 
                 $status == 0 && $k == $t;
382
 
                 $status = $db->seq($k,$value,R_NEXT)) {
383
 
                next if %$filter && !$filter->{$value};  # don't even bother
384
 
                push @results,$value;
385
 
            }
386
 
        }
387
 
    }
388
 
  }
389
 
  $self->update_filter($filter,\@results);
 
342
    my $self = shift;
 
343
    my ($types,$filter) = @_;
 
344
    my @types = ref $types eq 'ARRAY' ?  @$types : $types;
 
345
 
 
346
    my $index = $self->index_db('types');
 
347
    my $db    = tied(%$index);
 
348
 
 
349
    my @results;
 
350
 
 
351
    for my $type (@types) {
 
352
        my ($primary_tag,$source_tag);
 
353
        if (ref $type && $type->isa('Bio::DB::GFF::Typename')) {
 
354
            $primary_tag = $type->method;
 
355
            $source_tag  = $type->source;
 
356
        } else {
 
357
            ($primary_tag,$source_tag) = split ':',$type,2;
 
358
        }
 
359
        $source_tag ||= '';
 
360
        $primary_tag  = quotemeta($primary_tag);
 
361
        $source_tag    = quotemeta($source_tag);
 
362
        my $match = length $source_tag ? "^$primary_tag:$source_tag\$" : "^$primary_tag:";
 
363
        my $key      = lc "$primary_tag:$source_tag";
 
364
        my $value;
 
365
 
 
366
        # If filter is already provided, then it is usually faster to
 
367
        # fetch each object.
 
368
        if (%$filter) {  
 
369
            for my $id (keys %$filter) {
 
370
                my $obj = $self->_fetch($id) or next;
 
371
                push @results,$id if $obj->type =~ /$match/i;
 
372
            }
 
373
 
 
374
        }
 
375
 
 
376
        else {
 
377
            my $types   = $self->typeid_db;
 
378
            my @typeids = map {$types->{$_}} grep {/$match/} keys %$types;
 
379
            for my $t (@typeids) {
 
380
                my $k = $t;
 
381
                for (my $status = $db->seq($k,$value,R_CURSOR);
 
382
                    $status == 0 && $k == $t;
 
383
                    $status = $db->seq($k,$value,R_NEXT)) {
 
384
                    next if %$filter && !$filter->{$value};  # don't even bother
 
385
                    push @results,$value;
 
386
                }
 
387
            }
 
388
        }
 
389
    }
 
390
    $self->update_filter($filter,\@results);
390
391
}
391
392
 
392
393
sub filter_by_type_and_location {
393
 
  my $self = shift;
394
 
  my ($seq_id,$start,$end,$strand,$range_type,$typelist,$filter) = @_;
395
 
  $strand ||= 0;
396
 
 
397
 
  my $index    = $self->index_db('locations');
398
 
  my $db       = tied(%$index);
399
 
 
400
 
  my $binstart = defined $start ? int $start/BINSIZE : 0;
401
 
  my $binend   = defined $end   ? int $end/BINSIZE   : MAX_SEQUENCES-1;
402
 
 
403
 
  my %seenit;
404
 
  my @results;
405
 
 
406
 
  $start = MININT  if !defined $start;
407
 
  $end   = MAXINT  if !defined $end;
408
 
 
409
 
  my $seq_no = $self->seqid_id($seq_id);
410
 
  return unless defined $seq_no;
411
 
 
412
 
  if ($range_type eq 'overlaps' or $range_type eq 'contains') {
413
 
    my $keystart = $seq_no * MAX_SEQUENCES + $binstart;
414
 
    my $keystop  = $seq_no * MAX_SEQUENCES + $binend;
415
 
    my $value;
416
 
 
417
 
    for (my $status = $db->seq($keystart,$value,R_CURSOR);
418
 
         $status == 0 && $keystart <= $keystop;
419
 
         $status = $db->seq($keystart,$value,R_NEXT)) {
420
 
      my ($id,$fstart,$fend,$fstrand,$ftype) = unpack("i5",$value);
421
 
      next if $seenit{$id}++;
422
 
      next if $strand   && $fstrand != $strand;
423
 
      next if $typelist && !$typelist->{$ftype};
424
 
      if ($range_type eq 'overlaps') {
425
 
        next unless $fend >= $start && $fstart <= $end;
426
 
      }
427
 
      elsif ($range_type eq 'contains') {
428
 
        next unless $fstart >= $start && $fend <= $end;
429
 
      }
430
 
      next if %$filter && !$filter->{$id};  # don't bother
431
 
      push @results,$id;
432
 
    }
433
 
  }
434
 
 
435
 
  # for contained in, we look for features originating and terminating outside the specified range
436
 
  # this is incredibly inefficient, but fortunately the query is rare (?)
437
 
  elsif ($range_type eq 'contained_in') {
438
 
    my $keystart = $seq_no * MAX_SEQUENCES;
439
 
    my $keystop  = $seq_no * MAX_SEQUENCES + $binstart;
440
 
    my $value;
441
 
 
442
 
    # do the left part of the range
443
 
    for (my $status = $db->seq($keystart,$value,R_CURSOR);
444
 
         $status == 0 && $keystart <= $keystop;
445
 
         $status = $db->seq($keystart,$value,R_NEXT)) {
446
 
      my ($id,$fstart,$fend,$fstrand,$ftype) = unpack("i5",$value);
447
 
      next if $seenit{$id}++;
448
 
      next if $strand && $fstrand != $strand;
449
 
      next if $typelist && !$typelist->{$ftype};
450
 
      next unless $fstart <= $start && $fend >= $end;
451
 
      next if %$filter && !$filter->{$id};  # don't bother
452
 
      push @results,$id;
453
 
    }
454
 
 
455
 
    # do the right part of the range
456
 
    $keystart = $seq_no*MAX_SEQUENCES+$binend;
457
 
    for (my $status = $db->seq($keystart,$value,R_CURSOR);
458
 
         $status == 0;
459
 
         $status = $db->seq($keystart,$value,R_NEXT)) {
460
 
      my ($id,$fstart,$fend,$fstrand,$ftype) = unpack("i5",$value);
461
 
      next if $seenit{$id}++;
462
 
      next if $strand && $fstrand != $strand;
463
 
      next unless $fstart <= $start && $fend >= $end;
464
 
      next if $typelist && !$typelist->{$ftype};
465
 
      next if %$filter && !$filter->{$id};  # don't bother
466
 
      push @results,$id;
467
 
    }
468
 
 
469
 
  }
470
 
 
471
 
  $self->update_filter($filter,\@results);
 
394
    my $self = shift;
 
395
    my ($seq_id,$start,$end,$strand,$range_type,$typelist,$filter) = @_;
 
396
    $strand ||= 0;
 
397
 
 
398
    my $index    = $self->index_db('locations');
 
399
    my $db       = tied(%$index);
 
400
 
 
401
    my $binstart = defined $start ? int $start/BINSIZE : 0;
 
402
    my $binend   = defined $end   ? int $end/BINSIZE   : MAX_SEQUENCES-1;
 
403
 
 
404
    my %seenit;
 
405
    my @results;
 
406
 
 
407
    $start = MININT  if !defined $start;
 
408
    $end   = MAXINT  if !defined $end;
 
409
 
 
410
    my $seq_no = $self->seqid_id($seq_id);
 
411
    return unless defined $seq_no;
 
412
 
 
413
    if ($range_type eq 'overlaps' or $range_type eq 'contains') {
 
414
        my $keystart = $seq_no * MAX_SEQUENCES + $binstart;
 
415
        my $keystop  = $seq_no * MAX_SEQUENCES + $binend;
 
416
        my $value;
 
417
 
 
418
        for (my $status = $db->seq($keystart,$value,R_CURSOR);
 
419
             $status == 0 && $keystart <= $keystop;
 
420
             $status = $db->seq($keystart,$value,R_NEXT)) {
 
421
            my ($id,$fstart,$fend,$fstrand,$ftype) = unpack("i5",$value);
 
422
            next if $seenit{$id}++;
 
423
            next if $strand   && $fstrand != $strand;
 
424
            next if $typelist && !$typelist->{$ftype};
 
425
            if ($range_type eq 'overlaps') {
 
426
                next unless $fend >= $start && $fstart <= $end;
 
427
            }
 
428
            elsif ($range_type eq 'contains') {
 
429
                next unless $fstart >= $start && $fend <= $end;
 
430
            }
 
431
            next if %$filter && !$filter->{$id};  # don't bother
 
432
            push @results,$id;
 
433
        }
 
434
    }
 
435
 
 
436
    # for contained in, we look for features originating and terminating outside the specified range
 
437
    # this is incredibly inefficient, but fortunately the query is rare (?)
 
438
    elsif ($range_type eq 'contained_in') {
 
439
        my $keystart = $seq_no * MAX_SEQUENCES;
 
440
        my $keystop  = $seq_no * MAX_SEQUENCES + $binstart;
 
441
        my $value;
 
442
 
 
443
        # do the left part of the range
 
444
        for (my $status = $db->seq($keystart,$value,R_CURSOR);
 
445
             $status == 0 && $keystart <= $keystop;
 
446
             $status = $db->seq($keystart,$value,R_NEXT)) {
 
447
            my ($id,$fstart,$fend,$fstrand,$ftype) = unpack("i5",$value);
 
448
            next if $seenit{$id}++;
 
449
            next if $strand && $fstrand != $strand;
 
450
            next if $typelist && !$typelist->{$ftype};
 
451
            next unless $fstart <= $start && $fend >= $end;
 
452
            next if %$filter && !$filter->{$id};  # don't bother
 
453
            push @results,$id;
 
454
        }
 
455
 
 
456
        # do the right part of the range
 
457
        $keystart = $seq_no*MAX_SEQUENCES+$binend;
 
458
        for (my $status = $db->seq($keystart,$value,R_CURSOR);
 
459
             $status == 0;
 
460
             $status = $db->seq($keystart,$value,R_NEXT)) {
 
461
            my ($id,$fstart,$fend,$fstrand,$ftype) = unpack("i5",$value);
 
462
            next if $seenit{$id}++;
 
463
            next if $strand && $fstrand != $strand;
 
464
            next unless $fstart <= $start && $fend >= $end;
 
465
            next if $typelist && !$typelist->{$ftype};
 
466
            next if %$filter && !$filter->{$id};  # don't bother
 
467
            push @results,$id;
 
468
        }
 
469
 
 
470
    }
 
471
 
 
472
    $self->update_filter($filter,\@results);
472
473
}
473
474
 
474
475
sub build_summary_statistics {
491
492
    # features by typeid,seqid,start. In the second step, we read through
492
493
    # this sorted list. To avoid running out of memory, we use a db_file
493
494
    # temporary database
494
 
    my $fh   = File::Temp->new() or $self->throw("Couldn't create temporary file for sorting: $!");
 
495
    my $fh   = File::Temp->new() or $self->throw("Could not create temporary file '$name' for sorting: $!");
495
496
    my $name = $fh->filename;
496
497
    my %sort;
497
 
    my $numeric_cmp         = DB_File::BTREEINFO->new;
498
 
    $numeric_cmp->{compare} = sub { $_[0] <=> $_[1] };
499
 
    $numeric_cmp->{flags}   = R_DUP;
500
 
    my $s = tie %sort,'DB_File',$name,O_CREAT|O_RDWR,0666,$numeric_cmp 
501
 
        or $self->throw("Couldn't create temporary file for sorting: $!");
 
498
    my $num_cmp_tre         = DB_File::BTREEINFO->new;
 
499
    $num_cmp_tree->{compare} = sub { $_[0] <=> $_[1] };
 
500
    $num_cmp_tree->{flags}   = R_DUP;
 
501
    my $s = tie %sort, 'DB_File', $name, O_CREAT|O_RDWR, 0666, $num_cmp_tree 
 
502
        or $self->throw("Could not create Berkeley DB in temporary file '$name': $!");
502
503
 
503
504
    my $index    = $self->index_db('locations');
504
505
    my $db       = tied(%$index);
507
508
    my %seenit;
508
509
 
509
510
    for (my $status = $db->seq($keystart,$value,R_CURSOR);
510
 
         $status == 0;
511
 
         $status  = $db->seq($keystart,$value,R_NEXT)) {
512
 
        my ($id,$start,$end,$strand,$typeid) = unpack('i5',$value);
513
 
        next if $seenit{$id}++;
 
511
         $status == 0;
 
512
         $status  = $db->seq($keystart,$value,R_NEXT)) {
 
513
        my ($id,$start,$end,$strand,$typeid) = unpack('i5',$value);
 
514
        next if $seenit{$id}++;
514
515
 
515
 
        print STDERR $count," features sorted$le" if ++$count % 1000 == 0;
516
 
        my $seqid = int($keystart / MAX_SEQUENCES);
517
 
        my $key   = $self->_encode_summary_key($typeid,$seqid,$start-1);
518
 
        $sort{$key}=$end;
 
516
        print STDERR $count," features sorted$le" if ++$count % 1000 == 0;
 
517
        my $seqid = int($keystart / MAX_SEQUENCES);
 
518
        my $key   = $self->_encode_summary_key($typeid,$seqid,$start-1);
 
519
        $sort{$key}=$end;
519
520
    }
520
521
    print STDERR "COUNT = $count\n";
521
522
    
527
528
 
528
529
    # the second step allows us to iterate through this
529
530
    for (my $status = $s->seq($keystart,$end,R_CURSOR);
530
 
         $status == 0;
531
 
         $status  = $s->seq($keystart,$end,R_NEXT)) {
532
 
 
533
 
        print STDERR $count," features processed$le" if ++$count % 1000 == 0;
534
 
        my ($typeid,$seqid,$start) = $self->_decode_summary_key($keystart);
535
 
 
536
 
        my $bin   = int($start/$sbs);
537
 
        
538
 
        # because the input is sorted by start, no more features will contribute to the 
539
 
        # current bin so we can dispose of it
540
 
        if ($bin != $current_bin) {
541
 
            if ($seqid != $current_seqid or $typeid != $current_type) {
542
 
                # load all bins left over
543
 
                $self->_load_bins($insert,\%residuals,\$cum_count,$current_type,$current_seqid);
544
 
                %residuals = () ;
545
 
                $cum_count = 0;
546
 
            } else {
547
 
                # load all up to current one
548
 
                $self->_load_bins($insert,\%residuals,\$cum_count,$current_type,$current_seqid,$current_bin); 
549
 
            }
550
 
        }
551
 
 
552
 
        $last_bin = $current_bin;
553
 
        ($current_seqid,$current_type,$current_bin) = ($seqid,$typeid,$bin);
554
 
 
555
 
        # summarize across entire spanned region
556
 
        my $last_bin = int(($end-1)/$sbs);
557
 
        for (my $b=$bin;$b<=$last_bin;$b++) {
558
 
            $residuals{$b}++;
559
 
        }
 
531
         $status == 0;
 
532
         $status  = $s->seq($keystart,$end,R_NEXT)) {
 
533
 
 
534
        print STDERR $count," features processed$le" if ++$count % 1000 == 0;
 
535
        my ($typeid,$seqid,$start) = $self->_decode_summary_key($keystart);
 
536
 
 
537
        my $bin   = int($start/$sbs);
 
538
        
 
539
        # because the input is sorted by start, no more features will contribute to the 
 
540
        # current bin so we can dispose of it
 
541
        if ($bin != $current_bin) {
 
542
            if ($seqid != $current_seqid or $typeid != $current_type) {
 
543
                # load all bins left over
 
544
                $self->_load_bins($insert,\%residuals,\$cum_count,$current_type,$current_seqid);
 
545
                %residuals = () ;
 
546
                $cum_count = 0;
 
547
            } else {
 
548
                # load all up to current one
 
549
                $self->_load_bins($insert,\%residuals,\$cum_count,$current_type,$current_seqid,$current_bin); 
 
550
            }
 
551
        }
 
552
 
 
553
        $last_bin = $current_bin;
 
554
        ($current_seqid,$current_type,$current_bin) = ($seqid,$typeid,$bin);
 
555
 
 
556
        # summarize across entire spanned region
 
557
        my $last_bin = int(($end-1)/$sbs);
 
558
        for (my $b=$bin;$b<=$last_bin;$b++) {
 
559
            $residuals{$b}++;
 
560
        }
560
561
    }
561
562
 
562
563
    # handle tail case
563
 
    # load all bins left over   
 
564
    # load all bins left over
564
565
    $self->_load_bins($insert,\%residuals,\$cum_count,$current_type,$current_seqid);
565
566
 
566
567
    undef %sort;
571
572
    my $self = shift;
572
573
    my ($insert,$residuals,$cum_count,$typeid,$seqid,$stop_after) = @_;
573
574
    for my $b (sort {$a<=>$b} keys %$residuals) {
574
 
        last if defined $stop_after and $b > $stop_after;
575
 
        $$cum_count += $residuals->{$b};
576
 
        my $key         = $self->_encode_summary_key($typeid,$seqid,$b);
577
 
        $insert->{$key} = $$cum_count;
578
 
        delete $residuals->{$b}; # no longer needed
 
575
        last if defined $stop_after and $b > $stop_after;
 
576
        $$cum_count += $residuals->{$b};
 
577
        my $key         = $self->_encode_summary_key($typeid,$seqid,$b);
 
578
        $insert->{$key} = $$cum_count;
 
579
        delete $residuals->{$b}; # no longer needed
579
580
    }
580
581
}
581
582
 
582
583
sub coverage_array {
583
584
    my $self = shift;
584
585
    my ($seq_name,$start,$end,$types,$bins) = 
585
 
        rearrange([['SEQID','SEQ_ID','REF'],'START',['STOP','END'],
586
 
                   ['TYPES','TYPE','PRIMARY_TAG'],'BINS'],@_);
 
586
        rearrange([['SEQID','SEQ_ID','REF'],'START',['STOP','END'],
 
587
                   ['TYPES','TYPE','PRIMARY_TAG'],'BINS'],@_);
587
588
 
588
589
    $bins  ||= 1000;
589
590
    $start ||= 1;
590
591
    unless ($end) {
591
 
        my $segment = $self->segment($seq_name) or $self->throw("unknown seq_id $seq_name");
592
 
        $end        = $segment->end;
 
592
        my $segment = $self->segment($seq_name) or $self->throw("unknown seq_id $seq_name");
 
593
        $end        = $segment->end;
593
594
    }
594
595
 
595
596
    my $binsize = ($end-$start+1)/$bins;
607
608
 
608
609
    my (%bins,$report_tag);
609
610
    for my $typeid (sort keys %$t) {
610
 
        $report_tag ||= $typeid;
 
611
        $report_tag ||= $typeid;
611
612
 
612
 
        for (my $i=0;$i<@sum_bin_array;$i++) {
613
 
            my $cum_count;
614
 
            my $bin = $sum_bin_array[$i];
615
 
            my $key = $self->_encode_summary_key($typeid,$seqid,$bin);
616
 
            my $status = $db->seq($key,$cum_count,R_CURSOR);
617
 
            next unless $status == 0;
618
 
            push @{$bins{$typeid}},[$bin,$cum_count];
619
 
        }
 
613
        for (my $i=0;$i<@sum_bin_array;$i++) {
 
614
            my $cum_count;
 
615
            my $bin = $sum_bin_array[$i];
 
616
            my $key = $self->_encode_summary_key($typeid,$seqid,$bin);
 
617
            my $status = $db->seq($key,$cum_count,R_CURSOR);
 
618
            next unless $status == 0;
 
619
            push @{$bins{$typeid}},[$bin,$cum_count];
 
620
        }
620
621
    }
621
622
    
622
623
    my @merged_bins;
623
624
    my $firstbin = int(($start-1)/$binsize);
624
625
    for my $type (keys %bins) {
625
 
        my $arry       = $bins{$type};
626
 
        my $last_count = $arry->[0][1]-1;
627
 
        my $last_bin   = -1;
628
 
        my $i          = 0;
629
 
        my $delta;
630
 
        for my $b (@$arry) {
631
 
            my ($bin,$count) = @$b;
632
 
            $delta              = $count - $last_count if $bin > $last_bin;
633
 
            $merged_bins[$i++]  = $delta;
634
 
            $last_count         = $count;
635
 
            $last_bin           = $bin;
636
 
        }
 
626
        my $arry       = $bins{$type};
 
627
        my $last_count = $arry->[0][1]-1;
 
628
        my $last_bin   = -1;
 
629
        my $i          = 0;
 
630
        my $delta;
 
631
        for my $b (@$arry) {
 
632
            my ($bin,$count) = @$b;
 
633
            $delta              = $count - $last_count if $bin > $last_bin;
 
634
            $merged_bins[$i++]  = $delta;
 
635
            $last_count         = $count;
 
636
            $last_bin           = $bin;
 
637
        }
637
638
    }
638
639
    my $returned_type = $self->_id2type($report_tag);
639
640
    return wantarray ? (\@merged_bins,$returned_type) : \@merged_bins;
644
645
    my $self                 = shift;
645
646
    my ($typeid,$seqid,$bin) = @_;
646
647
    $self->throw('Cannot index chromosomes larger than '.C1*SUMMARY_BIN_SIZE/1e6.' megabases')
647
 
        if $bin > C1;
 
648
        if $bin > C1;
648
649
    return ($typeid-1)*C2 + ($seqid-1)*C1 + $bin;
649
650
}
650
651