DnaAlignFeatureAdaptor.pm 21.2 KB
Newer Older
1 2
=head1 LICENSE

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::DnaAlignFeatureAdaptor - Adaptor for DnaAlignFeatures

=head1 SYNOPSIS

36
  $dafa = $registry->get_adaptor( 'Human', 'Core', 'DnaAlignFeature' );
37

38
  my @features = @{ $dafa->fetch_all_by_Slice($slice) };
39

40
  $dafa->store(@features);
41 42 43

=head1 DESCRIPTION

44 45 46
This is an adaptor responsible for the retrieval and storage of
DnaDnaAlignFeatures from the database. This adaptor inherits most of its
functionality from the BaseAlignFeatureAdaptor and BaseFeatureAdaptor
47
superclasses.
48

49
=head1 METHODS
50 51 52 53 54 55 56

=cut


package Bio::EnsEMBL::DBSQL::DnaAlignFeatureAdaptor;
use vars qw(@ISA);
use strict;
57
use Bio::EnsEMBL::DnaDnaAlignFeature;
58
use Bio::EnsEMBL::DBSQL::BaseAlignFeatureAdaptor;
Laura Clarke's avatar
 
Laura Clarke committed
59
use Bio::EnsEMBL::Utils::Exception qw(throw warning);
60

61
@ISA = qw(Bio::EnsEMBL::DBSQL::BaseAlignFeatureAdaptor);
62 63


64
=head2 _tables
Graham McVicker's avatar
Graham McVicker committed
65 66

  Args       : none
67
  Example    : @tabs = $self->_tables
Graham McVicker's avatar
Graham McVicker committed
68
  Description: PROTECTED implementation of the abstract method inherited from
69 70
               BaseFeatureAdaptor.  Returns list of [tablename, alias] pairs
  Returntype : list of listrefs of strings
Graham McVicker's avatar
Graham McVicker committed
71 72
  Exceptions : none
  Caller     : Bio::EnsEMBL::DBSQL::BaseFeatureAdaptor::generic_fetch
73
  Status     : Stable
Graham McVicker's avatar
Graham McVicker committed
74 75 76

=cut

77
sub _tables {
78
  my $self = shift;
79

80
  return (['dna_align_feature', 'daf'],['external_db','exdb']);
81 82
}

Graham McVicker's avatar
Graham McVicker committed
83

84 85 86 87
sub _left_join{
    return (['external_db',"exdb.external_db_id = daf.external_db_id"]);
}

Graham McVicker's avatar
Graham McVicker committed
88 89 90 91 92 93 94 95 96
=head2 _columns

  Args       : none
  Example    : @columns = $self->_columns
  Description: PROTECTED implementation of abstract superclass method.  
               Returns a list of columns that are needed for object creation.
  Returntype : list of strings
  Exceptions : none
  Caller     : Bio::EnsEMBL::DBSQL::BaseFeatureAdaptor::generic_fetch
97
  Status     : Stable
Graham McVicker's avatar
Graham McVicker committed
98 99 100

=cut

101 102
sub _columns {
  my $self = shift;
103

104
  #warning, implementation of _objs_from_sth method depends on order of list
105 106 107 108 109 110 111 112 113 114 115 116 117
  return qw(daf.dna_align_feature_id
            daf.seq_region_id
            daf.analysis_id
            daf.seq_region_start
            daf.seq_region_end
            daf.seq_region_strand
            daf.hit_start
            daf.hit_end
            daf.hit_name
            daf.hit_strand
            daf.cigar_line
            daf.evalue
            daf.perc_ident
118 119
            daf.score
            daf.external_db_id
Eugene Kulesha's avatar
typo  
Eugene Kulesha committed
120
            daf.hcoverage
121 122 123
	    daf.external_data
	    exdb.db_name
	    exdb.db_display_name);
124
}
125

126

127
=head2 store
128

129
  Arg [1]    : list of Bio::EnsEMBL::DnaAlignFeatures @feats
Graham McVicker's avatar
Graham McVicker committed
130
               the features to store in the database
131
  Example    : $dna_align_feature_adaptor->store(@features);
Graham McVicker's avatar
Graham McVicker committed
132 133
  Description: Stores a list of DnaAlignFeatures in the database
  Returntype : none
134 135 136 137 138 139 140 141 142
  Exceptions : throw if any of the provided features cannot be stored
               which may occur if:
                 * The feature does not have an associate Slice
                 * The feature does not have an associated analysis
                 * The Slice the feature is associated with is on a seq_region
                   unknown to this database
               A warning is given if:
                 * The feature has already been stored in this db
  Caller     : Pipeline
