[Bio] / FigKernelPackages / FIGV.pm Repository:
ViewVC logotype

Annotation of /FigKernelPackages/FIGV.pm

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1.17 - (view) (download) (as text)

1 : overbeek 1.1 #
2 :     # Copyright (c) 2003-2006 University of Chicago and Fellowship
3 :     # for Interpretations of Genomes. All Rights Reserved.
4 :     #
5 :     # This file is part of the SEED Toolkit.
6 :     #
7 :     # The SEED Toolkit is free software. You can redistribute
8 :     # it and/or modify it under the terms of the SEED Toolkit
9 :     # Public License.
10 :     #
11 :     # You should have received a copy of the SEED Toolkit Public License
12 :     # along with this program; if not write to the University of Chicago
13 :     # at info@ci.uchicago.edu or the Fellowship for Interpretation of
14 :     # Genomes at veronika@thefig.info or download a copy from
15 :     # http://www.theseed.org/LICENSE.TXT.
16 :     #
17 :    
18 :     package FIGV;
19 :    
20 : overbeek 1.6 use Carp;
21 : overbeek 1.1 use strict;
22 :     use FIG;
23 :     use FIG_Config;
24 :     use SFXlate;
25 :     use SproutFIG;
26 :     use Tracer;
27 :     use Data::Dumper;
28 : olson 1.4 use vars qw($AUTOLOAD);
29 : olson 1.13 use DB_File;
30 :     use FileHandle;
31 : overbeek 1.1
32 :     sub new {
33 :     my($class,$org_dir,$low_level) = @_;
34 :    
35 :     my $fig;
36 :     if ($low_level && ($low_level =~ /sprout/i))
37 :     {
38 :     $fig = new SproutFIG($FIG_Config::sproutDB, $FIG_Config::sproutData);
39 :     }
40 : overbeek 1.17 elsif (! $org_dir)
41 :     {
42 :     $fig = new FIG;
43 :     return $fig;
44 :     }
45 : overbeek 1.1 else
46 :     {
47 :     $fig = new FIG;
48 :     }
49 :    
50 :     my $self = {};
51 :     $self->{_fig} = $fig;
52 :     $self->{_orgdir} = $org_dir;
53 :    
54 :     $org_dir =~ /(\d+\.\d+)$/;
55 :     $self->{_genome} = $1;
56 :     return bless $self, $class;
57 :     }
58 :    
59 : olson 1.4 #
60 :     # Redirect any method invocations that we don't implement out to the
61 :     # underlying FIG object.
62 :     #
63 :     sub AUTOLOAD
64 :     {
65 :     my($self, @args) = @_;
66 :    
67 : overbeek 1.6 if (ref($self) ne "FIGV") {
68 :     confess "BAD FIGV object passed to AUTOLOAD";
69 :     }
70 :    
71 : olson 1.4 my $meth = $AUTOLOAD;
72 :     $meth =~ s/.*:://;
73 :    
74 :     my $fig = $self->{_fig};
75 : olson 1.8 # my $args = Dumper(\@args);
76 : olson 1.4 if (wantarray)
77 :     {
78 :     my @res = $fig->$meth(@args);
79 : olson 1.5 # warn "FIGV invoke $meth($args) returns\n", Dumper(\@res);
80 : olson 1.4 return @res;
81 :     }
82 :     else
83 :     {
84 :     my $res = $fig->$meth(@args);
85 : olson 1.5 # warn "FIGV invoke $meth($args) returns\n", Dumper($res);
86 : olson 1.4 return $res;
87 :     }
88 :     }
89 :    
90 : olson 1.8 sub FIG
91 :     {
92 :     my($self) = @_;
93 :     return $self;
94 :     }
95 :    
96 :    
97 : olson 1.5 #
98 :     # To retrieve a subsystem in FIGV, we create the subsystem as normal via $fig->get_subsystem,
99 :     # then insert the row for the virtual org dir we are processing.
100 :     #
101 :    
102 :     sub get_subsystem
103 :     {
104 :     my($self,$ssa) = @_;
105 :    
106 :     my $fig = $self->{_fig};
107 :     my $newG = $self->{_genome};
108 :     my $newGdir = $self->{_orgdir};
109 :    
110 :     my $ss = $fig->get_subsystem($ssa);
111 :     return undef unless $ss;
112 :    
113 : olson 1.9 $self->load_ss_data();
114 : olson 1.8
115 : olson 1.9 my $bindings = $self->{_ss_bindings}->{$ssa};
116 :     my $variant = $self->{_ss_variants}->{$ssa};
117 :    
118 :     # warn "Adding virtual genome " . Dumper(\%bindings);
119 :     $ss->add_virtual_genome($self->genus_species(), $newG, $variant, $bindings);
120 :    
121 :     return $ss;
122 :     }
123 :    
124 : olson 1.14 sub active_subsystems
125 :     {
126 :     my($self, $genome, $all) = @_;
127 :    
128 :     my $fig = $self->{_fig};
129 :     my $newG = $self->{_genome};
130 :     my $newGdir = $self->{_orgdir};
131 :    
132 :     if ($genome ne $newG)
133 :     {
134 :     return $fig->active_subsystems($genome, $all);
135 :     }
136 :    
137 :     $self->load_ss_data();
138 :    
139 :     my $slist = {};
140 :     %{$slist} = %{$self->{_ss_variants}};
141 :    
142 :     if (not $all)
143 :     {
144 :     for my $ss (keys %$slist)
145 :     {
146 :     my $var = $slist->{$ss};
147 :     delete $slist->{$ss} if $var == 0 or $var == -1;
148 :     }
149 :     }
150 :     return $slist;
151 :     }
152 :    
153 : olson 1.9 sub peg_to_roles_in_subsystems
154 :     {
155 :     my($self,$peg) = @_;
156 :    
157 :     my $fig = $self->{_fig};
158 :     my $newG = $self->{_genome};
159 :     my $newGdir = $self->{_orgdir};
160 : olson 1.8
161 : olson 1.9 $self->load_ss_data();
162 : olson 1.8
163 : olson 1.9 my $ret = $self->{_ss_peg_index}->{$peg};
164 : olson 1.8
165 : olson 1.9 return $ret ? @$ret : ();
166 :     }
167 : olson 1.8
168 : olson 1.9 sub subsystems_for_peg
169 :     {
170 :     my($self,$peg) = @_;
171 :     return $self->peg_to_roles_in_subsystems($peg);
172 : olson 1.5 }
173 : olson 1.4
174 : overbeek 1.1 sub genomes {
175 :     my($self,$complete) = @_;
176 :     my $fig = $self->{_fig};
177 :    
178 :     return ($self->{_genome},$fig->genomes($complete));
179 :     }
180 :    
181 :     sub genus_species {
182 :     my($self,$genome) = @_;
183 :    
184 :     my $fig = $self->{_fig};
185 :     my $newG = $self->{_genome};
186 :     my $newGdir = $self->{_orgdir};
187 :    
188 :     if (($genome eq $newG) && open(GENOME,"<$newGdir/GENOME"))
189 :     {
190 :     my $x = <GENOME>;
191 :     close(GENOME);
192 :     chop $x;
193 :     return $x;
194 :     }
195 :     else
196 :     {
197 :     return $fig->genus_species($genome);
198 :     }
199 :     }
200 :    
201 : overbeek 1.7 sub get_genome_assignment_data {
202 :     my($self,$genome) = @_;
203 :    
204 :     my $fig = $self->{_fig};
205 :     my $newG = $self->{_genome};
206 :     my $newGdir = $self->{_orgdir};
207 :    
208 :     if ($genome eq $newG)
209 :     {
210 :     my %assign = map { ( $_ =~ /^(fig\|\d+\.\d+\.peg\.\d+)\t(\S.*\S)/) ? ($1 => $2) : () }
211 :     `cat $newGdir/proposed*functions`;
212 :     return [map { [$_,$assign{$_}] } sort { &FIG::by_fig_id($a,$b) } keys(%assign)];
213 :     }
214 :     else
215 :     {
216 :     return $fig->get_genome_assignment_data($genome);
217 :     }
218 :     }
219 :    
220 :     sub org_of {
221 :     my($self,$peg) = @_;
222 :    
223 :     if ($peg =~ /^fig\|(\d+\.\d+)\.peg\.\d+/)
224 :     {
225 :     return $self->genus_species($1);
226 :     }
227 :     return "";
228 :     }
229 :    
230 :     sub get_genome_subsystem_data {
231 :     my($self,$genome) = @_;
232 :    
233 :     my $fig = $self->{_fig};
234 :     my $newG = $self->{_genome};
235 :     my $newGdir = $self->{_orgdir};
236 :    
237 :     if ($genome eq $newG)
238 :     {
239 :     return [map { ($_ =~ /^(\S[^\t]+\S)\t(\S[^\t]*\S)\t(\S+)/) ? [$1,$2,$3] : () }
240 :     `cat $newGdir/Subsystems/bindings`];
241 :     }
242 :     else
243 :     {
244 :     return $fig->get_genome_subsystem_data($genome);
245 :     }
246 :     }
247 :    
248 :     sub orgname_of_orgid {
249 :     my($self,$genome) = @_;
250 :    
251 :     return $self->genus_species($genome);
252 :     }
253 :    
254 :     sub orgid_of_orgname {
255 :     my($self,$genome_name) = @_;
256 :    
257 :     my @genomes = $self->genomes('complete');
258 :     my $i;
259 :     for ($i=0; ($i < @genomes) && ($genome_name ne $self->genus_species($genomes[$i])); $i++) {}
260 :     return ($i == @genomes) ? undef : $genomes[$i];
261 :     }
262 :    
263 :     sub genus_species_domain {
264 :     my($self,$genome) = @_;
265 :    
266 :     return [$self->genus_species($genome),$self->genome_domain($genome)];
267 :     }
268 :    
269 :     sub protein_subsystem_to_roles {
270 :     my ($self,$peg,$subsystem) = @_;
271 :    
272 :     my $fig = $self->{_fig};
273 :     my $newG = $self->{_genome};
274 :     my $newGdir = $self->{_orgdir};
275 :    
276 :     if (&FIG::genome_of($peg) ne $newG)
277 :     {
278 :     return $fig->protein_subsystem_to_roles($peg,$subsystem);
279 :     }
280 :     else
281 :     {
282 :     my @roles = map { (($_ =~ /^([^\t]+)\t([^\t]+)\t(\S+)$/) && ($1 eq $subsystem) && ($3 eq $peg)) ?
283 :     $2 : () } `cat $newGdir/Subsystems/bindings`;
284 :     my %roles = map { $_ => 1 } @roles;
285 :     return [sort keys(%roles)];
286 :     }
287 :     }
288 :    
289 :     sub contig_lengths {
290 :     my ($self, $genome) = @_;
291 :    
292 :     my $fig = $self->{_fig};
293 :     my $newG = $self->{_genome};
294 :     my $newGdir = $self->{_orgdir};
295 :    
296 :     if ($genome ne $newG)
297 :     {
298 :     return $fig->contig_lengths($genome);
299 :     }
300 :     else
301 :     {
302 :     my $contig_lengths = {};
303 :     if (open(CONTIGS,"<$newGdir/contigs"))
304 :     {
305 :     $/ = "\n>";
306 :     while (defined(my $x = <CONTIGS>))
307 :     {
308 :     if ($x =~ />?(\S+)[^\n]*\n(.*)/s)
309 :     {
310 :     my $id = $1;
311 :     my $seq = $2;
312 :     $seq =~ s/\s//gs;
313 :     $contig_lengths->{$id} = length($seq);
314 :     }
315 :     }
316 :     close(CONTIGS);
317 :     $/ = "\n";
318 :     }
319 :     return $contig_lengths;
320 :     }
321 :     }
322 :    
323 : olson 1.10 sub contig_ln {
324 :     my ($self, $genome, $contig) = @_;
325 :    
326 :     my $fig = $self->{_fig};
327 :     my $newG = $self->{_genome};
328 :     my $newGdir = $self->{_orgdir};
329 :    
330 :     if ($genome ne $newG)
331 :     {
332 :     return $fig->contig_ln($genome, $contig);
333 :     }
334 :     else
335 :     {
336 :     if (open(CONTIGS,"<$newGdir/contigs"))
337 :     {
338 :     local $/ = "\n>";
339 :     while (defined(my $x = <CONTIGS>))
340 :     {
341 :     if ($x =~ />?(\S+)[^\n]*\n(.*)/s)
342 :     {
343 :     my $id = $1;
344 :     my $seq = $2;
345 :     $seq =~ s/\s//gs;
346 :     if ($id eq $contig)
347 :     {
348 :     return length($seq);
349 :     }
350 :     }
351 :     }
352 :     close(CONTIGS);
353 :     }
354 :     }
355 :     }
356 :    
357 :     sub contigs_of
358 :     {
359 :     my ($self, $genome) = @_;
360 :    
361 :     my $fig = $self->{_fig};
362 :     my $newG = $self->{_genome};
363 :     my $newGdir = $self->{_orgdir};
364 :    
365 :     if ($genome ne $newG)
366 :     {
367 :     return $fig->contigs_of($genome);
368 :     }
369 :     else
370 :     {
371 :     my @out;
372 :     if (open(CONTIGS,"<$newGdir/contigs"))
373 :     {
374 :     local $/ = "\n>";
375 :     while (defined(my $x = <CONTIGS>))
376 :     {
377 :     if ($x =~ />?(\S+)[^\n]*\n(.*)/s)
378 :     {
379 :     my $id = $1;
380 :     push(@out, $id);
381 :     }
382 :     }
383 :     close(CONTIGS);
384 :     }
385 :     return @out;
386 :     }
387 :     }
388 :    
389 :     =head3 dna_seq
390 :    
391 :     usage: $seq = dna_seq($genome,@locations)
392 :    
393 :     Returns the concatenated subsequences described by the list of locations. Each location
394 :     must be of the form
395 :    
396 :     Contig_Beg_End
397 :    
398 :     where Contig must be the ID of a contig for genome $genome. If Beg > End the location
399 :     describes a stretch of the complementary strand.
400 :    
401 :     =cut
402 :     #: Return Type $;
403 :     sub dna_seq {
404 :     my($self,$genome,@locations) = @_;
405 :    
406 :     my $fig = $self->{_fig};
407 :     my $newG = $self->{_genome};
408 :     my $newGdir = $self->{_orgdir};
409 :    
410 :     if ($genome ne $newG)
411 :     {
412 :     return $fig->dna_seq($genome, @locations);
413 :     }
414 :    
415 :     my %contigs;
416 :     if (open(CONTIGS,"<$newGdir/contigs"))
417 :     {
418 :     local $/ = "\n>";
419 :     while (defined(my $x = <CONTIGS>))
420 :     {
421 :     if ($x =~ />?(\S+)[^\n]*\n(.*)/s)
422 :     {
423 :     my $id = $1;
424 :     my $seq = $2;
425 :     $seq =~ s/\s//gs;
426 :     $contigs{$id} = $seq;
427 :     }
428 :     }
429 : paczian 1.15 close(CONTIGS);
430 : olson 1.10 }
431 :    
432 :     my(@pieces,$loc,$contig,$beg,$end,$ln,$rdbH);
433 :    
434 :     @locations = map { split(/,/,$_) } @locations;
435 :     @pieces = ();
436 :     foreach $loc (@locations)
437 :     {
438 :     if ($loc =~ /^(\S+)_(\d+)_(\d+)$/)
439 :     {
440 :     ($contig,$beg,$end) = ($1,$2,$3);
441 :     my $seq = $contigs{$contig};
442 :    
443 :     $ln = length($seq);
444 :    
445 :     if (! $ln) {
446 :     print STDERR "$genome/$contig: could not get length\n";
447 :     return "";
448 :     }
449 :    
450 :     if (&FIG::between(1,$beg,$ln) && &FIG::between(1,$end,$ln))
451 :     {
452 :     if ($beg < $end)
453 :     {
454 : paczian 1.15 push(@pieces, substr($seq, $beg - 1, ($end - $beg) + 1));
455 : olson 1.10 }
456 :     else
457 :     {
458 : paczian 1.15 push(@pieces, &FIG::reverse_comp(substr($seq, $end - 1, ($beg - $end) + 1)));
459 : olson 1.10 }
460 :     }
461 :     }
462 :     }
463 :     return lc(join("",@pieces));
464 :     }
465 :    
466 : overbeek 1.7 sub genome_szdna {
467 :     my ($self, $genome) = @_;
468 :    
469 :     my $fig = $self->{_fig};
470 :     my $newG = $self->{_genome};
471 :     my $newGdir = $self->{_orgdir};
472 :    
473 :     if ($genome ne $newG)
474 :     {
475 :     return $fig->genome_szdna($genome);
476 :     }
477 :     else
478 :     {
479 :     my $contig_lens = $self->contig_lengths($genome);
480 :     my $tot = 0;
481 :     while ( my($contig,$len) = each %$contig_lens)
482 :     {
483 :     $tot += $len;
484 :     }
485 :     return $tot;
486 :     }
487 :     }
488 :    
489 :     sub genome_version {
490 :     my ($self, $genome) = @_;
491 :    
492 :     my $fig = $self->{_fig};
493 :     my $newG = $self->{_genome};
494 :     my $newGdir = $self->{_orgdir};
495 :    
496 :     if ($genome ne $newG)
497 :     {
498 :     return $fig->genome_version($genome);
499 :     }
500 :     else
501 :     {
502 :     return "$genome.0";
503 :     }
504 :     }
505 :    
506 :     sub genome_pegs {
507 :     my ($self, $genome) = @_;
508 :    
509 :     my $fig = $self->{_fig};
510 :     my $newG = $self->{_genome};
511 :     my $newGdir = $self->{_orgdir};
512 :    
513 :     if ($genome ne $newG)
514 :     {
515 :     return $fig->genome_pegs($genome);
516 :     }
517 :     else
518 :     {
519 :     my @tmp = $self->all_features($genome,"peg");
520 :     my $n = @tmp;
521 :     return $n;
522 :     }
523 :     }
524 :    
525 :     sub genome_rnas {
526 :     my ($self, $genome) = @_;
527 :    
528 :     my $fig = $self->{_fig};
529 :     my $newG = $self->{_genome};
530 :     my $newGdir = $self->{_orgdir};
531 :    
532 :     if ($genome ne $newG)
533 :     {
534 :     return $fig->genome_rnas($genome);
535 :     }
536 :     else
537 :     {
538 :     my @tmp = $self->all_features($genome,"rna");
539 :     my $n = @tmp;
540 :     return $n;
541 :     }
542 :     }
543 :    
544 :     sub genome_domain {
545 :     my ($self, $genome) = @_;
546 :    
547 :     my $fig = $self->{_fig};
548 :     my $newG = $self->{_genome};
549 :     my $newGdir = $self->{_orgdir};
550 :    
551 :     if ($genome ne $newG)
552 :     {
553 :     return $fig->genome_domain($genome);
554 :     }
555 :     else
556 :     {
557 :     my $tax = $self->taxonomy_of($genome);
558 :     return ($tax =~ /^([^ \t;]+)/) ? $1 : "unknown";
559 :     }
560 :     }
561 : overbeek 1.2
562 :     sub genes_in_region {
563 :     my($self,$genome,$contig,$beg,$end) = @_;
564 :    
565 :     my $fig = $self->{_fig};
566 :     my $newG = $self->{_genome};
567 :     my $newGdir = $self->{_orgdir};
568 :    
569 :     if ($genome ne $newG)
570 :     {
571 :     return $fig->genes_in_region($genome,$contig,$beg,$end);
572 :     }
573 :     else
574 :     {
575 :     &load_tbl($self);
576 :     my $tblH = $self->{_tbl};
577 :     my $maxV = 0;
578 :     my $minV = 1000000000;
579 :     my $genes = [];
580 :     while ( my($fid,$tuple) = each %$tblH)
581 :     {
582 :     if (($tuple->[0]->[0] =~ /^(\S+)_(\d+)_\d+$/) && ($1 eq $contig))
583 :     {
584 :     my $beg1 = $2;
585 :     my $last = @{$tuple->[0]} - 1;
586 :     if (($tuple->[0]->[$last] =~ /^(\S+)_\d+_(\d+)$/) && ($1 eq $contig))
587 :     {
588 :     my $end1 = $2;
589 :     if (&overlaps($beg,$end,$beg1,$end1))
590 :     {
591 :     $minV = &FIG::min($minV,$beg1,$end1);
592 :     $maxV = &FIG::max($maxV,$beg1,$end1);
593 :     push(@$genes,$fid);
594 :     }
595 :     }
596 :     }
597 :     }
598 :     return ($genes,$minV,$maxV);
599 :     }
600 :     }
601 :    
602 :     sub overlaps {
603 :     my($b1,$e1,$b2,$e2) = @_;
604 :    
605 :     if ($b1 > $e1) { ($b1,$e1) = ($e1,$b1) }
606 :     if ($b2 > $e2) { ($b2,$e2) = ($e2,$b2) }
607 :     return &FIG::between($b1,$b2,$e1) || &FIG::between($b2,$b1,$e2);
608 :     }
609 :    
610 : overbeek 1.7 sub all_contigs {
611 :     my($self,$genome) = @_;
612 :    
613 :     my $fig = $self->{_fig};
614 :     my $newG = $self->{_genome};
615 :     my $newGdir = $self->{_orgdir};
616 :    
617 :     if ($genome ne $newG)
618 :     {
619 :     return $fig->all_contigs($genome);
620 :     }
621 :     else
622 :     {
623 :     &load_tbl($self);
624 :     my %contigs;
625 :     my $tblH = $self->{_tbl};
626 :     while ( my($fid,$tuple) = each %$tblH)
627 :     {
628 :     if ($tuple->[0]->[0] =~ /^(\S+)_\d+_\d+$/)
629 :     {
630 :     $contigs{$1} = 1;
631 :     }
632 :     }
633 :     return keys(%contigs);
634 :     }
635 :     }
636 :    
637 :     sub all_features {
638 :     my($self,$genome,$type) = @_;
639 :    
640 :     my $fig = $self->{_fig};
641 :     my $newG = $self->{_genome};
642 :     my $newGdir = $self->{_orgdir};
643 :    
644 :     if ($genome ne $newG)
645 :     {
646 :     return $fig->all_features($genome,$type);
647 :     }
648 :     else
649 :     {
650 :     &load_tbl($self);
651 :     my $tblH = $self->{_tbl};
652 :     return sort { &FIG::by_fig_id($a,$b) }
653 :     grep { ($_ =~ /^fig\|\d+\.\d+\.([^\.]+)/) && ($1 eq $type) }
654 :     keys(%$tblH);
655 :     }
656 :     }
657 :    
658 :     sub all_features_detailed_fast {
659 : olson 1.11 my($self,$genome, $regmin, $regmax, $contig) = @_;
660 : overbeek 1.7
661 :     my $fig = $self->{_fig};
662 :     my $newG = $self->{_genome};
663 :     my $newGdir = $self->{_orgdir};
664 :    
665 :     if ($genome ne $newG)
666 :     {
667 : olson 1.11 return $fig->all_features_detailed_fast($genome, $regmin, $regmax, $contig);
668 : overbeek 1.7 }
669 :     else
670 :     {
671 :     &load_tbl($self);
672 :     my $tblH = $self->{_tbl};
673 :     my $feat_details = [];
674 :     while ( my($fid,$tuple) = each %$tblH)
675 :     {
676 :     if ($fid =~ /^fig\|\d+\.\d+\.(\S+)\.\d+/)
677 :     {
678 :     my $type = $1;
679 : olson 1.11 if ($tuple->[0]->[0] =~ /^(\S+)_(\d+)_(\d+)$/)
680 : overbeek 1.7 {
681 : olson 1.11 my($ctg, $beg, $end) = ($1, $2, $3);
682 :     next if (defined($contig) and $contig ne $ctg);
683 :    
684 : overbeek 1.7 my($min,$max);
685 : olson 1.11 if ($beg < $end)
686 : overbeek 1.7 {
687 : olson 1.11 $min = $beg;
688 :     $max = $end;
689 : overbeek 1.7 }
690 :     else
691 :     {
692 : olson 1.11 $min = $end;
693 :     $max = $beg;
694 :     }
695 :    
696 :     if (not defined($regmin) or not defined($regmax) or
697 :     ($min < $regmin and $regmin < $max) or ($min < $regmax and $regmax < $max) or ($min > $regmin and $max < $regmax))
698 :     {
699 :     push(@$feat_details,[$fid,$tuple->[0]->[0],join(",",@{$tuple->[1]}),$type,$min,$max,$self->function_of($fid),'master','']);
700 : overbeek 1.7 }
701 :     }
702 :     }
703 :     }
704 :     return $feat_details;
705 :     }
706 :     }
707 :    
708 :     sub compute_clusters {
709 :     # Get the parameters.
710 :     my ($self, $pegList, $subsystem, $distance) = @_;
711 :     if (! defined $distance) {
712 :     $distance = 5000;
713 :     }
714 :    
715 :     my($peg,%by_contig);
716 :     foreach $peg (@$pegList)
717 :     {
718 :     my $loc;
719 :     if ($loc = $self->feature_location($peg))
720 :     {
721 :     my ($contig,$beg,$end) = &FIG::boundaries_of($loc);
722 :     my $genome = &FIG::genome_of($peg);
723 :     push(@{$by_contig{"$genome\t$contig"}},[($beg+$end)/2,$peg]);
724 :     }
725 :     }
726 :    
727 :     my @clusters = ();
728 :     foreach my $tuple (keys(%by_contig))
729 :     {
730 :     my $x = $by_contig{$tuple};
731 :     my @pegs = sort { $a->[0] <=> $b->[0] } @$x;
732 :     while ($x = shift @pegs)
733 :     {
734 :     my $clust = [$x->[1]];
735 :     while ((@pegs > 0) && (abs($pegs[0]->[0] - $x->[0]) <= $distance))
736 :     {
737 :     $x = shift @pegs;
738 :     push(@$clust,$x->[1]);
739 :     }
740 :    
741 :     if (@$clust > 1)
742 :     {
743 :     push(@clusters,$clust);
744 :     }
745 :     }
746 :     }
747 :     return sort { @$b <=> @$a } @clusters;
748 :     }
749 :    
750 :    
751 : overbeek 1.1 sub feature_location {
752 :     my($self,$fid) = @_;
753 :    
754 :     my $fig = $self->{_fig};
755 :     my $newG = $self->{_genome};
756 :     my $newGdir = $self->{_orgdir};
757 :    
758 :     if (($fid =~ /^fig\|(\d+\.\d+)/) && ($1 eq $newG))
759 :     {
760 :     &load_tbl($self);
761 :     if (my $x = $self->{_tbl}->{$fid})
762 :     {
763 :     return join(",",@{$x->[0]});
764 :     }
765 :     else
766 :     {
767 :     return undef;
768 :     }
769 :     }
770 :     else
771 :     {
772 :     return scalar $fig->feature_location($fid);
773 :     }
774 :     }
775 :    
776 :     sub function_of {
777 :     my($self,$fid) = @_;
778 :    
779 :     my $fig = $self->{_fig};
780 :     my $newG = $self->{_genome};
781 :     my $newGdir = $self->{_orgdir};
782 :    
783 :     if (($fid =~ /^fig\|(\d+\.\d+)/) && ($1 eq $newG))
784 :     {
785 :     &load_functions($self);
786 :     return $self->{_functions}->{$fid};
787 :     }
788 :     else
789 :     {
790 :     return scalar $fig->function_of($fid);
791 :     }
792 :     }
793 :    
794 :     sub feature_aliases {
795 :     my($self,$fid) = @_;
796 :    
797 :     my $fig = $self->{_fig};
798 :     my $newG = $self->{_genome};
799 :     my $newGdir = $self->{_orgdir};
800 :    
801 :     my @aliases;
802 :     if (($fid =~ /^fig\|(\d+\.\d+)/) && ($1 eq $newG))
803 :     {
804 :     &load_tbl($self);
805 :     if (my $x = $self->{_tbl}->{$fid})
806 :     {
807 :     @aliases = @{$x->[1]};
808 :     }
809 :     else
810 :     {
811 :     @aliases = ();
812 :     }
813 :     }
814 :     else
815 :     {
816 :     @aliases = $fig->feature_aliases($fid);
817 :     }
818 :     return wantarray() ? @aliases : join(",",@aliases);
819 :     }
820 :    
821 :     sub feature_annotations {
822 :     my($self,$fid,$rawtime) = @_;
823 :    
824 :     my $fig = $self->{_fig};
825 :     my $newG = $self->{_genome};
826 :     my $newGdir = $self->{_orgdir};
827 :    
828 :     my @annotations;
829 :     if (($fid =~ /^fig\|(\d+\.\d+)/) && ($1 eq $newG))
830 :     {
831 :     &load_ann($self);
832 :     if (my $x = $self->{_ann}->{$fid})
833 :     {
834 :     @annotations = @{$x};
835 :     }
836 :     else
837 :     {
838 :     @annotations = ();
839 :     }
840 :    
841 :     if ($rawtime)
842 :     {
843 :     return @annotations;
844 :     }
845 :     else
846 :     {
847 :     return map { $_->[1] = localtime($_->[1]); $_ } @annotations;
848 :     }
849 :     }
850 :     else
851 :     {
852 :     return $fig->feature_annotations($fid);
853 :     }
854 :     }
855 :    
856 :     sub get_translation {
857 :     my($self,$peg) = @_;
858 :    
859 :     my $fig = $self->{_fig};
860 :     my $newG = $self->{_genome};
861 :     my $newGdir = $self->{_orgdir};
862 :    
863 :     if (($peg =~ /^fig\|(\d+\.\d+)/) && ($1 eq $newG))
864 :     {
865 :     &load_pseq($self);
866 :     return $self->{_pseq}->{$peg};
867 :     }
868 :     else
869 :     {
870 :     return $fig->get_translation($peg);
871 :     }
872 :     }
873 :    
874 : olson 1.4 sub translation_length
875 :     {
876 :     my($self, $peg) = @_;
877 :    
878 :     my $fig = $self->{_fig};
879 :     my $newG = $self->{_genome};
880 :     my $newGdir = $self->{_orgdir};
881 :    
882 :     if (($peg =~ /^fig\|(\d+\.\d+)/) && ($1 eq $newG))
883 :     {
884 :     &load_pseq($self);
885 :     return length($self->{_pseq}->{$peg});
886 :     }
887 :     else
888 :     {
889 :     return $fig->translation_length($peg);
890 :     }
891 :     }
892 :    
893 :     sub translatable
894 :     {
895 :     my($self, $peg) = @_;
896 :    
897 :     my $fig = $self->{_fig};
898 :     my $newG = $self->{_genome};
899 :     my $newGdir = $self->{_orgdir};
900 :    
901 :     if (($peg =~ /^fig\|(\d+\.\d+)/) && ($1 eq $newG))
902 :     {
903 :     &load_pseq($self);
904 :     return length($self->{_pseq}->{$peg}) > 0 ? 1 : 0;;
905 :     }
906 :     else
907 :     {
908 :     return $fig->translatable($peg);
909 :     }
910 :     }
911 :    
912 :     sub is_real_feature
913 :     {
914 :     my($self, $fid) = @_;
915 :    
916 :     if ($self->is_virtual_feature($fid))
917 :     {
918 :     return 1;
919 :     }
920 :     else
921 :     {
922 :     return $self->{_fig}->is_real_feature($fid);
923 :     }
924 :    
925 :     }
926 :    
927 : olson 1.10 sub pegs_of
928 :     {
929 :     my($self, $genome) = @_;
930 :    
931 :     my $fig = $self->{_fig};
932 :     my $newG = $self->{_genome};
933 :     my $newGdir = $self->{_orgdir};
934 :    
935 :     if ($genome ne $newG)
936 :     {
937 :     return $fig->pegs_of($genome);
938 :     }
939 :    
940 :     $self->load_tbl();
941 :     return grep { /\.peg\./ } keys %{$self->{_tbl}};
942 :     }
943 :    
944 :    
945 :     sub rnas_of
946 :     {
947 :     my($self, $genome) = @_;
948 :    
949 :     my $fig = $self->{_fig};
950 :     my $newG = $self->{_genome};
951 :     my $newGdir = $self->{_orgdir};
952 :    
953 :     if ($genome ne $newG)
954 :     {
955 :     return $fig->pegs_of($genome);
956 :     }
957 :    
958 :     $self->load_tbl();
959 :     return grep { /\.rna\./ } keys %{$self->{_tbl}};
960 :     }
961 :    
962 : olson 1.4 sub is_virtual_feature
963 :     {
964 :     my($self, $peg) = @_;
965 :    
966 :     my $fig = $self->{_fig};
967 :     my $newG = $self->{_genome};
968 :     my $newGdir = $self->{_orgdir};
969 :    
970 :     if (($peg =~ /^fig\|(\d+\.\d+)/) && ($1 eq $newG))
971 :     {
972 :     &load_pseq($self);
973 :     return 1;
974 :     }
975 :     else
976 :     {
977 :     return 0;
978 :     }
979 :     }
980 :    
981 : overbeek 1.1 sub load_pseq {
982 :     my($self) = @_;
983 :    
984 :     if ($self->{_pseq}) { return };
985 :    
986 :     my $newGdir = $self->{_orgdir};
987 :     my $pseq = {};
988 :     if (open(FASTA,"<$newGdir/Features/peg/fasta"))
989 :     {
990 : olson 1.10 local $/ = "\n>";
991 : overbeek 1.1 my $x;
992 :     while (defined($x = <FASTA>))
993 :     {
994 : olson 1.10 if ($x =~ />?(\S+)[^\n]*\n(.*)/s)
995 : overbeek 1.1 {
996 :     my $peg = $1;
997 :     my $seq = $2;
998 :     $seq =~ s/\s//gs;
999 :     $pseq->{$peg} = $seq;
1000 :     }
1001 :     }
1002 :     close(FASTA);
1003 :     }
1004 :     $self->{_pseq} = $pseq;
1005 :     }
1006 :    
1007 : olson 1.9 sub load_ss_data
1008 :     {
1009 :     my($self) = @_;
1010 :    
1011 :     return if defined($self->{_ss_bindings});
1012 :    
1013 :     my $fig = $self->{_fig};
1014 :     my $newG = $self->{_genome};
1015 :     my $newGdir = $self->{_orgdir};
1016 :    
1017 :     open(S, "<$newGdir/Subsystems/bindings") or die "Cannot open $newGdir/Subsystems/bindings: $!";
1018 :    
1019 :     my $peg_index;
1020 :     my $bindings;
1021 : olson 1.14 my $genome_index;
1022 :     my %seen;
1023 : olson 1.9 while (<S>)
1024 :     {
1025 :     chomp;
1026 :     my($sname, $role, $peg) = split(/\t/);
1027 : olson 1.14
1028 :     my $genome = &FIG::genome_of($peg);
1029 : olson 1.9
1030 :     push(@{$bindings->{$sname}->{$role}}, $peg);
1031 :     push(@{$peg_index->{$peg}}, [$sname, $role]);
1032 : olson 1.14
1033 :     if (!$seen{$genome, $sname, $role})
1034 :     {
1035 :     push(@{$genome_index->{$genome}, [$sname, $role]});
1036 :     }
1037 :    
1038 : olson 1.9 }
1039 :     close(S);
1040 :    
1041 :     open(S, "<$newGdir/Subsystems/subsystems") or die "Cannot open $newGdir/Subsystems/subsystems: $!";
1042 :     my $variant;
1043 :     while (<S>)
1044 :     {
1045 :     chomp;
1046 :     my($sname, $v) = split(/\t/);
1047 :     $variant->{$sname} = $v;
1048 :     }
1049 :     close(S);
1050 :    
1051 :     $self->{_ss_bindings} = $bindings;
1052 :     $self->{_ss_variants} = $variant;
1053 :     $self->{_ss_peg_index} = $peg_index;
1054 : olson 1.14 $self->{_ss_genome_index} = $genome_index;
1055 : olson 1.9 }
1056 :    
1057 : overbeek 1.1 sub load_tbl {
1058 :     my($self) = @_;
1059 :    
1060 :     if ($self->{_tbl}) { return };
1061 :    
1062 :     my $newGdir = $self->{_orgdir};
1063 :     my $tbl = {};
1064 :     foreach my $x (`cat $newGdir/Features/*/tbl`)
1065 :     {
1066 :     if ($x =~ /^(\S+)\t(\S+)(\t(\S.*\S))?/)
1067 :     {
1068 :     my $fid = $1;
1069 :     my $loc = [split(/,/,$2)];
1070 : paczian 1.12 my $aliases = $4 ? [split(/\t/,$4)] : [];
1071 : overbeek 1.1 $tbl->{$fid} = [$loc,$aliases];
1072 :     }
1073 :     }
1074 :     $self->{_tbl} = $tbl;
1075 :     }
1076 :    
1077 :     sub load_functions {
1078 :     my($self) = @_;
1079 :    
1080 :     if ($self->{_functions}) { return };
1081 :    
1082 :     my $newG = $self->{_genome};
1083 :     my $newGdir = $self->{_orgdir};
1084 :    
1085 :     my $functions = {};
1086 :     foreach my $x (`cat $newGdir/*functions`)
1087 :     {
1088 :     if (($x =~ /^(fig\|(\d+\.\d+)\.\S+)\t(\S[^\t]*\S)/) && ($2 eq $newG))
1089 :     {
1090 :     $functions->{$1} = $3;
1091 :     }
1092 :     }
1093 :     $self->{_functions} = $functions;
1094 :     }
1095 :    
1096 : olson 1.16
1097 :     sub bbhs
1098 :     {
1099 :     my($self,$peg,$cutoff) = @_;
1100 :    
1101 :     my $fig = $self->{_fig};
1102 :     my $newG = $self->{_genome};
1103 :     my $newGdir = $self->{_orgdir};
1104 :    
1105 :     $self->load_bbh_index();
1106 :    
1107 :     $cutoff = 1.0e-10 unless defined($cutoff);
1108 :    
1109 :     my $tie = $self->{_bbh_tie};
1110 :    
1111 :     my @bbhs = $tie->get_dup($peg);
1112 :    
1113 :     # @bbhs now a list of comma-separated tuples (id2, score, bitscore)
1114 :     # print Dumper(\@bbhs);
1115 :    
1116 :     @bbhs = grep { $_->[1] < $cutoff } map { [split(/,/)] } @bbhs;
1117 :    
1118 :     if (not ($peg =~ /^fig\|(\d+\.\d+)/ && ($1 eq $newG)))
1119 :     {
1120 :     #
1121 :     # If this isn't one of our pegs, we retrieve the bbhs from the underlying
1122 :     # SEED and merge in the bbhs that we have here.
1123 :     #
1124 :    
1125 :     my @fbbhs = $fig->bbhs($peg, $cutoff);
1126 :     # print Dumper(\@fbbhs);
1127 :     push(@bbhs, @fbbhs);
1128 :     }
1129 :     return sort { $a->[1] <=> $b->[1] } @bbhs;
1130 :     }
1131 :    
1132 : olson 1.13 sub sims
1133 :     {
1134 :     my($self,$peg,$max,$maxP,$select) = @_;
1135 :    
1136 :     my $fig = $self->{_fig};
1137 :     my $newG = $self->{_genome};
1138 :     my $newGdir = $self->{_orgdir};
1139 :     my $max = $max ? $max : 10000;
1140 :     my $maxP = $maxP ? $maxP : 1.0e-5;
1141 :     $select = $select ? $select : "all";
1142 :    
1143 :     if (not ($peg =~ /^fig\|(\d+\.\d+)/) && ($1 eq $newG))
1144 :     {
1145 :     return $fig->sims($peg,$max,$maxP,$select);
1146 :     }
1147 :    
1148 :     $self->load_sims_index();
1149 :     my $sims_file = "$newGdir/expanded_similarities";
1150 :    
1151 :     my $info = $self->{_sims_index}->{$peg};
1152 :    
1153 :     if ($info !~ /^(\d+),(\d+)$/)
1154 :     {
1155 :     return;
1156 :     }
1157 :     my($seek, $len) = ($1, $2);
1158 :    
1159 :     my $sims_txt = &FIG::read_block($self->{_sims_fh}, $seek, $len);
1160 :    
1161 :     my @sims = map { bless $_, 'Sim' } grep { $_->[10] <= $maxP } map { [split(/\t/)] } @$sims_txt;
1162 :    
1163 :     if ($select =~ /^figx?$/)
1164 :     {
1165 :     @sims = grep { $_->[1] =~ /^fig/ } @sims;
1166 :     }
1167 :    
1168 :     if (@sims > $max)
1169 :     {
1170 :     $#sims = $max-1;
1171 :     }
1172 :    
1173 :     return @sims;
1174 :     }
1175 :    
1176 :    
1177 :     sub sims_old {
1178 : overbeek 1.1 my($self,$peg,$max,$maxP,$select) = @_;
1179 :    
1180 :     my $fig = $self->{_fig};
1181 :     my $newG = $self->{_genome};
1182 :     my $newGdir = $self->{_orgdir};
1183 :     my $max = $max ? $max : 10000;
1184 :     my $maxP = $maxP ? $maxP : 1.0e-5;
1185 :     $select = $select ? $select : "all";
1186 :    
1187 :     if (($peg =~ /^fig\|(\d+\.\d+)/) && ($1 eq $newG))
1188 :     {
1189 :     &load_sims($self);
1190 :     my @sims = ();
1191 :     my $raw_sims = $self->{_sims}->{$peg};
1192 :     if ($raw_sims)
1193 :     {
1194 :     foreach my $sim ( grep { $_->[10] <= $maxP } @$raw_sims )
1195 :     {
1196 :     my $id2 = $sim->id2;
1197 :     my $id1 = $sim->id1;
1198 :     my @relevant = ();
1199 :    
1200 :     my @maps_to = $fig->mapped_prot_ids( $id2 );
1201 :     my $ref_len = $maps_to[0]->[1];
1202 :    
1203 :     @maps_to = grep { $_->[0] !~ /^xxx\d+/ } @maps_to;
1204 :    
1205 :     if ( $select =~ /^figx?$/ ) # Only fig
1206 :     {
1207 :     @relevant = grep { $_->[0] =~ /^fig/ } @maps_to;
1208 :     }
1209 :     else # All
1210 :     {
1211 :     @relevant = @maps_to;
1212 :     }
1213 :    
1214 :     my $seen = {};
1215 :     foreach my $x ( @relevant )
1216 :     {
1217 :     my ( $x_id, $x_ln ) = @$x;
1218 :    
1219 :     next if $seen->{$x_id};
1220 :     $seen->{$x_id} = 1;
1221 :    
1222 :     my $delta2 = $ref_len - $x_ln; # Coordinate shift
1223 :     my $sim1 = [ @$sim ]; # Make a copy
1224 :     $sim1->[1] = $x_id;
1225 :     $sim1->[8] -= $delta2;
1226 :     $sim1->[9] -= $delta2;
1227 :     bless( $sim1, "Sim" );
1228 :     push( @sims, $sim1 );
1229 :     }
1230 :     }
1231 :     }
1232 :    
1233 :     if (@sims > $max) { $#sims = $max-1; }
1234 :     return @sims;
1235 :     }
1236 :     else
1237 :     {
1238 :     return $fig->sims($peg,$max,$maxP,$select);
1239 :     }
1240 :     }
1241 :    
1242 : olson 1.16 sub load_bbh_index
1243 :     {
1244 :     my($self) = @_;
1245 :    
1246 :     return if defined($self->{_bbh_index});
1247 :    
1248 :     my $fig = $self->{_fig};
1249 :     my $newG = $self->{_genome};
1250 :     my $newGdir = $self->{_orgdir};
1251 :    
1252 :     my $bbh_file = "$newGdir/bbhs";
1253 :     my $bbh_index_file = "$bbh_file.index";
1254 :     my $bbh_index = {};
1255 :    
1256 :     my $tied = tie %$bbh_index, 'DB_File', $bbh_index_file, O_RDONLY, 0666, $DB_BTREE;
1257 :    
1258 :     #
1259 :     # Set these even if failed so we don't keep trying to open and failing.
1260 :     #
1261 :     $self->{_bbh_index} = $bbh_index;
1262 :     $self->{_bbh_tie} = $tied;
1263 :    
1264 :     if (not $tied)
1265 :     {
1266 :     warn "Cannot tie bbh index $bbh_index_file: $!\n";
1267 :     }
1268 :     }
1269 :    
1270 : olson 1.13 sub load_sims_index
1271 :     {
1272 :     my($self) = @_;
1273 :    
1274 :     return if defined($self->{_sims_index});
1275 :    
1276 :     my $fig = $self->{_fig};
1277 :     my $newG = $self->{_genome};
1278 :     my $newGdir = $self->{_orgdir};
1279 :    
1280 :     my $sims_file = "$newGdir/expanded_similarities";
1281 :     my $sims_index_file = "$sims_file.index";
1282 :     my $sims_index = {};
1283 :    
1284 :     my $tied = tie %$sims_index, 'DB_File', $sims_index_file, O_RDONLY, 0666, $DB_BTREE;
1285 :    
1286 :     #
1287 :     # Set these even if failed so we don't keep trying to open and failing.
1288 :     #
1289 :     $self->{_sims_index} = $sims_index;
1290 :     $self->{_sims_tie} = $tied;
1291 :    
1292 :     if (not $tied)
1293 :     {
1294 :     warn "Cannot tie sims index $sims_index_file: $!\n";
1295 :     }
1296 :    
1297 :     #
1298 :     # open the sims file as well.
1299 :     #
1300 :    
1301 :     $self->{_sims_fh} = new FileHandle("<$sims_file");
1302 :    
1303 :     if (!$self->{_sims_fh})
1304 :     {
1305 :     warn "Cannot open sims file $sims_file: $!\n";
1306 :     }
1307 :     }
1308 :    
1309 : overbeek 1.1 sub load_sims {
1310 :     my($self) = @_;
1311 :    
1312 :     if ($self->{_sims}) { return };
1313 :    
1314 :     my $newGdir = $self->{_orgdir};
1315 :    
1316 :     my $sims = {};
1317 :     foreach my $x (`cat $newGdir/similarities`)
1318 :     {
1319 :     chop; $x;
1320 :     if ($x =~ /^(\S+)/)
1321 :     {
1322 :     push(@{$sims->{$1}},bless([split(/\t/,$x)],'Sim'));
1323 :     }
1324 :     }
1325 :     $self->{_sims} = $sims;
1326 :     }
1327 :    
1328 :     sub get_attributes {
1329 :     my($self,$peg) = @_;
1330 :    
1331 :     my $fig = $self->{_fig};
1332 :     my $newG = $self->{_genome};
1333 :     my $newGdir = $self->{_orgdir};
1334 :    
1335 :     if (($peg =~ /^fig\|(\d+\.\d+)\.peg\.\d+$/) && ($1 eq $newG))
1336 :     {
1337 :     &load_attr($self);
1338 :     if (my $x = $self->{_attr}->{$peg})
1339 :     {
1340 :     return @$x;
1341 :     }
1342 :     else
1343 :     {
1344 :     return ();
1345 :     }
1346 :     }
1347 :     else
1348 :     {
1349 :     return $fig->get_attributes($peg);
1350 :     }
1351 :     }
1352 :    
1353 :     sub load_attr {
1354 :     my($self) = @_;
1355 :    
1356 :     if ($self->{_attr}) { return };
1357 :    
1358 :     my $newGdir = $self->{_orgdir};
1359 :     my $attr = {};
1360 : overbeek 1.3 foreach my $x (`cat $newGdir/evidence.codes`)
1361 : overbeek 1.1 {
1362 :     if ($x =~ /^(\S+)\t(\S+)/)
1363 :     {
1364 :     push(@{$attr->{$1}},[$1,"evidence_code",$2,""]);
1365 :     }
1366 :     }
1367 :     $self->{_attr} = $attr;
1368 :     }
1369 :    
1370 :     sub load_ann {
1371 :     my($self) = @_;
1372 :    
1373 :     if ($self->{_ann}) { return };
1374 :    
1375 :     my $newGdir = $self->{_orgdir};
1376 :     my $ann = {};
1377 :     if (open(ANN,"<$newGdir/annotations"))
1378 :     {
1379 :     $/ = "\n//\n";
1380 :     while (defined(my $x = <ANN>))
1381 :     {
1382 :     chomp $x;
1383 :     if ($x =~ /^(\S+)\n([^\n]+)\n([^\n]+)\n(.*)/s)
1384 :     {
1385 :     push(@{$ann->{$1}},[$1,$2,$3,"$4\n"]);
1386 :     }
1387 :     }
1388 :     $/ = "\n";
1389 :     close(ANN);
1390 :     }
1391 :     $self->{_ann} = $ann;
1392 :     }
1393 : overbeek 1.6
1394 :     sub taxonomy_of {
1395 :     my($self,$genome) = @_;
1396 :     my $fig = $self->{_fig};
1397 :     my $newG = $self->{_genome};
1398 :     my $newGdir = $self->{_orgdir};
1399 :    
1400 :     if ($newG ne $genome)
1401 :     {
1402 :     return $fig->taxonomy_of($genome);
1403 :     }
1404 :    
1405 :     my $tax;
1406 :     if (open(TAX,"<$newGdir/TAXONOMY") && ($tax = <TAX>))
1407 :     {
1408 :     chop $tax;
1409 :     return $tax;
1410 :     }
1411 :     else
1412 :     {
1413 :     return "unknown";
1414 :     }
1415 :     }
1416 :    
1417 :     sub build_tree_of_complete {
1418 :     my($self,$min_for_label) = @_;
1419 :     return $self->build_tree_of_all($min_for_label, "complete");
1420 :     }
1421 :    
1422 :     sub build_tree_of_all {
1423 :     my($self, $min_for_label, $complete)=@_;
1424 :     my(@last,@tax,$i,$prefix,$lev,$genome,$tax);
1425 :    
1426 :     $min_for_label = $min_for_label ? $min_for_label : 10;
1427 :     open(TMP,">/tmp/tree$$") || die "could not open /tmp/tree$$";
1428 :     print TMP "1. root\n";
1429 :    
1430 :     @last = ();
1431 :    
1432 :    
1433 :     foreach $genome (grep { ! $self->is_environmental($_) } $self->sort_genomes_by_taxonomy($self->genomes($complete)))
1434 :     {
1435 :     $tax = $self->taxonomy_of($genome);
1436 :     @tax = split(/\s*;\s*/,$tax);
1437 :     push(@tax,$genome);
1438 :     for ($i=0; ((@last > $i) && (@tax > $i) && ($last[$i] eq $tax[$i])); $i++) {}
1439 :     while ($i < @tax)
1440 :     {
1441 :     $lev = $i+2;
1442 :     $prefix = " " x (4 * ($lev-1));
1443 :     print TMP "$prefix$lev\. $tax[$i]\n";
1444 :     $i++;
1445 :     }
1446 :     @last = @tax;
1447 :     }
1448 :     close(TMP);
1449 :     my $tree = &tree_utilities::build_tree_from_outline("/tmp/tree$$");
1450 :     $tree->[0] = 'All';
1451 :     &FIG::limit_labels($tree,$min_for_label);
1452 :     unlink("/tmp/tree$$");
1453 :     return ($tree,&tree_utilities::tips_of_tree($tree));
1454 :     }
1455 :    
1456 :     sub sort_genomes_by_taxonomy {
1457 :     my($self,@genomes) = @_;
1458 :    
1459 :     return map { $_->[0] }
1460 :     sort { $a->[1] cmp $b->[1] }
1461 :     map { [$_,$self->taxonomy_of($_)] }
1462 :     @genomes;
1463 :     }
1464 : overbeek 1.1
1465 : overbeek 1.6 sub taxonomic_groups_of_complete {
1466 :     my($self,$min_for_labels) = @_;
1467 :    
1468 :     my($tree,undef) = $self->build_tree_of_complete($min_for_labels);
1469 :     return &FIG::taxonomic_groups($tree);
1470 :     }
1471 :    
1472 : overbeek 1.1 1;

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3