t/Map/Map.t

   1 # -*-Perl-*- Test Harness script for Bioperl
   2 # $Id$
   3
   4 use strict;
   5
   6 BEGIN {
   7     use lib '.';
   8     use Bio::Root::Test;
   9
  10     test_begin(-tests => 267);
  11
  12     use_ok('Bio::Map::SimpleMap');
  13     use_ok('Bio::Map::Marker');
  14     use_ok('Bio::Map::Position');
  15     use_ok('Bio::Map::Relative');
  16     use_ok('Bio::Map::Mappable');
  17 }
  18
  19 ###
  20 # We explicitly test Bio::Map::SimpleMap, Bio::Map::Mappable, Bio::Map::Position,
  21 # Bio::Map::Marker and Bio::Map::Relative.
  22 #
  23 # We implicitly test Bio::Map::MapI, Bio::Map::MappableI, Bio::Map::PositionI,
  24 # and Bio::Map::PositionHandler.
  25 ###
  26
  27 # Test map basics
  28 my $map;
  29 {
  30     ok $map = Bio::Map::SimpleMap->new(-name  => 'my');
  31     ok $map->type('cyto');
  32     is $map->type, 'cyto';
  33     is $map->units, '';
  34     is $map->length, 0, "Length is ". $map->length;
  35     is $map->name, 'my';
  36     is $map->species('human'), 'human';
  37     is $map->species, 'human';
  38     is $map->unique_id, '1';
  39 }
  40
  41 # Test marker basics
  42 my $marker;
  43 {
  44     # make a plane one and add details after
  45     ok $marker = Bio::Map::Marker->new();
  46     is $marker->name('gene1'), 'gene1';
  47     ok $marker->position($map, 100);
  48     is $marker->position->value, 100;
  49     is $marker->map, $map;
  50
  51     # make positions a little easier to add by setting a default map first
  52     ok my $marker2 = Bio::Map::Marker->new(-name => 'gene3');
  53     ok $map->add_element($marker2); # one way of setting default
  54     is $marker2->default_map, $map;
  55     $marker2 = Bio::Map::Marker->new(-name => 'gene3');
  56     ok $marker2->default_map($map); # the other way of setting default
  57     is $marker2->default_map, $map;
  58     ok $marker2->position(300);
  59     is $marker2->position->value, 300;
  60     ok my $position = $marker2->position();
  61     is $position->value, 300;
  62
  63     # make one with details set in new()
  64     ok my $marker3 = Bio::Map::Marker->new(-name => 'gene2', -position => [$map, 200]);
  65     is $marker3->default_map, $map;
  66     is $marker3->position->value, 200;
  67
  68     # make one with multiple positions on multiple maps
  69     my $map2 = Bio::Map::SimpleMap->new();
  70     $marker2->positions([[$map, 150], [$map, 200], [$map2, 200], [$map2, 400]]);
  71     my @p = map($_->numeric, $marker2->each_position);
  72     is $p[0], 150;
  73     is $p[1], 200;
  74     is $p[2], 200;
  75     is $p[3], 300;
  76     is $p[4], 400;
  77     $marker2->purge_positions($map2);
  78     @p = map($_->numeric, $marker2->each_position);
  79     is $p[0], 150;
  80     is $p[1], 200;
  81     is $p[2], 300;
  82
  83     # make sure we can add positions with 0 value
  84     my $map3 = Bio::Map::SimpleMap->new();
  85     $marker->add_position($map3, 0);
  86     ok my @positions = $marker->get_positions($map3);
  87     is @positions, 1;
  88     is $positions[0]->value, 0;
  89 }
  90
  91 # Test position basics
  92 my $pos;
  93 {
  94     ok $pos = Bio::Map::Position->new();
  95     ok $pos->map($map);
  96     is $pos->map(), $map;
  97     ok $pos->element($marker);
  98     is $pos->element(), $marker;
  99
 100     ok $pos->value('10');
 101     is $pos->value(), '10';
 102     is $pos->numeric, 10;
 103     is $pos->sortable, 10;
 104     is $pos->start, 10;
 105     is $pos->end, 10;
 106
 107     # give a marker a single position with explicit position creation
 108     ok $pos = Bio::Map::Position->new(-map => $map, -value => 500);
 109     ok $marker->position($pos);
 110     ok my $got_pos = $marker->position();
 111     is $got_pos, $pos;
 112     is $marker->position->value, 500;
 113
 114     # add a position
 115     my $map2 = Bio::Map::SimpleMap->new(-name => 'genethon', -type => 'Genetic');
 116     my $pos2 = Bio::Map::Position->new(-map => $map2, -value => 100);
 117     $marker->add_position($pos2);
 118     ok my @positions = $marker->get_positions($map2);
 119     is @positions, 1;
 120     is $positions[0]->value, 100;
 121 }
 122
 123 # Test interaction of Markers and Maps via Positions
 124 {
 125     # markers know what maps they are on
 126     $marker->purge_positions;
 127     is $marker->known_maps, 0;
 128     $pos->element($marker);
 129     is $marker->known_maps, 1;
 130     ok $marker->in_map(1);
 131     ok $marker->in_map($map);
 132
 133     # maps know what markers are on themselves
 134     $map->purge_positions;
 135     my @els = $map->get_elements;
 136     is @els, 0;
 137     $pos->map($map);
 138     ok my @elements = $map->get_elements;
 139     is @elements, 1;
 140     is $elements[0], $marker;
 141
 142     # positions know what marker they are for and what map they are on
 143     is $pos->map, $map;
 144     is $pos->element, $marker;
 145 }
 146
 147 # We can compare Map objects to their own kind
 148 {
 149     # positions to positions
 150     {
 151         ok $pos->equals($pos);
 152         # these get tested properly when testing Relative, later
 153     }
 154
 155     # markers to markers
 156     {
 157         ok $marker->equals($marker);
 158         # these get tested properly when testing Mappables, later
 159     }
 160
 161     # maps to maps
 162     {
 163         my $human = Bio::Map::SimpleMap->new();
 164         my $mouse = Bio::Map::SimpleMap->new();
 165         my $chicken = Bio::Map::SimpleMap->new();
 166         my $aardvark = Bio::Map::SimpleMap->new();
 167
 168         # scenario 1: we have information about where some factors bind upstream
 169         # of a gene in 4 different species. Which factors are found in all the
 170         # species?
 171         my $fac1 = Bio::Map::Mappable->new();
 172         my $pos1 = Bio::Map::Position->new(-map => $human, -element => $fac1);
 173         my $pos2 = Bio::Map::Position->new(-map => $mouse, -element => $fac1);
 174         my $pos3 = Bio::Map::Position->new(-map => $chicken, -element => $fac1);
 175         my $pos4 = Bio::Map::Position->new(-map => $aardvark, -element => $fac1);
 176         my $fac2 = Bio::Map::Mappable->new();
 177         my $pos5 = Bio::Map::Position->new(-map => $human, -element => $fac2);
 178         my $pos6 = Bio::Map::Position->new(-map => $mouse, -element => $fac2);
 179         my $pos7 = Bio::Map::Position->new(-map => $chicken, -element => $fac2);
 180         my $fac3 = Bio::Map::Mappable->new();
 181         my $pos8 = Bio::Map::Position->new(-map => $human, -element => $fac3);
 182         my $pos9 = Bio::Map::Position->new(-map => $mouse, -element => $fac3);
 183
 184         # scenario 1 answer:
 185         ok my @factors = $human->common_elements([$mouse, $chicken, $aardvark]);
 186         is @factors, 1;
 187         ok @factors = $human->common_elements([$mouse, $chicken, $aardvark], -min_percent => 50);
 188         is @factors, 3;
 189         ok @factors = $human->common_elements([$mouse, $chicken, $aardvark], -min_percent => 50, -min_num => 3);
 190         is @factors, 2;
 191         ok @factors = $chicken->common_elements([$mouse, $human, $aardvark], -min_percent => 50, -require_self => 1);
 192         is @factors, 2;
 193         ok @factors = Bio::Map::SimpleMap->common_elements([$human, $mouse, $human, $aardvark], -min_percent => 50, -required => [$aardvark]);
 194         is @factors, 1;
 195     }
 196 }
 197
 198 # Test relative positions
 199 {
 200     my $map = Bio::Map::SimpleMap->new();
 201     my $pos1 = Bio::Map::Position->new(-map => $map, -start => 50, -length => 5);
 202     my $pos2 = Bio::Map::Position->new(-map => $map, -start => 100, -length => 5);
 203     ok my $relative = Bio::Map::Relative->new(-position => $pos2);
 204     ok $pos1->relative($relative);
 205     is $pos1->start, 50;
 206     is $pos1->absolute(1), 1;
 207     is $pos1->start, 150;
 208     is $pos1->absolute(0), 0;
 209     ok my $relative2 = Bio::Map::Relative->new(-map => 10);
 210     my $pos3 = Bio::Map::Position->new(-map => $map, -element => $marker, -start => -5, -length => 5);
 211     $pos3->relative($relative2);
 212     my $relative3 = Bio::Map::Relative->new(-position => $pos3);
 213     is $pos1->start($relative3), 145;
 214     is $pos1->numeric($relative3), 145;
 215     is $pos1->end($relative3), 149;
 216
 217     # Test the RangeI-related methods on relative positions
 218     {
 219         my $pos1 = Bio::Map::Position->new(-map => $map, -start => 50, -length => 10);
 220         my $pos2 = Bio::Map::Position->new(-map => $map, -start => 100, -length => 10);
 221         my $pos3 = Bio::Map::Position->new(-map => $map, -start => 45, -length => 1);
 222         my $pos4 = Bio::Map::Position->new(-map => $map, -start => 200, -length => 1);
 223         my $relative = Bio::Map::Relative->new(-position => $pos3);
 224         my $relative2 = Bio::Map::Relative->new(-position => $pos4);
 225         ok ! $pos1->overlaps($pos2);
 226         $pos1->relative($relative);
 227         ok $pos1->overlaps($pos2);
 228         ok $pos2->overlaps($pos1);
 229         ok $pos1->overlaps($pos2, undef, $relative2);
 230
 231         # Make sure it works with normal Ranges
 232         use Bio::Range;
 233         my $range = Bio::Range->new(-start => 100, -end => 109);
 234         ok $pos1->overlaps($range);
 235         ok ! $range->overlaps($pos1);
 236         $pos1->absolute(1);
 237         ok $range->overlaps($pos1);
 238         $pos1->absolute(0);
 239
 240         # Try the other methods briefly
 241         ok my $i = $pos1->intersection($pos2); # returns a mappable
 242         ($i) = $i->get_positions; # but we're just interested in the first (and only) position of mappable
 243         is $i->toString, '100..104';
 244         ok $i = $pos1->intersection($pos2, undef, $relative2);
 245         ($i) = $i->get_positions;
 246         is $i->toString, '-100..-96';
 247         is $i->map, $map;
 248         is $i->relative, $relative2;
 249         $i->absolute(1);
 250         is $i->toString, '100..104';
 251
 252         ok my $u = $pos1->union($pos2);
 253         ($u) = $u->get_positions;
 254         is $u->toString, '95..109';
 255         ok $u = $pos1->union($pos2, $relative2);
 256         ($u) = $u->get_positions;
 257         is $u->toString, '-105..-91';
 258         is $u->map, $map;
 259         is $u->relative, $relative2;
 260         $u->absolute(1);
 261         is $u->toString, '95..109';
 262
 263         ok ! $pos1->contains($pos2);
 264         $pos2->end(104);
 265         ok $pos1->contains($pos2);
 266         ok $pos1->contains(100);
 267
 268         ok ! $pos1->equals($pos2);
 269         $pos2->start(95);
 270         ok $pos1->equals($pos2);
 271     }
 272 }
 273
 274 # Test Mappables
 275 {
 276     ok my $gene = Bio::Map::Mappable->new();
 277     my $human = Bio::Map::SimpleMap->new();
 278     my $mouse = Bio::Map::SimpleMap->new();
 279     ok my $pos1 = Bio::Map::Position->new(-map => $human, -element => $gene, -start => 50, -length => 1000);
 280     my $pos2 = Bio::Map::Position->new(-map => $mouse, -start => 100, -length => 1000);
 281     $gene->add_position($pos2);
 282     my $gene_rel = Bio::Map::Relative->new(-element => $gene);
 283
 284     # scenario 1a: we know where a TF binds upstream of a gene in human.
 285     # we use four different programs to predict the site; how good were they?
 286     # scenaria 1b: to what extent do the predictions and known agree?
 287     my $factor = Bio::Map::Mappable->new();
 288     my $pos3 = Bio::Map::Position->new(-map => $human, -element => $factor, -start => -25, -length => 10, -relative => $gene_rel);
 289     my $perfect_prediction = Bio::Map::Mappable->new();
 290     my $pos4 = Bio::Map::Position->new(-map => $human, -element => $perfect_prediction, -start => 25, -length => 10);
 291     my $good_prediction = Bio::Map::Mappable->new();
 292     my $pos5 = Bio::Map::Position->new(-map => $human, -element => $good_prediction, -start => 24, -length => 10);
 293     my $ok_prediction = Bio::Map::Mappable->new();
 294     my $pos6 = Bio::Map::Position->new(-map => $human, -element => $ok_prediction, -start => 20, -length => 10);
 295     my $bad_prediction = Bio::Map::Mappable->new();
 296     my $pos7 = Bio::Map::Position->new(-map => $human, -element => $bad_prediction, -start => 10, -length => 10);
 297
 298     # scenario 2: we have the same program making a prediciton for a site
 299     # in two different species; is the predicted site conserved in terms of
 300     # its position relative to the gene?
 301     my $human_prediction = Bio::Map::Mappable->new();
 302     my $pos8 = Bio::Map::Position->new(-map => $human, -element => $human_prediction, -start => 25, -length => 10);
 303     my $mouse_prediction = Bio::Map::Mappable->new();
 304     my $pos9 = Bio::Map::Position->new(-map => $mouse, -element => $mouse_prediction, -start => 75, -length => 10);
 305
 306     # Test the RangeI-related methods
 307     {
 308         # scenario 1a answers:
 309         ok $perfect_prediction->equals($factor);
 310         ok $perfect_prediction->contains($factor);
 311         ok ! $ok_prediction->equals($factor);
 312         ok $ok_prediction->overlaps($factor);
 313         ok ! $bad_prediction->overlaps($factor);
 314         ok $bad_prediction->less_than($factor);
 315         ok ! $bad_prediction->greater_than($factor);
 316         ok $factor->greater_than($bad_prediction);
 317
 318         # scenario 1b answer:
 319         my $predictions = [$perfect_prediction, $good_prediction, $ok_prediction, $bad_prediction];
 320         ok my @groups = $factor->overlapping_groups($predictions, -relative => $gene_rel);
 321         is @groups, 2;
 322         is ${$groups[1]}[0], $pos7;
 323         is ${$groups[0]}[0], $pos6;
 324         is ${$groups[0]}[1], $pos5;
 325         is ${$groups[0]}[2]->toString($gene_rel), $pos4->toString($gene_rel);
 326         is ${$groups[0]}[3]->toString($gene_rel), $pos3->toString($gene_rel);
 327         ok my $di = $factor->disconnected_intersections($predictions, -relative => $gene_rel, -min_mappables_num => 3);
 328         my @di = $di->get_positions;
 329         is @di, 1;
 330         is $di[0]->toString, '-25..-21';
 331         ok my $du = $factor->disconnected_unions($predictions, -relative => $gene_rel, -min_mappables_num => 3);
 332         my @du = $du->get_positions;
 333         is @du, 1;
 334         is $du[0]->toString, '-30..-16';
 335
 336         # test the flags on overlapping_groups a bit more
 337         @groups = $factor->overlapping_groups($predictions, -relative => $gene_rel, -min_pos_num => 2);
 338         is @groups, 1;
 339         @groups = $factor->overlapping_groups($predictions, -relative => $gene_rel, -min_pos_num => 1, -min_mappables_num => 2);
 340         is @groups, 1;
 341         @groups = $factor->overlapping_groups($predictions, -relative => $gene_rel, -min_pos_num => 1, -min_mappables_num => 1, -min_mappables_percent => 50);
 342         is @groups, 1;
 343         @groups = $factor->overlapping_groups($predictions, -relative => $gene_rel, -min_pos_num => 1, -min_mappables_num => 1, -min_mappables_percent => 5);
 344         is @groups, 2;
 345         @groups = $factor->overlapping_groups($predictions, -relative => $gene_rel, -require_self => 1);
 346         is @groups, 1;
 347         @groups = $factor->overlapping_groups($predictions, -relative => $gene_rel, -required => [$factor]);
 348         is @groups, 1;
 349
 350         # scenario 2 answer:
 351         ok ! $human_prediction->overlaps($mouse_prediction);
 352         ok $human_prediction->overlaps($mouse_prediction, -relative => $gene_rel);
 353     }
 354 }
 355
 356 # complex (multi-mappable, multi-map) test of disconnected_*
 357 # and test Bio::Map::GeneMap, Bio::Map::Gene, Bio::Map::TranscriptionFactor,
 358 # Bio::Map::GeneRelative, Bio::Map::GenePosition and Bio::Map::Prediction
 359 use_ok('Bio::Map::Gene');
 360 use_ok('Bio::Map::GeneMap');
 361 use_ok('Bio::Map::TranscriptionFactor');
 362 use_ok('Bio::Map::GeneRelative');
 363 use_ok('Bio::Map::GenePosition');
 364 use_ok('Bio::Map::Prediction');
 365 {
 366     my @genes;
 367     my @predictions;
 368
 369     $genes[0] = Bio::Map::Gene->get(-universal_name => 'gene1');
 370     $genes[1] = Bio::Map::Gene->get(-universal_name => 'gene2');
 371     $genes[2] = Bio::Map::Gene->get(-universal_name => 'gene3');
 372     Bio::Map::GeneMap->get(-gene => 'gene1',
 373                            -species => 'species1',
 374                            -upstream => 1000);
 375     Bio::Map::GeneMap->get(-gene => 'gene1',
 376                            -species => 'species2',
 377                            -upstream => 2000);
 378     Bio::Map::GeneMap->get(-gene => 'gene2',
 379                            -species => 'species1',
 380                            -upstream => 1000);
 381     Bio::Map::GeneMap->get(-gene => 'gene2',
 382                            -species => 'species2',
 383                            -upstream => 2000);
 384     Bio::Map::GeneMap->get(-gene => 'gene3',
 385                            -species => 'species1',
 386                            -upstream => 1000);
 387     Bio::Map::GeneMap->get(-gene => 'gene3',
 388                            -species => 'species2',
 389                            -upstream => 2000);
 390
 391     $predictions[0] = Bio::Map::Prediction->new(-source => 'meme');
 392     Bio::Map::Position->new(-element => $predictions[0],
 393                             -map => Bio::Map::GeneMap->get(-gene => 'gene1', -species => 'species1'),
 394                             -start => 950,
 395                             -end => 960);
 396     Bio::Map::Position->new(-element => $predictions[0],
 397                             -map => Bio::Map::GeneMap->get(-gene => 'gene1', -species => 'species2'),
 398                             -start => 1950,
 399                             -end => 1960);
 400     Bio::Map::Position->new(-element => $predictions[0],
 401                             -map => Bio::Map::GeneMap->get(-gene => 'gene2', -species => 'species1'),
 402                             -start => 955,
 403                             -end => 965);
 404     Bio::Map::Position->new(-element => $predictions[0],
 405                             -map => Bio::Map::GeneMap->get(-gene => 'gene2', -species => 'species2'),
 406                             -start => 1955,
 407                             -end => 1965);
 408     $predictions[1] = Bio::Map::Prediction->new(-source => 'meme');
 409     Bio::Map::Position->new(-element => $predictions[1],
 410                             -map => Bio::Map::GeneMap->get(-gene => 'gene1', -species => 'species1'),
 411                             -start => 950,
 412                             -end => 960);
 413     Bio::Map::Position->new(-element => $predictions[1],
 414                             -map => Bio::Map::GeneMap->get(-gene => 'gene1', -species => 'species2'),
 415                             -start => 1950,
 416                             -end => 1960);
 417     Bio::Map::Position->new(-element => $predictions[1],
 418                             -map => Bio::Map::GeneMap->get(-gene => 'gene3', -species => 'species1'),
 419                             -start => 956,
 420                             -end => 966);
 421     Bio::Map::Position->new(-element => $predictions[1],
 422                             -map => Bio::Map::GeneMap->get(-gene => 'gene3', -species => 'species2'),
 423                             -start => 1956,
 424                             -end => 1966);
 425
 426     Bio::Map::Position->new(-element => $predictions[0],
 427                             -map => Bio::Map::GeneMap->get(-gene => 'gene2', -species => 'species2'),
 428                             -start => 19850,
 429                             -end => 19860);
 430     Bio::Map::Position->new(-element => $predictions[1],
 431                             -map => Bio::Map::GeneMap->get(-gene => 'gene2', -species => 'species2'),
 432                             -start => 19850,
 433                             -end => 19860);
 434
 435     my $rel = Bio::Map::GeneRelative->new(-gene => 0);
 436     my $di = Bio::Map::Mappable->disconnected_intersections(\@predictions,
 437                                                             -min_mappables_percent => 100,
 438                                                             -min_map_percent => 100,
 439                                                             -relative => $rel);
 440     my @positions = $di->get_positions;
 441     my $expected;
 442     $expected->{gene1}->{species1} = '-45..-41';
 443     $expected->{gene1}->{species2} = '-45..-41';
 444     $expected->{gene2}->{species1} = '-45..-41';
 445     $expected->{gene2}->{species2} = '-45..-41';
 446     $expected->{gene3}->{species1} = '-45..-41';
 447     $expected->{gene3}->{species2} = '-45..-41';
 448     foreach my $pos (@positions) {
 449         my $map = $pos->map;
 450         my $gname = $map ? $map->gene->universal_name : 'n/a';
 451         my $species = $map ? $map->species : 'n/a';
 452         if (defined $expected->{$gname}->{$species}) {
 453             is $expected->{$gname}->{$species}, $pos->toString;
 454             delete $expected->{$gname}->{$species};
 455         }
 456         unless (keys %{$expected->{$gname}} > 0) {
 457             delete $expected->{$gname};
 458         }
 459     }
 460     is scalar(keys %{$expected}), 0;
 461
 462     # check we don't have any extraneous positions
 463     $expected = 8;
 464     foreach my $map ($genes[0]->known_maps) {
 465         foreach my $pos ($map->get_positions) {
 466             $expected--;
 467         }
 468     }
 469     is $expected, 0;
 470     $expected = 8;
 471     foreach my $map ($genes[1]->known_maps) {
 472         foreach my $pos ($map->get_positions) {
 473             $expected--;
 474         }
 475     }
 476     is $expected, 0;
 477 }
 478
 479 {
 480         # make some maps that will represent an area around a particular gene in
 481         # particular species
 482     ok my $map1 = Bio::Map::GeneMap->get(-gene => 'BRCA2',
 483                                          -species => 'human',
 484                                          -description => 'breast cancer 2, early onset');
 485     ok my $gene = $map1->gene;
 486         my $map2 = Bio::Map::GeneMap->get(-gene => $gene,
 487                                       -species => 'mouse',
 488                                       -upstream => 500);
 489     is $map1->gene, $map2->gene;
 490     is $gene->universal_name, 'BRCA2';
 491     is $gene->description, 'breast cancer 2, early onset';
 492     is $map1->universal_name, 'BRCA2';
 493     is $map1->upstream, 1000;
 494     is $map2->upstream, 500;
 495     my $map3 = Bio::Map::GeneMap->get(-gene => 'BRCA2', -species => 'human');
 496     is $map3, $map1;
 497     is $map3->gene, $gene;
 498
 499         # model a TF that binds 500bp upstream of the BRCA2 gene in humans and
 500         # 250bp upstream of BRCA2 in mice
 501         ok my $rel = Bio::Map::GeneRelative->new(-description => "gene start");
 502     ok my $tf = Bio::Map::TranscriptionFactor->get(-universal_name => 'tf1');
 503     is $tf->universal_name, 'tf1';
 504         Bio::Map::Position->new(-map => $map1,
 505                             -element => $tf,
 506                             -start => -450,
 507                             -length => 10,
 508                             -relative => $rel);
 509         Bio::Map::Position->new(-map => $map2,
 510                             -element => $tf,
 511                             -start => -200,
 512                             -length => 10,
 513                             -relative => $rel);
 514
 515         # find out all the things that map near BRCA2 in all species
 516     my %answers = (human => ['human', 'tf1', -450, 'gene start', 551], mouse => ['mouse', 'tf1', -200, 'gene start', 301]);
 517         foreach my $map ($gene->known_maps) {
 518         my @answers = @{$answers{$map->species}};
 519                 foreach my $thing ($map->get_elements) {
 520             next if $thing eq $gene;
 521             foreach my $pos ($thing->get_positions($map)) {
 522                 is $map->species, shift @answers;
 523                 is $thing->universal_name, shift @answers;
 524                 is $pos->value, shift @answers;
 525                 is $pos->relative->description, shift @answers;
 526                 $pos->absolute(1);
 527                 is $pos->start, shift @answers;
 528             }
 529                 }
 530         is @answers, 0;
 531         delete $answers{$map->species};
 532         }
 533     is keys %answers, 0;
 534
 535     # test getting abolute positions of things relative to things that don't
 536     # even exist in the map yet: 1st exon of default transcript
 537     ok $rel = Bio::Map::GeneRelative->new(-exon => [1]);
 538     my $pos = Bio::Map::Position->new(-map => $map1,
 539                                       -element => $tf,
 540                                       -start => 5,
 541                                       -length => 5,
 542                                       -relative => $rel);
 543     is $pos->start, 5;
 544     is $pos->relative->absolute_conversion($pos), 1006;
 545     is $pos->start($pos->absolute_relative), 1006;
 546     $pos->absolute(1);
 547     is $pos->start, 1006;
 548     is $pos->end, 1010;
 549
 550     # now actually add some transcripts, exons, introns, coding etc. and retest
 551     ok my $trans = Bio::Map::GenePosition->new(-start => 0, -length => 700, -map => $map1, -type => 'transcript');
 552     ok $gene->add_transcript_position($trans);
 553     ok my $gene_pos = $gene->position($map1);
 554     is $gene_pos->start, 1001;
 555     is $gene_pos->end, 1700; # the gene position is big enough to hold the transcript
 556
 557     $trans = Bio::Map::GenePosition->new(-start => 100, -length => 800, -map => $map1, -type => 'transcript');
 558     $gene->add_transcript_position($trans);
 559     is $gene_pos->end, 1900;
 560     is $gene->active_transcript($map1), 2;
 561     my @t_pos = $gene->get_transcript_positions($map1);
 562     is @t_pos, 2;
 563
 564     # pos was relative to the active transcript, which has now changed...
 565     is $pos->start, 1106;
 566
 567     # make a new one relative to an explicit transcript
 568     $rel = Bio::Map::GeneRelative->new(-exon => [1, 2]);
 569     my $pos2 = Bio::Map::Position->new(-map => $map1,
 570                             -element => $tf,
 571                             -start => 15,
 572                             -length => 15,
 573                             -relative => $rel);
 574     is $pos2->start, 15;
 575     $pos2->absolute(1);
 576     is $pos2->start, 1116;
 577     is $pos2->end, 1130;
 578
 579     # which isn't affected by changing the active
 580     is $gene->active_transcript($map1, 1), 1;
 581     is $pos->start, 1006;
 582     is $pos2->start, 1116;
 583
 584     $map1->get_position_handler->purge_positions($pos2);
 585
 586     # add some exons to the first transcript
 587     ok my $exon = Bio::Map::GenePosition->new(-start => 0, -length => 100, -map => $map1, -type => 'exon');
 588     $gene->add_exon_position($exon, 1);
 589     is $pos->start, 1006;
 590     $exon->start(30); # not something you'd normally do; just for demo purposes
 591     is $exon->relative->absolute_conversion($exon), 1031;
 592     is $pos->start, 1036;
 593
 594     # add another exon before the previous one - this will be considered exon 1
 595     my $exon1 = Bio::Map::GenePosition->new(-start => 0, -length => 20, -map => $map1, -type => 'exon');
 596     $gene->add_exon_position($exon1, 1);
 597     is $gene->get_exon_position($map1, 2), $exon;
 598     ok my @exons = $gene->get_exon_positions($map1);
 599     is @exons, 2;
 600     is $exons[0], $exon1;
 601     is $exons[1], $exon;
 602     is $pos->start, 1006;
 603
 604     # add the intervening intron
 605     ok my $intron = Bio::Map::GenePosition->new(-start => 20, -length => 10, -map => $map1, -type => 'intron');
 606     ok ! $gene->get_intron_position($map1, 1);
 607     $gene->add_intron_position($intron, 1);
 608     is $gene->get_intron_position($map1, 1), $intron;
 609     ok ! $gene->get_intron_position($map1, 1, 2);
 610     ok ! $gene->get_intron_position($map1, 2);
 611     is $gene->get_intron_positions($map1), 1;
 612     is $intron->relative->absolute_conversion($intron), 1021;
 613
 614     # like for exon 1, we can be relative to the coding region without
 615     # defining it
 616     $rel = Bio::Map::GeneRelative->new(-coding => 0);
 617     my $pos3 = Bio::Map::Position->new(-map => $map1,
 618                             -element => $tf,
 619                             -start => -10,
 620                             -length => 5,
 621                             -relative => $rel);
 622     is $pos3->start, -10;
 623     $pos3->absolute(1);
 624     is $pos3->start, 991;
 625
 626     # add the coding region for transcript 1
 627     ok my $coding1a = Bio::Map::GenePosition->new(-start => 50, -length => 20, -map => $map1, -type => 'coding');
 628     $gene->coding_position($coding1a);
 629     is $pos3->start, 1041;
 630
 631     # try adding another coding region to the same transcript: we can't, so
 632     # the existing coding is replaced
 633     my $coding1b = Bio::Map::GenePosition->new(-start => 60, -length => 20, -map => $map1, -type => 'coding');
 634     $gene->coding_position($coding1b);
 635     is $pos3->start, 1051;
 636     ok ! $coding1a->element;
 637     # try adding things without using the add_x_position methods of Gene
 638     #...
 639
 640     # GenePositions can have sequence
 641     like $exon->seq, qr/N{70}/;
 642     my $pos4 = Bio::Map::GenePosition->new(-start => 200, -seq => 'ATGCCCAAAG', -map => $map1, -type => 'exon');
 643     is $pos4->seq, 'ATGCCCAAAG';
 644     $gene->add_exon_position($pos4, 1);
 645     is $gene->get_exon_positions($map1), 3;
 646     is $pos4->length, 10;
 647     $pos4->absolute(1);
 648     is $pos4->end, 1210;
 649     is $pos4->seq('ATGCC'), 'ATGCC';
 650     is $pos4->length, 5;
 651     is $pos4->end, 1205;
 652
 653     # so can GeneMaps
 654     my $map4 = Bio::Map::GeneMap->get(-gene => $gene,
 655                                       -species => 'chicken',
 656                                       -seq => 'ATGCCCAAAG');
 657     like $map4->seq, qr/ATGCCCAAAGN{991}/;
 658     is $map4->subseq(3,6), 'GCCC';
 659     is $map4->subseq(9,15), 'AGNNNNN'; # subseq is within map but beyond supplied seq, pads with Ns
 660
 661     # map sequence can be built from the sequence of PositionWithSequences on the map
 662     my $pos5 = Bio::Map::PositionWithSequence->new(-start => 11, -seq => 'ATG', -map => $map4);
 663     is $pos5->seq, 'ATG';
 664     is $map4->subseq(9,15), 'AGATGNN';
 665
 666     SKIP: {
 667         test_skip(-tests => 19,
 668                   -requires_modules => [qw(Bio::Tools::Run::Ensembl Bio::EnsEMBL::Registry XML::Twig)],
 669                   -requires_networking => 1);
 670
 671         # make transcript, coding, exon, intron positions on all maps for gene,
 672         # purging manually added GenePositions
 673         my $success = $gene->set_from_db;
 674
 675         skip('Failed to retreive anything from Ensembl; not sure why', 19) unless $success;
 676
 677         is $gene->get_transcript_position($map1)->toString($pos->absolute_relative), '1001..84737';
 678         is $gene->get_transcript_position($map2)->toString($pos->absolute_relative), '501..47617';
 679         is $gene->get_transcript_position($map4)->toString($pos->absolute_relative), '1373..37665';
 680         like $gene->description($map1), qr/Breast cancer/i;
 681         is $gene->display_id($map1), 'ENSG00000139618';
 682         is $gene->display_id($map2), 'ENSMUSG00000041147';
 683         is $gene->display_id($map4), 'ENSGALG00000017073';
 684         ok defined($gene->display_xref($map4));
 685         is $gene->external_name($map1), 'BRCA2';
 686         is $gene->biotype($map2), 'protein_coding';
 687         is $gene->source($map4), 'ensembl';
 688
 689         # we can add to a new map and the database info will be automatically there
 690         my $species = Bio::Species->new(-name => 'dog');
 691         $species->db_handle(Bio::DB::Taxonomy->new(-source => 'entrez'));
 692         my $map5 = Bio::Map::GeneMap->get(-gene => $gene, -species => $species);
 693         is $gene->display_id($map5), 'ENSCAFG00000006383';
 694
 695         # now the gene has a database connection, its maps and positions can get sequence
 696         ok my $seq = $map1->seq;
 697         is length($seq), 84737;
 698         is substr($seq, 0, 20), 'TGTTACAGAACCAACGAATT'; # start of upstream
 699         is substr($seq, -20, 20), 'CTACAAGTATTATTTTACAA'; # end of gene since no downstream
 700         ok substr($map1->subseq($gene->coding_position($map1)), 0, 3);
 701         #my $exon1_str = 'GGGCTTGTGGCGCGAGCTTCTGAAACTAGGCGGCAGAGGCGGAGCCGCTGTGGCACTGCTGCGCCTCTGCTGCGCCTCGGGTGTCTTTTGCGGCGGTGGGTCGCCGCCGGGAGAAGCGTGAGGGGACAGA';
 702         my $exon1_pos = $gene->get_exon_position($map1, 1);
 703         cmp_ok(length($map1->subseq($exon1_pos)), '>', 20);
 704         is $exon1_pos->seq, $map1->subseq($exon1_pos);
 705     }
 706
 707     # test a gene with multiple transcripts...
 708     #...
 709 }