diff --git a/modules/Bio/EnsEMBL/DBSQL/DensityFeatureAdaptor.pm b/modules/Bio/EnsEMBL/DBSQL/DensityFeatureAdaptor.pm index 91d76b4a1590176885a1af19c9416747e017c2e2..9f3361b5464900dff01dd0685dcb7afa10c8767c 100644 --- a/modules/Bio/EnsEMBL/DBSQL/DensityFeatureAdaptor.pm +++ b/modules/Bio/EnsEMBL/DBSQL/DensityFeatureAdaptor.pm @@ -162,6 +162,7 @@ sub fetch_all_by_Slice { my $density_type = undef; my $best_ratio_large = undef; my $density_type_large = undef; + my %dt_ratio_hash; foreach my $dt (@dtypes) { @@ -175,29 +176,19 @@ sub fetch_all_by_Slice { my $block_size = $slice->seq_region_length() / $dt->region_features(); $ratio = $wanted_block_size / $block_size; } - + # we prefer to use a block size that's smaller than the required one - # (better results on interpolation). remember larger block sizes though - # in case there is no smaller one in the database - if ($ratio < 1) { - if(!defined($best_ratio_large) || $ratio > $best_ratio_large) { - $best_ratio_large = $ratio; - $density_type_large = $dt; - } - } else { - if(!defined($best_ratio) || $ratio < $best_ratio) { - $best_ratio = $ratio; - $density_type = $dt; - } + # (better results on interpolation). + # give larger bits a disadvantage and make them comparable + if( $ratio < 1 ) { + $ratio = 5/$ratio; } - } - # fall back to larger block size if there is no smaller black size - # or it would require retrieving too many features - if( !$best_ratio || $best_ratio > 30 ) { - $best_ratio = $best_ratio_large; - $density_type = $density_type_large; + + $dt_ratio_hash{ $ratio } = $dt; } + $best_ratio = (sort {$a<=>$b} (keys %dt_ratio_hash))[0]; + #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) || @@ -205,6 +196,8 @@ sub fetch_all_by_Slice { return []; } + $density_type = $dt_ratio_hash{$best_ratio}; + my $constraint = "df.density_type_id = " . $density_type->dbID(); my @features = diff --git a/modules/t/densityFeatureAdaptor.t b/modules/t/densityFeatureAdaptor.t index a9f9c39c3205618c6c1cce2e384faf27e23fcb68..0fc98e3af3ba2b7ca22f4c3bae6027bb6c1db9f1 100644 --- a/modules/t/densityFeatureAdaptor.t +++ b/modules/t/densityFeatureAdaptor.t @@ -64,7 +64,6 @@ $dta->store($dt); ok($dt->is_stored($db)); - my $slice_adaptor = $db->get_SliceAdaptor(); my $slice = $slice_adaptor->fetch_by_region('chromosome', '20', 1, ($block_size*10)); @@ -131,16 +130,17 @@ ok($dt->is_stored($db)); @density_features = (); my $chr_20_slice = $slice_adaptor->fetch_by_region( "chromosome", 20 ); -my $step = $chr_20_slice->length() / 300; +my $step = POSIX::ceil( $chr_20_slice->length() / 300); $start = 1; while( $start < $chr_20_slice->length() ) { - my $end = int( $start + $step ); + my $end = $start+$step -1; + if( $end > $chr_20_slice->length ) { $end = $chr_20_slice->length();} push @density_features, Bio::EnsEMBL::DensityFeature->new(-seq_region => $chr_20_slice, - -start => int($start), - -end => int( $start+$step), + -start => $start, + -end => $end, -density_type => $dt, -density_value => 5 ); - $start += $step + 1; + $start += $step; } $dfa->store( @density_features ); @@ -155,6 +155,38 @@ ok( abs( $stored_features[0]->density_value() - 15) < 0.0001 ); debug( "Density value = ".$stored_features[0]->density_value() ); +# test the retreival of the right sized features +# first yet another density size +# comes to about 1000bp on chr 20 + +$dt = Bio::EnsEMBL::DensityType->new(-analysis => $analysis, + -region_features => 30_000, + -value_type => 'sum'); +$dta->store($dt); + +# need some features +@density_features = (); +$step = POSIX::ceil( $chr_20_slice->length() / 30_000); +for my $arr ( [1,1000,1],[1001,2000,2],[2001,3000,3],[31_000_000,31_000_999,31.0],[31_500_000, 31_500_999,31.5] ) { + push @density_features, Bio::EnsEMBL::DensityFeature->new(-seq_region => $chr_20_slice, + -start => $arr->[0], + -end => $arr->[1], + -density_type => $dt, + -density_value => $arr->[2] ); +} +$dfa->store( @density_features ); + +# now check for retrieval +my $sub_Slice = $chr_20_slice->sub_Slice( 1, 1000_000 ); +@stored_features = @{$dfa->fetch_all_by_Slice( $sub_Slice, 'GeneDensityTest', 2 )}; +ok( $stored_features[0]->length() > 10000 ); +@stored_features = @{$dfa->fetch_all_by_Slice( $sub_Slice, 'GeneDensityTest', 10 )}; +ok( $stored_features[0]->length() > 10000 ); +@stored_features = @{$dfa->fetch_all_by_Slice( $sub_Slice, 'GeneDensityTest', 100 )}; +ok( $stored_features[0]->length() == 1000 ); + + + $multi->restore('core', 'analysis');