DensityFeatureAdaptor.pm 24.9 KB
Newer Older
1
=head1 LICENSE
2

3
Copyright [1999-2013] Wellcome Trust Sanger Institute and the EMBL-European Bioinformatics Institute
4

5 6 7 8 9 10 11 12 13 14 15 16 17
Licensed under the Apache License, Version 2.0 (the "License");
you may not use this file except in compliance with the License.
You may obtain a copy of the License at

     http://www.apache.org/licenses/LICENSE-2.0

Unless required by applicable law or agreed to in writing, software
distributed under the License is distributed on an "AS IS" BASIS,
WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
See the License for the specific language governing permissions and
limitations under the License.

=cut
18 19 20 21 22


=head1 CONTACT

  Please email comments or questions to the public Ensembl
Magali Ruffier's avatar
Magali Ruffier committed
23
  developers list at <http://lists.ensembl.org/mailman/listinfo/dev>.
24 25

  Questions may also be sent to the Ensembl help desk at
Magali Ruffier's avatar
Magali Ruffier committed
26
  <http://www.ensembl.org/Help/Contact>.
27 28

=cut
29 30 31 32 33 34 35

=head1 NAME

Bio::EnsEMBL::DBSQL::DensityFeatureAdaptor

=head1 SYNOPSIS

36
  my $dfa = $database_adaptor->get_DensityFeatureAdaptor();
37

38 39
  my $interpolate   = 1;
  my $blocks_wanted = 50;
40

41 42 43 44
  @dense_feats = @{
    $dfa->fetch_all_by_Slice( $slice, 'SNPDensity', $blocks_wanted,
      $interpolate );
    }
45 46 47 48 49 50 51 52 53 54 55 56 57 58

=head1 DESCRIPTION

Density Feature Adaptor - An adaptor responsible for the creation of density
features from the database.

=head1 METHODS

=cut

package Bio::EnsEMBL::DBSQL::DensityFeatureAdaptor;
use vars qw(@ISA);
use strict;

59 60

use POSIX;
61
use Bio::EnsEMBL::DBSQL::BaseFeatureAdaptor;
62
use Bio::EnsEMBL::Utils::Cache;
63
use Bio::EnsEMBL::DensityFeature;
Web Admin's avatar
Web Admin committed
64
use Bio::EnsEMBL::DensityFeatureSet;
65
use Bio::EnsEMBL::DensityType;
66 67 68 69
use Bio::EnsEMBL::Utils::Exception qw(throw warning);

@ISA = qw(Bio::EnsEMBL::DBSQL::BaseFeatureAdaptor);

70
our $DENSITY_FEATURE_CACHE_SIZE = 20;
71

Magali Ruffier's avatar
Magali Ruffier committed
72

73 74 75 76 77 78 79 80 81
=head2 new

  Arg [1]    : list of args @args
               Superclass constructor arguments
  Example    : none
  Description: Constructor which just initializes internal cache structures
  Returntype : Bio::EnsEMBL::DBSQL::DensityFeatureAdaptor
  Exceptions : none
  Caller     : implementing subclass constructors
82
  Status     : Stable
83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99

=cut

sub new {
  my $caller = shift;

  my $class = ref($caller) || $caller;

  my $self = $class->SUPER::new(@_);

  #initialize an LRU cache
  my %cache;
  tie(%cache, 'Bio::EnsEMBL::Utils::Cache', $DENSITY_FEATURE_CACHE_SIZE);
  $self->{'_density_feature_cache'} = \%cache;

  return $self;
}
100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137

