PredictionTranscriptAdaptor.pm 14.9 KB
Newer Older
1
2
3
4
5
6
7
8
9
# EnsEMBL Exon reading writing adaptor for mySQL
#
# Author: Arne Stabenau
# 
# Date : 22.11.2001
#

=head1 NAME

10
11
Bio::EnsEMBL::DBSQL::PredictionTranscriptAdaptor -
Performs database interaction related to PredictionTranscripts
12
13
14

=head1 SYNOPSIS

15
16
17
18
19
#get a prediction transcript adaptor from the database
$pta = $database_adaptor->get_PredictionTranscriptAdaptor();

#get a slice on a region of chromosome 1
$sa = $database_adaptor->get_SliceAdaptor();
20
$slice = $sa->fetch_by_region('x', 100000, 200000);
21
22

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

25
26
=head1 CONTACT

27
Email questions to the EnsEMBL developer list: <ensembl-dev@ebi.ac.uk>
28
29
30
31
32
33
34
35

=cut

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

use vars qw( @ISA );
use strict;

36
use Bio::EnsEMBL::DBSQL::BaseFeatureAdaptor;
37
use Bio::EnsEMBL::DBSQL::AnalysisAdaptor;
38
use Bio::EnsEMBL::PredictionTranscript;
39
use Bio::EnsEMBL::Utils::Exception qw(deprecate throw warning);
40

41
@ISA = qw( Bio::EnsEMBL::DBSQL::BaseFeatureAdaptor );
42
43


44
45
46
47
48
49
50
51
52
# _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
53

54
sub _tables {
55
  my $self = shift;
56

57
  return ['prediction_transcript', 'pt'];
58
}
59

60

61
# _columns
62

63
64
65
66
67
68
69
70
71
#  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
#
72

73
sub _columns {
74
75
  my $self = shift;

76
77
78
79
80
81
  return qw( pt.prediction_transcript_id
             pt.seq_region_id
             pt.seq_region_start
             pt.seq_region_end
             pt.seq_region_strand
             pt.analysis_id);
82
83
}

84

85
86
87
88
=head2 fetch_by_stable_id

Arg [1]    : string $stable_id
             The stable id of the transcript to retrieve
89
Example    : $trans = $trans_adptr->fetch_by_stable_id('GENSCAN00000001234');
90
Description: Retrieves a prediction transcript via its stable id.  Note that
91
92
93
             the stable id is not actually stored in the database and is
             calculated upon retrieval as the logic_name of the transcripts
             analysis concatted with the '0'-left-padded internal identifier.
94
95
96
97
98
99
100
101
             prediction transcript
Returntype : Bio::EnsEMBL::PredictionTranscript
Caller     : general

=cut

sub fetch_by_stable_id {
  my $self = shift;
102
  my $stable_id = shift;
103

104
  throw('Stable_id argument expected') if(!$stable_id);
105

106
107
  #take the last 11 digits
  my $dbID = substr($stable_id, -11);
108

109
110
  #cast to an integer (remove lpadded 0s):
  $dbID += 0;
111

112
  return $self->fetch_by_dbID($dbID);
113
}
114

115

116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
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

=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

=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) {
    if(!defined($min_start) || $tr->start() < $min_start) {
      $min_start = $tr->start();
    }
    if(!defined($max_end) || $tr->end() > $max_end) {
      $max_end   = $tr->end();
    }
  }

  $min_start += $slice->start() - 1;
  $max_end   += $slice->start() - 1;

  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;
}





226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
=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
246

247
=cut
248

249
sub _objs_from_sth {
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
  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,
      $analysis_id );

  $sth->bind_columns(\$prediction_transcript_id,
                     \$seq_region_id,
                     \$seq_region_start,
                     \$seq_region_end,
                     \$seq_region_strand,
                     \$analysis_id );

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

296
297
298
299
300
301
302
303
304
305
  my $dest_slice_start;
  my $dest_slice_end;
  my $dest_slice_strand;
  my $dest_slice_length;
  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();
  }
306

307
 FEATURE: while($sth->fetch()) {
308

309
310
311
    #get the analysis object
    my $analysis = $analysis_hash{$analysis_id} ||=
      $aa->fetch_by_dbID($analysis_id);
312

313
314
315
316
317
318
319
    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();
320
    }
321

