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

Annotation of /FigKernelPackages/FIGV.pm

Parent Directory Parent Directory | Revision Log Revision Log


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

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3