=head2 fetch_all_by_Slice

  Arg [1]    : Bio::EnsEMBL::Slice $slice - The slice representing the region
               to retrieve density features from.
  Arg [2]    : string $logic_name - The logic name of the density features to
               retrieve.
  Arg [3]    : int $num_blocks (optional; default = 50) - The number of
               features that are desired. The ratio between the size of these
               features and the size of the features in the database will be
               used to determine which database features will be used.
  Arg [4]    : boolean $interpolate (optional; default = 0) - A flag indicating
               whether the features in the database should be interpolated to
               fit them to the requested number of features.  If true the
               features will be interpolated to provide $num_blocks features.
               This will not guarantee that exactly $num_blocks features are
               returned due to rounding etc. but it will be close.
  Arg [5]    : float $max_ratio - The maximum ratio between the size of the
               requested features (as determined by $num_blocks) and the actual
               size of the features in the database.  If this value is exceeded
               then an empty list will be returned.  This can be used to
               prevent excessive interpolation of the database values.
  Example    : #interpolate:
               $feats = $dfa->fetch_all_by_Slice($slice,'SNPDensity', 10, 1);
               #do not interpoloate, get what is in the database:
               $feats = $dfa->fetch_all_by_Slice($slice,'SNPDensity', 50);
               #interpolate, but not too much
               $feats = $dfa->fetch_all_by_Slice($slice,'SNPDensity',50,1,5.0);
  Description: Retrieves a set of density features which overlap the region
               of this slice. Density features are a discrete representation
               of a continuous value along a sequence, such as a density or
               percent coverage.  Density Features may be stored in chunks of
               different sizes in the database, and interpolated to fit the
               desired size for the requested region.  For example the database
               may store a single density value for each 1KB and also for each
               1MB.  When fetching for an entire chromosome the 1MB density
               chunks will be used if the requested number of blocks is not
               very high.
138 139
               Note that features which have been interpolated are not stored
               in the database and as such will have no dbID or adaptor set.
140 141 142 143 144
  Returntype : Bio::EnsEMBL::DensityFeature
  Exceptions : warning on invalid $num_blocks argument
               warning if no features with logic_name $logic_name exist
               warning if density_type table has invalid block_size value
  Caller     : general
145
  Status     : Stable
146 147 148 149 150 151 152 153 154 155 156 157 158 159

=cut

