BaseFeatureAdaptor.pm 31.5 KB
Newer Older
1
#
Simon Potter's avatar
Simon Potter committed
2
# EnsEMBL module for Bio::EnsEMBL::DBSQL::BaseFeatureAdaptor
3
#
Graham McVicker's avatar
Graham McVicker committed
4
# Copyright (c) 2003 Ensembl
5 6 7 8 9
#
# You may distribute this module under the same terms as perl itself

=head1 NAME

10
Bio::EnsEMBL::DBSQL::BaseFeatureAdaptor - An Abstract Base class for all
Graham McVicker's avatar
Graham McVicker committed
11
FeatureAdaptors
12 13 14

=head1 SYNOPSIS

15
Abstract class - should not be instantiated.  Implementation of
16 17 18 19 20
abstract methods must be performed by subclasses.

=head1 DESCRIPTION

This is a base adaptor for feature adaptors. This base class is simply a way
21
of eliminating code duplication through the implementation of methods
22 23
common to all feature adaptors.

24 25
=head1 CONTACT

Graham McVicker's avatar
Graham McVicker committed
26
Contact Ensembl development list for info: <ensembl-dev@ebi.ac.uk>
27 28 29

=head1 METHODS

30 31 32
=cut

package Bio::EnsEMBL::DBSQL::BaseFeatureAdaptor;
33
use vars qw(@ISA);
34 35 36
use strict;

use Bio::EnsEMBL::DBSQL::BaseAdaptor;
37
use Bio::EnsEMBL::Utils::Cache;
38 39
use Bio::EnsEMBL::Utils::Exception qw(warning throw deprecate);
use Bio::EnsEMBL::Utils::Argument qw(rearrange);
40 41 42

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

43 44
my $SLICE_FEATURE_CACHE_SIZE = 4;
my $MAX_SPLIT_QUERY_SEQ_REGIONS = 3;
45 46 47 48 49 50 51 52 53 54 55 56 57 58

=head2 new

  Arg [1]    : list of args @args
               Superclass constructor arguments
  Example    : none
  Description: Constructor which just initializes internal cache structures
  Returntype : Bio::EnsEMBL::BaseFeatureAdaptor
  Exceptions : none
  Caller     : implementing subclass constructors

=cut

sub new {
59
  my $caller = shift;
60

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

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

65 66 67 68
  #initialize an LRU cache
  my %cache;
  tie(%cache, 'Bio::EnsEMBL::Utils::Cache', $SLICE_FEATURE_CACHE_SIZE);
  $self->{'_slice_feature_cache'} = \%cache;
69 70 71

  return $self;
}
72

73

74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98

=head2 _straight_join

  Arg [1]    : (optional) boolean $new_val
  Example    : $self->_straight_join(1);
               $self->generic_fetch($constraint);
               $self->_straight_join(0);
  Description: Getter/Setter that turns on/off the use of a straight join
               in queries.
  Returntype : boolean
  Exceptions : none
  Caller     : general

=cut

sub _straight_join {
  my $self = shift;
  if(@_) {
    $self->{'_straight_join'} = shift;
  }

  return $self->{'_straight_join'};
}


99 100 101 102
=head2 generic_fetch

  Arg [1]    : (optional) string $constraint
               An SQL query constraint (i.e. part of the WHERE clause)
103 104 105 106 107 108 109 110
  Arg [2]    : (optional) Bio::EnsEMBL::AssemblyMapper $mapper
               A mapper object used to remap features
               as they are retrieved from the database
  Arg [3]    : (optional) Bio::EnsEMBL::Slice $slice
               A slice that features should be remapped to
  Arg [4]    : (optional) boolean $keep_all
               Set to 1 if all features, even ones entirely off slice,
               should be kept
111
  Example    : $fts = $a->generic_fetch('contig_id in (1234, 1235)', 'Swall');
112 113
  Description: Performs a database fetch and returns feature objects in
               contig coordinates.
114
  Returntype : listref of Bio::EnsEMBL::SeqFeature in contig coordinates
115 116 117 118
  Exceptions : none
  Caller     : BaseFeatureAdaptor, ProxyDnaAlignFeatureAdaptor::generic_fetch

=cut
119

