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

Annotation of /FigKernelPackages/FIGV.pm

Parent Directory Parent Directory | Revision Log Revision Log


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

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

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3