sub fetch_all_by_Slice {
  my ($self, $slice, $logic_name, $num_blocks, $interpolate, $max_ratio) = @_;

  if(defined($num_blocks) && $num_blocks < 1) {
    warning("Invalid number of density blocks [$num_blocks] requested.\n" .
	   "Returning empty list.");
    return [];
  }

  $num_blocks ||= 50;
  my $length = $slice->length();
160 161

  my $wanted_block_size = POSIX::ceil($length/$num_blocks);
162

163 164 165 166
  #
  # get out all of the density types and choose the one with the
  # block size closest to our desired block size
  #
167

168
  my $dta = $self->db()->get_DensityTypeAdaptor();
169

170
  my @dtypes = @{$dta->fetch_all_by_logic_name($logic_name)};
171
  if( ! @dtypes ){
172 173 174 175 176
    my @all_dtypes =  @{ $dta->fetch_all() };
    @all_dtypes or warning( "No DensityTypes in $dta" ) && return [];
    my $valid_list = join( ", ", map{$_->analysis->logic_name} @all_dtypes );
    warning( "No DensityTypes for logic name $logic_name. ".
	     "Select from $valid_list" );
177 178
    return [];
  }
179

180 181
  my $best_ratio   = undef;
  my $density_type = undef;
182 183
  my $best_ratio_large = undef;
  my $density_type_large = undef;
Arne Stabenau's avatar
Arne Stabenau committed
184
  my %dt_ratio_hash;
185

186
  foreach my $dt (@dtypes) {
187 188 189 190 191 192 193 194 195 196 197

    my $ratio;
    if( $dt->block_size() > 0 ) {
      $ratio = $wanted_block_size/$dt->block_size();
    } else {
      # This is only valid if the requested seq_region is the one the 
      # features are stored on. Please use sensibly. Or find better implementation.

      my $block_size = $slice->seq_region_length() / $dt->region_features();
      $ratio = $wanted_block_size / $block_size;
    }
Arne Stabenau's avatar
Arne Stabenau committed
198
    
199
    # we prefer to use a block size that's smaller than the required one
Arne Stabenau's avatar
Arne Stabenau committed
200 201 202 203
    # (better results on interpolation).
    # give larger bits a disadvantage and make them comparable
    if( $ratio < 1 ) {
      $ratio = 5/$ratio;
204
    }
Arne Stabenau's avatar
Arne Stabenau committed
205 206

    $dt_ratio_hash{ $ratio } = $dt;
207
  }
208

Arne Stabenau's avatar
Arne Stabenau committed
209 210
  $best_ratio = (sort {$a<=>$b} (keys %dt_ratio_hash))[0];
  
211 212 213 214 215 216 217
  #the ratio was not good enough, or this logic name was not in the
  #density type table, return an empty list
  if(!defined($best_ratio) ||
     (defined($max_ratio) && $best_ratio > $max_ratio)) {
    return [];
  }

Arne Stabenau's avatar
Arne Stabenau committed
218 219
  $density_type = $dt_ratio_hash{$best_ratio};

220
  my $constraint = "df.density_type_id = " . $density_type->dbID();
221 222

  my @features =
223
    @{$self->fetch_all_by_Slice_constraint($slice,$constraint)};
224 225 226 227 228 229

  return \@features if(!$interpolate);

  #interpolate the features into new features of a different size
  my @out;
  #sort the features on start position
230
  @features = sort( { $a->start() <=> $b->start() } @features );
231 232 233 234 235

  #resize the features that were returned
  my $start = 1;
  my $end   = $start+$wanted_block_size-1;

236 237 238 239 240 241 242 243 244
  my $value_type = $density_type->value_type();

  # create a new density type object for the interpolated features that
  # is not stored in the database
  $density_type = Bio::EnsEMBL::DensityType->new
    (-VALUE_TYPE => $value_type,
     -ANALYSIS   => $density_type->analysis(),
     -BLOCK_SIZE => $wanted_block_size);

245
  while($start < $length) {
246
#    $end = $length if($end > $length);
247 248 249 250 251 252 253 254 255 256 257 258 259 260

    my $density_value = 0.0;
    my ($f, $fstart, $fend, $portion);
    my @dvalues;

    #construct a new feature using all of the old density features that
    #we overlapped
    while(($f = shift(@features)) && $end > $f->{'start'}) {

      #what portion of this feature are we using to construct our new block?
      $fstart = ($f->{'start'} < $start) ? $start : $f->{'start'};
      $fend   = ($f->{'end'}   > $end  ) ? $end   : $f->{'end'};
      $fend   = $length if($fend > $length);

261
      if($value_type eq 'sum') {
262 263 264 265 266 267 268

        $portion = ($fend-$fstart+1)/$f->length();

        #take a percentage of density value, depending on how much of the
        #feature we overlapped
        $density_value += $portion * $f->{'density_value'};

269
      } elsif($value_type eq 'ratio') {
270 271 272 273 274 275

        #maintain a running total of the length used from each feature
        #and its value
        push(@dvalues,[$fend-$fstart+1,$f->{'density_value'}]);

      } else {
276
        throw("Unknown density value type [$value_type].");
277 278 279 280 281 282 283 284
      }

      #do not want to look at next feature if we only used part of this one:
      last if($fend < $f->{'end'});
    }

    #if we did not completely overlap the last feature, put it back on so
    #it can be partially used by the next block
Graham McVicker's avatar
Graham McVicker committed
285
    if(defined($f) && (!defined($fend) || $fend < $f->{'end'})) {
286 287 288
      unshift(@features, $f);
    }

289
    if($value_type eq 'ratio') {
290 291 292 293 294 295 296 297 298 299 300
      #take a weighted average of the all the density values of the features
      #used to construct this one
      my $total_len = $end - $start + 1;
      if($total_len > 0) {
        foreach my $pair (@dvalues) {
          my ($dlen, $dval) = @$pair;
          $density_value += $dval * ($dlen/$total_len);
        }
      }
    }

301 302
    # interpolated features are not stored in the db so do not set the dbID
    # or the adaptor
303 304 305 306
    push @out, Bio::EnsEMBL::DensityFeature->new
      (-seq_region    => $slice,
       -start         => $start,
       -end           => $end,
307
       -density_type  => $density_type,
308
       -density_value => $density_value);
309 310 311 312 313 314 315 316 317

    $start = $end + 1;
    $end  += $wanted_block_size;
  }

  return \@out;
}