120
sub generic_fetch {
121 122
  my ($self, $constraint, $mapper, $slice, $keep_all) = @_;

123
  my @tabs = $self->_tables;
124
  my $columns = join(', ', $self->_columns());
125

126
  my $db = $self->db();
127 128 129 130 131

  #
  # Construct a left join statement if one was defined, and remove the
  # left-joined table from the table list
  #
132
  my @left_join_list = $self->_left_join();
133 134
  my $left_join = '';
  my @tables;
135 136
  if(@left_join_list) {
    my %left_join_hash = map { $_->[0] => $_->[1] } @left_join_list;
137
    while(my $t = shift @tabs) {
138 139 140
      if( exists $left_join_hash{ $t->[0] } ) {
        my $condition = $left_join_hash{ $t->[0] };
        my $syn = $t->[1];
Web Admin's avatar
Web Admin committed
141
        $left_join .=  "LEFT JOIN\n       ".$t->[0]." $syn ON $condition ";
142
      } else {
143
        push @tables, $t;
144 145 146 147 148
      }
    }
  } else {
    @tables = @tabs;
  }
149

150 151 152 153 154 155 156

  my $straight_join = '';

  if($self->_straight_join()) {
    $straight_join = "STRAIGHT_JOIN";
  }

157 158 159
  #construct a nice table string like 'table1 t1, table2 t2'
  my $tablenames = join(', ', map({ join(' ', @$_) } @tables));

160
  my $sql = "SELECT $straight_join $columns\n  FROM $tablenames $left_join";
161 162 163 164 165

  my $default_where = $self->_default_where_clause;
  my $final_clause = $self->_final_clause;

  #append a where clause if it was defined
166
  if($constraint) {
Web Admin's avatar
Web Admin committed
167
    $sql .= "\n WHERE $constraint ";
168
    if($default_where) {
Web Admin's avatar
Web Admin committed
169
      $sql .= " AND\n       $default_where ";
170 171
    }
  } elsif($default_where) {
Web Admin's avatar
Web Admin committed
172
    $sql .= "\n WHERE $default_where ";
173 174 175
  }

  #append additional clauses which may have been defined
Web Admin's avatar
Web Admin committed
176
  $sql .= "\n$final_clause";
177

178
  ##print STDERR "\n\n$sql\n\n";
179

180 181 182 183 184
  my $sth = $db->prepare($sql);

  $sth->execute;

  my $res = $self->_objs_from_sth($sth, $mapper, $slice, $keep_all);
185

186
  return $res;
187 188 189 190
}


=head2 fetch_by_dbID
191 192

  Arg [1]    : int $id
193 194 195
               The unique database identifier for the feature to be obtained
  Example    : $feat = $adaptor->fetch_by_dbID(1234));
               $feat = $feat->transform('contig');
196
  Description: Returns the feature created from the database defined by the
197 198 199 200 201 202 203 204 205
               the id $id.  The feature will be returned in its native
               coordinate system.  That is, the coordinate system in which it
               is stored in the database.  In order to convert it to a
               particular coordinate system use the transfer() or transform()
               method.  If the feature is not found in the database then
               undef is returned instead
  Returntype : Bio::EnsEMBL::Feature or undef
  Exceptions : thrown if $id arg is not provided
               does not exist
206
  Caller     : general
207 208 209 210 211 212

=cut

sub fetch_by_dbID{
  my ($self,$id) = @_;

213
  throw("id argument is required") if(!defined $id);
214

215 216
  #construct a constraint like 't1.table1_id = 123'
  my @tabs = $self->_tables;
217 218
  my ($name, $syn) = @{$tabs[0]};
  my $constraint = "${syn}.${name}_id = $id";
219

220 221 222 223 224
  #Should only be one
  my ($feat) = @{$self->generic_fetch($constraint)};

  return undef if(!$feat);

225
  return $feat;
226 227 228
}


229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244
=head2 fetch_all_by_dbID_list

  Arg [1]    : listref of ints $id_list
               The unique database identifiers for the features to be obtained
  Example    : @feats = @{$adaptor->fetch_by_dbID_list([1234, 2131, 982]))};
  Description: Returns the features created from the database defined by the
               the ids in contained in the id list $id_list.  The features 
               will be returned in their native coordinate system. That is, 
               the coordinate system in which they are stored in the database.
               In order to convert the features to a particular coordinate 
               system use the transfer() or transform() method.  If none of the
               features are found in the database a reference to an empty 
               list is returned.
  Returntype : listref of Bio::EnsEMBL::Features
  Exceptions : thrown if $id arg is not provided
               does not exist
245
  Caller     : general
246 247 248

=cut

249 250
sub fetch_all_by_dbID_list {
  my ($self,$id_list) = @_;
251

252 253
  if(!defined($id_list) || ref($id_list) ne 'ARRAY') {
    throw("id_list list reference argument is required") if(!defined $id_list);
254 255
  }

256
  return [] if(!@$id_list);
257

258 259
  my @out;
  #construct a constraint like 't1.table1_id = 123'
260
  my @tabs = $self->_tables;
261
  my ($name, $syn) = @{$tabs[0]};
262

263 264 265
  #mysql is faster and we ensure that we do not exceed the max query size by
  #splitting large queries into smaller queries of 200 ids
  my $max_size = 200;
266

267 268 269 270 271 272 273
  while(@$id_list) {
    my @ids;
    if(@$id_list > $max_size) {
      @ids = splice(@$id_list, 0, $max_size);
    } else {
      @ids = splice(@$id_list, 0);
    }
274

275 276 277 278 279 280 281
    
    my $id_str;
    if(@ids > 1)  {
      $id_str = " IN (" . join(',', @ids). ")";
    } else {
      $id_str = " = " . $ids[0];
    }
282

283
    my $constraint = "${syn}.${name}_id $id_str";
284

285 286
    push @out, @{$self->generic_fetch($constraint)};
  }
Arne Stabenau's avatar
Arne Stabenau committed
287

288
  return \@out;
289 290 291 292
}



