|
0
|
1 =head1 LICENSE
|
|
|
2
|
|
|
3 Copyright (c) 1999-2012 The European Bioinformatics Institute and
|
|
|
4 Genome Research Limited. All rights reserved.
|
|
|
5
|
|
|
6 This software is distributed under a modified Apache license.
|
|
|
7 For license details, please see
|
|
|
8
|
|
|
9 http://www.ensembl.org/info/about/code_licence.html
|
|
|
10
|
|
|
11 =head1 CONTACT
|
|
|
12
|
|
|
13 Please email comments or questions to the public Ensembl
|
|
|
14 developers list at <dev@ensembl.org>.
|
|
|
15
|
|
|
16 Questions may also be sent to the Ensembl help desk at
|
|
|
17 <helpdesk@ensembl.org>.
|
|
|
18
|
|
|
19 =cut
|
|
|
20
|
|
|
21 =head1 NAME
|
|
|
22
|
|
|
23 Bio::EnsEMBL::Variation::TranscriptVariationAllele
|
|
|
24
|
|
|
25 =head1 SYNOPSIS
|
|
|
26
|
|
|
27 use Bio::EnsEMBL::Variation::TranscriptVariationAllele;
|
|
|
28
|
|
|
29 my $tva = Bio::EnsEMBL::Variation::TranscriptVariationAllele->new(
|
|
|
30 -transcript_variation => $tv,
|
|
|
31 -variation_feature_seq => 'A',
|
|
|
32 -is_reference => 0,
|
|
|
33 );
|
|
|
34
|
|
|
35 print "sequence with respect to the transcript: ", $tva->feature_seq, "\n";
|
|
|
36 print "sequence with respect to the variation feature: ", $tva->variation_feature_seq, "\n";
|
|
|
37 print "consequence SO terms: ", (join ",", map { $_->SO_term } @{ $tva->get_all_OverlapConsequences }), "\n";
|
|
|
38 print "amino acid change: ", $tva->peptide_allele_string, "\n";
|
|
|
39 print "resulting codon: ", $tva->codon, "\n";
|
|
|
40 print "reference codon: ", $tva->transcript_variation->get_reference_TranscriptVariationAllele->codon, "\n";
|
|
|
41 print "PolyPhen prediction: ", $tva->polyphen_prediction, "\n";
|
|
|
42 print "SIFT prediction: ", $tva->sift_prediction, "\n";
|
|
|
43
|
|
|
44 =head1 DESCRIPTION
|
|
|
45
|
|
|
46 A TranscriptVariationAllele object represents a single allele of a TranscriptVariation.
|
|
|
47 It provides methods that are specific to the sequence of the allele, such as codon,
|
|
|
48 peptide etc. Methods that depend only on position (e.g. CDS start) will be found in
|
|
|
49 the associated TranscriptVariation. Ordinarily you will not create these objects
|
|
|
50 yourself, but instead you would create a TranscriptVariation object which will then
|
|
|
51 construct TranscriptVariationAlleles based on the allele string of the associated
|
|
|
52 VariationFeature.
|
|
|
53
|
|
|
54 Note that any methods that are not specific to Transcripts will be found in the
|
|
|
55 VariationFeatureOverlapAllele superclass.
|
|
|
56
|
|
|
57 =cut
|
|
|
58
|
|
|
59 package Bio::EnsEMBL::Variation::TranscriptVariationAllele;
|
|
|
60
|
|
|
61 use strict;
|
|
|
62 use warnings;
|
|
|
63
|
|
|
64 use Bio::EnsEMBL::Variation::ProteinFunctionPredictionMatrix qw($AA_LOOKUP);
|
|
|
65 use Bio::EnsEMBL::Utils::Exception qw(throw warning deprecate);
|
|
|
66 use Bio::EnsEMBL::Variation::Utils::Sequence qw(hgvs_variant_notation format_hgvs_string);
|
|
|
67 use Bio::EnsEMBL::Utils::Sequence qw(reverse_comp);
|
|
|
68 use Bio::EnsEMBL::Variation::Utils::VariationEffect qw(within_cds within_intron stop_lost affects_start_codon);
|
|
|
69
|
|
|
70 use base qw(Bio::EnsEMBL::Variation::VariationFeatureOverlapAllele Bio::EnsEMBL::Variation::BaseTranscriptVariationAllele);
|
|
|
71
|
|
|
72
|
|
|
73 our $DEBUG = 0;
|
|
|
74
|
|
|
75 sub new_fast {
|
|
|
76 my ($class, $hashref) = @_;
|
|
|
77
|
|
|
78 # swap a transcript_variation argument for a variation_feature_overlap one
|
|
|
79
|
|
|
80 if ($hashref->{transcript_variation}) {
|
|
|
81 $hashref->{variation_feature_overlap} = delete $hashref->{transcript_variation};
|
|
|
82 }
|
|
|
83
|
|
|
84 # and call the superclass
|
|
|
85
|
|
|
86 return $class->SUPER::new_fast($hashref);
|
|
|
87 }
|
|
|
88
|
|
|
89 =head2 transcript_variation
|
|
|
90
|
|
|
91 Description: Get/set the associated TranscriptVariation
|
|
|
92 Returntype : Bio::EnsEMBL::Variation::TranscriptVariation
|
|
|
93 Exceptions : throws if the argument is the wrong type
|
|
|
94 Status : At Risk
|
|
|
95
|
|
|
96 =cut
|
|
|
97
|
|
|
98 sub transcript_variation {
|
|
|
99 my ($self, $tv) = @_;
|
|
|
100 assert_ref($tv, 'Bio::EnsEMBL::Variation::TranscriptVariation') if $tv;
|
|
|
101 return $self->variation_feature_overlap($tv);
|
|
|
102 }
|
|
|
103
|
|
|
104 =head2 variation_feature
|
|
|
105
|
|
|
106 Description: Get the associated VariationFeature
|
|
|
107 Returntype : Bio::EnsEMBL::Variation::VariationFeature
|
|
|
108 Exceptions : none
|
|
|
109 Status : At Risk
|
|
|
110
|
|
|
111 =cut
|
|
|
112
|
|
|
113 sub variation_feature {
|
|
|
114 my $self = shift;
|
|
|
115 return $self->transcript_variation->variation_feature;
|
|
|
116 }
|
|
|
117
|
|
|
118 =head2 pep_allele_string
|
|
|
119
|
|
|
120 Description: Return a '/' delimited string of the reference peptide and the
|
|
|
121 peptide resulting from this allele, or a single peptide if this
|
|
|
122 allele does not change the peptide (e.g. because it is synonymous)
|
|
|
123 Returntype : string or undef if this allele is not in the CDS
|
|
|
124 Exceptions : none
|
|
|
125 Status : At Risk
|
|
|
126
|
|
|
127 =cut
|
|
|
128
|
|
|
129 sub pep_allele_string {
|
|
|
130 my ($self) = @_;
|
|
|
131
|
|
|
132 my $pep = $self->peptide;
|
|
|
133
|
|
|
134 return undef unless $pep;
|
|
|
135
|
|
|
136 my $ref_pep = $self->transcript_variation->get_reference_TranscriptVariationAllele->peptide;
|
|
|
137
|
|
|
138 return undef unless $ref_pep;
|
|
|
139
|
|
|
140 return $ref_pep ne $pep ? $ref_pep.'/'.$pep : $pep;
|
|
|
141 }
|
|
|
142
|
|
|
143 =head2 codon_allele_string
|
|
|
144
|
|
|
145 Description: Return a '/' delimited string of the reference codon and the
|
|
|
146 codon resulting from this allele
|
|
|
147 Returntype : string or undef if this allele is not in the CDS
|
|
|
148 Exceptions : none
|
|
|
149 Status : At Risk
|
|
|
150
|
|
|
151 =cut
|
|
|
152
|
|
|
153 sub codon_allele_string {
|
|
|
154 my ($self) = @_;
|
|
|
155
|
|
|
156 my $codon = $self->codon;
|
|
|
157
|
|
|
158 return undef unless $codon;
|
|
|
159
|
|
|
160 my $ref_codon = $self->transcript_variation->get_reference_TranscriptVariationAllele->codon;
|
|
|
161
|
|
|
162 return $ref_codon.'/'.$codon;
|
|
|
163 }
|
|
|
164
|
|
|
165 =head2 display_codon_allele_string
|
|
|
166
|
|
|
167 Description: Return a '/' delimited string of the reference display_codon and the
|
|
|
168 display_codon resulting from this allele. The display_codon identifies
|
|
|
169 the nucleotides affected by this variant in UPPER CASE and other
|
|
|
170 nucleotides in lower case
|
|
|
171 Returntype : string or undef if this allele is not in the CDS
|
|
|
172 Exceptions : none
|
|
|
173 Status : At Risk
|
|
|
174
|
|
|
175 =cut
|
|
|
176
|
|
|
177 sub display_codon_allele_string {
|
|
|
178 my ($self) = @_;
|
|
|
179
|
|
|
180 my $display_codon = $self->display_codon;
|
|
|
181
|
|
|
182 return undef unless $display_codon;
|
|
|
183
|
|
|
184 my $ref_display_codon = $self->transcript_variation->get_reference_TranscriptVariationAllele->display_codon;
|
|
|
185
|
|
|
186 return $ref_display_codon.'/'.$display_codon;
|
|
|
187 }
|
|
|
188
|
|
|
189 =head2 peptide
|
|
|
190
|
|
|
191 Description: Return the amino acid sequence that this allele is predicted to result in
|
|
|
192 Returntype : string or undef if this allele is not in the CDS or is a frameshift
|
|
|
193 Exceptions : none
|
|
|
194 Status : At Risk
|
|
|
195
|
|
|
196 =cut
|
|
|
197
|
|
|
198 sub peptide {
|
|
|
199 my ($self, $peptide) = @_;
|
|
|
200
|
|
|
201 $self->{peptide} = $peptide if $peptide;
|
|
|
202
|
|
|
203 unless ($self->{peptide}) {
|
|
|
204
|
|
|
205 return undef unless $self->seq_is_unambiguous_dna;
|
|
|
206
|
|
|
207 if (my $codon = $self->codon) {
|
|
|
208
|
|
|
209 # the codon method can set the peptide in some circumstances
|
|
|
210 # so check here before we try an (expensive) translation
|
|
|
211 return $self->{peptide} if $self->{peptide};
|
|
|
212
|
|
|
213 # translate the codon sequence to establish the peptide allele
|
|
|
214
|
|
|
215 # allow for partial codons - split the sequence into whole and partial
|
|
|
216 # e.g. AAAGG split into AAA and GG
|
|
|
217 my $whole_codon = substr($codon, 0, int(length($codon) / 3) * 3);
|
|
|
218 my $partial_codon = substr($codon, int(length($codon) / 3) * 3);
|
|
|
219
|
|
|
220 my $pep = '';
|
|
|
221
|
|
|
222 if($whole_codon) {
|
|
|
223 # for mithocondrial dna we need to to use a different codon table
|
|
|
224 my $codon_table = $self->transcript_variation->_codon_table;
|
|
|
225
|
|
|
226 my $codon_seq = Bio::Seq->new(
|
|
|
227 -seq => $whole_codon,
|
|
|
228 -moltype => 'dna',
|
|
|
229 -alphabet => 'dna',
|
|
|
230 );
|
|
|
231
|
|
|
232 $pep .= $codon_seq->translate(undef, undef, undef, $codon_table)->seq;
|
|
|
233 }
|
|
|
234
|
|
|
235 if($partial_codon) {
|
|
|
236 $pep .= 'X';
|
|
|
237 }
|
|
|
238
|
|
|
239 $pep ||= '-';
|
|
|
240
|
|
|
241 $self->{peptide} = $pep;
|
|
|
242 }
|
|
|
243 }
|
|
|
244
|
|
|
245 return $self->{peptide};
|
|
|
246 }
|
|
|
247
|
|
|
248 =head2 codon
|
|
|
249
|
|
|
250 Description: Return the codon sequence that this allele is predicted to result in
|
|
|
251 Returntype : string or undef if this allele is not in the CDS or is a frameshift
|
|
|
252 Exceptions : none
|
|
|
253 Status : At Risk
|
|
|
254
|
|
|
255 =cut
|
|
|
256
|
|
|
257 sub codon {
|
|
|
258 my ($self, $codon) = @_;
|
|
|
259
|
|
|
260 $self->{codon} = $codon if defined $codon;
|
|
|
261
|
|
|
262 my $tv = $self->transcript_variation;
|
|
|
263
|
|
|
264 return undef unless $tv->translation_start && $tv->translation_end;
|
|
|
265
|
|
|
266 return undef unless $self->seq_is_dna;
|
|
|
267
|
|
|
268 unless ($self->{codon}) {
|
|
|
269
|
|
|
270 # try to calculate the codon sequence
|
|
|
271
|
|
|
272 my $seq = $self->feature_seq;
|
|
|
273
|
|
|
274 $seq = '' if $seq eq '-';
|
|
|
275
|
|
|
276 # calculate necessary coords and lengths
|
|
|
277
|
|
|
278 my $codon_cds_start = $tv->translation_start * 3 - 2;
|
|
|
279 my $codon_cds_end = $tv->translation_end * 3;
|
|
|
280 my $codon_len = $codon_cds_end - $codon_cds_start + 1;
|
|
|
281 my $vf_nt_len = $tv->cds_end - $tv->cds_start + 1;
|
|
|
282 my $allele_len = $self->seq_length;
|
|
|
283
|
|
|
284 if ($allele_len != $vf_nt_len) {
|
|
|
285 if (abs($allele_len - $vf_nt_len) % 3) {
|
|
|
286 # this is a frameshift variation, we don't attempt to
|
|
|
287 # calculate the resulting codon or peptide change as this
|
|
|
288 # could get quite complicated
|
|
|
289 return undef;
|
|
|
290 }
|
|
|
291 }
|
|
|
292
|
|
|
293 # splice the allele sequence into the CDS
|
|
|
294
|
|
|
295 my $cds = $tv->_translateable_seq;
|
|
|
296
|
|
|
297 substr($cds, $tv->cds_start-1, $vf_nt_len) = $seq;
|
|
|
298
|
|
|
299 # and extract the codon sequence
|
|
|
300
|
|
|
301 my $codon = substr($cds, $codon_cds_start-1, $codon_len + ($allele_len - $vf_nt_len));
|
|
|
302
|
|
|
303 if (length($codon) < 1) {
|
|
|
304 $self->{codon} = '-';
|
|
|
305 $self->{peptide} = '-';
|
|
|
306 }
|
|
|
307 else {
|
|
|
308 $self->{codon} = $codon;
|
|
|
309 }
|
|
|
310 }
|
|
|
311
|
|
|
312 return $self->{codon};
|
|
|
313 }
|
|
|
314
|
|
|
315 =head2 display_codon
|
|
|
316
|
|
|
317 Description: Return the codon sequence that this allele is predicted to result in
|
|
|
318 with the affected nucleotides identified in UPPER CASE and other
|
|
|
319 nucleotides in lower case
|
|
|
320 Returntype : string or undef if this allele is not in the CDS or is a frameshift
|
|
|
321 Exceptions : none
|
|
|
322 Status : At Risk
|
|
|
323
|
|
|
324 =cut
|
|
|
325
|
|
|
326 sub display_codon {
|
|
|
327 my $self = shift;
|
|
|
328
|
|
|
329 unless ($self->{_display_codon}) {
|
|
|
330
|
|
|
331 if ($self->codon && defined $self->transcript_variation->codon_position) {
|
|
|
332
|
|
|
333 my $display_codon = lc $self->codon;
|
|
|
334
|
|
|
335 # if this allele is an indel then just return all lowercase
|
|
|
336
|
|
|
337 if ($self->feature_seq ne '-') {
|
|
|
338
|
|
|
339 # codon_position is 1-based, while substr assumes the string starts at 0
|
|
|
340
|
|
|
341 my $pos = $self->transcript_variation->codon_position - 1;
|
|
|
342
|
|
|
343 my $len = length $self->feature_seq;
|
|
|
344
|
|
|
345 substr($display_codon, $pos, $len) = uc substr($display_codon, $pos, $len);
|
|
|
346 }
|
|
|
347
|
|
|
348 $self->{_display_codon} = $display_codon;
|
|
|
349 }
|
|
|
350 }
|
|
|
351
|
|
|
352 return $self->{_display_codon};
|
|
|
353 }
|
|
|
354
|
|
|
355 =head2 polyphen_prediction
|
|
|
356
|
|
|
357 Description: Return the qualitative PolyPhen-2 prediction for the effect of this allele.
|
|
|
358 (Note that we currently only have PolyPhen predictions for variants that
|
|
|
359 result in single amino acid substitutions in human)
|
|
|
360 Returntype : string (one of 'probably damaging', 'possibly damaging', 'benign', 'unknown')
|
|
|
361 if this is a non-synonymous mutation and a prediction is available, undef
|
|
|
362 otherwise
|
|
|
363 Exceptions : none
|
|
|
364 Status : At Risk
|
|
|
365
|
|
|
366 =cut
|
|
|
367
|
|
|
368 sub polyphen_prediction {
|
|
|
369 my ($self, $classifier, $polyphen_prediction) = @_;
|
|
|
370
|
|
|
371 $classifier ||= 'humvar';
|
|
|
372
|
|
|
373 my $analysis = "polyphen_${classifier}";
|
|
|
374
|
|
|
375 $self->{$analysis}->{prediction} = $polyphen_prediction if $polyphen_prediction;
|
|
|
376
|
|
|
377 unless ($self->{$analysis}->{prediction}) {
|
|
|
378 my ($prediction, $score) = $self->_protein_function_prediction($analysis);
|
|
|
379 $self->{$analysis}->{score} = $score;
|
|
|
380 $self->{$analysis}->{prediction} = $prediction;
|
|
|
381 }
|
|
|
382
|
|
|
383 return $self->{$analysis}->{prediction};
|
|
|
384 }
|
|
|
385
|
|
|
386 =head2 polyphen_score
|
|
|
387
|
|
|
388 Description: Return the PolyPhen-2 probability that this allele is deleterious (Note that we
|
|
|
389 currently only have PolyPhen predictions for variants that result in single
|
|
|
390 amino acid substitutions in human)
|
|
|
391 Returntype : float between 0 and 1 if this is a non-synonymous mutation and a prediction is
|
|
|
392 available, undef otherwise
|
|
|
393 Exceptions : none
|
|
|
394 Status : At Risk
|
|
|
395
|
|
|
396 =cut
|
|
|
397
|
|
|
398 sub polyphen_score {
|
|
|
399 my ($self, $classifier, $polyphen_score) = @_;
|
|
|
400
|
|
|
401 $classifier ||= 'humvar';
|
|
|
402
|
|
|
403 my $analysis = "polyphen_${classifier}";
|
|
|
404
|
|
|
405 $self->{$analysis}->{score} = $polyphen_score if defined $polyphen_score;
|
|
|
406
|
|
|
407 unless ($self->{$analysis}->{score}) {
|
|
|
408 my ($prediction, $score) = $self->_protein_function_prediction($analysis);
|
|
|
409 $self->{$analysis}->{score} = $score;
|
|
|
410 $self->{$analysis}->{prediction} = $prediction;
|
|
|
411 }
|
|
|
412
|
|
|
413 return $self->{$analysis}->{score};
|
|
|
414 }
|
|
|
415
|
|
|
416 =head2 sift_prediction
|
|
|
417
|
|
|
418 Description: Return the qualitative SIFT prediction for the effect of this allele.
|
|
|
419 (Note that we currently only have SIFT predictions for variants that
|
|
|
420 result in single amino acid substitutions in human)
|
|
|
421 Returntype : string (one of 'tolerated', 'deleterious') if this is a non-synonymous
|
|
|
422 mutation and a prediction is available, undef otherwise
|
|
|
423 Exceptions : none
|
|
|
424 Status : At Risk
|
|
|
425
|
|
|
426 =cut
|
|
|
427
|
|
|
428 sub sift_prediction {
|
|
|
429 my ($self, $sift_prediction) = @_;
|
|
|
430
|
|
|
431 $self->{sift_prediction} = $sift_prediction if $sift_prediction;
|
|
|
432
|
|
|
433 unless ($self->{sift_prediction}) {
|
|
|
434 my ($prediction, $score) = $self->_protein_function_prediction('sift');
|
|
|
435 $self->{sift_score} = $score;
|
|
|
436 $self->{sift_prediction} = $prediction unless $self->{sift_prediction};
|
|
|
437 }
|
|
|
438
|
|
|
439 return $self->{sift_prediction};
|
|
|
440 }
|
|
|
441
|
|
|
442 =head2 sift_score
|
|
|
443
|
|
|
444 Description: Return the SIFT score for this allele (Note that we currently only have SIFT
|
|
|
445 predictions for variants that result in single amino acid substitutions in human)
|
|
|
446 Returntype : float between 0 and 1 if this is a non-synonymous mutation and a prediction is
|
|
|
447 available, undef otherwise
|
|
|
448 Exceptions : none
|
|
|
449 Status : At Risk
|
|
|
450
|
|
|
451 =cut
|
|
|
452
|
|
|
453 sub sift_score {
|
|
|
454 my ($self, $sift_score) = @_;
|
|
|
455
|
|
|
456 $self->{sift_score} = $sift_score if defined $sift_score;
|
|
|
457
|
|
|
458 unless ($self->{sift_score}) {
|
|
|
459 my ($prediction, $score) = $self->_protein_function_prediction('sift');
|
|
|
460 $self->{sift_score} = $score;
|
|
|
461 $self->{sift_prediction} = $prediction;
|
|
|
462 }
|
|
|
463
|
|
|
464 return $self->{sift_score};
|
|
|
465 }
|
|
|
466
|
|
|
467
|
|
|
468 sub _protein_function_prediction {
|
|
|
469 my ($self, $analysis) = @_;
|
|
|
470
|
|
|
471 # we can only get results for variants that cause a single amino acid substitution,
|
|
|
472 # so check the peptide allele string first
|
|
|
473
|
|
|
474 if ($self->pep_allele_string && $self->pep_allele_string =~ /^[A-Z]\/[A-Z]$/ && defined $AA_LOOKUP->{$self->peptide}) {
|
|
|
475
|
|
|
476 if (my $matrix = $self->transcript_variation->_protein_function_predictions($analysis)) {
|
|
|
477
|
|
|
478 my ($prediction, $score) = $matrix->get_prediction(
|
|
|
479 $self->transcript_variation->translation_start,
|
|
|
480 $self->peptide,
|
|
|
481 );
|
|
|
482
|
|
|
483 return wantarray ? ($prediction, $score) : $prediction;
|
|
|
484 }
|
|
|
485 }
|
|
|
486
|
|
|
487 return undef;
|
|
|
488 }
|
|
|
489
|
|
|
490 =head2 hgvs_genomic
|
|
|
491
|
|
|
492 Description: Return a string representing the genomic-level effect of this allele in HGVS format
|
|
|
493 Returntype : string
|
|
|
494 Exceptions : none
|
|
|
495 Status : At Risk
|
|
|
496
|
|
|
497 =cut
|
|
|
498
|
|
|
499 sub hgvs_genomic {
|
|
|
500 return _hgvs_generic(@_,'genomic');
|
|
|
501 }
|
|
|
502
|
|
|
503 =head2 hgvs_coding
|
|
|
504
|
|
|
505 Description: Return a string representing the CDS-level effect of this allele in HGVS format
|
|
|
506 Returntype : string or undef if this allele is not in the CDS
|
|
|
507 Exceptions : none
|
|
|
508 Status : At Risk
|
|
|
509
|
|
|
510 =cut
|
|
|
511
|
|
|
512 sub hgvs_coding {
|
|
|
513
|
|
|
514 deprecate('HGVS coding support has been moved to hgvs_transcript. This method will be removed in the next release.');
|
|
|
515 return hgvs_transcript(@_);
|
|
|
516 }
|
|
|
517
|
|
|
518
|
|
|
519 =head2 hgvs_transcript
|
|
|
520
|
|
|
521 Description: Return a string representing the CDS-level effect of this allele in HGVS format
|
|
|
522 Returntype : string or undef if this allele is not in the CDS
|
|
|
523 Exceptions : none
|
|
|
524 Status : At Risk
|
|
|
525
|
|
|
526 =cut
|
|
|
527
|
|
|
528
|
|
|
529 sub hgvs_transcript {
|
|
|
530
|
|
|
531
|
|
|
532 my $self = shift;
|
|
|
533 my $notation = shift;
|
|
|
534
|
|
|
535 ##### set if string supplied
|
|
|
536 $self->{hgvs_transcript} = $notation if defined $notation;
|
|
|
537 ##### return if held
|
|
|
538 return $self->{hgvs_transcript} if defined $self->{hgvs_transcript} ;
|
|
|
539 return $self->{hgvs_coding} if defined $self->{hgvs_coding} ;
|
|
|
540
|
|
|
541 my $variation_feature_sequence = $self->variation_feature_seq() ;
|
|
|
542
|
|
|
543 ### don't try to handle odd characters
|
|
|
544 return undef if $variation_feature_sequence =~ m/[^ACGT\-]/ig;
|
|
|
545
|
|
|
546 ### no result for reference allele
|
|
|
547 return undef if $self->is_reference ==1;
|
|
|
548
|
|
|
549 ### else evaluate
|
|
|
550
|
|
|
551 ### get reference sequence strand
|
|
|
552 my $refseq_strand = $self->transcript_variation->transcript->strand();
|
|
|
553
|
|
|
554 if($DEBUG ==1){
|
|
|
555 my $var_name = $self->transcript_variation->variation_feature->variation_name();
|
|
|
556 print "\nHGVS transcript: Checking $var_name refseq strand => $refseq_strand seq name : " . $self->transcript_variation->transcript_stable_id() . " var strand " . $self->transcript_variation->variation_feature->strand() . " vf st " . $self->variation_feature->strand() ." seqname: " . $self->variation_feature->seqname() ." seq: " . $self->variation_feature_seq ."\n";
|
|
|
557 }
|
|
|
558
|
|
|
559 my $hgvs_notation ; ### store components of HGVS string in hash
|
|
|
560
|
|
|
561 ### vf strand is relative to transcript or transcript slice
|
|
|
562 if( $self->transcript_variation->variation_feature->strand() <0 && $refseq_strand >0 ||
|
|
|
563 $self->transcript_variation->variation_feature->strand() >0 && $refseq_strand < 0
|
|
|
564 ){
|
|
|
565 reverse_comp(\$variation_feature_sequence);
|
|
|
566 }
|
|
|
567
|
|
|
568 ### need to get ref seq from slice transcript is on for intron labelling
|
|
|
569 my ($slice_start, $slice_end, $slice) = $self->_var2transcript_slice_coords();
|
|
|
570
|
|
|
571
|
|
|
572 unless($slice_start){
|
|
|
573 #warn "TVA: VF not within transcript - no HGVS\n";
|
|
|
574 return undef;
|
|
|
575 }
|
|
|
576
|
|
|
577 ### decide event type from HGVS nomenclature
|
|
|
578 $hgvs_notation = hgvs_variant_notation( $variation_feature_sequence, ### alt_allele,
|
|
|
579 $slice->seq(), ### using this to extract ref allele
|
|
|
580 $slice_start,
|
|
|
581 $slice_end
|
|
|
582 );
|
|
|
583
|
|
|
584 ### This should not happen
|
|
|
585 unless($hgvs_notation->{'type'}){
|
|
|
586 #warn "Error - not continuing; no HGVS annotation\n";
|
|
|
587 return undef;
|
|
|
588 }
|
|
|
589
|
|
|
590 ### create reference name - transcript name & seq version
|
|
|
591 $hgvs_notation->{'ref_name'} = $self->transcript_variation->transcript_stable_id() . "." . $self->transcript_variation->transcript->version();
|
|
|
592
|
|
|
593
|
|
|
594 ### get position relative to transcript features [use HGVS coords not variation feature coords due to dups]
|
|
|
595 $hgvs_notation->{start} = $self->_get_cDNA_position( $hgvs_notation->{start} );
|
|
|
596 $hgvs_notation->{end} = $self->_get_cDNA_position( $hgvs_notation->{end} );
|
|
|
597
|
|
|
598
|
|
|
599 # Make sure that start is always less than end
|
|
|
600 my ($exon_start_coord, $intron_start_offset) = $hgvs_notation->{start} =~ m/(\-?[0-9]+)\+?(\-?[0-9]+)?/;
|
|
|
601 my ($exon_end_coord, $intron_end_offset) = $hgvs_notation->{end} =~ m/(\-?[0-9]+)\+?(\-?[0-9]+)?/;
|
|
|
602 $intron_start_offset ||= 0;
|
|
|
603 $intron_end_offset ||= 0;
|
|
|
604
|
|
|
605 ($hgvs_notation->{start},$hgvs_notation->{end}) = ($hgvs_notation->{end},$hgvs_notation->{start}) if
|
|
|
606 (($exon_start_coord + $intron_start_offset) > ($exon_end_coord + $intron_end_offset));
|
|
|
607
|
|
|
608
|
|
|
609
|
|
|
610 if($self->transcript->cdna_coding_start()){
|
|
|
611 $hgvs_notation->{'numbering'} = "c"; ### set 'c' if transcript is coding
|
|
|
612 }
|
|
|
613 else{
|
|
|
614 $hgvs_notation->{'numbering'} = "n"; ### set 'n' if transcript non-coding
|
|
|
615 }
|
|
|
616 ### generic formatting
|
|
|
617 $self->{hgvs_transcript} = format_hgvs_string( $hgvs_notation);
|
|
|
618
|
|
|
619 if($DEBUG ==1){print "HGVS notation for var " . $self->transcript_variation->variation_feature->variation_name() . " from hgvs transcript : " . $self->{hgvs_transcript} . " \n";}
|
|
|
620
|
|
|
621 return $self->{hgvs_transcript};
|
|
|
622
|
|
|
623 }
|
|
|
624
|
|
|
625
|
|
|
626 =head2 hgvs_protein
|
|
|
627
|
|
|
628 Description: Return a string representing the protein-level effect of this allele in HGVS format
|
|
|
629 Returntype : string or undef if this allele is not in the CDS
|
|
|
630 Exceptions : none
|
|
|
631 Status : At Risk
|
|
|
632
|
|
|
633 =cut
|
|
|
634
|
|
|
635 sub hgvs_protein {
|
|
|
636
|
|
|
637 my $self = shift;
|
|
|
638 my $notation = shift;
|
|
|
639 my $hgvs_notation;
|
|
|
640
|
|
|
641 if($DEBUG ==1){print "\nStarting hgvs_protein with ". $self->transcript_variation->variation_feature->variation_name() . "\n"; }
|
|
|
642
|
|
|
643 ### set if string supplied
|
|
|
644 $self->{hgvs_protein} = $notation if defined $notation;
|
|
|
645
|
|
|
646 ### return if set
|
|
|
647 return $self->{hgvs_protein} if defined $self->{hgvs_protein} ;
|
|
|
648
|
|
|
649 ### don't try to handle odd characters
|
|
|
650 return undef if $self->variation_feature_seq() =~ m/[^ACGT\-]/ig;
|
|
|
651
|
|
|
652 #### else, check viable input and create notation
|
|
|
653 return if $self->is_reference();
|
|
|
654
|
|
|
655 unless ($self->transcript_variation->translation_start()){
|
|
|
656 return undef ;
|
|
|
657 }
|
|
|
658
|
|
|
659 #### for debug
|
|
|
660 #my $var_name = $self->transcript_variation->variation_feature->variation_name();
|
|
|
661
|
|
|
662 ### get reference sequence [add seq version to transcript name]
|
|
|
663 $hgvs_notation->{ref_name} = $self->transcript_variation->transcript->translation()->display_id() . "." . $self->transcript_variation->transcript->translation()->version();
|
|
|
664
|
|
|
665 $hgvs_notation->{'numbering'} = 'p';
|
|
|
666
|
|
|
667 ### get default reference location [changed later in some cases eg. duplication]
|
|
|
668 $hgvs_notation->{start} = $self->transcript_variation->translation_start();
|
|
|
669 $hgvs_notation->{end} = $self->transcript_variation->translation_end();
|
|
|
670
|
|
|
671
|
|
|
672 ## get default reference & alt peptides [changed later to hgvs format]
|
|
|
673 $hgvs_notation->{alt} = $self->peptide;
|
|
|
674 $hgvs_notation->{ref} = $self->transcript_variation->get_reference_TranscriptVariationAllele->peptide;
|
|
|
675
|
|
|
676 if(defined $hgvs_notation->{alt} && defined $hgvs_notation->{ref}){
|
|
|
677 $hgvs_notation = _clip_alleles( $hgvs_notation);
|
|
|
678 }
|
|
|
679
|
|
|
680 #### define type - types are different for protein numbering
|
|
|
681 $hgvs_notation = $self->_get_hgvs_protein_type($hgvs_notation);
|
|
|
682
|
|
|
683 ##### Convert ref & alt peptides taking into account HGVS rules
|
|
|
684 $hgvs_notation = $self->_get_hgvs_peptides($hgvs_notation);
|
|
|
685
|
|
|
686 ##### String formatting
|
|
|
687 $self->{hgvs_protein} = $self->_get_hgvs_protein_format($hgvs_notation);
|
|
|
688 return $self->{hgvs_protein} ;
|
|
|
689 }
|
|
|
690
|
|
|
691 ### HGVS: format protein string
|
|
|
692 sub _get_hgvs_protein_format{
|
|
|
693
|
|
|
694 my $self = shift;
|
|
|
695 my $hgvs_notation = shift;
|
|
|
696
|
|
|
697 if ((defined $hgvs_notation->{ref} && defined $hgvs_notation->{alt} &&
|
|
|
698 $hgvs_notation->{ref} eq $hgvs_notation->{alt}) ||
|
|
|
699 (defined $hgvs_notation->{type} && $hgvs_notation->{type} eq "=")){
|
|
|
700
|
|
|
701 ### no protein change - return transcript nomenclature with flag for neutral protein consequence
|
|
|
702 $hgvs_notation->{'hgvs'} = $self->hgvs_transcript() . "(p.=)";
|
|
|
703 return $hgvs_notation->{'hgvs'} ;
|
|
|
704 }
|
|
|
705
|
|
|
706 ### all start with refseq name & numbering type
|
|
|
707 $hgvs_notation->{'hgvs'} = $hgvs_notation->{'ref_name'} . ":" . $hgvs_notation->{'numbering'} . ".";
|
|
|
708
|
|
|
709 ### handle stop_lost seperately regardless of cause by del/delins => p.XposAA1extnum_AA_to_stop
|
|
|
710 if(stop_lost($self)){ ### if deletion of stop add extX and number of new aa to alt
|
|
|
711 $hgvs_notation->{alt} = substr($hgvs_notation->{alt}, 0, 3);
|
|
|
712 if($hgvs_notation->{type} eq "del"){
|
|
|
713 my $aa_til_stop = _stop_loss_extra_AA($self,$hgvs_notation->{start}-1, "del");
|
|
|
714 if(defined $aa_til_stop){
|
|
|
715 $hgvs_notation->{alt} .= "extX" . $aa_til_stop;
|
|
|
716 }
|
|
|
717 }
|
|
|
718 elsif($hgvs_notation->{type} eq ">"){
|
|
|
719 my $aa_til_stop = _stop_loss_extra_AA($self,$hgvs_notation->{start}-1, "subs");
|
|
|
720 if(defined $aa_til_stop){
|
|
|
721 $hgvs_notation->{alt} .= "extX" . $aa_til_stop;
|
|
|
722 }
|
|
|
723 }
|
|
|
724 else{
|
|
|
725 # warn "TVA: stop loss for type $hgvs_notation->{type} not caught \n";
|
|
|
726 }
|
|
|
727 $hgvs_notation->{'hgvs'} .= $hgvs_notation->{ref} . $hgvs_notation->{start} . $hgvs_notation->{alt} ;
|
|
|
728 }
|
|
|
729
|
|
|
730 elsif( $hgvs_notation->{type} eq "dup"){
|
|
|
731
|
|
|
732 if($hgvs_notation->{start} < $hgvs_notation->{end}){
|
|
|
733 ### list only first and last peptides in long duplicated string
|
|
|
734 my $ref_pep_first = substr($hgvs_notation->{alt}, 0, 3);
|
|
|
735 my $ref_pep_last = substr($hgvs_notation->{alt}, -3, 3);
|
|
|
736 $hgvs_notation->{'hgvs'} .= $ref_pep_first . $hgvs_notation->{start} . "_" . $ref_pep_last . $hgvs_notation->{end} ."dup";
|
|
|
737
|
|
|
738 }
|
|
|
739 else{
|
|
|
740 $hgvs_notation->{'hgvs'} .= $hgvs_notation->{alt} . $hgvs_notation->{start} . "dup" ;
|
|
|
741 }
|
|
|
742 print "formating a dup $hgvs_notation->{hgvs} \n";
|
|
|
743 }
|
|
|
744
|
|
|
745 elsif($hgvs_notation->{type} eq ">"){
|
|
|
746 #### substitution
|
|
|
747
|
|
|
748 $hgvs_notation->{'hgvs'} .= $hgvs_notation->{ref}. $hgvs_notation->{start} . $hgvs_notation->{alt};
|
|
|
749 }
|
|
|
750
|
|
|
751 elsif( $hgvs_notation->{type} eq "delins" || $hgvs_notation->{type} eq "ins" ){
|
|
|
752
|
|
|
753 #### list first and last aa in reference only
|
|
|
754 my $ref_pep_first = substr($hgvs_notation->{ref}, 0, 3);
|
|
|
755 my $ref_pep_last;
|
|
|
756 if(substr($hgvs_notation->{ref}, -1, 1) eq "X"){
|
|
|
757 $ref_pep_last ="X";
|
|
|
758 }
|
|
|
759 else{
|
|
|
760 $ref_pep_last = substr($hgvs_notation->{ref}, -3, 3);
|
|
|
761 }
|
|
|
762 if($hgvs_notation->{ref} =~ /X$/){
|
|
|
763 ### For stops & add extX & distance to next stop to alt pep
|
|
|
764 my $aa_til_stop = _stop_loss_extra_AA($self,$hgvs_notation->{start}-1, "loss");
|
|
|
765 if(defined $aa_til_stop){
|
|
|
766 $hgvs_notation->{alt} .="extX" . $aa_til_stop;
|
|
|
767 }
|
|
|
768 }
|
|
|
769 if($hgvs_notation->{start} == $hgvs_notation->{end} && $hgvs_notation->{type} eq "delins"){
|
|
|
770 $hgvs_notation->{'hgvs'} .= $ref_pep_first . $hgvs_notation->{start} . $hgvs_notation->{end} . $hgvs_notation->{type} . $hgvs_notation->{alt} ;
|
|
|
771 }
|
|
|
772 else{
|
|
|
773 ### correct ordering if needed
|
|
|
774 if($hgvs_notation->{start} > $hgvs_notation->{end}){
|
|
|
775 ( $hgvs_notation->{start}, $hgvs_notation->{end}) = ($hgvs_notation->{end}, $hgvs_notation->{start} );
|
|
|
776 }
|
|
|
777
|
|
|
778 $hgvs_notation->{'hgvs'} .= $ref_pep_first . $hgvs_notation->{start} . "_" . $ref_pep_last . $hgvs_notation->{end} . $hgvs_notation->{type} . $hgvs_notation->{alt} ;
|
|
|
779 }
|
|
|
780 }
|
|
|
781
|
|
|
782 elsif($hgvs_notation->{type} eq "fs"){
|
|
|
783 if(defined $hgvs_notation->{alt} && $hgvs_notation->{alt} eq "X"){ ## stop gained
|
|
|
784 $hgvs_notation->{'hgvs'} .= $hgvs_notation->{ref} . $hgvs_notation->{start} . $hgvs_notation->{alt} ;
|
|
|
785
|
|
|
786 }
|
|
|
787 else{ ## not immediate stop - count aa until next
|
|
|
788
|
|
|
789 my $aa_til_stop = _stop_loss_extra_AA($self, $hgvs_notation->{start}-1, "fs");
|
|
|
790 if($aa_til_stop ){
|
|
|
791 ### use long form if new stop found
|
|
|
792 $hgvs_notation->{'hgvs'} .= $hgvs_notation->{ref} . $hgvs_notation->{start} . $hgvs_notation->{alt} . "fsX$aa_til_stop" ;
|
|
|
793 }
|
|
|
794 else{
|
|
|
795 ### otherwise use short form
|
|
|
796 $hgvs_notation->{'hgvs'} .= $hgvs_notation->{ref} . $hgvs_notation->{start} . "fs";
|
|
|
797 }
|
|
|
798 }
|
|
|
799 }
|
|
|
800
|
|
|
801 elsif( $hgvs_notation->{type} eq "del"){
|
|
|
802 $hgvs_notation->{alt} = "del";
|
|
|
803 if( length($hgvs_notation->{ref}) >3 ){
|
|
|
804 my $ref_pep_first = substr($hgvs_notation->{ref}, 0, 3);
|
|
|
805 my $ref_pep_last = substr($hgvs_notation->{ref}, -3, 3);
|
|
|
806 $hgvs_notation->{'hgvs'} .= $ref_pep_first . $hgvs_notation->{start} . "_" . $ref_pep_last . $hgvs_notation->{end} .$hgvs_notation->{alt} ;
|
|
|
807 }
|
|
|
808 else{
|
|
|
809 $hgvs_notation->{'hgvs'} .= $hgvs_notation->{ref} . $hgvs_notation->{start} . $hgvs_notation->{alt} ;
|
|
|
810 }
|
|
|
811 }
|
|
|
812
|
|
|
813 elsif($hgvs_notation->{start} ne $hgvs_notation->{end} ){
|
|
|
814 $hgvs_notation->{'hgvs'} .= $hgvs_notation->{ref} . $hgvs_notation->{start} . "_" . $hgvs_notation->{alt} . $hgvs_notation->{end} ;
|
|
|
815 }
|
|
|
816
|
|
|
817 else{
|
|
|
818 #### default to substitution
|
|
|
819 $hgvs_notation->{'hgvs'} .= $hgvs_notation->{ref}. $hgvs_notation->{start} . $hgvs_notation->{alt};
|
|
|
820 }
|
|
|
821
|
|
|
822 if($DEBUG==1){ print "Returning protein format: $hgvs_notation->{'hgvs'}\n";}
|
|
|
823 return $hgvs_notation->{'hgvs'};
|
|
|
824 }
|
|
|
825
|
|
|
826 ### HGVS: get type of variation event in protein terms
|
|
|
827 sub _get_hgvs_protein_type{
|
|
|
828
|
|
|
829 my $self = shift;
|
|
|
830 my $hgvs_notation = shift;
|
|
|
831
|
|
|
832
|
|
|
833 ### get allele length
|
|
|
834 my ($ref_length, $alt_length ) = $self->_get_allele_length();
|
|
|
835
|
|
|
836
|
|
|
837 if( defined $hgvs_notation->{ref} && defined $hgvs_notation->{alt} ){
|
|
|
838 ### Run type checks on peptides if available
|
|
|
839
|
|
|
840 if($hgvs_notation->{alt} eq $hgvs_notation->{ref} ){
|
|
|
841 #### synonymous indicated by =
|
|
|
842 $hgvs_notation->{type} = "=";
|
|
|
843 }
|
|
|
844 elsif( length($hgvs_notation->{ref}) ==1 && length($hgvs_notation->{alt}) ==1 ) {
|
|
|
845
|
|
|
846 $hgvs_notation->{type} = ">";
|
|
|
847 }
|
|
|
848 elsif($hgvs_notation->{ref} eq "-" || $hgvs_notation->{ref} eq "") {
|
|
|
849
|
|
|
850 $hgvs_notation->{type} = "ins";
|
|
|
851 }
|
|
|
852 elsif($hgvs_notation->{alt} eq "" ) {
|
|
|
853
|
|
|
854 $hgvs_notation->{type} = "del";
|
|
|
855 }
|
|
|
856 elsif( (length($hgvs_notation->{alt}) > length($hgvs_notation->{ref}) ) &&
|
|
|
857 $hgvs_notation->{alt} =~ / $hgvs_notation->{ref}/ ) {
|
|
|
858 ### capture duplication event described as TTT/TTTTT
|
|
|
859 $hgvs_notation->{type} = "dup";
|
|
|
860 }
|
|
|
861
|
|
|
862 elsif( (length($hgvs_notation->{alt}) >1 && length($hgvs_notation->{ref}) ==1) ||
|
|
|
863 (length($hgvs_notation->{alt}) ==1 && length($hgvs_notation->{ref}) >1) ) {
|
|
|
864 $hgvs_notation->{type} = "delins";
|
|
|
865 }
|
|
|
866 else{
|
|
|
867 $hgvs_notation->{type} = ">";
|
|
|
868 }
|
|
|
869 }
|
|
|
870
|
|
|
871
|
|
|
872 elsif($ref_length ne $alt_length && ($ref_length - $alt_length)%3 !=0 ){
|
|
|
873 $hgvs_notation->{type} = "fs";
|
|
|
874 }
|
|
|
875
|
|
|
876 elsif(length($self->variation_feature_seq()) >1 ){
|
|
|
877 if($hgvs_notation->{start} == ($hgvs_notation->{end} + 1) ){
|
|
|
878 ### convention for insertions - end one less than start
|
|
|
879 $hgvs_notation->{type} = "ins";
|
|
|
880 }
|
|
|
881 elsif( $hgvs_notation->{start} != $hgvs_notation->{end} ){
|
|
|
882 $hgvs_notation->{type} = "delins";
|
|
|
883 }
|
|
|
884 else{
|
|
|
885 $hgvs_notation->{type} = ">";
|
|
|
886 }
|
|
|
887 }
|
|
|
888 else{
|
|
|
889 # warn "Cannot define protein variant type [$ref_length - $alt_length]\n";
|
|
|
890 }
|
|
|
891 return $hgvs_notation ;
|
|
|
892
|
|
|
893 }
|
|
|
894
|
|
|
895 ### HGVS: get reference & alternative peptide
|
|
|
896 sub _get_hgvs_peptides{
|
|
|
897
|
|
|
898 my $self = shift;
|
|
|
899 my $hgvs_notation = shift;
|
|
|
900
|
|
|
901 if($hgvs_notation->{type} eq "fs"){
|
|
|
902 ### ensembl alt/ref peptides not the same as HGVS alt/ref - look up seperately
|
|
|
903 $hgvs_notation = $self->_get_fs_peptides($hgvs_notation);
|
|
|
904
|
|
|
905 }
|
|
|
906 elsif($hgvs_notation->{type} eq "ins" ){
|
|
|
907
|
|
|
908 ### HGVS ref are peptides flanking insertion
|
|
|
909 $hgvs_notation->{ref} = $self->_get_surrounding_peptides($hgvs_notation->{start});
|
|
|
910
|
|
|
911 if( $hgvs_notation->{alt} =~/\*/){
|
|
|
912 ## inserted bases after stop irrelevant; consider as substitution gaining stop MAINTAIN PREVIOUS BEHAVIOUR FOR NOW
|
|
|
913 #$hgvs_notation->{alt} =~ s/\*\w+$/\*/;
|
|
|
914 }
|
|
|
915 else{
|
|
|
916 ### Additional check that inserted bases do not duplicate 3' reference sequence [set to type = dup if so]
|
|
|
917 $hgvs_notation = $self->_check_for_peptide_duplication($hgvs_notation);
|
|
|
918 }
|
|
|
919 }
|
|
|
920
|
|
|
921
|
|
|
922 ### Convert peptide to 3 letter code as used in HGVS
|
|
|
923 $hgvs_notation->{ref} = Bio::SeqUtils->seq3(Bio::PrimarySeq->new(-seq => $hgvs_notation->{ref}, -id => 'ref', -alphabet => 'protein')) || "";
|
|
|
924 if( $hgvs_notation->{alt} eq "-"){
|
|
|
925 $hgvs_notation->{alt} = "del";
|
|
|
926 }
|
|
|
927 else{
|
|
|
928 $hgvs_notation->{alt} = Bio::SeqUtils->seq3(Bio::PrimarySeq->new(-seq => $hgvs_notation->{alt}, -id => 'ref', -alphabet => 'protein')) || "";
|
|
|
929 }
|
|
|
930
|
|
|
931 ### handle special cases
|
|
|
932 if( affects_start_codon($self) ){
|
|
|
933 #### handle initiator loss -> probably no translation => alt allele is '?'
|
|
|
934 $hgvs_notation->{alt} = "?";
|
|
|
935 $hgvs_notation->{type} = "";
|
|
|
936 }
|
|
|
937
|
|
|
938 elsif( $hgvs_notation->{type} eq "del"){
|
|
|
939 #### partial last codon
|
|
|
940 $hgvs_notation->{alt} = "del";
|
|
|
941 }
|
|
|
942 elsif($hgvs_notation->{type} eq "fs"){
|
|
|
943 ### only quote first ref peptide for frameshift
|
|
|
944 $hgvs_notation->{ref} = substr($hgvs_notation->{ref},0,3);
|
|
|
945 }
|
|
|
946
|
|
|
947 ### set stop to HGVS prefered "X"
|
|
|
948 if(defined $hgvs_notation->{ref}){ $hgvs_notation->{ref} =~ s/Ter|Xaa/X/g; }
|
|
|
949 if(defined $hgvs_notation->{alt}){ $hgvs_notation->{alt} =~ s/Ter|Xaa/X/g; }
|
|
|
950
|
|
|
951 return ($hgvs_notation);
|
|
|
952
|
|
|
953 }
|
|
|
954
|
|
|
955 ### HGVS: remove common peptides from alt and ref strings & shift start/end accordingly
|
|
|
956 sub _clip_alleles{
|
|
|
957
|
|
|
958 my $hgvs_notation = shift;
|
|
|
959
|
|
|
960 my $check_alt = $hgvs_notation->{alt} ;
|
|
|
961 my $check_ref = $hgvs_notation->{ref} ;
|
|
|
962 my $check_start = $hgvs_notation->{start};
|
|
|
963
|
|
|
964 ### strip same bases from start of string
|
|
|
965 for (my $p =0; $p <length ($hgvs_notation->{ref}); $p++){
|
|
|
966 my $check_next_ref = substr( $check_ref, 0, 1);
|
|
|
967 my $check_next_alt = substr( $check_alt, 0, 1);
|
|
|
968
|
|
|
969 if($check_next_ref eq "*" && $check_next_alt eq "*"){
|
|
|
970 ### stop re-created by variant - no protein change
|
|
|
971 $hgvs_notation->{type} = "=";
|
|
|
972
|
|
|
973 return($hgvs_notation);
|
|
|
974 }
|
|
|
975
|
|
|
976 if($check_next_ref eq $check_next_alt){
|
|
|
977 $check_start++;
|
|
|
978 $check_ref = substr( $check_ref, 1);
|
|
|
979 $check_alt = substr( $check_alt, 1);
|
|
|
980 }
|
|
|
981 else{
|
|
|
982 last;
|
|
|
983 }
|
|
|
984 }
|
|
|
985 #### strip same bases from end of string
|
|
|
986 for (my $q =0; $q < length ($check_ref); $q++){
|
|
|
987 my $check_next_ref = substr( $check_ref, -1, 1);
|
|
|
988 my $check_next_alt = substr( $check_alt, -1, 1);
|
|
|
989 if($check_next_ref eq $check_next_alt){
|
|
|
990 chop $check_ref;
|
|
|
991 chop $check_alt;
|
|
|
992 }
|
|
|
993 else{
|
|
|
994 last;
|
|
|
995 }
|
|
|
996 }
|
|
|
997
|
|
|
998 $hgvs_notation->{alt} = $check_alt;
|
|
|
999 $hgvs_notation->{ref} = $check_ref;
|
|
|
1000
|
|
|
1001 ### check if clipping force type change & adjust location
|
|
|
1002 if(length ($check_ref) == 0 && length ($check_alt) >= 1){
|
|
|
1003 ### re-set as ins not delins
|
|
|
1004 $hgvs_notation->{type} ="ins";
|
|
|
1005 ### insertion between ref base and next => adjust next
|
|
|
1006 if($hgvs_notation->{end} == $hgvs_notation->{start} ){$hgvs_notation->{end} = $check_start; }
|
|
|
1007 # $hgvs_notation->{start} = $check_start;
|
|
|
1008 }
|
|
|
1009 elsif(length ($check_ref) >=1 && length ($check_alt) == 0){
|
|
|
1010 ### re-set as del not delins
|
|
|
1011 $hgvs_notation->{type} ="del";
|
|
|
1012 $hgvs_notation->{start} = $check_start;
|
|
|
1013 }
|
|
|
1014 else{
|
|
|
1015 #### save trimmed peptide string & increased start position
|
|
|
1016 $hgvs_notation->{start} = $check_start;
|
|
|
1017
|
|
|
1018 }
|
|
|
1019
|
|
|
1020 return $hgvs_notation;
|
|
|
1021 }
|
|
|
1022
|
|
|
1023
|
|
|
1024
|
|
|
1025
|
|
|
1026
|
|
|
1027 #### HGVS: check allele lengths to look for frameshifts
|
|
|
1028 sub _get_allele_length{
|
|
|
1029
|
|
|
1030 my $self = shift;
|
|
|
1031 my $ref_length = 0;
|
|
|
1032 my $alt_length = 0;
|
|
|
1033
|
|
|
1034 my $al_string = $self->allele_string();
|
|
|
1035 my $ref_allele = (split/\//, $al_string)[0];
|
|
|
1036 $ref_allele =~ s/\-//;
|
|
|
1037 $ref_length = length $ref_allele;
|
|
|
1038
|
|
|
1039 my $alt_allele = $self->variation_feature_seq();
|
|
|
1040 $alt_allele =~ s/\-//;
|
|
|
1041 $alt_length = length $alt_allele;
|
|
|
1042
|
|
|
1043 return ($ref_length, $alt_length );
|
|
|
1044
|
|
|
1045 }
|
|
|
1046
|
|
|
1047 ### HGVS: list first different peptide [may not be first changed codon]
|
|
|
1048 sub _get_fs_peptides{
|
|
|
1049
|
|
|
1050 my $self = shift;
|
|
|
1051 my $hgvs_notation = shift;
|
|
|
1052
|
|
|
1053 ### get CDS with alt variant
|
|
|
1054 my $alt_cds = $self->_get_alternate_cds();
|
|
|
1055
|
|
|
1056 #### get new translation
|
|
|
1057 my $alt_trans = $alt_cds->translate()->seq();
|
|
|
1058
|
|
|
1059 ### get changed end (currently in single letter AA coding)
|
|
|
1060 my $ref_trans = $self->transcript->translate()->seq();
|
|
|
1061 $ref_trans .= "*"; ## appending ref stop for checking purposes
|
|
|
1062
|
|
|
1063 $hgvs_notation->{start} = $self->transcript_variation->translation_start() ;
|
|
|
1064
|
|
|
1065 if( $hgvs_notation->{start} > length $alt_trans){ ## deletion of stop, no further AA in alt seq
|
|
|
1066 $hgvs_notation->{alt} = "del";
|
|
|
1067 $hgvs_notation->{type} = "del";
|
|
|
1068 return $hgvs_notation;
|
|
|
1069 }
|
|
|
1070
|
|
|
1071 while ($hgvs_notation->{start} <= length $alt_trans){
|
|
|
1072 ### frame shift may result in the same AA#
|
|
|
1073
|
|
|
1074 $hgvs_notation->{ref} = substr($ref_trans, $hgvs_notation->{start}-1, 1);
|
|
|
1075 $hgvs_notation->{alt} = substr($alt_trans, $hgvs_notation->{start}-1, 1);
|
|
|
1076
|
|
|
1077 if($hgvs_notation->{ref} eq "*" && $hgvs_notation->{alt} eq "*"){
|
|
|
1078 ### variation at stop codon, but maintains stop codon
|
|
|
1079 return ($hgvs_notation);
|
|
|
1080 }
|
|
|
1081 last if $hgvs_notation->{ref} ne $hgvs_notation->{alt};
|
|
|
1082 $hgvs_notation->{start}++;
|
|
|
1083 }
|
|
|
1084 return ($hgvs_notation);
|
|
|
1085
|
|
|
1086 }
|
|
|
1087
|
|
|
1088 #### HGVS: if variant is an insertion, ref pep is initially "-", so seek AA before and after insertion
|
|
|
1089 sub _get_surrounding_peptides{
|
|
|
1090
|
|
|
1091 my $self = shift;
|
|
|
1092 my $ref_pos = shift;
|
|
|
1093 my $ref_trans = $self->transcript->translate()->seq();
|
|
|
1094
|
|
|
1095 my $end = substr($ref_trans, $ref_pos-1);
|
|
|
1096 my $ref_string = substr($ref_trans, $ref_pos-1, 2);
|
|
|
1097
|
|
|
1098 return ($ref_string);
|
|
|
1099
|
|
|
1100 }
|
|
|
1101
|
|
|
1102
|
|
|
1103 #### HGVS: alternate CDS needed to check for new stop when variant disrupts 'reference' stop
|
|
|
1104 sub _get_alternate_cds{
|
|
|
1105
|
|
|
1106 my $self = shift;
|
|
|
1107
|
|
|
1108 ### get reference sequence
|
|
|
1109 my $reference_cds_seq = $self->transcript->translateable_seq();
|
|
|
1110
|
|
|
1111 ### get sequences upstream and downstream of variant
|
|
|
1112 my $upstream_seq = substr($reference_cds_seq, 0, ($self->transcript_variation->cds_start() -1) );
|
|
|
1113 my $downstream_seq = substr($reference_cds_seq, ($self->transcript_variation->cds_end() ) );
|
|
|
1114
|
|
|
1115 ### fix alternate allele if deletion
|
|
|
1116 my $alt_allele = $self->variation_feature_seq();
|
|
|
1117 $alt_allele =~ s/\-//;
|
|
|
1118
|
|
|
1119 ### build alternate seq
|
|
|
1120 my $alternate_seq = $upstream_seq . $alt_allele . $downstream_seq ;
|
|
|
1121
|
|
|
1122 ### create seq obj with alternative allele in the CDS sequence
|
|
|
1123 my $alt_cds =Bio::PrimarySeq->new(-seq => $alternate_seq, -id => 'alt_cds', -alphabet => 'dna');
|
|
|
1124
|
|
|
1125 ### append UTR if available as stop may be disrupted
|
|
|
1126 my $utr = $self->transcript_variation->transcript->three_prime_utr();
|
|
|
1127
|
|
|
1128 if (defined $utr) {
|
|
|
1129 ### append the UTR to the alternative CDS
|
|
|
1130 $alt_cds->seq($alt_cds->seq() . $utr->seq());
|
|
|
1131 }
|
|
|
1132 else{
|
|
|
1133 ##warn "No UTR available for alternate CDS\n";
|
|
|
1134 }
|
|
|
1135
|
|
|
1136 return $alt_cds;
|
|
|
1137 }
|
|
|
1138
|
|
|
1139 ### HGVS: if inserted string is identical to 3' reference sequence, describe as duplication
|
|
|
1140 sub _check_for_peptide_duplication{
|
|
|
1141
|
|
|
1142 my $self = shift;
|
|
|
1143 my $hgvs_notation= shift;
|
|
|
1144
|
|
|
1145 ##### get reference sequence
|
|
|
1146 my $reference_cds_seq = $self->transcript->translateable_seq();
|
|
|
1147
|
|
|
1148 ##### get sequence upstream of variant
|
|
|
1149 my $upstream_seq = substr($reference_cds_seq, 0, ($self->transcript_variation->cds_start() -1) );
|
|
|
1150
|
|
|
1151 ##### create translation to check against inserted peptides
|
|
|
1152 my $upstream_cds =Bio::PrimarySeq->new(-seq => $upstream_seq, -id => 'alt_cds', -alphabet => 'dna');
|
|
|
1153 my $upstream_trans = $upstream_cds->translate()->seq();
|
|
|
1154
|
|
|
1155 ## Test whether alt peptide matches the reference sequence just before the variant
|
|
|
1156 my $test_new_start = $hgvs_notation->{'start'} - length($hgvs_notation->{'alt'}) -1 ;
|
|
|
1157 my $test_seq = substr($upstream_trans, $test_new_start, length($hgvs_notation->{'alt'}));
|
|
|
1158
|
|
|
1159 if ( $test_new_start >= 0 && $test_seq eq $hgvs_notation->{alt}) {
|
|
|
1160
|
|
|
1161 $hgvs_notation->{type} = 'dup';
|
|
|
1162 $hgvs_notation->{end} = $hgvs_notation->{start} -1;
|
|
|
1163 $hgvs_notation->{start} -= length($hgvs_notation->{alt});
|
|
|
1164
|
|
|
1165 }
|
|
|
1166
|
|
|
1167 return $hgvs_notation;
|
|
|
1168
|
|
|
1169 }
|
|
|
1170
|
|
|
1171 #### HGVS: if a stop is lost, seek the next in transcript & count number of extra AA's
|
|
|
1172 sub _stop_loss_extra_AA{
|
|
|
1173
|
|
|
1174 my $self = shift;
|
|
|
1175 my $ref_var_pos = shift; ### first effected AA - supply for frameshifts
|
|
|
1176 my $test = shift;
|
|
|
1177
|
|
|
1178 return undef unless $ref_var_pos;
|
|
|
1179
|
|
|
1180 my $extra_aa;
|
|
|
1181
|
|
|
1182 ### get the sequence with variant added
|
|
|
1183 my $alt_cds = $self->_get_alternate_cds();
|
|
|
1184 ### get new translation
|
|
|
1185 my $alt_trans = $alt_cds->translate();
|
|
|
1186
|
|
|
1187 my $ref_temp = $self->transcript->translate()->seq;
|
|
|
1188 my $ref_len = length($ref_temp);
|
|
|
1189
|
|
|
1190 if($DEBUG==1){
|
|
|
1191 print "alt translated:\n" . $alt_trans->seq() . "\n";
|
|
|
1192 print "ref translated:\n$ref_temp\n";;
|
|
|
1193 }
|
|
|
1194
|
|
|
1195 #### Find the number of residues that are translated until a termination codon is encountered
|
|
|
1196 if ($alt_trans->seq() =~ m/\*/) {
|
|
|
1197 if($DEBUG==1){print "Got $+[0] aa before stop, var event at $ref_var_pos \n";}
|
|
|
1198
|
|
|
1199 if($test eq "fs" ){
|
|
|
1200 ### frame shift - count from first AA effected by variant to stop
|
|
|
1201 $extra_aa = $+[0] - $ref_var_pos;
|
|
|
1202 if($DEBUG==1){ print "Stop change ($test): found $extra_aa amino acids before fs stop [ $+[0] - peptide ref_start: $ref_var_pos )]\n";}
|
|
|
1203 }
|
|
|
1204
|
|
|
1205 else{
|
|
|
1206 $extra_aa = $+[0] - 1 - $ref_len;
|
|
|
1207 if($DEBUG==1){ print "Stop change ($test): found $extra_aa amino acids before next stop [ $+[0] - 1 -normal stop $ref_len)]\n";}
|
|
|
1208 }
|
|
|
1209 }
|
|
|
1210
|
|
|
1211 # A special case is if the first aa is a stop codon => don't display the number of residues until the stop codon
|
|
|
1212 if(defined $extra_aa && $extra_aa >0){
|
|
|
1213 return $extra_aa ;
|
|
|
1214 }
|
|
|
1215 else{
|
|
|
1216 #warn "No stop found in alternate sequence\n";
|
|
|
1217 return undef;
|
|
|
1218 }
|
|
|
1219
|
|
|
1220 }
|
|
|
1221
|
|
|
1222 =head
|
|
|
1223 # We haven't implemented support for these methods yet
|
|
|
1224
|
|
|
1225 sub hgvs_rna {
|
|
|
1226 return _hgvs_generic(@_,'rna');
|
|
|
1227 }
|
|
|
1228
|
|
|
1229 sub hgvs_mitochondrial {
|
|
|
1230 return _hgvs_generic(@_,'mitochondrial');
|
|
|
1231 }
|
|
|
1232
|
|
|
1233 =cut
|
|
|
1234
|
|
|
1235 sub _hgvs_generic {
|
|
|
1236
|
|
|
1237 my $self = shift;
|
|
|
1238 my $reference = pop;
|
|
|
1239 my $notation = shift;
|
|
|
1240
|
|
|
1241 #ÊThe rna and mitochondrial modes have not yet been implemented, so return undef in case we get a call to these
|
|
|
1242 return undef if ($reference =~ m/rna|mitochondrial/);
|
|
|
1243
|
|
|
1244 my $sub = qq{hgvs_$reference};
|
|
|
1245
|
|
|
1246 $self->{$sub} = $notation if defined $notation;
|
|
|
1247
|
|
|
1248 unless ($self->{$sub}) {
|
|
|
1249
|
|
|
1250 # Use the transcript this VF is on as the reference feature
|
|
|
1251 my $reference_feature = $self->transcript;
|
|
|
1252 # If we want genomic coordinates, the reference_feature should actually be the slice for the underlying seq_region
|
|
|
1253 $reference_feature = $reference_feature->slice->seq_region_Slice if ($reference eq 'genomic');
|
|
|
1254
|
|
|
1255 # Calculate the HGVS notation on-the-fly and pass it to the TranscriptVariation in order to distribute the result to the other alleles
|
|
|
1256 $self->transcript_variation->$sub($self->variation_feature->get_all_hgvs_notations($reference_feature,substr($reference,0,1),undef,undef,$self->transcript_variation));
|
|
|
1257 }
|
|
|
1258
|
|
|
1259 return $self->{$sub};
|
|
|
1260 }
|
|
|
1261
|
|
|
1262
|
|
|
1263 ### HGVS: move variant to transcript slice
|
|
|
1264 sub _var2transcript_slice_coords{
|
|
|
1265
|
|
|
1266 my $self = shift;
|
|
|
1267
|
|
|
1268 my $ref_slice = $self->transcript->feature_Slice(); #### returns with strand same as feature
|
|
|
1269 my $tr_vf = $self->variation_feature->transfer($ref_slice);
|
|
|
1270
|
|
|
1271 # Return undef if this VariationFeature does not fall within the supplied feature.
|
|
|
1272 return undef if (!defined $tr_vf ||
|
|
|
1273 $tr_vf->start < 1 ||
|
|
|
1274 $tr_vf->end < 1 ||
|
|
|
1275 $tr_vf->start > ($self->transcript->end - $self->transcript->start + 1) ||
|
|
|
1276 $tr_vf->end > ($self->transcript->end - $self->transcript->start + 1));
|
|
|
1277
|
|
|
1278 return( $tr_vf->start() , $tr_vf->end(), $ref_slice);
|
|
|
1279 }
|
|
|
1280
|
|
|
1281
|
|
|
1282
|
|
|
1283 ### HGVS: get variant position in transcript
|
|
|
1284
|
|
|
1285 # intron:
|
|
|
1286 # If the position is in an intron, the boundary position of the closest exon and
|
|
|
1287 # a + or - offset into the intron is returned.
|
|
|
1288 # Ordered by genome forward not 5' -> 3'
|
|
|
1289
|
|
|
1290 # upstream:
|
|
|
1291 # If the position is 5' of the start codon, it is reported relative to the start codon
|
|
|
1292 # (-1 being the last nucleotide before the 'A' of ATG).
|
|
|
1293
|
|
|
1294 #downstream:
|
|
|
1295 # If the position is 3' pf the stop codon, it is reported with a '*' prefix and the offset
|
|
|
1296 # from the start codon (*1 being the first nucleotide after the last position of the stop codon)
|
|
|
1297
|
|
|
1298 sub _get_cDNA_position {
|
|
|
1299
|
|
|
1300 my $self = shift;
|
|
|
1301 my $position = shift; ### start or end of variation
|
|
|
1302
|
|
|
1303 my $transcript = $self->transcript();
|
|
|
1304 my $strand = $transcript->strand();
|
|
|
1305
|
|
|
1306 #### TranscriptVariation start/stop coord relative to transcript
|
|
|
1307 #### Switch to chromosome coordinates taking into account strand
|
|
|
1308 $position = ( $strand > 0 ?
|
|
|
1309 ( $self->transcript->start() + $position - 1 ) :
|
|
|
1310 ( $self->transcript->end() - $position + 1));
|
|
|
1311
|
|
|
1312
|
|
|
1313
|
|
|
1314 # Get all exons and sort them in positional order
|
|
|
1315 my @exons = sort {$a->start() <=> $b->start()} @{$transcript->get_all_Exons()};
|
|
|
1316 my $n_exons = scalar(@exons);
|
|
|
1317
|
|
|
1318 my $cdna_position;
|
|
|
1319 # Loop over the exons and get the coordinates of the variation in exon+intron notation
|
|
|
1320 for (my $i=0; $i<$n_exons; $i++) {
|
|
|
1321
|
|
|
1322 # Skip if the start point is beyond this exon
|
|
|
1323 next if ($position > $exons[$i]->end());
|
|
|
1324
|
|
|
1325
|
|
|
1326 # EXONIC: If the start coordinate is within this exon
|
|
|
1327 if ($position >= $exons[$i]->start()) {
|
|
|
1328 # Get the cDNA start coordinate of the exon and add the number of nucleotides from the exon boundary to the variation
|
|
|
1329 # If the transcript is in the opposite direction, count from the end instead
|
|
|
1330 $cdna_position = $exons[$i]->cdna_start($transcript) + ( $strand > 0 ?
|
|
|
1331 ( $position - $exons[$i]->start ) :
|
|
|
1332 ( $exons[$i]->end() - $position )
|
|
|
1333 );
|
|
|
1334 last; #### last exon checked
|
|
|
1335 }
|
|
|
1336 ## INTRONIC
|
|
|
1337 # Else the start coordinate is between this exon and the previous one, determine which one is closest and get coordinates relative to that one
|
|
|
1338 else {
|
|
|
1339
|
|
|
1340 my $updist = ($position - $exons[$i-1]->end());
|
|
|
1341 my $downdist = ($exons[$i]->start() - $position);
|
|
|
1342
|
|
|
1343 # If the distance to the upstream exon is the shortest, or equal and in the positive orientation, use that
|
|
|
1344 if ($updist < $downdist || ($updist == $downdist && $strand >= 0)) {
|
|
|
1345
|
|
|
1346 # If the orientation is reversed, we should use the cDNA start and a '-' offset
|
|
|
1347 $cdna_position = ($strand >= 0 ?
|
|
|
1348 $exons[$i-1]->cdna_end($transcript) . '+' :
|
|
|
1349 $exons[$i-1]->cdna_start($transcript) . '-')
|
|
|
1350 . $updist;
|
|
|
1351 }
|
|
|
1352 # Else if downstream is shortest...
|
|
|
1353 else {
|
|
|
1354 # If the orientation is reversed, we should use the cDNA end and a '+' offset
|
|
|
1355 $cdna_position = ($strand >= 0 ?
|
|
|
1356 $exons[$i]->cdna_start($transcript) . '-' :
|
|
|
1357 $exons[$i]->cdna_end($transcript) . '+')
|
|
|
1358 . $downdist;
|
|
|
1359 }
|
|
|
1360 last; ## last exon checked
|
|
|
1361 }
|
|
|
1362 }
|
|
|
1363
|
|
|
1364 # Shift the position to make it relative to the start & stop codons
|
|
|
1365 my $start_codon = $transcript->cdna_coding_start();
|
|
|
1366 my $stop_codon = $transcript->cdna_coding_end();
|
|
|
1367
|
|
|
1368 # Disassemble the cDNA coordinate into the exon and intron parts
|
|
|
1369 ### just built this now taking it appart again
|
|
|
1370 my ($cdna_coord, $intron_offset) = $cdna_position =~ m/([0-9]+)([\+\-][0-9]+)?/;
|
|
|
1371
|
|
|
1372
|
|
|
1373 # Start by correcting for the stop codon
|
|
|
1374 if (defined($stop_codon) && $cdna_coord > $stop_codon) {
|
|
|
1375 # Get the offset from the stop codon
|
|
|
1376 $cdna_coord -= $stop_codon;
|
|
|
1377 # Prepend a * to indicate the position is in the 3' UTR
|
|
|
1378 $cdna_coord = '*' . $cdna_coord;
|
|
|
1379 }
|
|
|
1380 elsif (defined($start_codon)) {
|
|
|
1381 # If the position is beyond the start codon, add 1 to get the correct offset
|
|
|
1382 $cdna_coord += ($cdna_coord >= $start_codon);
|
|
|
1383 # Subtract the position of the start codon
|
|
|
1384 $cdna_coord -= $start_codon;
|
|
|
1385 }
|
|
|
1386
|
|
|
1387 # Re-assemble the cDNA position [ return exon num & offset & direction for intron eg. 142+363]
|
|
|
1388 $cdna_position = $cdna_coord . (defined($intron_offset) ? $intron_offset : '');
|
|
|
1389 return $cdna_position;
|
|
|
1390 }
|
|
|
1391
|
|
|
1392
|
|
|
1393 1;
|
|
|
1394
|
|
|
1395
|