Magali Ruffier's avatar
Magali Ruffier committed
318 319 320 321 322 323 324 325 326 327 328 329 330 331 332 333 334 335 336 337 338 339 340 341 342 343 344 345 346 347 348 349 350 351 352 353 354 355 356 357 358 359
sub fetch_all {
  my $self = shift;
  my $logic_name = shift;

  my ($sth, $seq_region_id, $start, $end, $density_value, $density_type_id);
  my ($slice, $density_type, @out);
  my $sa = $self->db()->get_SliceAdaptor();
  my $dta = $self->db()->get_DensityTypeAdaptor();

  if ($logic_name) {
    $sth = $self->prepare("SELECT df.seq_region_id, df.seq_region_start, df.seq_region_end, df.density_value, df.density_type_id "
                         . "FROM density_feature df, density_type dt, analysis a "
                         . "WHERE df.density_type_id = dt.density_type_id "
                         . "AND dt.analysis_id = a.analysis_id "
                         . "AND logic_name = ?" );
    $sth->bind_param(1, $logic_name, SQL_VARCHAR);
    $sth->execute();
  } else {
    $sth = $self->prepare("SELECT df.seq_region_id, df.seq_region_start, df.seq_region_end, df.density_value, df.density_type_id "
                         . "FROM density_feature df, density_type dt "
                         . "WHERE df.density_type_id = dt.density_type_id" );
    $sth->execute();
  }


  $sth->bind_columns(\($seq_region_id, $start, $end, $density_value, $density_type_id));


  while ($sth->fetch()) {
    $slice = $sa->fetch_by_seq_region_id($seq_region_id);
    $density_type = $dta->fetch_by_dbID($density_type_id);
    push (@out, Bio::EnsEMBL::DensityFeature->new
      (-seq_region    => $slice,
       -start         => $start,
       -end           => $end,
       -density_type  => $density_type,
       -density_value => $density_value));
  }
  return \@out;
}


360 361 362 363 364 365 366 367 368 369 370 371 372 373 374 375 376 377 378 379 380 381 382 383 384 385 386 387 388 389 390 391 392 393 394 395 396 397 398 399 400 401 402 403 404 405 406 407 408 409 410 411 412 413 414 415 416 417 418 419 420 421 422 423 424 425 426 427 428 429 430 431 432 433 434 435 436 437 438 439 440 441 442 443 444 445 446 447 448 449 450
=head2 fetch_all_by_Slice_constraint

  Arg [1]    : Bio::EnsEMBL::Slice $slice
               the slice from which to obtain features
  Arg [2]    : (optional) string $constraint
               An SQL query constraint (i.e. part of the WHERE clause)
  Example    : $fs = $a->fetch_all_by_Slice_constraint($slc, 'density_type_id = 88');
  Description: Returns a listref of features created from the database which 
               are on the Slice defined by $slice and fulfill the SQL 
               constraint defined by $constraint. 
               Note that this is a re-implementation of a method with the same name
               in the BaseFeatureAdaptor.
               This was necessary to remove the use of symliked sequences for density features
  Returntype : listref of Bio::EnsEMBL::SeqFeatures in Slice coordinates
  Exceptions : thrown if $slice is not defined
  Caller     : Bio::EnsEMBL::Slice
  Status     : Stable

=cut


sub fetch_all_by_Slice_constraint {
  my $self = shift;
  my $slice = shift;
  my $constraint = shift;

  my @result;

  #Pull out here as we need to remember them & reset accordingly
  my $bind_params = $self->bind_param_generic_fetch();

  if ( !ref($slice)
       || !(    $slice->isa('Bio::EnsEMBL::Slice')
             or $slice->isa('Bio::EnsEMBL::LRGSlice') ) )
  {
    throw("Bio::EnsEMBL::Slice argument expected.");
  }

  $constraint ||= '';

  # If the logic name was invalid, undef was returned
  if ( !defined($constraint) ) { return [] }

  my $key;
  my $cache;

  # Will only use feature_cache if hasn't been no_cache attribute set
  if (
    !( defined( $self->db()->no_cache() ) && $self->db()->no_cache() ) )
  {

    #strain test and add to constraint if so to stop caching.
    if ( $slice->isa('Bio::EnsEMBL::StrainSlice') ) {
      my $string =
        $self->dbc()->db_handle()->quote( $slice->strain_name() );

      if ( $constraint ne "" ) {
        $constraint .= " AND $string = $string ";
      } else {
        $constraint .= " $string = $string ";
      }
    }

    # Check the cache and return the cached results if we have already
    # done this query.  The cache key is the made up from the slice
    # name, the constraint, and the bound parameters (if there are any).
    $key = uc( join( ':', $slice->name(), $constraint ) );

    if ( defined($bind_params) ) {
      $key .= ':'
        . join( ':', map { $_->[0] . '/' . $_->[1] } @{$bind_params} );
    }

    $cache = $self->_slice_feature_cache();
    if ( exists( $cache->{$key} ) ) {
      # Clear the bound parameters and return the cached data.
      $self->{'_bind_param_generic_fetch'} = ();
      return $cache->{$key};
    }
  } ## end if ( !( defined( $self...)))

  $self->_bind_param_generic_fetch($bind_params);
  my $features = $self->_slice_fetch( $slice, $constraint );

  if ( defined($key) ) {
    $cache->{$key} = $features;
  }

  return $features;
}