293
=head2 fetch_all_by_Slice
294 295 296 297 298

  Arg [1]    : Bio::EnsEMBL::Slice $slice
               the slice from which to obtain features
  Arg [2]    : (optional) string $logic_name
               the logic name of the type of features to obtain
299 300 301 302 303 304
  Example    : $fts = $a->fetch_all_by_Slice($slice, 'Swall');
  Description: Returns a listref of features created from the database 
               which are on the Slice defined by $slice. If $logic_name is 
               defined only features with an analysis of type $logic_name 
               will be returned. 
  Returntype : listref of Bio::EnsEMBL::SeqFeatures in Slice coordinates
305 306 307 308 309
  Exceptions : none
  Caller     : Bio::EnsEMBL::Slice

=cut

310
sub fetch_all_by_Slice {
311
  my ($self, $slice, $logic_name) = @_;
312

313
  #fetch by constraint with empty constraint
314
  return $self->fetch_all_by_Slice_constraint($slice, '', $logic_name);
315 316 317
}


318
=head2 fetch_all_by_Slice_and_score
319 320 321

  Arg [1]    : Bio::EnsEMBL::Slice $slice
               the slice from which to obtain features
322
  Arg [2]    : (optional) float $score
323 324 325
               lower bound of the the score of the features retrieved
  Arg [3]    : (optional) string $logic_name
               the logic name of the type of features to obtain
326
  Example    : $fts = $a->fetch_all_by_Slice($slice, 'Swall');
327 328 329 330 331
  Description: Returns a list of features created from the database which are 
               are on the Slice defined by $slice and which have a score 
               greated than $score. If $logic_name is defined, 
               only features with an analysis of type $logic_name will be 
               returned. 
332
  Returntype : listref of Bio::EnsEMBL::SeqFeatures in Slice coordinates
333 334 335 336 337
  Exceptions : none
  Caller     : Bio::EnsEMBL::Slice

=cut

338
sub fetch_all_by_Slice_and_score {
339 340 341
  my ($self, $slice, $score, $logic_name) = @_;
  my $constraint;

342
  if(defined $score) {
343 344 345 346
    #get the synonym of the primary_table
    my @tabs = $self->_tables;
    my $syn = $tabs[0]->[1];
    $constraint = "${syn}.score > $score";
347 348
  }

349 350
  return $self->fetch_all_by_Slice_constraint($slice, $constraint, 
					      $logic_name);
351
}
352 353


354
=head2 fetch_all_by_Slice_constraint
355

356 357 358 359 360
  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)
  Arg [3]    : (optional) string $logic_name
361
               the logic name of the type of features to obtain
362 363
  Example    : $fs = $a->fetch_all_by_Slice_constraint($slc, 'perc_ident > 5');
  Description: Returns a listref of features created from the database which 
364 365 366 367
               are on the Slice defined by $slice and fulfill the SQL 
               constraint defined by $constraint. If logic name is defined, 
               only features with an analysis of type $logic_name will be 
               returned. 
368
  Returntype : listref of Bio::EnsEMBL::SeqFeatures in Slice coordinates
369 370
  Exceptions : thrown if $slice is not defined
  Caller     : Bio::EnsEMBL::Slice
371 372 373

=cut

