PredictionTranscriptAdaptor.pm 17.2 KB
Newer Older
1
2
=head1 LICENSE

Andy Yates's avatar
Andy Yates committed
3
  Copyright (c) 1999-2012 The European Bioinformatics Institute and
4
5
6
7
8
9
10
11
12
13
  Genome Research Limited.  All rights reserved.

  This software is distributed under a modified Apache license.
  For license details, please see

    http://www.ensembl.org/info/about/code_licence.html

=head1 CONTACT

  Please email comments or questions to the public Ensembl
14
  developers list at <dev@ensembl.org>.
15
16
17
18
19

  Questions may also be sent to the Ensembl help desk at
  <helpdesk@ensembl.org>.

=cut
20
21
22

=head1 NAME

23
24
Bio::EnsEMBL::DBSQL::PredictionTranscriptAdaptor -
Performs database interaction related to PredictionTranscripts
25
26
27

=head1 SYNOPSIS

28
29
  # get a prediction transcript adaptor from the database
  $pta = $database_adaptor->get_PredictionTranscriptAdaptor();
30

31
32
  # get a slice on a region of chromosome 1
  $sa = $database_adaptor->get_SliceAdaptor();
Jan-hinnerk Vogel's avatar
 
Jan-hinnerk Vogel committed
33

34
  $slice = $sa->fetch_by_region( 'chromosome', 'x', 100000, 200000 );
35

36
37
  # get all the prediction transcripts from the slice region
  $prediction_transcripts = @{ $pta->fetch_all_by_Slice($slice) };
38

39
=head1 METHODS
40
41
42
43
44
45
46
47

=cut

package Bio::EnsEMBL::DBSQL::PredictionTranscriptAdaptor;

use vars qw( @ISA );
use strict;

48
use Bio::EnsEMBL::DBSQL::BaseFeatureAdaptor;
49
use Bio::EnsEMBL::DBSQL::AnalysisAdaptor;
50
use Bio::EnsEMBL::PredictionTranscript;
51
use Bio::EnsEMBL::Utils::Exception qw(deprecate throw warning);
52

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


56
57
58
59
60
61
62
63
64
# _tables
#
#  Arg [1]    : none
#  Example    : none
#  Description: Implements abstract superclass method to define the table used
#               to retrieve prediction transcripts from the database
#  Returntype : string
#  Exceptions : none
#  Caller     : generic_fetch
65

66
sub _tables {
67
  my $self = shift;
68

69
  return ['prediction_transcript', 'pt'];
70
}
71

72

73
# _columns
74

75
76
77
78
79
80
81
82
83
#  Arg [1]    : none
#  Example    : none
#  Description: Implements abstract superclass method to define the columns
#               retrieved in database queries used to create prediction 
#               transcripts.
#  Returntype : list of strings
#  Exceptions : none
#  Caller     : generic_fetch
#
84

85
sub _columns {
86
87
  my $self = shift;

88
89
90
91
92
  return qw( pt.prediction_transcript_id
             pt.seq_region_id
             pt.seq_region_start
             pt.seq_region_end
             pt.seq_region_strand
93
94
             pt.analysis_id
             pt.display_label);
95
96
}

97

98
99
=head2 fetch_by_stable_id

100
101
102
103
104
105
106
107
108
109
110
111
  Arg [1]    : string $stable_id
               The stable id of the transcript to retrieve
  Example    : $trans = $trans_adptr->fetch_by_stable_id('GENSCAN00000001234');
  Description: Retrieves a prediction transcript via its display_label.
               This method is called fetch_by_stable_id for polymorphism with
               the TranscriptAdaptor.  Prediction transcript display_labels are
               not necessarily stable in that the same identifier may be reused
               for a completely different prediction transcript in a subsequent
               database release.
  Returntype : Bio::EnsEMBL::PredictionTranscript
  Caller     : general
  Status     : Stable
112
113
114
115
116

=cut

sub fetch_by_stable_id {
  my $self = shift;
117
  my $stable_id = shift;
118

119
  throw('Stable_id argument expected') if(!$stable_id);
120

121
  my $syn = $self->_tables()->[1];
122

123
124
  $self->bind_param_generic_fetch($stable_id,SQL_VARCHAR);
  my $pts = $self->generic_fetch("$syn.display_label = ?");
125

126
  return (@$pts) ? $pts->[0] : undef;
127
}
128

129

130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145

=head2 fetch_all_by_Slice

  Arg [1]    : Bio::EnsEMBL::Slice $slice
               The slice to fetch transcripts on.
  Arg [3]    : (optional) boolean $load_exons
               if true, exons will be loaded immediately rather than
               lazy loaded later.
  Example    : $transcripts = $
  Description: Overrides superclass method to optionally load exons
               immediately rather than lazy-loading them later.  This
               is more efficient when there are a lot of transcripts whose
               exons are going to be used.
  Returntype : reference to list of transcripts
  Exceptions : thrown if exon cannot be placed on transcript slice
  Caller     : Slice::get_all_Transcripts
