← Back to branch summary

~ubuntu-branches/ubuntu/trusty/bioperl/trusty

« back to all changes in this revision

Viewing changes to Bio/DB/Taxonomy/list.pm

  • Committer: Package Import Robot
  • Author(s): Charles Plessy
  • Date: 2013-09-22 13:39:48 UTC
  • mfrom: (3.1.11 sid)
  • Revision ID: package-import@ubuntu.com-20130922133948-c6z62zegjyp7ztou
Tags: 1.6.922-1
http://bugs.debian.org/722910
* New upstream release.
* Replaces and Breaks grinder (<< 0.5.3-3~) because of overlaping contents.
  Closes: #722910
* Stop Replacing and Breaking bioperl ( << 1.6.9 ): not needed anymore. 

Show diffs side-by-side

added added

removed removed

Lines of Context:
Bio/DB/Taxonomy/list.pm
11
11
 
12
12
# POD documentation - main docs before the code
13
13
 
 
14
 
14
15
=head1 NAME
15
16
 
16
17
Bio::DB::Taxonomy::list - An implementation of Bio::DB::Taxonomy
20
21
 
21
22
  use Bio::DB::Taxonomy;
22
23
 
 
24
  my $db = Bio::DB::Taxonomy->new( -source => 'list' );
 
25
 
 
26
  my @ranks = ('superkingdom', 'class', 'genus', 'species');
23
27
  my @names = ('Eukaryota', 'Mammalia', 'Homo', 'Homo sapiens');
24
 
  my @ranks = qw(superkingdom class genus species);
25
 
  my $db = Bio::DB::Taxonomy->new(-source => 'list', -names => \@names,
26
 
                                                    -ranks => \@ranks);
 
28
  $db->add_lineage(-names => \@names, -ranks => \@ranks);
27
29
 
28
30
  @names = ('Eukaryota', 'Mammalia', 'Mus', 'Mus musculus');
29
31
  $db->add_lineage(-names => \@names, -ranks => \@ranks);
83
85
 
84
86
# Let the code begin...
85
87
 
 
88
 
86
89
package Bio::DB::Taxonomy::list;
 
90
 
87
91
use strict;
88
92
use Bio::Taxon;
89
93
 
90
94
use base qw(Bio::DB::Taxonomy);
91
95
 
 
96
our $prefix = 'list';
 
97
 
 
98
 
92
99
=head2 new
93
100
 
94
101
 Title   : new
102
109
sub new {
103
110
    my ($class, @args) = @_;
104
111
    my $self = $class->SUPER::new(@args);
105
 
    $self->{db} = {};
106
 
    $self->add_lineage(@args) if @args;
 
112
    my %args = @args;
 
113
    delete $args{'-source'};
 
114
 
 
115
    $self->add_lineage(%args) if %args;
107
116
    
108
117
    return $self;
109
118
}
110
119
 
 
120
 
111
121
=head2 add_lineage
112
122
 
113
123
 Title   : add_lineage
114
 
 Usage   : $db->add_lineage(-names => \@names)
 
124
 Usage   : my @names = ('Eukaryota', 'Mammalia', 'Homo', 'Homo sapiens');
 
125
           my @ranks = ('superkingdom', 'class', 'genus', 'species');
 
126
           $db->add_lineage( -names => \@names, -ranks => \@ranks );
115
127
 Function: Add a lineage to the database, where the lineage is described by
116
128
           a list of scientific names in the order root->leaf. The rank of each
117
129
           name can optionally be described by supplying an additional list
118
130
           of rank names in the same order (eg. superkingdom->species).
119
 
 Returns : n/a
 
131
 Returns : 1 for success
120
132
 Args    : -names => [] : array ref of scientific names, REQUIRED
121
133
           -ranks => [] : array ref of rank names, same order as above, OPTIONAL
122
134
 
123
135
=cut
124
136
 