374
sub fetch_all_by_Slice_constraint {
375 376
  my($self, $orig_slice, $original_constraint, $logic_name) = @_;
  my @result_features;
377

378
  if(!ref($orig_slice) || !$orig_slice->isa("Bio::EnsEMBL::Slice")) {
Graham McVicker's avatar
Graham McVicker committed
379
    throw("Bio::EnsEMBL::Slice argument expected.");
380 381
  }

382 383 384 385 386 387
  $original_constraint ||= '';
  $original_constraint = $self->_logic_name_to_constraint($original_constraint,
                                                          $logic_name);

  #if the logic name was invalid, undef was returned
  return [] if(!defined($original_constraint));
388

389
  #check the cache and return if we have already done this query
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 451 452 453 454 455 456 457 458 459 460 461 462 463 464 465 466 467 468 469 470 471 472 473 474 475 476 477 478 479 480 481 482 483 484 485 486 487 488 489 490 491 492 493 494 495 496 497 498 499 500 501 502 503 504 505 506 507 508 509 510 511 512 513 514 515 516 517 518 519
  my $key = uc(join(':', $orig_slice->name, $original_constraint));

  if(exists($self->{'_slice_feature_cache'}->{$key})) {
    return $self->{'_slice_feature_cache'}->{$key};
  }

  my $slice_adaptor = $orig_slice->adaptor();

  #retrieve normalized 'non-symlinked' slices
  #this allows us to support haplotypes and PARs
  my @projection =
    @{$slice_adaptor->fetch_normalized_slice_projection($orig_slice)};

  if(@projection == 0) {
    throw('Could not retrieve normalized Slices. Database contains ' .
          'incorrect assembly_exception information.');
  }

  #we want to retrieve all features calculated on the FULL original slice 
  #as well as any symlinked slices.  

  #Filter out any partial slices from the normalized projection that are on 
  #the same seq region as the original slice
  my @new_projection = 
    grep { $_->[2]->seq_region_name() ne $orig_slice->seq_region_name() } 
    @projection;

  push( @new_projection, [ 1, $orig_slice->length(), $orig_slice ] );

  #fetch features for the primary slice AND all symlinked slices
  foreach my $segment (@new_projection) {
    my ($offset, $slice);
    ($offset, undef, $slice) = @$segment;

    my $slice_start  = $slice->start();
    my $slice_end    = $slice->end();
    my $slice_strand = $slice->strand();
    my $slice_cs     = $slice->coord_system();
    my $slice_seq_region = $slice->seq_region_name();

    #get the synonym and name of the primary_table
    my @tabs = $self->_tables;
    my ($tab_name, $tab_syn) = @{$tabs[0]};

    #find out what coordinate systems the features are in
    my $csa = $self->db->get_CoordSystemAdaptor();
    my @feat_css = @{$csa->fetch_all_by_feature_table($tab_name)};

    my $asma = $self->db->get_AssemblyMapperAdaptor();
    my @features;

    # fetch the features from each coordinate system they are stored in
  COORD_SYSTEM: foreach my $feat_cs (@feat_css) {
      my $mapper;
      my @coords;
      my @ids;
      
      if($feat_cs->equals($slice_cs)) {
        #no mapping is required if this is the same coord system
        my $constraint = $original_constraint;
        
        # obtain seq_region_id of this slice from db
        my $seq_region_id = 
          $self->db->get_SliceAdaptor->get_seq_region_id($slice);
        $constraint .= " AND " if($constraint);
        $constraint .=
          "${tab_syn}.seq_region_id = $seq_region_id AND " .
          "${tab_syn}.seq_region_start <= $slice_end AND " .
          "${tab_syn}.seq_region_end >= $slice_start";
        my $fs = $self->generic_fetch($constraint,undef,$slice);

        #features may still have to have coordinates made relative to slice 
        #start
        $fs = $self->_remap($fs, $mapper, $slice);

        push @features, @$fs;
      } else {
        $mapper = $asma->fetch_by_CoordSystems($slice_cs, $feat_cs);
        
        # Get a list of coordinates and corresponding internal ids for the
        # regions we are interested in
        @coords = $mapper->map($slice_seq_region, $slice_start, $slice_end,
                               $slice_strand, $slice_cs);

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

        next COORD_SYSTEM if(!@coords);

        @ids = map {$_->id()} @coords;
        @ids = @{$asma->seq_regions_to_ids($feat_cs, \@ids)};

        #if the regions are large and only partially spanned
        #it is faster to to limit the query with start and end constraints
        #however, it is difficult to tell if a region is large and only 
        #partially wanted. The easy approach is just to limit the queries if 
        #there are less than a certain number of regions. As well seperate 
        #queries are needed otherwise the indices will not be useful
        if(@coords > $MAX_SPLIT_QUERY_SEQ_REGIONS) {
          #do one query, and do not limit with start / end constraints
          my $constraint = $original_constraint;
          my $id_str = join(',', @ids);
          $constraint .= " AND " if($constraint);
          $constraint .= "${tab_syn}.seq_region_id IN ($id_str)";
          
          my $fs = 
            $self->generic_fetch($constraint, $mapper, $slice);
          
          $fs = $self->_remap($fs, $mapper, $slice);

          push @features, @$fs;

        } else {
          #do multiple split queries using start / end constraints
          my $len = @coords;
          for(my $i = 0; $i < $len; $i++) {
            my $constraint = $original_constraint;
            $constraint .= " AND " if($constraint);
            $constraint .=
              "${tab_syn}.seq_region_id = "     . $ids[$i] . " AND " .
              "${tab_syn}.seq_region_start <= " . $coords[$i]->end() . " AND ".
              "${tab_syn}.seq_region_end >= "   . $coords[$i]->start();
            my $fs = $self->generic_fetch($constraint,$mapper,$slice);

            $fs = $self->_remap($fs, $mapper, $slice);
            
            push @features, @$fs;
          }
        }
      }
    } #COORD system loop
Web Admin's avatar
Web Admin committed
520
    
521 522 523 524 525 526 527 528 529 530 531 532 533 534 535 536
    #if this was a symlinked slice offset the feature coordinates as needed
    if($slice != $orig_slice) {
      foreach my $f (@features) {
        #function calls are slow!
        if($offset != 1) {
          $f->{'start'} += $offset-1;
          $f->{'end'}   += $offset-1;
        }
        $f->{'slice'} = $orig_slice;
        push @result_features, $f;
      }
    } else {
      push @result_features, @features;
    }
    
  } #slice & symmlinked slice loop
537

538
  $self->{'_slice_feature_cache'}->{$key} = \@result_features;
539

540 541
  return \@result_features;
}
542

