1
# $Id: Gene.pm,v 1.6 2006/07/17 14:16:53 sendu Exp $
3
# BioPerl module for Bio::Map::Gene
5
# Cared for by Sendu Bala <bix@sendu.me.uk>
9
# You may distribute this module under the same terms as perl itself
11
# POD documentation - main docs before the code
15
Bio::Map::Gene - An gene modelled as a mappable element.
21
my $gene = Bio::Map::Gene->get(-universal_name => 'BRCA2',
22
-description => 'breast cancer 2, early onset');
24
# Normally you get Gene objects from GeneMaps
25
use Bio::Map::GeneMap;
27
# Model a gene with its orthologous versions found in different species,
28
# but at abstract locations within each genome
29
my $map1 = Bio::Map::GeneMap->get(-universal_name => 'BRCA2', -species => $human);
30
my $map2 = Bio::Map::GeneMap->get(-universal_name => 'BRCA2', -species => $mouse);
34
# Genes can have special kinds of positions (Bio::Map::GenePosition) that
35
# define where various sub-regions of the gene are, relative to one of the
36
# normal Positions the gene has placing it on a map.
37
my $trans = Bio::Map::GenePosition->new(-start => 0, -length => 700,
38
-map => $map1, -type => 'transcript');
39
$gene->add_transcript_position($trans);
40
my $exon = Bio::Map::GenePosition->new(-start => 0, -length => 100,
41
-map => $map1, -type => 'exon');
42
$gene->add_exon_position($exon, 1);
43
# (so now the gene has 1 transcript 700bp long which starts at the beginning
44
# of the gene, and we've defined the first of many exons which starts at the
45
# start of the transcript and is 100bp long)
49
Model a gene as an abstract mappable element. This is for when you don't care
50
exactly where a gene is in a genome, but just want to model other things (like
51
transcription factor binding sites) that are near it so you can answer questions
52
like "what binds near this gene?", or "which genes does this bind near?".
54
See t/Map/Map.t for more example usage.
60
User feedback is an integral part of the evolution of this and other
61
Bioperl modules. Send your comments and suggestions preferably to the
62
Bioperl mailing list. Your participation is much appreciated.
64
bioperl-l@bioperl.org - General discussion
65
http://bioperl.org/wiki/Mailing_lists - About the mailing lists
69
Report bugs to the Bioperl bug tracking system to help us keep track
70
of the bugs and their resolution. Bug reports can be submitted via the
73
http://bugzilla.open-bio.org/
75
=head1 AUTHOR - Sendu Bala
81
The rest of the documentation details each of the object methods.
82
Internal methods are usually preceded with a _
86
# Let the code begin...
88
package Bio::Map::Gene;
91
use Bio::Map::GenePosition;
93
use base qw(Bio::Map::Mappable);
100
# Bio::Tools::Run::Ensembl is in bioperl-run package which may not be
101
# installed, but its functionality is only optional here
102
eval {require Bio::Tools::Run::Ensembl;};
109
Usage : my $gene = Bio::Map::Gene->new();
110
Function: Builds a new Bio::Map::Gene object
111
Returns : Bio::Map::Gene
112
Args : -universal_name => string : name of the gene (in a form common to all
113
species that have the gene, but unique
114
amongst non-orthologous genes), REQUIRED
115
-description => string : free text description of the gene
120
my ($class, @args) = @_;
121
my $self = $class->SUPER::new(@args);
123
my ($u_name, $desc) = $self->_rearrange([qw(UNIVERSAL_NAME DESCRIPTION)], @args);
124
$u_name || $self->throw("You must supply a -universal_name");
125
$self->universal_name($u_name);
127
defined $desc && $self->description($desc);
135
Usage : my $gene = Bio::Map::Gene->get();
136
Function: Builds a new Bio::Map::Gene object (like new()), or gets a
137
pre-existing one that shares the same universal_name.
138
Returns : Bio::Map::Gene
139
Args : -universal_name => string, name of the gene (in a form common to all
140
species that have the gene, but unique amongst
141
non-orthologous genes), REQUIRED
142
-description => string, free text description of the gene
147
my ($class, @args) = @_;
148
my ($u_name, $desc) = Bio::Root::Root->_rearrange([qw(UNIVERSAL_NAME DESCRIPTION)], @args);
150
if ($u_name && defined $GENES->{$u_name}) {
151
$GENES->{$u_name}->description($desc) if $desc;
152
return $GENES->{$u_name};
155
return $class->new(@args);
158
=head2 universal_name
160
Title : universal_name
161
Usage : my $name = $gene->universal_name
162
Function: Get/Set Mappable name, corresponding to the name of the gene in a
163
form shared by orthologous versions of the gene in different species,
164
but otherwise unique.
166
Args : none to get, OR string to set
171
my ($self, $value) = @_;
172
if (defined $value) {
173
delete $GENES->{$self->{'_uname'}} if $self->{'_uname'};
174
$self->{'_uname'} = $value;
175
$GENES->{$value} = $self;
177
return $self->{'_uname'};
183
Usage : my $description = $gene->description();
184
$gene->description($description);
185
Function: Get/set information relating to the gene, in this case the
186
description (eg. 'full name of gene')
187
Returns : string (empty string if not defined)
188
Args : none to get general version, OR Bio::Map::GeneMap to get map-specific
190
string to set general version, optionally AND Bio::Map::GeneMap to
191
set map-specific version
197
return $self->_gene_data('description', @_);
203
Usage : my $display_id = $gene->display_id();
204
$gene->display_id($display_id);
205
Function: Get/set information relating to the gene, in this case the
206
display_id (eg. 'ENSG00000155287')
207
Returns : string (empty string if not defined)
208
Args : none to get general version, OR Bio::Map::GeneMap to get map-specific
210
string to set general version, optionally AND Bio::Map::GeneMap to
211
set map-specific version
217
return $self->_gene_data('display_id', @_);
223
Usage : my $display_xref = $gene->display_xref();
224
$gene->display_xref($display_xref);
225
Function: Get/set information relating to the gene, in this case the
226
display_xref (eg. 'HUGO:23472').
227
Returns : string (empty string if not defined)
228
Args : none to get general version, OR Bio::Map::GeneMap to get map-specific
230
string to set general version, optionally AND Bio::Map::GeneMap to
231
set map-specific version
237
return $self->_gene_data('display_xref', @_);
243
Usage : my $external_db = $gene->external_db();
244
$gene->external_db($external_db);
245
Function: Get/set information relating to the gene, in this case the
246
external_db (eg. 'HUGO').
247
Returns : string (empty string if not defined)
248
Args : none to get general version, OR Bio::Map::GeneMap to get map-specific
250
string to set general version, optionally AND Bio::Map::GeneMap to
251
set map-specific version
257
return $self->_gene_data('external_db', @_);
262
Title : external_name
263
Usage : my $external_name = $gene->external_name();
264
$gene->external_name($external_name);
265
Function: Get/set information relating to the gene, in this case the (eg.
266
'gene_name', probably the same as or similar to what you set
267
universal_name() to, but could be a species-specific alternative).
268
Returns : string (empty string if not defined)
269
Args : none to get general version, OR Bio::Map::GeneMap to get map-specific
271
string to set general version, optionally AND Bio::Map::GeneMap to
272
set map-specific version
278
return $self->_gene_data('external_name', @_);
284
Usage : my $biotype = $gene->biotype();
285
$gene->biotype($biotype);
286
Function: Get/set information relating to the gene, in this case the biotype
287
(eg. 'protein_coding').
288
Returns : string (empty string if not defined)
289
Args : none to get general version, OR Bio::Map::GeneMap to get map-specific
291
string to set general version, optionally AND Bio::Map::GeneMap to
292
set map-specific version
298
return $self->_gene_data('biotype', @_);
304
Usage : my $source = $gene->source();
305
$gene->source($source);
306
Function: Get/set information relating to the gene, in this case the source
308
Returns : string (empty string if not defined)
309
Args : none to get general version, OR Bio::Map::GeneMap to get map-specific
311
string to set general version, optionally AND Bio::Map::GeneMap to
312
set map-specific version
318
return $self->_gene_data('source', @_);
324
Usage : my $position = $mappable->position($map);
325
Function: Get the main Position of this Mappable on a given map. (A gene may
326
have many positions on a map, but all but one of them are
327
Bio::Map::GenePosition objects that describe sub-regions of the gene
328
which are relative to the 'main' Bio::Map::Position position, which
329
is the only one that is directly relative to the map - this is the
330
Position returned by this method.)
331
Returns : Bio::Map::Position
332
Args : L<Bio::Map::MapI> object.
337
my ($self, $map) = @_;
338
($map && $self->in_map($map)) || return;
340
foreach my $pos ($self->get_positions($map, 1)) {
341
next if $pos->isa('Bio::Map::GenePosition');
343
#*** could do sanity checking; there should only be 1 non-GenePosition
344
# object here, and it should have a relative of type 'map', and it
345
# should sort before or equal to all other positions
349
=head2 add_transcript_position
351
Title : add_transcript_position
352
Usage : $gene->add_transcript_position($position);
353
Function: Set the bounds of a transcript on a map (that of the supplied
354
position). All transcript positions added this way must have
355
coordinates relative to the main position of the 'gene' mappable on
356
this transcript's map. The first position added using this method
357
must have a start of 0. The supplied Position will be given a type of
358
'transcript' and relative of (gene => 0). The active_transcript for
359
the Position's map will be set to this one.
361
Args : Bio::Map::GenePosition (which must have its map() defined, and be for
362
a map this gene is on)
366
sub add_transcript_position {
367
my ($self, $pos) = @_;
368
($pos && $pos->isa('Bio::Map::GenePosition')) || return;
370
my $map = $pos->map || $self->throw("Supplied GenePosition has no map");
371
$self->in_map($map) || $self->throw("Supplied GenePosition is not on a map that this gene belong to");
372
my @transcripts = $self->get_transcript_positions($map);
373
if (@transcripts == 0) {
374
# first transcript needs start of 0
375
if ($pos->start != 0) {
376
$self->warn("The first transcript position added to a map needs a start of 0, not adding");
381
$pos->type('transcript');
382
$pos->relative->gene(0);
383
$self->SUPER::add_position($pos);
385
# need to remember the order these were added, but remember what we store
386
# here could become invalid if positions are purged outside of this class
387
push(@{$self->{t_order}->{$map}}, $pos);
389
# adjust main position's length to hold this transcript
390
my $main_pos = $self->position($map);
391
my $increase = ($pos->length + $pos->start($pos->absolute_relative)) - ($main_pos->end + 1);
393
$main_pos->end($main_pos->end + $increase);
396
# make this new transcript the active one
397
$self->active_transcript($map, scalar(@transcripts) + 1);
400
=head2 active_transcript
402
Title : active_transcript
403
Usage : my $active = $gene->active_transcript($map);
404
$gene->active_transcript($map, $int);
405
Function: Get/set the active transcript number (an int of 1 would mean the 1st
406
transcript position added to the object for the given map, ie. would
407
correspond to the the 1st Position object in the list returned by
408
get_transcript_positions($map)). The active transcript is the one
409
considered by other methods and objects when dealing with positions
410
relative to 'the' transcript.
411
Returns : int, 0 means there were no transcript positions on the given map,
412
undef is some other problem
413
Args : Just Bio::Map::GeneMap to get
414
Bio::Map::GeneMap AND int to set
418
sub active_transcript {
419
my ($self, $map, $int) = @_;
422
my @transcripts = $self->get_transcript_positions($map);
423
if (@transcripts > 0) {
425
if ($int > 0 && $int <= @transcripts) {
426
$self->{active_transcript}->{$map} = $int;
430
$self->warn("Supplied int '$int' not a good number (higher than the number of transcripts on the map?)");
435
if (defined $self->{active_transcript}->{$map}) {
436
return $self->{active_transcript}->{$map};
439
# default to the total number of transcripts on the map, ie. the
440
# most recently added
441
$self->{active_transcript}->{$map} = @transcripts;
442
return $self->{active_transcript}->{$map};
449
=head2 get_transcript_positions
451
Title : get_transcript_positions
452
Usage : my @transcript_positions = $gene->get_transcript_positions($map);
453
Function: Get all the transcript positions of this gene on the given map, in
454
the order they were added to the map.
455
Returns : list of Bio::Map::GenePosition
456
Args : Bio::Map::GeneMap
460
sub get_transcript_positions {
461
my ($self, $map) = @_;
463
$map->isa('Bio::Map::GeneMap') or return;
464
return $self->_get_typed_positions($map, 'transcript');
467
=head2 get_transcript_position
469
Title : get_transcript_position
470
Usage : my $position = $gene->get_transcript_position($map, $int);
471
Function: Get the $int'th transcript position added to the map. If no
472
transcripts have been added to the map, and the default transcript
473
was requested, $gene->position is returned, as that will have the
474
same start and end as the first transcript.
475
Returns : Bio::Map::GenePosition
476
Args : Bio::Map::GeneMap AND int (if int not supplied, or 0, returns
477
the currently active transcript position)
481
sub get_transcript_position {
482
my ($self, $map, $value) = @_;
484
$value ||= $self->active_transcript($map);
485
my @transcripts = $self->get_transcript_positions($map);
486
if (@transcripts == 0 && $value == 0) {
487
return $self->position($map);
489
return $self->_get_list_element($value, @transcripts);
492
=head2 coding_position
494
Title : coding_position
495
Usage : $gene->coding_position($position, $transcript_number);
496
$gene->coding_position($map, $transcript_number);
497
Function: Get/set the bounds of a coding region of a given transcript on a map
498
(that of the supplied position).
500
When setting, coordinates must be relative to the transcript start.
501
The supplied position will be given a type 'coding' and a relative
502
(-transcript => $transcript_number). There can be only one coding
503
position per transcript (hence this is a get/set).
505
When getting, if a coding region has not been defined for the
506
requested transcript, $gene->get_transcript_position($map,
507
$transcript_number) is returned, as if assuming the entirety of the
508
transcript is coding.
510
Returns : Bio::Map::GenePosition
511
Args : Bio::Map::GeneMap AND int (the transcript number) to get, OR to set:
512
Bio::Map::GenePosition (which must have its map() defined, and be for
513
a map this gene is on) AND int (the transcript number)
514
In both cases, if transcript number not supplied or 0 this will be
515
resolved to the current active transcript number - there must be at
516
least one transcript on the map
520
sub coding_position {
521
my ($self, $thing, $transcript_num) = @_;
522
ref($thing) || return;
523
$transcript_num ||= 0;
525
# deliberate test for PositionI so _add_type_position can do nothing if
526
# its not a GenePosition
527
if ($thing->isa('Bio::Map::PositionI')) {
528
my $map = $thing->map || return;
529
my ($existing_pos) = $self->_get_typed_positions($map, 'coding', $transcript_num);
532
$self->purge_positions($existing_pos);
534
$self->_add_type_position('coding', $thing, $transcript_num);
538
my ($pos) = $self->_get_typed_positions($thing, 'coding', $transcript_num);
539
return $pos || $self->get_transcript_position($thing, $transcript_num);
542
=head2 add_exon_position
544
Title : add_exon_position
545
Usage : $gene->add_exon_position($position, $transcript_number);
546
Function: Set the bounds of an exon of a given transcript on a map (that of the
547
supplied position). Coordinates must be relative to the transcript
548
start. The supplied position will be given a type 'exon' and a
549
relative (-transcript => $transcript_number).
551
Args : Bio::Map::GenePosition (which must have its map() defined, and be for
552
a map this gene is on) AND int (the transcript number; if not
553
supplied or 0 this will be resolved to the current active transcript
554
number - there must be at least one transcript on the map)
558
sub add_exon_position {
560
$self->_add_type_position('exon', @_);
563
=head2 get_exon_positions
565
Title : get_exon_positions
566
Usage : my @positions = $gene->get_exon_positions($map, $int);
567
Function: Get all the exon positions that are relative to the $int'th
568
transcript position added to the map. Exons are returned sorted by
569
their start positions.
570
Returns : array of Bio::Map::GenePosition
571
Args : Bio::Map::GeneMap AND int (the transcript number; if second int not
572
supplied, or 0, considers the currently active transcript)
576
sub get_exon_positions {
577
my ($self, $map, $value) = @_;
580
return $self->_get_typed_positions($map, 'exon', $value);
583
=head2 get_exon_position
585
Title : get_exon_position
586
Usage : my $position = $gene->get_exon_position($map, $exon_num, $int);
587
Function: Get the $exon_num'th exon position that is relative to the $int'th
588
transcript position added to the map. Exons are numbered in Position
589
order, not the order they were added to the map. If no exons have
590
been added to the map, and the first exon was requested,
591
$gene->get_transcript_position($map, $int) is returned, as that will
592
have the same start as the first exon, and could have the same end
593
for a single exon gene.
594
Returns : Bio::Map::GenePosition
595
Args : Bio::Map::GeneMap AND int (the exon you want) AND int (the transcript
596
number; if second int not supplied, or 0, considers the currently
601
sub get_exon_position {
602
my ($self, $map, $exon_num, $value) = @_;
603
my @exons = $self->get_exon_positions($map, $value);
604
if (@exons == 0 && $exon_num == 1) {
605
return $self->get_transcript_position($map, $value);
607
return $self->_get_list_element($exon_num, @exons);
610
=head2 add_intron_position
612
Title : add_intron_position
613
Usage : $gene->add_intron_position($position, $transcript_number);
614
Function: Set the bounds of an intron of a given transcript on a map (that of
615
the supplied position). Coordinates must be relative to the
616
transcript start. The supplied position will be given a type 'intron'
617
and a relative (-transcript => $transcript_number).
619
Args : Bio::Map::GenePosition (which must have its map() defined, and be for
620
a map this gene is on) AND int (the transcript number; if not
621
supplied or 0 this will be resolved to the current active transcript
622
number - there must be at least one transcript on the map)
626
sub add_intron_position {
628
$self->_add_type_position('intron', @_);
631
=head2 get_intron_positions
633
Title : get_intron_positions
634
Usage : my @positions = $gene->get_intron_positions($map, $int);
635
Function: Get all the intron positions that are relative to the $int'th
636
transcript position added to the map. Introns are returned sorted by
637
their start positions.
638
Returns : array of Bio::Map::GenePosition
639
Args : Bio::Map::GeneMap AND int (the transcript number; if second int not
640
supplied, or 0, considers the currently active transcript)
644
sub get_intron_positions {
645
my ($self, $map, $value) = @_;
648
return $self->_get_typed_positions($map, 'intron', $value);
651
=head2 get_intron_position
653
Title : get_intron_position
654
Usage : my $position = $gene->get_intron_position($map, $intron_num, $int);
655
Function: Get the $intron_num'th intron position that is relative to the
656
$int'th transcript position added to the map. Introns are numbered in
657
Position order, not the order they were added to the map.
658
Returns : Bio::Map::GenePosition
659
Args : Bio::Map::GeneMap AND int (the intron you want) AND int (the
660
transcript number; if second int not supplied, or 0, considers the
661
currently active transcript)
665
sub get_intron_position {
666
my ($self, $map, $intron_num, $value) = @_;
667
my @introns = $self->get_intron_positions($map, $value);
668
return $self->_get_list_element($intron_num, @introns);
674
Usage : $gene->set_from_db(); # for an instance only
675
Bio::Map::Gene->set_from_db(); # decide that all future genes added
676
# to maps will be set from db
677
Function: Creates all the various types of positions (transcripts, coding,
678
exons, introns) for this gene on all its maps. The information comes
679
from an Ensembl database via Bio::Tools::Run::Ensembl. NB: will
680
purge any existing Bio::Map::GenePosition objects that were
681
previously on the maps this gene is one.
682
Returns : undef on failure, otherwise the number of maps that successfully
683
had positions added to them
684
Args : boolean (no argument/undef is treated as 1, ie. do set from db;
685
supply 0 to turn off)
687
NB: Bio::Tools::Run::Ensembl is available in the bioperl-run package;
688
see it for details on setting up a database to use.
690
Once set, any new maps (species) this gene is added to will
691
automatically also have their positions set_from_db
696
my ($self, $bool) = @_;
697
return unless $USE_ENSEMBL;
698
return unless Bio::Tools::Run::Ensembl->registry_setup();
699
defined($bool) || ($bool = 1);
701
unless (ref($self)) {
702
$SET_FROM_DB = $bool;
706
$self->{_set_from_db} = $bool;
709
foreach my $map ($self->known_maps) {
710
$success += $self->_set_from_db($map);
716
# set from db for a particular map (species)
718
my ($self, $map) = @_;
719
my $gene_name = $self->universal_name || return 0;
720
$SET_FROM_DB || $self->{_set_from_db} || return;
722
my $species = $map->species;
724
my $slice_adaptor = Bio::Tools::Run::Ensembl->get_adaptor($species, 'Slice') || return 0;
725
my $gene = Bio::Tools::Run::Ensembl->get_gene_by_name(-species => $species,
727
-use_orthologues => 'Homo sapiens',
728
-use_swiss_lookup => 1,
729
-use_entrez_lookup => 1) || return 0;
731
# attach species(map)-specific gene info to self
732
$self->description($gene->description, $map);
733
$self->display_id($gene->display_id, $map);
734
$self->display_xref($gene->display_xref->display_id, $map);
735
$self->external_db($gene->external_db, $map);
736
$self->external_name($gene->external_name, $map);
737
$self->biotype($gene->biotype, $map);
738
$self->source($gene->source, $map);
740
# get the transcripts for this map
741
my $trans_ref = $gene->get_all_Transcripts;
742
unless ($trans_ref && @{$trans_ref} > 0) {
746
# purge all existing GenePositions from the map
747
my $handler = $map->get_position_handler();
748
foreach my $pos ($map->get_positions) {
749
if ($pos->isa('Bio::Map::GenePosition')) {
750
$handler->purge_positions($pos);
754
# assume all transcripts on the same strand, sort them
755
my $strand = ${$trans_ref}[0]->strand;
756
my @transcripts = sort { $strand == -1 ? ($b->end <=> $a->end) : ($a->start <=> $b->start) } @{$trans_ref};
758
# store slice of first transcript so we can use it to get seq data, and
759
# add chromosome info to our map if not set
760
my $primary_slice = $slice_adaptor->fetch_by_transcript_stable_id($transcripts[0]->stable_id, 0);
761
my $uid = $map->unique_id;
762
@{$self->{_ensembl}->{$uid}} = ($slice_adaptor, $primary_slice, $strand);
764
#my $cyto = $map->location || Bio::Map::CytoPosition->new();
765
#unless ($cyto->chr) {
766
# $cyto->chr($primary_slice->seq_region_name);
768
#$map->location($cyto);
770
# adjustment needed to make all transcript coords relative to the start of
771
# the first transcript which must start at 0
772
my $adjust = $strand == -1 ? $transcripts[0]->end : $transcripts[0]->start;
773
my $orig_adjust = $adjust;
774
my $adjustment = sub { return $strand == -1 ? $adjust - shift() : shift() - $adjust; };
776
# go through all the transcripts, remembering the longest
777
my $longest_trans = 0;
780
foreach my $transcript (@transcripts) {
781
# length is the total number of bases the exons cover, not genomic span
782
my $length = $transcript->length();
783
if ($length > $longest_trans) {
784
$longest_trans = $length;
788
# make positions for this transcript
789
my $slice = $slice_adaptor->fetch_by_transcript_stable_id($transcript->stable_id, 0);
790
my $start = &$adjustment($slice->start());
791
my $end = &$adjustment($slice->end());
792
($start, $end) = ($end, $start) if $start > $end;
794
my $trans_pos = Bio::Map::GenePosition->new(-map => $map, -start => $start, -end => $end, -type => 'transcript');
795
$self->add_transcript_position($trans_pos);
797
# all subsequent coordinates need to be relative to the start of this
799
$adjust = $strand == -1 ? $slice->end : $slice->start;
801
# there may not be a coding region
802
if (defined($transcript->coding_region_start)) {
803
my $atg = &$adjustment($transcript->coding_region_start());
804
my $stop = &$adjustment($transcript->coding_region_end());
805
($atg, $stop) = ($stop, $atg) if $atg > $stop;
807
my $cod_pos = Bio::Map::GenePosition->new(-map => $map, -start => $atg, -end => $stop, -type => 'coding');
808
$self->coding_position($cod_pos);
812
foreach my $exon (@{$transcript->get_all_Exons}) {
813
my $start = &$adjustment($exon->start());
814
my $end = &$adjustment($exon->end());
815
($start, $end) = ($end, $start) if $start > $end;
817
my $throw_species = ref($species) ? $species->scientific_name : $species;
818
defined($end) || $self->throw("gene $gene_name in species $throw_species (".$gene->display_id.") had exon $start with no end");
819
my $pos = Bio::Map::GenePosition->new(-map => $map, -start => $start, -end => $end, -type => 'exon');
820
$self->add_exon_position($pos);
824
foreach my $intron (@{$transcript->get_all_Introns}) {
825
my $start = &$adjustment($intron->start());
826
my $end = &$adjustment($intron->end());
827
($start, $end) = ($end, $start) if $start > $end;
829
my $pos = Bio::Map::GenePosition->new(-map => $map, -start => $start, -end => $end, -type => 'intron');
830
$self->add_intron_position($pos);
833
$adjust = $orig_adjust;
834
} continue { $count++ };
836
$self->active_transcript($map, $longest);
841
# get safely sorted positions of a certain type
842
sub _get_typed_positions {
843
my ($self, $map, $type, $transcript_number) = @_;
844
if (defined $transcript_number && $transcript_number == 0) {
845
$transcript_number = $self->active_transcript($map);
849
foreach my $pos ($self->get_positions($map, 1)) {
850
$pos->isa('Bio::Map::GenePosition') || next;
851
$pos->type eq $type || next;
853
if (defined $transcript_number) {
854
my $rel = $pos->relative || next;
855
$rel->type eq 'transcript' || next;
856
my $rel_transcript_num = $rel->transcript || $self->active_transcript($map);
857
$rel_transcript_num == $transcript_number || next;
860
push(@positions, $pos);
863
# avoid sorting using $pos->sortable since we would go infinite from the
864
# call to absolute_conversion - we don't need absolute_conversion here
865
# since we know the raw starts are all relative to the same thing, or in
866
# the case of transcripts, we want them sorted in the way they were added
867
if (defined $transcript_number) {
868
# ensure we get raw start; ask for starts relative to the things
869
# the positions are relative to. Precompute answer for efficiency
870
my @sort = map { $_->[1] }
871
sort { $a->[0] <=> $b->[0] }
872
map { [$_->start($_->relative), $_] }
877
my @known_order = @{$self->{t_order}->{$map} || []};
878
@known_order || return;
880
# transcripts might have been removed, so known_order could be invalid
881
return @known_order if @known_order == @positions; #*** dangerous assumption?
882
my %exists = map { $_ => $_ } @positions;
884
foreach my $pos (@known_order) {
885
exists $exists{$pos} || next;
886
push(@new_order, $pos);
888
@{$self->{t_order}->{$map}} = @new_order;
893
# get a certain element from an array, checking the array has that element
894
sub _get_list_element {
895
my ($self, $wanted, @list) = @_;
896
($wanted && $wanted > 0) || return;
898
my $index = $wanted - 1;
899
if ($index >= 0 && $index <= $#list) {
900
return $list[$index];
905
# add a certain type of posiiton
906
sub _add_type_position {
907
my ($self, $type, $pos, $transcript_num) = @_;
908
($pos && $pos->isa('Bio::Map::GenePosition')) || return;
910
my $map = $pos->map || $self->throw("Supplied GenePosition has no map");
911
$self->in_map($map) || $self->throw("Supplied GenePosition is not on a map that this gene belong to");
913
$transcript_num ||= $self->active_transcript($map) || $self->throw("Asked to be relative to the active transcript, but there is no transcript");
915
# sanity check - must be within the transcript
916
my $transcript_pos = $self->get_transcript_position($map, $transcript_num) || $self->throw("Asked to be relative to transcript $transcript_num, but there is no such transcript");
917
$transcript_pos->end || ($self->warn("no transcript pos end for pos for gene ".$self->universal_name." and species ".$pos->map->species."!") && exit);
918
$pos->end || ($self->warn("no pos end for pos for gene ".$self->universal_name." and species ".$pos->map->species."!") && exit);
919
unless ($transcript_pos->contains($pos)) {
920
$self->warn("$type coordinates must lie within those of the transcript, not adding $type");
925
$pos->relative->transcript($transcript_num);
926
$self->SUPER::add_position($pos);
929
# get/setter for general/map-specific data
931
my ($self, $type, $thing, $map) = @_;
932
$thing or return ($self->{$type}->{general} || '');
934
if (ref($thing) && $thing->isa('Bio::Map::GeneMap')) {
935
return $self->{$type}->{$thing} || '';
938
if ($map && $map->isa('Bio::Map::GeneMap')) {
939
$self->{$type}->{$map} = $thing;
942
$self->{$type}->{general} = $thing;
947
# for exclusive use by GeneMap so it can get sequence data
949
my ($self, $map) = @_;
951
my $uid = $map->unique_id || return;
952
if (defined $self->{_ensembl}->{$uid}) {
953
return @{$self->{_ensembl}->{$uid}};
added
removed
Bio/Map/Gene.pm