143
  Status     : Stable
144 145 146

=cut

147
sub store {
148
  my ( $self, @feats ) = @_;
149

150
  throw("Must call store with features") if ( scalar(@feats) == 0 );
151 152

  my @tabs = $self->_tables;
153
  my ($tablename) = @{ $tabs[0] };
154

155
  my $db               = $self->db();
156 157 158
  my $analysis_adaptor = $db->get_AnalysisAdaptor();

  my $sth = $self->prepare(
159 160 161 162
    "INSERT INTO $tablename (seq_region_id, seq_region_start,
                             seq_region_end, seq_region_strand,
                             hit_start, hit_end, hit_strand, hit_name,
                             cigar_line, analysis_id, score, evalue,
163 164
                             perc_ident, external_db_id, hcoverage)
     VALUES (?,?,?,?,?,?,?,?,?,?,?,?,?,?,?)"    # 15 arguments
165 166 167 168 169 170
  );

FEATURE:
  foreach my $feat (@feats) {
    if ( !ref $feat || !$feat->isa("Bio::EnsEMBL::DnaDnaAlignFeature") )
    {
171
      throw("feature must be a Bio::EnsEMBL::DnaDnaAlignFeature,"
172 173 174
          . " not a ["
          . ref($feat)
          . "]." );
175 176
    }

177 178 179 180
    if ( $feat->is_stored($db) ) {
      warning( "DnaDnaAlignFeature ["
          . $feat->dbID()
          . "] is already stored in this database." );
181 182 183
      next FEATURE;
    }

184 185
    my $hstart  = $feat->hstart();
    my $hend    = $feat->hend();
186
    my $hstrand = $feat->hstrand();
187
    $self->_check_start_end_strand( $hstart, $hend, $hstrand );
188 189

    my $cigar_string = $feat->cigar_string();
190
    if ( !$cigar_string ) {
191
      $cigar_string = $feat->length() . 'M';
192 193
      warning( "DnaDnaAlignFeature does not define a cigar_string.\n"
          . "Assuming ungapped block with cigar_line=$cigar_string ." );
194
    }
195 196

    my $hseqname = $feat->hseqname();
197
    if ( !$hseqname ) {
198 199 200
      throw("DnaDnaAlignFeature must define an hseqname.");
    }

201 202 203
    if ( !defined( $feat->analysis ) ) {
      throw(
        "An analysis must be attached to the features to be stored.");
204 205
    }

206
    #store the analysis if it has not been stored yet
207 208
    if ( !$feat->analysis->is_stored($db) ) {
      $analysis_adaptor->store( $feat->analysis() );
209
    }
210

211 212
    my $original = $feat;
    my $seq_region_id;
213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229
    ( $feat, $seq_region_id ) = $self->_pre_store($feat);

    $sth->bind_param( 1,  $seq_region_id,        SQL_INTEGER );
    $sth->bind_param( 2,  $feat->start,          SQL_INTEGER );
    $sth->bind_param( 3,  $feat->end,            SQL_INTEGER );
    $sth->bind_param( 4,  $feat->strand,         SQL_TINYINT );
    $sth->bind_param( 5,  $hstart,               SQL_INTEGER );
    $sth->bind_param( 6,  $hend,                 SQL_INTEGER );
    $sth->bind_param( 7,  $hstrand,              SQL_TINYINT );
    $sth->bind_param( 8,  $hseqname,             SQL_VARCHAR );
    $sth->bind_param( 9,  $cigar_string,         SQL_LONGVARCHAR );
    $sth->bind_param( 10, $feat->analysis->dbID, SQL_INTEGER );
    $sth->bind_param( 11, $feat->score,          SQL_DOUBLE );
    $sth->bind_param( 12, $feat->p_value,        SQL_DOUBLE );
    $sth->bind_param( 13, $feat->percent_id,     SQL_FLOAT );
    $sth->bind_param( 14, $feat->external_db_id, SQL_INTEGER );
    $sth->bind_param( 15, $feat->hcoverage,      SQL_DOUBLE );
230 231

    $sth->execute();
232 233

    $original->dbID( $sth->{'mysql_insertid'} );
234
    $original->adaptor($self);
235
  } ## end foreach my $feat (@feats)
236 237

  $sth->finish();