543

544 545 546 547 548 549 550 551 552 553 554 555 556 557 558 559 560 561
#
# Helper function containing some common feature storing functionality
#
# Given a Feature this will return a copy (or the same feature if no changes 
# to the feature are needed) of the feature which is relative to the start
# of the seq_region it is on. The seq_region_id of the seq_region it is on
# is also returned.
#
# This method will also ensure that the database knows which coordinate
# systems that this feature is stored in.
#

sub _pre_store {
  my $self    = shift;
  my $feature = shift;

  if(!ref($feature) || !$feature->isa('Bio::EnsEMBL::Feature')) {
    throw('Expected Feature argument.');
562
  }
563

564 565 566
  $self->_check_start_end_strand($feature->start(),$feature->end(),
                                 $feature->strand());

567 568 569 570

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

  my $slice_adaptor = $db->get_SliceAdaptor();
571 572 573 574 575 576 577 578 579 580 581 582 583 584 585 586 587 588 589 590 591 592 593 594 595 596 597 598 599
  my $slice = $feature->slice();

  if(!ref($slice) || !$slice->isa('Bio::EnsEMBL::Slice')) {
    throw('Feature must be attached to Slice to be stored.');
  }

  # make sure that the feature coordinates are relative to
  # the start of the entire seq_region
  if($slice->start != 1 || $slice->strand != 1) {
    #move the feature onto a slice of the entire seq_region
    $slice = $slice_adaptor->fetch_by_region($slice->coord_system->name(),
                                             $slice->seq_region_name(),
                                             undef, #start
                                             undef, #end
                                             undef, #strand
                                             $slice->coord_system->version());

    $feature = $feature->transfer($slice);

    if(!$feature) {
      throw('Could not transfer Feature to slice of ' .
            'entire seq_region prior to storing');
    }
  }

  #
  # Ensure that this type of feature is known to be stored in this coord
  # system.
  #
600
  my $csa = $db->get_CoordSystemAdaptor();
601 602 603 604 605 606 607
  my $cs = $slice->coord_system;

  my ($tab) = $self->_tables();
  my $tabname = $tab->[0];

  $csa->add_feature_table($cs, $tabname);

608 609 610 611 612 613 614 615 616 617 618 619 620 621
  # we have to update the meta coord table in both the dna db and the feature
  # db so that the feature db can be used independently later

  if($db->dnadb() != $db) {
    my $dnadb = $db->dnadb();
    $db->dnadb(undef); # unset the dnadb temporarily

    #get a coord system adaptor from the feature database
    $csa = $db->get_CoordSystemAdaptor();
    $csa->add_feature_table($cs, $tabname);

    $db->dnadb($dnadb); # reinstate the dnadb
  }

622 623 624 625 626 627 628 629 630 631 632 633 634 635 636 637 638 639 640 641 642 643 644 645 646 647 648 649 650 651 652 653 654 655 656 657 658 659 660 661 662 663 664 665 666 667 668 669 670 671 672 673 674
  my $seq_region_id = $slice_adaptor->get_seq_region_id($slice);

  if(!$seq_region_id) {
    throw('Feature is associated with seq_region which is not in this DB.');
  }

  return ($feature, $seq_region_id);
}


#
# helper function used to validate start/end/strand and 
# hstart/hend/hstrand etc.
#
sub _check_start_end_strand {
  my $self = shift;
  my $start = shift;
  my $end   = shift;
  my $strand = shift;

  #
  # Make sure that the start, end, strand are valid
  #
  if(int($start) != $start) {
    throw("Invalid Feature start [$start].  Must be integer.");
  }
  if(int($end) != $end) {
    throw("Invalid Feature end [$end]. Must be integer.");
  }
  if(int($strand) != $strand || $strand < -1 || $strand > 1) {
    throw("Invalid Feature strand [$strand]. Must be -1, 0 or 1.");
  }
  if($end < $start) {
    throw("Invalid Feature start/end [$start/$end]. Start must be less " .
          "than or equal to end.");
  }

  return 1;
}


