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

Annotation of /FigKernelPackages/FIGV.pm

Parent Directory Parent Directory | Revision Log Revision Log


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

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

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3