146
  Status     : Stable
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168

=cut

sub fetch_all_by_Slice {
  my $self  = shift;
  my $slice = shift;
  my $logic_name = shift;
  my $load_exons = shift;

  my $transcripts = $self->SUPER::fetch_all_by_Slice($slice,$logic_name);

  # if there are 0 or 1 transcripts still do lazy-loading
  if(!$load_exons || @$transcripts < 2) {
    return $transcripts;
  }

  # preload all of the exons now, instead of lazy loading later
  # faster than 1 query per transcript

  # get extent of region spanned by transcripts
  my ($min_start, $max_end);
  foreach my $tr (@$transcripts) {
169
170
    if(!defined($min_start) || $tr->seq_region_start() < $min_start) {
      $min_start = $tr->seq_region_start();
171
    }
172
173
    if(!defined($max_end) || $tr->seq_region_end() > $max_end) {
      $max_end   = $tr->seq_region_end();
174
175
176
    }
  }

177
178
179
#  mades no sense, the limit for the slice will be defined by the transcripts
#  $min_start += $slice->start() - 1;
#  $max_end   += $slice->start() - 1;
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241

  my $ext_slice;

  if($min_start >= $slice->start() && $max_end <= $slice->end()) {
    $ext_slice = $slice;
  } else {
    my $sa = $self->db()->get_SliceAdaptor();
    $ext_slice = $sa->fetch_by_region
      ($slice->coord_system->name(), $slice->seq_region_name(),
       $min_start,$max_end, $slice->strand(), $slice->coord_system->version());
  }

  # associate exon identifiers with transcripts

  my %tr_hash = map {$_->dbID => $_} @$transcripts;

  my $tr_id_str = '(' . join(',', keys %tr_hash) . ')';

  my $sth = $self->prepare
    ("SELECT prediction_transcript_id, prediction_exon_id, exon_rank " .
     "FROM   prediction_exon " .
     "WHERE  prediction_transcript_id IN $tr_id_str");

  $sth->execute();

  my ($ex_id, $tr_id, $rank);
  $sth->bind_columns(\$tr_id, \$ex_id, \$rank);

  my %ex_tr_hash;

  while($sth->fetch()) {
    $ex_tr_hash{$ex_id} ||= [];
    push @{$ex_tr_hash{$ex_id}}, [$tr_hash{$tr_id}, $rank];
  }

  $sth->finish();

  my $ea = $self->db()->get_PredictionExonAdaptor();
  my $exons = $ea->fetch_all_by_Slice($ext_slice);

  # move exons onto transcript slice, and add them to transcripts
  foreach my $ex (@$exons) {
    $ex = $ex->transfer($slice) if($slice != $ext_slice);

    if(!$ex) {
      throw("Unexpected. PredictionExon could not be transfered onto " .
            "PredictionTranscript slice.");
    }

    foreach my $row (@{$ex_tr_hash{$ex->dbID()}}) {
      my ($tr, $rank) = @$row;
      $tr->add_Exon($ex, $rank);
    }
  }

  return $transcripts;
}





242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
=head2 _objs_from_sth

  Arg [1]    : DBI:st $sth 
               An executed DBI statement handle
  Arg [2]    : (optional) Bio::EnsEMBL::Mapper $mapper 
               An mapper to be used to convert contig coordinates
               to assembly coordinates.
  Arg [3]    : (optional) Bio::EnsEMBL::Slice $slice
               A slice to map the prediction transcript to.   
  Example    : $p_transcripts = $self->_objs_from_sth($sth);
  Description: Creates a list of Prediction transcripts from an executed DBI
               statement handle.  The columns retrieved via the statement 
               handle must be in the same order as the columns defined by the
               _columns method.  If the slice argument is provided then the
               the prediction transcripts will be in returned in the coordinate
               system of the $slice argument.  Otherwise the prediction 
               transcripts will be returned in the RawContig coordinate system.
  Returntype : reference to a list of Bio::EnsEMBL::PredictionTranscripts
  Exceptions : none
  Caller     : superclass generic_fetch
262
  Status     : Stable
263

264
=cut
265

