t/GeneCoordinateMapper.t

   1 # -*-Perl-*- Test Harness script for Bioperl
   2 # $Id$
   3
   4 use strict;
   5
   6 BEGIN {
   7     use lib 't/lib';
   8     use BioperlTest;
   9
  10     test_begin(-tests => 116);
  11
  12     use_ok('Bio::Location::Simple');
  13     use_ok('Bio::Coordinate::Pair');
  14     use_ok('Bio::Coordinate::ExtrapolatingPair');
  15         use_ok('Bio::Coordinate::GeneMapper');
  16 }
  17
  18 #
  19 # Extrapolating pairs
  20 #
  21 #    No gaps returned, matches extrapolated
  22 #     returns always a match or undef
  23 #     -strict
  24 #
  25
  26
  27 # the  reverse strand pair
  28 my $inr = Bio::Location::Simple->new(-start=>2, -end=>5, -strand=>1);
  29 my $outr = Bio::Location::Simple->new(-start=>10, -end=>13, -strand=>-1);
  30 ok my $pairr = Bio::Coordinate::ExtrapolatingPair->
  31     new(-in => $inr,
  32         -out => $outr
  33        );
  34
  35 my $posr = Bio::Location::Simple->new
  36     (-start => 3, -end => 4, -strand=> 1 );
  37 my $resr = $pairr->map($posr);
  38 is $resr->start, 11;
  39 is $resr->end, 12;
  40 is $resr->strand, -1;
  41
  42
  43
  44 # propepide
  45 my $match1 = Bio::Location::Simple->new
  46     (-seq_id => 'propeptide', -start => 21, -end => 40, -strand=>1 );
  47 # peptide
  48 my $match2 = Bio::Location::Simple->new
  49     (-seq_id => 'peptide', -start => 1, -end => 20, -strand=>1 );
  50
  51 ok my $pair = Bio::Coordinate::ExtrapolatingPair->
  52     new(-in => $match1,
  53         -out => $match2,
  54         -strict => 1
  55        );
  56
  57 ok $pair->test;
  58 is $pair->strand(), 1; #  = in->strand * out->strand
  59 is $pair->in->seq_id(), 'propeptide';
  60 is $pair->strict(), 1;
  61
  62 my ($count, $pos, $pos2, $res, $match, $res2);
  63
  64 # match within
  65 $pos = Bio::Location::Simple->new
  66     (-start => 25, -end => 25, -strand=> -1 );
  67 $res = $pair->map($pos);
  68
  69 isa_ok $res, 'Bio::Location::Simple';
  70 is $res->start, 5;
  71 is $res->end, 5;
  72 is $res->strand, -1;
  73 is $res->seq_id, 'peptide';
  74
  75
  76 # match outside = undef
  77 $pos = Bio::Location::Simple->new (-start => 5, -end => 5 );
  78 $res = $pair->map($pos);
  79
  80 is $res, undef;
  81
  82 #
  83 # partial match = match
  84 #
  85 $pos2 = Bio::Location::Simple->new
  86     (-start => 20, -end => 22, -strand=> -1 );
  87
  88 ok $res = $pair->map($pos2);
  89
  90 is $res->start, 0;
  91 is $res->end, 2;
  92 is $res->seq_id, 'peptide';
  93 is $res->strand, -1;
  94
  95
  96 #
  97 # partial match2 =  match & gap
  98 #
  99 $pos2 = Bio::Location::Simple->new (-start => 40, -end => 41, -strand=> 1 );
 100 ok $res = $pair->map($pos2);
 101 is $res->start, 20;
 102 is $res->end, 20;
 103
 104 #
 105 #enveloping
 106 #
 107 $pos2 = Bio::Location::Simple->new (-start => 19, -end => 41, -strand=> 1 );
 108 ok $res = $pair->map($pos2);
 109 is $res->start, 1;
 110 is $res->end, 20;
 111
 112 #
 113 # testing the changing the strand
 114 #
 115
 116 # chr
 117 $match1 = Bio::Location::Simple->new
 118     (-seq_id => 'chr', -start => 21, -end => 40, -strand=>1 );
 119 # gene
 120 $match2 = Bio::Location::Simple->new
 121     (-seq_id => 'gene', -start => 1, -end => 20, -strand=>-1 );
 122
 123  $pair = Bio::Coordinate::ExtrapolatingPair->
 124 #my $pair = Bio::Coordinate::Pair->
 125     new(-in => $match1,
 126         -out => $match2,
 127         -strict => 0
 128        );
 129
 130 $pos = Bio::Location::Simple->new
 131     (-start => 38, -end => 40, -strand=> 1 );
 132 $res = $pair->map($pos);
 133 is $res->start, 1;
 134 is $res->end, 3;
 135 is $res->strand, -1;
 136
 137 $pos = Bio::Location::Simple->new
 138     (-start => 1, -end => 3, -strand=> 1 );
 139 $res = $pair->map($pos);
 140 is $res->start, 38;
 141 is $res->end, 40;
 142 is $res->strand, -1;
 143
 144
 145 #
 146 #
 147 # Gene Mapper
 148 #
 149 #
 150
 151 ok my $m = Bio::Coordinate::GeneMapper->new(-in => 'propeptide',
 152                                            -out => 'peptide');
 153 #$m->verbose(2);
 154
 155 is $m->peptide_offset(5), 5;
 156
 157
 158 # match within
 159 $pos = Bio::Location::Simple->new
 160     (-start => 25, -end => 25, -strand=> 1 );
 161 $res = $m->map($pos);
 162
 163 is $res->start, 20;
 164 is $res->end, 20;
 165 is $res->strand, 1;
 166 is $res->seq_id, 'peptide';
 167
 168
 169 #
 170 # nozero
 171 #
 172
 173 # match within
 174 $pos = Bio::Location::Simple->new
 175     (-start => 4, -end => 5, -strand=> 1 );
 176 $res = $m->map($pos);
 177 is $res->start, -1;
 178 is $res->end, 0;
 179
 180 is $m->nozero('in&out'), 'in&out';
 181 $res = $m->map($pos);
 182 is $res->start, -2;
 183 is $res->end, -1;
 184 is $m->nozero(0), 0;
 185
 186
 187
 188 ok $m->swap;
 189 $pos = Bio::Location::Simple->new
 190     (-start => 5, -end => 5, -strand=> 1 );
 191 $res = $m->map($pos);
 192 is $res->start, 10;
 193
 194 # cds -> propeptide
 195 is $m->in('cds'), 'cds';
 196 is $m->out('propeptide'), 'propeptide';
 197
 198 $res = $m->map($pos);
 199 is $res->start, 2;
 200 ok $res = $m->_translate($pos);
 201 is $res->start, 2;
 202 ok $res = $m->_reverse_translate($pos);
 203 is $res->start, 13;
 204 is $res->end, 15;
 205
 206 $pos = Bio::Location::Simple->new
 207     (-start => 26, -end => 26, -strand=> 1 );
 208 $m->out('peptide');
 209 $res = $m->map($pos);
 210 is $res->start, 4;
 211
 212
 213 #
 214 # frame
 215 #
 216
 217 $pos = Bio::Location::Simple->new
 218     (-start => 1, -end => 3, -strand=> 1 );
 219 $res = $m->_frame($pos);
 220 is $res->start, 1;
 221 is $res->end, 3;
 222
 223
 224 # Collection representing exons
 225 #
 226 #  cds    1   5     6   10    11  15
 227 #  exon   1   5     1   5     1   5
 228 #  gene   1   5    11   15   21   25
 229 #         |---|     |---|     |---|
 230 #-----|-----------------------|---|--
 231 # chr 1   5   9    15   19   25   29
 232 #         pair1     pair2     pair3
 233
 234 # gene
 235 my $e1 = Bio::Location::Simple->new
 236     (-seq_id => 'gene', -start => 5, -end => 9, -strand=>1 );
 237 my $e2 = Bio::Location::Simple->new
 238     (-seq_id => 'gene', -start => 15, -end => 19, -strand=>1 );
 239 my $e3 = Bio::Location::Simple->new
 240     (-seq_id => 'gene', -start => 25, -end => 29, -strand=>1 );
 241 my @cexons = ($e1, $e2, $e3);
 242
 243 $m= Bio::Coordinate::GeneMapper->new();
 244
 245 $m->in('chr');
 246 $m->out('gene');
 247 my $off = $m->cds(5);
 248 is $off->start, 5; # start of the coding region
 249 is $m->exons(@cexons), 3;
 250
 251 $m->out('exon');
 252 $pos = Bio::Location::Simple->new
 253     (-start => 6, -end => 7, -strand=> 1 );
 254 $res = $m->map($pos);
 255
 256 is $res->start, 2;
 257 is $res->end, 3;
 258
 259 $m->out('negative_intron');
 260 $pos = Bio::Location::Simple->new
 261     (-start => 12, -end => 14, -strand=> 1 );
 262 $res = $m->map($pos);
 263 is $res->start, -3;
 264 is $res->end, -1;
 265 is $res->seq_id, 'intron1';
 266
 267
 268 # cds
 269 $m->out('cds');
 270 $pos = Bio::Location::Simple->new
 271     (-start => 5, -end => 9, -strand=> 1 );
 272 $res = $m->map($pos);
 273 is $res->start, 1;
 274 is $res->end, 5;
 275
 276 $pos = Bio::Location::Simple->new
 277     (-start => 15, -end => 25, -strand=> 1 );
 278 $res = $m->map($pos);
 279 is $res->start, 6;
 280 is $res->end, 11;
 281
 282 $pos = Bio::Location::Simple->new
 283     (-start => 5, -end => 19, -strand=> 1 );
 284 $res = $m->map($pos);
 285 is $res->start, 1;
 286 is $res->end, 10;
 287
 288
 289 #
 290 # chr to cds ; ranges into one
 291 #
 292 my $exons = Bio::Location::Split->new(-seq_id => 'gene');
 293 $exons->add_sub_Location($e1);
 294 $exons->add_sub_Location($e2);
 295 $exons->add_sub_Location($e3);
 296
 297 $res = $m->map($exons);
 298 isa_ok $res,'Bio::Location::Simple';
 299 is $res->start, 1;
 300 is $res->end, 15;
 301
 302 #
 303 # cds to chr; single range into two
 304 #
 305 $m->in('cds');
 306 $m->out('gene');
 307
 308 $pos = Bio::Location::Simple->new
 309     (-start => 4, -end => 7, -strand=> 1 );
 310 $res = $m->map($pos);
 311 is $res->start, 4;
 312 is $res->end, 12;
 313
 314
 315
 316 # Collection representing exons
 317 #
 318 #  cds  -11  -7    -6  -2    -1   3  :27
 319 #  cds   -6  -2    -1 1 3     4   8  :17
 320 #  exon   1   5     1   5     1   5
 321 #  gene -21  -17  -11  -7    -1 1 3  :27
 322 #  gene -11  -7    -1 1 3     9   13 :17
 323 #         |---|     |---|     |---|
 324 #-----|-----------------------|---|--
 325 # chr 1   5   9    15   19   25   29
 326 #         pair1     pair2     pair3
 327
 328 $m= Bio::Coordinate::GeneMapper->new();
 329
 330 $m->in('chr');
 331 $m->out('gene');
 332 $off = $m->cds(17);
 333 is $off->start, 17; # start of the coding region
 334 is $m->exons(@cexons), 3;
 335
 336 # testing parameter handling in the constructor
 337 ok $m = Bio::Coordinate::GeneMapper->new(-in => 'gene',
 338                                         -out => 'peptide',
 339                                         -cds => 3,
 340                                         -exons => @cexons,
 341                                         -utr => 7,
 342                                         -peptide_offset => 5
 343                                        );
 344
 345
 346 #
 347 # Real life data
 348 # Mapping SNPs into  human serum protein MSE55 and
 349 # human galecting LGALS2 from Ensembl:
 350 #
 351
 352 #Ensembl Gene ID        Exon Start (Chr bp)     Exon End (Chr bp)       Exon Coding Start (Chr bp)
 353 #       Exon Coding End (Chr bp)        Strand
 354
 355 my @gene1_dump = split ( /\n/, qq {
 356 ENSG00000128283 34571058        34571126                        1
 357 ENSG00000128283 34576610        34577350        34576888        34577350        1
 358 ENSG00000128283 34578646        34579858        34578646        34579355        1
 359 });
 360
 361
 362 my @gene2_dump = split ( /\n/, qq {
 363 ENSG00000100079 34590438        34590464                        -1
 364 ENSG00000100079 34582387        34582469        34582387        34582469        -1
 365 ENSG00000100079 34581114        34581273        34581114        34581273        -1
 366 ENSG00000100079 34580784        34580950        34580804        34580950        -1
 367 }); # exon start should be less than end or is this intentional?
 368
 369 #Chromosome Name        Location (bp)   Strand  Reference ID
 370 my @snp_dump = split ( /\n/, qq {
 371 22      34572694        1       2235335
 372 22      34572799        1       2235336
 373 22      34572843        1       2235337
 374 22      34574896        1       2076087
 375 22      34575256        1       2076088
 376 22      34578830        1       2281098
 377 22      34579111        1       2281099
 378 22      34580411        1       2235338
 379 22      34580591        1       2281097
 380 22      34580845        1       2235339
 381 22      34581963        1       2281100
 382 22      34583722        1       140057
 383 22      34585003        1       140058
 384 22      34587726        1       968725
 385 22      34588207        1       2284055
 386 22      34591507        1       1969639
 387 22      34591949        1       140059
 388 });
 389 shift @snp_dump;
 390
 391 my ($cdsr, @exons) = read_gene_data(@gene1_dump);
 392
 393 ok my $g1 = Bio::Coordinate::GeneMapper->new(-in=>'chr', -out=>'gene');
 394 $g1->cds($cdsr);
 395
 396 #$pos = Bio::Location::Simple->new
 397 #    (-start => 34576888, -end => 34579355, -strand=> 1 );
 398 $res = $g1->map($cdsr);
 399 is $res->start, 1;
 400 is $res->end, 2468;
 401
 402 $g1->exons(@exons);
 403 $g1->in('gene');
 404 $g1->out('cds');
 405 $res = $g1->map($res);
 406 is $res->start, 1;
 407 is $res->end, 1173;
 408
 409 #map_snps($g1, @snp_dump);
 410
 411
 412 #gene 2 in reverse strand
 413 ($cdsr, @exons) = read_gene_data(@gene2_dump);
 414 ok my $g2 = Bio::Coordinate::GeneMapper->new(-in=>'chr', -out=>'gene');
 415 $g2->cds($cdsr);
 416
 417 $pos = Bio::Location::Simple->new
 418     (-start => $cdsr->end-2, -end => $cdsr->end, -strand=> 1 );
 419 $res = $g2->map($pos);
 420 is $res->start, 1;
 421 is $res->end, 3;
 422 is $res->strand, -1;
 423
 424
 425 $g2->exons(@exons);
 426
 427 #map_snps($g2, @snp_dump);
 428
 429
 430 $match1 = Bio::Location::Simple->new
 431     (-seq_id => 'a', -start => 5, -end => 17, -strand=>1 );
 432 $match2 = Bio::Location::Simple->new
 433     (-seq_id => 'b', -start => 1, -end => 13, -strand=>-1 );
 434 ok $pair = Bio::Coordinate::Pair->new(-in => $match1,
 435                                          -out => $match2,
 436                                         );
 437
 438 #
 439 # split location
 440 #
 441
 442 ok my $split = Bio::Location::Split->new();
 443 ok $split->add_sub_Location(Bio::Location::Simple->new(-start=>6,
 444                                                       -end=>8,
 445                                                       -strand=>1));
 446 $split->add_sub_Location(Bio::Location::Simple->new(-start=>15,
 447                                                    -end=>16,
 448                                                    -strand=>1));
 449
 450 $res=$pair->map($split);
 451 ok my @sublocs = $res->each_Location(1);
 452 is @sublocs, 2;
 453
 454 #print Dumper \@sublocs;
 455 is $sublocs[0]->start, 2;
 456 is $sublocs[0]->end, 3;
 457 is $sublocs[1]->start, 10;
 458 is $sublocs[1]->end, 12;
 459
 460 # testing  cds -> gene/chr which generates a split location from a simple one
 461 # exons in reverse strand!
 462 #
 463 #  pept   33222     111
 464 #  cds    8   4     3 1-1
 465 #  exon   5   1     5   1
 466 #  gene  13   9     3 1-2
 467 #         |---|     |---|
 468 #-----|-------------------
 469 # chr 1   5   9    15   19
 470 #           e1        e2
 471
 472 # gene
 473 $e1 = Bio::Location::Simple->new
 474     (-seq_id => 'gene', -start => 5, -end => 9, -strand=>-1 );
 475 $e2 = Bio::Location::Simple->new
 476     (-seq_id => 'gene', -start => 15, -end => 19, -strand=>-1 );
 477 @cexons = ($e1, $e2);
 478 my $cds= Bio::Location::Simple->new
 479     (-seq_id => 'gene', -start => 5, -end => 17, -strand=>-1 );
 480
 481 $m = Bio::Coordinate::GeneMapper->new(-in=>'cds', -out=>'chr');
 482
 483 $m->cds($cds); # this has to be set first!?
 484 is $m->exons(@cexons), 2;
 485
 486
 487 my $cds_f= Bio::Location::Simple->new
 488     (-start => 2, -end => 7, );
 489 $res = $m->map($cds_f);
 490
 491 ok @sublocs = $res->each_Location(1);
 492 is @sublocs, 2;
 493
 494 is $sublocs[0]->start, 6;
 495 is $sublocs[0]->end, 9;
 496 is $sublocs[1]->start, 15;
 497 is $sublocs[1]->end, 16;
 498
 499
 500 # test inex, exon & negative_intron
 501
 502 $m->in('gene');
 503 $m->out('inex');
 504
 505 $pos = Bio::Location::Simple->new
 506     (-seq_id => 'gene', -start => 2, -end => 10, -strand=> 1 );
 507
 508 $res = $m->map($pos);
 509 is $res->each_Location, 3;
 510
 511
 512 $m->out('intron');
 513 $res = $m->map($pos);
 514 is $res->start, 1;
 515 is $res->end, 5;
 516 is $res->strand, 1;
 517
 518 $m->out('negative_intron');
 519 $res = $m->map($pos);
 520 is $res->start, -5;
 521 is $res->end, -1;
 522 is $res->strand, 1;
 523
 524 is $m->_mapper_code2string('1-2'), 'chr-gene';
 525 is $m->_mapper_string2code('chr-gene'), '1-2';
 526
 527
 528 #todo:
 529 #  strict mapping mode
 530 #  extrapolating pair code into Bio::Coordinate::Pair ?
 531
 532
 533
 534
 535
 536
 537 sub read_gene_data {
 538     my ($self,@gene_dump) = @_;
 539     my ($cds_start, $cds_end, $strand, @exons);
 540
 541     #one line per exon
 542     my ($first, $first_line);
 543     for my $line ( @gene_dump ) {
 544
 545         my ($geneid, $exon_start, $exon_end, $exon_cstart,
 546             $exon_cend, $exon_strand) = split /\t/, $line;
 547
 548         $strand = $exon_strand if $exon_strand;
 549         #print join (' ', $geneid, $exon_start, $exon_strand), "\n";
 550
 551         # CDS location in chromosome coordinates
 552         $cds_start = $exon_cstart if !$cds_start and $exon_cstart;
 553         $cds_end = $exon_cend if $exon_cend;
 554
 555
 556         if ($exon_start > $exon_end) {
 557             ($exon_start, $exon_end) = ($exon_end, $exon_start);
 558         }
 559
 560         my $exon = Bio::Location::Simple->new
 561             (-seq_id => 'gene', -start => $exon_start,
 562              -end => $exon_end, -strand=>$strand, -verbose=>2);
 563         push @exons, $exon;
 564     }
 565
 566     if ($cds_start > $cds_end) {
 567         ($cds_start, $cds_end) = ($cds_end, $cds_start);
 568     }
 569
 570     my $cdsr = Bio::Location::Simple->new (-start => $cds_start,
 571                                            -end => $cds_end,
 572                                            -strand=> $strand);
 573
 574     return ($cdsr, @exons);
 575 }
 576
 577
 578 sub map_snps {
 579     my ($mapper, @snps) =@_;
 580     $mapper->in('chr');
 581     $mapper->out('cds');
 582     foreach my $line (@snps) {
 583         $mapper->out('cds');
 584
 585         my ($chr, $start, $strand, $id) = split /\t/, $line;
 586         my $loc = Bio::Location::Simple->new
 587             ( -start => $start,
 588              -end => $start, -strand=>$strand );
 589
 590         my $res = $mapper->map($loc);
 591         my $cds_start = 0;
 592         $cds_start = $res->start if defined $res;#defined $res->start;
 593         print $id, "\t", $cds_start, "\n";
 594
 595         # coding
 596         if ($cds_start) {
 597             $mapper->out('propeptide');
 598             my $frame_obj = $mapper->_frame($res);
 599             my $res = $mapper->map($loc);
 600             my $cds_start = 0;
 601             $cds_start = $res->start if defined $res;#defined $res->start;
 602             print  "\t\t", $cds_start, " (", $frame_obj->start, ")\n";
 603         }
 604     }
 605 }