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

Annotation of /FigKernelPackages/FIGV.pm

Parent Directory Parent Directory | Revision Log Revision Log


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

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3