266
sub _objs_from_sth {
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
  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 is used.
  #

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

  my @ptranscripts;
  my %analysis_hash;
  my %slice_hash;
  my %sr_name_hash;
  my %sr_cs_hash;

  my ($prediction_transcript_id,
      $seq_region_id,
      $seq_region_start,
      $seq_region_end,
      $seq_region_strand,
289
290
      $analysis_id,
      $display_label);
291
292
293
294
295
296

  $sth->bind_columns(\$prediction_transcript_id,
                     \$seq_region_id,
                     \$seq_region_start,
                     \$seq_region_end,
                     \$seq_region_strand,
297
298
                     \$analysis_id,
                     \$display_label);
299
300
301
302
303
304
305
306
307
308
309
310
311
312

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

315
316
317
318
  my $dest_slice_start;
  my $dest_slice_end;
  my $dest_slice_strand;
  my $dest_slice_length;
319
  my $dest_slice_sr_name;
320
  my $dest_slice_sr_id;
321
322
323
324
325
  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();
326
    $dest_slice_sr_name = $dest_slice->seq_region_name();
327
    $dest_slice_sr_id =  $dest_slice->get_seq_region_id();
328
  }
329

330
 FEATURE: while($sth->fetch()) {
331

332
333
334
    #get the analysis object
    my $analysis = $analysis_hash{$analysis_id} ||=
      $aa->fetch_by_dbID($analysis_id);
335
336
    #need to get the internal_seq_region, if present
    $seq_region_id = $self->get_seq_region_id_internal($seq_region_id);
337
338
339
340
341
342
343
    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();
344
    }
345

346
347
    my $sr_name = $sr_name_hash{$seq_region_id};
    my $sr_cs   = $sr_cs_hash{$seq_region_id};
348
349
350
351
352
353
    #
    # remap the feature coordinates to another coord system 
    # if a mapper was provided
    #
    if($mapper) {

354
      ($seq_region_id,$seq_region_start,$seq_region_end,$seq_region_strand) =
355
356
357
358
        $mapper->fastmap($sr_name, $seq_region_start, $seq_region_end,
			 $seq_region_strand, $sr_cs);

      #skip features that map to gaps or coord system boundaries
359
      next FEATURE if(!defined($seq_region_id));
360
361

      #get a slice in the coord system we just mapped to
362
363
364
365
366
367
368
369
#      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);
#      }
370
371
    }

372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
    #
    # 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;
	  $seq_region_strand *= -1;
	}
      }

389
390
      #throw away features off the end of the requested slice
      if($seq_region_end < 1 || $seq_region_start > $dest_slice_length ||
391
	( $dest_slice_sr_id ne $seq_region_id )) {
392
393
394
395
	next FEATURE;
      }


396
      $slice = $dest_slice;
397
    }
398

399
400
401
402
403
404
405
406
407
408
409
410
411
    # Finally, create the new PredictionTranscript.
    push( @ptranscripts,
          $self->_create_feature('Bio::EnsEMBL::PredictionTranscript', {
                                   '-start'    => $seq_region_start,
                                   '-end'      => $seq_region_end,
                                   '-strand'   => $seq_region_strand,
                                   '-adaptor'  => $self,
                                   '-slice'    => $slice,
                                   '-analysis' => $analysis,
                                   '-dbID' => $prediction_transcript_id,
                                   '-display_label' => $display_label
                                 } ) );

412
  }
413

414
  return \@ptranscripts;
415
416
417
418
419
420
}



=head2 store

421
422
423
424
425
  Arg [1]    : list of Bio::EnsEMBL::PredictionTranscript @pre_transcripts 
  Example    : $prediction_transcript_adaptor->store(@pre_transcripts);
  Description: Stores a list of given prediction transcripts in database. 
               Puts dbID and Adaptor into each object stored object.
  Returntype : none
426
427
  Exceptions : on wrong argument type 
  Caller     : general 
428
  Status     : Stable
429
430
431
432

=cut

