1
# $Id: GeneMap.pm,v 1.17 2006/07/17 14:16:53 sendu Exp $
3
# BioPerl module for Bio::Map::GeneMap
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::GeneMap - A MapI implementation to represent the area around a gene
19
use Bio::Map::GeneMap;
21
use Bio::Map::TranscriptionFactor;
22
use Bio::Map::GeneRelative;
24
# make some maps that will represent an area around a particular gene in
25
# particular species (by default, the map represents the area in the genome
26
# 1000bp upstream of the gene)
27
my $map1 = Bio::Map::GeneMap->get(-gene => 'BRCA2',
29
-description => 'breast cancer 2, early onset');
30
my $map2 = Bio::Map::GeneMap->get(-gene => 'BRCA2',
33
# model a TF that binds 500bp upstream of the BRCA2 gene in humans and
34
# 250bp upstream of BRCA2 in mice
35
my $rel = Bio::Map::GeneRelative->new(-description => "gene start");
36
my $tf = Bio::Map::TranscriptionFactor->get(-universal_name => 'tf1');
37
Bio::Map::Position->new(-map => $map1,
42
Bio::Map::Position->new(-map => $map2,
48
# find out all the things that map near BRCA2 in all species
49
foreach my $map ($gene->known_maps) {
50
foreach my $thing ($map->get_elements) {
51
next if $thing eq $gene;
52
foreach my $pos ($thing->get_positions($map)) {
53
print "In species ", $map->species, ", ",
54
$thing->universal_name, " maps at ", $pos->value,
55
" relative to ", $pos->relative->description, " of gene ",
56
$gene->universal_name, "\n";
61
# a GeneMap isa PrimarySeq and so can have sequence associated with it
63
my $subseq = $map1->subseq(2,3); # TG
67
Model the abstract notion of the area around a gene - you don't care exactly
68
where this area is in the genome, you just want to be able to say "something
69
binds upstream of gene X" and "something else binds 20bp upstream of the first
72
It's useful for modelling transcription factor bindings sites, letting you find
73
out which transcription factors bind near a gene of interest, or which genes
74
are bound by a transcription factor of interest.
76
See t/Map/Map.t for more example usage.
82
User feedback is an integral part of the evolution of this and other
83
Bioperl modules. Send your comments and suggestions preferably to
84
the Bioperl mailing list. Your participation is much appreciated.
86
bioperl-l@bioperl.org - General discussion
87
http://bioperl.org/wiki/Mailing_lists - About the mailing lists
91
Report bugs to the Bioperl bug tracking system to help us keep track
92
of the bugs and their resolution. Bug reports can be submitted via the
95
http://bugzilla.open-bio.org/
97
=head1 AUTHOR - Sendu Bala
103
The rest of the documentation details each of the object methods.
104
Internal methods are usually preceded with a _
108
# Let the code begin...
110
package Bio::Map::GeneMap;
114
use Bio::Map::Position;
116
use base qw(Bio::Map::SimpleMap Bio::PrimarySeq);
123
Usage : my $obj = Bio::Map::GeneMap->new();
124
Function: Builds a new Bio::Map::GeneMap object (that has placed on it a
125
mappable element (Bio::Map::Gene) representing a gene).
126
Returns : Bio::Map::GeneMap
127
Args : -gene => string name of the gene this map will be for
128
(in a form common to all species that have the gene,
129
but unique amongst non-orthologous genes) or a
130
Bio::Map::Gene object, REQUIRED
131
-species => Bio::Taxon or string representing species, REQUIRED
132
-uid => string, unique identifier for this map (must be
133
unique amongst all gene/species combinations)
134
-description => string, free text description of the gene
135
-upstream => int, the number of bases the map extends before the
136
start of the gene element (default 1000).
137
-downstream => int, the number of bases the map extends beyond the
138
end of the gene element (default 0).
139
-seq => string, the sequence of the map, presumably the
140
genomic sequence -upstream bases of the gene,
141
including the gene, and -downstream bases of the gene
146
my ($class, @args) = @_;
147
my $self = $class->SUPER::new(@args);
149
my ($uid, $gene, $species, $desc, $up, $down, $seq) = $self->_rearrange([qw(UID
157
unless (defined $gene && defined $species) {
158
$self->throw("You must supply both -species and -gene");
161
$self->gene(-gene => $gene, -description => $desc, -upstream => $up, -downstream => $down);
162
$self->seq($seq) if $seq;
164
unless (defined($uid)) {
165
# trigger the special behaviour in our unique_id method by supplying it
166
# the unique_id we got from our parent class
167
$self->unique_id($self->unique_id);
176
Usage : my $map = Bio::Map::GeneMap->get();
177
Function: Builds a new Bio::Map::GeneMap object (like new()), or gets a
178
pre-existing one that corresponds to your arguements.
179
Returns : Bio::Map::GeneMap
180
Args : -gene => string name of the gene this map will be for
181
(in a form common to all species that have the gene,
182
but unique amongst non-orthologous genes) or a
183
Bio::Map::Gene object, REQUIRED
184
-species => Bio::Taxon or string representing species, REQUIRED
185
-uid => string, unique identifier for this map (must be
186
unique amongst all gene/species combinations)
187
-description => string, free text description of the gene
188
-upstream => int, the number of bases the map extends before the
189
start of the gene element (default 1000).
190
-downstream => int, the number of bases the map extends beyond the
191
end of the gene element (default 0).
192
-seq => string, the sequence of the map, presumably the
193
genomic sequence -upstream bases of the gene,
194
including the gene, and -downstream bases of the gene
196
If you supply a -uid, and a map had previously been created and
197
given that uid, that same map object will be returned. Otherwise, the
198
combination of -gene and -species will be used to determine
199
if the same map had previously been made. If a corresponding map
200
hadn't previously been made, a new map object will be created and
206
my ($class, @args) = @_;
207
my ($uid, $gene, $species, $desc, $up, $down, $seq) = Bio::Root::Root->_rearrange([qw(UID
216
if ($uid && defined $GENEMAPS->{by_uid}->{$uid}) {
217
$gene_map = $GENEMAPS->{by_uid}->{$uid};
219
elsif ($gene && $species) {
220
my $name = ref($gene) ? $gene->universal_name : $gene;
221
if (defined $GENEMAPS->{by_ns}->{$name}->{$species}) {
222
$gene_map = $GENEMAPS->{by_ns}->{$name}->{$species};
226
$gene_map->gene->description($desc) if $desc;
227
$gene_map->upstream($up) if defined($up);
228
$gene_map->downstream($down) if defined($down);
229
$gene_map->seq($seq) if $seq;
233
return $class->new(@args);
239
Usage : my $id = $map->unique_id;
240
Function: Get/set the unique ID for this map
242
Args : none to get, OR string to set
247
my ($self, $id) = @_;
249
delete $GENEMAPS->{by_uid}->{$self->{'_uid'}};
250
$self->{'_uid'} = $id;
251
$GENEMAPS->{by_uid}->{$id} = $self;
253
return $self->{'_uid'};
259
Usage : my $species = $map->species;
260
Function: Get/set Species for a map. It is not recommended to change this once
262
Returns : Bio::Taxon object or string
263
Args : none to get, OR Bio::Taxon or string to set
268
my ($self, $value) = @_;
270
my $old_species = $self->{_species};
271
$self->{'_species'} = $value;
272
my $name = $self->universal_name || return $value;
274
delete $GENEMAPS->{by_ns}->{$name}->{$old_species};
276
$GENEMAPS->{by_ns}->{$name}->{$value} = $self;
278
return $self->{'_species'};
284
Usage : my $type = $map->type
285
Function: Get Map type
286
Returns : string 'gene'
298
Usage : my $gene = $map->gene;
299
$map->gene(-gene => $gene);
300
Function: Get/set the mappable element on this map that represents the gene
301
this map is for. Once set, it is not recommended to re-set the gene
302
to something else. Behaviour in that case is undefined.
303
Returns : Bio::Map::Gene
304
Args : none to get, OR to set:
305
-gene => Bio::Map::Gene or string of the universal name (see
306
Bio::Map::Gene docs), REQUIRED
307
-description => string, applied to the Bio::Map::Gene
308
-upstream => int, the number of bases the map extends before the
309
start of the gene element (default 1000).
310
-downstream => int, the number of bases the map extends beyond the
311
end of the gene element (default 0).
316
my ($self, @args) = @_;
319
my ($gene, $desc, $up, $down) = $self->_rearrange([qw(GENE
322
DOWNSTREAM)], @args);
323
$self->throw("You must supply -gene") unless $gene;
325
my $gene_obj = ref($gene) ? $gene : Bio::Map::Gene->get(-universal_name => $gene, -description => $desc);
326
if (defined $self->{gene}) {
327
if ($self->{gene} ne $gene_obj) {
328
$self->warn("Changing the gene that this map is for, which could be bad");
329
$self->purge_positions($self->{gene});
330
delete $GENEMAPS->{by_ns}->{$self->universal_name}->{$self->species};
331
$self->{gene} = $gene_obj;
334
# change the gene's position on us if necessary
335
$self->upstream($up) if defined $up;
336
$self->downstream($down) if defined $down;
339
# give the gene object a position on us
341
$up >= 0 || $self->throw("-upstream must be a positive integer");
342
Bio::Map::Position->new(-map => $self, -start => ($up + 1), -element => $gene_obj);
343
$self->{gene} = $gene_obj;
344
$self->downstream($down || 0);
346
# set other gene positions from db if already user-requested
347
$gene_obj->_set_from_db($self);
350
$GENEMAPS->{by_ns}->{$self->universal_name}->{$self->species} = $self;
353
return $self->{gene};
356
=head2 universal_name
358
Title : universal_name
359
Usage : my $name = $map->universal_name
360
Function: Get/set the name of Bio::Map::Gene object associated with this map.
361
It is not recommended to change this once set.
363
Args : none to get, OR string to set
368
my ($self, $value) = @_;
369
$self->gene || return;
371
my $species = $self->species;
372
delete $GENEMAPS->{by_ns}->{$self->gene->universal_name}->{$species};
373
$self->gene->universal_name($value);
374
$GENEMAPS->{by_ns}->{$value}->{$species} = $self;
376
return $self->gene->universal_name;
382
Usage : my $distance = $map->upstream;
383
$map->upstream($distance);
384
Function: Get/set how long the map is before the start of the Bio::Map::Gene
387
Args : none to get, OR int to set (the number of bases the map extends
388
before the start of the gene)
393
my ($self, $value) = @_;
395
my $pos = $self->gene->position($self);
396
if (defined($value)) {
397
$value >= 0 || $self->throw("Supplied value must be a positive integer");
398
$pos->start($value + 1);
401
return $pos->start - 1;
407
Usage : my $distance = $map->downstream;
408
$map->downstream($distance);
409
Function: Get/set the nominal end of the map relative to the end of the
410
Bio::Map::Gene object on this map.
412
Args : none to get, OR int to set (the number of bases the map extends
413
beyond the end of the gene)
419
if (@_) { $self->{_downstream} = shift }
420
return $self->{_downstream} || 0;
426
Usage : my $length = $map->length();
427
Function: Retrieves the length of the map. This is normally the length of the
428
upstream region + length of the gene + length of the downstream
429
region, but may be longer if positions have been placed on the map
430
beyond the end of the nominal downstream region.
438
my $expected_length = $self->gene->position($self)->length + $self->upstream + $self->downstream;
439
my $actual_length = $self->SUPER::length;
440
return $actual_length > $expected_length ? $actual_length : $expected_length;
446
Usage : $string = $obj->seq()
447
Function: Get/set the sequence as a string of letters. When getting, If the
448
GeneMap object didn't have sequence attached directly to it for the
449
region requested, the map's gene's database will be asked for the
450
sequence, and failing that, the map's gene's positions will be asked
451
for their sequences. Areas for which no sequence could be found will
452
be filled with Ns, unless no sequence was found anywhere, in which
453
case undef is returned.
455
Args : Optionally on set the new value (a string). An optional second
456
argument presets the alphabet (otherwise it will be guessed).
461
my ($self, @args) = @_;
462
my $seq = $self->SUPER::seq(@args);
463
my $expected_length = $self->length;
464
if (! $seq || CORE::length($seq) < $expected_length) {
465
my @have = split('', $seq || '');
467
for (0..($expected_length - 1)) {
468
$result[$_] = shift(@have) || 'N';
471
# build map sequence by asking gene or positions
472
my @slice_stuff = $self->gene->_get_slice($self);
474
my ($slice_adaptor, $slice, $strand) = @slice_stuff;
475
my ($start, $end, $gene_start) = (CORE::length($seq || '') + 1, $expected_length, $self->upstream + 1);
477
# convert map coords to genomic coords
478
my $adjust = $strand == -1 ? $slice->end : $slice->start;
479
my $adjustment = sub { return $strand == -1 ? $adjust - shift() : shift() + $adjust; };
480
my $converted_start = &$adjustment($start - $gene_start);
481
my $converted_end = &$adjustment($end - $gene_start);
482
($converted_start, $converted_end) = ($converted_end, $converted_start) if $converted_start > $converted_end;
484
# get sequence from a new slice of desired region
485
#*** what happens if desired region starts or ends off end of chromo?...
486
my $new_slice = $slice_adaptor->fetch_by_region($slice->coord_system_name, $slice->seq_region_name, $converted_start, $converted_end);
487
if ($new_slice && (my $seq_str = $new_slice->seq)) {
489
$seq_str = $self->_revcom($seq_str);
491
splice(@result, CORE::length($seq || ''), CORE::length($seq_str), split('', $seq_str));
495
foreach my $pos ($self->get_positions) {
496
next unless $pos->can('seq');
497
my @pos_seq = split('', $pos->seq(undef, undef, 1) || next);
498
for my $i ($pos->start($pos->absolute_relative)..$pos->end($pos->absolute_relative)) {
500
my $base = shift(@pos_seq);
501
if ($result[$i] eq 'N') {
508
$seq = join('', @result);
516
Usage : $substring = $obj->subseq(10, 40);
517
Function: Returns the subseq from start to end, where the first base
518
is 1 and the number is inclusive, ie 1-2 are the first two
519
bases of the sequence. If the GeneMap object didn't have sequence
520
attached directly to it for the region requested, the map's gene's
521
database will be asked for the sequence, and failing that, the map's
522
gene's positions will be asked for their sequences. Areas for which
523
no sequence could be found will be filled with Ns, unless no
524
sequence was found anywhere, in which case undef is returned. subseq
525
requests that extend beyond the end of the map will throw.
527
Args : integer for start position AND integer for end position
529
Bio::LocationI location for subseq (strand honored)
531
Bio::RangeI (eg. a Bio::Map::PositionI)
536
my ($self, $start, $end) = @_;
538
if ($start && ref($start) && $start->isa('Bio::RangeI')) {
540
if ($start->isa('Bio::Map::Position')) {
541
($start, $end) = ($thing->start($thing->absolute_relative), $thing->end($thing->absolute_relative));
544
($start, $end) = ($thing->start, $thing->end);
548
# *** this implementation potentially wastefull? Should duplicate code
549
# from seq() to do this just for the desired region??
550
my $orig_seq = $self->{seq};
551
$self->{seq} = $self->seq();
552
my $subseq = $self->{seq} ? $self->SUPER::subseq($start, $end) : '';
553
$self->{seq} = $orig_seq;
558
# quick revcom for strings (silly to create a PrimarySeq just to revcom and then
559
# return a string again)
561
my ($self, $seq) = @_;
563
$seq = reverse($seq);
564
$seq =~ tr/acgtrymkswhbvdnxACGTRYMKSWHBVDNX/tgcayrkmswdvbhnxTGCAYRKMSWDVBHNX/;
added
removed
Bio/Map/GeneMap.pm