Magali Ruffier's avatar
Magali Ruffier committed
451

452 453 454
sub _tables {
  my $self = shift;

455
  return (['density_feature', 'df']);
456 457 458 459 460 461 462 463
}


sub _columns {
  my $self = shift;

  return qw( df.density_feature_id
             df.seq_region_id df.seq_region_start df.seq_region_end
464
             df.density_value df.density_type_id );
465 466 467 468 469 470 471 472 473 474 475 476 477 478
}



sub _objs_from_sth {
  my ($self, $sth, $mapper, $dest_slice) = @_;

  #
  # This code is ugly because an attempt has been made to remove as many
  # function calls as possible for speed purposes.  Thus many caches,
  # and a fair bit of gymnastics have been used.
  #

  my $sa = $self->db()->get_SliceAdaptor();
479
  my $dta = $self->db()->get_DensityTypeAdaptor();
480 481

  my @features;
482
  my %dtype_hash;
483 484 485 486 487
  my %slice_hash;
  my %sr_name_hash;
  my %sr_cs_hash;

  my($density_feature_id, $seq_region_id, $seq_region_start, $seq_region_end,
488
     $density_value, $density_type_id );
489 490

  $sth->bind_columns(\$density_feature_id, \$seq_region_id, \$seq_region_start,
491
                     \$seq_region_end, \$density_value, \$density_type_id);
492 493 494 495 496 497 498 499 500 501 502 503 504 505 506 507 508 509 510 511

  my $asm_cs;
  my $cmp_cs;
  my $asm_cs_vers;
  my $asm_cs_name;
  my $cmp_cs_vers;
  my $cmp_cs_name;
  if($mapper) {
    $asm_cs = $mapper->assembled_CoordSystem();
    $cmp_cs = $mapper->component_CoordSystem();
    $asm_cs_name = $asm_cs->name();
    $asm_cs_vers = $asm_cs->version();
    $cmp_cs_name = $cmp_cs->name();
    $cmp_cs_vers = $cmp_cs->version();
  }

  my $dest_slice_start;
  my $dest_slice_end;
  my $dest_slice_strand;
  my $dest_slice_length;
512
  my $dest_slice_sr_name;
513
  my $dest_slice_sr_id;
514

515 516 517 518 519
  if($dest_slice) {
    $dest_slice_start  = $dest_slice->start();
    $dest_slice_end    = $dest_slice->end();
    $dest_slice_strand = $dest_slice->strand();
    $dest_slice_length = $dest_slice->length();
520
    $dest_slice_sr_name = $dest_slice->seq_region_name();
521
    $dest_slice_sr_id  = $dest_slice->get_seq_region_id();
522 523 524
  }

  FEATURE: while($sth->fetch()) {
525 526 527
    #get the density type object
    my $density_type = $dtype_hash{$density_type_id} ||=
      $dta->fetch_by_dbID($density_type_id);
528 529

    #get the slice object
530 531
    #need to get the internal_seq_region, if present
    $seq_region_id = $self->get_seq_region_id_internal($seq_region_id);
532 533 534 535 536 537 538 539
    my $slice = $slice_hash{"ID:".$seq_region_id};
    if(!$slice) {
      $slice = $sa->fetch_by_seq_region_id($seq_region_id);
      $slice_hash{"ID:".$seq_region_id} = $slice;
      $sr_name_hash{$seq_region_id} = $slice->seq_region_name();
      $sr_cs_hash{$seq_region_id} = $slice->coord_system();
    }

540 541 542 543
    my $sr_name = $sr_name_hash{$seq_region_id};
    my $sr_cs   = $sr_cs_hash{$seq_region_id};
 
   #
544 545 546 547 548 549 550
    # remap the feature coordinates to another coord system
    # if a mapper was provided
    #
    if($mapper) {

      my $len = $seq_region_end - $seq_region_start + 1;

551 552 553 554 555 556 557 558 559 560
      my @coords;

      if (defined $dest_slice && $mapper->isa('Bio::EnsEMBL::ChainedAssemblyMapper')  ) {
	    
	  @coords = $mapper->map( $sr_name, $seq_region_start, $seq_region_end,
                          1, $sr_cs, 0, $dest_slice);

      } else {
	  @coords = $mapper->map($sr_name, $seq_region_start, $seq_region_end,1, $sr_cs);
      }
561 562 563 564 565 566 567 568 569

      #filter out gaps
      @coords = grep {!$_->isa('Bio::EnsEMBL::Mapper::Gap')} @coords;

      #throw out density features mapped to gaps, or split
      next FEATURE if(@coords != 1);

      $seq_region_start  = $coords[0]->{'start'};
      $seq_region_end    = $coords[0]->{'end'};
570
      $seq_region_id     = $coords[0]->{'id'};
571

572
      if($density_type->value_type() eq 'sum') {
573 574 575 576 577 578
        #adjust density value so it reflects length of feature actually used
        my $newlen = $seq_region_end - $seq_region_start +1;
        $density_value *= $newlen/$len if($newlen != $len);
      }

      #get a slice in the coord system we just mapped to
579 580 581 582 583 584 585 586
#      if($asm_cs == $sr_cs || ($cmp_cs != $sr_cs && $asm_cs->equals($sr_cs))) {
        $slice = $slice_hash{"ID:".$seq_region_id} ||=
          $sa->fetch_by_seq_region_id($seq_region_id);
#      } else {
#        $slice = $slice_hash{"NAME:$sr_name:$asm_cs_name:$asm_cs_vers"} ||=
#          $sa->fetch_by_region($asm_cs_name, $sr_name, undef, undef, undef,
#                               $asm_cs_vers);
#      }
587 588 589 590 591 592 593 594 595 596 597 598 599 600 601 602
    }

    #
    # If a destination slice was provided convert the coords
    # If the dest_slice starts at 1 and is foward strand, nothing needs doing
    #
    if($dest_slice) {
      if($dest_slice_start != 1 || $dest_slice_strand != 1) {
        if($dest_slice_strand == 1) {
          $seq_region_start = $seq_region_start - $dest_slice_start + 1;
          $seq_region_end   = $seq_region_end   - $dest_slice_start + 1;
        } else {
          my $tmp_seq_region_start = $seq_region_start;
          $seq_region_start = $dest_slice_end - $seq_region_end + 1;
          $seq_region_end   = $dest_slice_end - $tmp_seq_region_start + 1;
        }
603
      }
604

605 606
      #throw away features entirely off the end of the requested slice
      if($seq_region_end < 1 || $seq_region_start > $dest_slice_length ||
607
	 ( $dest_slice_sr_id ne $seq_region_id )) {
608
	next FEATURE;
609 610 611 612
      }
      $slice = $dest_slice;
    }

613 614 615 616 617 618 619 620 621 622 623
    push( @features,
          $self->_create_feature( 'Bio::EnsEMBL::DensityFeature', {
                                    -dbID       => $density_feature_id,
                                    -adaptor    => $self,
                                    -start      => $seq_region_start,
                                    -end        => $seq_region_end,
                                    -seq_region => $slice,
                                    -density_value => $density_value,
                                    -density_type  => $density_type
                                  } ) );

624
  }
625

626 627 628 629 630 631 632 633 634 635
  return \@features;
}


