Translation.pm 17.4 KB
Newer Older
1
=head1 LICENSE
2

3
  Copyright (c) 1999-2011 The European Bioinformatics Institute and
4
5
6
7
8
9
10
11
12
  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
  developers list at <dev@ensembl.org>.
  Questions may also be sent to the Ensembl help desk at
  <helpdesk@ensembl.org>.
13

14
15
16
17
18
19
=cut
=head1 NAME
=head1 SYNOPSIS
=head1 DESCRIPTION
=head1 METHODS
=cut
20

21
22
23
24
25
26
27
package Bio::EnsEMBL::Utils::VegaCuration::Translation;
use strict;
use warnings;
use vars qw(@ISA);
use Data::Dumper;
use Bio::EnsEMBL::Utils::VegaCuration::Transcript;
@ISA = qw(Bio::EnsEMBL::Utils::VegaCuration::Transcript);
28

29
=head2 check_CDS_start_end_remarks
30

31
32
33
34
35
   Args       : B::E::Transcript
   Example    : my $results = $support->check_CDS_end_remarks($transcript)
   Description: identifies incorrect 'CDS end...' transcript remarks in a
                otter-derived Vega database
   Returntype : hashref
36

37
=cut
38

39
sub check_CDS_start_end_remarks {
40
  my $self = shift; 
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
  my $trans = shift;
  # info for checking
  my @remarks = @{$trans->get_all_Attributes('remark')};
  my $coding_end   = $trans->cdna_coding_end;
  my $coding_start = $trans->cdna_coding_start;
  my $trans_end    = $trans->length;
  my $trans_seq    = $trans->seq->seq;
  my $stop_codon   = substr($trans_seq, $coding_end-3, 3);
  my $start_codon  = substr($trans_seq, $coding_start-1, 3);
  #hashref to return results
  my $results;
  
  #extra CDS end not found remarks
  if (grep {$_->value eq 'CDS end not found'} @remarks) {
    if (   ($coding_end != $trans_end) 
	     && ( grep {$_ eq $stop_codon} qw(TGA TAA TAG) ) ) {
      $results->{'END_EXTRA'} = 1;
    }
  }
  #missing CDS end not found remark
  if ( $coding_end == $trans_end ) {
    if (! grep {$_->value eq 'CDS end not found'} @remarks) {
      if (grep {$_ eq $stop_codon} qw(TGA TAA TAG)) {
	$results->{'END_MISSING_2'} = 1;
      }
      else {
	$results->{'END_MISSING_1'} = $stop_codon;
      }
    }
  }
  #extra CDS start not found remark
  if (grep {$_->value eq 'CDS start not found'} @remarks) {
    if (   ($coding_start != 1) 
	     && ($start_codon eq 'ATG') ) {
      $results->{'START_EXTRA'} = 1;
    }
  }
  #missing CDS start not found remark
  if ( $coding_start == 1) {
    if ( ! grep {$_->value eq 'CDS start not found'} @remarks) {
      if ($start_codon eq 'ATG') {
	$results->{'START_MISSING_2'} = 1;
      } else {
	$results->{'START_MISSING_1'} = $start_codon;
      }
    }
  }
  return $results;
}
90

91
=head2 check_CDS_start_end_remarks_loutre
92

93
94
95
96
97
   Args       : B::E::Transcript
   Example    : my $results = $support->check_CDS_end_remarks($transcript)
   Description: identifies incorrect 'CDS end...' transcript attribs in a loutre
                of a loutre-derived Vega database.
   Returntype : hashref
98

99
=cut
100

