Mercurial > repos > willmclaren > ensembl_vep

comparison variant_effect_predictor/Bio/LiveSeq/Mutator.pm @ 0:21066c0abaf5 draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Uploaded
author willmclaren
date Fri, 03 Aug 2012 10:04:48 -0400
parents
children
comparison
equal deleted inserted replaced
-1:000000000000 0:21066c0abaf5
1 # $Id: Mutator.pm,v 1.26 2002/10/22 07:38:34 lapp Exp $
2 #
3 # bioperl module for Bio::LiveSeq::Mutator
4 #
5 # Cared for by Heikki Lehvaslaiho <heikki@ebi.ac.uk>
6 #
7 # Copyright Joseph Insana
8 #
9 # You may distribute this module under the same terms as perl itself
10 #
11 # POD documentation - main docs before the code
12
13 =head1 NAME
14
15 Bio::LiveSeq::Mutator - Package mutating LiveSequences
16
17 =head1 SYNOPSIS
18
19 # $gene is a Bio::LiveSeq::Gene object
20 my $mutate = Bio::LiveSeq::Mutator->new('-gene' => $gene,
21 '-numbering' => "coding"
22 );
23 # $mut is a Bio::LiveSeq::Mutation object
24 $mutate->add_Mutation($mut);
25 # $results is a Bio::Variation::SeqDiff object
26 my $results=$mutate->change_gene();
27 if ($results) {
28 my $out = Bio::Variation::IO->new( '-format' => 'flat');
29 $out->write($results);
30 }
31
32 =head1 DESCRIPTION
33
34 This class mutates Bio::LiveSeq::Gene objects and returns a
35 Bio::Variation::SeqDiff object. Mutations are described as
36 Bio::LiveSeq::Mutation objects. See L<Bio::LiveSeq::Gene>,
37 L<Bio::Variation::SeqDiff>, and L<Bio::LiveSeq::Mutation> for details.
38
39 =head1 FEEDBACK
40
41
42 User feedback is an integral part of the evolution of this and other
43 Bioperl modules. Send your comments and suggestions preferably to the
44 Bioperl mailing lists Your participation is much appreciated.
45
46 bioperl-l@bioperl.org - General discussion
47 http://bio.perl.org/MailList.html - About the mailing lists
48
49 =head2 Reporting Bugs
50
51 report bugs to the Bioperl bug tracking system to help us keep track
52 the bugs and their resolution. Bug reports can be submitted via
53 email or the web:
54
55 bioperl-bugs@bio.perl.org
56 http://bugzilla.bioperl.org/
57
58 =head1 AUTHOR - Heikki Lehvaslaiho & Joseph A.L. Insana
59
60 Email: heikki@ebi.ac.uk
61 insana@ebi.ac.uk, jinsana@gmx.net
62
63 Address:
64
65 EMBL Outstation, European Bioinformatics Institute
66 Wellcome Trust Genome Campus, Hinxton
67 Cambs. CB10 1SD, United Kingdom
68
69 =head1 APPENDIX
70
71 The rest of the documentation details each of the object
72 methods. Internal methods are usually preceded with a _
73
74 =cut
75
76 # Let the code begin...
77
78 package Bio::LiveSeq::Mutator;
79 use vars qw(@ISA);
80 use strict;
81
82 use vars qw($VERSION @ISA);
83 use Bio::Variation::SeqDiff;
84 use Bio::Variation::DNAMutation;
85 use Bio::Variation::RNAChange;
86 use Bio::Variation::AAChange;
87 use Bio::Variation::Allele;
88 use Bio::LiveSeq::Mutation;
89
90 #use integer;
91 # Object preamble - inheritance
92
93 use Bio::Root::Root;
94
95 @ISA = qw( Bio::Root::Root );
96
97 sub new {
98 my($class,@args) = @_;
99 my $self;
100 $self = {};
101 bless $self, $class;
102
103 my ($gene, $numbering) =
104 $self->_rearrange([qw(GENE
105 NUMBERING
106 )],
107 @args);
108
109 $self->{ 'mutations' } = [];
110
111 $gene && $self->gene($gene);
112 $numbering && $self->numbering($numbering);
113
114 #class constant;
115 $self->{'flanklen'} = 25;
116 return $self; # success - we hope!
117 }
118
119 =head2 gene
120
121 Title : gene
122 Usage : $mutobj = $obj->gene;
123 : $mutobj = $obj->gene($objref);
124 Function:
125
126 Returns or sets the link-reference to a
127 Bio::LiveSeq::Gene object. If no value has ben set, it
128 will return undef
129
130 Returns : an object reference or undef
131 Args : a Bio::LiveSeq::Gene
132
133 See L<Bio::LiveSeq::Gene> for more information.
134
135 =cut
136
137 sub gene {
138 my ($self,$value) = @_;
139 if (defined $value) {
140 if( ! $value->isa('Bio::LiveSeq::Gene') ) {
141 $self->throw("Is not a Bio::LiveSeq::Gene object but a [$value]");
142 return undef;
143 }
144 else {
145 $self->{'gene'} = $value;
146 }
147 }
148 unless (exists $self->{'gene'}) {
149 return (undef);
150 } else {
151 return $self->{'gene'};
152 }
153 }
154
155
156 =head2 numbering
157
158 Title : numbering
159 Usage : $obj->numbering();
160 Function:
161
162 Sets and returns coordinate system used in positioning the
163 mutations. See L<change_gene> for details.
164
165 Example :
166 Returns : string
167 Args : string (coding [transcript number] | gene | entry)
168
169 =cut
170
171
172 sub numbering {
173 my ($self,$value) = @_;
174 if( defined $value) {
175 if ($value =~ /(coding)( )?(\d+)?/ or $value eq 'entry' or $value eq 'gene') {
176 $self->{'numbering'} = $value;
177 } else { # defaulting to 'coding'
178 $self->{'numbering'} = 'coding';
179 }
180 }
181 unless (exists $self->{'numbering'}) {
182 return 'coding';
183 } else {
184 return $self->{'numbering'};
185 }
186 }
187
188 =head2 add_Mutation
189
190 Title : add_Mutation
191 Usage : $self->add_Mutation($ref)
192 Function: adds a Bio::LiveSeq::Mutation object
193 Example :
194 Returns :
195 Args : a Bio::LiveSeq::Mutation
196
197 See L<Bio::LiveSeq::Mutation> for more information.
198
199 =cut
200
201 sub add_Mutation{
202 my ($self,$value) = @_;
203 if( $value->isa('Bio::Liveseq::Mutation') ) {
204 my $com = ref $value;
205 $self->throw("Is not a Mutation object but a [$com]" );
206 return undef;
207 }
208 if (! $value->pos) {
209 $self->warn("No value for mutation position in the sequence!");
210 return undef;
211 }
212 if (! $value->seq && ! $value->len) {
213 $self->warn("Either mutated sequence or length of the deletion must be given!");
214 return undef;
215 }
216 push(@{$self->{'mutations'}},$value);
217 }
218
219 =head2 each_Mutation
220
221 Title : each_Mutation
222 Usage : foreach $ref ( $a->each_Mutation )
223 Function: gets an array of Bio::LiveSeq::Mutation objects
224 Example :
225 Returns : array of Mutations
226 Args :
227
228 See L<Bio::LiveSeq::Mutation> for more information.
229
230 =cut
231
232 sub each_Mutation{
233 my ($self) = @_;
234 return @{$self->{'mutations'}};
235 }
236
237
238 =head2 mutation
239
240 Title : mutation
241 Usage : $mutobj = $obj->mutation;
242 : $mutobj = $obj->mutation($objref);
243 Function:
244
245 Returns or sets the link-reference to the current mutation
246 object. If the value is not set, it will return undef.
247 Internal method.
248
249 Returns : an object reference or undef
250
251 =cut
252
253
254 sub mutation {
255 my ($self,$value) = @_;
256 if (defined $value) {
257 if( ! $value->isa('Bio::LiveSeq::Mutation') ) {
258 $self->throw("Is not a Bio::LiveSeq::Mutation object but a [$value]");
259 return undef;
260 }
261 else {
262 $self->{'mutation'} = $value;
263 }
264 }
265 unless (exists $self->{'mutation'}) {
266 return (undef);
267 } else {
268 return $self->{'mutation'};
269 }
270 }
271
272 =head2 DNA
273
274 Title : DNA
275 Usage : $mutobj = $obj->DNA;
276 : $mutobj = $obj->DNA($objref);
277 Function:
278
279 Returns or sets the reference to the LiveSeq object holding
280 the reference sequence. If there is no link, it will return
281 undef.
282 Internal method.
283
284 Returns : an object reference or undef
285
286 =cut
287
288 sub DNA {
289 my ($self,$value) = @_;
290 if (defined $value) {
291 if( ! $value->isa('Bio::LiveSeq::DNA') and ! $value->isa('Bio::LiveSeq::Transcript') ) {
292 $self->throw("Is not a Bio::LiveSeq::DNA/Transcript object but a [$value]");
293 return undef;
294 }
295 else {
296 $self->{'DNA'} = $value;
297 }
298 }
299 unless (exists $self->{'DNA'}) {
300 return (undef);
301 } else {
302 return $self->{'DNA'};
303 }
304 }
305
306
307 =head2 RNA
308
309 Title : RNA
310 Usage : $mutobj = $obj->RNA;
311 : $mutobj = $obj->RNA($objref);
312 Function:
313
314 Returns or sets the reference to the LiveSeq object holding
315 the reference sequence. If the value is not set, it will return
316 undef.
317 Internal method.
318
319 Returns : an object reference or undef
320
321 =cut
322
323
324 sub RNA {
325 my ($self,$value) = @_;
326 if (defined $value) {
327 if( ! $value->isa('Bio::LiveSeq::Transcript') ) {
328 $self->throw("Is not a Bio::LiveSeq::RNA/Transcript object but a [$value]");
329 return undef;
330 }
331 else {
332 $self->{'RNA'} = $value;
333 }
334 }
335 unless (exists $self->{'RNA'}) {
336 return (undef);
337 } else {
338 return $self->{'RNA'};
339 }
340 }
341
342
343 =head2 dnamut
344
345 Title : dnamut
346 Usage : $mutobj = $obj->dnamut;
347 : $mutobj = $obj->dnamut($objref);
348 Function:
349
350 Returns or sets the reference to the current DNAMutation object.
351 If the value is not set, it will return undef.
352 Internal method.
353
354 Returns : a Bio::Variation::DNAMutation object or undef
355
356 See L<Bio::Variation::DNAMutation> for more information.
357
358 =cut
359
360
361 sub dnamut {
362 my ($self,$value) = @_;
363 if (defined $value) {
364 if( ! $value->isa('Bio::Variation::DNAMutation') ) {
365 $self->throw("Is not a Bio::Variation::DNAMutation object but a [$value]");
366 return undef;
367 }
368 else {
369 $self->{'dnamut'} = $value;
370 }
371 }
372 unless (exists $self->{'dnamut'}) {
373 return (undef);
374 } else {
375 return $self->{'dnamut'};
376 }
377 }
378
379
380 =head2 rnachange
381
382 Title : rnachange
383 Usage : $mutobj = $obj->rnachange;
384 : $mutobj = $obj->rnachange($objref);
385 Function:
386
387 Returns or sets the reference to the current RNAChange object.
388 If the value is not set, it will return undef.
389 Internal method.
390
391 Returns : a Bio::Variation::RNAChange object or undef
392
393 See L<Bio::Variation::RNAChange> for more information.
394
395 =cut
396
397
398 sub rnachange {
399 my ($self,$value) = @_;
400 if (defined $value) {
401 if( ! $value->isa('Bio::Variation::RNAChange') ) {
402 $self->throw("Is not a Bio::Variation::RNAChange object but a [$value]");
403 return undef;
404 }
405 else {
406 $self->{'rnachange'} = $value;
407 }
408 }
409 unless (exists $self->{'rnachange'}) {
410 return (undef);
411 } else {
412 return $self->{'rnachange'};
413 }
414 }
415
416
417 =head2 aachange
418
419 Title : aachange
420 Usage : $mutobj = $obj->aachange;
421 : $mutobj = $obj->aachange($objref);
422 Function:
423
424 Returns or sets the reference to the current AAChange object.
425 If the value is not set, it will return undef.
426 Internal method.
427
428 Returns : a Bio::Variation::AAChange object or undef
429
430 See L<Bio::Variation::AAChange> for more information.
431
432 =cut
433
434
435 sub aachange {
436 my ($self,$value) = @_;
437 if (defined $value) {
438 if( ! $value->isa('Bio::Variation::AAChange') ) {
439 $self->throw("Is not a Bio::Variation::AAChange object but a [$value]");
440 return undef;
441 }
442 else {
443 $self->{'aachange'} = $value;
444 }
445 }
446 unless (exists $self->{'aachange'}) {
447 return (undef);
448 } else {
449 return $self->{'aachange'};
450 }
451 }
452
453
454 =head2 exons
455
456 Title : exons
457 Usage : $mutobj = $obj->exons;
458 : $mutobj = $obj->exons($objref);
459 Function:
460
461 Returns or sets the reference to a current array of Exons.
462 If the value is not set, it will return undef.
463 Internal method.
464
465 Returns : an array of Bio::LiveSeq::Exon objects or undef
466
467 See L<Bio::LiveSeq::Exon> for more information.
468
469 =cut
470
471
472 sub exons {
473 my ($self,$value) = @_;
474 if (defined $value) {
475 $self->{'exons'} = $value;
476 }
477 unless (exists $self->{'exons'}) {
478 return (undef);
479 } else {
480 return $self->{'exons'};
481 }
482 }
483
484 =head2 change_gene_with_alignment
485
486 Title : change_gene_with_alignment
487 Usage : $results=$mutate->change_gene_with_alignment($aln);
488
489 Function:
490
491 Returns a Bio::Variation::SeqDiff object containing the
492 results of the changes in the alignment. The alignment has
493 to be pairwise and have one sequence named 'QUERY', the
494 other one is assumed to be a part of the sequence from
495 $gene.
496
497 This method offers a shortcut to change_gene and
498 automates the creation of Bio::LiveSeq::Mutation objects.
499 Use it with almost identical sequnces, e.g. to locate a SNP.
500
501 Args : Bio::SimpleAlign object representing a short local alignment
502 Returns : Bio::Variation::SeqDiff object or 0 on error
503
504 See L<Bio::LiveSeq::Mutation>, L<Bio::SimpleAlign>, and
505 L<Bio::Variation::SeqDiff> for more information.
506
507 =cut
508
509 sub change_gene_with_alignment {
510 my ($self, $aln) = @_;
511
512 #
513 # Sanity checks
514 #
515
516 $self->throw("Argument is not a Bio::SimpleAlign object but a [$aln]")
517 unless $aln->isa('Bio::SimpleAlign');
518 $self->throw("'Pairwise alignments only, please")
519 if $aln->no_sequences != 2;
520
521 # find out the order the two sequences are given
522 my $queryseq_pos = 1; #default
523 my $refseq_pos = 2;
524 unless ($aln->get_seq_by_pos(1)->id eq 'QUERY') {
525 carp('Query sequence has to be named QUERY')
526 if $aln->get_seq_by_pos(2)->id ne 'QUERY';
527 $queryseq_pos = 2; # alternative
528 $refseq_pos = 1;
529 }
530
531 # trim the alignment
532 my $start = $aln->column_from_residue_number('QUERY', 1);
533 my $end = $aln->column_from_residue_number('QUERY',
534 $aln->get_seq_by_pos($queryseq_pos)->end );
535
536 my $aln2 = $aln->slice($start, $end);
537
538 #
539 # extracting mutations
540 #
541
542 my $cs = $aln2->consensus_string(51);
543 my $queryseq = $aln2->get_seq_by_pos($queryseq_pos);
544 my $refseq = $aln2->get_seq_by_pos($refseq_pos);
545
546 while ( $cs =~ /(\?+)/g) {
547 # pos in local coordinates
548 my $pos = pos($cs) - length($1) + 1;
549 my $mutation = create_mutation($self,
550 $refseq,
551 $queryseq,
552 $pos,
553 CORE::length($1)
554 );
555 # reset pos to refseq coordinates
556 $pos += $refseq->start - 1;
557 $mutation->pos($pos);
558
559 $self->add_Mutation($mutation);
560 }
561 return $self->change_gene();
562 }
563
564 =head2 create_mutation
565
566 Title : create_mutation
567 Usage :
568 Function:
569
570 Formats sequence differences from two sequences into
571 Bio::LiveSeq::Mutation objects which can be applied to a
572 gene.
573
574 To keep it generic, sequence arguments need not to be
575 Bio::LocatableSeq. Coordinate change to parent sequence
576 numbering needs to be done by the calling code.
577
578 Called from change_gene_with_alignment
579
580 Args : Bio::PrimarySeqI inheriting object for the reference sequence
581 Bio::PrimarySeqI inheriting object for the query sequence
582 integer for the start position of the local sequence difference
583 integer for the length of the sequence difference
584 Returns : Bio::LiveSeq::Mutation object
585
586 =cut
587
588 sub create_mutation {
589 my ($self, $refseq, $queryseq, $pos, $len) = @_;
590
591 $self->throw("Is not a Bio::PrimarySeqI object but a [$refseq]")
592 unless $refseq->isa('Bio::PrimarySeqI');
593 $self->throw("Is not a Bio::PrimarySeqI object but a [$queryseq]")
594 unless $queryseq->isa('Bio::PrimarySeqI');
595 $self->throw("Position is not a positive integer but [$pos]")
596 unless $pos =~ /^\+?\d+$/;
597 $self->throw("Length is not a positive integer but [$len]")
598 unless $len =~ /^\+?\d+$/;
599
600 my $mutation;
601 my $refstring = $refseq->subseq($pos, $pos + $len - 1);
602 my $varstring = $queryseq->subseq($pos, $pos + $len - 1);
603
604 if ($len == 1 and $refstring =~ /[^\.\-\*\?]/ and
605 $varstring =~ /[^\.\-\*\?]/ ) { #point
606 $mutation = new Bio::LiveSeq::Mutation (-seq => $varstring,
607 -pos => $pos,
608 );
609 }
610 elsif ( $refstring =~ /^[^\.\-\*\?]+$/ and
611 $varstring !~ /^[^\.\-\*\?]+$/ ) { # deletion
612 $mutation = new Bio::LiveSeq::Mutation (-pos => $pos,
613 -len => $len
614 );
615 }
616 elsif ( $refstring !~ /^[^\.\-\*\?]+$/ and
617 $varstring =~ /^[^\.\-\*\?]+$/ ) { # insertion
618 $mutation = new Bio::LiveSeq::Mutation (-seq => $varstring,
619 -pos => $pos,
620 -len => 0
621 );
622 } else { # complex
623 $mutation = new Bio::LiveSeq::Mutation (-seq => $varstring,
624 -pos => $pos,
625 -len => $len
626 );
627 }
628
629 return $mutation;
630 }
631
632 =head2 change_gene
633
634 Title : change_gene
635 Usage : my $mutate = Bio::LiveSeq::Mutator->new(-gene => $gene,
636 numbering => "coding"
637 );
638 # $mut is Bio::LiveSeq::Mutation object
639 $mutate->add_Mutation($mut);
640 my $results=$mutate->change_gene();
641
642 Function:
643
644 Returns a Bio::Variation::SeqDiff object containing the
645 results of the changes performed according to the
646 instructions present in Mutation(s). The -numbering
647 argument decides what molecule is being changed and what
648 numbering scheme being used:
649
650 -numbering => "entry"
651
652 determines the DNA level, using the numbering from the
653 beginning of the sequence
654
655 -numbering => "coding"
656
657 determines the RNA level, using the numbering from the
658 beginning of the 1st transcript
659
660 Alternative transcripts can be used by specifying
661 "coding 2" or "coding 3" ...
662
663 -numbering => "gene"
664
665 determines the DNA level, using the numbering from the
666 beginning of the 1st transcript and inluding introns.
667 The meaning equals 'coding' if the reference molecule
668 is cDNA.
669
670 Args : Bio::LiveSeq::Gene object
671 Bio::LiveSeq::Mutation object(s)
672 string specifying a numbering scheme (defaults to 'coding')
673 Returns : Bio::Variation::SeqDiff object or 0 on error
674
675 =cut
676
677 sub change_gene {
678 my ($self) = @_;
679
680 #
681 # Sanity check
682 #
683 unless ($self->gene) {
684 $self->warn("Input object Bio::LiveSeq::Gene is not given");
685 return 0;
686 }
687 #
688 # Setting the reference sequence based on -numbering
689 #
690 my @transcripts=@{$self->gene->get_Transcripts};
691 my $refseq; # will hold Bio::LiveSeq:Transcript object or Bio::LiveSeq::DNA
692
693 # 'gene' eq 'coding' if reference sequence is cDNA
694 $self->numbering ('coding') if $self->gene->get_DNA->alphabet eq 'rna' and $self->numbering eq 'gene';
695
696 if ($self->numbering =~ /(coding)( )?(\d+)?/ ) {
697 $self->numbering($1);
698 my $transnumber = $3;
699 $transnumber-- if $3; # 1 -> 0, 2 -> 1
700 if ($transnumber && $transnumber >= 0 && $transnumber <= $#transcripts) {
701 $refseq=$transcripts[$transnumber];
702 } else {
703 $transnumber && $self->warn("The alternative transcript number ". $transnumber+1 .
704 "- does not exist. Reverting to the 1st transcript\n");
705 $refseq=$transcripts[0];
706 }
707 } else {
708 $refseq=$transcripts[0]->{'seq'};
709 }
710 #
711 # Recording the state: SeqDiff object creation ?? transcript no.??
712 #
713 my $seqDiff = Bio::Variation::SeqDiff->new();
714 $seqDiff->alphabet($self->gene->get_DNA->alphabet);
715 $seqDiff->numbering($self->numbering);
716 my ($DNAobj, $RNAobj);
717 if ($refseq->isa("Bio::LiveSeq::Transcript")) {
718 $self->RNA($refseq);
719 $self->DNA($refseq->{'seq'});
720 $seqDiff->rna_ori($refseq->seq );
721 $seqDiff->aa_ori($refseq->get_Translation->seq);
722 } else {
723 $self->DNA($refseq);
724 $self->RNA($transcripts[0]);
725 $seqDiff->rna_ori($self->RNA->seq);
726 $seqDiff->aa_ori($self->RNA->get_Translation->seq);
727 }
728 $seqDiff->dna_ori($self->DNA->seq);
729 # put the accession number into the SeqDiff object ID
730 $seqDiff->id($self->DNA->accession_number);
731
732 # the atg_offset takes in account that DNA object could be a subset of the
733 # whole entry (via the light_weight loader)
734 my $atg_label=$self->RNA->start;
735 my $atg_offset=$self->DNA->position($atg_label)+($self->DNA->start)-1;
736 $seqDiff->offset($atg_offset - 1);
737 $self->DNA->coordinate_start($atg_label);
738
739 my @exons = $self->RNA->all_Exons;
740 $seqDiff->cds_end($exons[$#exons]->end);
741
742 #
743 # Converting mutation positions to labels
744 #
745 $self->warn("no mutations"), return 0
746 unless $self->_mutationpos2label($refseq, $seqDiff);
747
748 # need to add more than one rna & aa
749 #foreach $transcript (@transcripts) {
750 # $seqDiff{"ori_transcript_${i}_seq"}=$transcript->seq;
751 # $seqDiff{"ori_translation_${i}_seq"}=$transcript->get_Translation->seq;
752 #}
753
754 # do changes
755 my $k;
756 foreach my $mutation ($self->each_Mutation) {
757 next unless $mutation->label > 0;
758 $self->mutation($mutation);
759
760 $mutation->issue(++$k);
761 #
762 # current position on the transcript
763 #
764 if ($self->numbering =~ /coding/) {
765 $mutation->transpos($mutation->pos); # transpos given by user
766 } else {
767 #transpos of label / It will be 0 if mutating an intron, negative if upstream of ATG
768 $mutation->transpos($self->RNA->position($mutation->label,$atg_label));
769 }
770 #
771 # Calculate adjacent labels based on the position on the current sequence
772 #
773 $mutation->prelabel($self->DNA->label(-1, $mutation->label)); # 1 before label
774 if ($mutation->len == 0) {
775 $mutation->postlabel($mutation->label);
776 $mutation->lastlabel($mutation->label);
777 } elsif ($mutation->len == 1) {
778 $mutation->lastlabel($mutation->label); # last nucleotide affected
779 $mutation->postlabel($self->DNA->label(2,$mutation->lastlabel)); # $len after label
780 } else {
781 $mutation->lastlabel($self->DNA->label($mutation->len,$mutation->label));
782 $mutation->postlabel($self->DNA->label(2,$mutation->lastlabel));
783 }
784 my $dnamut = $self->_set_DNAMutation($seqDiff);
785 #
786 #
787 #
788 if ($self->_rnaAffected) {
789 $self->_set_effects($seqDiff, $dnamut);
790 }
791 elsif ($seqDiff->offset != 0 and $dnamut->region ne 'intron') {
792 $self->_untranslated ($seqDiff, $dnamut);
793 } else {
794 #$self->warn("Mutation starts outside coding region, RNAChange object not created");
795 }
796
797 #########################################################################
798 # Mutations are done here! #
799 $refseq->labelchange($mutation->seq, $mutation->label, $mutation->len); #
800 #########################################################################
801
802 $self->_post_mutation ($seqDiff);
803
804 $self->dnamut(undef);
805 $self->rnachange(undef);
806 $self->aachange(undef);
807 $self->exons(undef);
808 }
809 # record the final state of all three sequences
810 $seqDiff->dna_mut($self->DNA->seq);
811 $seqDiff->rna_mut($self->RNA->seq);
812 if ($refseq->isa("Bio::LiveSeq::Transcript")) {
813 $seqDiff->aa_mut($refseq->get_Translation->seq);
814 } else {
815 $seqDiff->aa_mut($self->RNA->get_Translation->seq);
816 }
817
818 #$seqDiff{mut_dna_seq}=$gene->get_DNA->seq;
819 #my $i=1;
820 #foreach $transcript (@transcripts) {
821 # $seqDiff{"mut_transcript_${i}_seq"}=$transcript->seq;
822 # $seqDiff{"mut_translation_${i}_seq"}=$transcript->get_Translation->seq;
823 #}
824 return $seqDiff;
825 }
826
827 =head2 _mutationpos2label
828
829 Title : _mutationpos2label
830 Usage :
831 Function: converts mutation positions into labels
832 Example :
833 Returns : number of valid mutations
834 Args : LiveSeq sequence object
835
836 =cut
837
838 sub _mutationpos2label {
839 my ($self, $refseq, $SeqDiff) = @_;
840 my $count;
841 my @bb = @{$self->{'mutations'}};
842 my $cc = scalar @bb;
843 foreach my $mut (@{$self->{'mutations'}}) {
844 # if ($self->numbering eq 'gene' and $mut->pos < 1) {
845 # my $tmp = $mut->pos;
846 # print STDERR "pos: ", "$tmp\n";
847 # $tmp++ if $tmp < 1;
848 # $tmp += $SeqDiff->offset;
849 # print STDERR "pos2: ", "$tmp\n";
850 # $mut->pos($tmp);
851 # }
852 # elsif ($self->numbering eq 'entry') {
853 if ($self->numbering eq 'entry') {
854 my $tmp = $mut->pos;
855 $tmp -= $SeqDiff->offset;
856 $tmp-- if $tmp < 1;
857 $mut->pos($tmp);
858 }
859
860 my $label = $refseq->label($mut->pos); # get the label for the position
861 $mut->label($label), $count++ if $label > 0 ;
862 #print STDERR "x", $mut->pos,'|' ,$mut->label, "\n";
863 }
864 return $count;
865 }
866
867 #
868 # Calculate labels around mutated nucleotide
869 #
870
871 =head2 _set_DNAMutation
872
873 Title : _set_DNAMutation
874 Usage :
875 Function:
876
877 Stores DNA level mutation attributes before mutation into
878 Bio::Variation::DNAMutation object. Links it to SeqDiff
879 object.
880
881 Example :
882 Returns : Bio::Variation::DNAMutation object
883 Args : Bio::Variation::SeqDiff object
884
885 See L<Bio::Variation::DNAMutation> and L<Bio::Variation::SeqDiff>.
886
887 =cut
888
889 sub _set_DNAMutation {
890 my ($self, $seqDiff) = @_;
891
892 my $dnamut_start = $self->mutation->label - $seqDiff->offset;
893 # if negative DNA positions (before ATG)
894 $dnamut_start-- if $dnamut_start <= 0;
895 my $dnamut_end;
896 ($self->mutation->len == 0 or $self->mutation->len == 1) ?
897 ($dnamut_end = $dnamut_start) :
898 ($dnamut_end = $dnamut_start+$self->mutation->len);
899 #print "start:$dnamut_start, end:$dnamut_end\n";
900 my $dnamut = Bio::Variation::DNAMutation->new(-start => $dnamut_start,
901 -end => $dnamut_end,
902 );
903 $dnamut->mut_number($self->mutation->issue);
904 $dnamut->isMutation(1);
905 my $da_m = Bio::Variation::Allele->new;
906 $da_m->seq($self->mutation->seq) if $self->mutation->seq;
907 $dnamut->allele_mut($da_m);
908 $dnamut->add_Allele($da_m);
909 # allele_ori
910 my $allele_ori = $self->DNA->labelsubseq($self->mutation->prelabel,
911 undef,
912 $self->mutation->postlabel); # get seq
913 chop $allele_ori; # chop the postlabel nucleotide
914 $allele_ori=substr($allele_ori,1); # away the prelabel nucleotide
915 my $da_o = Bio::Variation::Allele->new;
916 $da_o->seq($allele_ori) if $allele_ori;
917 $dnamut->allele_ori($da_o);
918 ($self->mutation->len == 0) ?
919 ($dnamut->length($self->mutation->len)) : ($dnamut->length(CORE::length $allele_ori));
920 #print " --------------- $dnamut_start -$len- $dnamut_end -\n";
921 $seqDiff->add_Variant($dnamut);
922 $self->dnamut($dnamut);
923 $dnamut->mut_number($self->mutation->issue);
924 # setting proof
925 if ($seqDiff->numbering eq "entry" or $seqDiff->numbering eq "gene") {
926 $dnamut->proof('experimental');
927 } else {
928 $dnamut->proof('computed');
929 }
930 # how many nucleotides to store upstream and downstream of the change
931 my $flanklen = $self->{'flanklen'};
932 #print `date`, " flanking sequences start\n";
933 my $uplabel = $self->DNA->label(1-$flanklen,$self->mutation->prelabel); # this could be unavailable!
934
935 my $upstreamseq;
936 if ($uplabel > 0) {
937 $upstreamseq =
938 $self->DNA->labelsubseq($uplabel, undef, $self->mutation->prelabel);
939 } else { # from start (less than $flanklen nucleotides)
940 $upstreamseq =
941 $self->DNA->labelsubseq($self->DNA->start, undef, $self->mutation->prelabel);
942 }
943 $dnamut->upStreamSeq($upstreamseq);
944 my $dnstreamseq = $self->DNA->labelsubseq($self->mutation->postlabel, $flanklen);
945 $dnamut->dnStreamSeq($dnstreamseq); # $flanklen or less nucleotides
946 return $dnamut;
947 }
948
949
950 #
951 ### Check if mutation propagates to RNA (and AA) level
952 #
953 # side effect: sets intron/exon information
954 # returns a boolean value
955 #
956
957 sub _rnaAffected {
958 my ($self) = @_;
959 my @exons=$self->RNA->all_Exons;
960 my $RNAstart=$self->RNA->start;
961 my $RNAend=$self->RNA->end;
962 my ($firstexon,$before,$after,$i);
963 my ($rnaAffected) = 0;
964
965 # check for inserted labels (that require follows instead of <,>)
966 my $DNAend=$self->RNA->{'seq'}->end;
967 if ($self->mutation->prelabel > $DNAend or $self->mutation->postlabel > $DNAend) {
968 #this means one of the two labels is an inserted one
969 #(coming from a previous mutation. This would falsify all <,>
970 #checks, so the follow() has to be used
971 $self->warn("Attention, workaround not fully tested yet! Expect unpredictable results.\n");
972 if (($self->mutation->postlabel==$RNAstart) or (follows($self->mutation->postlabel,$RNAstart))) {
973 $self->warn("RNA not affected because change occurs before RNAstart");
974 }
975 elsif (($RNAend==$self->mutation->prelabel) or (follows($RNAend,$self->mutation->prelabel))) {
976 $self->warn("RNA not affected because change occurs after RNAend");
977 }
978 elsif (scalar @exons == 1) {
979 #no introns, just one exon
980 $rnaAffected = 1; # then RNA is affected!
981 } else {
982 # otherwise check for change occurring inside an intron
983 $firstexon=shift(@exons);
984 $before=$firstexon->end;
985
986 foreach $i (0..$#exons) {
987 $after=$exons[$i]->start;
988 if (follows($self->mutation->prelabel,$before) or
989 ($after==$self->mutation->prelabel) or
990 follows($after,$self->mutation->prelabel) or
991 follows($after,$self->mutation->postlabel)) {
992
993 $rnaAffected = 1;
994 # $i is number of exon and can be used for proximity check
995 }
996 $before=$exons[$i]->end;
997 }
998
999 }
1000 } else {
1001 my $strand = $exons[0]->strand;
1002 if (($strand == 1 and $self->mutation->postlabel <= $RNAstart) or
1003 ($strand != 1 and $self->mutation->postlabel >= $RNAstart)) {
1004 #$self->warn("RNA not affected because change occurs before RNAstart");
1005 $rnaAffected = 0;
1006 }
1007 elsif (($strand == 1 and $self->mutation->prelabel >= $RNAend) or
1008 ($strand != 1 and $self->mutation->prelabel <= $RNAend)) {
1009 #$self->warn("RNA not affected because change occurs after RNAend");
1010 $rnaAffected = 0;
1011 my $dist;
1012 if ($strand == 1){
1013 $dist = $self->mutation->prelabel - $RNAend;
1014 } else {
1015 $dist = $RNAend - $self->mutation->prelabel;
1016 }
1017 $self->dnamut->region_dist($dist);
1018 }
1019 elsif (scalar @exons == 1) {
1020 #if just one exon -> no introns,
1021 $rnaAffected = 1; # then RNA is affected!
1022 } else {
1023 # otherwise check for mutation occurring inside an intron
1024 $firstexon=shift(@exons);
1025 $before=$firstexon->end;
1026 if ( ($strand == 1 and $self->mutation->prelabel < $before) or
1027 ($strand == -1 and $self->mutation->prelabel > $before)
1028 ) {
1029 $rnaAffected = 1 ;
1030
1031 #print "Exon 1 : ", $firstexon->start, " - ", $firstexon->end, "<br>\n";
1032 my $afterdist = $self->mutation->prelabel - $firstexon->start;
1033 my $beforedist = $firstexon->end - $self->mutation->postlabel;
1034 my $exonvalue = $i + 1;
1035 $self->dnamut->region('exon');
1036 $self->dnamut->region_value($exonvalue);
1037 if ($afterdist < $beforedist) {
1038 $afterdist++;
1039 $afterdist++;
1040 $self->dnamut->region_dist($afterdist);
1041 #print "splice site $afterdist nt upstream!<br>";
1042 } else {
1043 $self->dnamut->region_dist($beforedist);
1044 #print "splice site $beforedist nt downstream!<br>";
1045 }
1046 } else {
1047 #print "first exon : ", $firstexon->start, " - ", $firstexon->end, "<br>\n";
1048 foreach $i (0..$#exons) {
1049 $after=$exons[$i]->start;
1050 #proximity test for intronic mutations
1051 if ( ($strand == 1 and
1052 $self->mutation->prelabel >= $before and
1053 $self->mutation->postlabel <= $after)
1054 or
1055 ($strand == -1 and
1056 $self->mutation->prelabel <= $before and
1057 $self->mutation->postlabel >= $after) ) {
1058 $self->dnamut->region('intron');
1059 #$self->dnamut->region_value($i);
1060 my $afterdist = $self->mutation->prelabel - $before;
1061 my $beforedist = $after - $self->mutation->postlabel;
1062 my $intronvalue = $i + 1;
1063 if ($afterdist < $beforedist) {
1064 $afterdist++;
1065 $self->dnamut->region_value($intronvalue);
1066 $self->dnamut->region_dist($afterdist);
1067 #print "splice site $afterdist nt upstream!<br>";
1068 } else {
1069 $self->dnamut->region_value($intronvalue);
1070 $self->dnamut->region_dist($beforedist * -1);
1071 #print "splice site $beforedist nt downstream!<br>";
1072 }
1073 $self->rnachange(undef);
1074 last;
1075 }
1076 #proximity test for exon mutations
1077 elsif ( ( $strand == 1 and
1078 $exons[$i]->start <= $self->mutation->prelabel and
1079 $exons[$i]->end >= $self->mutation->postlabel) or
1080 ( $strand == -1 and
1081 $exons[$i]->start >= $self->mutation->prelabel and
1082 $exons[$i]->end <= $self->mutation->postlabel) ) {
1083 $rnaAffected = 1;
1084
1085 my $afterdist = $self->mutation->prelabel - $exons[$i]->start;
1086 my $beforedist = $exons[$i]->end - $self->mutation->postlabel;
1087 my $exonvalue = $i + 1;
1088 $self->dnamut->region('exon');
1089 if ($afterdist < $beforedist) {
1090 $afterdist++;
1091 $self->dnamut->region_value($exonvalue+1);
1092 $self->dnamut->region_dist($afterdist);
1093 #print "splice site $afterdist nt upstream!<br>";
1094 } else {
1095 #$beforedist;
1096 $self->dnamut->region_value($exonvalue+1);
1097 $self->dnamut->region_dist($beforedist * -1);
1098 #print "splice site $beforedist nt downstream!<br>";
1099 }
1100 last;
1101 }
1102 $before=$exons[$i]->end;
1103 }
1104 }
1105 }
1106 }
1107 #$self->warn("RNA not affected because change occurs inside an intron");
1108 #return(0); # if still not returned, then not affected, return 0
1109 return $rnaAffected;
1110 }
1111
1112 #
1113 # ### Creation of RNA and AA variation objects
1114 #
1115
1116 =head2 _set_effects
1117
1118 Title : _set_effects
1119 Usage :
1120 Function:
1121
1122 Stores RNA and AA level mutation attributes before mutation
1123 into Bio::Variation::RNAChange and
1124 Bio::Variation::AAChange objects. Links them to
1125 SeqDiff object.
1126
1127 Example :
1128 Returns :
1129 Args : Bio::Variation::SeqDiff object
1130 Bio::Variation::DNAMutation object
1131
1132 See L<Bio::Variation::RNAChange>, L<Bio::Variation::RNAChange>,
1133 L<Bio::Variation::SeqDiff>, and L<Bio::Variation::DNAMutation>.
1134
1135 =cut
1136
1137 sub _set_effects {
1138 my ($self, $seqDiff, $dnamut) = @_;
1139 my ($rnapos_end, $upstreamseq, $dnstreamseq);
1140 my $flanklen = $self->{'flanklen'};
1141
1142 ($self->mutation->len == 0) ?
1143 ($rnapos_end = $self->mutation->transpos) :
1144 ($rnapos_end = $self->mutation->transpos + $self->mutation->len -1);
1145 my $rnachange = Bio::Variation::RNAChange->new(-start => $self->mutation->transpos,
1146 -end => $rnapos_end
1147 );
1148 $rnachange->isMutation(1);
1149
1150 # setting proof
1151 if ($seqDiff->numbering eq "coding") {
1152 $rnachange->proof('experimental');
1153 } else {
1154 $rnachange->proof('computed');
1155 }
1156
1157 $seqDiff->add_Variant($rnachange);
1158 $self->rnachange($rnachange);
1159 $rnachange->DNAMutation($dnamut);
1160 $dnamut->RNAChange($rnachange);
1161 $rnachange->mut_number($self->mutation->issue);
1162
1163 # setting the codon_position of the "start" nucleotide of the change
1164 $rnachange->codon_pos(($self->RNA->frame($self->mutation->label))+1); # codon_pos=frame+1
1165
1166 my @exons=$self->RNA->all_Exons;
1167 $self->exons(\@exons);
1168 #print `date`, " before flank, after exons. RNAObj query\n";
1169 # if cannot retrieve from Transcript, Transcript::upstream_seq will be used
1170 # before "fac7 g 65" bug discovered
1171 # $uplabel=$self->RNA->label(1-$flanklen,$prelabel);
1172 my $RNAprelabel=$self->RNA->label(-1,$self->mutation->label); # to fix fac7g65 bug
1173 # for the fix, all prelabel used in the next block have been changed to RNAprelabel
1174 my $uplabel=$self->RNA->label(1-$flanklen,$RNAprelabel);
1175 if ($self->RNA->valid($uplabel)) {
1176 $upstreamseq = $self->RNA->labelsubseq($uplabel, undef, $RNAprelabel);
1177 } else {
1178 $upstreamseq = $self->RNA->labelsubseq($self->RNA->start, undef, $RNAprelabel)
1179 if $self->RNA->valid($RNAprelabel);
1180 my $lacking=$flanklen-length($upstreamseq); # how many missing
1181 my $upstream_atg=$exons[0]->subseq(-$lacking,-1);
1182 $upstreamseq=$upstream_atg . $upstreamseq;
1183 }
1184
1185 $rnachange->upStreamSeq($upstreamseq);
1186
1187 # won't work OK if postlabel NOT in Transcript
1188 # now added RNApostlabel but this has to be /fully tested/
1189 # for the fix, all postlabel used in the next block have been changed to RNApostlabel
1190 my $RNApostlabel; # to fix fac7g64 bug
1191 if ($self->mutation->len == 0) {
1192 $RNApostlabel=$self->mutation->label;
1193 } else {
1194 my $mutlen = 1 + $self->mutation->len;
1195 $RNApostlabel=$self->RNA->label($mutlen,$self->mutation->label);
1196 }
1197 $dnstreamseq=$self->RNA->labelsubseq($RNApostlabel, $flanklen);
1198 if ($dnstreamseq eq '-1') { # if out of transcript was requested
1199 my $lastexon=$exons[-1];
1200 my $lastexonlength=$lastexon->length;
1201 $dnstreamseq=$self->RNA->labelsubseq($RNApostlabel); # retrieves till RNAend
1202 my $lacking=$flanklen-length($dnstreamseq); # how many missing
1203 my $downstream_stop=$lastexon->subseq($lastexonlength+1,undef,$lacking);
1204 $dnstreamseq .= $downstream_stop;
1205 } else {
1206 $rnachange->dnStreamSeq($dnstreamseq);
1207 }
1208 # AAChange creation
1209 my $AAobj=$self->RNA->get_Translation;
1210 # storage of prelabel here, to be used in create_mut_objs_after
1211 my $aachange = Bio::Variation::AAChange->new(-start => $RNAprelabel
1212 );
1213 $aachange->isMutation(1);
1214 $aachange->proof('computed');
1215
1216 $seqDiff->add_Variant($aachange);
1217 $self->aachange($aachange);
1218 $rnachange->AAChange($aachange);
1219 $aachange->RNAChange($rnachange);
1220
1221 $aachange->mut_number($self->mutation->issue);
1222 # $before_mutation{aachange}=$aachange;
1223
1224 my $ra_o = Bio::Variation::Allele->new;
1225 $ra_o->seq($dnamut->allele_ori->seq) if $dnamut->allele_ori->seq;
1226 $rnachange->allele_ori($ra_o);
1227
1228 $rnachange->length(CORE::length $rnachange->allele_ori->seq);
1229
1230 my $ra_m = Bio::Variation::Allele->new;
1231 $ra_m->seq($self->mutation->seq) if $self->mutation->seq;
1232 $rnachange->allele_mut($ra_m);
1233 $rnachange->add_Allele($ra_m);
1234
1235 #$rnachange->allele_mut($seq);
1236 $rnachange->end($rnachange->start) if $rnachange->length == 0;
1237
1238 # this holds the aminoacid sequence that will be affected by the mutation
1239 my $aa_allele_ori=$AAobj->labelsubseq($self->mutation->label,undef,
1240 $self->mutation->lastlabel);
1241
1242 my $aa_o = Bio::Variation::Allele->new;
1243 $aa_o->seq($aa_allele_ori) if $aa_allele_ori;
1244 $aachange->allele_ori($aa_o);
1245 #$aachange->allele_ori($aa_allele_ori);
1246
1247 my $aa_length_ori = length($aa_allele_ori);
1248 $aachange->length($aa_length_ori); #print "==========$aa_length_ori\n";
1249 $aachange->end($aachange->start + $aa_length_ori - 1 );
1250 }
1251
1252 =head2 _untranslated
1253
1254 Title : _untranslated
1255 Usage :
1256 Function:
1257
1258 Stores RNA change attributes before mutation
1259 into Bio::Variation::RNAChange object. Links it to
1260 SeqDiff object.
1261
1262 Example :
1263 Returns :
1264 Args : Bio::Variation::SeqDiff object
1265 Bio::Variation::DNAMutation object
1266
1267 See L<Bio::Variation::RNAChange>, L<Bio::Variation::SeqDiff> and
1268 L<Bio::Variation::DNAMutation> for details.
1269
1270 =cut
1271
1272 sub _untranslated {
1273 my ($self, $seqDiff, $dnamut) = @_;
1274 my $rnapos_end;
1275 ($self->mutation->len == 0) ?
1276 ($rnapos_end = $self->mutation->transpos) :
1277 ($rnapos_end = $self->mutation->transpos + $self->mutation->len -1);
1278 my $rnachange = Bio::Variation::RNAChange->new(-start => $self->mutation->transpos,
1279 -end => $rnapos_end
1280 );
1281 #my $rnachange = Bio::Variation::RNAChange->new;
1282
1283 $rnachange->isMutation(1);
1284 my $ra_o = Bio::Variation::Allele->new;
1285 $ra_o->seq($dnamut->allele_ori->seq) if $dnamut->allele_ori->seq;
1286 $rnachange->allele_ori($ra_o);
1287 my $ra_m = Bio::Variation::Allele->new;
1288 $ra_m->seq($dnamut->allele_mut->seq) if $dnamut->allele_mut->seq;
1289 $rnachange->allele_mut($ra_m);
1290 $rnachange->add_Allele($ra_m);
1291 $rnachange->upStreamSeq($dnamut->upStreamSeq);
1292 $rnachange->dnStreamSeq($dnamut->dnStreamSeq);
1293 $rnachange->length($dnamut->length);
1294 $rnachange->mut_number($dnamut->mut_number);
1295 # setting proof
1296 if ($seqDiff->numbering eq "coding") {
1297 $rnachange->proof('experimental');
1298 } else {
1299 $rnachange->proof('computed');
1300 }
1301
1302 my $dist;
1303 if ($rnachange->end < 0) {
1304 $rnachange->region('5\'UTR');
1305 $dnamut->region('5\'UTR');
1306 my $dist = $dnamut->end ;
1307 $dnamut->region_dist($dist);
1308 $dist = $seqDiff->offset - $self->gene->maxtranscript->start + 1 + $dist;
1309 $rnachange->region_dist($dist);
1310 return if $dist < 1; # if mutation is not in mRNA
1311 } else {
1312 $rnachange->region('3\'UTR');
1313 $dnamut->region('3\'UTR');
1314 my $dist = $dnamut->start - $seqDiff->cds_end + $seqDiff->offset;
1315 $dnamut->region_dist($dist);
1316 $dist = $seqDiff->cds_end - $self->gene->maxtranscript->end -1 + $dist;
1317 $rnachange->region_dist($dist);
1318 return if $dist > 0; # if mutation is not in mRNA
1319 }
1320 $seqDiff->add_Variant($rnachange);
1321 $self->rnachange($rnachange);
1322 $rnachange->DNAMutation($dnamut);
1323 $dnamut->RNAChange($rnachange);
1324 }
1325
1326 # args: reference to label changearray, reference to position changearray
1327 # Function: take care of the creation of mutation objects, with
1328 # information AFTER the change takes place
1329 sub _post_mutation {
1330 my ($self, $seqDiff) = @_;
1331
1332 if ($self->rnachange and $self->rnachange->region eq 'coding') {
1333
1334 #$seqDiff->add_Variant($self->rnachange);
1335
1336 my $aachange=$self->aachange;
1337 my ($AAobj,$aa_start_prelabel,$aa_start,$mut_translation);
1338 $AAobj=$self->RNA->get_Translation;
1339 $aa_start_prelabel=$aachange->start;
1340 $aa_start=$AAobj->position($self->RNA->label(2,$aa_start_prelabel));
1341 $aachange->start($aa_start);
1342 $mut_translation=$AAobj->seq;
1343
1344 # this now takes in account possible preinsertions
1345 my $aa_m = Bio::Variation::Allele->new;
1346 $aa_m->seq(substr($mut_translation,$aa_start-1)) if substr($mut_translation,$aa_start-1);
1347 $aachange->allele_mut($aa_m);
1348 $aachange->add_Allele($aa_m);
1349 #$aachange->allele_mut(substr($mut_translation,$aa_start-1));
1350 #$aachange->allele_mut($mut_translation);
1351 my ($rlenori, $rlenmut);
1352 $rlenori = CORE::length($aachange->RNAChange->allele_ori->seq);
1353 $rlenmut = CORE::length($aachange->RNAChange->allele_mut->seq);
1354 #point mutation
1355
1356 if ($rlenori == 1 and $rlenmut == 1 and $aachange->allele_ori->seq ne '*') {
1357 my $alleleseq;
1358 if ($aachange->allele_mut->seq) {
1359 $alleleseq = substr($aachange->allele_mut->seq, 0, 1);
1360 $aachange->allele_mut->seq($alleleseq);
1361 }
1362 $aachange->end($aachange->start);
1363 $aachange->length(1);
1364 }
1365 elsif ( $rlenori == $rlenmut and
1366 $aachange->allele_ori->seq ne '*' ) { #complex inframe mutation
1367 $aachange->allele_mut->seq(substr $aachange->allele_mut->seq,
1368 0,
1369 length($aachange->allele_ori->seq));
1370 }
1371 #inframe mutation
1372 elsif ((int($rlenori-$rlenmut))%3 == 0) {
1373 if ($aachange->RNAChange->allele_mut->seq and
1374 $aachange->RNAChange->allele_ori->seq ) {
1375 # complex
1376 my $rna_len = length ($aachange->RNAChange->allele_mut->seq);
1377 my $len = $rna_len/3;
1378 $len++ unless $rna_len%3 == 0;
1379 $aachange->allele_mut->seq(substr $aachange->allele_mut->seq, 0, $len );
1380 }
1381 elsif ($aachange->RNAChange->codon_pos == 1){
1382 # deletion
1383 if ($aachange->RNAChange->allele_mut->seq eq '') {
1384 $aachange->allele_mut->seq('');
1385 $aachange->end($aachange->start + $aachange->length - 1 );
1386 }
1387 # insertion
1388 elsif ($aachange->RNAChange->allele_ori->seq eq '' ) {
1389 $aachange->allele_mut->seq(substr $aachange->allele_mut->seq, 0,
1390 length ($aachange->RNAChange->allele_mut->seq) / 3);
1391 $aachange->allele_ori->seq('');
1392 $aachange->end($aachange->start + $aachange->length - 1 );
1393 $aachange->length(0);
1394 }
1395 } else {
1396 #elsif ($aachange->RNAChange->codon_pos == 2){
1397 # deletion
1398 if (not $aachange->RNAChange->allele_mut->seq ) {
1399 $aachange->allele_mut->seq(substr $aachange->allele_mut->seq, 0, 1);
1400 }
1401 # insertion
1402 elsif (not $aachange->RNAChange->allele_ori->seq) {
1403 $aachange->allele_mut->seq(substr $aachange->allele_mut->seq, 0,
1404 length ($aachange->RNAChange->allele_mut->seq) / 3 +1);
1405 }
1406 }
1407 } else {
1408 #frameshift
1409 #my $pos = index $aachange->allele_mut
1410 #$aachange->allele_mut(substr($aachange->allele_mut, 0, 1));
1411 $aachange->length(CORE::length($aachange->allele_ori->seq));
1412 my $aaend = $aachange->start + $aachange->length -1;
1413 $aachange->end($aachange->start);
1414 }
1415
1416 # splicing site deletion check
1417 my @beforeexons=@{$self->exons};
1418 my @afterexons=$self->RNA->all_Exons;
1419 my $i;
1420 if (scalar(@beforeexons) ne scalar(@afterexons)) {
1421 my $mut_number = $self->mutation->issue;
1422 $self->warn("Exons have been modified at mutation n.$mut_number!");
1423 $self->rnachange->exons_modified(1);
1424 } else {
1425 EXONCHECK:
1426 foreach $i (0..$#beforeexons) {
1427 if ($beforeexons[$i] ne $afterexons[$i]) {
1428 my $mut_number = $self->mutation->issue;
1429 $self->warn("Exons have been modified at mutation n.$mut_number!");
1430 $self->rnachange->exons_modified(1);
1431 last EXONCHECK;
1432 }
1433 }
1434 }
1435 } else {
1436 #$seqDiff->rnachange(undef);
1437 #print "getting here?";
1438 }
1439 return 1;
1440 }
1441
1442 1;