#
# Given a list of features checks if they are in the correct coord system
# by looking at the first features slice.  If they are not then they are
# converted and placed on the slice.
#
sub _remap {
  my ($self, $features, $mapper, $slice) = @_;

  #check if any remapping is actually needed
  if(@$features && (!$features->[0]->isa('Bio::EnsEMBL::Feature') ||
                    $features->[0]->slice == $slice)) {
    return $features;
675
  }
676

677
  #remapping has not been done, we have to do our own conversion from
678 679 680 681 682 683 684 685 686 687
  #to slice coords

  my @out;

  my $slice_start = $slice->start();
  my $slice_end   = $slice->end();
  my $slice_strand = $slice->strand();
  my $slice_cs    = $slice->coord_system();

  my ($seq_region, $start, $end, $strand);
688

689
  my $slice_seq_region = $slice->seq_region_name();
690

691
  foreach my $f (@$features) {
692
    #since feats were obtained in contig coords, attached seq is a contig
693 694 695 696 697 698
    my $fslice = $f->slice();
    if(!$fslice) {
      throw("Feature does not have attached slice.\n");
    }
    my $fseq_region = $fslice->seq_region_name();
    my $fcs = $fslice->coord_system();
699

700 701
    if(!$slice_cs->equals($fcs)) {
      #slice of feature in different coord system, mapping required
702

703 704 705 706 707 708 709 710 711 712 713 714 715 716
      ($seq_region, $start, $end, $strand) =
        $mapper->fastmap($fseq_region,$f->start(),$f->end(),$f->strand(),$fcs);

      # undefined start means gap
      next if(!defined $start);
    } else {
      $start      = $f->start();
      $end        = $f->end();
      $strand     = $f->strand();
      $seq_region = $f->slice->seq_region_name();
    }

    # maps to region outside desired area
    next if ($start > $slice_end) || ($end < $slice_start) || ($slice_seq_region ne $seq_region);
717

718
    #shift the feature start, end and strand in one call
719
    if($slice_strand == -1) {
720
      $f->move( $slice_end - $end + 1, $slice_end - $start + 1, $strand * -1 );
721 722 723
    } else {
      $f->move( $start - $slice_start + 1, $end - $slice_start + 1, $strand );
    }
724 725 726

    $f->slice($slice);

Web Admin's avatar
Web Admin committed
727
    push @out,$f;
728
  }
729 730 731 732 733 734 735 736 737 738 739 740 741 742 743 744 745 746 747 748 749 750 751 752 753 754 755 756 757 758 759 760 761 762 763 764 765 766 767 768 769 770 771 772 773 774 775 776 777 778 779 780

  return \@out;
}


#
# Given a logic name and an existing constraint this will
# add an analysis table constraint to the feature.  Note that if no
# analysis_id exists in the columns of the primary table then no
# constraint is added at all
#
sub _logic_name_to_constraint {
  my $self = shift;
  my $constraint = shift;
  my $logic_name = shift;

  return $constraint if(!$logic_name);

  #make sure that an analysis_id exists in the primary table
  my ($prim_tab) = $self->_tables();
  my $prim_synonym = $prim_tab->[1];

  my $found_analysis=0;
  foreach my $col ($self->_columns) {
    my ($syn,$col_name) = split(/\./,$col);
    next if($syn ne $prim_synonym);
    if($col_name eq 'analysis_id') {
      $found_analysis = 1;
      last;
    }
  }

  if(!$found_analysis) {
    warning("This feature is not associated with an analysis.\n" .
            "Ignoring logic_name argument = [$logic_name].\n");
    return $constraint;
  }

  my $aa = $self->db->get_AnalysisAdaptor();
  my $an = $aa->fetch_by_logic_name($logic_name);

  if(!$an) {
    warning("No analysis exists with logic_name = [$logic_name].\n" .
            "Returning empty list.\n");
    return undef;
  }

  my $an_id = $an->dbID();

  $constraint .= ' AND' if($constraint);
  $constraint .= " ${prim_synonym}.analysis_id = $an_id";
  return $constraint;
781 782 783
}


784
=head2 store
785

786
  Arg [1]    : list of Bio::EnsEMBL::SeqFeature
787 788 789 790 791 792 793
  Example    : $adaptor->store(@feats);
  Description: ABSTRACT  Subclasses are responsible for implementing this 
               method.  It should take a list of features and store them in 
               the database.
  Returntype : none
  Exceptions : thrown method is not implemented by subclass
  Caller     : general
794 795 796 797 798 799

=cut

