Bio/Map/PositionI.pm

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