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

Annotation of /FigKernelPackages/FIGV.pm

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1.39 - (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 : olson 1.27 sub is_complete
66 :     {
67 :     return 1;
68 :     }
69 :    
70 : olson 1.4 #
71 :     # Redirect any method invocations that we don't implement out to the
72 :     # underlying FIG object.
73 :     #
74 :     sub AUTOLOAD
75 :     {
76 :     my($self, @args) = @_;
77 :    
78 : overbeek 1.6 if (ref($self) ne "FIGV") {
79 :     confess "BAD FIGV object passed to AUTOLOAD";
80 :     }
81 :    
82 : olson 1.27 no strict 'refs';
83 :    
84 : olson 1.4 my $meth = $AUTOLOAD;
85 :     $meth =~ s/.*:://;
86 : olson 1.27 my $fmeth = "FIG::$meth";
87 : olson 1.4
88 :     my $fig = $self->{_fig};
89 : olson 1.8 # my $args = Dumper(\@args);
90 : olson 1.4 if (wantarray)
91 :     {
92 :     my @res = $fig->$meth(@args);
93 : olson 1.27 # my @res = &$fmeth($self, @args);
94 : olson 1.5 # warn "FIGV invoke $meth($args) returns\n", Dumper(\@res);
95 : olson 1.4 return @res;
96 :     }
97 :     else
98 :     {
99 :     my $res = $fig->$meth(@args);
100 : olson 1.27 # my $res = &$fmeth($self, @args);
101 : olson 1.5 # warn "FIGV invoke $meth($args) returns\n", Dumper($res);
102 : olson 1.4 return $res;
103 :     }
104 :     }
105 :    
106 : olson 1.8 sub FIG
107 :     {
108 :     my($self) = @_;
109 :     return $self;
110 :     }
111 :    
112 : olson 1.27 sub sort_fids_by_taxonomy
113 :     {
114 :     my($self,@fids) = @_;
115 :    
116 :     return map { $_->[0] }
117 :     sort { $a->[1] cmp $b->[1] }
118 :     map { [$_,$self->taxonomy_of($self->genome_of($_))] }
119 :     @fids;
120 :     }
121 : olson 1.8
122 : olson 1.34 sub get_basic_statistics
123 :     {
124 :     my($self, $genome) = @_;
125 :    
126 :     my $fig = $self->{_fig};
127 :     my $newG = $self->{_genome};
128 :     my $newGdir = $self->{_orgdir};
129 :    
130 :     if ($genome ne $newG)
131 :     {
132 :     return $fig->get_basic_statistics($genome);
133 :     }
134 :    
135 :     #
136 :     # Check cache.
137 :     #
138 :    
139 :     my $cache = "$newGdir/cache.basic_statistics";
140 :     my $fh = new FileHandle($cache);
141 :     if ($fh)
142 :     {
143 :     my $stats = {};
144 :     while (<$fh>)
145 :     {
146 :     chomp;
147 :     my($k, $v) = split(/\t/);
148 :     $stats->{$k} = $v;
149 :     }
150 :     close($fh);
151 :     return $stats;
152 :     }
153 :    
154 :    
155 :     my $subsystem_data = $self->get_genome_subsystem_data($genome);
156 :    
157 :     my %sscount = map { $_->[0] => 1 } @$subsystem_data;
158 :     my $nss=scalar(keys(%sscount));
159 :    
160 :     my $statistics = {
161 :     num_subsystems => $nss,
162 :     num_contigs => scalar($self->all_contigs($genome)),
163 :     num_basepairs => $self->genome_szdna($genome),
164 :     genome_name => $self->genus_species($genome),
165 :     genome_domain => $self->genome_domain($genome),
166 :     genome_pegs => $self->genome_pegs($genome),
167 :     genome_rnas => $self->genome_rnas($genome),
168 :     genome_version => $self->genome_version($genome)
169 :     };
170 :    
171 :     $fh = new FileHandle(">$cache");
172 :     if ($fh)
173 :     {
174 :     while (my($k, $v) = each %$statistics)
175 :     {
176 :     print $fh join("\t", $k, $v), "\n";
177 :     }
178 :     close($fh);
179 :     }
180 :    
181 :     return $statistics;
182 :     }
183 :    
184 :    
185 :     sub get_peg_statistics {
186 :     my ($self, $genome) = @_;
187 :    
188 :     my $fig = $self->{_fig};
189 :     my $newG = $self->{_genome};
190 :     my $newGdir = $self->{_orgdir};
191 :    
192 :     if ($genome ne $newG)
193 :     {
194 :     return $fig->get_peg_statistics($genome);
195 :     }
196 :    
197 :     #
198 :     # Check cache.
199 :     #
200 :    
201 :     my $cache = "$newGdir/cache.peg_statistics";
202 :     my $fh = new FileHandle($cache);
203 :     if ($fh)
204 :     {
205 :     my $stats = {};
206 :     while (<$fh>)
207 :     {
208 :     chomp;
209 :     my($k, $v) = split(/\t/);
210 :     $stats->{$k} = $v;
211 :     }
212 :     close($fh);
213 :     return $stats;
214 :     }
215 :    
216 :    
217 :     my $subsystem_data = $self->get_genome_subsystem_data($genome);
218 :     my $assignment_data = $self->get_genome_assignment_data($genome);
219 :    
220 :     my $hypo_sub = 0;
221 :     my $hypo_nosub = 0;
222 :     my $nothypo_sub = 0;
223 :     my $nothypo_nosub = 0;
224 :     my %in = map { $_->[2] => 1 } @$subsystem_data;
225 :     my $in = keys(%in);
226 :    
227 :     my %sscount = map { $_->[0] => 1 } @$subsystem_data;
228 :    
229 :     foreach $_ (@$assignment_data)
230 :     {
231 :     my($peg,$func) = @$_;
232 :     my $is_hypo = &FIG::hypo($func);
233 :    
234 :     if ($is_hypo && $in{$peg}) { $hypo_sub++ }
235 :     elsif ($is_hypo && ! $in{$peg}) { $hypo_nosub++ }
236 :     elsif ((! $is_hypo) && (! $in{$peg})) { $nothypo_nosub++ }
237 :     elsif ((! $is_hypo) && $in{$peg}) { $nothypo_sub++ }
238 :     }
239 :     my $tot = $hypo_sub + $nothypo_sub + $hypo_nosub + $nothypo_nosub;
240 :    
241 :     my ($fracHS, $fracNHS, $fracHNS, $fracNHNS);
242 :    
243 :     if ($tot == 0) {
244 :     $fracHS = sprintf "%.2f", 0.0;
245 :     $fracNHS = sprintf "%.2f", 0.0;
246 :     $fracHNS = sprintf "%.2f", 0.0;
247 :     $fracNHNS = sprintf "%.2f", 0.0;
248 :     } else {
249 :     $fracHS = sprintf "%.2f", $hypo_sub / $tot * 100;
250 :     $fracNHS = sprintf "%.2f", $nothypo_sub / $tot * 100;
251 :     $fracHNS = sprintf "%.2f", $hypo_nosub / $tot * 100;
252 :     $fracNHNS = sprintf "%.2f", $nothypo_nosub / $tot * 100;
253 :     }
254 :    
255 :     my $statistics = {
256 :     hypothetical_in_subsystem => $hypo_sub,
257 :     hypothetical_not_in_subsystem => $hypo_nosub,
258 :     non_hypothetical_in_subsystem => $nothypo_sub,
259 :     non_hypothetical_not_in_subsystem => $nothypo_nosub,
260 :     hypothetical_in_subsystem_percent => $fracHS,
261 :     hypothetical_not_in_subsystem_percent => $fracHNS,
262 :     non_hypothetical_in_subsystem_percent => $fracNHS,
263 :     non_hypothetical_not_in_subsystem_percent => $fracNHNS
264 :     };
265 :    
266 :     $fh = new FileHandle(">$cache");
267 :     if ($fh)
268 :     {
269 :     while (my($k, $v) = each %$statistics)
270 :     {
271 :     print $fh join("\t", $k, $v), "\n";
272 :     }
273 :     close($fh);
274 :     }
275 :    
276 :     return $statistics;
277 :     }
278 :    
279 : olson 1.5 #
280 :     # To retrieve a subsystem in FIGV, we create the subsystem as normal via $fig->get_subsystem,
281 :     # then insert the row for the virtual org dir we are processing.
282 :     #
283 :    
284 :     sub get_subsystem
285 :     {
286 :     my($self,$ssa) = @_;
287 :    
288 :     my $fig = $self->{_fig};
289 :     my $newG = $self->{_genome};
290 :     my $newGdir = $self->{_orgdir};
291 :    
292 :     my $ss = $fig->get_subsystem($ssa);
293 :     return undef unless $ss;
294 :    
295 : olson 1.9 $self->load_ss_data();
296 : olson 1.8
297 : olson 1.9 my $bindings = $self->{_ss_bindings}->{$ssa};
298 :     my $variant = $self->{_ss_variants}->{$ssa};
299 :    
300 :     # warn "Adding virtual genome " . Dumper(\%bindings);
301 :     $ss->add_virtual_genome($self->genus_species(), $newG, $variant, $bindings);
302 :    
303 :     return $ss;
304 :     }
305 :    
306 : olson 1.14 sub active_subsystems
307 :     {
308 :     my($self, $genome, $all) = @_;
309 :    
310 :     my $fig = $self->{_fig};
311 :     my $newG = $self->{_genome};
312 :     my $newGdir = $self->{_orgdir};
313 :    
314 :     if ($genome ne $newG)
315 :     {
316 :     return $fig->active_subsystems($genome, $all);
317 :     }
318 :    
319 :     $self->load_ss_data();
320 :    
321 :     my $slist = {};
322 : olson 1.25
323 :     if ($self->{_ss_variants})
324 :     {
325 :     %{$slist} = %{$self->{_ss_variants}};
326 :     }
327 : olson 1.14
328 :     if (not $all)
329 :     {
330 :     for my $ss (keys %$slist)
331 :     {
332 :     my $var = $slist->{$ss};
333 : paarmann 1.32 delete $slist->{$ss} if $var eq 0 or $var eq -1;
334 : olson 1.14 }
335 :     }
336 :     return $slist;
337 :     }
338 :    
339 : overbeek 1.19 sub peg_to_subsystems
340 :     {
341 :     my($self, $peg) = @_;
342 :    
343 :     my $variant = $self->{_ss_variants};
344 :     my %sub = map { $_ => 1 }
345 :     grep { $variant->{$_} !~ /^(-1)|0$/ }
346 :     map { $_->[0] }
347 :     $self->peg_to_roles_in_subsystems($peg);
348 :     return sort keys(%sub);
349 :     }
350 :    
351 : olson 1.9 sub peg_to_roles_in_subsystems
352 :     {
353 :     my($self,$peg) = @_;
354 :    
355 :     my $fig = $self->{_fig};
356 :     my $newG = $self->{_genome};
357 :     my $newGdir = $self->{_orgdir};
358 : olson 1.8
359 : olson 1.22 if ($peg !~ /fig\|$newG\.peg/)
360 :     {
361 :     return $fig->peg_to_roles_in_subsystems($peg);
362 :     }
363 :    
364 : olson 1.9 $self->load_ss_data();
365 : olson 1.8
366 : olson 1.9 my $ret = $self->{_ss_peg_index}->{$peg};
367 : olson 1.8
368 : olson 1.9 return $ret ? @$ret : ();
369 :     }
370 : olson 1.8
371 : olson 1.9 sub subsystems_for_peg
372 :     {
373 :     my($self,$peg) = @_;
374 :     return $self->peg_to_roles_in_subsystems($peg);
375 : olson 1.5 }
376 : olson 1.4
377 : overbeek 1.1 sub genomes {
378 :     my($self,$complete) = @_;
379 :     my $fig = $self->{_fig};
380 :    
381 :     return ($self->{_genome},$fig->genomes($complete));
382 :     }
383 :    
384 :     sub genus_species {
385 :     my($self,$genome) = @_;
386 :    
387 :     my $fig = $self->{_fig};
388 :     my $newG = $self->{_genome};
389 :     my $newGdir = $self->{_orgdir};
390 :    
391 :     if (($genome eq $newG) && open(GENOME,"<$newGdir/GENOME"))
392 :     {
393 :     my $x = <GENOME>;
394 :     close(GENOME);
395 :     chop $x;
396 :     return $x;
397 :     }
398 :     else
399 :     {
400 :     return $fig->genus_species($genome);
401 :     }
402 :     }
403 :    
404 : overbeek 1.7 sub get_genome_assignment_data {
405 :     my($self,$genome) = @_;
406 :    
407 :     my $fig = $self->{_fig};
408 :     my $newG = $self->{_genome};
409 :     my $newGdir = $self->{_orgdir};
410 :    
411 :     if ($genome eq $newG)
412 :     {
413 : olson 1.31 my @fnfiles = <$newGdir/proposed*functions>;
414 :     if (@fnfiles == 0)
415 :     {
416 :     @fnfiles = <$newGdir/assigned*functions>;
417 :     }
418 :    
419 :     if (@fnfiles == 0)
420 :     {
421 :     return [];
422 :     }
423 : overbeek 1.7 my %assign = map { ( $_ =~ /^(fig\|\d+\.\d+\.peg\.\d+)\t(\S.*\S)/) ? ($1 => $2) : () }
424 : olson 1.31 `cat @fnfiles`;
425 : overbeek 1.7 return [map { [$_,$assign{$_}] } sort { &FIG::by_fig_id($a,$b) } keys(%assign)];
426 :     }
427 :     else
428 :     {
429 :     return $fig->get_genome_assignment_data($genome);
430 :     }
431 :     }
432 :    
433 :     sub org_of {
434 :     my($self,$peg) = @_;
435 :    
436 :     if ($peg =~ /^fig\|(\d+\.\d+)\.peg\.\d+/)
437 :     {
438 :     return $self->genus_species($1);
439 :     }
440 :     return "";
441 :     }
442 :    
443 :     sub get_genome_subsystem_data {
444 :     my($self,$genome) = @_;
445 :    
446 :     my $fig = $self->{_fig};
447 :     my $newG = $self->{_genome};
448 :     my $newGdir = $self->{_orgdir};
449 :    
450 :     if ($genome eq $newG)
451 :     {
452 :     return [map { ($_ =~ /^(\S[^\t]+\S)\t(\S[^\t]*\S)\t(\S+)/) ? [$1,$2,$3] : () }
453 :     `cat $newGdir/Subsystems/bindings`];
454 :     }
455 :     else
456 :     {
457 :     return $fig->get_genome_subsystem_data($genome);
458 :     }
459 :     }
460 :    
461 :     sub orgname_of_orgid {
462 :     my($self,$genome) = @_;
463 :    
464 :     return $self->genus_species($genome);
465 :     }
466 :    
467 :     sub orgid_of_orgname {
468 :     my($self,$genome_name) = @_;
469 :    
470 :     my @genomes = $self->genomes('complete');
471 :     my $i;
472 :     for ($i=0; ($i < @genomes) && ($genome_name ne $self->genus_species($genomes[$i])); $i++) {}
473 :     return ($i == @genomes) ? undef : $genomes[$i];
474 :     }
475 :    
476 :     sub genus_species_domain {
477 :     my($self,$genome) = @_;
478 :    
479 :     return [$self->genus_species($genome),$self->genome_domain($genome)];
480 :     }
481 :    
482 :     sub protein_subsystem_to_roles {
483 :     my ($self,$peg,$subsystem) = @_;
484 :    
485 :     my $fig = $self->{_fig};
486 :     my $newG = $self->{_genome};
487 :     my $newGdir = $self->{_orgdir};
488 :    
489 :     if (&FIG::genome_of($peg) ne $newG)
490 :     {
491 :     return $fig->protein_subsystem_to_roles($peg,$subsystem);
492 :     }
493 :     else
494 :     {
495 :     my @roles = map { (($_ =~ /^([^\t]+)\t([^\t]+)\t(\S+)$/) && ($1 eq $subsystem) && ($3 eq $peg)) ?
496 :     $2 : () } `cat $newGdir/Subsystems/bindings`;
497 :     my %roles = map { $_ => 1 } @roles;
498 :     return [sort keys(%roles)];
499 :     }
500 :     }
501 :    
502 :     sub contig_lengths {
503 :     my ($self, $genome) = @_;
504 :    
505 :     my $fig = $self->{_fig};
506 :     my $newG = $self->{_genome};
507 :     my $newGdir = $self->{_orgdir};
508 :    
509 :     if ($genome ne $newG)
510 :     {
511 :     return $fig->contig_lengths($genome);
512 :     }
513 :     else
514 :     {
515 : olson 1.37 my $contig_lengths = $self->{_contig_len_index};
516 :    
517 : olson 1.39 if (!tied(%$contig_lengths))
518 : overbeek 1.7 {
519 : olson 1.37 $contig_lengths = {};
520 :     if (open(CONTIGS,"<$newGdir/contigs"))
521 :     {
522 :     $/ = "\n>";
523 :     while (defined(my $x = <CONTIGS>))
524 :     {
525 :     if ($x =~ />?(\S+)[^\n]*\n(.*)/s)
526 :     {
527 :     my $id = $1;
528 :     my $seq = $2;
529 :     $seq =~ s/\s//gs;
530 :     $contig_lengths->{$id} = length($seq);
531 :     }
532 :     }
533 :     close(CONTIGS);
534 :     $/ = "\n";
535 :     }
536 : overbeek 1.7 }
537 :     return $contig_lengths;
538 :     }
539 :     }
540 :    
541 : olson 1.10 sub contig_ln {
542 :     my ($self, $genome, $contig) = @_;
543 :    
544 :     my $fig = $self->{_fig};
545 :     my $newG = $self->{_genome};
546 :     my $newGdir = $self->{_orgdir};
547 :    
548 :     if ($genome ne $newG)
549 :     {
550 :     return $fig->contig_ln($genome, $contig);
551 :     }
552 :     else
553 :     {
554 : olson 1.37 my $contig_len = $self->{_contig_len_index};
555 :    
556 : olson 1.39 if (tied(%$contig_len))
557 : olson 1.37 {
558 :     return $contig_len->{$contig};
559 :     }
560 :    
561 : olson 1.10 if (open(CONTIGS,"<$newGdir/contigs"))
562 :     {
563 :     local $/ = "\n>";
564 :     while (defined(my $x = <CONTIGS>))
565 :     {
566 : olson 1.33 chomp $x;
567 : olson 1.10 if ($x =~ />?(\S+)[^\n]*\n(.*)/s)
568 :     {
569 :     my $id = $1;
570 :     my $seq = $2;
571 :     $seq =~ s/\s//gs;
572 :     if ($id eq $contig)
573 :     {
574 :     return length($seq);
575 :     }
576 :     }
577 :     }
578 :     close(CONTIGS);
579 :     }
580 :     }
581 :     }
582 :    
583 :     sub contigs_of
584 :     {
585 :     my ($self, $genome) = @_;
586 :    
587 :     my $fig = $self->{_fig};
588 :     my $newG = $self->{_genome};
589 :     my $newGdir = $self->{_orgdir};
590 :    
591 :     if ($genome ne $newG)
592 :     {
593 :     return $fig->contigs_of($genome);
594 :     }
595 :     else
596 :     {
597 :     my @out;
598 : olson 1.37 $self->load_contigs_index();
599 : olson 1.36
600 :     my $contigs = $self->{_contigs_index};
601 : olson 1.39 if (tied(%$contigs))
602 : olson 1.36 {
603 :     return keys %$contigs;
604 :     }
605 :    
606 : olson 1.10 if (open(CONTIGS,"<$newGdir/contigs"))
607 :     {
608 :     local $/ = "\n>";
609 :     while (defined(my $x = <CONTIGS>))
610 :     {
611 : olson 1.33 chomp $x;
612 : olson 1.10 if ($x =~ />?(\S+)[^\n]*\n(.*)/s)
613 :     {
614 :     my $id = $1;
615 :     push(@out, $id);
616 :     }
617 :     }
618 :     close(CONTIGS);
619 :     }
620 :     return @out;
621 :     }
622 :     }
623 :    
624 :     =head3 dna_seq
625 :    
626 :     usage: $seq = dna_seq($genome,@locations)
627 :    
628 :     Returns the concatenated subsequences described by the list of locations. Each location
629 :     must be of the form
630 :    
631 :     Contig_Beg_End
632 :    
633 :     where Contig must be the ID of a contig for genome $genome. If Beg > End the location
634 :     describes a stretch of the complementary strand.
635 :    
636 :     =cut
637 :     #: Return Type $;
638 :     sub dna_seq {
639 :     my($self,$genome,@locations) = @_;
640 :    
641 :     my $fig = $self->{_fig};
642 :     my $newG = $self->{_genome};
643 :     my $newGdir = $self->{_orgdir};
644 :    
645 :     if ($genome ne $newG)
646 :     {
647 :     return $fig->dna_seq($genome, @locations);
648 :     }
649 :    
650 : olson 1.36 my $contigs = $self->{_contigs_index};
651 :     if (!$contigs)
652 : olson 1.10 {
653 : olson 1.36 $contigs = {};
654 :     if (open(CONTIGS,"<$newGdir/contigs"))
655 : olson 1.10 {
656 : olson 1.36 local $/ = "\n>";
657 :     while (defined(my $x = <CONTIGS>))
658 : olson 1.10 {
659 : olson 1.36 chomp $x;
660 :     if ($x =~ />?(\S+)[^\n]*\n(.*)/s)
661 :     {
662 :     my $id = $1;
663 :     my $seq = $2;
664 :     $seq =~ s/\s//gs;
665 :     $contigs->{$id} = $seq;
666 :     }
667 : olson 1.10 }
668 : olson 1.36 close(CONTIGS);
669 : olson 1.10 }
670 :     }
671 :    
672 :     my(@pieces,$loc,$contig,$beg,$end,$ln,$rdbH);
673 :    
674 :     @locations = map { split(/,/,$_) } @locations;
675 :     @pieces = ();
676 :     foreach $loc (@locations)
677 :     {
678 :     if ($loc =~ /^(\S+)_(\d+)_(\d+)$/)
679 :     {
680 :     ($contig,$beg,$end) = ($1,$2,$3);
681 : olson 1.36 my $seq = $contigs->{$contig};
682 : olson 1.10
683 :     $ln = length($seq);
684 :    
685 :     if (! $ln) {
686 :     print STDERR "$genome/$contig: could not get length\n";
687 :     return "";
688 :     }
689 :    
690 :     if (&FIG::between(1,$beg,$ln) && &FIG::between(1,$end,$ln))
691 :     {
692 :     if ($beg < $end)
693 :     {
694 : paczian 1.15 push(@pieces, substr($seq, $beg - 1, ($end - $beg) + 1));
695 : olson 1.10 }
696 :     else
697 :     {
698 : paczian 1.15 push(@pieces, &FIG::reverse_comp(substr($seq, $end - 1, ($beg - $end) + 1)));
699 : olson 1.10 }
700 :     }
701 :     }
702 :     }
703 :     return lc(join("",@pieces));
704 :     }
705 :    
706 : overbeek 1.7 sub genome_szdna {
707 :     my ($self, $genome) = @_;
708 :    
709 :     my $fig = $self->{_fig};
710 :     my $newG = $self->{_genome};
711 :     my $newGdir = $self->{_orgdir};
712 :    
713 :     if ($genome ne $newG)
714 :     {
715 :     return $fig->genome_szdna($genome);
716 :     }
717 :     else
718 :     {
719 :     my $contig_lens = $self->contig_lengths($genome);
720 :     my $tot = 0;
721 :     while ( my($contig,$len) = each %$contig_lens)
722 :     {
723 :     $tot += $len;
724 :     }
725 :     return $tot;
726 :     }
727 :     }
728 :    
729 :     sub genome_version {
730 :     my ($self, $genome) = @_;
731 :    
732 :     my $fig = $self->{_fig};
733 :     my $newG = $self->{_genome};
734 :     my $newGdir = $self->{_orgdir};
735 :    
736 :     if ($genome ne $newG)
737 :     {
738 :     return $fig->genome_version($genome);
739 :     }
740 :     else
741 :     {
742 :     return "$genome.0";
743 :     }
744 :     }
745 :    
746 :     sub genome_pegs {
747 :     my ($self, $genome) = @_;
748 :    
749 :     my $fig = $self->{_fig};
750 :     my $newG = $self->{_genome};
751 :     my $newGdir = $self->{_orgdir};
752 :    
753 :     if ($genome ne $newG)
754 :     {
755 :     return $fig->genome_pegs($genome);
756 :     }
757 :     else
758 :     {
759 :     my @tmp = $self->all_features($genome,"peg");
760 :     my $n = @tmp;
761 :     return $n;
762 :     }
763 :     }
764 :    
765 :     sub genome_rnas {
766 :     my ($self, $genome) = @_;
767 :    
768 :     my $fig = $self->{_fig};
769 :     my $newG = $self->{_genome};
770 :     my $newGdir = $self->{_orgdir};
771 :    
772 :     if ($genome ne $newG)
773 :     {
774 :     return $fig->genome_rnas($genome);
775 :     }
776 :     else
777 :     {
778 :     my @tmp = $self->all_features($genome,"rna");
779 :     my $n = @tmp;
780 :     return $n;
781 :     }
782 :     }
783 :    
784 :     sub genome_domain {
785 :     my ($self, $genome) = @_;
786 :    
787 :     my $fig = $self->{_fig};
788 :     my $newG = $self->{_genome};
789 :     my $newGdir = $self->{_orgdir};
790 :    
791 :     if ($genome ne $newG)
792 :     {
793 :     return $fig->genome_domain($genome);
794 :     }
795 :     else
796 :     {
797 :     my $tax = $self->taxonomy_of($genome);
798 :     return ($tax =~ /^([^ \t;]+)/) ? $1 : "unknown";
799 :     }
800 :     }
801 : overbeek 1.2
802 :     sub genes_in_region {
803 :     my($self,$genome,$contig,$beg,$end) = @_;
804 :    
805 :     my $fig = $self->{_fig};
806 :     my $newG = $self->{_genome};
807 :     my $newGdir = $self->{_orgdir};
808 :    
809 :     if ($genome ne $newG)
810 :     {
811 :     return $fig->genes_in_region($genome,$contig,$beg,$end);
812 :     }
813 :     else
814 :     {
815 :     &load_tbl($self);
816 :     my $tblH = $self->{_tbl};
817 :     my $maxV = 0;
818 :     my $minV = 1000000000;
819 :     my $genes = [];
820 :     while ( my($fid,$tuple) = each %$tblH)
821 :     {
822 :     if (($tuple->[0]->[0] =~ /^(\S+)_(\d+)_\d+$/) && ($1 eq $contig))
823 :     {
824 :     my $beg1 = $2;
825 :     my $last = @{$tuple->[0]} - 1;
826 :     if (($tuple->[0]->[$last] =~ /^(\S+)_\d+_(\d+)$/) && ($1 eq $contig))
827 :     {
828 :     my $end1 = $2;
829 :     if (&overlaps($beg,$end,$beg1,$end1))
830 :     {
831 :     $minV = &FIG::min($minV,$beg1,$end1);
832 :     $maxV = &FIG::max($maxV,$beg1,$end1);
833 :     push(@$genes,$fid);
834 :     }
835 :     }
836 :     }
837 :     }
838 :     return ($genes,$minV,$maxV);
839 :     }
840 :     }
841 :    
842 :     sub overlaps {
843 :     my($b1,$e1,$b2,$e2) = @_;
844 :    
845 :     if ($b1 > $e1) { ($b1,$e1) = ($e1,$b1) }
846 :     if ($b2 > $e2) { ($b2,$e2) = ($e2,$b2) }
847 :     return &FIG::between($b1,$b2,$e1) || &FIG::between($b2,$b1,$e2);
848 :     }
849 :    
850 : overbeek 1.7 sub all_contigs {
851 :     my($self,$genome) = @_;
852 :    
853 :     my $fig = $self->{_fig};
854 :     my $newG = $self->{_genome};
855 :     my $newGdir = $self->{_orgdir};
856 :    
857 :     if ($genome ne $newG)
858 :     {
859 :     return $fig->all_contigs($genome);
860 :     }
861 :     else
862 :     {
863 :     &load_tbl($self);
864 :     my %contigs;
865 :     my $tblH = $self->{_tbl};
866 :     while ( my($fid,$tuple) = each %$tblH)
867 :     {
868 :     if ($tuple->[0]->[0] =~ /^(\S+)_\d+_\d+$/)
869 :     {
870 :     $contigs{$1} = 1;
871 :     }
872 :     }
873 :     return keys(%contigs);
874 :     }
875 :     }
876 :    
877 :     sub all_features {
878 :     my($self,$genome,$type) = @_;
879 :    
880 :     my $fig = $self->{_fig};
881 :     my $newG = $self->{_genome};
882 :     my $newGdir = $self->{_orgdir};
883 :    
884 :     if ($genome ne $newG)
885 :     {
886 :     return $fig->all_features($genome,$type);
887 :     }
888 :     else
889 :     {
890 :     &load_tbl($self);
891 :     my $tblH = $self->{_tbl};
892 :     return sort { &FIG::by_fig_id($a,$b) }
893 :     grep { ($_ =~ /^fig\|\d+\.\d+\.([^\.]+)/) && ($1 eq $type) }
894 :     keys(%$tblH);
895 :     }
896 :     }
897 :    
898 :     sub all_features_detailed_fast {
899 : olson 1.11 my($self,$genome, $regmin, $regmax, $contig) = @_;
900 : overbeek 1.7
901 :     my $fig = $self->{_fig};
902 :     my $newG = $self->{_genome};
903 :     my $newGdir = $self->{_orgdir};
904 :    
905 :     if ($genome ne $newG)
906 :     {
907 : olson 1.11 return $fig->all_features_detailed_fast($genome, $regmin, $regmax, $contig);
908 : overbeek 1.7 }
909 :     else
910 :     {
911 :     &load_tbl($self);
912 :     my $tblH = $self->{_tbl};
913 :     my $feat_details = [];
914 :     while ( my($fid,$tuple) = each %$tblH)
915 :     {
916 :     if ($fid =~ /^fig\|\d+\.\d+\.(\S+)\.\d+/)
917 :     {
918 :     my $type = $1;
919 : olson 1.11 if ($tuple->[0]->[0] =~ /^(\S+)_(\d+)_(\d+)$/)
920 : overbeek 1.7 {
921 : olson 1.11 my($ctg, $beg, $end) = ($1, $2, $3);
922 :     next if (defined($contig) and $contig ne $ctg);
923 :    
924 : overbeek 1.7 my($min,$max);
925 : olson 1.11 if ($beg < $end)
926 : overbeek 1.7 {
927 : olson 1.11 $min = $beg;
928 :     $max = $end;
929 : overbeek 1.7 }
930 :     else
931 :     {
932 : olson 1.11 $min = $end;
933 :     $max = $beg;
934 :     }
935 :    
936 :     if (not defined($regmin) or not defined($regmax) or
937 :     ($min < $regmin and $regmin < $max) or ($min < $regmax and $regmax < $max) or ($min > $regmin and $max < $regmax))
938 :     {
939 :     push(@$feat_details,[$fid,$tuple->[0]->[0],join(",",@{$tuple->[1]}),$type,$min,$max,$self->function_of($fid),'master','']);
940 : overbeek 1.7 }
941 :     }
942 :     }
943 :     }
944 :     return $feat_details;
945 :     }
946 :     }
947 :    
948 :     sub compute_clusters {
949 :     # Get the parameters.
950 :     my ($self, $pegList, $subsystem, $distance) = @_;
951 :     if (! defined $distance) {
952 :     $distance = 5000;
953 :     }
954 :    
955 :     my($peg,%by_contig);
956 :     foreach $peg (@$pegList)
957 :     {
958 :     my $loc;
959 :     if ($loc = $self->feature_location($peg))
960 :     {
961 :     my ($contig,$beg,$end) = &FIG::boundaries_of($loc);
962 :     my $genome = &FIG::genome_of($peg);
963 :     push(@{$by_contig{"$genome\t$contig"}},[($beg+$end)/2,$peg]);
964 :     }
965 :     }
966 :    
967 :     my @clusters = ();
968 :     foreach my $tuple (keys(%by_contig))
969 :     {
970 :     my $x = $by_contig{$tuple};
971 :     my @pegs = sort { $a->[0] <=> $b->[0] } @$x;
972 :     while ($x = shift @pegs)
973 :     {
974 :     my $clust = [$x->[1]];
975 :     while ((@pegs > 0) && (abs($pegs[0]->[0] - $x->[0]) <= $distance))
976 :     {
977 :     $x = shift @pegs;
978 :     push(@$clust,$x->[1]);
979 :     }
980 :    
981 :     if (@$clust > 1)
982 :     {
983 :     push(@clusters,$clust);
984 :     }
985 :     }
986 :     }
987 :     return sort { @$b <=> @$a } @clusters;
988 :     }
989 :    
990 :    
991 : overbeek 1.1 sub feature_location {
992 :     my($self,$fid) = @_;
993 :    
994 :     my $fig = $self->{_fig};
995 :     my $newG = $self->{_genome};
996 :     my $newGdir = $self->{_orgdir};
997 :    
998 :     if (($fid =~ /^fig\|(\d+\.\d+)/) && ($1 eq $newG))
999 :     {
1000 :     &load_tbl($self);
1001 :     if (my $x = $self->{_tbl}->{$fid})
1002 :     {
1003 :     return join(",",@{$x->[0]});
1004 :     }
1005 :     else
1006 :     {
1007 :     return undef;
1008 :     }
1009 :     }
1010 :     else
1011 :     {
1012 :     return scalar $fig->feature_location($fid);
1013 :     }
1014 :     }
1015 :    
1016 :     sub function_of {
1017 :     my($self,$fid) = @_;
1018 :    
1019 :     my $fig = $self->{_fig};
1020 :     my $newG = $self->{_genome};
1021 :     my $newGdir = $self->{_orgdir};
1022 :    
1023 :     if (($fid =~ /^fig\|(\d+\.\d+)/) && ($1 eq $newG))
1024 :     {
1025 :     &load_functions($self);
1026 :     return $self->{_functions}->{$fid};
1027 :     }
1028 :     else
1029 :     {
1030 :     return scalar $fig->function_of($fid);
1031 :     }
1032 :     }
1033 :    
1034 :     sub feature_aliases {
1035 :     my($self,$fid) = @_;
1036 :    
1037 :     my $fig = $self->{_fig};
1038 :     my $newG = $self->{_genome};
1039 :     my $newGdir = $self->{_orgdir};
1040 :    
1041 :     my @aliases;
1042 :     if (($fid =~ /^fig\|(\d+\.\d+)/) && ($1 eq $newG))
1043 :     {
1044 :     &load_tbl($self);
1045 :     if (my $x = $self->{_tbl}->{$fid})
1046 :     {
1047 :     @aliases = @{$x->[1]};
1048 :     }
1049 :     else
1050 :     {
1051 :     @aliases = ();
1052 :     }
1053 :     }
1054 :     else
1055 :     {
1056 :     @aliases = $fig->feature_aliases($fid);
1057 :     }
1058 :     return wantarray() ? @aliases : join(",",@aliases);
1059 :     }
1060 :    
1061 :     sub feature_annotations {
1062 :     my($self,$fid,$rawtime) = @_;
1063 :    
1064 :     my $fig = $self->{_fig};
1065 :     my $newG = $self->{_genome};
1066 :     my $newGdir = $self->{_orgdir};
1067 :    
1068 :     my @annotations;
1069 :     if (($fid =~ /^fig\|(\d+\.\d+)/) && ($1 eq $newG))
1070 :     {
1071 :     &load_ann($self);
1072 :     if (my $x = $self->{_ann}->{$fid})
1073 :     {
1074 :     @annotations = @{$x};
1075 :     }
1076 :     else
1077 :     {
1078 :     @annotations = ();
1079 :     }
1080 :    
1081 :     if ($rawtime)
1082 :     {
1083 :     return @annotations;
1084 :     }
1085 :     else
1086 :     {
1087 :     return map { $_->[1] = localtime($_->[1]); $_ } @annotations;
1088 :     }
1089 :     }
1090 :     else
1091 :     {
1092 :     return $fig->feature_annotations($fid);
1093 :     }
1094 :     }
1095 :    
1096 :     sub get_translation {
1097 :     my($self,$peg) = @_;
1098 :    
1099 :     my $fig = $self->{_fig};
1100 :     my $newG = $self->{_genome};
1101 :     my $newGdir = $self->{_orgdir};
1102 :    
1103 :     if (($peg =~ /^fig\|(\d+\.\d+)/) && ($1 eq $newG))
1104 :     {
1105 :     &load_pseq($self);
1106 :     return $self->{_pseq}->{$peg};
1107 :     }
1108 :     else
1109 :     {
1110 :     return $fig->get_translation($peg);
1111 :     }
1112 :     }
1113 :    
1114 : olson 1.4 sub translation_length
1115 :     {
1116 :     my($self, $peg) = @_;
1117 :    
1118 :     my $fig = $self->{_fig};
1119 :     my $newG = $self->{_genome};
1120 :     my $newGdir = $self->{_orgdir};
1121 :    
1122 :     if (($peg =~ /^fig\|(\d+\.\d+)/) && ($1 eq $newG))
1123 :     {
1124 :     &load_pseq($self);
1125 :     return length($self->{_pseq}->{$peg});
1126 :     }
1127 :     else
1128 :     {
1129 :     return $fig->translation_length($peg);
1130 :     }
1131 :     }
1132 :    
1133 :     sub translatable
1134 :     {
1135 :     my($self, $peg) = @_;
1136 :    
1137 :     my $fig = $self->{_fig};
1138 :     my $newG = $self->{_genome};
1139 :     my $newGdir = $self->{_orgdir};
1140 :    
1141 :     if (($peg =~ /^fig\|(\d+\.\d+)/) && ($1 eq $newG))
1142 :     {
1143 :     &load_pseq($self);
1144 :     return length($self->{_pseq}->{$peg}) > 0 ? 1 : 0;;
1145 :     }
1146 :     else
1147 :     {
1148 :     return $fig->translatable($peg);
1149 :     }
1150 :     }
1151 :    
1152 :     sub is_real_feature
1153 :     {
1154 :     my($self, $fid) = @_;
1155 :    
1156 :     if ($self->is_virtual_feature($fid))
1157 :     {
1158 :     return 1;
1159 :     }
1160 :     else
1161 :     {
1162 :     return $self->{_fig}->is_real_feature($fid);
1163 :     }
1164 :    
1165 :     }
1166 :    
1167 : olson 1.10 sub pegs_of
1168 :     {
1169 :     my($self, $genome) = @_;
1170 :    
1171 :     my $fig = $self->{_fig};
1172 :     my $newG = $self->{_genome};
1173 :     my $newGdir = $self->{_orgdir};
1174 :    
1175 :     if ($genome ne $newG)
1176 :     {
1177 :     return $fig->pegs_of($genome);
1178 :     }
1179 :    
1180 :     $self->load_tbl();
1181 :     return grep { /\.peg\./ } keys %{$self->{_tbl}};
1182 :     }
1183 :    
1184 :    
1185 :     sub rnas_of
1186 :     {
1187 :     my($self, $genome) = @_;
1188 :    
1189 :     my $fig = $self->{_fig};
1190 :     my $newG = $self->{_genome};
1191 :     my $newGdir = $self->{_orgdir};
1192 :    
1193 :     if ($genome ne $newG)
1194 :     {
1195 :     return $fig->pegs_of($genome);
1196 :     }
1197 :    
1198 :     $self->load_tbl();
1199 :     return grep { /\.rna\./ } keys %{$self->{_tbl}};
1200 :     }
1201 :    
1202 : olson 1.4 sub is_virtual_feature
1203 :     {
1204 :     my($self, $peg) = @_;
1205 :    
1206 :     my $fig = $self->{_fig};
1207 :     my $newG = $self->{_genome};
1208 :     my $newGdir = $self->{_orgdir};
1209 :    
1210 :     if (($peg =~ /^fig\|(\d+\.\d+)/) && ($1 eq $newG))
1211 :     {
1212 :     &load_pseq($self);
1213 :     return 1;
1214 :     }
1215 :     else
1216 :     {
1217 :     return 0;
1218 :     }
1219 :     }
1220 :    
1221 : olson 1.35 ####################################################
1222 :     #
1223 :     # Following are some MG-RAST specific features. FIGV seems a good a place as any to include them.#
1224 :     #
1225 :     #
1226 :    
1227 :     =head3 taxa_to_seed_ids
1228 :    
1229 :     Given a prefix of a taxonomy string, return the list of metagenome fids that
1230 :     mapped to SEED organisms in that taxonomy.
1231 :    
1232 :     =cut
1233 :    
1234 :     sub taxa_to_seed_ids
1235 :     {
1236 :     my($self, $tax) = @_;
1237 :    
1238 :     my $btree_file = "$self->{_orgdir}/taxa_summary_by_blast.btree";
1239 :     my %btree;
1240 :    
1241 :     my $btree_tie = tie %btree, 'DB_File', $btree_file, O_RDONLY, 0666, $DB_BTREE;
1242 :    
1243 :     if (!$btree_tie)
1244 :     {
1245 :     warn "Cannot tie $btree_file: $!\n";
1246 :     return undef;
1247 :     }
1248 :    
1249 :     my $key = $tax;
1250 :     my ($res, $value);
1251 :     my @vals;
1252 :     my %seen;
1253 :    
1254 :     for ($res = $btree_tie->seq($key, $value, R_CURSOR);
1255 :     $res == 0 and $key =~ /^$tax/;
1256 :     $res = $btree_tie->seq($key, $value, R_NEXT))
1257 :     {
1258 :     print "$key\n";
1259 :     push(@vals, map { [$key, $_] } grep { not $seen{$_}++ } split(/$;/, $value));
1260 :     }
1261 :     return @vals;
1262 :     }
1263 :    
1264 : overbeek 1.1 sub load_pseq {
1265 :     my($self) = @_;
1266 :    
1267 :     if ($self->{_pseq}) { return };
1268 :    
1269 :     my $newGdir = $self->{_orgdir};
1270 :     my $pseq = {};
1271 :     if (open(FASTA,"<$newGdir/Features/peg/fasta"))
1272 :     {
1273 : olson 1.10 local $/ = "\n>";
1274 : overbeek 1.1 my $x;
1275 :     while (defined($x = <FASTA>))
1276 :     {
1277 : olson 1.21 chomp $x;
1278 : olson 1.10 if ($x =~ />?(\S+)[^\n]*\n(.*)/s)
1279 : overbeek 1.1 {
1280 :     my $peg = $1;
1281 :     my $seq = $2;
1282 :     $seq =~ s/\s//gs;
1283 :     $pseq->{$peg} = $seq;
1284 :     }
1285 :     }
1286 :     close(FASTA);
1287 :     }
1288 :     $self->{_pseq} = $pseq;
1289 :     }
1290 :    
1291 : olson 1.9 sub load_ss_data
1292 :     {
1293 :     my($self) = @_;
1294 :    
1295 :     return if defined($self->{_ss_bindings});
1296 :    
1297 :     my $fig = $self->{_fig};
1298 :     my $newG = $self->{_genome};
1299 :     my $newGdir = $self->{_orgdir};
1300 :    
1301 :     open(S, "<$newGdir/Subsystems/bindings") or die "Cannot open $newGdir/Subsystems/bindings: $!";
1302 :    
1303 :     my $peg_index;
1304 :     my $bindings;
1305 : olson 1.14 my $genome_index;
1306 :     my %seen;
1307 : olson 1.9 while (<S>)
1308 :     {
1309 :     chomp;
1310 :     my($sname, $role, $peg) = split(/\t/);
1311 : olson 1.14
1312 : olson 1.34 my($genome) = ($peg =~ /fig\|(\d+\.\d+)/);
1313 :    
1314 :     # my $genome = &FIG::genome_of($peg);
1315 : olson 1.9
1316 :     push(@{$bindings->{$sname}->{$role}}, $peg);
1317 :     push(@{$peg_index->{$peg}}, [$sname, $role]);
1318 : olson 1.14
1319 :     if (!$seen{$genome, $sname, $role})
1320 :     {
1321 : parrello 1.24 push(@{$genome_index->{$genome}}, [$sname, $role]);
1322 : olson 1.14 }
1323 :    
1324 : olson 1.9 }
1325 :     close(S);
1326 :    
1327 :     open(S, "<$newGdir/Subsystems/subsystems") or die "Cannot open $newGdir/Subsystems/subsystems: $!";
1328 :     my $variant;
1329 :     while (<S>)
1330 :     {
1331 :     chomp;
1332 :     my($sname, $v) = split(/\t/);
1333 :     $variant->{$sname} = $v;
1334 :     }
1335 :     close(S);
1336 :    
1337 :     $self->{_ss_bindings} = $bindings;
1338 :     $self->{_ss_variants} = $variant;
1339 :     $self->{_ss_peg_index} = $peg_index;
1340 : olson 1.14 $self->{_ss_genome_index} = $genome_index;
1341 : olson 1.9 }
1342 :    
1343 : overbeek 1.1 sub load_tbl {
1344 :     my($self) = @_;
1345 :    
1346 :     if ($self->{_tbl}) { return };
1347 :    
1348 :     my $newGdir = $self->{_orgdir};
1349 :     my $tbl = {};
1350 :     foreach my $x (`cat $newGdir/Features/*/tbl`)
1351 :     {
1352 :     if ($x =~ /^(\S+)\t(\S+)(\t(\S.*\S))?/)
1353 :     {
1354 :     my $fid = $1;
1355 :     my $loc = [split(/,/,$2)];
1356 : paczian 1.12 my $aliases = $4 ? [split(/\t/,$4)] : [];
1357 : overbeek 1.1 $tbl->{$fid} = [$loc,$aliases];
1358 :     }
1359 :     }
1360 :     $self->{_tbl} = $tbl;
1361 :     }
1362 :    
1363 :     sub load_functions {
1364 :     my($self) = @_;
1365 :    
1366 :     if ($self->{_functions}) { return };
1367 :    
1368 :     my $newG = $self->{_genome};
1369 :     my $newGdir = $self->{_orgdir};
1370 :    
1371 :     my $functions = {};
1372 :     foreach my $x (`cat $newGdir/*functions`)
1373 :     {
1374 :     if (($x =~ /^(fig\|(\d+\.\d+)\.\S+)\t(\S[^\t]*\S)/) && ($2 eq $newG))
1375 :     {
1376 :     $functions->{$1} = $3;
1377 :     }
1378 :     }
1379 :     $self->{_functions} = $functions;
1380 :     }
1381 :    
1382 : olson 1.16
1383 :     sub bbhs
1384 :     {
1385 :     my($self,$peg,$cutoff) = @_;
1386 :    
1387 :     my $fig = $self->{_fig};
1388 :     my $newG = $self->{_genome};
1389 :     my $newGdir = $self->{_orgdir};
1390 :    
1391 :     $self->load_bbh_index();
1392 :    
1393 :     $cutoff = 1.0e-10 unless defined($cutoff);
1394 :    
1395 :     my $tie = $self->{_bbh_tie};
1396 :    
1397 :     my @bbhs = $tie->get_dup($peg);
1398 :    
1399 :     # @bbhs now a list of comma-separated tuples (id2, score, bitscore)
1400 :     # print Dumper(\@bbhs);
1401 :    
1402 :     @bbhs = grep { $_->[1] < $cutoff } map { [split(/,/)] } @bbhs;
1403 :    
1404 :     if (not ($peg =~ /^fig\|(\d+\.\d+)/ && ($1 eq $newG)))
1405 :     {
1406 :     #
1407 :     # If this isn't one of our pegs, we retrieve the bbhs from the underlying
1408 :     # SEED and merge in the bbhs that we have here.
1409 :     #
1410 :    
1411 :     my @fbbhs = $fig->bbhs($peg, $cutoff);
1412 :     # print Dumper(\@fbbhs);
1413 :     push(@bbhs, @fbbhs);
1414 :     }
1415 :     return sort { $a->[1] <=> $b->[1] } @bbhs;
1416 :     }
1417 :    
1418 : olson 1.13 sub sims
1419 :     {
1420 : olson 1.27 my($self,$pegarg,$max,$maxP,$select, $max_expand, $filters) = @_;
1421 : olson 1.13
1422 :     my $fig = $self->{_fig};
1423 :     my $newG = $self->{_genome};
1424 :     my $newGdir = $self->{_orgdir};
1425 : parrello 1.24 $max = $max ? $max : 10000;
1426 :     $maxP = $maxP ? $maxP : 1.0e-5;
1427 : olson 1.13 $select = $select ? $select : "all";
1428 :    
1429 : olson 1.27 my $prefix = "_sims";
1430 :     my $flip_prefix = "_flips";
1431 :    
1432 :     if ($select eq 'raw')
1433 :     {
1434 :     $prefix = "_raw_sims";
1435 :     $flip_prefix = "_raw_flips";
1436 :     }
1437 :    
1438 :     #
1439 :     # Partition pegs into one set that is part of this organism
1440 :     # and another set that is not.
1441 :     #
1442 :     my @pegs = ref($pegarg) ? @$pegarg : ($pegarg);
1443 :     my(@part, @not_part, %part);
1444 :     my %sims;
1445 :    
1446 :     for my $peg (@pegs)
1447 :     {
1448 :     if ($peg =~ /^fig\|(\d+\.\d+)/ and $1 eq $newG)
1449 :     {
1450 :     push(@part, $peg);
1451 :     $part{$peg}++;
1452 :     }
1453 :     else
1454 :     {
1455 :     push(@not_part, $peg);
1456 :     $sims{$peg} = [];
1457 :     }
1458 :     }
1459 :    
1460 :     #
1461 :     # Retrieve a batch of the sims for the not-part pegs, and partition
1462 :     # into %sims.
1463 :     #
1464 :     for my $sim ($fig->sims(\@not_part, $max, $maxP, $select, $max_expand, $filters))
1465 :     {
1466 :     push(@{$sims{$sim->id1}}, $sim);
1467 :     }
1468 :    
1469 :     my @out;
1470 :     my $start_len;
1471 :    
1472 :     for my $peg (@pegs)
1473 : olson 1.13 {
1474 : olson 1.27 $start_len = @out;
1475 :     if (not $part{$peg})
1476 : overbeek 1.23 {
1477 : olson 1.27 my @flips = $self->retrieve_sims($peg, $flip_prefix, $maxP, $select);
1478 :    
1479 :     @flips = sort { $b->bsc <=> $a->bsc } @flips;
1480 :    
1481 :     # my @old = $fig->sims($peg,$max,$maxP,$select);
1482 :    
1483 :     my @old = sort { $b->bsc <=> $a->bsc } @{$sims{$peg}};
1484 :    
1485 :     # my @merged = ();
1486 :     my $i1 = 0;
1487 :     my $i2 = 0;
1488 :     while (($i1 < @flips) || ($i2 < @old))
1489 : overbeek 1.23 {
1490 : olson 1.27 if (($i1 == @flips) || (($i2 < @old) && ($flips[$i1]->[11] < $old[$i2]->[11])))
1491 :     {
1492 :     # push(@merged,$old[$i2]);
1493 :     push(@out,$old[$i2]);
1494 :     $i2++;
1495 :     }
1496 :     else
1497 :     {
1498 :     # push(@merged,$flips[$i1]);
1499 :     push(@out,$flips[$i1]);
1500 :     $i1++;
1501 :     }
1502 : overbeek 1.23 }
1503 :     }
1504 : olson 1.27 else
1505 :     {
1506 :     my @sims = $self->retrieve_sims($peg, $prefix, $maxP, $select);
1507 :     push(@out, @sims);
1508 :     }
1509 :    
1510 :     if (@out - $start_len > $max)
1511 :     {
1512 :     $#out = $start_len + $max - 1;
1513 :     }
1514 : olson 1.13 }
1515 : olson 1.27
1516 :     return @out;
1517 :     }
1518 : olson 1.13
1519 : olson 1.27 sub retrieve_sims
1520 :     {
1521 :     my($self, $peg, $prefix, $maxP, $select) = @_;
1522 :    
1523 : olson 1.13 $self->load_sims_index();
1524 :    
1525 : olson 1.27 my $info = $self->{"${prefix}_index"}->{$peg};
1526 : olson 1.13
1527 :     if ($info !~ /^(\d+),(\d+)$/)
1528 :     {
1529 :     return;
1530 :     }
1531 :     my($seek, $len) = ($1, $2);
1532 :    
1533 : olson 1.27 my $sims_txt = &FIG::read_block($self->{"${prefix}_fh"}, $seek, $len);
1534 : olson 1.13
1535 :     my @sims = map { bless $_, 'Sim' } grep { $_->[10] <= $maxP } map { [split(/\t/)] } @$sims_txt;
1536 :    
1537 :     if ($select =~ /^figx?$/)
1538 :     {
1539 :     @sims = grep { $_->[1] =~ /^fig/ } @sims;
1540 :     }
1541 :    
1542 :     return @sims;
1543 :     }
1544 :    
1545 :     sub sims_old {
1546 : overbeek 1.1 my($self,$peg,$max,$maxP,$select) = @_;
1547 :    
1548 :     my $fig = $self->{_fig};
1549 :     my $newG = $self->{_genome};
1550 :     my $newGdir = $self->{_orgdir};
1551 : parrello 1.24 $max = $max ? $max : 10000;
1552 :     $maxP = $maxP ? $maxP : 1.0e-5;
1553 : overbeek 1.1 $select = $select ? $select : "all";
1554 :    
1555 :     if (($peg =~ /^fig\|(\d+\.\d+)/) && ($1 eq $newG))
1556 :     {
1557 :     &load_sims($self);
1558 :     my @sims = ();
1559 :     my $raw_sims = $self->{_sims}->{$peg};
1560 :     if ($raw_sims)
1561 :     {
1562 :     foreach my $sim ( grep { $_->[10] <= $maxP } @$raw_sims )
1563 :     {
1564 :     my $id2 = $sim->id2;
1565 :     my $id1 = $sim->id1;
1566 :     my @relevant = ();
1567 :    
1568 :     my @maps_to = $fig->mapped_prot_ids( $id2 );
1569 :     my $ref_len = $maps_to[0]->[1];
1570 :    
1571 :     @maps_to = grep { $_->[0] !~ /^xxx\d+/ } @maps_to;
1572 :    
1573 :     if ( $select =~ /^figx?$/ ) # Only fig
1574 :     {
1575 :     @relevant = grep { $_->[0] =~ /^fig/ } @maps_to;
1576 :     }
1577 :     else # All
1578 :     {
1579 :     @relevant = @maps_to;
1580 :     }
1581 :    
1582 :     my $seen = {};
1583 :     foreach my $x ( @relevant )
1584 :     {
1585 :     my ( $x_id, $x_ln ) = @$x;
1586 :    
1587 :     next if $seen->{$x_id};
1588 :     $seen->{$x_id} = 1;
1589 :    
1590 :     my $delta2 = $ref_len - $x_ln; # Coordinate shift
1591 :     my $sim1 = [ @$sim ]; # Make a copy
1592 :     $sim1->[1] = $x_id;
1593 :     $sim1->[8] -= $delta2;
1594 :     $sim1->[9] -= $delta2;
1595 :     bless( $sim1, "Sim" );
1596 :     push( @sims, $sim1 );
1597 :     }
1598 :     }
1599 :     }
1600 :    
1601 :     if (@sims > $max) { $#sims = $max-1; }
1602 :     return @sims;
1603 :     }
1604 :     else
1605 :     {
1606 :     return $fig->sims($peg,$max,$maxP,$select);
1607 :     }
1608 :     }
1609 :    
1610 : olson 1.27
1611 :     sub coupled_to
1612 :     {
1613 :     my($self,$peg) = @_;
1614 :    
1615 :     my $fig = $self->{_fig};
1616 :     my $newG = $self->{_genome};
1617 :     my $newGdir = $self->{_orgdir};
1618 :    
1619 :     if ($peg =~ /^fig\|$newG\.peg/)
1620 :     {
1621 :    
1622 :     $self->load_coupling_index();
1623 :    
1624 :     my $tie = $self->{_pch_tie};
1625 : olson 1.30
1626 :     if (not defined($tie))
1627 :     {
1628 :     return;
1629 :     }
1630 : olson 1.27
1631 :     my @dat = $tie->get_dup($peg);
1632 :    
1633 :     return map { [split(/$;/, $_)] } @dat;
1634 :     }
1635 :     else
1636 :     {
1637 :     return $fig->coupled_to($peg);
1638 :     }
1639 :    
1640 :     }
1641 :    
1642 :     sub coupling_evidence
1643 :     {
1644 :     my($self,$peg1, $peg2) = @_;
1645 :    
1646 :     my $fig = $self->{_fig};
1647 :     my $newG = $self->{_genome};
1648 :     my $newGdir = $self->{_orgdir};
1649 :    
1650 :     if ($peg1 =~ /^fig\|$newG\.peg/)
1651 :     {
1652 :     $self->load_coupling_index();
1653 :    
1654 :     my $tie = $self->{_pch_ev_tie};
1655 :    
1656 : olson 1.30 if (not defined($tie))
1657 :     {
1658 :     return;
1659 :     }
1660 :    
1661 : olson 1.27 my @dat = $tie->get_dup("$peg1$;$peg2");
1662 :    
1663 :     my @a;
1664 :     return map { @a = split(/$;/, $_); [@a[0,1,4]] } @dat;
1665 :     }
1666 :     else
1667 :     {
1668 :     return $fig->coupling_evidence($peg1, $peg2);
1669 :     }
1670 :    
1671 :     }
1672 :    
1673 :     sub coupling_and_evidence
1674 :     {
1675 :     my($self,$peg1) = @_;
1676 :    
1677 :     my $fig = $self->{_fig};
1678 :     my $newG = $self->{_genome};
1679 :     my $newGdir = $self->{_orgdir};
1680 :    
1681 :     if ($peg1 =~ /^fig\|$newG\.peg/)
1682 :     {
1683 :     $self->load_coupling_index();
1684 :    
1685 :     my $tie = $self->{_pch_tie};
1686 :     my $evtie = $self->{_pch_ev_tie};
1687 :    
1688 : olson 1.30 if (not(defined($tie) and defined($evtie)))
1689 :     {
1690 :     return;
1691 :     }
1692 :    
1693 : olson 1.27 my @out;
1694 :     my @coupled_to = $tie->get_dup($peg1);
1695 :     for my $ent (@coupled_to)
1696 :     {
1697 :     my ($peg2, $score) = split(/$;/, $ent);
1698 :    
1699 :     my @ev = $evtie->get_dup("$peg1$;$peg2");
1700 :    
1701 :     my $l = [];
1702 :     for my $event (@ev)
1703 :     {
1704 :     my($peg3, $peg4, $iden3, $iden4, $rep) = split(/$;/, $event);
1705 :     push(@$l, [$peg3, $peg4]);
1706 :     }
1707 :     push(@out, [$score, $peg2, $l]);
1708 :     }
1709 :     return @out;
1710 :     }
1711 :     else
1712 :     {
1713 :     return $fig->coupling_and_evidence($peg1);
1714 :     }
1715 :    
1716 :     }
1717 :    
1718 :     sub load_coupling_index
1719 :     {
1720 :     my($self) = @_;
1721 :    
1722 :     return if defined($self->{_pch_index});
1723 :    
1724 :     my $fig = $self->{_fig};
1725 :     my $newG = $self->{_genome};
1726 :     my $newGdir = $self->{_orgdir};
1727 :    
1728 :     my $pch_btree = "$newGdir/pchs.btree";
1729 :     my $ev_btree = "$newGdir/pchs.evidence.btree";
1730 :    
1731 :     my $pch_index = {};
1732 :     my $ev_index = {};
1733 :    
1734 :     my $tied = tie %$pch_index, 'DB_File', $pch_btree, O_RDONLY, 0666, $DB_BTREE;
1735 :     my $evtied = tie %$ev_index, 'DB_File', $ev_btree, O_RDONLY, 0666, $DB_BTREE;
1736 :    
1737 :     #
1738 :     # Set these even if failed so we don't keep trying to open and failing.
1739 :     #
1740 :     $self->{_pch_index} = $pch_index;
1741 :     $self->{_pch_tie} = $tied;
1742 :     $self->{_pch_ev_index} = $ev_index;
1743 :     $self->{_pch_ev_tie} = $evtied;
1744 :    
1745 :     if (not $tied)
1746 :     {
1747 :     warn "Cannot tie pch index $pch_btree: $!\n";
1748 :     }
1749 :     }
1750 :    
1751 : olson 1.16 sub load_bbh_index
1752 :     {
1753 :     my($self) = @_;
1754 :    
1755 :     return if defined($self->{_bbh_index});
1756 :    
1757 :     my $fig = $self->{_fig};
1758 :     my $newG = $self->{_genome};
1759 :     my $newGdir = $self->{_orgdir};
1760 :    
1761 :     my $bbh_file = "$newGdir/bbhs";
1762 :     my $bbh_index_file = "$bbh_file.index";
1763 :     my $bbh_index = {};
1764 :    
1765 :     my $tied = tie %$bbh_index, 'DB_File', $bbh_index_file, O_RDONLY, 0666, $DB_BTREE;
1766 :    
1767 :     #
1768 :     # Set these even if failed so we don't keep trying to open and failing.
1769 :     #
1770 :     $self->{_bbh_index} = $bbh_index;
1771 :     $self->{_bbh_tie} = $tied;
1772 :    
1773 :     if (not $tied)
1774 :     {
1775 :     warn "Cannot tie bbh index $bbh_index_file: $!\n";
1776 :     }
1777 :     }
1778 :    
1779 : olson 1.36 sub load_contigs_index
1780 :     {
1781 :     my($self) = @_;
1782 :    
1783 :     return if defined($self->{_contigs_index});
1784 :    
1785 :     my $fig = $self->{_fig};
1786 :     my $newG = $self->{_genome};
1787 :     my $newGdir = $self->{_orgdir};
1788 :    
1789 :     my $contig_index_file = "$newGdir/contigs.btree";
1790 :     my $contig_index = {};
1791 : olson 1.37 my $contig_len_index_file = "$newGdir/contig_len.btree";
1792 :     my $contig_len_index = {};
1793 : olson 1.36
1794 :     my $tied = tie %$contig_index, 'DB_File', $contig_index_file, O_RDONLY, 0666, $DB_BTREE;
1795 : olson 1.37 if (not $tied)
1796 :     {
1797 :     warn "Cannot tie contig index $contig_index_file: $!\n";
1798 :     }
1799 :    
1800 :     my $ltied = tie %$contig_len_index, 'DB_File', $contig_len_index_file, O_RDONLY, 0666, $DB_BTREE;
1801 :     if (not $ltied)
1802 :     {
1803 :     warn "Cannot tie contig length index $contig_len_index_file: $!\n";
1804 :     }
1805 : olson 1.36
1806 :     #
1807 :     # Set these even if failed so we don't keep trying to open and failing.
1808 :     #
1809 : olson 1.37 $self->{_contigs_index} = $contig_index;
1810 :     $self->{_contigs_tie} = $tied;
1811 :     $self->{_contig_len_index} = $contig_len_index;
1812 :     $self->{_contig_len_tie} = $ltied;
1813 : olson 1.36
1814 :     }
1815 :    
1816 : olson 1.13 sub load_sims_index
1817 :     {
1818 :     my($self) = @_;
1819 :    
1820 :     return if defined($self->{_sims_index});
1821 :    
1822 :     my $fig = $self->{_fig};
1823 :     my $newG = $self->{_genome};
1824 :     my $newGdir = $self->{_orgdir};
1825 :    
1826 :     my $sims_file = "$newGdir/expanded_similarities";
1827 :     my $sims_index_file = "$sims_file.index";
1828 :     my $sims_index = {};
1829 :    
1830 : olson 1.27 my $flips_file = "$newGdir/expanded_similarities.flips";
1831 :     my $flips_index_file = "$flips_file.index";
1832 :     my $flips_index = {};
1833 :    
1834 :     my $raw_sims_file = "$newGdir/similarities";
1835 :     my $raw_sims_index_file = "$raw_sims_file.index";
1836 :     my $raw_sims_index = {};
1837 :    
1838 :     my $raw_flips_file = "$newGdir/similarities.flips";
1839 :     my $raw_flips_index_file = "$raw_flips_file.index";
1840 :     my $raw_flips_index = {};
1841 :    
1842 :     $self->open_sims($sims_index, $sims_file, $sims_index_file, "_sims");
1843 :     $self->open_sims($flips_index, $flips_file, $flips_index_file, "_flips");
1844 :     $self->open_sims($raw_sims_index, $raw_sims_file, $raw_sims_index_file, "_raw_sims");
1845 :     $self->open_sims($raw_flips_index, $raw_flips_file, $raw_flips_index_file, "_raw_flips");
1846 :     }
1847 :    
1848 :     #
1849 :     # Open a sims file, tie it to a hash, and store the info into the $self obj.
1850 :     #
1851 :     sub open_sims
1852 :     {
1853 :     my($self, $hash, $sims_file, $index_file, $prefix) = @_;
1854 :    
1855 :     my $tied = tie %$hash, 'DB_File', $index_file, O_RDONLY, 0666, $DB_BTREE;
1856 : olson 1.13
1857 :     #
1858 :     # Set these even if failed so we don't keep trying to open and failing.
1859 :     #
1860 : olson 1.27 $self->{"${prefix}_index"} = $hash;
1861 :     $self->{"${prefix}_tie"} = $tied;
1862 : olson 1.13
1863 :     if (not $tied)
1864 :     {
1865 : olson 1.27 warn "Cannot tie sims index $index_file: $!\n";
1866 : olson 1.13 }
1867 :    
1868 :     #
1869 :     # open the sims file as well.
1870 :     #
1871 :    
1872 : olson 1.27 $self->{"${prefix}_fh"} = new FileHandle("<$sims_file");
1873 : olson 1.13
1874 : olson 1.27 if (!$self->{"${prefix}_fh"})
1875 : olson 1.13 {
1876 :     warn "Cannot open sims file $sims_file: $!\n";
1877 :     }
1878 :     }
1879 :    
1880 : overbeek 1.1 sub load_sims {
1881 :     my($self) = @_;
1882 :    
1883 :     if ($self->{_sims}) { return };
1884 :    
1885 :     my $newGdir = $self->{_orgdir};
1886 :    
1887 :     my $sims = {};
1888 :     foreach my $x (`cat $newGdir/similarities`)
1889 :     {
1890 :     chop; $x;
1891 :     if ($x =~ /^(\S+)/)
1892 :     {
1893 :     push(@{$sims->{$1}},bless([split(/\t/,$x)],'Sim'));
1894 :     }
1895 :     }
1896 :     $self->{_sims} = $sims;
1897 :     }
1898 :    
1899 :     sub get_attributes {
1900 :     my($self,$peg) = @_;
1901 :    
1902 :     my $fig = $self->{_fig};
1903 :     my $newG = $self->{_genome};
1904 :     my $newGdir = $self->{_orgdir};
1905 :    
1906 :     if (($peg =~ /^fig\|(\d+\.\d+)\.peg\.\d+$/) && ($1 eq $newG))
1907 :     {
1908 :     &load_attr($self);
1909 :     if (my $x = $self->{_attr}->{$peg})
1910 :     {
1911 :     return @$x;
1912 :     }
1913 :     else
1914 :     {
1915 :     return ();
1916 :     }
1917 :     }
1918 :     else
1919 :     {
1920 :     return $fig->get_attributes($peg);
1921 :     }
1922 :     }
1923 :    
1924 :     sub load_attr {
1925 :     my($self) = @_;
1926 :    
1927 :     if ($self->{_attr}) { return };
1928 :    
1929 :     my $newGdir = $self->{_orgdir};
1930 :     my $attr = {};
1931 : overbeek 1.3 foreach my $x (`cat $newGdir/evidence.codes`)
1932 : overbeek 1.1 {
1933 :     if ($x =~ /^(\S+)\t(\S+)/)
1934 :     {
1935 :     push(@{$attr->{$1}},[$1,"evidence_code",$2,""]);
1936 :     }
1937 :     }
1938 :     $self->{_attr} = $attr;
1939 :     }
1940 :    
1941 :     sub load_ann {
1942 :     my($self) = @_;
1943 :    
1944 :     if ($self->{_ann}) { return };
1945 :    
1946 :     my $newGdir = $self->{_orgdir};
1947 :     my $ann = {};
1948 :     if (open(ANN,"<$newGdir/annotations"))
1949 :     {
1950 :     $/ = "\n//\n";
1951 :     while (defined(my $x = <ANN>))
1952 :     {
1953 :     chomp $x;
1954 :     if ($x =~ /^(\S+)\n([^\n]+)\n([^\n]+)\n(.*)/s)
1955 :     {
1956 :     push(@{$ann->{$1}},[$1,$2,$3,"$4\n"]);
1957 :     }
1958 :     }
1959 :     $/ = "\n";
1960 :     close(ANN);
1961 :     }
1962 :     $self->{_ann} = $ann;
1963 :     }
1964 : overbeek 1.6
1965 :     sub taxonomy_of {
1966 :     my($self,$genome) = @_;
1967 :     my $fig = $self->{_fig};
1968 :     my $newG = $self->{_genome};
1969 :     my $newGdir = $self->{_orgdir};
1970 :    
1971 :     if ($newG ne $genome)
1972 :     {
1973 :     return $fig->taxonomy_of($genome);
1974 :     }
1975 :    
1976 :     my $tax;
1977 :     if (open(TAX,"<$newGdir/TAXONOMY") && ($tax = <TAX>))
1978 :     {
1979 :     chop $tax;
1980 :     return $tax;
1981 :     }
1982 :     else
1983 :     {
1984 :     return "unknown";
1985 :     }
1986 :     }
1987 :    
1988 :     sub build_tree_of_complete {
1989 :     my($self,$min_for_label) = @_;
1990 :     return $self->build_tree_of_all($min_for_label, "complete");
1991 :     }
1992 :    
1993 :     sub build_tree_of_all {
1994 :     my($self, $min_for_label, $complete)=@_;
1995 :     my(@last,@tax,$i,$prefix,$lev,$genome,$tax);
1996 :    
1997 :     $min_for_label = $min_for_label ? $min_for_label : 10;
1998 :     open(TMP,">/tmp/tree$$") || die "could not open /tmp/tree$$";
1999 :     print TMP "1. root\n";
2000 :    
2001 :     @last = ();
2002 :    
2003 :    
2004 :     foreach $genome (grep { ! $self->is_environmental($_) } $self->sort_genomes_by_taxonomy($self->genomes($complete)))
2005 :     {
2006 :     $tax = $self->taxonomy_of($genome);
2007 :     @tax = split(/\s*;\s*/,$tax);
2008 :     push(@tax,$genome);
2009 :     for ($i=0; ((@last > $i) && (@tax > $i) && ($last[$i] eq $tax[$i])); $i++) {}
2010 :     while ($i < @tax)
2011 :     {
2012 :     $lev = $i+2;
2013 :     $prefix = " " x (4 * ($lev-1));
2014 :     print TMP "$prefix$lev\. $tax[$i]\n";
2015 :     $i++;
2016 :     }
2017 :     @last = @tax;
2018 :     }
2019 :     close(TMP);
2020 :     my $tree = &tree_utilities::build_tree_from_outline("/tmp/tree$$");
2021 :     $tree->[0] = 'All';
2022 :     &FIG::limit_labels($tree,$min_for_label);
2023 :     unlink("/tmp/tree$$");
2024 :     return ($tree,&tree_utilities::tips_of_tree($tree));
2025 :     }
2026 :    
2027 :     sub sort_genomes_by_taxonomy {
2028 :     my($self,@genomes) = @_;
2029 :    
2030 :     return map { $_->[0] }
2031 :     sort { $a->[1] cmp $b->[1] }
2032 :     map { [$_,$self->taxonomy_of($_)] }
2033 :     @genomes;
2034 :     }
2035 : overbeek 1.1
2036 : overbeek 1.6 sub taxonomic_groups_of_complete {
2037 :     my($self,$min_for_labels) = @_;
2038 :    
2039 :     my($tree,undef) = $self->build_tree_of_complete($min_for_labels);
2040 :     return &FIG::taxonomic_groups($tree);
2041 :     }
2042 :    
2043 : olson 1.31 =head2 Search Database
2044 :    
2045 :     Searches the database for objects that match the query string in some way.
2046 :    
2047 :     Returns a list of results if the query is ambiguous or an unique identifier
2048 :     otherwise.
2049 :    
2050 :     =cut
2051 :    
2052 :     sub search_database {
2053 :     # get parameters
2054 :     my ($self, $query, $options) = @_;
2055 :    
2056 :     my $fig = $self->{_fig};
2057 :     my $newG = $self->{_genome};
2058 :     my $newGdir = $self->{_orgdir};
2059 :    
2060 :     #
2061 :     # Check for local ids and genome.
2062 :     #
2063 :    
2064 :    
2065 :     if ($query eq $newG or lc($query) eq lc($self->genus_species($newG)))
2066 :     {
2067 :     return { type => 'organism', result => $newG };
2068 :     }
2069 :    
2070 :     if ($query =~ /^fig\|(\d+\.\d+)\.peg/ and $1 eq $newG)
2071 :     {
2072 :     return { type => 'feature', result => $query };
2073 :     }
2074 :    
2075 :     #
2076 :     # Match on functions
2077 :     #
2078 :    
2079 :     &load_functions($self);
2080 :    
2081 :     my @protlist;
2082 :     my $fn = $self->{_functions};
2083 :     for my $id (keys %$fn)
2084 :     {
2085 :     if ($fn->{$id} =~ /$query/i)
2086 :     {
2087 :     push @protlist, [$id, $fn->{$id}, $newG];
2088 :     }
2089 :     }
2090 :    
2091 :     #
2092 :     # Pull FIG lookups.
2093 :     #
2094 :    
2095 :     my $res = $fig->search_database($query, $options);
2096 :    
2097 :     #
2098 :     # And integrate our protein list
2099 :     #
2100 :    
2101 :     my $res_prot;
2102 :     if (ref($res) eq 'ARRAY')
2103 :     {
2104 :     for my $ent (@$res)
2105 :     {
2106 :     if ($ent->{type} eq 'proteins')
2107 :     {
2108 :     $res_prot = $ent->{result};
2109 :     }
2110 :     }
2111 :     if (!$res_prot)
2112 :     {
2113 :     $res_prot = {type => 'proteins', result => \@protlist}
2114 :     }
2115 :     else
2116 :     {
2117 :     push(@$res_prot, @protlist);
2118 :     }
2119 :     }
2120 :    
2121 :     return $res;
2122 :    
2123 :     }
2124 :    
2125 :    
2126 : paczian 1.29 =head3 model_directory
2127 : paczian 1.26
2128 : paczian 1.29 C<< FIG->model_directory($organism); >>
2129 : paczian 1.26
2130 : paczian 1.29 Returns the model directory of an organism.
2131 : paczian 1.26
2132 :     =over 4
2133 :    
2134 :     =item $organism
2135 :    
2136 :     The seed-taxonomy id of the organism, e.g. 83333.1. If you pass the
2137 :     string 'All', you will get the directory for the global model.
2138 :    
2139 :     =back
2140 :    
2141 :     =cut
2142 :    
2143 : paczian 1.29 sub model_directory {
2144 : paczian 1.26 my ($self, $organism) = @_;
2145 :    
2146 : olson 1.27 my $directory;
2147 : paczian 1.26
2148 : olson 1.27 if ($organism eq $self->{_genome})
2149 :     {
2150 :     $directory = $self->{_orgdir} . "/Models";
2151 :     $directory .= "/$organism" if defined($organism);
2152 : formsma 1.38 }
2153 :     elsif(!defined $organism && defined $self->{_genome})
2154 :     {
2155 :     $directory = $self->{_orgdir} . "/Models";
2156 :     }
2157 :     else {
2158 : paczian 1.29 $directory = $self->{_fig}->model_directory($organism);
2159 : paczian 1.26 }
2160 :    
2161 :     return $directory;
2162 :     }
2163 :    
2164 : overbeek 1.1 1;

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3