238
} ## end sub store
239 240 241 242 243 244 245 246 247 248 249 250 251 252


sub save {
  my ($self, $features) = @_;

  my @feats = @$features;
  throw("Must call store with features") if( scalar(@feats) == 0 );

  my @tabs = $self->_tables;
  my ($tablename) = @{$tabs[0]};

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

253
  my $sql = qq{INSERT INTO $tablename (seq_region_id, seq_region_start, seq_region_end, seq_region_strand, hit_start, hit_end, hit_strand, hit_name, cigar_line, analysis_id, score, evalue, perc_ident, external_db_id, hcoverage, external_data) VALUES (?, ?, ?, ?, ?, ?, ?, ?, ?, ?, ?, ?, ?, ?, ?, ?)};
254

255 256
  my %analyses = ();

257 258 259 260 261 262 263 264 265 266 267 268 269 270
  my $sth = $self->prepare($sql);
     
 FEATURE: foreach my $feat ( @feats ) {
    if( !ref $feat || !$feat->isa("Bio::EnsEMBL::DnaDnaAlignFeature") ) {
      throw("feature must be a Bio::EnsEMBL::DnaDnaAlignFeature,"
            . " not a [".ref($feat)."].");
    }

    if($feat->is_stored($db)) {
      warning("DnaDnaAlignFeature [".$feat->dbID."] is already stored" .
              " in this database.");
      next FEATURE;
    }

271 272 273 274 275 276 277 278
    my $hstart  = $feat->hstart || 0; # defined $feat->hstart  ? $feat->hstart : $feat->start ;
    my $hend    = $feat->hend   || 0; # defined $feat->hend    ? $feat->hend : $feat->end;
    my $hstrand = $feat->hstrand|| 0; # defined $feat->hstrand ? $feat->hstrand : $feat->strand;
    if( $hstart && $hend ) {
      if($hend < $hstart) {
        throw("Invalid Feature start/end [$hstart/$hend]. Start must be less than or equal to end.");
      }
    }
279 280 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295 296 297 298 299
    my $cigar_string = $feat->cigar_string();
    if(!$cigar_string) {
      $cigar_string = $feat->length() . 'M';
      warning("DnaDnaAlignFeature does not define a cigar_string.\n" .
              "Assuming ungapped block with cigar_line=$cigar_string .");
    }

    my $hseqname = $feat->hseqname();
    if(!$hseqname) {
      throw("DnaDnaAlignFeature must define an hseqname.");
    }

    if(!defined($feat->analysis)) {
      throw("An analysis must be attached to the features to be stored.");
    }

    #store the analysis if it has not been stored yet
    if(!$feat->analysis->is_stored($db)) {
      $analysis_adaptor->store($feat->analysis());
    }

300 301
    $analyses{ $feat->analysis->dbID }++;

302 303 304 305
    my $original = $feat;
    my $seq_region_id;
    ($feat, $seq_region_id) = $self->_pre_store_userdata($feat);

306 307
    my $extra_data = $feat->extra_data ? $self->dump_data($feat->extra_data) : '';

308 309 310 311 312 313 314 315 316 317
    $sth->bind_param(1,$seq_region_id,SQL_INTEGER);
    $sth->bind_param(2,$feat->start,SQL_INTEGER);
    $sth->bind_param(3,$feat->end,SQL_INTEGER);
    $sth->bind_param(4,$feat->strand,SQL_TINYINT);
    $sth->bind_param(5,$hstart,SQL_INTEGER);
    $sth->bind_param(6,$hend,SQL_INTEGER);
    $sth->bind_param(7,$hstrand,SQL_TINYINT);
    $sth->bind_param(8,$hseqname,SQL_VARCHAR);
    $sth->bind_param(9,$cigar_string,SQL_LONGVARCHAR);
    $sth->bind_param(10,$feat->analysis->dbID,SQL_INTEGER);
318
    $sth->bind_param(11,$feat->score,SQL_DOUBLE);
319
#    $sth->bind_param(11,$feat->score); # if the above statement does not work it means you need to upgrade DBD::mysql, meantime you can replace it with this line
320 321
    $sth->bind_param(12,$feat->p_value,SQL_DOUBLE);
    $sth->bind_param(13,$feat->percent_id,SQL_FLOAT);
322 323
    $sth->bind_param(14,$feat->external_db_id,SQL_INTEGER);
    $sth->bind_param(15,$feat->hcoverage,SQL_DOUBLE);
324
    $sth->bind_param(16,$extra_data,SQL_LONGVARCHAR);
325

326 327

    $sth->execute();
328 329
    $original->dbID($sth->{'mysql_insertid'});
    $original->adaptor($self);
330
  }
331 332

  $sth->finish();
333 334 335

## js5 hack to update meta_coord table... 
  if( keys %analyses ) {
336

337 338
    my $sth = $self->prepare( 'select sr.coord_system_id, max(daf.seq_region_end-daf.seq_region_start) from seq_region as sr, dna_align_feature as daf where daf.seq_region_id=sr.seq_region_id and analysis_id in ('.join(',',keys %analyses).') group by coord_system_id' );
    $sth->execute;
339

340 341 342 343
    foreach( @{ $sth->fetchall_arrayref } ) {
      my $sth2 = $self->prepare( qq(insert ignore into meta_coord values("dna_align_feature",$_->[0],$_->[1])) );
      $sth2->execute;
      $sth2->finish;
344 345

      $sth2 = $self->prepare( qq(update meta_coord set max_length = $_->[1] where coord_system_id = $_->[0] and table_name="dna_align_feature" and max_length < $_->[1]) );
346 347 348
      $sth2->execute;
      $sth2->finish;
    }
349

350 351 352
    $sth->finish;
  }

353 354 355
}


