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

Annotation of /FigKernelPackages/FIGV.pm

Parent Directory Parent Directory | Revision Log Revision Log


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

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3