ProteinAlignFeatureAdaptor.pm 14.5 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

=head1 NAME

32 33
Bio::EnsEMBL::DBSQL::ProteinAlignFeatureAdaptor - 
Adaptor for ProteinAlignFeatures
34 35 36

=head1 SYNOPSIS

37 38
  $pafa =
    $registry->get_adaptor( 'Human', 'Core', 'ProteinAlignFeature' );
39

40
  my @features = @{ $pafa->fetch_all_by_Slice($slice) };
41

42
  $pafa->store(@features);
43

44
=head1 METHODS
45 46 47 48 49 50 51 52

=cut


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

53
use Bio::EnsEMBL::DBSQL::BaseAlignFeatureAdaptor;
54
use Bio::EnsEMBL::DnaPepAlignFeature;
Laura Clarke's avatar
 
Laura Clarke committed
55
use Bio::EnsEMBL::Utils::Exception qw(throw warning);
56
@ISA = qw(Bio::EnsEMBL::DBSQL::BaseAlignFeatureAdaptor);
57 58


59 60
=head2 store

61 62 63
  Arg [1]    : list of Bio::EnsEMBL::DnaPepAlignFeature @feats
  Example    : $protein_align_feature_adaptor->store(@feats);
  Description: stores a list of ProteinAlignFeatures in the database
64
  Returntype : none
65 66 67 68 69 70 71 72 73
  Exceptions : throw if any of the provided features cannot be stored
               which may occur if:
                 * The feature does not have an associated 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
74
  Status     : Stable
75 76 77

=cut

78

79
sub store{
80 81 82 83 84 85 86 87 88 89 90 91 92 93 94
 my ($self, @feats) = @_;

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

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

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

  my $sth = $self->prepare(
     "INSERT INTO $tablename (seq_region_id, seq_region_start, seq_region_end,
                             seq_region_strand, hit_start, hit_end,
                             hit_name, cigar_line,
95 96
                             analysis_id, score, evalue, perc_ident, external_db_id, hcoverage)
     VALUES (?,?,?,?,?,?,?,?,?,?, ?, ?, ?, ?)");
97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120

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

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

   #sanity check the hstart and hend
   my $hstart  = $feat->hstart();
   my $hend    = $feat->hend();
   $self->_check_start_end_strand($hstart,$hend,1);

   my $cigar_string = $feat->cigar_string();
   if(!$cigar_string) {
     $cigar_string = $feat->length() . 'M';
     warning("DnaPepAlignFeature does not define a cigar_string.\n" .
             "Assuming ungapped block with cigar_string=$cigar_string\n");
   }
121

122 123 124
   my $hseqname = $feat->hseqname();
   if(!$hseqname) {
     throw("DnaPepAlignFeature must define an hseqname.");
125 126
   }

127 128
   if(!defined($feat->analysis)) {
     throw("An analysis must be attached to the features to be stored.");
129
   }
130 131 132 133 134 135 136

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

   my $slice = $feat->slice();
Ian Longden's avatar
Ian Longden committed
137
   if(!defined($slice) || !($slice->isa("Bio::EnsEMBL::Slice") or $slice->isa('Bio::EnsEMBL::LRGSlice')) ) {
138 139 140 141 142 143 144
     throw("A slice must be attached to the features to be stored.");
   }

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

145 146 147 148 149 150 151 152 153 154 155 156
   $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,$feat->hstart,SQL_INTEGER);
   $sth->bind_param(6,$feat->hend,SQL_INTEGER);
   $sth->bind_param(7,$feat->hseqname,SQL_VARCHAR);
   $sth->bind_param(8,$feat->cigar_string,SQL_LONGVARCHAR);
   $sth->bind_param(9,$feat->analysis->dbID,SQL_INTEGER);
   $sth->bind_param(10,$feat->score,SQL_DOUBLE);
   $sth->bind_param(11,$feat->p_value,SQL_DOUBLE);
   $sth->bind_param(12,$feat->percent_id,SQL_REAL);
157 158
   $sth->bind_param(13,$feat->external_db_id,SQL_INTEGER);
   $sth->bind_param(14,$feat->hcoverage,SQL_DOUBLE);
159 160

   $sth->execute();
161 162
   $original->dbID($sth->{'mysql_insertid'});
   $original->adaptor($self);
163
 }
164

165 166 167
 $sth->finish();
}

168

169
=head2 _objs_from_sth
170

171 172 173 174 175 176 177 178 179 180
  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()
  Example    : @dna_dna_align_feats = $self->_obj_from_hashref
  Description: PROTECTED implementation of superclass abstract method. 
               Creates DnaDnaAlignFeature objects from a DBI hashref
  Returntype : listref of Bio::EnsEMBL::ProteinAlignFeatures
  Exceptions : none
  Caller     : Bio::EnsEMBL::BaseFeatureAdaptor::generic_fetch
181
  Status     : Stable
182

183
=cut
184

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

188 189 190 191 192
  #
  # 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.
  #
193

194 195
  my $sa = $self->db()->get_SliceAdaptor();
  my $aa = $self->db->get_AnalysisAdaptor();
196

197 198 199 200 201 202 203 204
  my @features;
  my %analysis_hash;
  my %slice_hash;
  my %sr_name_hash;
  my %sr_cs_hash;

  my ($protein_align_feature_id, $seq_region_id, $seq_region_start,
      $seq_region_end, $analysis_id, $seq_region_strand, $hit_start,
205
      $hit_end, $hit_name, $cigar_line, $evalue, $perc_ident, $score,
206
      $external_db_id, $hcoverage, $external_db_name, $external_display_db_name );