=head2 list_dbIDs

  Arg [1]    : none
  Example    : @feature_ids = @{$density_feature_adaptor->list_dbIDs()};
  Description: Gets an array of internal ids for all density features in the
               current db
636
  Arg[1]     : <optional> int. not set to 0 for the ids to be sorted by the seq_region.
637 638 639
  Returntype : list of ints
  Exceptions : none
  Caller     : ?
640
  Status     : Stable
641 642 643 644

=cut

sub list_dbIDs {
645
   my ($self, $ordered) = @_;
646

647
   return $self->_list_dbIDs("density_feature",undef, $ordered);
648 649
}

650 651 652 653 654 655 656 657 658 659 660 661

=head2 store

  Arg [1]    : list of Bio::EnsEMBL::DensityFeatures @df
               the simple features to store in the database
  Example    : $density_feature_adaptor->store(1234, @density_feats);
  Description: Stores a list of density feature objects in the database
  Returntype : none
  Exceptions : thrown if @df is not defined, if any of the features do not
               have an attached slice.
               or if any elements of @df are not Bio::EnsEMBL::SeqFeatures 
  Caller     : general
662
  Status     : Stable
663 664 665 666 667 668 669 670 671 672 673 674 675 676 677 678 679 680 681 682 683 684 685 686 687 688 689 690 691 692 693 694 695 696 697 698

=cut