sub store {
433
  my ( $self, @pre_transcripts ) = @_;
434

435
  my $ptstore_sth = $self->prepare
436
437
438
439
440
441
442
443
    (qq{INSERT INTO prediction_transcript (seq_region_id, seq_region_start,
                                           seq_region_end, seq_region_strand, 
                                           analysis_id, display_label)
        VALUES( ?, ?, ?, ?, ?, ?)});

  my $ptupdate_sth = $self->prepare
    (qq{UPDATE prediction_transcript SET display_label = ?
        WHERE  prediction_transcript_id = ?});
444

445
446
447
  my $db = $self->db();
  my $analysis_adaptor = $db->get_AnalysisAdaptor();
  my $pexon_adaptor = $db->get_PredictionExonAdaptor();
448

449
450
451
  FEATURE: foreach my $pt (@pre_transcripts) {
    if(!ref($pt) || !$pt->isa('Bio::EnsEMBL::PredictionTranscript')) {
      throw('Expected PredictionTranscript argument not [' . ref($pt).']');
452
    }
453
454
455
456
457

    #skip prediction transcripts that have already been stored
    if($pt->is_stored($db)) {
      warning('Not storing already stored prediction transcript '. $pt->dbID);
      next FEATURE;
458
    }
459
460
461
462
463
464
465
466

    #get analysis and store it if it is not in the db
    my $analysis = $pt->analysis();
    if(!$analysis) {
      throw('Prediction transcript must have analysis to be stored.');
    }
    if(!$analysis->is_stored($db)) {
      $analysis_adaptor->store($analysis);
467
468
    }

469
470
471
472
473
474
475
476
477
478
    #ensure that the transcript coordinates are correct, they may not be,
    #if somebody has done some exon coordinate juggling and not recalculated
    #the transcript coords.
    $pt->recalculate_coordinates();

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

    #store the prediction transcript
479
480
481
482
483
484
485
486
    $ptstore_sth->bind_param(1,$seq_region_id,SQL_INTEGER);
    $ptstore_sth->bind_param(2,$pt->start,SQL_INTEGER);
    $ptstore_sth->bind_param(3,$pt->end,SQL_INTEGER);
    $ptstore_sth->bind_param(4,$pt->strand,SQL_TINYINT);
    $ptstore_sth->bind_param(5,$analysis->dbID,SQL_INTEGER);
    $ptstore_sth->bind_param(6,$pt->display_label,SQL_VARCHAR);

    $ptstore_sth->execute();
487
488
489
490
491
492
493
494
495

    my $pt_id = $ptstore_sth->{'mysql_insertid'};
    $original->dbID($pt_id);
    $original->adaptor($self);

    #store the exons
    my $rank = 1;
    foreach my $pexon (@{$original->get_all_Exons}) {
      $pexon_adaptor->store($pexon, $pt_id, $rank++);
496
    }
497
498
499
500
501

    # if a display label was not defined autogenerate one
    if(!defined($pt->display_label())) {
      my $zeros = '0' x (11 - length($pt_id));
      my $display_label = uc($analysis->logic_name()) . $zeros . $pt_id;
502
503
504
      $ptupdate_sth->bind_param(1,$display_label,SQL_VARCHAR);
      $ptupdate_sth->bind_param(2,$pt_id,SQL_INTEGER);
      $ptupdate_sth->execute();
505
506
      $original->display_label($display_label);
    }
507
508
509
510
  }
}


511
512
513

=head2 remove

514
515
516
517
  Arg [1]    : Bio::EnsEMBL::PredictionTranscript $pt 
  Example    : $prediction_transcript_adaptor->remove($pt);
  Description: removes given prediction transcript $pt from database. 
  Returntype : none
518
  Exceptions : throws if argument not a  Bio::EnsEMBL::PredictionTranscript
519
  Caller     : general
520
  Status     : Stable
521
522
523

=cut

524
525
526
sub remove {
  my $self = shift;
  my $pre_trans = shift;
527
528
529
530
531
532
533

  if(!ref($pre_trans)||!$pre_trans->isa('Bio::EnsEMBL::PredictionTranscript')){
    throw('Expected PredictionTranscript argument.');
  }

  if(!$pre_trans->is_stored($self->db())) {
    warning('PredictionTranscript is not stored in this DB - not removing.');
534
535
536
    return;
  }

537
  #remove all associated prediction exons
Jan-hinnerk Vogel's avatar
Jan-hinnerk Vogel committed
538
  my $pexon_adaptor = $self->db()->get_PredictionExonAdaptor();
539
540
541
542
543
544
  foreach my $pexon (@{$pre_trans->get_all_Exons}) {
    $pexon_adaptor->remove($pexon);
  }

  #remove the prediction transcript
  my $sth = $self->prepare( "DELETE FROM prediction_transcript
545
                             WHERE prediction_transcript_id = ?" );
546
547
  $sth->bind_param(1,$pre_trans->dbID,SQL_INTEGER);
  $sth->execute();
548

549
550
551
  #unset the adaptor and internal id
  $pre_trans->dbID(undef);
  $pre_trans->adaptor(undef);
552
553
554
}


Glenn Proctor's avatar
Glenn Proctor committed
555
556
557
558
=head2 list_dbIDs

  Arg [1]    : none
  Example    : @feature_ids = @{$prediction_transcript_adaptor->list_dbIDs()};
559
560
  Description: Gets an array of internal ids for all prediction transcript
               features in the current db
561
  Arg[1]     : <optional> int. not 0 for the ids to be sorted by the seq_region.
Glenn Proctor's avatar
Glenn Proctor committed
562
563
564
  Returntype : list of ints
  Exceptions : none
  Caller     : ?
565
  Status     : Stable
Glenn Proctor's avatar
Glenn Proctor committed
566
567
568
569

=cut

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

572
   return $self->_list_dbIDs("prediction_transcript", undef, $ordered);
Glenn Proctor's avatar
Glenn Proctor committed
573
}
574

575
1;