101
102
103
sub check_CDS_start_end_remarks_loutre {
  my $self = shift;
  my $trans = shift;
104

105
106
107
108
  # info for checking
  my @stops = qw(TGA TAA TAG);
  my %attributes;
  foreach my $attribute (@{$trans->get_all_Attributes()}) {
109
    push @{$attributes{$attribute->code}}, $attribute;
110
111
112
113
114
115
116
117
  }
#  warn $trans->stable_id;
#  warn Data::Dumper::Dumper(\%attributes);
  my $coding_end   = $trans->cdna_coding_end;
  my $coding_start = $trans->cdna_coding_start;
  my $trans_end    = $trans->length;
  my $trans_seq    = $trans->seq->seq;
  my $stop_codon_offset = 3 + $trans->translation->end_Exon->end_phase;
118
  my $initial_exon_phase = $trans->translation->start_Exon->phase;
119
120
  my $stop_codon   = substr($trans_seq, $coding_end-3, 3);
  my $start_codon  = substr($trans_seq, $coding_start-1, 3);
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135

  my $start_codon_incorrect = 1;
  if ($start_codon eq 'ATG' ) {
    $start_codon_incorrect = 0;
  }
  elsif ($start_codon eq 'CTG') {
    foreach my $attrib (@{$attributes{'remark'}}) {
      if ($attrib->value =~ /non[- ]ATG start/) {
	$start_codon_incorrect = 0;
      }
    }
  }

#  warn "$start_codon -- $initial_exon_phase -- $coding_start -- $start_codon_incorrect";

136
137
  #hashref to return results
  my $results;
138

139
140
  #extra CDS end not found remarks
  if ($attributes{'cds_end_NF'}) {
141
    if ( ($attributes{'cds_end_NF'}->[0]->value == 1)
142
143
144
145
146
147
148
149
150
151
	   && ($coding_end != $trans_end) 
	   && ( grep {$_ eq $stop_codon} @stops) ) {
#      warn $trans->stable_id.": $coding_end--$trans_end--$stop_codon";
#      warn $trans->translation->end_Exon->end_phase;
      $results->{'END_EXTRA'} = $stop_codon;
    }
  }
  #missing CDS end not found remark
  if ( $coding_end == $trans_end ) {
    if ($attributes{'cds_end_NF'}) {
152
      if ($attributes{'cds_end_NF'}->[0]->value == 0 ) {
153
154
155
156
157
158
159
160
161
162
163
164
	if (! grep {$_ eq $stop_codon} @stops) {
#	  warn $trans->translation->end_Exon->end_phase;
	  $results->{'END_MISSING'}{'WRONG'} = $stop_codon;
	}
      }
    }
    elsif (! grep {$_ eq $stop_codon} @stops) {
      $results->{'END_MISSING'}{'ABSENT'} = $stop_codon;
    }
  }
  #extra CDS start not found remark 
  if ( $attributes{'cds_start_NF'}) {
165
166
167
168
169
    if (   ($attributes{'cds_start_NF'}->[0]->value == 1 )
        && (! $start_codon_incorrect)) {
      unless ( ($coding_start == 1) && ( $initial_exon_phase > 0)) {
	$results->{'START_EXTRA'} = $start_codon;
      }
170
171
172
173
174
    }
  }
  #missing CDS start not found remark
  if ( $coding_start == 1) {
    if ( $attributes{'cds_start_NF'} ) {
175
176
177
178
179
180
      if ( $attributes{'cds_start_NF'}->[0]->value == 0 ) {
	if ($start_codon_incorrect) {
	  $results->{'START_MISSING'}{'ABSENT'} = $start_codon;
	}
	elsif ($initial_exon_phase > 0) {
	  $results->{'START_MISSING'}{'INITIAL_PHASE'} = $initial_exon_phase;
181
182
183
184
185
186
	}
      }
    }
    elsif ($start_codon ne 'ATG') {
      $results->{'START_MISSING'}{'ABSENT'} = $start_codon;
    }
187

188
189
190
  }
  return $results;
}
191

192
=head2 get_havana_seleno_comments
193

194
195
196
197
   Args       : none
   Example    : my $results = $support->get_havana_seleno_comments
   Description: parses the HEREDOC containing Havana comments in this module
   Returntype : hashref
198

199
=cut
200

201
202
203
204
205
206
207
sub get_havana_seleno_comments {
  my $seen_translations;
  while (<DATA>) {
    next if /^\s+$/ or /#+/;
    my ($obj,$comment) = split /=/;
    $obj =~ s/^\s+|\s+$//g;
    $comment =~ s/^\s+|\s+$//g;
208
209
210
    # We add the origin as now "seen" can come from a number of places, and have
    # a number of consequences in different cases, not just discounted Secs from this method. -- ds23
    $seen_translations->{$obj} = [ $comment,"notsec-havana" ];
211
212
213
  }
  return $seen_translations;
}
Maurice Hendrix's avatar
bug fix    
Maurice Hendrix committed
214

