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

Annotation of /FigKernelPackages/FIGV.pm

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1.59 - (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 : 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 :     # Just to make this a little less fragile, convert user of form
1260 :     # master:name to master. -- GJO
1261 :    
1262 :     $user = 'master'; # from this point on all assignments are treated as master assignments
1263 :    
1264 :     $function =~ s/\s+/ /sg;
1265 :     $function =~ s/^\s+//;
1266 :     $function =~ s/\s+$//;
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 :     if ( ! open( TMP, ">>$file" ) )
1307 :     {
1308 :     print STDERR "FAILED ASSIGNMENT: $fid\t$function\t$confidence\n";
1309 :     return 0;
1310 :     }
1311 :    
1312 :     flock(TMP,LOCK_EX) || confess "cannot lock assigned_functions";
1313 :     seek(TMP,0,2) || confess "failed to seek to the end of the file";
1314 :     print TMP "$fid\t$function\t$confidence\n";
1315 :     close(TMP);
1316 :     chmod(0777,$file);
1317 :    
1318 :     return 1;
1319 :     }
1320 :    
1321 :     sub add_annotation {
1322 :     my($self,$feature_id,$user,$annotation, $time_made) = @_;
1323 :     my($genome);
1324 :    
1325 :     if (!$self->is_virtual_feature($feature_id))
1326 :     {
1327 :     return $self->{_fig}->add_annotation($feature_id,$user,$annotation, $time_made);
1328 :     }
1329 :    
1330 :     $time_made = time unless $time_made =~ /^\d+$/;
1331 :    
1332 :     if ($self->is_deleted_fid($feature_id)) { return 0 }
1333 :    
1334 :     # print STDERR "add: fid=$feature_id user=$user annotation=$annotation\n";
1335 :     if ($genome = $self->genome_of($feature_id))
1336 :     {
1337 :     my $file = "$self->{_orgdir}/annotations";
1338 :     my $ma = ($annotation =~ /^Set master function to/);
1339 :    
1340 :     if (open(TMP,">>$file"))
1341 :     {
1342 :     flock(TMP,LOCK_EX) || confess "cannot lock annotations";
1343 :     seek(TMP,0,2) || confess "failed to seek to the end of the file";
1344 :    
1345 :     my $dataLine = "$feature_id\n$time_made\n$user\n$annotation" . ((substr($annotation,-1) eq "\n") ? "" : "\n");
1346 :     print TMP $dataLine . "//\n";
1347 :     close(TMP);
1348 :     chmod 0777, $file;
1349 :    
1350 :     #
1351 :     # Update local cache.
1352 :     #
1353 :     my $ann = $self->{_ann};
1354 :     push(@{$ann->{$feature_id}}, [$feature_id, $time_made, $user, $annotation . "\n"]);
1355 :     }
1356 :     }
1357 :     return 0;
1358 :     }
1359 :    
1360 :    
1361 : overbeek 1.1 sub function_of {
1362 :     my($self,$fid) = @_;
1363 :    
1364 :     my $fig = $self->{_fig};
1365 :     my $newG = $self->{_genome};
1366 :     my $newGdir = $self->{_orgdir};
1367 :    
1368 :     if (($fid =~ /^fig\|(\d+\.\d+)/) && ($1 eq $newG))
1369 :     {
1370 :     &load_functions($self);
1371 : olson 1.56
1372 :     my $fn = $self->{_functions}->{$fid};
1373 :     if (wantarray)
1374 :     {
1375 :     return ['master', $fn];
1376 :     }
1377 :     else
1378 :     {
1379 :     return $fn;
1380 :     }
1381 : overbeek 1.1 }
1382 :     else
1383 :     {
1384 : olson 1.56 return $fig->function_of($fid);
1385 : overbeek 1.1 }
1386 :     }
1387 :    
1388 : redwards 1.54
1389 :     =pod
1390 :    
1391 :     find_features_by_annotation
1392 :    
1393 :     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.
1394 :    
1395 :     If the case boolean is set the search is case insensitive. Otherwise case sensitive.
1396 :    
1397 :     =cut
1398 :    
1399 :     sub find_features_by_annotation {
1400 :     my($self,$anno_hash, $case)=@_;
1401 :     $self->load_functions;
1402 :    
1403 :     if ($case) {map {$anno_hash->{uc($_)}=1} keys %$anno_hash}
1404 :    
1405 :     my $res={};
1406 :     foreach my $id (keys %{$self->{_functions}})
1407 :     {
1408 :     my $fn = $self->{_functions}->{$id};
1409 :     $case ? $fn = uc($fn) : 1;
1410 :     if ($anno_hash->{$fn}) {push @{$res->{$fn}}, $id}
1411 :     }
1412 :    
1413 :     return $res;
1414 :     }
1415 :    
1416 :    
1417 :     =pod
1418 :    
1419 :     search_features_by_annotation
1420 :    
1421 :     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.
1422 :    
1423 :     If the case boolean is set the search is case insensitive. Otherwise case sensitive.
1424 :    
1425 :     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.
1426 :    
1427 :    
1428 :     =cut
1429 :    
1430 :     sub search_features_by_annotation {
1431 :     my($self,$term, $case)=@_;
1432 :     $self->load_functions;
1433 :    
1434 :     # to make case insensitive convert everything to uppercase
1435 :     # alternative is to use two regexps, one for case insens and one for not case insens
1436 :     # but the bad thing about that approach is that if you have a case insensitive search
1437 :     # you do two regexps for each failed match
1438 :    
1439 :     $case ? $term = uc($term) : 1;
1440 :    
1441 :     my $res={};
1442 :     foreach my $id (keys %{$self->{_functions}})
1443 :     {
1444 :     # we set two variables, one that has case changed for case insensitive searches
1445 :     my $fn = my $fnc = $self->{_functions}->{$id};
1446 :     $case ? $fn = uc($fn) : 1;
1447 :     if ($fn =~ m/$term/) {push @{$res->{$fnc}}, $id}
1448 :     }
1449 :    
1450 :     return $res;
1451 :     }
1452 :    
1453 :    
1454 : overbeek 1.1 sub feature_aliases {
1455 :     my($self,$fid) = @_;
1456 :    
1457 :     my $fig = $self->{_fig};
1458 :     my $newG = $self->{_genome};
1459 :     my $newGdir = $self->{_orgdir};
1460 :    
1461 :     my @aliases;
1462 :     if (($fid =~ /^fig\|(\d+\.\d+)/) && ($1 eq $newG))
1463 :     {
1464 :     &load_tbl($self);
1465 :     if (my $x = $self->{_tbl}->{$fid})
1466 :     {
1467 :     @aliases = @{$x->[1]};
1468 :     }
1469 :     else
1470 :     {
1471 :     @aliases = ();
1472 :     }
1473 :     }
1474 :     else
1475 :     {
1476 :     @aliases = $fig->feature_aliases($fid);
1477 :     }
1478 :     return wantarray() ? @aliases : join(",",@aliases);
1479 :     }
1480 :    
1481 :     sub feature_annotations {
1482 :     my($self,$fid,$rawtime) = @_;
1483 :    
1484 :     my $fig = $self->{_fig};
1485 :     my $newG = $self->{_genome};
1486 :     my $newGdir = $self->{_orgdir};
1487 :    
1488 :     my @annotations;
1489 :     if (($fid =~ /^fig\|(\d+\.\d+)/) && ($1 eq $newG))
1490 :     {
1491 :     &load_ann($self);
1492 :     if (my $x = $self->{_ann}->{$fid})
1493 :     {
1494 :     @annotations = @{$x};
1495 :     }
1496 :     else
1497 :     {
1498 :     @annotations = ();
1499 :     }
1500 :    
1501 :     if ($rawtime)
1502 :     {
1503 :     return @annotations;
1504 :     }
1505 :     else
1506 :     {
1507 :     return map { $_->[1] = localtime($_->[1]); $_ } @annotations;
1508 :     }
1509 :     }
1510 :     else
1511 :     {
1512 :     return $fig->feature_annotations($fid);
1513 :     }
1514 :     }
1515 :    
1516 :     sub get_translation {
1517 :     my($self,$peg) = @_;
1518 :    
1519 :     my $fig = $self->{_fig};
1520 :     my $newG = $self->{_genome};
1521 :     my $newGdir = $self->{_orgdir};
1522 :    
1523 :     if (($peg =~ /^fig\|(\d+\.\d+)/) && ($1 eq $newG))
1524 :     {
1525 : olson 1.55 &load_feature_indexes($self);
1526 :    
1527 :     my $out = $self->{_feat}->{peg}->{$peg};
1528 :     if (defined($out))
1529 :     {
1530 :     return $out;
1531 :     }
1532 :    
1533 :     #
1534 :     # If we have a blast-formatted fasta, use fastacmd to
1535 :     # do the lookup, and cache the output for later use.
1536 :     #
1537 :    
1538 :     if ($self->{_feat_fasta}->{peg})
1539 :     {
1540 :     my $id = "gnl|$peg";
1541 :     my $cmd = "$FIG_Config::ext_bin/fastacmd -d $self->{_feat_fasta}->{peg} -s '$id'";
1542 :     open(P, "$cmd|") or die "get_translation: cmd failed with $!: $cmd";
1543 :     $_ = <P>;
1544 :     my $out;
1545 :     while (<P>)
1546 :     {
1547 :     s/\s+//g;
1548 :     $out .= $_;
1549 :     }
1550 :     close(P);
1551 :     $self->{_feat}->{$peg} = $out;
1552 :     return $out;
1553 :     }
1554 :     else
1555 :     {
1556 :     return $self->{_feat}->{$peg};
1557 :     }
1558 : overbeek 1.1 }
1559 :     else
1560 :     {
1561 :     return $fig->get_translation($peg);
1562 :     }
1563 :     }
1564 :    
1565 : olson 1.4 sub translation_length
1566 :     {
1567 :     my($self, $peg) = @_;
1568 :    
1569 :     my $fig = $self->{_fig};
1570 :     my $newG = $self->{_genome};
1571 :     my $newGdir = $self->{_orgdir};
1572 :    
1573 :     if (($peg =~ /^fig\|(\d+\.\d+)/) && ($1 eq $newG))
1574 :     {
1575 : olson 1.55 my $t = $self->get_translation($peg);
1576 :     return length($t);
1577 : olson 1.4 }
1578 :     else
1579 :     {
1580 :     return $fig->translation_length($peg);
1581 :     }
1582 :     }
1583 :    
1584 :     sub translatable
1585 :     {
1586 :     my($self, $peg) = @_;
1587 :    
1588 :     my $fig = $self->{_fig};
1589 :     my $newG = $self->{_genome};
1590 :     my $newGdir = $self->{_orgdir};
1591 :    
1592 :     if (($peg =~ /^fig\|(\d+\.\d+)/) && ($1 eq $newG))
1593 :     {
1594 : olson 1.55 return $self->translation_length($peg) > 0 ? 1 : 0;
1595 : olson 1.4 }
1596 :     else
1597 :     {
1598 :     return $fig->translatable($peg);
1599 :     }
1600 :     }
1601 :    
1602 : gdpusch 1.48 sub pick_gene_boundaries {
1603 :     return &FIG::pick_gene_boundaries(@_);
1604 :     }
1605 :    
1606 :     sub call_start {
1607 :     return &FIG::call_start(@_);
1608 :     }
1609 :    
1610 : olson 1.4 sub is_real_feature
1611 :     {
1612 :     my($self, $fid) = @_;
1613 :    
1614 :     if ($self->is_virtual_feature($fid))
1615 :     {
1616 :     return 1;
1617 :     }
1618 :     else
1619 :     {
1620 :     return $self->{_fig}->is_real_feature($fid);
1621 :     }
1622 : parrello 1.49
1623 : olson 1.4 }
1624 :    
1625 : olson 1.56 sub is_deleted_fid
1626 :     {
1627 :     my($self, $fid) = @_;
1628 :    
1629 :     if ($self->is_virtual_feature($fid))
1630 :     {
1631 :     return 0;
1632 :     }
1633 :     else
1634 :     {
1635 :     return $self->{_fig}->is_real_feature($fid);
1636 :     }
1637 :    
1638 :     }
1639 :    
1640 : olson 1.10 sub pegs_of
1641 :     {
1642 :     my($self, $genome) = @_;
1643 :    
1644 :     my $fig = $self->{_fig};
1645 :     my $newG = $self->{_genome};
1646 :     my $newGdir = $self->{_orgdir};
1647 :    
1648 :     if ($genome ne $newG)
1649 :     {
1650 :     return $fig->pegs_of($genome);
1651 :     }
1652 :    
1653 : olson 1.55 $self->load_feature_indexes();
1654 :    
1655 :     my $t = $self->{_feat_btree}->{peg};
1656 :     if ($t)
1657 :     {
1658 :     return keys %$t;
1659 :     }
1660 :     else
1661 :     {
1662 :     warn "Reverting to load_tbl for $newG\n";
1663 :     $self->load_tbl();
1664 :     return grep { /\.peg\./ } keys %{$self->{_tbl}};
1665 :     }
1666 : olson 1.10 }
1667 :    
1668 :    
1669 :     sub rnas_of
1670 :     {
1671 :     my($self, $genome) = @_;
1672 :    
1673 :     my $fig = $self->{_fig};
1674 :     my $newG = $self->{_genome};
1675 :     my $newGdir = $self->{_orgdir};
1676 :    
1677 :     if ($genome ne $newG)
1678 :     {
1679 : olson 1.55 return $fig->rna_of($genome);
1680 :     }
1681 :    
1682 :     $self->load_feature_indexes();
1683 :    
1684 :     my $t = $self->{_feat_btree}->{rna};
1685 :     if ($t)
1686 :     {
1687 :     return keys %$t;
1688 :     }
1689 :     else
1690 :     {
1691 :     warn "Reverting to load_tbl for $newG\n";
1692 :     $self->load_tbl();
1693 :     return grep { /\.rna\./ } keys %{$self->{_tbl}};
1694 : olson 1.10 }
1695 :    
1696 :     }
1697 :    
1698 : olson 1.4 sub is_virtual_feature
1699 :     {
1700 :     my($self, $peg) = @_;
1701 :    
1702 :     my $fig = $self->{_fig};
1703 :     my $newG = $self->{_genome};
1704 :     my $newGdir = $self->{_orgdir};
1705 :    
1706 :     if (($peg =~ /^fig\|(\d+\.\d+)/) && ($1 eq $newG))
1707 :     {
1708 :     return 1;
1709 :     }
1710 :     else
1711 :     {
1712 :     return 0;
1713 :     }
1714 :     }
1715 :    
1716 : olson 1.35 ####################################################
1717 :     #
1718 :     # Following are some MG-RAST specific features. FIGV seems a good a place as any to include them.#
1719 :     #
1720 :     #
1721 :    
1722 :     =head3 taxa_to_seed_ids
1723 :    
1724 :     Given a prefix of a taxonomy string, return the list of metagenome fids that
1725 :     mapped to SEED organisms in that taxonomy.
1726 :    
1727 :     =cut
1728 :    
1729 :     sub taxa_to_seed_ids
1730 :     {
1731 :     my($self, $tax) = @_;
1732 :    
1733 :     my $btree_file = "$self->{_orgdir}/taxa_summary_by_blast.btree";
1734 :     my %btree;
1735 :    
1736 :     my $btree_tie = tie %btree, 'DB_File', $btree_file, O_RDONLY, 0666, $DB_BTREE;
1737 :    
1738 :     if (!$btree_tie)
1739 :     {
1740 :     warn "Cannot tie $btree_file: $!\n";
1741 :     return undef;
1742 :     }
1743 :    
1744 :     my $key = $tax;
1745 :     my ($res, $value);
1746 :     my @vals;
1747 :     my %seen;
1748 :    
1749 :     for ($res = $btree_tie->seq($key, $value, R_CURSOR);
1750 :     $res == 0 and $key =~ /^$tax/;
1751 :     $res = $btree_tie->seq($key, $value, R_NEXT))
1752 :     {
1753 :     print "$key\n";
1754 :     push(@vals, map { [$key, $_] } grep { not $seen{$_}++ } split(/$;/, $value));
1755 :     }
1756 :     return @vals;
1757 :     }
1758 :    
1759 : olson 1.55 sub load_feature_indexes
1760 :     {
1761 :     my($self) = @_;
1762 :    
1763 :     if ($self->{_feat}) { return };
1764 :    
1765 :     my $newGdir = $self->{_orgdir};
1766 :    
1767 :     for my $fdir (<$newGdir/Features/*>)
1768 :     {
1769 :     my $ftype = basename($fdir);
1770 :    
1771 :     #
1772 :     # If we have a tbl.btree, tie that for our use.
1773 :     #
1774 :    
1775 :     my $tbl_idx = {};
1776 :     my $tie = tie %$tbl_idx, 'DB_File', "$fdir/tbl.btree", O_RDONLY, 0666, $DB_BTREE;
1777 :     if ($tie)
1778 :     {
1779 :     $self->{_feat_tie}->{$ftype} = $tie;
1780 :     $self->{_feat_btree}->{$ftype} = $tbl_idx;
1781 :     }
1782 :    
1783 :     my $tbl_list = [];
1784 :     my $ltie = tie @$tbl_list, 'DB_File', "$fdir/tbl.recno", O_RDONLY, 0666, $DB_RECNO;
1785 :     if ($tie)
1786 :     {
1787 :     $self->{_feat_ltie}->{$ftype} = $ltie;
1788 :     $self->{_feat_recno}->{$ftype} = $tbl_list;
1789 :     }
1790 :    
1791 :     #
1792 :     # If we have fasta.norm.phr, set _pseq_fasta to the fasta file to use with fastacmd.
1793 :     #
1794 :    
1795 :     my $pseq = {};
1796 :    
1797 :     if (-f "$fdir/fasta.norm.phr")
1798 :     {
1799 :     $self->{_feat_fasta}->{$ftype} = "$fdir/fasta.norm";
1800 :     }
1801 :     else
1802 :     {
1803 :     #
1804 :     # Otherwise, we need to load the data.
1805 :     #
1806 :    
1807 :     if (open(FASTA,"<$newGdir/Features/peg/fasta"))
1808 :     {
1809 :     local $/ = "\n>";
1810 :     my $x;
1811 :     while (defined($x = <FASTA>))
1812 :     {
1813 :     chomp $x;
1814 :     if ($x =~ />?(\S+)[^\n]*\n(.*)/s)
1815 :     {
1816 :     my $peg = $1;
1817 :     my $seq = $2;
1818 :     $seq =~ s/\s//gs;
1819 :     $pseq->{$peg} = $seq;
1820 :     }
1821 :     }
1822 :     close(FASTA);
1823 :     }
1824 :     }
1825 :     $self->{_feat}->{$ftype} = $pseq;
1826 :     }
1827 :     }
1828 :    
1829 : overbeek 1.1 sub load_pseq {
1830 :     my($self) = @_;
1831 :    
1832 :     if ($self->{_pseq}) { return };
1833 :    
1834 :     my $newGdir = $self->{_orgdir};
1835 : olson 1.55 my $fdir = "$newGdir/Features/peg";
1836 :    
1837 :     #
1838 :     # If we have a tbl.btree, tie that for our use.
1839 :     #
1840 :    
1841 :     my $tbl_idx = {};
1842 :     my $tie = tie %$tbl_idx, 'DB_File', "$fdir/tbl.btree", O_RDONLY, 0666, $DB_BTREE;
1843 :     if ($tie)
1844 :     {
1845 :     $self->{_tbl_tie} = $tie;
1846 :     $self->{_tbl_btree} = $tbl_idx;
1847 :     }
1848 :    
1849 :     #
1850 :     # If we have fasta.norm.phr, set _pseq_fasta to the fasta file to use with fastacmd.
1851 :     #
1852 :    
1853 : overbeek 1.1 my $pseq = {};
1854 : olson 1.55
1855 :     if (-f "$fdir/fasta.norm.phr")
1856 : overbeek 1.1 {
1857 : olson 1.55 $self->{_pseq_fasta} = "$fdir/fasta.norm";
1858 :     }
1859 :     else
1860 :     {
1861 :     #
1862 :     # Otherwise, we need to load the data.
1863 :     #
1864 :    
1865 :     if (open(FASTA,"<$newGdir/Features/peg/fasta"))
1866 : overbeek 1.1 {
1867 : olson 1.55 local $/ = "\n>";
1868 :     my $x;
1869 :     while (defined($x = <FASTA>))
1870 : overbeek 1.1 {
1871 : olson 1.55 chomp $x;
1872 :     if ($x =~ />?(\S+)[^\n]*\n(.*)/s)
1873 :     {
1874 :     my $peg = $1;
1875 :     my $seq = $2;
1876 :     $seq =~ s/\s//gs;
1877 :     $pseq->{$peg} = $seq;
1878 :     }
1879 : overbeek 1.1 }
1880 : olson 1.55 close(FASTA);
1881 : overbeek 1.1 }
1882 :     }
1883 :     $self->{_pseq} = $pseq;
1884 :     }
1885 :    
1886 : olson 1.9 sub load_ss_data
1887 :     {
1888 :     my($self) = @_;
1889 :    
1890 :     return if defined($self->{_ss_bindings});
1891 :    
1892 :     my $fig = $self->{_fig};
1893 :     my $newG = $self->{_genome};
1894 :     my $newGdir = $self->{_orgdir};
1895 :    
1896 :     open(S, "<$newGdir/Subsystems/bindings") or die "Cannot open $newGdir/Subsystems/bindings: $!";
1897 :    
1898 :     my $peg_index;
1899 :     my $bindings;
1900 : olson 1.14 my $genome_index;
1901 :     my %seen;
1902 : olson 1.9 while (<S>)
1903 :     {
1904 :     chomp;
1905 :     my($sname, $role, $peg) = split(/\t/);
1906 : olson 1.14
1907 : olson 1.34 my($genome) = ($peg =~ /fig\|(\d+\.\d+)/);
1908 :    
1909 :     # my $genome = &FIG::genome_of($peg);
1910 : parrello 1.49
1911 : olson 1.9 push(@{$bindings->{$sname}->{$role}}, $peg);
1912 :     push(@{$peg_index->{$peg}}, [$sname, $role]);
1913 : olson 1.14
1914 :     if (!$seen{$genome, $sname, $role})
1915 :     {
1916 : parrello 1.24 push(@{$genome_index->{$genome}}, [$sname, $role]);
1917 : olson 1.14 }
1918 : parrello 1.49
1919 : olson 1.9 }
1920 :     close(S);
1921 :    
1922 :     open(S, "<$newGdir/Subsystems/subsystems") or die "Cannot open $newGdir/Subsystems/subsystems: $!";
1923 :     my $variant;
1924 :     while (<S>)
1925 :     {
1926 :     chomp;
1927 :     my($sname, $v) = split(/\t/);
1928 :     $variant->{$sname} = $v;
1929 :     }
1930 :     close(S);
1931 :    
1932 : gdpusch 1.59 $self->{_ss_bindings} = $bindings;
1933 :     $self->{_ss_variants} = $variant;
1934 :     $self->{_ss_peg_index} = $peg_index;
1935 : olson 1.14 $self->{_ss_genome_index} = $genome_index;
1936 : olson 1.9 }
1937 :    
1938 : overbeek 1.1 sub load_tbl {
1939 :     my($self) = @_;
1940 :    
1941 :     if ($self->{_tbl}) { return };
1942 :    
1943 :     my $newGdir = $self->{_orgdir};
1944 :     my $tbl = {};
1945 :     foreach my $x (`cat $newGdir/Features/*/tbl`)
1946 :     {
1947 : gdpusch 1.59 chomp $x;
1948 : overbeek 1.1 if ($x =~ /^(\S+)\t(\S+)(\t(\S.*\S))?/)
1949 :     {
1950 :     my $fid = $1;
1951 :     my $loc = [split(/,/,$2)];
1952 : paczian 1.12 my $aliases = $4 ? [split(/\t/,$4)] : [];
1953 : overbeek 1.1 $tbl->{$fid} = [$loc,$aliases];
1954 :     }
1955 : gdpusch 1.59 else {
1956 :     warn "Bad feature line in $newGdir: $x";
1957 :     }
1958 : overbeek 1.1 }
1959 : gdpusch 1.59 print STDERR ("Loaded ", (scalar keys %$tbl), " features from $newGdir\n") if $ENV{FIG_VERBOSE};
1960 :    
1961 : overbeek 1.1 $self->{_tbl} = $tbl;
1962 :     }
1963 :    
1964 :     sub load_functions {
1965 :     my($self) = @_;
1966 :    
1967 :     if ($self->{_functions}) { return };
1968 :    
1969 :     my $newG = $self->{_genome};
1970 :     my $newGdir = $self->{_orgdir};
1971 :    
1972 :     my $functions = {};
1973 : olson 1.55
1974 :     my $bfile = "$newGdir/assigned_functions.btree";
1975 :     my $tie;
1976 :     if (-f $bfile)
1977 :     {
1978 :     $tie = tie %$functions, 'DB_File', $bfile, O_RDONLY, 0666, $DB_BTREE;
1979 :     }
1980 :    
1981 :     if (!$tie)
1982 : overbeek 1.1 {
1983 : olson 1.55 foreach my $x (`cat $newGdir/*functions`)
1984 : overbeek 1.1 {
1985 : olson 1.55 if (($x =~ /^(fig\|(\d+\.\d+)\.\S+)\t(\S[^\t]*\S)/) && ($2 eq $newG))
1986 :     {
1987 :     $functions->{$1} = $3;
1988 :     }
1989 : overbeek 1.1 }
1990 :     }
1991 : olson 1.55
1992 : overbeek 1.1 $self->{_functions} = $functions;
1993 :     }
1994 :    
1995 : olson 1.16
1996 :     sub bbhs
1997 :     {
1998 :     my($self,$peg,$cutoff) = @_;
1999 :    
2000 :     my $fig = $self->{_fig};
2001 :     my $newG = $self->{_genome};
2002 :     my $newGdir = $self->{_orgdir};
2003 :    
2004 :     $self->load_bbh_index();
2005 :    
2006 :     $cutoff = 1.0e-10 unless defined($cutoff);
2007 :    
2008 :     my $tie = $self->{_bbh_tie};
2009 :    
2010 :     my @bbhs = $tie->get_dup($peg);
2011 :    
2012 :     # @bbhs now a list of comma-separated tuples (id2, score, bitscore)
2013 :     # print Dumper(\@bbhs);
2014 :    
2015 :     @bbhs = grep { $_->[1] < $cutoff } map { [split(/,/)] } @bbhs;
2016 :    
2017 :     if (not ($peg =~ /^fig\|(\d+\.\d+)/ && ($1 eq $newG)))
2018 :     {
2019 :     #
2020 :     # If this isn't one of our pegs, we retrieve the bbhs from the underlying
2021 :     # SEED and merge in the bbhs that we have here.
2022 :     #
2023 :    
2024 :     my @fbbhs = $fig->bbhs($peg, $cutoff);
2025 :     # print Dumper(\@fbbhs);
2026 :     push(@bbhs, @fbbhs);
2027 :     }
2028 :     return sort { $a->[1] <=> $b->[1] } @bbhs;
2029 : parrello 1.49 }
2030 : olson 1.16
2031 : olson 1.13 sub sims
2032 :     {
2033 : olson 1.27 my($self,$pegarg,$max,$maxP,$select, $max_expand, $filters) = @_;
2034 : olson 1.13
2035 :     my $fig = $self->{_fig};
2036 :     my $newG = $self->{_genome};
2037 :     my $newGdir = $self->{_orgdir};
2038 : parrello 1.24 $max = $max ? $max : 10000;
2039 :     $maxP = $maxP ? $maxP : 1.0e-5;
2040 : olson 1.13 $select = $select ? $select : "all";
2041 :    
2042 : olson 1.27 my $prefix = "_sims";
2043 :     my $flip_prefix = "_flips";
2044 :    
2045 :     if ($select eq 'raw')
2046 :     {
2047 :     $prefix = "_raw_sims";
2048 :     $flip_prefix = "_raw_flips";
2049 :     }
2050 :    
2051 :     #
2052 :     # Partition pegs into one set that is part of this organism
2053 :     # and another set that is not.
2054 :     #
2055 :     my @pegs = ref($pegarg) ? @$pegarg : ($pegarg);
2056 :     my(@part, @not_part, %part);
2057 :     my %sims;
2058 : parrello 1.49
2059 : olson 1.27 for my $peg (@pegs)
2060 :     {
2061 :     if ($peg =~ /^fig\|(\d+\.\d+)/ and $1 eq $newG)
2062 :     {
2063 :     push(@part, $peg);
2064 :     $part{$peg}++;
2065 :     }
2066 :     else
2067 :     {
2068 :     push(@not_part, $peg);
2069 :     $sims{$peg} = [];
2070 :     }
2071 :     }
2072 :    
2073 :     #
2074 : olson 1.58 # Retrieve a batch of the sims from the normal SEED sims for the not-part pegs, and partition
2075 : olson 1.27 # into %sims.
2076 :     #
2077 : olson 1.58
2078 :     if (@not_part)
2079 : olson 1.27 {
2080 : olson 1.58 for my $sim ($fig->sims(\@not_part, $max, $maxP, $select, $max_expand, $filters))
2081 :     {
2082 :     push(@{$sims{$sim->id1}}, $sim);
2083 :     }
2084 : olson 1.27 }
2085 :    
2086 :     my @out;
2087 :     my $start_len;
2088 :    
2089 :     for my $peg (@pegs)
2090 : olson 1.13 {
2091 : olson 1.27 $start_len = @out;
2092 : olson 1.58
2093 : olson 1.27 if (not $part{$peg})
2094 : overbeek 1.23 {
2095 : olson 1.58 #
2096 :     # For the ids that are not part of this organism, pull the flips that we
2097 :     # computed - these are the sims between this org and the one in the requested
2098 :     # peg.
2099 :     #
2100 : olson 1.27 my @flips = $self->retrieve_sims($peg, $flip_prefix, $maxP, $select);
2101 : parrello 1.49
2102 : olson 1.27 @flips = sort { $b->bsc <=> $a->bsc } @flips;
2103 : parrello 1.49
2104 : olson 1.27 # my @old = $fig->sims($peg,$max,$maxP,$select);
2105 : parrello 1.49
2106 : olson 1.58 #
2107 :     # Merge these sims together with any that we retrieved
2108 :     # earlier - recall sims are returned in bsc order.
2109 :     #
2110 :    
2111 : olson 1.27 my @old = sort { $b->bsc <=> $a->bsc } @{$sims{$peg}};
2112 :    
2113 :     # my @merged = ();
2114 :     my $i1 = 0;
2115 :     my $i2 = 0;
2116 :     while (($i1 < @flips) || ($i2 < @old))
2117 : overbeek 1.23 {
2118 : olson 1.27 if (($i1 == @flips) || (($i2 < @old) && ($flips[$i1]->[11] < $old[$i2]->[11])))
2119 :     {
2120 :     # push(@merged,$old[$i2]);
2121 :     push(@out,$old[$i2]);
2122 :     $i2++;
2123 :     }
2124 :     else
2125 :     {
2126 :     # push(@merged,$flips[$i1]);
2127 :     push(@out,$flips[$i1]);
2128 :     $i1++;
2129 :     }
2130 : overbeek 1.23 }
2131 :     }
2132 : olson 1.27 else
2133 :     {
2134 : olson 1.58 #
2135 :     # For the ids that are in this organism, we just pull the sims directly.
2136 :     #
2137 : olson 1.27 my @sims = $self->retrieve_sims($peg, $prefix, $maxP, $select);
2138 :     push(@out, @sims);
2139 :     }
2140 :    
2141 :     if (@out - $start_len > $max)
2142 :     {
2143 :     $#out = $start_len + $max - 1;
2144 :     }
2145 : olson 1.13 }
2146 : parrello 1.49
2147 : olson 1.27 return @out;
2148 :     }
2149 : olson 1.13
2150 : olson 1.27 sub retrieve_sims
2151 :     {
2152 :     my($self, $peg, $prefix, $maxP, $select) = @_;
2153 : parrello 1.49
2154 : olson 1.13 $self->load_sims_index();
2155 :    
2156 : olson 1.27 my $info = $self->{"${prefix}_index"}->{$peg};
2157 : olson 1.13
2158 :     if ($info !~ /^(\d+),(\d+)$/)
2159 :     {
2160 :     return;
2161 :     }
2162 :     my($seek, $len) = ($1, $2);
2163 :    
2164 : olson 1.27 my $sims_txt = &FIG::read_block($self->{"${prefix}_fh"}, $seek, $len);
2165 : olson 1.13
2166 :     my @sims = map { bless $_, 'Sim' } grep { $_->[10] <= $maxP } map { [split(/\t/)] } @$sims_txt;
2167 :    
2168 :     if ($select =~ /^figx?$/)
2169 :     {
2170 :     @sims = grep { $_->[1] =~ /^fig/ } @sims;
2171 :     }
2172 :    
2173 :     return @sims;
2174 :     }
2175 :    
2176 :     sub sims_old {
2177 : overbeek 1.1 my($self,$peg,$max,$maxP,$select) = @_;
2178 :    
2179 :     my $fig = $self->{_fig};
2180 :     my $newG = $self->{_genome};
2181 :     my $newGdir = $self->{_orgdir};
2182 : parrello 1.24 $max = $max ? $max : 10000;
2183 :     $maxP = $maxP ? $maxP : 1.0e-5;
2184 : overbeek 1.1 $select = $select ? $select : "all";
2185 :    
2186 :     if (($peg =~ /^fig\|(\d+\.\d+)/) && ($1 eq $newG))
2187 :     {
2188 :     &load_sims($self);
2189 :     my @sims = ();
2190 :     my $raw_sims = $self->{_sims}->{$peg};
2191 :     if ($raw_sims)
2192 :     {
2193 :     foreach my $sim ( grep { $_->[10] <= $maxP } @$raw_sims )
2194 :     {
2195 :     my $id2 = $sim->id2;
2196 :     my $id1 = $sim->id1;
2197 :     my @relevant = ();
2198 :    
2199 :     my @maps_to = $fig->mapped_prot_ids( $id2 );
2200 :     my $ref_len = $maps_to[0]->[1];
2201 :    
2202 :     @maps_to = grep { $_->[0] !~ /^xxx\d+/ } @maps_to;
2203 : parrello 1.49
2204 : overbeek 1.1 if ( $select =~ /^figx?$/ ) # Only fig
2205 :     {
2206 :     @relevant = grep { $_->[0] =~ /^fig/ } @maps_to;
2207 :     }
2208 :     else # All
2209 :     {
2210 :     @relevant = @maps_to;
2211 :     }
2212 :    
2213 :     my $seen = {};
2214 :     foreach my $x ( @relevant )
2215 :     {
2216 :     my ( $x_id, $x_ln ) = @$x;
2217 : parrello 1.49
2218 : overbeek 1.1 next if $seen->{$x_id};
2219 :     $seen->{$x_id} = 1;
2220 :    
2221 :     my $delta2 = $ref_len - $x_ln; # Coordinate shift
2222 :     my $sim1 = [ @$sim ]; # Make a copy
2223 :     $sim1->[1] = $x_id;
2224 :     $sim1->[8] -= $delta2;
2225 :     $sim1->[9] -= $delta2;
2226 :     bless( $sim1, "Sim" );
2227 :     push( @sims, $sim1 );
2228 :     }
2229 :     }
2230 :     }
2231 :    
2232 :     if (@sims > $max) { $#sims = $max-1; }
2233 :     return @sims;
2234 :     }
2235 :     else
2236 :     {
2237 :     return $fig->sims($peg,$max,$maxP,$select);
2238 :     }
2239 :     }
2240 :    
2241 : olson 1.27
2242 :     sub coupled_to
2243 :     {
2244 :     my($self,$peg) = @_;
2245 :    
2246 :     my $fig = $self->{_fig};
2247 :     my $newG = $self->{_genome};
2248 :     my $newGdir = $self->{_orgdir};
2249 :    
2250 :     if ($peg =~ /^fig\|$newG\.peg/)
2251 :     {
2252 : parrello 1.49
2253 : olson 1.27 $self->load_coupling_index();
2254 : parrello 1.49
2255 : olson 1.27 my $tie = $self->{_pch_tie};
2256 : olson 1.30
2257 :     if (not defined($tie))
2258 :     {
2259 :     return;
2260 :     }
2261 : parrello 1.49
2262 : olson 1.27 my @dat = $tie->get_dup($peg);
2263 : parrello 1.49
2264 : olson 1.27 return map { [split(/$;/, $_)] } @dat;
2265 :     }
2266 :     else
2267 :     {
2268 :     return $fig->coupled_to($peg);
2269 :     }
2270 :    
2271 :     }
2272 :    
2273 :     sub coupling_evidence
2274 :     {
2275 :     my($self,$peg1, $peg2) = @_;
2276 :    
2277 :     my $fig = $self->{_fig};
2278 :     my $newG = $self->{_genome};
2279 :     my $newGdir = $self->{_orgdir};
2280 :    
2281 :     if ($peg1 =~ /^fig\|$newG\.peg/)
2282 :     {
2283 :     $self->load_coupling_index();
2284 : parrello 1.49
2285 : olson 1.27 my $tie = $self->{_pch_ev_tie};
2286 : parrello 1.49
2287 : olson 1.30 if (not defined($tie))
2288 :     {
2289 :     return;
2290 :     }
2291 :    
2292 : olson 1.27 my @dat = $tie->get_dup("$peg1$;$peg2");
2293 :    
2294 :     my @a;
2295 :     return map { @a = split(/$;/, $_); [@a[0,1,4]] } @dat;
2296 :     }
2297 :     else
2298 :     {
2299 :     return $fig->coupling_evidence($peg1, $peg2);
2300 :     }
2301 :    
2302 :     }
2303 :    
2304 :     sub coupling_and_evidence
2305 :     {
2306 :     my($self,$peg1) = @_;
2307 :    
2308 :     my $fig = $self->{_fig};
2309 :     my $newG = $self->{_genome};
2310 :     my $newGdir = $self->{_orgdir};
2311 :    
2312 :     if ($peg1 =~ /^fig\|$newG\.peg/)
2313 :     {
2314 :     $self->load_coupling_index();
2315 : parrello 1.49
2316 : olson 1.27 my $tie = $self->{_pch_tie};
2317 :     my $evtie = $self->{_pch_ev_tie};
2318 :    
2319 : olson 1.30 if (not(defined($tie) and defined($evtie)))
2320 :     {
2321 :     return;
2322 :     }
2323 :    
2324 : olson 1.27 my @out;
2325 :     my @coupled_to = $tie->get_dup($peg1);
2326 :     for my $ent (@coupled_to)
2327 :     {
2328 :     my ($peg2, $score) = split(/$;/, $ent);
2329 : parrello 1.49
2330 : olson 1.27 my @ev = $evtie->get_dup("$peg1$;$peg2");
2331 :    
2332 :     my $l = [];
2333 :     for my $event (@ev)
2334 :     {
2335 :     my($peg3, $peg4, $iden3, $iden4, $rep) = split(/$;/, $event);
2336 :     push(@$l, [$peg3, $peg4]);
2337 :     }
2338 :     push(@out, [$score, $peg2, $l]);
2339 :     }
2340 :     return @out;
2341 :     }
2342 :     else
2343 :     {
2344 :     return $fig->coupling_and_evidence($peg1);
2345 :     }
2346 :    
2347 :     }
2348 :    
2349 : olson 1.45 sub in_pch_pin_with
2350 :     {
2351 :     my($self, $peg1, $diverse) = @_;
2352 :    
2353 :     my @all = $self->in_pch_pin_with_and_evidence($peg1);
2354 :    
2355 :     if ($diverse)
2356 :     {
2357 :     return map { $_->[0] } grep { $_->[1] == 1 } @all;
2358 :     }
2359 :     else
2360 :     {
2361 :     return map { $_->[0] } @all;
2362 :     }
2363 :     }
2364 :    
2365 :     sub in_pch_pin_with_and_evidence
2366 :     {
2367 :     my($self,$peg1) = @_;
2368 :    
2369 :     my $fig = $self->{_fig};
2370 :     my $newG = $self->{_genome};
2371 :     my $newGdir = $self->{_orgdir};
2372 :    
2373 :     if ($peg1 !~ /^fig\|$newG\.peg/)
2374 :     {
2375 :     return $fig->in_pch_pin_with_and_evidence($peg1);
2376 :     }
2377 :    
2378 :     $self->load_coupling_index();
2379 :    
2380 :     my $evtie = $self->{_pch_ev_tie};
2381 :    
2382 :     my($key, $value, $st);
2383 :    
2384 :     my %max;
2385 : parrello 1.49
2386 : olson 1.45 $key = "$peg1$;";
2387 :     my $qkey = quotemeta($key);
2388 :     for ($st = $evtie->seq($key, $value, R_CURSOR); $st == 0 and $key =~ /^$qkey/; $st = $evtie->seq($key, $value, R_NEXT))
2389 :     {
2390 :     # print "key=$key value=$value\n";
2391 :    
2392 :     my($peg3, $peg4, $iden3, $iden4, $rep) = split(/$;/, $value);
2393 :    
2394 :     if (exists($max{$peg3}))
2395 :     {
2396 :     $max{$peg3} = $rep if $rep > $max{$peg3};
2397 :     }
2398 :     else
2399 :     {
2400 :     $max{$peg3} = $rep;
2401 :     }
2402 :     }
2403 :    
2404 :     return map { [$_, $max{$_}] } keys %max;
2405 :     }
2406 :    
2407 :    
2408 : olson 1.27 sub load_coupling_index
2409 :     {
2410 :     my($self) = @_;
2411 :    
2412 :     return if defined($self->{_pch_index});
2413 :    
2414 :     my $fig = $self->{_fig};
2415 :     my $newG = $self->{_genome};
2416 :     my $newGdir = $self->{_orgdir};
2417 :    
2418 :     my $pch_btree = "$newGdir/pchs.btree";
2419 :     my $ev_btree = "$newGdir/pchs.evidence.btree";
2420 :    
2421 :     my $pch_index = {};
2422 :     my $ev_index = {};
2423 :    
2424 :     my $tied = tie %$pch_index, 'DB_File', $pch_btree, O_RDONLY, 0666, $DB_BTREE;
2425 :     my $evtied = tie %$ev_index, 'DB_File', $ev_btree, O_RDONLY, 0666, $DB_BTREE;
2426 :    
2427 :     #
2428 :     # Set these even if failed so we don't keep trying to open and failing.
2429 :     #
2430 :     $self->{_pch_index} = $pch_index;
2431 :     $self->{_pch_tie} = $tied;
2432 :     $self->{_pch_ev_index} = $ev_index;
2433 :     $self->{_pch_ev_tie} = $evtied;
2434 : parrello 1.49
2435 : olson 1.27 if (not $tied)
2436 :     {
2437 :     warn "Cannot tie pch index $pch_btree: $!\n";
2438 :     }
2439 :     }
2440 :    
2441 : olson 1.16 sub load_bbh_index
2442 :     {
2443 :     my($self) = @_;
2444 :    
2445 :     return if defined($self->{_bbh_index});
2446 :    
2447 :     my $fig = $self->{_fig};
2448 :     my $newG = $self->{_genome};
2449 :     my $newGdir = $self->{_orgdir};
2450 :    
2451 :     my $bbh_file = "$newGdir/bbhs";
2452 :     my $bbh_index_file = "$bbh_file.index";
2453 :     my $bbh_index = {};
2454 :    
2455 :     my $tied = tie %$bbh_index, 'DB_File', $bbh_index_file, O_RDONLY, 0666, $DB_BTREE;
2456 :    
2457 :     #
2458 :     # Set these even if failed so we don't keep trying to open and failing.
2459 :     #
2460 :     $self->{_bbh_index} = $bbh_index;
2461 :     $self->{_bbh_tie} = $tied;
2462 : parrello 1.49
2463 : olson 1.16 if (not $tied)
2464 :     {
2465 :     warn "Cannot tie bbh index $bbh_index_file: $!\n";
2466 :     }
2467 :     }
2468 :    
2469 : olson 1.36 sub load_contigs_index
2470 :     {
2471 :     my($self) = @_;
2472 :    
2473 :     return if defined($self->{_contigs_index});
2474 :    
2475 :     my $fig = $self->{_fig};
2476 :     my $newG = $self->{_genome};
2477 :     my $newGdir = $self->{_orgdir};
2478 :    
2479 :     my $contig_index_file = "$newGdir/contigs.btree";
2480 :     my $contig_index = {};
2481 : olson 1.37 my $contig_len_index_file = "$newGdir/contig_len.btree";
2482 :     my $contig_len_index = {};
2483 : olson 1.36
2484 :     my $tied = tie %$contig_index, 'DB_File', $contig_index_file, O_RDONLY, 0666, $DB_BTREE;
2485 : olson 1.37 if (not $tied)
2486 :     {
2487 : olson 1.47 # warn "Cannot tie contig index $contig_index_file: $!\n";
2488 : olson 1.37 }
2489 :    
2490 :     my $ltied = tie %$contig_len_index, 'DB_File', $contig_len_index_file, O_RDONLY, 0666, $DB_BTREE;
2491 :     if (not $ltied)
2492 :     {
2493 : olson 1.47 # warn "Cannot tie contig length index $contig_len_index_file: $!\n";
2494 : olson 1.37 }
2495 : olson 1.36
2496 :     #
2497 :     # Set these even if failed so we don't keep trying to open and failing.
2498 :     #
2499 : olson 1.37 $self->{_contigs_index} = $contig_index;
2500 :     $self->{_contigs_tie} = $tied;
2501 :     $self->{_contig_len_index} = $contig_len_index;
2502 :     $self->{_contig_len_tie} = $ltied;
2503 : parrello 1.49
2504 : olson 1.36 }
2505 :    
2506 : olson 1.13 sub load_sims_index
2507 :     {
2508 :     my($self) = @_;
2509 :    
2510 :     return if defined($self->{_sims_index});
2511 :    
2512 :     my $fig = $self->{_fig};
2513 :     my $newG = $self->{_genome};
2514 :     my $newGdir = $self->{_orgdir};
2515 :    
2516 :     my $sims_file = "$newGdir/expanded_similarities";
2517 :     my $sims_index_file = "$sims_file.index";
2518 :     my $sims_index = {};
2519 :    
2520 : olson 1.27 my $flips_file = "$newGdir/expanded_similarities.flips";
2521 :     my $flips_index_file = "$flips_file.index";
2522 :     my $flips_index = {};
2523 :    
2524 :     my $raw_sims_file = "$newGdir/similarities";
2525 :     my $raw_sims_index_file = "$raw_sims_file.index";
2526 :     my $raw_sims_index = {};
2527 :    
2528 :     my $raw_flips_file = "$newGdir/similarities.flips";
2529 :     my $raw_flips_index_file = "$raw_flips_file.index";
2530 :     my $raw_flips_index = {};
2531 :    
2532 :     $self->open_sims($sims_index, $sims_file, $sims_index_file, "_sims");
2533 :     $self->open_sims($flips_index, $flips_file, $flips_index_file, "_flips");
2534 :     $self->open_sims($raw_sims_index, $raw_sims_file, $raw_sims_index_file, "_raw_sims");
2535 :     $self->open_sims($raw_flips_index, $raw_flips_file, $raw_flips_index_file, "_raw_flips");
2536 :     }
2537 :    
2538 :     #
2539 :     # Open a sims file, tie it to a hash, and store the info into the $self obj.
2540 :     #
2541 :     sub open_sims
2542 :     {
2543 :     my($self, $hash, $sims_file, $index_file, $prefix) = @_;
2544 :    
2545 :     my $tied = tie %$hash, 'DB_File', $index_file, O_RDONLY, 0666, $DB_BTREE;
2546 : olson 1.13
2547 :     #
2548 :     # Set these even if failed so we don't keep trying to open and failing.
2549 :     #
2550 : olson 1.27 $self->{"${prefix}_index"} = $hash;
2551 :     $self->{"${prefix}_tie"} = $tied;
2552 : parrello 1.49
2553 : olson 1.13 if (not $tied)
2554 :     {
2555 : olson 1.27 warn "Cannot tie sims index $index_file: $!\n";
2556 : olson 1.13 }
2557 :    
2558 :     #
2559 :     # open the sims file as well.
2560 :     #
2561 :    
2562 : olson 1.27 $self->{"${prefix}_fh"} = new FileHandle("<$sims_file");
2563 : olson 1.13
2564 : olson 1.27 if (!$self->{"${prefix}_fh"})
2565 : olson 1.13 {
2566 :     warn "Cannot open sims file $sims_file: $!\n";
2567 :     }
2568 :     }
2569 :    
2570 : overbeek 1.1 sub load_sims {
2571 :     my($self) = @_;
2572 :    
2573 :     if ($self->{_sims}) { return };
2574 :    
2575 :     my $newGdir = $self->{_orgdir};
2576 : parrello 1.49
2577 : overbeek 1.1 my $sims = {};
2578 :     foreach my $x (`cat $newGdir/similarities`)
2579 :     {
2580 : olson 1.46 chomp $x;
2581 : overbeek 1.1 if ($x =~ /^(\S+)/)
2582 :     {
2583 :     push(@{$sims->{$1}},bless([split(/\t/,$x)],'Sim'));
2584 :     }
2585 :     }
2586 :     $self->{_sims} = $sims;
2587 :     }
2588 :    
2589 :     sub get_attributes {
2590 : olson 1.44 my($self, $id, $attr) = @_;
2591 : overbeek 1.1
2592 :     my $fig = $self->{_fig};
2593 :     my $newG = $self->{_genome};
2594 :     my $newGdir = $self->{_orgdir};
2595 :    
2596 : olson 1.44 if ((($id =~ /^fig\|(\d+\.\d+)\.peg\.\d+$/) and ($1 eq $newG)) or
2597 :     (($id =~ /^(\d+\.\d+)/ and $1 eq $newG)))
2598 : overbeek 1.1 {
2599 :     &load_attr($self);
2600 : olson 1.44
2601 :     my $t = $self->{_attr_id_tie};
2602 :     if (!$t)
2603 : overbeek 1.1 {
2604 : olson 1.44 #
2605 :     # Tied index not present, bail back to simple attributes.
2606 :     #
2607 :    
2608 :     my $l = $self->{_attr_id}->{$id};
2609 :     return $l ? @$l : ();
2610 : overbeek 1.1 }
2611 : olson 1.44
2612 :     my @v = $t->get_dup($id);
2613 :    
2614 :     my @out;
2615 :     for my $v (@v)
2616 : overbeek 1.1 {
2617 : olson 1.44 my @a = split(/$;/, $v);
2618 :    
2619 :     if (!defined($attr) or $attr eq $a[1])
2620 :     {
2621 :     push(@out, [@a]);
2622 :     }
2623 : overbeek 1.1 }
2624 : olson 1.44 return @out;
2625 : overbeek 1.1 }
2626 :     else
2627 :     {
2628 : olson 1.44 my @f = $fig->get_attributes($id, $attr);
2629 :     if (!defined($id) and defined(my $t = $self->{_attr_key_tie}))
2630 :     {
2631 :     #
2632 :     # lookup locally for $attr matches and merge with other output.
2633 :     #
2634 :    
2635 :     my @mine = $t->get_dup($attr);
2636 :     @mine = map { [split(/$;/, $_)] } @mine;
2637 :    
2638 :     return (@mine, @f);
2639 :     }
2640 :     else
2641 :     {
2642 :     return @f;
2643 :     }
2644 : overbeek 1.1 }
2645 :     }
2646 :    
2647 :     sub load_attr {
2648 :     my($self) = @_;
2649 :    
2650 : olson 1.44 if ($self->{_attr_id}) { return };
2651 : overbeek 1.1
2652 :     my $newGdir = $self->{_orgdir};
2653 : olson 1.44
2654 :     my $id = {};
2655 :     my $key = {};
2656 :    
2657 :     $self->{_attr_id_tie} = tie %$id, 'DB_File', "$newGdir/attr_id.btree", O_RDONLY, 0, $DB_BTREE;
2658 :     $self->{_attr_key_tie} = tie %$key, 'DB_File', "$newGdir/attr_key.btree", O_RDONLY, 0, $DB_BTREE;
2659 :    
2660 :     $self->{_attr_id} = $id;
2661 :     $self->{_attr_key} = $key;
2662 :    
2663 :     #
2664 :     # If the tie didn't work for ids, at least load up the evidence codes.
2665 :     #
2666 :    
2667 :     if (!$self->{_attr_id_tie})
2668 : overbeek 1.1 {
2669 : olson 1.44 foreach my $x (`cat $newGdir/evidence.codes`)
2670 : overbeek 1.1 {
2671 : olson 1.44 if ($x =~ /^(\S+)\t(\S+)/)
2672 :     {
2673 :     push(@{$id->{$1}},[$1,"evidence_code",$2,""]);
2674 :     }
2675 : overbeek 1.1 }
2676 :     }
2677 : olson 1.44
2678 : overbeek 1.1 }
2679 : parrello 1.49
2680 : overbeek 1.1 sub load_ann {
2681 :     my($self) = @_;
2682 :    
2683 :     if ($self->{_ann}) { return };
2684 :    
2685 :     my $newGdir = $self->{_orgdir};
2686 :     my $ann = {};
2687 :     if (open(ANN,"<$newGdir/annotations"))
2688 :     {
2689 :     $/ = "\n//\n";
2690 :     while (defined(my $x = <ANN>))
2691 :     {
2692 :     chomp $x;
2693 :     if ($x =~ /^(\S+)\n([^\n]+)\n([^\n]+)\n(.*)/s)
2694 :     {
2695 :     push(@{$ann->{$1}},[$1,$2,$3,"$4\n"]);
2696 :     }
2697 :     }
2698 :     $/ = "\n";
2699 :     close(ANN);
2700 :     }
2701 :     $self->{_ann} = $ann;
2702 :     }
2703 : overbeek 1.6
2704 :     sub taxonomy_of {
2705 :     my($self,$genome) = @_;
2706 :     my $fig = $self->{_fig};
2707 :     my $newG = $self->{_genome};
2708 :     my $newGdir = $self->{_orgdir};
2709 :    
2710 :     if ($newG ne $genome)
2711 :     {
2712 :     return $fig->taxonomy_of($genome);
2713 :     }
2714 :    
2715 :     my $tax;
2716 :     if (open(TAX,"<$newGdir/TAXONOMY") && ($tax = <TAX>))
2717 :     {
2718 :     chop $tax;
2719 :     return $tax;
2720 :     }
2721 :     else
2722 :     {
2723 :     return "unknown";
2724 :     }
2725 :     }
2726 :    
2727 :     sub build_tree_of_complete {
2728 :     my($self,$min_for_label) = @_;
2729 :     return $self->build_tree_of_all($min_for_label, "complete");
2730 :     }
2731 :    
2732 :     sub build_tree_of_all {
2733 :     my($self, $min_for_label, $complete)=@_;
2734 :     my(@last,@tax,$i,$prefix,$lev,$genome,$tax);
2735 :    
2736 :     $min_for_label = $min_for_label ? $min_for_label : 10;
2737 :     open(TMP,">/tmp/tree$$") || die "could not open /tmp/tree$$";
2738 :     print TMP "1. root\n";
2739 :    
2740 :     @last = ();
2741 :    
2742 :    
2743 :     foreach $genome (grep { ! $self->is_environmental($_) } $self->sort_genomes_by_taxonomy($self->genomes($complete)))
2744 :     {
2745 :     $tax = $self->taxonomy_of($genome);
2746 :     @tax = split(/\s*;\s*/,$tax);
2747 :     push(@tax,$genome);
2748 :     for ($i=0; ((@last > $i) && (@tax > $i) && ($last[$i] eq $tax[$i])); $i++) {}
2749 :     while ($i < @tax)
2750 :     {
2751 :     $lev = $i+2;
2752 :     $prefix = " " x (4 * ($lev-1));
2753 :     print TMP "$prefix$lev\. $tax[$i]\n";
2754 :     $i++;
2755 :     }
2756 :     @last = @tax;
2757 :     }
2758 :     close(TMP);
2759 :     my $tree = &tree_utilities::build_tree_from_outline("/tmp/tree$$");
2760 :     $tree->[0] = 'All';
2761 :     &FIG::limit_labels($tree,$min_for_label);
2762 :     unlink("/tmp/tree$$");
2763 :     return ($tree,&tree_utilities::tips_of_tree($tree));
2764 :     }
2765 :    
2766 :     sub sort_genomes_by_taxonomy {
2767 :     my($self,@genomes) = @_;
2768 :    
2769 :     return map { $_->[0] }
2770 :     sort { $a->[1] cmp $b->[1] }
2771 :     map { [$_,$self->taxonomy_of($_)] }
2772 :     @genomes;
2773 :     }
2774 : parrello 1.49
2775 : overbeek 1.6 sub taxonomic_groups_of_complete {
2776 :     my($self,$min_for_labels) = @_;
2777 :    
2778 :     my($tree,undef) = $self->build_tree_of_complete($min_for_labels);
2779 :     return &FIG::taxonomic_groups($tree);
2780 :     }
2781 :    
2782 : olson 1.31 =head2 Search Database
2783 :    
2784 :     Searches the database for objects that match the query string in some way.
2785 :    
2786 :     Returns a list of results if the query is ambiguous or an unique identifier
2787 :     otherwise.
2788 :    
2789 :     =cut
2790 :    
2791 :     sub search_database {
2792 :     # get parameters
2793 :     my ($self, $query, $options) = @_;
2794 :    
2795 :     my $fig = $self->{_fig};
2796 :     my $newG = $self->{_genome};
2797 :     my $newGdir = $self->{_orgdir};
2798 : parrello 1.49
2799 : olson 1.31 #
2800 :     # Check for local ids and genome.
2801 :     #
2802 :    
2803 : parrello 1.49
2804 : olson 1.31 if ($query eq $newG or lc($query) eq lc($self->genus_species($newG)))
2805 :     {
2806 :     return { type => 'organism', result => $newG };
2807 :     }
2808 :    
2809 :     if ($query =~ /^fig\|(\d+\.\d+)\.peg/ and $1 eq $newG)
2810 :     {
2811 :     return { type => 'feature', result => $query };
2812 :     }
2813 :    
2814 :     #
2815 :     # Match on functions
2816 :     #
2817 :    
2818 :     &load_functions($self);
2819 :    
2820 :     my @protlist;
2821 :     my $fn = $self->{_functions};
2822 :     for my $id (keys %$fn)
2823 :     {
2824 :     if ($fn->{$id} =~ /$query/i)
2825 :     {
2826 :     push @protlist, [$id, $fn->{$id}, $newG];
2827 :     }
2828 :     }
2829 :    
2830 :     #
2831 :     # Pull FIG lookups.
2832 :     #
2833 :    
2834 :     my $res = $fig->search_database($query, $options);
2835 :    
2836 :     #
2837 :     # And integrate our protein list
2838 :     #
2839 :    
2840 :     my $res_prot;
2841 :     if (ref($res) eq 'ARRAY')
2842 :     {
2843 :     for my $ent (@$res)
2844 :     {
2845 :     if ($ent->{type} eq 'proteins')
2846 :     {
2847 :     $res_prot = $ent->{result};
2848 :     }
2849 :     }
2850 :     if (!$res_prot)
2851 :     {
2852 :     $res_prot = {type => 'proteins', result => \@protlist}
2853 :     }
2854 :     else
2855 :     {
2856 :     push(@$res_prot, @protlist);
2857 :     }
2858 :     }
2859 : parrello 1.49
2860 : olson 1.31 return $res;
2861 : parrello 1.49
2862 : olson 1.31 }
2863 :    
2864 :    
2865 : paczian 1.29 =head3 model_directory
2866 : paczian 1.26
2867 : parrello 1.49 FIG->model_directory($organism);
2868 : paczian 1.26
2869 : paczian 1.29 Returns the model directory of an organism.
2870 : paczian 1.26
2871 :     =over 4
2872 :    
2873 :     =item $organism
2874 :    
2875 :     The seed-taxonomy id of the organism, e.g. 83333.1. If you pass the
2876 :     string 'All', you will get the directory for the global model.
2877 :    
2878 :     =back
2879 :    
2880 :     =cut
2881 :    
2882 : paczian 1.29 sub model_directory {
2883 : paczian 1.26 my ($self, $organism) = @_;
2884 :    
2885 : olson 1.27 my $directory;
2886 : paczian 1.26
2887 : olson 1.27 if ($organism eq $self->{_genome})
2888 :     {
2889 :     $directory = $self->{_orgdir} . "/Models";
2890 :     $directory .= "/$organism" if defined($organism);
2891 : formsma 1.38 }
2892 :     elsif(!defined $organism && defined $self->{_genome})
2893 :     {
2894 :     $directory = $self->{_orgdir} . "/Models";
2895 :     }
2896 :     else {
2897 : paczian 1.29 $directory = $self->{_fig}->model_directory($organism);
2898 : paczian 1.26 }
2899 :    
2900 :     return $directory;
2901 :     }
2902 :    
2903 : overbeek 1.1 1;

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3