sub store{
  my $self = @_;

800
  throw("Abstract method store not defined by implementing subclass\n");
801 802 803
}


804 805 806 807 808 809 810 811 812
=head2 remove

  Arg [1]    : A feature $feature 
  Example    : $feature_adaptor->remove($feature);
  Description: This removes a feature from the database.  The table the
               feature is removed from is defined by the abstract method
               _tablename, and the primary key of the table is assumed
               to be _tablename() . '_id'.  The feature argument must 
               be an object implementing the dbID method, and for the
813
               feature to be removed from the database a dbID value must
814 815
               be returned.
  Returntype : none
816 817
  Exceptions : thrown if $feature arg does not implement dbID(), or if
               $feature->dbID is not a true value
818 819 820 821
  Caller     : general

=cut

822

823 824 825
sub remove {
  my ($self, $feature) = @_;

826 827
  if(!$feature || !ref($feature) || !$feature->isa('Bio::EnsEMBL::Feature')) {
    throw('Feature argument is required');
828 829
  }

830 831
  if(!$feature->is_stored($self->db)) {
    throw("This feature is not stored in this database");
832 833
  }

834 835
  my @tabs = $self->_tables;
  my ($table) = @{$tabs[0]};
836 837 838 839

  my $sth = $self->prepare("DELETE FROM $table WHERE ${table}_id = ?");
  $sth->execute($feature->dbID());

840 841 842 843
  #unset the feature dbID ad adaptor
  $feature->dbID(undef);
  $feature->adaptor(undef);

844 845 846 847
  return;
}


848
=head2 remove_by_Slice
849

850 851 852 853 854 855 856 857
  Arg [1]    : Bio::Ensembl::Slice $slice
  Example    : $feature_adaptor->remove_by_RawContig($slice);
  Description: This removes features from the database which lie on a region
               represented by the passed in slice.  Only features which are
               fully contained by the slice are deleted; features which overlap
               the edge of the slice are not removed.
               The table the features are removed from is defined by
               the abstract method_tablename.
858
  Returntype : none
859
  Exceptions : thrown if no slice is supplied
860 861 862 863
  Caller     : general

=cut

864 865
sub remove_by_Slice {
  my ($self, $slice) = @_;
866

867 868
  if(!$slice || !ref($slice) || !$slice->isa('Bio::EnsEMBL::Slice')) {
    throw("Slice argument is required");
869 870
  }

871 872
  my @tabs = $self->_tables;
  my ($table_name) = @{$tabs[0]};
873

874 875 876
  my $seq_region_id = $self->db->get_SliceAdaptor->get_seq_region_id($slice);
  my $start = $slice->start();
  my $end   = $slice->end();
877

878 879 880 881 882 883 884
  #
  # Delete only features fully on the slice, not overlapping ones
  #
  my $sth = $self->prepare("DELETE FROM $table_name " .
                           "WHERE seq_region_id = ? " .
                           "AND   seq_region_start >= ? " .
                           "AND   seq_region_end <= ?");
885

886 887
  $sth->execute($seq_region_id, $start, $end);
  $sth->finish();
888 889 890 891
}



892 893


894 895 896 897 898 899 900 901 902 903 904 905 906 907 908 909
#_tables
#
#  Args       : none
#  Example    : $tablename = $self->_table_name()
#  Description: ABSTRACT PROTECTED Subclasses are responsible for implementing
#               this method.  It should list of [tablename, alias] pairs.  
#               Additionally the primary table (with the dbID, analysis_id, and
#               score) should be the first table in the list.
#               e.g:
#               ( ['repeat_feature',   'rf'],
#                 ['repeat_consensus', 'rc']);
#               used to obtain features.  
#  Returntype : list of [tablename, alias] pairs
#  Exceptions : thrown if not implemented by subclass
#  Caller     : BaseFeatureAdaptor::generic_fetch
#
910

911
sub _tables {
912 913
  my $self = shift;

914
  throw("abstract method _tables not defined by implementing" .
Simon Potter's avatar
Simon Potter committed
915
               " subclass of BaseFeatureAdaptor");
916 917 918
  return undef;
}

919

920 921 922 923 924 925 926 927 928 929 930
#_columns
#
#  Args       : none
#  Example    : $tablename = $self->_columns()
#  Description: ABSTRACT PROTECTED Subclasses are responsible for implementing
#               this method.  It should return a list of columns to be used
#               for feature creation
#  Returntype : list of strings
#  Exceptions : thrown if not implemented by subclass
#  Caller     : BaseFeatureAdaptor::generic_fetch
#
931 932 933 934

sub _columns {
  my $self = shift;

935
  throw("abstract method _columns not defined by implementing" .
Simon Potter's avatar
Simon Potter committed
936
               " subclass of BaseFeatureAdaptor");
937 938 939
}


940