sub store{
  my ($self,@df) = @_;

  if( scalar(@df) == 0 ) {
    throw("Must call store with list of DensityFeatures");
  }
#mysql> desc density_feature;
#+--------------------+---------+------+-----+---------+----------------+
#| Field              | Type    | Null | Key | Default | Extra          |
#+--------------------+---------+------+-----+---------+----------------+
#| density_feature_id | int(11) |      | PRI | NULL    | auto_increment |
#| density_type_id    | int(11) |      | MUL | 0       |                |
#| seq_region_id      | int(11) |      |     | 0       |                |
#| seq_region_start   | int(11) |      |     | 0       |                |
#| seq_region_end     | int(11) |      |     | 0       |                |
#| density_value      | float   |      |     | 0       |                |
#+--------------------+---------+------+-----+---------+----------------+

  my $sth = $self->prepare
    ("INSERT INTO density_feature (seq_region_id, seq_region_start, " .
                                  "seq_region_end, density_type_id, " .
                                  "density_value) " .
     "VALUES (?,?,?,?,?)");

  my $db = $self->db();
  my $analysis_adaptor = $db->get_AnalysisAdaptor();

 FEATURE: foreach my $df ( @df ) {

    if( !ref $df || !$df->isa("Bio::EnsEMBL::DensityFeature") ) {
      throw("DensityFeature must be an Ensembl DensityFeature, " .
            "not a [".ref($df)."]");
    }
699 700 701
    
    # we dont store 0 value density features
    next if( $df->density_value == 0 );
702 703 704 705 706 707 708 709 710 711 712
    if($df->is_stored($db)) {
      warning("DensityFeature [".$df->dbID."] is already stored" .
              " in this database.");
      next FEATURE;
    }

    if(!defined($df->density_type)) {
      throw("A density type must be attached to the features to be stored.");
    }

    #store the density_type if it has not been stored yet
713

714
    if(!$df->density_type->is_stored($db)) {
715 716
      my $dta = $db->get_DensityTypeAdaptor();
      $dta->store($df->density_type());
717 718 719 720 721 722
    }

    my $original = $df;
    my $seq_region_id;
    ($df, $seq_region_id) = $self->_pre_store($df);

723 724 725 726 727 728
    $sth->bind_param(1,$seq_region_id,SQL_INTEGER);
    $sth->bind_param(2,$df->start,SQL_INTEGER);
    $sth->bind_param(3,$df->end,SQL_INTEGER);
    $sth->bind_param(4,$df->density_type->dbID,SQL_INTEGER);
    $sth->bind_param(5,$df->density_value,SQL_FLOAT);
    $sth->execute();
729 730 731 732 733 734

    $original->dbID($sth->{'mysql_insertid'});
    $original->adaptor($self);
  }
}

735 736 737 738 739 740 741 742 743 744 745 746 747
=head2 fetch_Featureset_by_Slice 

  Arg [1-5]  : see
               Bio::EnsEMBL::DBSQL::DensityFeatureAdaptor::fetch_all_by_Slice()
               for argument documentation
  Example    : $featureset = $dfa->fetch_FeatureSet_by_Slice($slice,'SNPDensity', 10, 1);
  Description: wrapper around
               Bio::EnsEMBL::DBSQL::DensityFeatureAdaptor::fetch_all_by_Slice()
               which returns a Bio::EnsEMBL::DensityFeatureSet and also caches
               results
  Returntype : Bio::EnsEMBL::DensityFeatureSet
  Exceptions : none
  Caller     : general
748
  Status     : Stable
749 750

=cut
751

Web Admin's avatar
Web Admin committed
752
sub fetch_Featureset_by_Slice {
753 754 755 756 757 758 759 760 761 762 763 764 765 766 767 768 769
  my ($self, $slice, $logic_name, $num_blocks, $interpolate, $max_ratio) = @_;

  my $key = join(":", $slice->name,
                      $logic_name,
                      $num_blocks || 50,
                      $interpolate || 0,
                      $max_ratio);
                      
  unless ($self->{'_density_feature_cache'}->{$key}) {
    my $dfeats = $self->fetch_all_by_Slice($slice, $logic_name, $num_blocks,
                                           $interpolate, $max_ratio);
    $self->{'_density_feature_cache'}->{$key} =
      new Bio::EnsEMBL::DensityFeatureSet($dfeats);
  }
  return $self->{'_density_feature_cache'}->{$key};
}

770

771
1;