322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
    #
    # remap the feature coordinates to another coord system 
    # if a mapper was provided
    #
    if($mapper) {
      my $sr_name = $sr_name_hash{$seq_region_id};
      my $sr_cs   = $sr_cs_hash{$seq_region_id};

      ($sr_name,$seq_region_start,$seq_region_end,$seq_region_strand) =
        $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
      next FEATURE if(!defined($sr_name));

      #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{"NAME:$sr_name:$cmp_cs_name:$cmp_cs_vers"} ||=
          $sa->fetch_by_region($cmp_cs_name, $sr_name,undef, undef, undef,
                               $cmp_cs_vers);
342
      } else {
343
344
345
        $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);
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
    #
    # 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;
	}
	
	#throw away features off the end of the requested slice
	if($seq_region_end < 1 || $seq_region_start > $dest_slice_length) {
	  next FEATURE;
	}
      }

      $slice = $dest_slice;
372
    }
373
374
375
376
377
378
379
380
381
382

    #finally, create the new repeat feature
    push @ptranscripts, Bio::EnsEMBL::PredictionTranscript->new
      ( '-start'         =>  $seq_region_start,
        '-end'           =>  $seq_region_end,
        '-strand'        =>  $seq_region_strand,
        '-adaptor'       =>  $self,
        '-slice'         =>  $slice,
        '-analysis'      =>  $analysis,
        '-dbID'          =>  $prediction_transcript_id );
383
  }
384

385
  return \@ptranscripts;
386
387
388
389
390
391
}



=head2 store

392
393
394
395
396
  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
397
398
  Exceptions : on wrong argument type 
  Caller     : general 
399
400
401
402

=cut

sub store {
403
  my ( $self, @pre_transcripts ) = @_;
404

405
406
407
408
  my $ptstore_sth = $self->prepare
    ("INSERT INTO prediction_transcript (seq_region_id, seq_region_start, " .
     "                    seq_region_end, seq_region_strand, analysis_id) " .
     "VALUES( ?, ?, ?, ?, ?)");
409

410
411
412
  my $db = $self->db();
  my $analysis_adaptor = $db->get_AnalysisAdaptor();
  my $pexon_adaptor = $db->get_PredictionExonAdaptor();
413

414
415
416
  FEATURE: foreach my $pt (@pre_transcripts) {
    if(!ref($pt) || !$pt->isa('Bio::EnsEMBL::PredictionTranscript')) {
      throw('Expected PredictionTranscript argument not [' . ref($pt).']');
417
    }
418
419
420
421
422

    #skip prediction transcripts that have already been stored
    if($pt->is_stored($db)) {
      warning('Not storing already stored prediction transcript '. $pt->dbID);
      next FEATURE;
423
    }
424
425
426
427
428
429
430
431

    #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);
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
    #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
    $ptstore_sth->execute($seq_region_id,
                          $pt->start(),
                          $pt->end(),
                          $pt->strand(),
                          $analysis->dbID());

    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++);
458
459
460
461
462
    }
  }
}


463
464
465

=head2 remove

466
467
468
469
470
471
  Arg [1]    : Bio::EnsEMBL::PredictionTranscript $pt 
  Example    : $prediction_transcript_adaptor->remove($pt);
  Description: removes given prediction transcript $pt from database. 
  Returntype : none
  Exceptions : none 
  Caller     : general
472
473
474

=cut

475
476
477
sub remove {
  my $self = shift;
  my $pre_trans = shift;
478
479
480
481
482
483
484

  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.');
485
486
487
    return;
  }

488
489
490
491
492
493
494
495
  #remove all associated prediction exons
  my $pexon_adaptor = $self->get_PredictionExonAdaptor();
  foreach my $pexon (@{$pre_trans->get_all_Exons}) {
    $pexon_adaptor->remove($pexon);
  }

  #remove the prediction transcript
  my $sth = $self->prepare( "DELETE FROM prediction_transcript
496
                             WHERE prediction_transcript_id = ?" );
497
498
  $sth->execute( $pre_trans->dbID );

499
500
501
  #unset the adaptor and internal id
  $pre_trans->dbID(undef);
  $pre_trans->adaptor(undef);
502
503
504
}


Glenn Proctor's avatar
Glenn Proctor committed
505
506
507
508
=head2 list_dbIDs

  Arg [1]    : none
  Example    : @feature_ids = @{$prediction_transcript_adaptor->list_dbIDs()};
509
510
  Description: Gets an array of internal ids for all prediction transcript
               features in the current db
Glenn Proctor's avatar
Glenn Proctor committed
511
512
513
514
515
516
517
518
519
520
521
  Returntype : list of ints
  Exceptions : none
  Caller     : ?

=cut

sub list_dbIDs {
   my ($self) = @_;

   return $self->_list_dbIDs("prediction_transcript");
}
522

523
1;