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

Annotation of /FigKernelPackages/FIGV.pm

Parent Directory Parent Directory | Revision Log Revision Log


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

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

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3