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

Annotation of /FigKernelPackages/FIGV.pm

Parent Directory Parent Directory | Revision Log Revision Log


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

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3