Bio/Map/PositionI.pm

   1 #
   2 # BioPerl module for Bio::Map::PositionI
   3 #
   4 # Please direct questions and support issues to <bioperl-l@bioperl.org>
   5 #
   6 # Cared for by Sendu Bala <bix@sendu.me.uk>
   7 #
   8 # Copyright Jason Stajich
   9 #
  10 # You may distribute this module under the same terms as perl itself
  11
  12 # POD documentation - main docs before the code
  13
  14 =head1 NAME
  15
  16 Bio::Map::PositionI - Abstracts the notion of a position having a value in the context of a marker and a Map
  17
  18 =head1 SYNOPSIS
  19
  20     # do not use this module directly
  21     # See Bio::Map::Position for an example of
  22     # implementation.
  23
  24 =head1 DESCRIPTION
  25
  26 This object stores one of the positions that a mappable object
  27 (e.g. Marker) may have in a map.
  28
  29 Positions can have non-numeric values or other methods to store the locations,
  30 so they have a method numeric() which does the conversion. numeric()
  31 returns the position in a form that can be compared between other positions of
  32 the same type. It is not necessarily a value suitable for sorting positions (it
  33 may be the distance from the previous position); for that purpose the result of
  34 sortable() should be used.
  35
  36 A 'position', in addition to being a single point, can also be an area and so
  37 can be imagined as a range and compared with other positions on the basis of
  38 overlap, intersection etc.
  39
  40 =head1 FEEDBACK
  41
  42 =head2 Mailing Lists
  43
  44 User feedback is an integral part of the evolution of this and other
  45 Bioperl modules. Send your comments and suggestions preferably to
  46 the Bioperl mailing list.  Your participation is much appreciated.
  47
  48   bioperl-l@bioperl.org                  - General discussion
  49   http://bioperl.org/wiki/Mailing_lists  - About the mailing lists
  50
  51 =head2 Support
  52
  53 Please direct usage questions or support issues to the mailing list:
  54
  55 I<bioperl-l@bioperl.org>
  56
  57 rather than to the module maintainer directly. Many experienced and
  58 reponsive experts will be able look at the problem and quickly
  59 address it. Please include a thorough description of the problem
  60 with code and data examples if at all possible.
  61
  62 =head2 Reporting Bugs
  63
  64 Report bugs to the Bioperl bug tracking system to help us keep track
  65 of the bugs and their resolution. Bug reports can be submitted via the
  66 web:
  67
  68   https://github.com/bioperl/bioperl-live/issues
  69
  70 =head1 AUTHOR - Jason Stajich
  71
  72 Email jason-at-bioperl.org
  73
  74 =head1 CONTRIBUTORS
  75
  76 Lincoln Stein, lstein-at-cshl.org
  77 Heikki Lehvaslaiho, heikki-at-bioperl-dot-org
  78 Sendu Bala, bix@sendu.me.uk
  79
  80 =head1 APPENDIX
  81
  82 The rest of the documentation details each of the object methods.
  83 Internal methods are usually preceded with a _
  84
  85 =cut
  86
  87 # Let the code begin...
  88
  89 package Bio::Map::PositionI;
  90 use strict;
  91 use Bio::Map::PositionHandler;
  92 use Bio::Map::Mappable;
  93 use Scalar::Util qw(looks_like_number);
  94
  95 use base qw(Bio::Map::EntityI Bio::RangeI);
  96
  97 =head2 EntityI methods
  98
  99  These are fundamental to coordination of Positions and other entities, so are
 100  implemented at the interface level
 101
 102 =cut
 103
 104 =head2 get_position_handler
 105
 106  Title   : get_position_handler
 107  Usage   : my $position_handler = $entity->get_position_handler();
 108  Function: Gets a PositionHandlerI that $entity is registered with.
 109  Returns : Bio::Map::PositionHandlerI object
 110  Args    : none
 111
 112 =cut
 113
 114 sub get_position_handler {
 115     my $self = shift;
 116     unless (defined $self->{_eh}) {
 117         my $ph = Bio::Map::PositionHandler->new(-self => $self);
 118         $self->{_eh} = $ph;
 119         $ph->register;
 120     }
 121     return $self->{_eh};
 122 }
 123
 124 =head2 PositionHandlerI-related methods
 125
 126  These are fundamental to coordination of Positions and other entities, so are
 127  implemented at the interface level
 128
 129 =cut
 130
 131 =head2 map
 132
 133  Title   : map
 134  Usage   : my $map = $position->map();
 135            $position->map($map);
 136  Function: Get/Set the map the position is in.
 137  Returns : L<Bio::Map::MapI>
 138  Args    : none to get
 139            new L<Bio::Map::MapI> to set
 140
 141 =cut
 142
 143 sub map {
 144     my ($self, $map) = @_;
 145     return $self->get_position_handler->map($map);
 146 }
 147
 148 =head2 element
 149
 150  Title   : element
 151  Usage   : my $element = $position->element();
 152            $position->element($element);
 153  Function: Get/Set the element the position is for.
 154  Returns : L<Bio::Map::MappableI>
 155  Args    : none to get
 156            new L<Bio::Map::MappableI> to set
 157
 158 =cut
 159
 160 sub element {
 161     my ($self, $element) = @_;
 162     return $self->get_position_handler->element($element);
 163 }
 164
 165 =head2 marker
 166
 167  Title   : marker
 168  Function: This is a synonym of the element() method
 169  Status  : deprecated, will be removed in the next version
 170
 171 =cut
 172
 173 *marker = \&element;
 174
 175 =head2 PositionI-specific methods
 176
 177 =cut
 178
 179 =head2 value
 180
 181  Title   : value
 182  Usage   : my $pos = $position->value();
 183  Function: Get/Set the value for this position
 184  Returns : scalar, value
 185  Args    : [optional] new value to set
 186
 187 =cut
 188
 189 sub value {
 190     my $self = shift;
 191     $self->throw_not_implemented();
 192 }
 193
 194 =head2 numeric
 195
 196  Title   : numeric
 197  Usage   : my $num = $position->numeric;
 198  Function: Read-only method that is guaranteed to return a numeric
 199            representation of the start of this position.
 200  Returns : scalar numeric
 201  Args    : none to get the co-ordinate normally (see absolute() method), OR
 202            Bio::Map::RelativeI to get the co-ordinate converted to be
 203            relative to what this Relative describes.
 204
 205 =cut
 206
 207 sub numeric {
 208     my $self = shift;
 209     $self->throw_not_implemented();
 210 }
 211
 212 =head2 sortable
 213
 214  Title   : sortable
 215  Usage   : my $num = $position->sortable();
 216  Function: Read-only method that is guaranteed to return a value suitable
 217            for correctly sorting this kind of position amongst other positions
 218            of the same kind on the same map. Note that sorting different kinds
 219            of position together is unlikely to give sane results.
 220  Returns : numeric
 221  Args    : none
 222
 223 =cut
 224
 225 sub sortable {
 226     my $self = shift;
 227     $self->throw_not_implemented();
 228 }
 229
 230 =head2 relative
 231
 232   Title   : relative
 233   Usage   : my $relative = $position->relative();
 234             $position->relative($relative);
 235   Function: Get/set the thing this Position's coordinates (numerical(), start(),
 236             end()) are relative to, as described by a Relative object.
 237   Returns : Bio::Map::RelativeI (default is one describing "relative to the
 238             start of the Position's map")
 239   Args    : none to get, OR
 240             Bio::Map::RelativeI to set
 241
 242 =cut
 243
 244 sub relative {
 245     my $self = shift;
 246     $self->throw_not_implemented();
 247 }
 248
 249 =head2 absolute
 250
 251   Title   : absolute
 252   Usage   : my $absolute = $position->absolute();
 253             $position->absolute($absolute);
 254   Function: Get/set how this Position's co-ordinates (numerical(), start(),
 255             end()) are reported. When absolute is off, co-ordinates are
 256             relative to the thing described by relative(). Ie. the value
 257             returned by start() will be the same as the value you set start()
 258             to. When absolute is on, co-ordinates are converted to be relative
 259             to the start of the map.
 260
 261             So if relative() currently points to a Relative object describing
 262             "relative to another position which is 100 bp from the start of
 263             the map", this Position's start() had been set to 50 and absolute()
 264             returns 1, $position->start() will return 150. If absolute() returns
 265             0 in the same situation, $position->start() would return 50.
 266
 267   Returns : boolean (default 0)
 268   Args    : none to get, OR
 269             boolean to set
 270
 271 =cut
 272
 273 sub absolute {
 274     my $self = shift;
 275     $self->throw_not_implemented();
 276 }
 277
 278 =head2 RangeI-based methods
 279
 280 =cut
 281
 282 =head2 start
 283
 284   Title   : start
 285   Usage   : my $start = $position->start();
 286             $position->start($start);
 287   Function: Get/set the start co-ordinate of this position.
 288   Returns : the start of this position
 289   Args    : scalar numeric to set, OR
 290             none to get the co-ordinate normally (see absolute() method), OR
 291             Bio::Map::RelativeI to get the co-ordinate converted to be
 292             relative to what this Relative describes.
 293
 294 =cut
 295
 296 =head2 end
 297
 298   Title   : end
 299   Usage   : my $end = $position->end();
 300             $position->end($end);
 301   Function: Get/set the end co-ordinate of this position.
 302   Returns : the end of this position
 303   Args    : scalar numeric to set, OR
 304             none to get the co-ordinate normally (see absolute() method), OR
 305             Bio::Map::RelativeI to get the co-ordinate converted to be
 306             relative to what this Relative describes.
 307
 308 =cut
 309
 310 =head2 length
 311
 312   Title   : length
 313   Usage   : $length = $position->length();
 314   Function: Get the length of this position.
 315   Returns : the length of this position
 316   Args    : none
 317
 318 =cut
 319
 320 =head2 strand
 321
 322   Title   : strand
 323   Usage   : $strand = $position->strand();
 324   Function: Get the strand of this position; it is always 1 since maps to not
 325             have strands.
 326   Returns : 1
 327   Args    : none
 328
 329 =cut
 330
 331 sub strand {
 332     return 1;
 333 }
 334
 335 =head2 toString
 336
 337   Title   : toString
 338   Usage   : print $position->toString(), "\n";
 339   Function: stringifies this range
 340   Returns : a string representation of the range of this Position
 341   Args    : optional Bio::Map::RelativeI to have the co-ordinates reported
 342             relative to the thing described by that Relative
 343
 344 =cut
 345
 346 sub toString {
 347     my $self = shift;
 348     $self->throw_not_implemented();
 349 }
 350
 351 =head1 RangeI-related methods
 352
 353 These methods work by considering only the values of start() and end(), as
 354 modified by considering every such co-ordinate relative to the start of the map
 355 (ie. absolute(1) is set temporarily during the calculation), or any supplied
 356 Relative. For the boolean methods, when the comparison Position is on the same
 357 map as the calling Position, there is no point supplying a Relative since the
 358 answer will be the same as without. Relative is most useful when comparing
 359 Positions on different maps and you have a Relative that describes some special
 360 place on each map like 'the start of the gene', where the actual start of the
 361 gene relative to the start of the map is different for each map.
 362
 363 The methods do not consider maps during their calculations - things on different
 364 maps can overlap/contain/intersect/etc. each other.
 365
 366 The geometrical methods (intersect, union etc.) do things to the geometry of
 367 ranges, and return Bio::Map::PositionI compliant objects or triplets (start,
 368 stop, strand) from which new positions could be built. When a PositionI is made
 369 it will have a map transferred to it if all the arguments shared the same map.
 370 If a Relative was supplied the result will have that same Relative.
 371
 372 Note that the strand-testing args are there for compatibility with the RangeI
 373 interface. They have no meaning when only using PositionI objects since maps do
 374 not have strands. Typically you will just set the argument to undef if you want
 375 to supply the argument after it.
 376
 377 =head2 equals
 378
 379   Title   : equals
 380   Usage   : if ($p1->equals($p2)) {...}
 381   Function: Test whether $p1 has the same start, end, length as $p2.
 382   Returns : true if they are describing the same position (regardless of map)
 383   Args    : arg #1 = a Bio::RangeI (eg. a Bio::Map::Position) to compare this
 384                      one to (mandatory)
 385             arg #2 = optional strand-testing arg ('strong', 'weak', 'ignore')
 386             arg #3 = optional Bio::Map::RelativeI to ask if the Positions
 387                      equal in terms of their relative position to the thing
 388                      described by that Relative
 389
 390 =cut
 391
 392 sub equals {
 393     # overriding the RangeI implementation so we can handle Relative
 394     my ($self, $other, $so, $rel) = @_;
 395
 396     my ($own_start, $own_end) = $self->_pre_rangei($self, $rel);
 397     my ($other_start, $other_end) = $self->_pre_rangei($other, $rel);
 398
 399     return ($self->_testStrand($other, $so) and
 400             $own_start == $other_start and $own_end == $other_end);
 401 }
 402
 403
 404 =head2 less_than
 405
 406  Title   : less_than
 407  Usage   : if ($position->less_than($other_position)) {...}
 408  Function: Ask if this Position ends before another starts.
 409  Returns : boolean
 410  Args    : arg #1 = a Bio::RangeI (eg. a Bio::Map::Position) to compare this
 411                     one to (mandatory)
 412            arg #2 = optional Bio::Map::RelativeI to ask if the Position is less
 413                     in terms of their relative position to the thing described
 414                     by that Relative
 415
 416 =cut
 417
 418 sub less_than {
 419     my ($self, $other, $rel) = @_;
 420
 421     my ($own_start, $own_end) = $self->_pre_rangei($self, $rel);
 422     my ($other_start, $other_end) = $self->_pre_rangei($other, $rel);
 423
 424     return $own_end < $other_start;
 425 }
 426
 427 =head2 greater_than
 428
 429  Title   : greater_than
 430  Usage   : if ($position->greater_than($other_position)) {...}
 431  Function: Ask if this Position starts after another ends.
 432  Returns : boolean
 433  Args    : arg #1 = a Bio::RangeI (eg. a Bio::Map::Position) to compare this
 434                     one to (mandatory)
 435            arg #2 = optional Bio::Map::RelativeI to ask if the Position is
 436                     greater in terms of their relative position to the thing
 437                     described by that Relative
 438
 439 =cut
 440
 441 sub greater_than {
 442     my ($self, $other, $rel) = @_;
 443
 444     my ($own_start, $own_end) = $self->_pre_rangei($self, $rel);
 445     my ($other_start, $other_end) = $self->_pre_rangei($other, $rel);
 446
 447     return $own_start > $other_end;
 448 }
 449
 450 =head2 overlaps
 451
 452   Title   : overlaps
 453   Usage   : if ($p1->overlaps($p2)) {...}
 454   Function: Tests if $p1 overlaps $p2.
 455   Returns : True if the positions overlap (regardless of map), false otherwise
 456   Args    : arg #1 = a Bio::RangeI (eg. a Bio::Map::Position) to compare this
 457                      one to (mandatory)
 458             arg #2 = optional strand-testing arg ('strong', 'weak', 'ignore')
 459             arg #3 = optional Bio::Map::RelativeI to ask if the Positions
 460                      overlap in terms of their relative position to the thing
 461                      described by that Relative
 462             arg #4 = optional minimum percentage length of the overlap before
 463                      reporting an overlap exists (default 0)
 464
 465 =cut
 466
 467 sub overlaps {
 468     # overriding the RangeI implementation so we can handle Relative
 469     my ($self, $other, $so, $rel, $min_percent) = @_;
 470     $min_percent ||= 0;
 471
 472     my ($own_min, $other_min) = (0, 0);
 473     if ($min_percent > 0) {
 474         $own_min = (($self->length / 100) * $min_percent) - 1;
 475         $other_min = (($other->length / 100) * $min_percent) - 1;
 476     }
 477
 478     my ($own_start, $own_end) = $self->_pre_rangei($self, $rel);
 479     my ($other_start, $other_end) = $self->_pre_rangei($other, $rel);
 480
 481     return ($self->_testStrand($other, $so) and not
 482             (($own_start + $own_min > $other_end or $own_end - $own_min < $other_start) ||
 483              ($own_start > $other_end - $other_min or $own_end < $other_start + $other_min)));
 484 }
 485
 486 =head2 contains
 487
 488   Title   : contains
 489   Usage   : if ($p1->contains($p2)) {...}
 490   Function: Tests whether $p1 totally contains $p2.
 491   Returns : true if the argument is totally contained within this position
 492             (regardless of map), false otherwise
 493   Args    : arg #1 = a Bio::RangeI (eg. a Bio::Map::Position) to compare this
 494                      one to, or scalar number (mandatory)
 495             arg #2 = optional strand-testing arg ('strong', 'weak', 'ignore')
 496             arg #3 = optional Bio::Map::RelativeI to ask if the Position
 497                      is contained in terms of their relative position to the
 498                      thing described by that Relative
 499
 500 =cut
 501
 502 sub contains {
 503     # overriding the RangeI implementation so we can handle Relative
 504     my ($self, $other, $so, $rel) = @_;
 505
 506     my ($own_start, $own_end) = $self->_pre_rangei($self, $rel);
 507     my ($other_start, $other_end) = $self->_pre_rangei($other, $rel);
 508
 509     return ($self->_testStrand($other, $so) and
 510                         $other_start >= $own_start and $other_end <= $own_end);
 511 }
 512
 513 =head2 intersection
 514
 515  Title   : intersection
 516  Usage   : ($start, $stop, $strand) = $p1->intersection($p2)
 517            ($start, $stop, $strand) = Bio::Map::Position->intersection(\@positions);
 518            $mappable = $p1->intersection($p2, undef, $relative);
 519            $mappable = Bio::Map::Position->intersection(\@positions);
 520  Function: gives the range that is contained by all ranges
 521  Returns : undef if they do not overlap, OR
 522            Bio::Map::Mappable object who's positions are the
 523            cross-map-calculated intersection of the input positions on all the
 524            maps that the input positions belong to, OR, in list context, a three
 525            element array (start, end, strand)
 526  Args    : arg #1 = [REQUIRED] a Bio::RangeI (eg. a Bio::Map::Position) to
 527                     compare this one to, or an array ref of Bio::RangeI
 528            arg #2 = optional strand-testing arg ('strong', 'weak', 'ignore')
 529            arg #3 = optional Bio::Map::RelativeI to ask how the Positions
 530                     intersect in terms of their relative position to the thing
 531                     described by that Relative
 532
 533 =cut
 534
 535 sub intersection {
 536     # overriding the RangeI implementation so we can transfer map and handle
 537     # Relative
 538     my ($self, $given, $so, $rel) = @_;
 539         $self->throw("missing arg: you need to pass in another argument") unless $given;
 540
 541     my @positions;
 542     if ($self eq "Bio::Map::PositionI") {
 543                 $self = "Bio::Map::Position";
 544                 $self->warn("calling static methods of an interface is deprecated; use $self instead");
 545         }
 546         if (ref $self) {
 547                 push(@positions, $self);
 548         }
 549     ref($given) eq 'ARRAY' ? push(@positions, @{$given}) : push(@positions, $given);
 550     $self->throw("Need at least 2 Positions") unless @positions >= 2;
 551
 552     my ($intersect, $i_start, $i_end, $c_start, $c_end, %known_maps);
 553     while (@positions > 0) {
 554         unless ($intersect) {
 555             $intersect = shift(@positions);
 556             ($i_start, $i_end) = $self->_pre_rangei($intersect, $rel);
 557             my $map = $intersect->map;
 558             $known_maps{$map->unique_id} = $map;
 559         }
 560
 561         my $compare = shift(@positions);
 562         ($c_start, $c_end) = $self->_pre_rangei($compare, $rel);
 563         return unless $compare->_testStrand($intersect, $so);
 564         if ($compare->isa('Bio::Map::PositionI')) {
 565             my $this_map = $compare->map;
 566             if ($this_map) {
 567                 $known_maps{$this_map->unique_id} = $this_map;
 568             }
 569         }
 570         else {
 571             $self->throw("Only Bio::Map::PositionI objects are supported, not [$compare]");
 572         }
 573
 574         my @starts = sort {$a <=> $b} ($i_start, $c_start);
 575         my @ends   = sort {$a <=> $b} ($i_end, $c_end);
 576
 577         my $start = pop @starts; # larger of the 2 starts
 578         my $end = shift @ends;   # smaller of the 2 ends
 579
 580         my $intersect_strand;    # strand for the intersection
 581         if (defined($intersect->strand) && defined($compare->strand) && $intersect->strand == $compare->strand) {
 582             $intersect_strand = $compare->strand;
 583         }
 584         else {
 585             $intersect_strand = 0;
 586         }
 587
 588         if ($start > $end) {
 589             return;
 590         }
 591         else {
 592             $intersect = $self->new(-start  => $start,
 593                                     -end    => $end,
 594                                     -strand => $intersect_strand);
 595         }
 596     }
 597
 598     $intersect || return;
 599     my ($start, $end, $strand) = ($intersect->start, $intersect->end, $intersect->strand);
 600
 601     my @intersects;
 602     foreach my $known_map (values %known_maps) {
 603         my $new_intersect = $intersect->new(-start => $start,
 604                                             -end => $end,
 605                                             -strand => $strand,
 606                                             -map => $known_map);
 607         $new_intersect->relative($rel) if $rel;
 608         push(@intersects, $new_intersect);
 609     }
 610     unless (@intersects) {
 611         $intersect->relative($rel) if $rel;
 612         @intersects = ($intersect);
 613     }
 614
 615     my $result = Bio::Map::Mappable->new();
 616     $result->add_position(@intersects); # sneaky, add_position can take a list of positions
 617     return $result;
 618 }
 619
 620 =head2 union
 621
 622  Title   : union
 623  Usage   : ($start, $stop, $strand) = $p1->union($p2);
 624            ($start, $stop, $strand) = Bio::Map::Position->union(@positions);
 625            my $mappable = $p1->union($p2);
 626            my $mappable = Bio::Map::Position->union(@positions);
 627  Function: finds the minimal position/range that contains all of the positions
 628  Returns : Bio::Map::Mappable object who's positions are the
 629            cross-map-calculated union of the input positions on all the maps
 630            that the input positions belong to, OR, in list context, a three
 631            element array (start, end, strand)
 632  Args    : a Bio::Map::PositionI to compare this one to, or a list of such
 633            OR
 634            a single Bio::Map::PositionI or array ref of such AND a
 635            Bio::Map::RelativeI to ask for the Position's union in terms of their
 636            relative position to the thing described by that Relative
 637
 638 =cut
 639
 640 sub union {
 641     # overriding the RangeI implementation so we can transfer map and handle
 642     # Relative
 643     my ($self, @args) = @_;
 644     $self->throw("Not enough arguments") unless @args >= 1;
 645
 646     my @positions;
 647     my $rel;
 648     if ($self eq "Bio::Map::PositionI") {
 649                 $self = "Bio::Map::Position";
 650                 $self->warn("calling static methods of an interface is deprecated; use $self instead");
 651         }
 652         if (ref $self) {
 653                 push(@positions, $self);
 654         }
 655     if (ref $args[0] eq 'ARRAY') {
 656         push(@positions, @{shift(@args)});
 657     }
 658     else {
 659         push(@positions, shift(@args));
 660     }
 661     if ($args[0] && $args[0]->isa('Bio::Map::RelativeI')) {
 662         $rel = shift(@args);
 663     }
 664     foreach my $arg (@args) {
 665         # avoid pushing undefined values into @positions
 666         push(@positions, $arg) if $arg;
 667     }
 668     $self->throw("Need at least 2 Positions") unless @positions >= 2;
 669
 670     my (@starts, @ends, %known_maps, $union_strand);
 671     foreach my $compare (@positions) {
 672         # RangeI union allows start or end to be undefined; however _pre_rangei
 673         # will throw
 674         my ($start, $end) = $self->_pre_rangei($compare, $rel);
 675
 676         if ($compare->isa('Bio::Map::PositionI')) {
 677             my $this_map = $compare->map;
 678             if ($this_map) {
 679                 $known_maps{$this_map->unique_id} = $this_map;
 680             }
 681         }
 682         else {
 683             $self->throw("Only Bio::Map::PositionI objects are supported, not [$compare]");
 684         }
 685
 686         if (! defined $union_strand) {
 687                         $union_strand = $compare->strand;
 688                 }
 689         else {
 690                         if (! defined $compare->strand or $union_strand ne $compare->strand) {
 691                                 $union_strand = 0;
 692                         }
 693                 }
 694
 695         push(@starts, $start);
 696         push(@ends, $end);
 697     }
 698
 699         @starts = sort { $a <=> $b } @starts;
 700         @ends   = sort { $a <=> $b } @ends;
 701         my $start = shift @starts;
 702         my $end = pop @ends;
 703
 704     my @unions;
 705     foreach my $known_map (values %known_maps) {
 706         my $new_union = $self->new(-start => $start,
 707                                    -end => $end,
 708                                    -strand => $union_strand,
 709                                    -map => $known_map);
 710         $new_union->relative($rel) if $rel;
 711         push(@unions, $new_union);
 712     }
 713     unless (@unions) {
 714         @unions = ($self->new(-start => $start,
 715                          -end => $end,
 716                          -strand => $union_strand));
 717         $unions[0]->relative($rel) if $rel;
 718     }
 719
 720     my $result = Bio::Map::Mappable->new();
 721     $result->add_position(@unions); # sneaky, add_position can take a list of positions
 722     return $result;
 723 }
 724
 725 =head2 overlap_extent
 726
 727  Title   : overlap_extent
 728  Usage   : ($a_unique,$common,$b_unique) = $a->overlap_extent($b)
 729  Function: Provides actual amount of overlap between two different
 730            positions
 731  Example :
 732  Returns : array of values containing the length unique to the calling
 733            position, the length common to both, and the length unique to
 734            the argument position
 735  Args    : a position
 736
 737 =cut
 738
 739 #*** should this be overridden from RangeI?
 740
 741 =head2 disconnected_ranges
 742
 743  Title   : disconnected_ranges
 744  Usage   : my @disc_ranges = Bio::Map::Position->disconnected_ranges(@ranges);
 745  Function: Creates the minimal set of positions such that each input position is
 746            fully contained by at least one output position, and none of the
 747            output positions overlap.
 748  Returns : Bio::Map::Mappable with the calculated disconnected ranges
 749  Args    : a Bio::Map::PositionI to compare this one to, or a list of such,
 750            OR
 751            a single Bio::Map::PositionI or array ref of such AND a
 752            Bio::Map::RelativeI to consider all Position's co-ordinates in terms
 753            of their relative position to the thing described by that Relative,
 754            AND, optionally, an int for the minimum percentage of overlap that
 755            must be present before considering two ranges to be overlapping
 756            (default 0)
 757
 758 =cut
 759
 760 sub disconnected_ranges {
 761     # overriding the RangeI implementation so we can transfer map and handle
 762     # Relative
 763     my ($self, @args) = @_;
 764     $self->throw("Not enough arguments") unless @args >= 1;
 765
 766     my @positions;
 767     my $rel;
 768     my $overlap = 0;
 769     if ($self eq "Bio::Map::PositionI") {
 770                 $self = "Bio::Map::Position";
 771                 $self->warn("calling static methods of an interface is deprecated; use $self instead");
 772         }
 773         if (ref $self) {
 774                 push(@positions, $self);
 775         }
 776     if (ref $args[0] eq 'ARRAY') {
 777         push(@positions, @{shift(@args)});
 778     }
 779     else {
 780         push(@positions, shift(@args));
 781     }
 782     if ($args[0] && $args[0]->isa('Bio::Map::RelativeI')) {
 783         $rel = shift(@args);
 784         $overlap = shift(@args);
 785     }
 786     foreach my $arg (@args) {
 787         push(@positions, $arg) if $arg;
 788     }
 789     $self->throw("Need at least 2 Positions") unless @positions >= 2;
 790
 791     my %known_maps;
 792     foreach my $pos (@positions) {
 793         $pos->isa('Bio::Map::PositionI') || $self->throw("Must supply only Bio::Map::PositionI objects, not [$pos]");
 794         my $map = $pos->map || next;
 795         $known_maps{$map->unique_id} = $map;
 796     }
 797     my %prior_positions;
 798     foreach my $map (values %known_maps) {
 799         foreach my $pos ($map->get_positions) {
 800             $prior_positions{$pos} = 1;
 801         }
 802     }
 803
 804     my @outranges = ();
 805     foreach my $inrange (@positions) {
 806         my @outranges_new = ();
 807         my %overlapping_ranges = ();
 808
 809         for (my $i=0; $i<@outranges; $i++) {
 810             my $outrange = $outranges[$i];
 811             if ($inrange->overlaps($outrange, undef, $rel, $overlap)) {
 812                 my $union_able = $inrange->union($outrange, $rel); # using $inrange->union($outrange, $rel); gives >6x speedup,
 813                                                                    # but different answer, not necessarily incorrect...
 814                 foreach my $pos ($union_able->get_positions) {
 815                     $overlapping_ranges{$pos->toString} = $pos; # we flatten down to a result on a single map
 816                                                                 # to avoid creating 10s of thousands of positions during this process;
 817                                                                 # we then apply the final answer to all maps at the very end
 818                     last;
 819                 }
 820             }
 821             else {
 822                 push(@outranges_new, $outrange);
 823             }
 824         }
 825
 826         @outranges = @outranges_new;
 827
 828         my @overlappers = values %overlapping_ranges;
 829         if (@overlappers) {
 830             if (@overlappers > 1) {
 831                 my $merged_range_able = shift(@overlappers)->union(\@overlappers, $rel);
 832                 push(@outranges, $merged_range_able->get_positions);
 833             }
 834             else {
 835                 push(@outranges, @overlappers);
 836             }
 837         }
 838         else {
 839             push(@outranges, $self->new(-start => $inrange->start($rel), -end => $inrange->end($rel), -strand => $inrange->strand, -map => $inrange->map, -relative => $rel));
 840         }
 841     }
 842
 843     # purge positions that were created whilst calculating the answer, but
 844     # aren't the final answer and weren't there previously
 845     my %answers = map { $_ => 1 } @outranges;
 846     foreach my $map (values %known_maps) {
 847         foreach my $pos ($map->get_positions) {
 848             if (! exists $prior_positions{$pos} && ! exists $answers{$pos}) {
 849                 $map->purge_positions($pos);
 850             }
 851         }
 852     }
 853
 854     my %post_positions;
 855     foreach my $map (values %known_maps) {
 856         foreach my $pos ($map->get_positions) {
 857             $post_positions{$pos} = 1;
 858         }
 859     }
 860
 861     @outranges || return;
 862
 863     # make an outrange on all known maps
 864     my @final_positions;
 865     foreach my $map (values %known_maps) {
 866         foreach my $pos (@outranges) {
 867             if ($pos->map eq $map) {
 868                 push(@final_positions, $pos);
 869             }
 870             else {
 871                 push(@final_positions, $pos->new(-start => $pos->start,
 872                                                  -end => $pos->end,
 873                                                  -relative => $pos->relative,
 874                                                  -map => $map));
 875             }
 876         }
 877     }
 878
 879     # assign the positions to a result mappable
 880     my $result = Bio::Map::Mappable->new();
 881     $result->add_position(@final_positions); # sneaky, add_position can take a list of positions
 882     return $result;
 883 }
 884
 885 # get start & end suitable for rangeI methods, taking relative into account
 886 sub _pre_rangei {
 887     my ($self, $other, $rel) = @_;
 888     $self->throw("Must supply an object") unless $other;
 889     if ($rel) {
 890         $self->throw("Must supply an object for the Relative argument") unless ref($rel);
 891         $self->throw("This is [$rel], not a Bio::Map::RelativeI") unless $rel->isa('Bio::Map::RelativeI');
 892     }
 893
 894     my ($other_start, $other_end);
 895     if (ref($other)) {
 896         if (ref($other) eq 'ARRAY') {
 897             $self->throw("_pre_rangei got an array");
 898         }
 899         $self->throw("This is [$other], not a Bio::RangeI object") unless defined $other && $other->isa('Bio::RangeI');
 900
 901         if ($other->isa('Bio::Map::PositionI')) {
 902             # to get the desired start/end we need the position to be on a map;
 903             # if it isn't on one temporarily place it on self's map
 904             # - this lets us have 'generic' positions that aren't on any map
 905             # but have a relative defined and can thus be usefully compared to
 906             # positions that /are/ on maps
 907             my $other_map = $other->map;
 908             unless ($other_map) {
 909                 my $self_map = $self->map || $self->throw("Trying to compare two positions but neither had been placed on a map");
 910                 $other->map($self_map);
 911             }
 912
 913             # want start and end positions relative to the supplied rel or map start
 914             $rel ||= $other->absolute_relative;
 915             $other_start = $other->start($rel);
 916             $other_end = $other->end($rel);
 917
 918             unless ($other_map) {
 919                 $self->map->purge_positions($other);
 920             }
 921         }
 922         else {
 923             $other_start = $other->start;
 924             $other_end = $other->end;
 925         }
 926     }
 927     else {
 928         $self->throw("not a number") unless looks_like_number($other);
 929         $other_start = $other_end = $other;
 930     }
 931
 932         $other->throw("start is undefined") unless defined $other_start;
 933         $other->throw("end is undefined") unless defined $other_end;
 934
 935     return ($other_start, $other_end);
 936 }
 937
 938 1;