356
=head2 _objs_from_sth
Graham McVicker's avatar
Graham McVicker committed
357

358 359 360 361
  Arg [1]    : DBI statement handle $sth
               an exectuted DBI statement handle generated by selecting 
               the columns specified by _columns() from the table specified 
               by _table()
Graham McVicker's avatar
Graham McVicker committed
362 363 364
  Example    : @dna_dna_align_feats = $self->_obj_from_hashref
  Description: PROTECTED implementation of superclass abstract method. 
               Creates DnaDnaAlignFeature objects from a DBI hashref
365
  Returntype : listref of Bio::EnsEMBL::DnaDnaAlignFeatures
Graham McVicker's avatar
Graham McVicker committed
366 367
  Exceptions : none
  Caller     : Bio::EnsEMBL::BaseFeatureAdaptor::generic_fetch
368
  Status     : Stable
Graham McVicker's avatar
Graham McVicker committed
369 370

=cut
371

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

375 376 377 378 379
  #
  # 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 is used.
  #
380

381 382 383 384 385
  # In case of userdata we need the features on the dest_slice.  In case
  # of get_all_supporting_features dest_slice is not provided.
  my $sa = (   $dest_slice
             ? $dest_slice->adaptor()
             : $self->db()->get_SliceAdaptor() );
386
  my $aa = $self->db->get_AnalysisAdaptor();
387

388 389 390 391 392 393
  my @features;
  my %analysis_hash;
  my %slice_hash;
  my %sr_name_hash;
  my %sr_cs_hash;

394 395 396 397 398 399 400 401
  my ( $dna_align_feature_id, $seq_region_id,
       $analysis_id,          $seq_region_start,
       $seq_region_end,       $seq_region_strand,
       $hit_start,            $hit_end,
       $hit_name,             $hit_strand,
       $cigar_line,           $evalue,
       $perc_ident,           $score,
       $external_db_id,       $hcoverage,
402
       $extra_data,
403 404 405 406 407 408 409 410 411 412
       $external_db_name,     $external_display_db_name );

  $sth->bind_columns( \( $dna_align_feature_id, $seq_region_id,
                         $analysis_id,          $seq_region_start,
                         $seq_region_end,       $seq_region_strand,
                         $hit_start,            $hit_end,
                         $hit_name,             $hit_strand,
                         $cigar_line,           $evalue,
                         $perc_ident,           $score,
                         $external_db_id,       $hcoverage,
413
                         $extra_data,
414 415
                         $external_db_name, $external_display_db_name )
  );
416 417 418 419 420 421 422

  my $asm_cs;
  my $cmp_cs;
  my $asm_cs_vers;
  my $asm_cs_name;
  my $cmp_cs_vers;
  my $cmp_cs_name;
423 424 425 426

  if ( defined($mapper) ) {
    $asm_cs      = $mapper->assembled_CoordSystem();
    $cmp_cs      = $mapper->component_CoordSystem();
427 428 429 430 431
    $asm_cs_name = $asm_cs->name();
    $asm_cs_vers = $asm_cs->version();
    $cmp_cs_name = $cmp_cs->name();
    $cmp_cs_vers = $cmp_cs->version();
  }
