Mercurial > repos > mahtabm > ensembl

comparison variant_effect_predictor/Bio/Coordinate/GeneMapper.pm @ 0:1f6dce3d34e0

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Uploaded
author mahtabm
date Thu, 11 Apr 2013 02:01:53 -0400
parents
children
comparison
equal deleted inserted replaced
-1:000000000000 0:1f6dce3d34e0
1 # $Id: GeneMapper.pm,v 1.13.2.2 2003/03/13 11:56:30 heikki Exp $
2 #
3 # bioperl module for Bio::Coordinate::GeneMapper
4 #
5 # Cared for by Heikki Lehvaslaiho <heikki@ebi.ac.uk>
6 #
7 # Copyright Heikki Lehvaslaiho
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::Coordinate::GeneMapper - transformations between gene related coordinate systems
16
17 =head1 SYNOPSIS
18
19 use Bio::Coordinate::GeneMapper;
20
21 # get a Bio::RangeI representing the start, end and strand of the CDS
22 # in chromosomal (or entry) coordinates
23 my $cds;
24
25 # get a Bio::Location::Split or an array of Bio::LocationI objects
26 # holding the start, end and strand of all the exons in chromosomal
27 # (or entry) coordinates
28 my $exons;
29
30 # create a gene mapper and set it to map from chromosomal to cds coordinates
31 my $gene = Bio::Coordinate::GeneMapper->new(-in=>'chr',
32 -out=>'cds',
33 -cds=>$cds,
34 -exons=>$exons
35 );
36
37 # get a a Bio::Location or sequence feature in input (chr) coordinates
38 my $loc;
39
40 # map the location into output coordinates and get a new location object
41 $newloc = $gene->map($loc);
42
43
44 =head1 DESCRIPTION
45
46 Bio::Coordinate::GeneMapper is a module for simplifying the mappings
47 of coodinate locations between various gene related locations in human
48 genetics. It also adds a special human genetics twist to coordinate
49 systems by making it possible to disable the use of zero
50 (0). Locations before position one start from -1. See method
51 L<nozero>.
52
53 It understands by name the following coordinate systems and mapping
54 between them:
55
56 peptide (peptide length)
57 ^
58 | -peptide_offset
59 |
60 frame propeptide (propeptide length)
61 ^ ^
62 \ |
63 translate \ |
64 \ |
65 cds (transcript start and end)
66 ^
67 negative_intron | \
68 ^ | \ transcribe
69 \ | \
70 intron exon \
71 ^ ^ ^ /
72 splice \ \ / | /
73 \ \ / | /
74 \ inex | /
75 \ ^ | /
76 \ \ |/
77 ----- gene (gene_length)
78 ^
79 | - gene_offset
80 |
81 chr (or entry)
82
83
84 This structure is kept in the global variable $DAG which is a
85 representation of a Directed Acyclic Graph. The path calculations
86 traversing this graph are done in a helper class. See
87 L<Bio::Coordinate::Graph>.
88
89 Of these, two operations are special cases, translate and splice.
90 Translating and reverse translating are implemented as internal
91 methods that do the simple 1E<lt>-E<gt>3 conversion. Splicing needs
92 additional information that is provided by method L<exons> which takes
93 in an array of Bio::LocationI objects.
94
95 Most of the coordinate system names should be selfexplanatory to
96 anyone familiar with genes. Negative intron coordinate system is
97 starts counting backwards from -1 as the last nucleotide in the
98 intron. This used when only exon and a few flanking intron nucleotides
99 are known.
100
101
102 This class models coordinates within one transcript of a gene, so to
103 tackle multiple transcripts you need several instances of the
104 class. It is therefore valid to argue that the name of the class
105 should be TranscriptMapper. GeneMapper is a catchier name, so it
106 stuck.
107
108
109 =head1 FEEDBACK
110
111 =head2 Mailing Lists
112
113 User feedback is an integral part of the evolution of this and other
114 Bioperl modules. Send your comments and suggestions preferably to the
115 Bioperl mailing lists Your participation is much appreciated.
116
117 bioperl-l@bioperl.org - General discussion
118 http://bio.perl.org/MailList.html - About the mailing lists
119
120 =head2 Reporting Bugs
121
122 Report bugs to the Bioperl bug tracking system to help us keep track
123 the bugs and their resolution. Bug reports can be submitted via email
124 or the web:
125
126 bioperl-bugs@bio.perl.org
127 http://bugzilla.bioperl.org/
128
129 =head1 AUTHOR - Heikki Lehvaslaiho
130
131 Email: heikki@ebi.ac.uk
132 Address:
133
134 EMBL Outstation, European Bioinformatics Institute
135 Wellcome Trust Genome Campus, Hinxton
136 Cambs. CB10 1SD, United Kingdom
137
138 =head1 APPENDIX
139
140 The rest of the documentation details each of the object
141 methods. Internal methods are usually preceded with a _
142
143 =cut
144
145
146 # Let the code begin...
147
148 package Bio::Coordinate::GeneMapper;
149 use vars qw(@ISA %COORDINATE_SYSTEMS %COORDINATE_INTS $TRANSLATION $DAG
150 $NOZERO_VALUES $NOZERO_KEYS);
151 use strict;
152
153 # Object preamble - inherits from Bio::Root::Root
154
155 use Bio::Root::Root;
156 use Bio::Coordinate::Result;
157 use Bio::Location::Simple;
158 use Bio::Coordinate::Graph;
159 use Bio::Coordinate::Collection;
160 use Bio::Coordinate::Pair;
161 use Bio::Coordinate::ExtrapolatingPair;
162
163 @ISA = qw(Bio::Root::Root Bio::Coordinate::MapperI);
164
165 # first set internal values for all translation tables
166
167 %COORDINATE_SYSTEMS = (
168 peptide => 10,
169 propeptide => 9,
170 frame => 8,
171 cds => 7,
172 negative_intron => 6,
173 intron => 5,
174 exon => 4,
175 inex => 3,
176 gene => 2,
177 chr => 1
178 );
179
180 %COORDINATE_INTS = (
181 10 => 'peptide',
182 9 => 'propeptide',
183 8 => 'frame',
184 7 => 'cds',
185 6 => 'negative_intron',
186 5 => 'intron',
187 4 => 'exon',
188 3 => 'inex',
189 2 => 'gene',
190 1 => 'chr'
191 );
192
193 $TRANSLATION = $COORDINATE_SYSTEMS{'cds'}. "-".
194 $COORDINATE_SYSTEMS{'propeptide'};
195
196 $DAG = {
197 10 => [],
198 9 => [10],
199 8 => [],
200 7 => [8, 9],
201 6 => [],
202 5 => [6],
203 4 => [7],
204 3 => [4, 5],
205 2 => [3, 4, 5, 7],
206 1 => [2]
207 };
208
209 $NOZERO_VALUES = {0 => 0, 'in' => 1, 'out' => 2, 'in&out' => 3 };
210 $NOZERO_KEYS = { 0 => 0, 1 => 'in', 2 => 'out', 3 => 'in&out' };
211
212
213 sub new {
214 my($class,@args) = @_;
215 my $self = $class->SUPER::new(@args);
216
217 # prime the graph
218 my $graph = new Bio::Coordinate::Graph;
219 $graph->hash_of_arrays($DAG);
220 $self->graph($graph);
221
222 my($in, $out, $peptide_offset, $exons,
223 $cds, $nozero, $strict) =
224 $self->_rearrange([qw(IN
225 OUT
226 PEPTIDE_OFFSET
227 EXONS
228 CDS
229 NOZERO
230 STRICT
231 )],
232 @args);
233
234 # direction of mapping when going chr to protein
235 $self->{_direction} = 1;
236
237 $in && $self->in($in);
238 $out && $self->out($out);
239 $cds && $self->cds($cds);
240 $exons && ref($exons) =~ /ARRAY/i && $self->exons(@$exons);
241 $peptide_offset && $self->peptide_offset($peptide_offset);
242 $nozero && $self->nozero($nozero);
243 $strict && $self->strict($strict);
244
245 return $self; # success - we hope!
246 }
247
248 =head2 in
249
250 Title : in
251 Usage : $obj->in('peptide');
252 Function: Set and read the input coordinate system.
253 Example :
254 Returns : value of input system
255 Args : new value (optional)
256
257 =cut
258
259 sub in {
260 my ($self,$value) = @_;
261 if( defined $value) {
262 $self->throw("Not a valid input coordinate system name [$value]\n".
263 "Valid values are ". join(", ", keys %COORDINATE_SYSTEMS ))
264 unless defined $COORDINATE_SYSTEMS{$value};
265
266 $self->{'_in'} = $COORDINATE_SYSTEMS{$value};
267 }
268 return $COORDINATE_INTS{ $self->{'_in'} };
269 }
270
271
272 =head2 out
273
274 Title : out
275 Usage : $obj->out('peptide');
276 Function: Set and read the output coordinate system.
277 Example :
278 Returns : value of output system
279 Args : new value (optional)
280
281 =cut
282
283 sub out {
284 my ($self,$value) = @_;
285 if( defined $value) {
286 $self->throw("Not a valid input coordinate system name [$value]\n".
287 "Valid values are ". join(", ", keys %COORDINATE_SYSTEMS ))
288 unless defined $COORDINATE_SYSTEMS{$value};
289
290 $self->{'_out'} = $COORDINATE_SYSTEMS{$value};
291 }
292 return $COORDINATE_INTS{ $self->{'_out'} };
293 }
294
295 =head2 strict
296
297 Title : strict
298 Usage : $obj->strict('peptide');
299 Function: Set and read weather strict boundaried of coordinate
300 systems are enforced.
301 When strict is on, the end of the coordinate range must be defined.
302 Example :
303 Returns : boolean
304 Args : boolean (optional)
305
306 =cut
307
308 sub strict {
309 my ($self,$value) = @_;
310 if( defined $value) {
311 $value ? ( $self->{'_strict'} = 1 ) : ( $self->{'_strict'} = 0 );
312 ## update in each mapper !!
313 }
314 return $self->{'_strict'} || 0 ;
315 }
316
317
318 =head2 nozero
319
320 Title : nozero
321 Usage : $obj->nozero(1);
322 Function: Flag to disable the use of zero in the input,
323 output or both coordinate systems. Use of coordinate
324 systems without zero is a peculiarity common in
325 human genetics community.
326 Example :
327 Returns : 0 (default), or 'in', 'out', 'in&out'
328 Args : 0 (default), or 'in', 'out', 'in&out'
329
330 =cut
331
332 sub nozero {
333 my ($self,$value) = @_;
334
335 if (defined $value) {
336 $self->throw("Not a valid value for nozero [$value]\n".
337 "Valid values are ". join(", ", keys %{$NOZERO_VALUES} ))
338 unless defined $NOZERO_VALUES->{$value};
339 $self->{'_nozero'} = $NOZERO_VALUES->{$value};
340 }
341
342 my $res = $self->{'_nozero'} || 0;
343 return $NOZERO_KEYS->{$res};
344 }
345
346 =head2 graph
347
348 Title : graph
349 Usage : $obj->graph($new_graph);
350 Function: Set and read the graph object representing relationships
351 between coordinate systems
352 Example :
353 Returns : Bio::Coordinate::Graph object
354 Args : new Bio::Coordinate::Graph object (optional)
355
356 =cut
357
358 sub graph {
359 my ($self,$value) = @_;
360 if( defined $value) {
361 $self->throw("Not a valid graph [$value]\n")
362 unless $value->isa('Bio::Coordinate::Graph');
363 $self->{'_graph'} = $value;
364 }
365 return $self->{'_graph'};
366 }
367
368 =head2 peptide
369
370 Title : peptide
371 Usage : $obj->peptide_offset($peptide_coord);
372 Function: Read and write the offset of peptide from the start of propeptide
373 and peptide length
374 Returns : a Bio::Location::Simple object
375 Args : a Bio::LocationI object
376
377 =cut
378
379 sub peptide {
380 my ($self, $value) = @_;
381 if( defined $value) {
382 $self->throw("I need a Bio::LocationI, not [". $value. "]")
383 unless $value->isa('Bio::LocationI');
384
385 $self->throw("Peptide start not defined")
386 unless defined $value->start;
387 $self->{'_peptide_offset'} = $value->start - 1;
388
389 $self->throw("Peptide end not defined")
390 unless defined $value->end;
391 $self->{'_peptide_length'} = $value->end - $self->{'_peptide_offset'};
392
393
394 my $a = $self->_create_pair
395 ('propeptide', 'peptide', $self->strict,
396 $self->{'_peptide_offset'}, $self->{'_peptide_length'} );
397 my $mapper = $COORDINATE_SYSTEMS{'propeptide'}. "-". $COORDINATE_SYSTEMS{'peptide'};
398 $self->{'_mappers'}->{$mapper} = $a;
399 }
400 return Bio::Location::Simple->new
401 (-seq_id => 'propeptide',
402 -start => $self->{'_peptide_offset'} + 1 ,
403 -end => $self->{'_peptide_length'} + $self->{'_peptide_offset'},
404 -strand => 1
405 );
406 }
407
408 =head2 peptide_offset
409
410 Title : peptide_offset
411 Usage : $obj->peptide_offset(20);
412 Function: Set and read the offset of peptide from the start of propeptide
413 Returns : set value or 0
414 Args : new value (optional)
415
416 =cut
417
418 sub peptide_offset {
419 my ($self,$offset, $len) = @_;
420 if( defined $offset) {
421 $self->throw("I need an integer, not [$offset]")
422 unless $offset =~ /^[+-]?\d+$/;
423 $self->{'_peptide_offset'} = $offset;
424
425 if (defined $len) {
426 $self->throw("I need an integer, not [$len]")
427 unless $len =~ /^[+-]?\d+$/;
428 $self->{'_peptide_length'} = $len;
429 }
430
431 my $a = $self->_create_pair
432 ('propeptide', 'peptide', $self->strict, $offset, $self->{'_peptide_length'} );
433 my $mapper = $COORDINATE_SYSTEMS{'propeptide'}. "-". $COORDINATE_SYSTEMS{'peptide'};
434 $self->{'_mappers'}->{$mapper} = $a;
435 }
436 return $self->{'_peptide_offset'} || 0;
437 }
438
439 =head2 peptide_length
440
441 Title : peptide_length
442 Usage : $obj->peptide_length(20);
443 Function: Set and read the offset of peptide from the start of propeptide
444 Returns : set value or 0
445 Args : new value (optional)
446
447 =cut
448
449
450 sub peptide_length {
451 my ($self, $len) = @_;
452 if( defined $len) {
453 $self->throw("I need an integer, not [$len]")
454 if defined $len && $len !~ /^[+-]?\d+$/;
455 $self->{'_peptide_length'} = $len;
456 }
457 return $self->{'_peptide_length'};
458 }
459
460
461 =head2 exons
462
463 Title : exons
464 Usage : $obj->exons(@exons);
465 Function: Set and read the offset of CDS from the start of transcipt
466 You do not have to sort the exons before calling this method as
467 they will be sorted automatically.
468 If you have not defined the CDS, is will be set to span all
469 exons here.
470 Returns : array of Bio::LocationI exons in genome coordinates or 0
471 Args : array of Bio::LocationI exons in genome (or entry) coordinates
472
473 =cut
474
475 sub exons {
476 my ($self,@value) = @_;
477 my $cds_mapper = $COORDINATE_SYSTEMS{'gene'}. "-". $COORDINATE_SYSTEMS{'cds'};
478 my $inex_mapper =
479 $COORDINATE_SYSTEMS{'gene'}. "-". $COORDINATE_SYSTEMS{'inex'};
480 my $exon_mapper =
481 $COORDINATE_SYSTEMS{'gene'}. "-". $COORDINATE_SYSTEMS{'exon'};
482 my $intron_mapper =
483 $COORDINATE_SYSTEMS{'gene'}. "-". $COORDINATE_SYSTEMS{'intron'};
484 my $negative_intron_mapper =
485 $COORDINATE_SYSTEMS{'intron'}. "-". $COORDINATE_SYSTEMS{'negative_intron'};
486 my $exon_cds_mapper = $COORDINATE_SYSTEMS{'exon'}. "-". $COORDINATE_SYSTEMS{'cds'};
487
488 if(@value) {
489 if (ref($value[0]) &&
490 $value[0]->isa('Bio::SeqFeatureI') and
491 $value[0]->location->isa('Bio::Location::SplitLocationI')) {
492 @value = $value[0]->location->each_Location;
493 } else {
494 $self->throw("I need an array , not [@value]")
495 unless ref \@value eq 'ARRAY';
496 $self->throw("I need a reference to an array of Bio::LocationIs, not to [".
497 $value[0]. "]")
498 unless ref $value[0] and $value[0]->isa('Bio::LocationI');
499 }
500
501 #
502 # sort the input array
503 #
504 # and if the used has not defined CDS assume it is the complete exonic range
505 if (defined $value[0]->strand && $value[0]->strand == - 1) { #reverse strand
506 @value = map { $_->[0] }
507 sort { $b->[1] <=> $a->[1] }
508 map { [ $_, $_->start] }
509 @value;
510
511 unless ($self->cds) {
512 $self->cds(new Bio::Location::Simple(-start => $value[-1]->start,
513 -end => $value[0]->end,
514 -strand=> $value[0]->strand,
515 -seq_id=> $value[0]->seq_id,
516 )
517 );
518 }
519 } else { #undef or forward strand
520 @value = map { $_->[0] }
521 sort { $a->[1] <=> $b->[1] }
522 map { [ $_, $_->start] }
523 @value;
524 unless ($self->cds) {
525 $self->cds(new Bio::Location::Simple(-start => $value[0]->start,
526 -end => $value[-1]->end,
527 -strand=> $value[0]->strand,
528 -seq_id=> $value[0]->seq_id,
529 )
530 );
531 }
532
533 }
534
535 $self->{'_chr_exons'} = \@value;
536
537 # transform exons from chromosome to gene coordinates
538 # but only if gene coordinate system has been set
539 my @exons ;
540 #my $gene_mapper = $self->$COORDINATE_SYSTEMS{'chr'}. "-". $COORDINATE_SYSTEMS{'gene'};
541 my $gene_mapper = "1-2";
542 if (defined $self->{'_mappers'}->{$gene_mapper} ) {
543
544 my $tmp_in = $self->{'_in'};
545 my $tmp_out = $self->{'_out'};
546 my $tmp_verb = $self->verbose;
547 $self->verbose(0);
548
549 $self->in('chr');
550 $self->out('gene');
551
552 @exons = map {$self->map($_)} @value;
553
554 $self->{'_in'} = ($tmp_in);
555 $self->{'_out'} = ($tmp_out);
556 $self->verbose($tmp_verb);
557 } else {
558 @exons = @value;
559 }
560
561 my $cds_map = Bio::Coordinate::Collection->new;
562 my $inex_map = Bio::Coordinate::Collection->new;
563 my $exon_map = Bio::Coordinate::Collection->new;
564 my $exon_cds_map = Bio::Coordinate::Collection->new;
565 my $intron_map = Bio::Coordinate::Collection->new;
566 my $negative_intron_map = Bio::Coordinate::Collection->new;
567
568 my $tr_end = 0;
569 my $coffset;
570 my $exon_counter;
571 my $prev_exon_end;
572
573 for my $exon ( @exons ) {
574
575 $exon_counter++;
576
577 #
578 # gene -> cds
579 #
580
581 my $match1 = Bio::Location::Simple->new
582 (-seq_id =>'gene' ,
583 -start => $exon->start,
584 -end => $exon->end, -strand=>1 );
585 my $match2 = Bio::Location::Simple->new
586 (-seq_id => 'cds',
587 -start => $tr_end + 1,
588 -end => $tr_end + $exon->end - $exon->start +1,
589 -strand=>$exon->strand );
590
591 $cds_map->add_mapper(Bio::Coordinate::Pair->new
592 (-in => $match1,
593 -out => $match2,
594 )
595 );
596
597 if ($exon->start <= 1 and $exon->end >= 1) {
598 $coffset = $tr_end - $exon->start + 1;
599 }
600 $tr_end = $tr_end + $exon->end - $exon->start + 1;
601
602 #
603 # gene -> intron
604 #
605
606 if (defined $prev_exon_end) {
607 my $match3 = Bio::Location::Simple->new
608 (-seq_id =>'gene',
609 -start => $prev_exon_end + 1,
610 -end => $exon->start -1, -strand=>$exon->strand );
611
612 my $match4 = Bio::Location::Simple->new
613 (-seq_id => 'intron'. ($exon_counter -1),
614 -start => 1,
615 -end => $exon->start - 1 - $prev_exon_end,
616 -strand=>$exon->strand );
617
618 # negative intron coordinates
619 my $match5 = Bio::Location::Simple->new
620 (-seq_id => 'intron'. ($exon_counter -1),
621 -start => -1 * ($exon->start - 2 - $prev_exon_end) -1,
622 -end => -1,
623 -strand=>$exon->strand );
624
625 $inex_map->add_mapper(Bio::Coordinate::Pair->new
626 (-in => $match3,
627 -out => $match4
628 )
629 );
630 $intron_map->add_mapper(Bio::Coordinate::Pair->new
631 (-in => $self->_clone_loc($match3),
632 -out => $self->_clone_loc($match4)
633 )
634 );
635 $negative_intron_map->add_mapper(Bio::Coordinate::Pair->new
636 (-in => $self->_clone_loc($match4),
637 -out => $match5
638 ));
639
640 }
641
642 # store the value
643 $prev_exon_end = $exon->end;
644
645 #
646 # gene -> exon
647 #
648 my $match6 = Bio::Location::Simple->new
649 (-seq_id => 'exon'. $exon_counter,
650 -start => 1,
651 -end => $exon->end - $exon->start +1,
652 -strand=> $exon->strand );
653
654 my $pair2 = Bio::Coordinate::Pair->new(-in => $self->_clone_loc($match1),
655 -out => $match6
656 );
657 my $pair3 = Bio::Coordinate::Pair->new(-in => $self->_clone_loc($match6),
658 -out => $self->_clone_loc($match2)
659 );
660 $inex_map->add_mapper(Bio::Coordinate::Pair->new
661 (-in => $self->_clone_loc($match1),
662 -out => $match6
663 )
664 );
665 $exon_map->add_mapper(Bio::Coordinate::Pair->new
666 (-in => $self->_clone_loc($match1),
667 -out => $self->_clone_loc($match6)
668 )
669 );
670 $exon_cds_map->add_mapper(Bio::Coordinate::Pair->new
671 (-in => $self->_clone_loc($match6),
672 -out => $self->_clone_loc($match2)
673 )
674 );
675
676 }
677
678 # move coordinate start if exons have negative values
679 if ($coffset) {
680 foreach my $m ($cds_map->each_mapper) {
681 $m->out->start($m->out->start - $coffset);
682 $m->out->end($m->out->end - $coffset);
683 }
684
685 }
686
687 $self->{'_mappers'}->{$cds_mapper} = $cds_map;
688 $self->{'_mappers'}->{$exon_cds_mapper} = $exon_cds_map;
689 $self->{'_mappers'}->{$inex_mapper} = $inex_map;
690 $self->{'_mappers'}->{$exon_mapper} = $exon_map;
691 $self->{'_mappers'}->{$intron_mapper} = $intron_map;
692 $self->{'_mappers'}->{$negative_intron_mapper} = $negative_intron_map;
693 }
694 return @{$self->{'_chr_exons'}} || 0;
695 }
696
697 =head2 _clone_loc
698
699 Title : _clone_loc
700 Usage : $copy_of_loc = $obj->_clone_loc($loc);
701 Function: Make a deep copy of a simple location
702 Returns : a Bio::Location::Simple object
703 Args : a Bio::Location::Simple object to be cloned
704
705 =cut
706
707
708 sub _clone_loc { # clone a simple location
709 my ($self,$loc) = @_;
710
711 $self->throw("I need a Bio::Location::Simple , not [". ref $loc. "]")
712 unless $loc->isa('Bio::Location::Simple');
713
714 return Bio::Location::Simple->new
715 (-seq_id => $loc->seq_id,
716 -start => $loc->start,
717 -end => $loc->end,
718 -strand=> $loc->strand,
719 -location_type => $loc->location_type
720 );
721 }
722
723
724 =head2 cds
725
726 Title : cds
727 Usage : $obj->cds(20);
728 Function: Set and read the offset of CDS from the start of transcipt
729
730 Simple input can be an integer which gives the start of the
731 coding region in genomic coordinate. If you want to provide
732 the end of the coding region or indicate the use of the
733 opposite strand, you have to pass a Bio::RangeI
734 (e.g. Bio::Location::Simple or Bio::SegFeature::Generic)
735 object to this method.
736
737 Returns : set value or 0
738 Args : new value (optional)
739
740 =cut
741
742 sub cds {
743 my ($self,$value) = @_;
744 if( defined $value) {
745 if ($value =~ /^[+-]?\d+$/ ) {
746 my $loc = Bio::Location::Simple->new(-start=>$value);
747 $self->{'_cds'} = $loc;
748 }
749 elsif (ref $value && $value->isa('Bio::RangeI') ) {
750 $self->{'_cds'} = $value;
751 } else {
752 $self->throw("I need an integer or Bio::RangeI, not [$value]")
753 }
754 # strand !!
755 my $len;
756
757 $len = $self->{'_cds'}->end - $self->{'_cds'}->start +1
758 if defined $self->{'_cds'}->end;
759
760 my $a = $self->_create_pair
761 ('chr', 'gene', 0,
762 $self->{'_cds'}->start-1,
763 $len,
764 $self->{'_cds'}->strand);
765 my $mapper = $COORDINATE_SYSTEMS{'chr'}. "-". $COORDINATE_SYSTEMS{'gene'};
766 $self->{'_mappers'}->{$mapper} = $a;
767
768 # recalculate exon-based mappers
769 if ( defined $self->{'_chr_exons'} ) {
770 $self->exons(@{$self->{'_chr_exons'}});
771 }
772
773 }
774 return $self->{'_cds'} || 0;
775 }
776
777
778 =head2 map
779
780 Title : map
781 Usage : $newpos = $obj->map(5);
782 Function: Map the location from the input coordinate system
783 to a new value in the output coordinate system.
784 Example :
785 Returns : new value in the output coordiante system
786 Args : a Bio::Location::Simple
787
788 =cut
789
790 sub map {
791 my ($self,$value) = @_;
792 my ($res);
793
794 $self->throw("Need to pass me a Bio::Location::Simple or ".
795 "Bio::Location::Simple or Bio::SeqFeatureI, not [".
796 ref($value). "]")
797 unless ref($value) && ($value->isa('Bio::Location::Simple') or
798 $value->isa('Bio::Location::SplitLocationI') or
799 $value->isa('Bio::SeqFeatureI'));
800 $self->throw("Input coordinate system not set")
801 unless $self->{'_in'};
802 $self->throw("Output coordinate system not set")
803 unless $self->{'_out'};
804 $self->throw("Do not be silly. Input and output coordinate ".
805 "systems are the same!")
806 unless $self->{'_in'} != $self->{'_out'};
807
808 $self->_check_direction();
809
810 $value = $value->location if $value->isa('Bio::SeqFeatureI');
811 print STDERR "=== Start location: ". $value->start. ",".
812 $value->end. " (". $value->strand. ")\n" if $self->verbose > 0;
813
814
815 # if nozero coordinate system is used in the input values
816 if ( defined $self->{'_nozero'} &&
817 ( $self->{'_nozero'} == 1 || $self->{'_nozero'} == 3 ) ) {
818 $value->start($value->start + 1)
819 if defined $value->start && $value->start < 1;
820 $value->end($value->end + 1)
821 if defined $value->end && $value->end < 1;
822 }
823
824 my @steps = $self->_get_path();
825 print "mapping ", $self->{'_in'}, "->", $self->{'_out'},
826 " Mappers: ", join(", ", @steps), "\n" if $self->verbose > 0;
827
828 foreach my $mapper (@steps) {
829 if ($mapper eq $TRANSLATION) {
830 if ($self->direction == 1) {
831 $value = $self->_translate($value);
832 print STDERR "+ $TRANSLATION cds -> propeptide (translate) \n"
833 if $self->verbose > 0;
834 } else {
835 $value = $self->_reverse_translate($value);
836 print STDERR "+ $TRANSLATION propeptide -> cds (reverse translate) \n"
837 if $self->verbose > 0;
838 }
839 }
840 # keep the start and end values, and go on to next iteration
841 # if this mapper is not set
842 elsif ( ! defined $self->{'_mappers'}->{$mapper} ) {
843 # update mapper name
844 $mapper =~ /\d+-(\d+)/; my ($counter) = $1;
845 $value->seq_id($COORDINATE_INTS{$counter});
846 print STDERR "- $mapper\n" if $self->verbose > 0;
847 } else {
848
849 #
850 # the DEFAULT : generic mapping
851 #
852 $value = $self->{'_mappers'}->{$mapper}->map($value);
853 $value->purge_gaps
854 if ($value && $value->isa('Bio::Location::SplitLocationI') && $value->can('gap'));
855 print STDERR "+ $mapper (", $self->direction, "): start ",
856 $value->start, " end ", $value->end, "\n"
857 if $value && $self->verbose > 0;
858 }
859 }
860
861 # if nozero coordinate system is asked to be used in the output values
862 if ( defined $value && defined $self->{'_nozero'} &&
863 ( $self->{'_nozero'} == 2 || $self->{'_nozero'} == 3 ) ) {
864
865 $value->start($value->start - 1)
866 if defined $value->start && $value->start < 1;
867 $value->end($value->end - 1)
868 if defined $value->end && $value->end < 1;
869 }
870
871 # handle merging of adjacent split locations!
872
873 if (ref $value eq "Bio::Coordinate::Result" && $value->each_match > 1 ) {
874 my $prevloc;
875 my $merging = 0;
876 my $newvalue;
877 my @matches;
878 foreach my $loc ( $value->each_Location(1) ) {
879 unless ($prevloc) {
880 $prevloc = $loc;
881 push @matches, $prevloc;
882 next;
883 }
884 if ($prevloc->end == ($loc->start - 1) && $prevloc->seq_id eq $loc->seq_id) {
885 $prevloc->end($loc->end);
886 $merging = 1;
887 } else {
888 push @matches, $loc;
889 $prevloc = $loc;
890 }
891 }
892 if ($merging) {
893 if (@matches > 1 ) {
894 $newvalue = Bio::Coordinate::Result->new;
895 map {$newvalue->add_sub_Location} @matches;
896 } else {
897 $newvalue = Bio::Coordinate::Result::Match->new
898 (-seq_id => $matches[0]->seq_id,
899 -start => $matches[0]->start,
900 -end => $matches[0]->end,
901 -strand=> $matches[0]->strand );
902 }
903 $value = $newvalue;
904 }
905 }
906 elsif (ref $value eq "Bio::Coordinate::Result" && $value->each_match == 1 ){
907 $value = $value->match;
908 }
909
910
911 return $value;
912 }
913
914 =head2 direction
915
916 Title : direction
917 Usage : $obj->direction('peptide');
918 Function: Read-only method for the direction of mapping deduced from
919 predefined input and output coordinate names.
920 Example :
921 Returns : 1 or -1, mapping direction
922 Args : new value (optional)
923
924 =cut
925
926 sub direction {
927 my ($self) = @_;
928 return $self->{'_direction'};
929 }
930
931
932 =head2 swap
933
934 Title : swap
935 Usage : $obj->swap;
936 Function: Swap the direction of transformation
937 (input <-> output)
938 Example :
939 Returns : 1
940 Args :
941
942 =cut
943
944 sub swap {
945 my ($self,$value) = @_;
946
947 ($self->{'_in'}, $self->{'_out'}) = ($self->{'_out'}, $self->{'_in'});
948 map { $self->{'_mappers'}->{$_}->swap } keys %{$self->{'_mappers'}};
949
950 # record the changed direction;
951 $self->{_direction} *= -1;
952
953 return 1;
954 }
955
956
957 =head2 to_string
958
959 Title : to_string
960 Usage : $newpos = $obj->to_string(5);
961 Function: Dump the internal mapper values into a human readable format
962 Example :
963 Returns : string
964 Args :
965
966 =cut
967
968 sub to_string {
969 my ($self) = shift;
970
971 print "-" x 40, "\n";
972
973 # chr-gene
974 my $mapper_str = 'chr-gene';
975 my $mapper = $self->_mapper_string2code($mapper_str);
976
977 printf "\n %-12s (%s)\n", $mapper_str, $mapper ;
978 if (defined $self->cds) {
979 my $end = $self->cds->end -1 if defined $self->cds->end;
980 printf "%16s%s: %s (%s)\n", ' ', 'gene offset', $self->cds->start-1 , $end || '';
981 printf "%16s%s: %s\n", ' ', 'gene strand', $self->cds->strand || 0;
982 }
983
984 # gene-intron
985 $mapper_str = 'gene-intron';
986 $mapper = $self->_mapper_string2code($mapper_str);
987 printf "\n %-12s (%s)\n", $mapper_str, $mapper ;
988
989 my $i = 1;
990 foreach my $pair ( $self->{'_mappers'}->{$mapper}->each_mapper ) {
991 printf "%8s :%8s -> %-12s\n", $i, $pair->in->start, $pair->out->start ;
992 printf "%8s :%8s -> %-12s\n", '', $pair->in->end, $pair->out->end ;
993 $i++;
994 }
995
996 # intron-negative_intron
997 $mapper_str = 'intron-negative_intron';
998 $mapper = $self->_mapper_string2code($mapper_str);
999 printf "\n %-12s (%s)\n", $mapper_str, $mapper ;
1000
1001 $i = 1;
1002 foreach my $pair ( $self->{'_mappers'}->{$mapper}->each_mapper ) {
1003 printf "%8s :%8s -> %-12s\n", $i, $pair->in->start, $pair->out->start ;
1004 printf "%8s :%8s -> %-12s\n", '', $pair->in->end, $pair->out->end ;
1005 $i++;
1006 }
1007
1008
1009 # gene-exon
1010 $mapper_str = 'gene-exon';
1011 $mapper = $self->_mapper_string2code($mapper_str);
1012 printf "\n %-12s (%s)\n", $mapper_str, $mapper ;
1013
1014 $i = 1;
1015 foreach my $pair ( $self->{'_mappers'}->{$mapper}->each_mapper ) {
1016 printf "%8s :%8s -> %-12s\n", $i, $pair->in->start, $pair->out->start ;
1017 printf "%8s :%8s -> %-12s\n", '', $pair->in->end, $pair->out->end ;
1018 $i++;
1019 }
1020
1021
1022 # gene-cds
1023 $mapper_str = 'gene-cds';
1024 $mapper = $self->_mapper_string2code($mapper_str);
1025 printf "\n %-12s (%s)\n", $mapper_str, $mapper ;
1026
1027 $i = 1;
1028 foreach my $pair ( $self->{'_mappers'}->{$mapper}->each_mapper ) {
1029 printf "%8s :%8s -> %-12s\n", $i, $pair->in->start, $pair->out->start ;
1030 printf "%8s :%8s -> %-12s\n", '', $pair->in->end, $pair->out->end ;
1031 $i++;
1032 }
1033
1034 # cds-propeptide
1035 $mapper_str = 'cds-propeptide';
1036 $mapper = $self->_mapper_string2code($mapper_str);
1037 printf "\n %-12s (%s)\n", $mapper_str, $mapper ;
1038 printf "%9s%-12s\n", "", '"translate"';
1039
1040
1041 # propeptide-peptide
1042 $mapper_str = 'propeptide-peptide';
1043 $mapper = $self->_mapper_string2code($mapper_str);
1044 printf "\n %-12s (%s)\n", $mapper_str, $mapper ;
1045 printf "%16s%s: %s\n", ' ', "peptide offset", $self->peptide_offset;
1046
1047
1048
1049 print "\nin : ", $self->in, "\n";
1050 print "out: ", $self->out, "\n";
1051 my $dir;
1052 $self->direction ? ($dir='forward') : ($dir='reverse');
1053 printf "direction: %-8s(%s)\n", $dir, $self->direction;
1054 print "\n", "-" x 40, "\n";
1055
1056 1;
1057 }
1058
1059 sub _mapper_code2string {
1060 my ($self, $code) = @_;
1061 my ($a, $b) = $code =~ /(\d+)-(\d+)/;
1062 return $COORDINATE_INTS{$a}. '-'. $COORDINATE_INTS{$b};
1063
1064 }
1065
1066 sub _mapper_string2code {
1067 my ($self, $string) =@_;
1068 my ($a, $b) = $string =~ /([^-]+)-(.*)/;
1069 return $COORDINATE_SYSTEMS{$a}. '-'. $COORDINATE_SYSTEMS{$b};
1070 }
1071
1072
1073 =head2 _create_pair
1074
1075 Title : _create_pair
1076 Usage : $mapper = $obj->_create_pair('chr', 'gene', 0, 2555, 10000, -1);
1077 Function: Internal helper method to create a mapper between
1078 two coordinate systems
1079 Returns : a Bio::Coordinate::Pair object
1080 Args : string, input coordinate system name,
1081 string, output coordinate system name,
1082 boolean, strict mapping
1083 positive integer, offset
1084 positive integer, length
1085 1 || -1 , strand
1086
1087 =cut
1088
1089 sub _create_pair {
1090 my ($self, $in, $out, $strict, $offset, $length, $strand ) = @_;
1091 $strict ||= 0;
1092 $strand ||= 1;
1093 $length ||= 20;
1094
1095 my $match1 = Bio::Location::Simple->new
1096 (-seq_id => $in,
1097 -start => $offset+1,
1098 -end => $offset+$length, -strand=>1 );
1099
1100 my $match2 = Bio::Location::Simple->new
1101 (-seq_id => $out,
1102 -start => 1,
1103 -end => $length, -strand=>$strand );
1104
1105 my $pair = Bio::Coordinate::ExtrapolatingPair->new
1106 (-in => $match1,
1107 -out => $match2,
1108 -strict => $strict
1109 );
1110
1111 return $pair;
1112
1113 }
1114
1115
1116 =head2 _translate
1117
1118 Title : _translate
1119 Usage : $newpos = $obj->_translate($loc);
1120 Function: Translate the location from the CDS coordinate system
1121 to a new value in the propeptide coordinate system.
1122 Example :
1123 Returns : new location
1124 Args : a Bio::Location::Simple or Bio::Location::SplitLocationI
1125
1126 =cut
1127
1128 sub _translate {
1129 my ($self,$value) = @_;
1130
1131 $self->throw("Need to pass me a Bio::Location::Simple or ".
1132 "Bio::Location::SplitLocationI, not [". ref($value). "]")
1133 unless defined $value &&
1134 ($value->isa('Bio::Location::Simple') || $value->isa('Bio::Location::SplitLocationI'));
1135
1136 my $seqid = 'propeptide';
1137
1138 if ($value->isa("Bio::Location::SplitLocationI")) {
1139 my $split = new Bio::Location::Split(-seq_id=>$seqid);
1140 foreach my $loc ( $value->each_Location(1) ) {
1141
1142 my $match = new Bio::Location::Simple(-start => int($loc->start / 3 )+1,
1143 -end => int($loc->end / 3 )+1,
1144 -seq_id => $seqid,
1145 -strand => 1
1146 );
1147 $split->add_sub_Location($match);
1148 }
1149 return $split;
1150
1151 } else {
1152 return new Bio::Location::Simple(-start => int($value->start / 3 )+1,
1153 -end => int($value->end / 3 )+1,
1154 -seq_id => $seqid,
1155 -strand => 1
1156 );
1157 }
1158 }
1159
1160 sub _frame {
1161 my ($self,$value) = @_;
1162
1163 $self->throw("Need to pass me a Bio::Location::Simple or ".
1164 "Bio::Location::SplitLocationI, not [". ref($value). "]")
1165 unless defined $value &&
1166 ($value->isa('Bio::Location::Simple') || $value->isa('Bio::Location::SplitLocationI'));
1167
1168 my $seqid = 'propeptide';
1169
1170 if ($value->isa("Bio::Location::SplitLocationI")) {
1171 my $split = new Bio::Location::Split(-seq_id=>$seqid);
1172 foreach my $loc ( $value->each_Location(1) ) {
1173
1174 my $match = new Bio::Location::Simple(-start => ($value->start-1) % 3 +1,
1175 -end => ($value->end-1) % 3 +1,
1176 -seq_id => 'frame',
1177 -strand => 1
1178 );
1179 $split->add_sub_Location($match);
1180 }
1181 return $split;
1182 } else {
1183 return new Bio::Location::Simple(-start => ($value->start-1) % 3 +1,
1184 -end => ($value->end-1) % 3 +1,
1185 -seq_id => 'frame',
1186 -strand => 1
1187 );
1188 }
1189 }
1190
1191
1192 =head2 _reverse_translate
1193
1194 Title : _reverse_translate
1195 Usage : $newpos = $obj->_reverse_translate(5);
1196 Function: Reverse translate the location from the propeptide
1197 coordinate system to a new value in the CSD.
1198 Note that a single peptide location expands to cover
1199 the codon triplet
1200 Example :
1201 Returns : new location in the CDS coordinate system
1202 Args : a Bio::Location::Simple or Bio::Location::SplitLocationI
1203
1204 =cut
1205
1206 sub _reverse_translate {
1207 my ($self,$value) = @_;
1208
1209
1210 $self->throw("Need to pass me a Bio::Location::Simple or ".
1211 "Bio::Location::SplitLocationI, not [". ref($value). "]")
1212 unless defined $value &&
1213 ($value->isa('Bio::Location::Simple') || $value->isa('Bio::Location::SplitLocationI'));
1214
1215 my $seqid = 'cds';
1216
1217 if ($value->isa("Bio::Location::SplitLocationI")) {
1218 my $split = new Bio::Location::Split(-seq_id=>$seqid);
1219 foreach my $loc ( $value->each_Location(1) ) {
1220
1221 my $match = new Bio::Location::Simple(-start => $value->start * 3 - 2,
1222 -end => $value->end * 3,
1223 -seq_id => $seqid,
1224 -strand => 1
1225 );
1226 $split->add_sub_Location($match);
1227 }
1228 return $split;
1229
1230 } else {
1231 return new Bio::Location::Simple(-start => $value->start * 3 - 2,
1232 -end => $value->end * 3,
1233 -seq_id => $seqid,
1234 -strand => 1
1235 );
1236 }
1237 }
1238
1239
1240 =head2 _check_direction
1241
1242 Title : _check_direction
1243 Usage : $obj->_check_direction();
1244 Function: Check and swap when needed the direction the location
1245 mapping Pairs based on input and output values
1246 Example :
1247 Returns : new location
1248 Args : a Bio::Location::Simple
1249
1250 =cut
1251
1252 sub _check_direction {
1253 my ($self) = @_;
1254
1255 my $new_direction = 1;
1256 $new_direction = -1 if $self->{'_in'} > $self->{'_out'};
1257
1258 unless ($new_direction == $self->{_direction} ) {
1259 map { $self->{'_mappers'}->{$_}->swap } keys %{$self->{'_mappers'}};
1260 # record the changed direction;
1261 $self->{_direction} *= -1;
1262 }
1263 1;
1264 }
1265
1266
1267 =head2 _get_path
1268
1269 Title : _get_path
1270 Usage : $obj->_get_path('peptide');
1271 Function: internal method for finding that shortest path between
1272 input and output coordinate systems.
1273 Calculations and caching are handled by the graph class.
1274 See L<Bio::Coordinate::Graph>.
1275 Example :
1276 Returns : array of the mappers
1277 Args : none
1278
1279 =cut
1280
1281 sub _get_path {
1282 my ($self) = @_;
1283
1284 my $start = $self->{'_in'} || 0;
1285 my $end = $self->{'_out'} || 0;
1286
1287 # note the order
1288 # always go from smaller to bigger: it makes caching more efficient
1289 my $reverse;
1290 if ($start > $end) {
1291 ($start, $end) = ($end, $start );
1292 $reverse++;
1293 }
1294
1295 my @mappers;
1296 if (exists $self->{'_previous_path'} and
1297 $self->{'_previous_path'} eq "$start$end" ) {
1298 # use cache
1299 @mappers = @{$self->{'_mapper_path'}};
1300 } else {
1301 my $mapper;
1302 my $prev_node = '';
1303 @mappers =
1304 map { $mapper = "$prev_node-$_"; $prev_node = $_; $mapper; }
1305 $self->{'_graph'}->shortest_path($start, $end);
1306 shift @mappers;
1307
1308 $self->{'_previous_path'} = "$start$end";
1309 $self->{'_mapper_path'} = \@mappers;
1310 }
1311
1312 $reverse ? return reverse @mappers : return @mappers;
1313 }
1314
1315
1316 1;
1317