207 208 209 210

  $sth->bind_columns(\$protein_align_feature_id, \$seq_region_id,
           \$seq_region_start,\$seq_region_end, \$analysis_id,
           \$seq_region_strand, \$hit_start,\$hit_end, \$hit_name,
211
           \$cigar_line, \$evalue, \$perc_ident, \$score,
212
           \$external_db_id, \$hcoverage, \$external_db_name, \$external_display_db_name );
213 214 215 216 217 218 219 220 221 222 223 224 225 226 227

  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();
  }
228

229 230 231 232
  my $dest_slice_start;
  my $dest_slice_end;
  my $dest_slice_strand;
  my $dest_slice_length;
233
  my $dest_slice_sr_name;
234
  my $dest_slice_sr_id;
235 236 237 238 239
  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();
240
    $dest_slice_sr_name = $dest_slice->seq_region_name();
241
    $dest_slice_sr_id  = $dest_slice->get_seq_region_id();
242
  }
243

244 245 246 247
  FEATURE: while($sth->fetch()) {
    #get the analysis object
    my $analysis = $analysis_hash{$analysis_id} ||=
      $aa->fetch_by_dbID($analysis_id);
248

249 250
    #get the slice object
    my $slice = $slice_hash{"ID:".$seq_region_id};
251

252 253 254 255 256 257
    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();
    }
258

259 260
    my $sr_name = $sr_name_hash{$seq_region_id};
    my $sr_cs   = $sr_cs_hash{$seq_region_id};
261 262 263 264 265 266
    #
    # remap the feature coordinates to another coord system
    # if a mapper was provided
    #
    if($mapper) {

267 268 269 270 271 272 273 274 275 276 277 278 279 280 281
	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,
                          $seq_region_strand, $sr_cs );
	}
282 283

      #skip features that map to gaps or coord system boundaries
284 285 286 287 288 289 290 291 292 293 294
      next FEATURE if(!defined($seq_region_id));

 #     #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} ||=
          $sa->fetch_by_seq_region_id($seq_region_id);
#      } else {
#        $slice = $slice_hash{"ID:".$seq_region_id} ||=
#          $sa->fetch_by_seq_region_id($asm_cs_name, $sr_name, undef, undef, undef,
#                               $asm_cs_vers);
#      }
295
    }
296

297 298 299 300
    #
    # If a destination slice was provided convert the coords
    #
    if($dest_slice) {
301 302 303 304 305 306 307 308 309 310 311 312 313 314 315 316 317 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 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
      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;

	    } 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...)]
398 399 400
       
      #throw away features off the end of the requested slice
      if($seq_region_end < 1 || $seq_region_start > $dest_slice_length ||
401
	( $dest_slice_sr_id ne $seq_region_id )) {
402
	next FEATURE;
403
      }
404

405
      $slice = $dest_slice;
406
    }
407

408 409 410 411 412 413 414 415 416 417 418 419 420 421 422 423 424 425 426 427 428 429
    # Finally, create the new ProteinAlignFeature.
    push(
      @features,
      $self->_create_feature_fast(
        'Bio::EnsEMBL::DnaPepAlignFeature', {
          'slice'        => $slice,
          'start'        => $seq_region_start,
          'end'          => $seq_region_end,
          'strand'       => $seq_region_strand,
          'hseqname'     => $hit_name,
          'hstart'       => $hit_start,
          'hend'         => $hit_end,
          'hstrand'      => 1,                  # dna_pep_align features
                                                # are always hstrand 1
          'score'        => $score,
          'p_value'      => $evalue,
          'percent_id'   => $perc_ident,
          'cigar_string' => $cigar_line,
          'analysis'     => $analysis,
          'adaptor'      => $self,
          'dbID'           => $protein_align_feature_id,
          'external_db_id' => $external_db_id,
430 431 432
          'hcoverage'      => $hcoverage,
	  'dbname'         => $external_db_name,
	  'db_display_name' => $external_display_db_name
433 434
        } ) );

435
  }
436

437
  return \@features;
438
}
439 440


441

442
sub _tables {
443 444
  my $self = shift;

445
  return (['protein_align_feature', 'paf'], ['external_db','exdb']);
446 447
}

448

449 450
sub _columns {
  my $self = shift;
451 452

  #warning _objs_from_hashref method depends on ordering of this list 
453 454 455 456 457 458 459 460 461 462 463 464
  return qw( paf.protein_align_feature_id
             paf.seq_region_id
             paf.seq_region_start
             paf.seq_region_end
             paf.analysis_id
             paf.seq_region_strand
             paf.hit_start
             paf.hit_end
             paf.hit_name
             paf.cigar_line
             paf.evalue
             paf.perc_ident
465 466
             paf.score
             paf.external_db_id
467 468 469 470 471 472 473
             paf.hcoverage
	     exdb.db_name
	     exdb.db_display_name);
}

sub _left_join{
    return (['external_db',"exdb.external_db_id = paf.external_db_id"]);
474 475
}

Glenn Proctor's avatar
Glenn Proctor committed
476 477 478 479
=head2 list_dbIDs

  Arg [1]    : none
  Example    : @feature_ids = @{$protein_align_feature_adaptor->list_dbIDs()};
480 481
  Description: Gets an array of internal ids for all protein align 
               features in the current db
482
  Arg[1]     : <optional> int. not 0 for the ids to be sorted by the seq_region.
483
  Returntype : listref of ints
Glenn Proctor's avatar
Glenn Proctor committed
484 485
  Exceptions : none
  Caller     : ?
486
  Status     : Stable
Glenn Proctor's avatar
Glenn Proctor committed
487 488 489 490

=cut

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

493
   return $self->_list_dbIDs("protein_align_feature", undef, $ordered);
Glenn Proctor's avatar
Glenn Proctor committed
494
}
495 496 497 498 499 500






501
1;