432

433 434 435 436
  my $dest_slice_start;
  my $dest_slice_end;
  my $dest_slice_strand;
  my $dest_slice_length;
437
  my $dest_slice_sr_name;
438
  my $dest_slice_seq_region_id;
439

440 441 442 443 444 445
  if ( defined($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();
    $dest_slice_sr_name       = $dest_slice->seq_region_name();
446
    $dest_slice_seq_region_id = $dest_slice->get_seq_region_id();
447
  }
448

449 450 451
FEATURE:
  while ( $sth->fetch() ) {
    # Get the analysis object.
452 453
    my $analysis = $analysis_hash{$analysis_id} ||=
      $aa->fetch_by_dbID($analysis_id);
454

455 456
    # Get the slice object.
    my $slice = $slice_hash{ "ID:" . $seq_region_id };
457

458
    if ( !defined($slice) ) {
459
      $slice = $sa->fetch_by_seq_region_id($seq_region_id);
460 461
      if ( defined($slice) ) {
        $slice_hash{ "ID:" . $seq_region_id } = $slice;
Anne Parker's avatar
 
Anne Parker committed
462
        $sr_name_hash{$seq_region_id} = $slice->seq_region_name();
463
        $sr_cs_hash{$seq_region_id}   = $slice->coord_system();
Anne Parker's avatar
 
Anne Parker committed
464
      }
465
    }
466

467 468
    my $sr_name = $sr_name_hash{$seq_region_id};
    my $sr_cs   = $sr_cs_hash{$seq_region_id};
469 470 471 472

    # Remap the feature coordinates to another coord system
    # if a mapper was provided.
    if ( defined($mapper) ) {
473 474 475 476 477 478 479 480 481 482 483 484 485 486

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

	} else {

	    ( $seq_region_id,  $seq_region_start,
	      $seq_region_end, $seq_region_strand )
		=
		$mapper->fastmap( $sr_name, $seq_region_start, $seq_region_end,
487
                          $seq_region_strand, $sr_cs );
488
	}
489 490 491 492 493 494 495 496 497

      # Skip features that map to gaps or coord system boundaries.
      if ( !defined($seq_region_id) ) { next FEATURE }

      # Get a slice in the coord system we just mapped to.
      if ( $asm_cs == $sr_cs
           || ( $cmp_cs != $sr_cs && $asm_cs->equals($sr_cs) ) )
      {
        $slice = $slice_hash{ "ID:" . $seq_region_id } ||=
498
          $sa->fetch_by_seq_region_id($seq_region_id);
499
      } else {
500
        $slice = $slice_hash{ "ID:" . $seq_region_id } ||=
501
          $sa->fetch_by_seq_region_id($seq_region_id);
502 503
      }
    }
504

505
    # If a destination slice was provided, convert the coords.  If the
506
    # dest_slice starts at 1 and is forward strand, nothing needs doing.
507
    if ( defined($dest_slice) ) {
508 509 510 511 512 513 514 515 516 517 518 519 520 521 522 523 524 525 526 527 528 529 530 531 532 533 534 535 536 537 538 539 540 541 542 543 544 545 546 547 548 549 550 551 552 553 554 555 556 557 558 559
      my $seq_region_len = $dest_slice->seq_region_length();

      if ($dest_slice_strand == 1) { # Positive strand
		
	$seq_region_start = $seq_region_start - $dest_slice_start + 1;
	$seq_region_end   = $seq_region_end - $dest_slice_start + 1;

	if ($dest_slice->is_circular()) {
	  # Handle cicular chromosomes.

	  if ($seq_region_start > $seq_region_end) {
	    # Looking at a feature overlapping the chromsome origin.

	    if ($seq_region_end > $dest_slice_start) {

	      # Looking at the region in the beginning of the
	      # chromosome.
	      $seq_region_start -= $seq_region_len;
	    }

	    if ($seq_region_end < 0) {
	      $seq_region_end += $seq_region_len;
	    }

	  } else {

	    if (   $dest_slice_start > $dest_slice_end
		   && $seq_region_end < 0) {
	      # Looking at the region overlapping the chromosome
	      # origin and a feature which is at the beginning of the
	      # chromosome.
	      $seq_region_start += $seq_region_len;
	      $seq_region_end   += $seq_region_len;
	    }
	  }

	}		       ## end if ($dest_slice->is_circular...)

      } else {			# Negative strand

	my $start = $dest_slice_end - $seq_region_end + 1;
	my $end = $dest_slice_end - $seq_region_start + 1;

	if ($dest_slice->is_circular()) {

	  if ($dest_slice_start > $dest_slice_end) { 
	    # slice spans origin or replication

	    if ($seq_region_start >= $dest_slice_start) {
	      $end += $seq_region_len;
	      $start += $seq_region_len 
		if $seq_region_end > $dest_slice_start;
560

561 562 563 564 565 566 567 568 569 570 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 600 601 602 603 604 605 606 607 608 609
	    } elsif ($seq_region_start <= $dest_slice_end) {
	      # do nothing
	    } elsif ($seq_region_end >= $dest_slice_start) {
	      $start += $seq_region_len;
	      $end += $seq_region_len;

	    } elsif ($seq_region_end <= $dest_slice_end) {

	      $end += $seq_region_len
		if $end < 0;

	    } elsif ($seq_region_start > $seq_region_end) {
		  
	      $end += $seq_region_len;

	    } else {
		  
	    }
      
	  } else {

	    if ($seq_region_start <= $dest_slice_end and $seq_region_end >= $dest_slice_start) {
	      # do nothing
	    } elsif ($seq_region_start > $seq_region_end) {
	      if ($seq_region_start <= $dest_slice_end) {
	  
		$start -= $seq_region_len;

	      } elsif ($seq_region_end >= $dest_slice_start) {
		$end += $seq_region_len;

	      } else {
		    
	      }
	    }
	  }

	}

	$seq_region_start = $start;
	$seq_region_end = $end;
	$seq_region_strand *= -1;

      }	## end else [ if ($dest_slice_strand...)]

      # Throw away features off the end of the requested slice.
      if (    $seq_region_end < 1
	      || $seq_region_start > $dest_slice_length
	      || ( $dest_slice_seq_region_id ne $seq_region_id ) )
610
        {
611 612
          next FEATURE;
        }
613

614
      $slice = $dest_slice;
615
    }
616

617 618 619
    # Inlining the following in the hash causes major issues with 5.16 and messes up the hash 
    my $evalled_extra_data = $extra_data ? $self->get_dumped_data($extra_data) : '';

620 621 622
    # Finally, create the new DnaAlignFeature.
    push( @features,
          $self->_create_feature_fast(
623 624 625 626 627 628 629 630 631 632 633 634 635 636 637 638 639 640
             'Bio::EnsEMBL::DnaDnaAlignFeature', {
               'slice'          => $slice,
               'start'          => $seq_region_start,
               'end'            => $seq_region_end,
               'strand'         => $seq_region_strand,
               'hseqname'       => $hit_name,
               'hstart'         => $hit_start,
               'hend'           => $hit_end,
               'hstrand'        => $hit_strand,
               'score'          => $score,
               'p_value'        => $evalue,
               'percent_id'     => $perc_ident,
               'cigar_string'   => $cigar_line,
               'analysis'       => $analysis,
               'adaptor'        => $self,
               'dbID'           => $dna_align_feature_id,
               'external_db_id' => $external_db_id,
               'hcoverage'      => $hcoverage,
641
               'extra_data'     => $evalled_extra_data,
642
               'dbname'                    => $external_db_name,
643
               'db_display_name'           => $external_display_db_name
644 645 646
             } ) );

  } ## end while ( $sth->fetch() )
647

648
  return \@features;
649
} ## end sub _objs_from_sth
650

Glenn Proctor's avatar
Glenn Proctor committed
651 652 653 654
=head2 list_dbIDs

  Arg [1]    : none
  Example    : @feature_ids = @{$dna_align_feature_adaptor->list_dbIDs()};
655 656
  Description: Gets an array of internal ids for all dna align features in 
               the current db
657
  Arg[1]     : <optional> int. not 0 for the ids to be sorted by the seq_region.
Glenn Proctor's avatar
Glenn Proctor committed
658 659 660
  Returntype : list of ints
  Exceptions : none
  Caller     : ?
661
  Status     : Stable
Glenn Proctor's avatar
Glenn Proctor committed
662 663 664 665

=cut

sub list_dbIDs {
666
   my ($self, $ordered) = @_;
Glenn Proctor's avatar
Glenn Proctor committed
667

668
   return $self->_list_dbIDs("dna_align_feature",undef, $ordered);
Glenn Proctor's avatar
Glenn Proctor committed
669 670
}

671 672 673
1;