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

Annotation of /FigKernelPackages/FIGV.pm

Parent Directory Parent Directory | Revision Log Revision Log


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

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

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3