215
216
217
218
219
220
sub check_for_stops {
  my $support = shift;
  my ($gene,$seen_transcripts,$log_object) = @_;
  my $transcripts;
  my $has_log_object=defined($log_object);
  if($has_log_object){
221
222
    my @help = $log_object->species_params->get_trans($gene->stable_id);
    $transcripts=\@help;
223
224
  }else{
    $log_object=$support;
225
    $transcripts=$gene->get_all_Transcripts;
226
  }
227

228
229
230
231
232
233
  my $gname = $gene->get_all_Attributes('name')->[0]->value;
  my $gsi = $gene->stable_id;
  my $scodon = 'TGA';
  my $mod_date = $support->date_format( $gene->modified_date,'%d/%m/%y' );
  my $hidden_remak_ttributes;
 TRANS:
234
  foreach my $trans (@$transcripts) {
235
236
237
238
239
240
241
242
243
244
    my $tsi = $trans->stable_id;
    my $tID = $trans->dbID;
    my $tname = $trans->get_all_Attributes('name')->[0]->value;
    if($has_log_object){
      $hidden_remak_ttributes=$log_object->species_params->get_attributes->{$tsi}->{'hidden_remark'};
    }else{
      $hidden_remak_ttributes=$trans->get_all_Attributes('hidden_remark');
    }
    foreach my $rem (@$hidden_remak_ttributes) {
      if ($rem->value =~ /not_for_Vega/) {
245
        #$support->log_verbose("Skipping transcript $tname ($tsi) since 'not_for_Vega'\n",1);
246
        $log_object->_save_log('log_verbose', '', $gsi, '', $tsi, '', "Skipping transcript $tname ($tsi) since 'not_for_Vega'");
247
248
249
250
251
252
253
254
255
        next TRANS;
      }
    }

    # go no further if there is a ribosomal framshift attribute
    foreach my $attrib (@{$trans->get_all_Attributes('_rib_frameshift')}) {
      if ($attrib->value) {
        $log_object->_save_log('log', '', $gsi, '', $tsi, '', "Skipping $tsi ($tname) since it has a ribosomal frameshift attribute");
        next TRANS;
256
257
      }
    }
258

259
    #$support->log_verbose("Studying transcript $tsi ($tname, $tID)\n",1);
260
    $log_object->_save_log('log_verbose', '', $gsi, '', $tsi, '', "Studying transcript $tsi ($tname, $tID)");
261
262
263
264
265
266
267
    my $peptide;
		
    # go no further if the transcript doesn't translate or if there are no stops
    next TRANS unless ($peptide = $trans->translate);
    my $pseq = $peptide->seq;
    my $orig_seq = $pseq;
    # (translate method trims stops from sequence end)
268

269
    next TRANS unless ($pseq =~ /\*/);
270
    # warn sprintf("Stop codon is '%s'\n",substr($trans->translateable_seq,-3));
271
    #$support->log_verbose("Stops found in $tsi ($tname)\n",1);
272
    $log_object->_save_log('log_verbose', '', $gsi, '', $tsi, '', "Stops found in $tsi ($tname)");
273

274
275
276
277
278
    # find out where and how many stops there are
    my @found_stops;
    my $mrna   = $trans->translateable_seq;
    my $offset = 0;
    my $tstop;
279
    while ($pseq =~ /^([^\*]*)\*(.*)/) {
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
      my $pseq1_f = $1;
      $pseq = $2;
      my $seq_flag = 0;
      $offset += length($pseq1_f) * 3;
      my $stop = substr($mrna, $offset, 3);
      my $aaoffset = int($offset / 3)+1;
      push(@found_stops, [ $stop, $aaoffset ]);
      $tstop .= "$aaoffset ";
      $offset += 3;
    }
		
    # are all stops TGA...?
    my $num_stops = scalar(@found_stops);
    my $num_tga = 0;
    my $positions;
    foreach my $stop (@found_stops) {
      $positions .= $stop->[0]."(".$stop->[1].") ";
      if ($stop->[0] eq $scodon) {
	$num_tga++;
      }
    }
    my $source = $gene->source;
    #...no - an internal stop codon error in the database...
    if ($num_tga < $num_stops) {
      if ($source eq 'havana') {
305
        #$support->log_warning("INTERNAL STOPS HAVANA: Transcript $tsi ($tname) from gene $gname has non \'$scodon\' stop codons [$mod_date]:\nSequence = $orig_seq\nStops at $positions)\n\n");
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
        $log_object->_save_log('log_warning', '', $gsi, 'TRANSCRIPT', $tsi, 'VQCT_internal_stop', "INTERNAL STOPS HAVANA: Transcript $tsi ($tname) from gene $gname has non \'$scodon\' stop codons [$mod_date]: Sequence = $orig_seq Stops at $positions)");
      }
      else {
        #$support->log_warning("INTERNAL STOPS EXTERNAL: Transcript $tsi ($tname) from gene $gname has non \'$scodon\' stop codons[$mod_date]:\nSequence = $orig_seq\nStops at $positions)\n\n");
        $log_object->_save_log('log_warning', '', $gsi, 'TRANSCRIPT', $tsi, 'VQCT_internal_stop', "INTERNAL STOPS EXTERNAL: Transcript $tsi ($tname) from gene $gname has non \'$scodon\' stop codons[$mod_date]: Sequence = $orig_seq Stops at $positions)");  
      }
    }
    #...yes - check remarks
    else {
      my $flag_remark  = 0; # 1 if word seleno has been used
      my $flag_remark2 = 0; # 1 if existing remark has correct numbering
      my $alabel       = 'Annotation_remark- selenocysteine ';
      my $alabel2      = 'selenocysteine ';
      my $annot_stops;
      my $remarks;
      my $att;
      #get both hidden_remarks and remarks
      foreach my $remark_type ('remark','hidden_remark') {
324
        if($has_log_object){
325
326
327
328
          $att=$log_object->species_params->get_attributes->{$trans->stable_id}->{$remark_type};
        }else{
          $att=$trans->get_all_Attributes($remark_type)
        }
329
330
331
        foreach my $attrib ( @$att) {
          push @{$remarks->{$remark_type}}, $attrib->value;
        }
332
333
334
      }
      #parse remarks to check syntax for location of edits
      while (my ($attrib,$remarks)= each %$remarks) {
335
336
337
        foreach my $text (@{$remarks}) {					
          if ( ($attrib eq 'remark') && ($text=~/^$alabel(.*)/) ){
            #$support->log_warning("seleno remark for $tsi stored as Annotation_remark not hidden remark) [$mod_date]\n");
Maurice Hendrix's avatar
bug fix    
Maurice Hendrix committed
338
            $log_object->_save_log('log_warning', '', $gsi, '', $tsi, 'VQCT_wrong_selC_coord', "seleno remark for $tsi stored as Annotation_remark not hidden remark) [$mod_date]");        
339
340
341
342
343
344
            $annot_stops=$1;
          }
          elsif ($text =~ /^$alabel2(.*)/) {
            $annot_stops=$1;
          }
        }
345
      }
Maurice Hendrix's avatar
bug fix    
Maurice Hendrix committed
346

347
348
349
      #check the location of the annotated edits matches actual stops in the sequence
      my @annotated_stops;
      if ($annot_stops){
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
        my $i = 0;
        foreach my $offset (split(/\s+/, $annot_stops)) {
          if ($offset !~ /^\d+$/) {
            $log_object->_save_log('log', '', $gene->stable_id, '', $tsi, '', "Annotated stop for transcript tsi ($tname) is not an integer \"$offset\" - might need to investigate [$mod_date]");
          }
          #OK if it matches a known stop
          elsif (
            defined($found_stops[$i]) && defined($found_stops[$i]->[1]) && ($found_stops[$i]->[1] == $offset)) {
            push  @annotated_stops, $offset;
          }
          # catch old annotations where number was in DNA not peptide coordinates
          elsif (defined($found_stops[$i]) && defined($found_stops[$i]->[1]) && (($found_stops[$i]->[1] * 3) == $offset)) {
            $log_object->_save_log('log_warning', '', $gene->stable_id, 'DNA', $tsi, 'VQCT_wrong_selC_coord', "DNA: Annotated stop for transcript tsi ($tname) is in DNA not peptide coordinates) [$mod_date]");
          }
          # catch old annotations where number off by one
          elsif (defined($found_stops[$i]) && defined($found_stops[$i]->[1]) && (($found_stops[$i]->[1]) == $offset+1)) {
            $log_object->_save_log('log_warning', '', $gene->stable_id, 'PEPTIDE', $tsi, 'VQCT_wrong_selC_coord', "PEPTIDE: Annotated stop for transcript $tsi ($tname) is out by one) [$mod_date]");
          }
Maurice Hendrix's avatar
bug fix    
Maurice Hendrix committed
368
          elsif (defined($offset)  && ($offset=~/^\d+$/)){
369
            if ($offset == length($orig_seq)+1) {
370
              if($seen_transcripts->{$tsi} && $seen_transcripts->{$tsi}->[1] eq 'known-tga-stop') {
371
                $log_object->_save_log('log', '', $gene->stable_id, 'TRANSCRIPT', $tsi, '', "Annotated stop for transcript $tsi ($tname) known to be a stop codon. Ok. [$mod_date]");
372
373
              } elsif($seen_transcripts->{$tsi} && $seen_transcripts->{$tsi}->[1] eq 'known-terminal-sec') {
                $log_object->_save_log('log', '', $gene->stable_id, 'TRANSCRIPT', $tsi, '', "Annotated stop for transcript $tsi ($tname) known to be a terminal Sec. Ok. [$mod_date]");
374
              } else {
375
                $log_object->_save_log('log_warning', '', $gene->stable_id, 'TRANSCRIPT', $tsi, '', "Annotated stop for transcript $tsi ($tname) \"$offset\" matches actual stop codon yet has no entry in script config to disambiguate it. Please investigate and add appropriate entry to config arrays in add_selcys.pl. [$mod_date]");
376
              }
377
378
379
380
381
382
383
384
            }
            else {
              $log_object->_save_log('log_warning', '', $gene->stable_id, 'TRANSCRIPT', $tsi, 'VQCT_wrong_selC_coord', "Annotated stop for transcript $tsi ($tname) \"$offset\" does not match a TGA codon) [$mod_date]");
              push  @annotated_stops, $offset;
            }
          }						
          $i++;
        }
385
      }
Maurice Hendrix's avatar
bug fix    
Maurice Hendrix committed
386

387
388
      #check location of found stops matches annotated ones - any new ones are reported
      foreach my $stop ( @found_stops ) {
389
390
391
        my $pos = $stop->[1];
        my $seq = $stop->[0];
        unless ( grep { $pos == $_} @annotated_stops) {
392
          if ($seen_transcripts->{$tsi} && $seen_transcripts->{$tsi}->[1] eq 'notsec-havana') {
393
            #$support->log_verbose("Transcript $tsi ($tname) has potential selenocysteines but has been discounted by annotators:\n\t".$seen_transcripts->{$tsi}.") [$mod_date]\n");
394
            $log_object->_save_log('log_verbose', '', $gene->stable_id, '', $tsi, 'VQCT_pot_selC', "Transcript $tsi ($tname) has potential selenocysteines but has been discounted by annotators: ".$seen_transcripts->{$tsi}->[0].") [$mod_date]");
395
396
397
          }
          else {
            #$support->log("POTENTIAL SELENO ($seq) in $tsi ($tname, gene $gname) found at $pos [$mod_date]\n");
Maurice Hendrix's avatar
bug fix    
Maurice Hendrix committed
398
            $log_object->_save_log('log', '', $gene->stable_id, '', $tsi, 'VQCT_pot_selC', "POTENTIAL SELENO ($seq) in $tsi ($tname, gene $gname) found at $pos [$mod_date]");
399
400
          }
        }
401
402
403
404
405
406
407
408
409
410
411
412
413
      }
    }
  }
}
sub _save_log{
  my $self=shift;
  my $log_type = shift;
  my $chrom_name=shift || '';
  my $gsi=shift || '';
  my $type=shift || '';  
  my $tsi=shift || '';  
  my $tag=shift || '';
  my $txt=shift || '';
414
  $self->$log_type($txt."\n");
415
}
Maurice Hendrix's avatar
bug fix    
Maurice Hendrix committed
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
#details of annotators comments
__DATA__
OTTHUMT00000144659 = FIXED- changed to transcript
OTTHUMT00000276377 = FIXED- changed to transcript
OTTHUMT00000257741 = FIXED- changed to nmd
OTTHUMT00000155694 = NOT_FIXED- should be nmd but external annotation but cannot be fixed
OTTHUMT00000155695 = NOT_FIXED- should be nmd but external annotation but cannot be fixed
OTTHUMT00000282573 = FIXED- changed to unprocessed pseudogene
OTTHUMT00000285227 = FIXED- changed start site
OTTHUMT00000151008 = FIXED- incorrect trimming of CDS, removed extra stop codon
OTTHUMT00000157999 = FIXED- changed incorrect stop
OTTHUMT00000150523 = FIXED- incorrect trimming of CDS
OTTHUMT00000150525 = FIXED- incorrect trimming of CDS
OTTHUMT00000150522 = FIXED- incorrect trimming of CDS
OTTHUMT00000150521 = FIXED- incorrect trimming of CDS
OTTHUMT00000246819 = FIXED- corrected frame
OTTHUMT00000314078 = FIXED- corrected frame
OTTHUMT00000080133 = FIXED- corrected frame
OTTHUMT00000286423 = FIXED- changed to transcript
OTTMUST00000055509 = FIXED- error
OTTMUST00000038729 = FIXED- corrected frame
OTTMUST00000021760 = FIXED- corrected frame
OTTMUST00000023057 = FIXED- corrected frame
OTTMUST00000015207 = FIXED- corrected frame
OTTMUST00000056646 = FIXED- error
OTTMUST00000059686 = FIXED- corrected frame
OTTMUST00000013426 = FIXED- corrected frame
444
OTTMUST00000044412 = FIXED- error