941 942 943 944 945 946 947 948 949 950 951 952
# _default_where_clause
#
#  Arg [1]    : none
#  Example    : none
#  Description: May be overridden to provide an additional where constraint to 
#               the SQL query which is generated to fetch feature records.
#               This constraint is always appended to the end of the generated
#               where clause
#  Returntype : string
#  Exceptions : none
#  Caller     : generic_fetch
#
953

954 955
sub _default_where_clause {
  my $self = shift;
956

957 958
  return '';
}
959 960 961



962
# _left_join
963

964 965 966 967 968 969 970 971 972
#  Arg [1]    : none
#  Example    : none
#  Description: Can be overridden by a subclass to specify any left joins
#               which should occur. The table name specigfied in the join
#               must still be present in the return values of
#  Returntype : a {'tablename' => 'join condition'} pair
#  Exceptions : none
#  Caller     : general
#
973

974 975
sub _left_join {
  my $self = shift;
976

977 978
  return ();
}
979 980


981 982 983 984 985 986 987 988 989 990 991 992 993 994

#_final_clause

#  Arg [1]    : none
#  Example    : none
#  Description: May be overriden to provide an additional clause to the end
#               of the SQL query used to fetch feature records.  
#               This is useful to add a required ORDER BY clause to the 
#               query for example.
#  Returntype : string
#  Exceptions : none
#  Caller     : generic_fetch

sub _final_clause {
995 996 997 998 999 1000
  my $self = shift;

  return '';
}


1001
#_objs_from_sth
1002

1003 1004 1005 1006 1007 1008 1009 1010 1011
#  Arg [1]    : DBI::row_hashref $hashref containing key-value pairs 
#               for each of the columns specified by the _columns method
#  Example    : my @feats = $self->_obj_from_hashref
#  Description: ABSTRACT PROTECTED The subclass is responsible for implementing
#               this method.  It should take in a DBI row hash reference and
#               return a list of created features in contig coordinates.
#  Returntype : list of Bio::EnsEMBL::*Features in contig coordinates
#  Exceptions : thrown if not implemented by subclass
#  Caller     : BaseFeatureAdaptor::generic_fetch
1012

1013 1014
sub _objs_from_sth {
  my $self = shift;
1015

1016 1017 1018
  throw("abstract method _obj_from_sth not defined by implementing"
             . " subclass of BaseFeatureAdaptor");
}
1019

1020 1021 1022 1023 1024 1025 1026 1027 1028 1029 1030 1031
# deleteObj
#
#  Arg [1]    : none
#  Example    : none
#  Description: Cleans up internal caches and references to other objects so
#               that correct garbage collection may occur.
#  Returntype : none
#  Exceptions : none
#  Caller     : Bio::EnsEMBL::DBConnection::deleteObj


sub deleteObj {
1032 1033
  my $self = shift;

1034 1035
  #flush feature cache
  %{$self->{'_slice_feature_cache'}} = ();
1036 1037 1038
}


Graham McVicker's avatar
Graham McVicker committed
1039
=head1 DEPRECATED METHODS
1040

Graham McVicker's avatar
Graham McVicker committed
1041
=cut
1042

1043 1044 1045 1046

=head2 fetch_all_by_RawContig_constraint

  Description: DEPRECATED use fetch_all_by_RawContig_constraint instead
1047 1048 1049

=cut

1050
sub fetch_all_by_RawContig_constraint {
1051
  my $self = shift;
1052 1053
  deprecate('Use fetch_all_by_Slice_constraint() instead.');
  return $self->fetch_all_by_slice_constraint(@_);
1054 1055
}

1056
=head2 fetch_all_by_RawContig
1057

1058
  Description: DEPRECATED use fetch_all_by_Slice instead
1059 1060 1061

=cut

1062
sub fetch_all_by_RawContig {
1063
  my $self = shift;
1064 1065 1066
  deprecate('Use fetch_all_by_Slice() instead.');
  return $self->fetch_all_by_Slice(@_);
}
1067

1068
=head2 fetch_all_by_RawContig_and_score
1069

1070
  Description: DEPRECATED use fetch_all_by_Slice_and_score instead
1071

1072
=cut
1073

1074 1075 1076 1077 1078 1079 1080 1081 1082
sub fetch_all_by_RawContig_and_score{
  my $self = shift;
  deprecate('Use fetch_all_by_Slice_and_score() instead.');
  return $self->fetch_all_by_Slice_and_score(@_);
}

=head2 remove_by_RawContig

  Description: DEPRECATED use remove_by_Slice instead
1083 1084 1085

=cut

1086
sub remove_by_RawContig {
1087
  my $self = shift;
1088 1089
  deprecate("Use remove_by_Slice instead");
  return $self->remove_by_Slice(@_);
1090 1091
}

1092

1093 1094 1095
1;