125
137
sub add_lineage {
126
 
    my $self = shift;
127
 
    my ($names, $ranks) = $self->_rearrange([qw (NAMES RANKS)], @_);
128
 
    $self->throw("-names must be supplied and its value must be an array reference") unless $names && ref($names) eq 'ARRAY';
129
 
    my @names = @{$names};
130
 
    
131
 
    my @ranks;
 
138
    my ($self, @args) = @_;
 
139
    my ($names, $ranks) = $self->_rearrange([qw (NAMES RANKS)], @args);
 
140
    $self->throw("-names must be supplied and its value must be an array reference")
 
141
        unless $names && ref($names) eq 'ARRAY';
 
142
 
 
143
    my $names_idx = scalar @$names - 1;
 
144
 
132
145
    if ($ranks) {
133
 
        $self->throw("-ranks must be an array reference") unless ref($ranks) eq 'ARRAY';
134
 
        $self->throw("The -names and -ranks lists must be of equal length") unless @{$names} == @{$ranks};
135
 
        @ranks = @{$ranks};
136
 
    }
137
 
    else {
138
 
        for (0..$#names) {
139
 
            push(@ranks, 'no rank');
140
 
        }
 
146
        $self->throw("-ranks must be an array reference")
 
147
            unless ref($ranks) eq 'ARRAY';
 
148
        $self->throw("The -names and -ranks lists must be of equal length")
 
149
            unless $names_idx == scalar @$ranks - 1;
141
150
    }
142
151
    
143
152
    # This is non-trivial because names are not guaranteed unique in a taxonomy,
153
162
    # ('Mammalia', 'Hominidae', 'Homo', 'Homo sapiens')
154
163
    #
155
164
    # Clearly with limited information we can't do a perfect job, but we can try
156
 
    # and do a reasonable one.
157
 
    
158
 
    
159
 
    #...
160
 
    
161
 
    
162
 
    # All that said, let's just do the trivial implementation now and see how
163
 
    # bad it is! (assumes ranks are unique except for 'no rank')
164
 
    
165
 
    
166
 
    my $first_lineage = $self->{db}->{node_ids} ? 0 : 1;
167
 
    
168
 
    my $ancestor_node_id;
 
165
    # and do a reasonable one. So, let's just do the trivial implementation now
 
166
    # and see how bad it is! (assumes ranks are unique except for 'no rank')
 
167
    
 
168
    my $ancestors  = $self->{ancestors};
 
169
    my $node_data  = $self->{node_data};
 
170
    my $name_to_id = $self->{name_to_id};
 
171
    my $children   = $self->{children};
 
172
 
 
173
    my $my_ancestor_id = '';
169
174
    my @node_ids;
170
 
    for my $i (0..$#names) {
171
 
        my $name = $names[$i];
172
 
        my $rank = $ranks[$i];
173
 
        
 
175
    for my $i (0 .. $names_idx) {
 
176
        my $name = $names->[$i];
 
177
        my $rank = $ranks->[$i]; # if undef, this node has 'no rank'
 
178
 
174
179
        # This is a new node with a new id if we haven't seen this name before.
175
180
        # It's also always a new node if this is the first lineage going into
176
181
        # the db.
177
182
        #
178
183
        # We need to handle, however, situations in the future where we try to
179
184
        # merge in a new lineage but we have non-unique names in the lineage
180
 
        # and possible missing classes in some lineages
181
 
        # (eg.
182
 
        # '... Anophelinae, Anopheles, Anopheles, Angusticorn, Anopheles...'
 
185
        # and possible missing classes in some lineages, e.g.
 
186
        #    '... Anophelinae, Anopheles, Anopheles, Angusticorn, Anopheles...'
183
187
        # merged with
184
 
        # '... Anophelinae, Anopheles, Angusticorn, Anopheles...'),
 
188
        #    '... Anophelinae, Anopheles, Angusticorn, Anopheles...'),
185
189
        # but still need the 'tree' to be correct
186
 
        
187
 
        my $is_new = 0;
188
 
        if ($first_lineage || ! exists $self->{db}->{name_to_id}->{$name}) {
189
 
            $is_new = 1;
190
 
        }
191
 
        
 
190
 
 
191
        # Look for a node that is consistent with this lineage
192
192
        my $node_id;
193
 
        unless ($is_new) {
194
 
            my @same_named = @{$self->{db}->{name_to_id}->{$name}};
 
193
        SAME_NAMED: for my $same_id (@{$name_to_id->{$name}}) {
 
194
 
 
195
            # Taxa are the same if it they have the same ancestor or none
 
196
            my $this_ancestor_id = $ancestors->{$same_id} || '';
 
197
            if ($my_ancestor_id eq $this_ancestor_id) {
 
198
                $node_id = $same_id;
 
199
                last SAME_NAMED;
 
200
            }
195
201
            
196
 
            # look for the node that is consistent with this lineage
197
 
            SAME_NAMED: foreach my $s_id (@same_named) {
198
 
                my $this_ancestor_id;
199
 
                if ($ancestor_node_id) {
200
 
                    $this_ancestor_id = $self->{db}->{ancestors}->{$s_id};
201
 
                    if ($ancestor_node_id eq $this_ancestor_id) {
202
 
                        $node_id = $s_id;
 
202
            # Compare children
 
203
            next if $i >= $names_idx; # this taxon has no child
 
204
            my $my_child_name = $names->[$i + 1];
 
205
            #while ( my ($this_child_id, undef) = each %{$children->{$same_id}} ) {
 
206
            for my $this_child_id (keys %{$children->{$same_id}}) {
 
207
                if ($my_child_name eq $node_data->{$this_child_id}->[0]) { # both children have same name
 
208
                    if ($my_ancestor_id) {
 
209
                        my @s_ancestors;
 
210
                        while ($this_ancestor_id = $ancestors->{$this_ancestor_id}) {
 
211
                            if ($my_ancestor_id eq $this_ancestor_id) {
 
212
                                $my_ancestor_id = $ancestors->{$same_id};
 
213
                                push @node_ids, @s_ancestors, $my_ancestor_id;
 
214
                                $node_id = $same_id;
 
215
                                last SAME_NAMED;
 
216
                            }
 
217
                            unshift @s_ancestors, $this_ancestor_id;
 
218
                        }
 
219
                    } else {
 
220
                        # This new lineage (@$names) doesn't start at the
 
221
                        # same root as the existing lineages. Assuming
 
222
                        # '$name' corresponds to node $same_id");
 
223
                        $node_id = $same_id;
203
224
                        last SAME_NAMED;
204
225
                    }
205
226
                }
206
 
                
207
 
                if ($names[$i + 1]) {
208
 
                    my $my_child_name = $names[$i + 1];
209
 
                    my @children_ids = keys %{$self->{db}->{children}->{$s_id} || {}};
210
 
                    foreach my $c_id (@children_ids) {
211
 
                        my $this_child_name = $self->{db}->{node_data}->{$c_id}->[0];
212
 
                        if ($my_child_name eq $this_child_name) {
213
 
                            
214
 
                            if ($ancestor_node_id) {
215
 
                                my @s_ancestors;
216
 
                                while ($this_ancestor_id = $self->{db}->{ancestors}->{$this_ancestor_id}) {
217
 
                                    if ($ancestor_node_id eq $this_ancestor_id) {
218
 
                                        $node_id = $s_id;
219
 
                                        $ancestor_node_id = $self->{db}->{ancestors}->{$s_id};
220
 
                                        push(@node_ids, @s_ancestors, $ancestor_node_id);
221
 
                                        last SAME_NAMED;
222
 
                                    }
223
 
                                    unshift(@s_ancestors, $this_ancestor_id);
224
 
                                }
225
 
                            }
226
 
                            else {
227
 
                                #$self->warn("This new lineage (@names) doesn't start at the same root as the existing lineages.".
228
 
                                #            "\nI'm assuming '$name' corresponds to node $s_id");
229
 
                                $node_id = $s_id;
230
 
                                last SAME_NAMED;
231
 
                            }
232
 
                        }
233
 
                    }
234
 
                }
235
 
            }
236
 
            
237
 
            $node_id || $is_new++;
238
 
        }
239
 
        
240
 
        if ($is_new) {
241
 
            my $next_num = ++$self->{db}->{node_ids};
242
 
            # 'list' so definitely not confused with ncbi taxonomy ids
243
 
            $node_id = 'list'.$next_num;
244
 
            push(@{$self->{db}->{name_to_id}->{$name}}, $node_id);
245
 
        }
246
 
        
247
 
        unless (exists $self->{db}->{node_data}->{$node_id}) {
248
 
            $self->{db}->{node_data}->{$node_id} = [($name, '')];
249
 
        }
250
 
        my $node_data = $self->{db}->{node_data}->{$node_id};
251
 
        
252
 
        if (!$node_data->[1] || ($node_data->[1] eq 'no rank' && $rank ne 'no rank')) {
253
 
            $node_data->[1] = $rank;
254
 
        }
255
 
        
256
 
        if ($ancestor_node_id) {
257
 
            if ($self->{db}->{ancestors}->{$node_id} && $self->{db}->{ancestors}->{$node_id} ne $ancestor_node_id) {
258
 
                $self->throw("This lineage (".join(', ', @names).") and a previously computed lineage share a node name but have different ancestries for that node. Can't cope!");
259
 
            }
260
 
            $self->{db}->{ancestors}->{$node_id} = $ancestor_node_id;
261
 
        }
262
 
        
263
 
        $ancestor_node_id = $node_id;
264
 
        push(@node_ids, $node_id);
 
227
            }
 
228
        }
 
229
        
 
230
        if (not defined $node_id) {
 
231
            # This is a new node. Add it to the database, using the prefix 'list'
 
232
            # for its ID to distinguish it from the IDs from other taxonomies.
 
233
            my $next_num = ++$self->{node_ids};
 
234
            $node_id = $prefix.$next_num;
 
235
            push @{$self->{name_to_id}->{$name}}, $node_id;
 
236
            $self->{node_data}->{$node_id}->[0] = $name;
 
237
        }
 
238
 
 
239
        if ( (defined $rank) && (not defined $node_data->{$node_id}->[1]) ) {
 
240
            # Save rank if node in database has no rank but the current node has one
 
241
            $self->{node_data}->{$node_id}->[1] = $rank;
 
242
        }
 
243
 
 
244
        if ($my_ancestor_id) {
 
245
            if ($self->{ancestors}->{$node_id} && $self->{ancestors}->{$node_id} ne $my_ancestor_id) {
 
246
                $self->throw("The lineage '".join(', ', @$names)."' and a ".
 
247
                    "previously stored lineage share a node name but have ".
 
248
                    "different ancestries for that node. Can't cope!");
 
249
            }
 
250
            $self->{ancestors}->{$node_id} = $my_ancestor_id;
 
251
        }
 
252
        
 
253
        $my_ancestor_id = $node_id;
 
254
        push @node_ids, $node_id;
265
255
    }
266
256
    
267
 
    # go through the lineage in reverse so we can remember the children
268
 
    my $child_id;
269
 
    foreach my $node_id (reverse @node_ids) {
270
 
        unless ($child_id) {
271
 
            $child_id = $node_id;
272
 
            next;
273
 
        }
274
 
        
275
 
        $self->{db}->{children}->{$node_id}->{$child_id} = 1;
276
 
        $child_id = $node_id;
 
257
    # Go through the lineage in reverse so we can remember the children
 
258
    for (my $i = $names_idx - 1; $i >= 0; $i--) {
 
259
        $self->{children}->{$node_ids[$i]}->{$node_ids[$i+1]} = undef;
277
260
    }
 
261
    return 1;
278
262
}
279
263
 
 
264
 
280
265
=head2 Bio::DB::Taxonomy Interface implementation
281
266
 
 
267
=head2 get_num_taxa
 
268
 
 
269
 Title   : get_num_taxa
 
270
 Usage   : my $num = $db->get_num_taxa();
 
271
 Function: Get the number of taxa stored in the database.
 
272
 Returns : A number
 
273
 Args    : None
 
274
 
282
275
=cut
283
276
 
 
277
sub get_num_taxa {
 
278
    my ($self) = @_;
 
279
    return $self->{node_ids} || 0;
 
280
}
 
281
 
 
282
 
284
283
=head2 get_taxon
285
284
 
286
285
 Title   : get_taxon
287
286
 Usage   : my $taxon = $db->get_taxon(-taxonid => $taxonid)
288
287
 Function: Get a Bio::Taxon object from the database.
289
288
 Returns : Bio::Taxon object
290
 
 Args    : just a single value which is the database id, OR named args:
291
 
           -taxonid => taxonomy id (to query by taxonid; NB: these are not
292
 
                       NCBI taxonomy ids but 'list' pre-fixed ids unique to the
293
 
                       list database)
294
 
            OR
295
 
           -name    => string (to query by a taxonomy name)
 
289
 Args    : A single value which is the ID of the taxon to retrieve
 
290
             OR named args, as follows:
 
291
           -taxonid => Taxonomy ID (NB: these are not NCBI taxonomy ids but
 
292
                       'list' pre-fixed ids unique to the list database).
 
293
             OR
 
294
           -name    => String (to query by a taxonomy name). A given taxon name
 
295
                       can match different taxonomy objects. When that is the
 
296
                       case, a warning is displayed and the first matching taxon
 
297
                       is reported. See get_taxonids() to get all matching taxon
 
298
                       IDs.
 
299
             OR
 
300
           -names   => Array ref of lineage names, like in add_lineage(). To
 
301
                       overcome the limitations of -name, you can use -names to
 
302
                       provide the full lineage of the taxon you want and get a
 
303
                       unique, unambiguous taxon object.
296
304
 
297
305
=cut
298
306
 
299
307
sub get_taxon {
300
 
    my $self = shift;
301
 
    my ($taxonid, $name);
302
 
    
303
 
    if (@_ > 1) {
304
 
        ($taxonid, $name) = $self->_rearrange([qw(TAXONID NAME)],@_);
 
308
    my ($self, @args) = @_;
 
309
 
 
310
    my $taxonid;
 
311
    if (scalar @args == 1) {
 
312
        # Argument is a taxon ID
 
313
        $taxonid = $args[0];
 
314
    } else {
 
315
        # Got named arguments
 
316
        my ($name, $names);
 
317
        ($taxonid, $name, $names) = $self->_rearrange([qw(TAXONID NAME NAMES)], @args);
305
318
        if ($name) {
306
 
            ($taxonid, my @others) = $self->get_taxonids($name);
307
 
            $self->warn("There were multiple ids ($taxonid @others) matching '$name', using '$taxonid'") if @others > 0;
308
 
        }
309
 
    }
310
 
    else {
311
 
        $taxonid = shift;
312
 
    }
313
 
    
314
 
    my $node = $self->{db}->{node_data}->{$taxonid} || return;
315
 
    my ($sci_name, $rank) = @{$node};
316
 
    
317
 
    my $taxon = Bio::Taxon->new(
318
 
                        -name         => $sci_name,
319
 
                        -object_id    => $taxonid, # since this is NOT a real ncbi taxid, set it as simply the object id
320
 
                        -rank         => $rank );
321
 
    # we can't use -dbh or the db_handle() method ourselves or we'll go
322
 
    # infinite on the merge attempt
323
 
    $taxon->{'db_handle'} = $self;
324
 
    
325
 
    $self->_handle_internal_id($taxon, 1);
326
 
    
 
319
            $names = [$name];
 
320
        }
 
321
        if ($names) {
 
322
            $name = $names->[-1];
 
323
 
 
324
            my @taxonids = $self->get_taxonids($name);
 
325
            $taxonid = $taxonids[0];
 
326
 
 
327
            # Use provided lineage to find correct ID amongst several matching IDs
 
328
            if ( (scalar @taxonids > 1) && (scalar @$names > 1) ) {
 
329
                for my $query_taxonid (@taxonids) {
 
330
                    my $matched = 1;
 
331
                    my $db_ancestor = $self->get_taxon($query_taxonid);
 
332
                    for (my $i = $#$names-1; $i >= 0; $i--) {
 
333
                        my $query_ancestor_name = $names->[$i];
 
334
                        $db_ancestor = $db_ancestor->ancestor;
 
335
                        my $db_ancestor_name = '';
 
336
                        if ($db_ancestor) {
 
337
                            $db_ancestor_name = $db_ancestor->node_name;
 
338
                        }
 
339
                        if (not ($query_ancestor_name eq $db_ancestor_name) ) {
 
340
                            $matched = 0;
 
341
                            last; # done testing this taxonid
 
342
                        }
 
343
                    }
 
344
                    if ($matched == 1) {
 
345
                       @taxonids = [$query_taxonid];
 
346
                       $taxonid  =  $query_taxonid;
 
347
                       last; # done testing all taxonids
 
348
                    }
 
349
                }
 
350
            }
 
351
 
 
352
            # Warn if several taxon IDs matched
 
353
            if (scalar @taxonids > 1) {
 
354
                $self->warn("There were multiple ids (@taxonids) matching '$name',".
 
355
                    " using '$taxonid'") if scalar @taxonids > 1;
 
356
            }
 
357
 
 
358
        }
 
359
    }
 
360
    
 
361
    # Now that we have the taxon ID, retrieve the corresponding Taxon object
 
362
    my $taxon;
 
363
    my $node = $self->{node_data}->{$taxonid};
 
364
    if ($node) {
 
365
        my ($sci_name, $rank) = @$node;
 
366
        $taxon = Bio::Taxon->new(
 
367
            -name      => $sci_name,
 
368
            -object_id => $taxonid, # not an ncbi taxid, simply an object id
 
369
        );
 
370
 
 
371
        if ($rank) {
 
372
            $taxon->rank($rank);
 
373
        }
 
374
 
 
375
        # we can't use -dbh or the db_handle() method ourselves or we'll go
 
376
        # infinite on the merge attempt
 
377
        $taxon->{'db_handle'} = $self;
 
378
    
 
379
        $self->_handle_internal_id($taxon, 1);
 
380
    }
 
381
 
327
382
    return $taxon;
328
383
}
329
384
 
330
385
*get_Taxonomy_Node = \&get_taxon;
331
386
 
 
387
 
332
388
=head2 get_taxonids
333
389
 
334
390
 Title   : get_taxonids
342
398
=cut
343
399
 
344
400
sub get_taxonids {
345
 
    my ($self, $query) = @_;
346
 
    return @{$self->{db}->{name_to_id}->{$query} || []};
 
401
    my ($self, $name) = @_;
 
402
    return wantarray() ? @{$self->{name_to_id}->{$name} || []} : $self->{name_to_id}->{$name}->[0];
347
403
}
348
404
 
349
405
*get_taxonid = \&get_taxonids;
350
406
 
 
407
 
351
408
=head2 ancestor
352
409
 
353
410
 Title   : ancestor
363
420
    my ($self, $taxon) = @_;
364
421
    $taxon || return; # for bug 2092, or something similar to it at least: shouldn't need this!
365
422
    $self->throw("Must supply a Bio::Taxon") unless ref($taxon) && $taxon->isa('Bio::Taxon');
366
 
    $self->throw("The supplied Taxon must belong to this database") unless $taxon->db_handle && $taxon->db_handle eq $self;
 
423
    $self->throw("The supplied Taxon must belong to this database")
 
424
        unless $taxon->db_handle && $taxon->db_handle eq $self;
367
425
    my $id = $taxon->id || $self->throw("The supplied Taxon is missing its id!");
368
426
    
369
 
    my $ancestor_id = $self->{db}->{ancestors}->{$id} || return;
 
427
    my $ancestor_id = $self->{ancestors}->{$id} || return;
370
428
    return $self->get_taxon($ancestor_id);
371
429
}
372
430
 
 
431
 
373
432
=head2 each_Descendent
374
433
 
375
434
 Title   : each_Descendent
384
443
sub each_Descendent {
385
444
    my ($self, $taxon) = @_;
386
445
    $self->throw("Must supply a Bio::Taxon") unless ref($taxon) && $taxon->isa('Bio::Taxon');
387
 
    $self->throw("The supplied Taxon must belong to this database") unless $taxon->db_handle && $taxon->db_handle eq $self;
 
446
    $self->throw("The supplied Taxon must belong to this database")
 
447
        unless $taxon->db_handle && $taxon->db_handle eq $self;
388
448
    my $id = $taxon->id || $self->throw("The supplied Taxon is missing its id!");
389
449
    
390
 
    my @children_ids = keys %{$self->{db}->{children}->{$id} || {}};
391
450
    my @children;
392
 
    foreach my $child_id (@children_ids) {
393
 
        push(@children, $self->get_taxon($child_id) || next);
 
451
    while ( my ($child_id, undef) = each %{$self->{children}->{$id}} ) {
 
452
        push @children, ($self->get_taxon($child_id) || next);
394
453
    }
395
454
    
396
455
    return @children;
397
456
}
398
457
 
 
458
 
399
459
1;