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

Annotation of /FigKernelPackages/FIGV.pm

Parent Directory Parent Directory | Revision Log Revision Log


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

1 : overbeek 1.1 #
2 :     # Copyright (c) 2003-2006 University of Chicago and Fellowship
3 :     # for Interpretations of Genomes. All Rights Reserved.
4 :     #
5 :     # This file is part of the SEED Toolkit.
6 :     #
7 :     # The SEED Toolkit is free software. You can redistribute
8 :     # it and/or modify it under the terms of the SEED Toolkit
9 :     # Public License.
10 :     #
11 :     # You should have received a copy of the SEED Toolkit Public License
12 :     # along with this program; if not write to the University of Chicago
13 :     # at info@ci.uchicago.edu or the Fellowship for Interpretation of
14 :     # Genomes at veronika@thefig.info or download a copy from
15 :     # http://www.theseed.org/LICENSE.TXT.
16 :     #
17 :    
18 :     package FIGV;
19 :    
20 : overbeek 1.6 use Carp;
21 : overbeek 1.1 use strict;
22 :     use FIG;
23 :     use FIG_Config;
24 :     use SFXlate;
25 :     use SproutFIG;
26 :     use Tracer;
27 :     use Data::Dumper;
28 : olson 1.4 use vars qw($AUTOLOAD);
29 : olson 1.13 use DB_File;
30 :     use FileHandle;
31 : overbeek 1.1
32 :     sub new {
33 :     my($class,$org_dir,$low_level) = @_;
34 :    
35 :     my $fig;
36 :     if ($low_level && ($low_level =~ /sprout/i))
37 :     {
38 :     $fig = new SproutFIG($FIG_Config::sproutDB, $FIG_Config::sproutData);
39 :     }
40 : overbeek 1.17 elsif (! $org_dir)
41 :     {
42 :     $fig = new FIG;
43 :     return $fig;
44 :     }
45 : overbeek 1.1 else
46 :     {
47 :     $fig = new FIG;
48 :     }
49 :    
50 :     my $self = {};
51 :     $self->{_fig} = $fig;
52 :     $self->{_orgdir} = $org_dir;
53 :    
54 :     $org_dir =~ /(\d+\.\d+)$/;
55 :     $self->{_genome} = $1;
56 :     return bless $self, $class;
57 :     }
58 :    
59 : paarmann 1.20
60 :     # return the path of the organism directory
61 :     sub organism_directory {
62 :     return $_[0]->{_orgdir};
63 :     }
64 :    
65 : olson 1.27 sub is_complete
66 :     {
67 :     return 1;
68 :     }
69 :    
70 : olson 1.4 #
71 :     # Redirect any method invocations that we don't implement out to the
72 :     # underlying FIG object.
73 :     #
74 :     sub AUTOLOAD
75 :     {
76 :     my($self, @args) = @_;
77 :    
78 : overbeek 1.6 if (ref($self) ne "FIGV") {
79 :     confess "BAD FIGV object passed to AUTOLOAD";
80 :     }
81 :    
82 : olson 1.27 no strict 'refs';
83 :    
84 : olson 1.4 my $meth = $AUTOLOAD;
85 :     $meth =~ s/.*:://;
86 : olson 1.27 my $fmeth = "FIG::$meth";
87 : olson 1.4
88 :     my $fig = $self->{_fig};
89 : olson 1.8 # my $args = Dumper(\@args);
90 : olson 1.4 if (wantarray)
91 :     {
92 :     my @res = $fig->$meth(@args);
93 : olson 1.27 # my @res = &$fmeth($self, @args);
94 : olson 1.5 # warn "FIGV invoke $meth($args) returns\n", Dumper(\@res);
95 : olson 1.4 return @res;
96 :     }
97 :     else
98 :     {
99 :     my $res = $fig->$meth(@args);
100 : olson 1.27 # my $res = &$fmeth($self, @args);
101 : olson 1.5 # warn "FIGV invoke $meth($args) returns\n", Dumper($res);
102 : olson 1.4 return $res;
103 :     }
104 :     }
105 :    
106 : olson 1.8 sub FIG
107 :     {
108 :     my($self) = @_;
109 :     return $self;
110 :     }
111 :    
112 : olson 1.27 sub sort_fids_by_taxonomy
113 :     {
114 :     my($self,@fids) = @_;
115 :    
116 :     return map { $_->[0] }
117 :     sort { $a->[1] cmp $b->[1] }
118 :     map { [$_,$self->taxonomy_of($self->genome_of($_))] }
119 :     @fids;
120 :     }
121 : olson 1.8
122 : olson 1.34 sub get_basic_statistics
123 :     {
124 :     my($self, $genome) = @_;
125 :    
126 :     my $fig = $self->{_fig};
127 :     my $newG = $self->{_genome};
128 :     my $newGdir = $self->{_orgdir};
129 :    
130 :     if ($genome ne $newG)
131 :     {
132 :     return $fig->get_basic_statistics($genome);
133 :     }
134 :    
135 :     #
136 :     # Check cache.
137 :     #
138 :    
139 :     my $cache = "$newGdir/cache.basic_statistics";
140 :     my $fh = new FileHandle($cache);
141 :     if ($fh)
142 :     {
143 :     my $stats = {};
144 :     while (<$fh>)
145 :     {
146 :     chomp;
147 :     my($k, $v) = split(/\t/);
148 :     $stats->{$k} = $v;
149 :     }
150 :     close($fh);
151 :     return $stats;
152 :     }
153 :    
154 :    
155 :     my $subsystem_data = $self->get_genome_subsystem_data($genome);
156 :    
157 :     my %sscount = map { $_->[0] => 1 } @$subsystem_data;
158 :     my $nss=scalar(keys(%sscount));
159 :    
160 :     my $statistics = {
161 :     num_subsystems => $nss,
162 :     num_contigs => scalar($self->all_contigs($genome)),
163 :     num_basepairs => $self->genome_szdna($genome),
164 :     genome_name => $self->genus_species($genome),
165 :     genome_domain => $self->genome_domain($genome),
166 :     genome_pegs => $self->genome_pegs($genome),
167 :     genome_rnas => $self->genome_rnas($genome),
168 :     genome_version => $self->genome_version($genome)
169 :     };
170 :    
171 :     $fh = new FileHandle(">$cache");
172 :     if ($fh)
173 :     {
174 :     while (my($k, $v) = each %$statistics)
175 :     {
176 :     print $fh join("\t", $k, $v), "\n";
177 :     }
178 :     close($fh);
179 :     }
180 :    
181 :     return $statistics;
182 :     }
183 :    
184 :    
185 :     sub get_peg_statistics {
186 :     my ($self, $genome) = @_;
187 :    
188 :     my $fig = $self->{_fig};
189 :     my $newG = $self->{_genome};
190 :     my $newGdir = $self->{_orgdir};
191 :    
192 :     if ($genome ne $newG)
193 :     {
194 :     return $fig->get_peg_statistics($genome);
195 :     }
196 :    
197 :     #
198 :     # Check cache.
199 :     #
200 :    
201 :     my $cache = "$newGdir/cache.peg_statistics";
202 :     my $fh = new FileHandle($cache);
203 :     if ($fh)
204 :     {
205 :     my $stats = {};
206 :     while (<$fh>)
207 :     {
208 :     chomp;
209 :     my($k, $v) = split(/\t/);
210 :     $stats->{$k} = $v;
211 :     }
212 :     close($fh);
213 :     return $stats;
214 :     }
215 :    
216 :    
217 :     my $subsystem_data = $self->get_genome_subsystem_data($genome);
218 :     my $assignment_data = $self->get_genome_assignment_data($genome);
219 :    
220 :     my $hypo_sub = 0;
221 :     my $hypo_nosub = 0;
222 :     my $nothypo_sub = 0;
223 :     my $nothypo_nosub = 0;
224 :     my %in = map { $_->[2] => 1 } @$subsystem_data;
225 :     my $in = keys(%in);
226 :    
227 :     my %sscount = map { $_->[0] => 1 } @$subsystem_data;
228 :    
229 :     foreach $_ (@$assignment_data)
230 :     {
231 :     my($peg,$func) = @$_;
232 :     my $is_hypo = &FIG::hypo($func);
233 :    
234 :     if ($is_hypo && $in{$peg}) { $hypo_sub++ }
235 :     elsif ($is_hypo && ! $in{$peg}) { $hypo_nosub++ }
236 :     elsif ((! $is_hypo) && (! $in{$peg})) { $nothypo_nosub++ }
237 :     elsif ((! $is_hypo) && $in{$peg}) { $nothypo_sub++ }
238 :     }
239 :     my $tot = $hypo_sub + $nothypo_sub + $hypo_nosub + $nothypo_nosub;
240 :    
241 :     my ($fracHS, $fracNHS, $fracHNS, $fracNHNS);
242 :    
243 :     if ($tot == 0) {
244 :     $fracHS = sprintf "%.2f", 0.0;
245 :     $fracNHS = sprintf "%.2f", 0.0;
246 :     $fracHNS = sprintf "%.2f", 0.0;
247 :     $fracNHNS = sprintf "%.2f", 0.0;
248 :     } else {
249 :     $fracHS = sprintf "%.2f", $hypo_sub / $tot * 100;
250 :     $fracNHS = sprintf "%.2f", $nothypo_sub / $tot * 100;
251 :     $fracHNS = sprintf "%.2f", $hypo_nosub / $tot * 100;
252 :     $fracNHNS = sprintf "%.2f", $nothypo_nosub / $tot * 100;
253 :     }
254 :    
255 :     my $statistics = {
256 :     hypothetical_in_subsystem => $hypo_sub,
257 :     hypothetical_not_in_subsystem => $hypo_nosub,
258 :     non_hypothetical_in_subsystem => $nothypo_sub,
259 :     non_hypothetical_not_in_subsystem => $nothypo_nosub,
260 :     hypothetical_in_subsystem_percent => $fracHS,
261 :     hypothetical_not_in_subsystem_percent => $fracHNS,
262 :     non_hypothetical_in_subsystem_percent => $fracNHS,
263 :     non_hypothetical_not_in_subsystem_percent => $fracNHNS
264 :     };
265 :    
266 :     $fh = new FileHandle(">$cache");
267 :     if ($fh)
268 :     {
269 :     while (my($k, $v) = each %$statistics)
270 :     {
271 :     print $fh join("\t", $k, $v), "\n";
272 :     }
273 :     close($fh);
274 :     }
275 :    
276 :     return $statistics;
277 :     }
278 :    
279 : olson 1.5 #
280 :     # To retrieve a subsystem in FIGV, we create the subsystem as normal via $fig->get_subsystem,
281 :     # then insert the row for the virtual org dir we are processing.
282 :     #
283 :    
284 :     sub get_subsystem
285 :     {
286 :     my($self,$ssa) = @_;
287 :    
288 :     my $fig = $self->{_fig};
289 :     my $newG = $self->{_genome};
290 :     my $newGdir = $self->{_orgdir};
291 :    
292 :     my $ss = $fig->get_subsystem($ssa);
293 :     return undef unless $ss;
294 :    
295 : olson 1.9 $self->load_ss_data();
296 : olson 1.8
297 : olson 1.9 my $bindings = $self->{_ss_bindings}->{$ssa};
298 :     my $variant = $self->{_ss_variants}->{$ssa};
299 :    
300 :     # warn "Adding virtual genome " . Dumper(\%bindings);
301 :     $ss->add_virtual_genome($self->genus_species(), $newG, $variant, $bindings);
302 :    
303 :     return $ss;
304 :     }
305 :    
306 : olson 1.14 sub active_subsystems
307 :     {
308 :     my($self, $genome, $all) = @_;
309 :    
310 :     my $fig = $self->{_fig};
311 :     my $newG = $self->{_genome};
312 :     my $newGdir = $self->{_orgdir};
313 :    
314 :     if ($genome ne $newG)
315 :     {
316 :     return $fig->active_subsystems($genome, $all);
317 :     }
318 :    
319 :     $self->load_ss_data();
320 :    
321 :     my $slist = {};
322 : olson 1.25
323 :     if ($self->{_ss_variants})
324 :     {
325 :     %{$slist} = %{$self->{_ss_variants}};
326 :     }
327 : olson 1.14
328 :     if (not $all)
329 :     {
330 :     for my $ss (keys %$slist)
331 :     {
332 :     my $var = $slist->{$ss};
333 : paarmann 1.32 delete $slist->{$ss} if $var eq 0 or $var eq -1;
334 : olson 1.14 }
335 :     }
336 :     return $slist;
337 :     }
338 :    
339 : overbeek 1.19 sub peg_to_subsystems
340 :     {
341 :     my($self, $peg) = @_;
342 :    
343 :     my $variant = $self->{_ss_variants};
344 :     my %sub = map { $_ => 1 }
345 :     grep { $variant->{$_} !~ /^(-1)|0$/ }
346 :     map { $_->[0] }
347 :     $self->peg_to_roles_in_subsystems($peg);
348 :     return sort keys(%sub);
349 :     }
350 :    
351 : olson 1.9 sub peg_to_roles_in_subsystems
352 :     {
353 :     my($self,$peg) = @_;
354 :    
355 :     my $fig = $self->{_fig};
356 :     my $newG = $self->{_genome};
357 :     my $newGdir = $self->{_orgdir};
358 : olson 1.8
359 : olson 1.22 if ($peg !~ /fig\|$newG\.peg/)
360 :     {
361 :     return $fig->peg_to_roles_in_subsystems($peg);
362 :     }
363 :    
364 : olson 1.9 $self->load_ss_data();
365 : olson 1.8
366 : olson 1.9 my $ret = $self->{_ss_peg_index}->{$peg};
367 : olson 1.8
368 : olson 1.9 return $ret ? @$ret : ();
369 :     }
370 : olson 1.8
371 : olson 1.9 sub subsystems_for_peg
372 :     {
373 :     my($self,$peg) = @_;
374 :     return $self->peg_to_roles_in_subsystems($peg);
375 : olson 1.5 }
376 : olson 1.4
377 : overbeek 1.1 sub genomes {
378 :     my($self,$complete) = @_;
379 :     my $fig = $self->{_fig};
380 :    
381 :     return ($self->{_genome},$fig->genomes($complete));
382 :     }
383 :    
384 :     sub genus_species {
385 :     my($self,$genome) = @_;
386 :    
387 :     my $fig = $self->{_fig};
388 :     my $newG = $self->{_genome};
389 :     my $newGdir = $self->{_orgdir};
390 :    
391 :     if (($genome eq $newG) && open(GENOME,"<$newGdir/GENOME"))
392 :     {
393 :     my $x = <GENOME>;
394 :     close(GENOME);
395 :     chop $x;
396 :     return $x;
397 :     }
398 :     else
399 :     {
400 :     return $fig->genus_species($genome);
401 :     }
402 :     }
403 :    
404 : overbeek 1.7 sub get_genome_assignment_data {
405 :     my($self,$genome) = @_;
406 :    
407 :     my $fig = $self->{_fig};
408 :     my $newG = $self->{_genome};
409 :     my $newGdir = $self->{_orgdir};
410 :    
411 :     if ($genome eq $newG)
412 :     {
413 : olson 1.31 my @fnfiles = <$newGdir/proposed*functions>;
414 :     if (@fnfiles == 0)
415 :     {
416 :     @fnfiles = <$newGdir/assigned*functions>;
417 :     }
418 :    
419 :     if (@fnfiles == 0)
420 :     {
421 :     return [];
422 :     }
423 : overbeek 1.7 my %assign = map { ( $_ =~ /^(fig\|\d+\.\d+\.peg\.\d+)\t(\S.*\S)/) ? ($1 => $2) : () }
424 : olson 1.31 `cat @fnfiles`;
425 : overbeek 1.7 return [map { [$_,$assign{$_}] } sort { &FIG::by_fig_id($a,$b) } keys(%assign)];
426 :     }
427 :     else
428 :     {
429 :     return $fig->get_genome_assignment_data($genome);
430 :     }
431 :     }
432 :    
433 :     sub org_of {
434 :     my($self,$peg) = @_;
435 :    
436 :     if ($peg =~ /^fig\|(\d+\.\d+)\.peg\.\d+/)
437 :     {
438 :     return $self->genus_species($1);
439 :     }
440 :     return "";
441 :     }
442 :    
443 :     sub get_genome_subsystem_data {
444 :     my($self,$genome) = @_;
445 :    
446 :     my $fig = $self->{_fig};
447 :     my $newG = $self->{_genome};
448 :     my $newGdir = $self->{_orgdir};
449 :    
450 :     if ($genome eq $newG)
451 :     {
452 :     return [map { ($_ =~ /^(\S[^\t]+\S)\t(\S[^\t]*\S)\t(\S+)/) ? [$1,$2,$3] : () }
453 :     `cat $newGdir/Subsystems/bindings`];
454 :     }
455 :     else
456 :     {
457 :     return $fig->get_genome_subsystem_data($genome);
458 :     }
459 :     }
460 :    
461 :     sub orgname_of_orgid {
462 :     my($self,$genome) = @_;
463 :    
464 :     return $self->genus_species($genome);
465 :     }
466 :    
467 :     sub orgid_of_orgname {
468 :     my($self,$genome_name) = @_;
469 :    
470 :     my @genomes = $self->genomes('complete');
471 :     my $i;
472 :     for ($i=0; ($i < @genomes) && ($genome_name ne $self->genus_species($genomes[$i])); $i++) {}
473 :     return ($i == @genomes) ? undef : $genomes[$i];
474 :     }
475 :    
476 :     sub genus_species_domain {
477 :     my($self,$genome) = @_;
478 :    
479 :     return [$self->genus_species($genome),$self->genome_domain($genome)];
480 :     }
481 :    
482 :     sub protein_subsystem_to_roles {
483 :     my ($self,$peg,$subsystem) = @_;
484 :    
485 :     my $fig = $self->{_fig};
486 :     my $newG = $self->{_genome};
487 :     my $newGdir = $self->{_orgdir};
488 :    
489 :     if (&FIG::genome_of($peg) ne $newG)
490 :     {
491 :     return $fig->protein_subsystem_to_roles($peg,$subsystem);
492 :     }
493 :     else
494 :     {
495 :     my @roles = map { (($_ =~ /^([^\t]+)\t([^\t]+)\t(\S+)$/) && ($1 eq $subsystem) && ($3 eq $peg)) ?
496 :     $2 : () } `cat $newGdir/Subsystems/bindings`;
497 :     my %roles = map { $_ => 1 } @roles;
498 :     return [sort keys(%roles)];
499 :     }
500 :     }
501 :    
502 :     sub contig_lengths {
503 :     my ($self, $genome) = @_;
504 :    
505 :     my $fig = $self->{_fig};
506 :     my $newG = $self->{_genome};
507 :     my $newGdir = $self->{_orgdir};
508 :    
509 :     if ($genome ne $newG)
510 :     {
511 :     return $fig->contig_lengths($genome);
512 :     }
513 :     else
514 :     {
515 : olson 1.37 my $contig_lengths = $self->{_contig_len_index};
516 :    
517 : olson 1.39 if (!tied(%$contig_lengths))
518 : overbeek 1.7 {
519 : olson 1.37 $contig_lengths = {};
520 :     if (open(CONTIGS,"<$newGdir/contigs"))
521 :     {
522 :     $/ = "\n>";
523 :     while (defined(my $x = <CONTIGS>))
524 :     {
525 : olson 1.40 chomp $x;
526 : olson 1.37 if ($x =~ />?(\S+)[^\n]*\n(.*)/s)
527 :     {
528 :     my $id = $1;
529 :     my $seq = $2;
530 :     $seq =~ s/\s//gs;
531 :     $contig_lengths->{$id} = length($seq);
532 :     }
533 :     }
534 :     close(CONTIGS);
535 :     $/ = "\n";
536 :     }
537 : overbeek 1.7 }
538 :     return $contig_lengths;
539 :     }
540 :     }
541 :    
542 : olson 1.10 sub contig_ln {
543 :     my ($self, $genome, $contig) = @_;
544 :    
545 :     my $fig = $self->{_fig};
546 :     my $newG = $self->{_genome};
547 :     my $newGdir = $self->{_orgdir};
548 :    
549 :     if ($genome ne $newG)
550 :     {
551 :     return $fig->contig_ln($genome, $contig);
552 :     }
553 :     else
554 :     {
555 : olson 1.37 my $contig_len = $self->{_contig_len_index};
556 :    
557 : olson 1.39 if (tied(%$contig_len))
558 : olson 1.37 {
559 :     return $contig_len->{$contig};
560 :     }
561 :    
562 : olson 1.10 if (open(CONTIGS,"<$newGdir/contigs"))
563 :     {
564 :     local $/ = "\n>";
565 :     while (defined(my $x = <CONTIGS>))
566 :     {
567 : olson 1.33 chomp $x;
568 : olson 1.10 if ($x =~ />?(\S+)[^\n]*\n(.*)/s)
569 :     {
570 :     my $id = $1;
571 :     my $seq = $2;
572 :     $seq =~ s/\s//gs;
573 :     if ($id eq $contig)
574 :     {
575 :     return length($seq);
576 :     }
577 :     }
578 :     }
579 :     close(CONTIGS);
580 :     }
581 :     }
582 :     }
583 :    
584 :     sub contigs_of
585 :     {
586 :     my ($self, $genome) = @_;
587 :    
588 :     my $fig = $self->{_fig};
589 :     my $newG = $self->{_genome};
590 :     my $newGdir = $self->{_orgdir};
591 :    
592 :     if ($genome ne $newG)
593 :     {
594 :     return $fig->contigs_of($genome);
595 :     }
596 :     else
597 :     {
598 :     my @out;
599 : olson 1.37 $self->load_contigs_index();
600 : olson 1.36
601 :     my $contigs = $self->{_contigs_index};
602 : olson 1.39 if (tied(%$contigs))
603 : olson 1.36 {
604 :     return keys %$contigs;
605 :     }
606 :    
607 : olson 1.10 if (open(CONTIGS,"<$newGdir/contigs"))
608 :     {
609 :     local $/ = "\n>";
610 :     while (defined(my $x = <CONTIGS>))
611 :     {
612 : olson 1.33 chomp $x;
613 : olson 1.10 if ($x =~ />?(\S+)[^\n]*\n(.*)/s)
614 :     {
615 :     my $id = $1;
616 :     push(@out, $id);
617 :     }
618 :     }
619 :     close(CONTIGS);
620 :     }
621 :     return @out;
622 :     }
623 :     }
624 :    
625 :     =head3 dna_seq
626 :    
627 :     usage: $seq = dna_seq($genome,@locations)
628 :    
629 :     Returns the concatenated subsequences described by the list of locations. Each location
630 :     must be of the form
631 :    
632 :     Contig_Beg_End
633 :    
634 :     where Contig must be the ID of a contig for genome $genome. If Beg > End the location
635 :     describes a stretch of the complementary strand.
636 :    
637 :     =cut
638 :     #: Return Type $;
639 :     sub dna_seq {
640 :     my($self,$genome,@locations) = @_;
641 :    
642 :     my $fig = $self->{_fig};
643 :     my $newG = $self->{_genome};
644 :     my $newGdir = $self->{_orgdir};
645 :    
646 :     if ($genome ne $newG)
647 :     {
648 :     return $fig->dna_seq($genome, @locations);
649 :     }
650 :    
651 : olson 1.36 my $contigs = $self->{_contigs_index};
652 : olson 1.41 if (!tied %$contigs)
653 : olson 1.10 {
654 : olson 1.36 $contigs = {};
655 :     if (open(CONTIGS,"<$newGdir/contigs"))
656 : olson 1.10 {
657 : olson 1.36 local $/ = "\n>";
658 :     while (defined(my $x = <CONTIGS>))
659 : olson 1.10 {
660 : olson 1.36 chomp $x;
661 :     if ($x =~ />?(\S+)[^\n]*\n(.*)/s)
662 :     {
663 :     my $id = $1;
664 :     my $seq = $2;
665 :     $seq =~ s/\s//gs;
666 :     $contigs->{$id} = $seq;
667 :     }
668 : olson 1.10 }
669 : olson 1.36 close(CONTIGS);
670 : olson 1.10 }
671 :     }
672 :    
673 :     my(@pieces,$loc,$contig,$beg,$end,$ln,$rdbH);
674 :    
675 :     @locations = map { split(/,/,$_) } @locations;
676 :     @pieces = ();
677 :     foreach $loc (@locations)
678 :     {
679 :     if ($loc =~ /^(\S+)_(\d+)_(\d+)$/)
680 :     {
681 :     ($contig,$beg,$end) = ($1,$2,$3);
682 : olson 1.36 my $seq = $contigs->{$contig};
683 : olson 1.10
684 :     $ln = length($seq);
685 :    
686 :     if (! $ln) {
687 :     print STDERR "$genome/$contig: could not get length\n";
688 :     return "";
689 :     }
690 :    
691 :     if (&FIG::between(1,$beg,$ln) && &FIG::between(1,$end,$ln))
692 :     {
693 :     if ($beg < $end)
694 :     {
695 : paczian 1.15 push(@pieces, substr($seq, $beg - 1, ($end - $beg) + 1));
696 : olson 1.10 }
697 :     else
698 :     {
699 : paczian 1.15 push(@pieces, &FIG::reverse_comp(substr($seq, $end - 1, ($beg - $end) + 1)));
700 : olson 1.10 }
701 :     }
702 :     }
703 :     }
704 :     return lc(join("",@pieces));
705 :     }
706 :    
707 : overbeek 1.7 sub genome_szdna {
708 :     my ($self, $genome) = @_;
709 :    
710 :     my $fig = $self->{_fig};
711 :     my $newG = $self->{_genome};
712 :     my $newGdir = $self->{_orgdir};
713 :    
714 :     if ($genome ne $newG)
715 :     {
716 :     return $fig->genome_szdna($genome);
717 :     }
718 :     else
719 :     {
720 :     my $contig_lens = $self->contig_lengths($genome);
721 :     my $tot = 0;
722 :     while ( my($contig,$len) = each %$contig_lens)
723 :     {
724 :     $tot += $len;
725 :     }
726 :     return $tot;
727 :     }
728 :     }
729 :    
730 :     sub genome_version {
731 :     my ($self, $genome) = @_;
732 :    
733 :     my $fig = $self->{_fig};
734 :     my $newG = $self->{_genome};
735 :     my $newGdir = $self->{_orgdir};
736 :    
737 :     if ($genome ne $newG)
738 :     {
739 :     return $fig->genome_version($genome);
740 :     }
741 :     else
742 :     {
743 :     return "$genome.0";
744 :     }
745 :     }
746 :    
747 :     sub genome_pegs {
748 :     my ($self, $genome) = @_;
749 :    
750 :     my $fig = $self->{_fig};
751 :     my $newG = $self->{_genome};
752 :     my $newGdir = $self->{_orgdir};
753 :    
754 :     if ($genome ne $newG)
755 :     {
756 :     return $fig->genome_pegs($genome);
757 :     }
758 :     else
759 :     {
760 :     my @tmp = $self->all_features($genome,"peg");
761 :     my $n = @tmp;
762 :     return $n;
763 :     }
764 :     }
765 :    
766 :     sub genome_rnas {
767 :     my ($self, $genome) = @_;
768 :    
769 :     my $fig = $self->{_fig};
770 :     my $newG = $self->{_genome};
771 :     my $newGdir = $self->{_orgdir};
772 :    
773 :     if ($genome ne $newG)
774 :     {
775 :     return $fig->genome_rnas($genome);
776 :     }
777 :     else
778 :     {
779 :     my @tmp = $self->all_features($genome,"rna");
780 :     my $n = @tmp;
781 :     return $n;
782 :     }
783 :     }
784 :    
785 :     sub genome_domain {
786 :     my ($self, $genome) = @_;
787 :    
788 :     my $fig = $self->{_fig};
789 :     my $newG = $self->{_genome};
790 :     my $newGdir = $self->{_orgdir};
791 :    
792 :     if ($genome ne $newG)
793 :     {
794 :     return $fig->genome_domain($genome);
795 :     }
796 :     else
797 :     {
798 :     my $tax = $self->taxonomy_of($genome);
799 :     return ($tax =~ /^([^ \t;]+)/) ? $1 : "unknown";
800 :     }
801 :     }
802 : overbeek 1.2
803 :     sub genes_in_region {
804 :     my($self,$genome,$contig,$beg,$end) = @_;
805 :    
806 :     my $fig = $self->{_fig};
807 :     my $newG = $self->{_genome};
808 :     my $newGdir = $self->{_orgdir};
809 :    
810 :     if ($genome ne $newG)
811 :     {
812 :     return $fig->genes_in_region($genome,$contig,$beg,$end);
813 :     }
814 :     else
815 :     {
816 :     &load_tbl($self);
817 :     my $tblH = $self->{_tbl};
818 :     my $maxV = 0;
819 :     my $minV = 1000000000;
820 :     my $genes = [];
821 :     while ( my($fid,$tuple) = each %$tblH)
822 :     {
823 :     if (($tuple->[0]->[0] =~ /^(\S+)_(\d+)_\d+$/) && ($1 eq $contig))
824 :     {
825 :     my $beg1 = $2;
826 :     my $last = @{$tuple->[0]} - 1;
827 :     if (($tuple->[0]->[$last] =~ /^(\S+)_\d+_(\d+)$/) && ($1 eq $contig))
828 :     {
829 :     my $end1 = $2;
830 :     if (&overlaps($beg,$end,$beg1,$end1))
831 :     {
832 :     $minV = &FIG::min($minV,$beg1,$end1);
833 :     $maxV = &FIG::max($maxV,$beg1,$end1);
834 :     push(@$genes,$fid);
835 :     }
836 :     }
837 :     }
838 :     }
839 :     return ($genes,$minV,$maxV);
840 :     }
841 :     }
842 :    
843 :     sub overlaps {
844 :     my($b1,$e1,$b2,$e2) = @_;
845 :    
846 :     if ($b1 > $e1) { ($b1,$e1) = ($e1,$b1) }
847 :     if ($b2 > $e2) { ($b2,$e2) = ($e2,$b2) }
848 :     return &FIG::between($b1,$b2,$e1) || &FIG::between($b2,$b1,$e2);
849 :     }
850 :    
851 : overbeek 1.7 sub all_contigs {
852 :     my($self,$genome) = @_;
853 :    
854 :     my $fig = $self->{_fig};
855 :     my $newG = $self->{_genome};
856 :     my $newGdir = $self->{_orgdir};
857 :    
858 :     if ($genome ne $newG)
859 :     {
860 :     return $fig->all_contigs($genome);
861 :     }
862 :     else
863 :     {
864 :     &load_tbl($self);
865 :     my %contigs;
866 :     my $tblH = $self->{_tbl};
867 :     while ( my($fid,$tuple) = each %$tblH)
868 :     {
869 :     if ($tuple->[0]->[0] =~ /^(\S+)_\d+_\d+$/)
870 :     {
871 :     $contigs{$1} = 1;
872 :     }
873 :     }
874 :     return keys(%contigs);
875 :     }
876 :     }
877 :    
878 :     sub all_features {
879 :     my($self,$genome,$type) = @_;
880 :    
881 :     my $fig = $self->{_fig};
882 :     my $newG = $self->{_genome};
883 :     my $newGdir = $self->{_orgdir};
884 :    
885 :     if ($genome ne $newG)
886 :     {
887 :     return $fig->all_features($genome,$type);
888 :     }
889 :     else
890 :     {
891 :     &load_tbl($self);
892 :     my $tblH = $self->{_tbl};
893 :     return sort { &FIG::by_fig_id($a,$b) }
894 :     grep { ($_ =~ /^fig\|\d+\.\d+\.([^\.]+)/) && ($1 eq $type) }
895 :     keys(%$tblH);
896 :     }
897 :     }
898 :    
899 :     sub all_features_detailed_fast {
900 : olson 1.11 my($self,$genome, $regmin, $regmax, $contig) = @_;
901 : overbeek 1.7
902 :     my $fig = $self->{_fig};
903 :     my $newG = $self->{_genome};
904 :     my $newGdir = $self->{_orgdir};
905 :    
906 :     if ($genome ne $newG)
907 :     {
908 : olson 1.11 return $fig->all_features_detailed_fast($genome, $regmin, $regmax, $contig);
909 : overbeek 1.7 }
910 :     else
911 :     {
912 :     &load_tbl($self);
913 :     my $tblH = $self->{_tbl};
914 :     my $feat_details = [];
915 :     while ( my($fid,$tuple) = each %$tblH)
916 :     {
917 :     if ($fid =~ /^fig\|\d+\.\d+\.(\S+)\.\d+/)
918 :     {
919 :     my $type = $1;
920 : olson 1.11 if ($tuple->[0]->[0] =~ /^(\S+)_(\d+)_(\d+)$/)
921 : overbeek 1.7 {
922 : olson 1.11 my($ctg, $beg, $end) = ($1, $2, $3);
923 :     next if (defined($contig) and $contig ne $ctg);
924 :    
925 : overbeek 1.7 my($min,$max);
926 : olson 1.11 if ($beg < $end)
927 : overbeek 1.7 {
928 : olson 1.11 $min = $beg;
929 :     $max = $end;
930 : overbeek 1.7 }
931 :     else
932 :     {
933 : olson 1.11 $min = $end;
934 :     $max = $beg;
935 :     }
936 :    
937 :     if (not defined($regmin) or not defined($regmax) or
938 :     ($min < $regmin and $regmin < $max) or ($min < $regmax and $regmax < $max) or ($min > $regmin and $max < $regmax))
939 :     {
940 :     push(@$feat_details,[$fid,$tuple->[0]->[0],join(",",@{$tuple->[1]}),$type,$min,$max,$self->function_of($fid),'master','']);
941 : overbeek 1.7 }
942 :     }
943 :     }
944 :     }
945 :     return $feat_details;
946 :     }
947 :     }
948 :    
949 :     sub compute_clusters {
950 :     # Get the parameters.
951 :     my ($self, $pegList, $subsystem, $distance) = @_;
952 :     if (! defined $distance) {
953 :     $distance = 5000;
954 :     }
955 :    
956 :     my($peg,%by_contig);
957 :     foreach $peg (@$pegList)
958 :     {
959 :     my $loc;
960 :     if ($loc = $self->feature_location($peg))
961 :     {
962 :     my ($contig,$beg,$end) = &FIG::boundaries_of($loc);
963 :     my $genome = &FIG::genome_of($peg);
964 :     push(@{$by_contig{"$genome\t$contig"}},[($beg+$end)/2,$peg]);
965 :     }
966 :     }
967 :    
968 :     my @clusters = ();
969 :     foreach my $tuple (keys(%by_contig))
970 :     {
971 :     my $x = $by_contig{$tuple};
972 :     my @pegs = sort { $a->[0] <=> $b->[0] } @$x;
973 :     while ($x = shift @pegs)
974 :     {
975 :     my $clust = [$x->[1]];
976 :     while ((@pegs > 0) && (abs($pegs[0]->[0] - $x->[0]) <= $distance))
977 :     {
978 :     $x = shift @pegs;
979 :     push(@$clust,$x->[1]);
980 :     }
981 :    
982 :     if (@$clust > 1)
983 :     {
984 :     push(@clusters,$clust);
985 :     }
986 :     }
987 :     }
988 :     return sort { @$b <=> @$a } @clusters;
989 :     }
990 :    
991 : overbeek 1.42 sub boundaries_of {
992 :     my($self,@args) = @_;
993 :    
994 :     my $fig = $self->{_fig};
995 : olson 1.43 return $fig->boundaries_of(@args);
996 : overbeek 1.42 }
997 :    
998 :    
999 : overbeek 1.7
1000 : overbeek 1.1 sub feature_location {
1001 :     my($self,$fid) = @_;
1002 :    
1003 :     my $fig = $self->{_fig};
1004 :     my $newG = $self->{_genome};
1005 :     my $newGdir = $self->{_orgdir};
1006 :    
1007 :     if (($fid =~ /^fig\|(\d+\.\d+)/) && ($1 eq $newG))
1008 :     {
1009 :     &load_tbl($self);
1010 :     if (my $x = $self->{_tbl}->{$fid})
1011 :     {
1012 :     return join(",",@{$x->[0]});
1013 :     }
1014 :     else
1015 :     {
1016 :     return undef;
1017 :     }
1018 :     }
1019 :     else
1020 :     {
1021 :     return scalar $fig->feature_location($fid);
1022 :     }
1023 :     }
1024 :    
1025 :     sub function_of {
1026 :     my($self,$fid) = @_;
1027 :    
1028 :     my $fig = $self->{_fig};
1029 :     my $newG = $self->{_genome};
1030 :     my $newGdir = $self->{_orgdir};
1031 :    
1032 :     if (($fid =~ /^fig\|(\d+\.\d+)/) && ($1 eq $newG))
1033 :     {
1034 :     &load_functions($self);
1035 :     return $self->{_functions}->{$fid};
1036 :     }
1037 :     else
1038 :     {
1039 :     return scalar $fig->function_of($fid);
1040 :     }
1041 :     }
1042 :    
1043 :     sub feature_aliases {
1044 :     my($self,$fid) = @_;
1045 :    
1046 :     my $fig = $self->{_fig};
1047 :     my $newG = $self->{_genome};
1048 :     my $newGdir = $self->{_orgdir};
1049 :    
1050 :     my @aliases;
1051 :     if (($fid =~ /^fig\|(\d+\.\d+)/) && ($1 eq $newG))
1052 :     {
1053 :     &load_tbl($self);
1054 :     if (my $x = $self->{_tbl}->{$fid})
1055 :     {
1056 :     @aliases = @{$x->[1]};
1057 :     }
1058 :     else
1059 :     {
1060 :     @aliases = ();
1061 :     }
1062 :     }
1063 :     else
1064 :     {
1065 :     @aliases = $fig->feature_aliases($fid);
1066 :     }
1067 :     return wantarray() ? @aliases : join(",",@aliases);
1068 :     }
1069 :    
1070 :     sub feature_annotations {
1071 :     my($self,$fid,$rawtime) = @_;
1072 :    
1073 :     my $fig = $self->{_fig};
1074 :     my $newG = $self->{_genome};
1075 :     my $newGdir = $self->{_orgdir};
1076 :    
1077 :     my @annotations;
1078 :     if (($fid =~ /^fig\|(\d+\.\d+)/) && ($1 eq $newG))
1079 :     {
1080 :     &load_ann($self);
1081 :     if (my $x = $self->{_ann}->{$fid})
1082 :     {
1083 :     @annotations = @{$x};
1084 :     }
1085 :     else
1086 :     {
1087 :     @annotations = ();
1088 :     }
1089 :    
1090 :     if ($rawtime)
1091 :     {
1092 :     return @annotations;
1093 :     }
1094 :     else
1095 :     {
1096 :     return map { $_->[1] = localtime($_->[1]); $_ } @annotations;
1097 :     }
1098 :     }
1099 :     else
1100 :     {
1101 :     return $fig->feature_annotations($fid);
1102 :     }
1103 :     }
1104 :    
1105 :     sub get_translation {
1106 :     my($self,$peg) = @_;
1107 :    
1108 :     my $fig = $self->{_fig};
1109 :     my $newG = $self->{_genome};
1110 :     my $newGdir = $self->{_orgdir};
1111 :    
1112 :     if (($peg =~ /^fig\|(\d+\.\d+)/) && ($1 eq $newG))
1113 :     {
1114 :     &load_pseq($self);
1115 :     return $self->{_pseq}->{$peg};
1116 :     }
1117 :     else
1118 :     {
1119 :     return $fig->get_translation($peg);
1120 :     }
1121 :     }
1122 :    
1123 : olson 1.4 sub translation_length
1124 :     {
1125 :     my($self, $peg) = @_;
1126 :    
1127 :     my $fig = $self->{_fig};
1128 :     my $newG = $self->{_genome};
1129 :     my $newGdir = $self->{_orgdir};
1130 :    
1131 :     if (($peg =~ /^fig\|(\d+\.\d+)/) && ($1 eq $newG))
1132 :     {
1133 :     &load_pseq($self);
1134 :     return length($self->{_pseq}->{$peg});
1135 :     }
1136 :     else
1137 :     {
1138 :     return $fig->translation_length($peg);
1139 :     }
1140 :     }
1141 :    
1142 :     sub translatable
1143 :     {
1144 :     my($self, $peg) = @_;
1145 :    
1146 :     my $fig = $self->{_fig};
1147 :     my $newG = $self->{_genome};
1148 :     my $newGdir = $self->{_orgdir};
1149 :    
1150 :     if (($peg =~ /^fig\|(\d+\.\d+)/) && ($1 eq $newG))
1151 :     {
1152 :     &load_pseq($self);
1153 :     return length($self->{_pseq}->{$peg}) > 0 ? 1 : 0;;
1154 :     }
1155 :     else
1156 :     {
1157 :     return $fig->translatable($peg);
1158 :     }
1159 :     }
1160 :    
1161 : gdpusch 1.48 sub pick_gene_boundaries {
1162 :     return &FIG::pick_gene_boundaries(@_);
1163 :     }
1164 :    
1165 :     sub call_start {
1166 :     return &FIG::call_start(@_);
1167 :     }
1168 :    
1169 : olson 1.4 sub is_real_feature
1170 :     {
1171 :     my($self, $fid) = @_;
1172 :    
1173 :     if ($self->is_virtual_feature($fid))
1174 :     {
1175 :     return 1;
1176 :     }
1177 :     else
1178 :     {
1179 :     return $self->{_fig}->is_real_feature($fid);
1180 :     }
1181 :    
1182 :     }
1183 :    
1184 : olson 1.10 sub pegs_of
1185 :     {
1186 :     my($self, $genome) = @_;
1187 :    
1188 :     my $fig = $self->{_fig};
1189 :     my $newG = $self->{_genome};
1190 :     my $newGdir = $self->{_orgdir};
1191 :    
1192 :     if ($genome ne $newG)
1193 :     {
1194 :     return $fig->pegs_of($genome);
1195 :     }
1196 :    
1197 :     $self->load_tbl();
1198 :     return grep { /\.peg\./ } keys %{$self->{_tbl}};
1199 :     }
1200 :    
1201 :    
1202 :     sub rnas_of
1203 :     {
1204 :     my($self, $genome) = @_;
1205 :    
1206 :     my $fig = $self->{_fig};
1207 :     my $newG = $self->{_genome};
1208 :     my $newGdir = $self->{_orgdir};
1209 :    
1210 :     if ($genome ne $newG)
1211 :     {
1212 :     return $fig->pegs_of($genome);
1213 :     }
1214 :    
1215 :     $self->load_tbl();
1216 :     return grep { /\.rna\./ } keys %{$self->{_tbl}};
1217 :     }
1218 :    
1219 : olson 1.4 sub is_virtual_feature
1220 :     {
1221 :     my($self, $peg) = @_;
1222 :    
1223 :     my $fig = $self->{_fig};
1224 :     my $newG = $self->{_genome};
1225 :     my $newGdir = $self->{_orgdir};
1226 :    
1227 :     if (($peg =~ /^fig\|(\d+\.\d+)/) && ($1 eq $newG))
1228 :     {
1229 :     &load_pseq($self);
1230 :     return 1;
1231 :     }
1232 :     else
1233 :     {
1234 :     return 0;
1235 :     }
1236 :     }
1237 :    
1238 : olson 1.35 ####################################################
1239 :     #
1240 :     # Following are some MG-RAST specific features. FIGV seems a good a place as any to include them.#
1241 :     #
1242 :     #
1243 :    
1244 :     =head3 taxa_to_seed_ids
1245 :    
1246 :     Given a prefix of a taxonomy string, return the list of metagenome fids that
1247 :     mapped to SEED organisms in that taxonomy.
1248 :    
1249 :     =cut
1250 :    
1251 :     sub taxa_to_seed_ids
1252 :     {
1253 :     my($self, $tax) = @_;
1254 :    
1255 :     my $btree_file = "$self->{_orgdir}/taxa_summary_by_blast.btree";
1256 :     my %btree;
1257 :    
1258 :     my $btree_tie = tie %btree, 'DB_File', $btree_file, O_RDONLY, 0666, $DB_BTREE;
1259 :    
1260 :     if (!$btree_tie)
1261 :     {
1262 :     warn "Cannot tie $btree_file: $!\n";
1263 :     return undef;
1264 :     }
1265 :    
1266 :     my $key = $tax;
1267 :     my ($res, $value);
1268 :     my @vals;
1269 :     my %seen;
1270 :    
1271 :     for ($res = $btree_tie->seq($key, $value, R_CURSOR);
1272 :     $res == 0 and $key =~ /^$tax/;
1273 :     $res = $btree_tie->seq($key, $value, R_NEXT))
1274 :     {
1275 :     print "$key\n";
1276 :     push(@vals, map { [$key, $_] } grep { not $seen{$_}++ } split(/$;/, $value));
1277 :     }
1278 :     return @vals;
1279 :     }
1280 :    
1281 : overbeek 1.1 sub load_pseq {
1282 :     my($self) = @_;
1283 :    
1284 :     if ($self->{_pseq}) { return };
1285 :    
1286 :     my $newGdir = $self->{_orgdir};
1287 :     my $pseq = {};
1288 :     if (open(FASTA,"<$newGdir/Features/peg/fasta"))
1289 :     {
1290 : olson 1.10 local $/ = "\n>";
1291 : overbeek 1.1 my $x;
1292 :     while (defined($x = <FASTA>))
1293 :     {
1294 : olson 1.21 chomp $x;
1295 : olson 1.10 if ($x =~ />?(\S+)[^\n]*\n(.*)/s)
1296 : overbeek 1.1 {
1297 :     my $peg = $1;
1298 :     my $seq = $2;
1299 :     $seq =~ s/\s//gs;
1300 :     $pseq->{$peg} = $seq;
1301 :     }
1302 :     }
1303 :     close(FASTA);
1304 :     }
1305 :     $self->{_pseq} = $pseq;
1306 :     }
1307 :    
1308 : olson 1.9 sub load_ss_data
1309 :     {
1310 :     my($self) = @_;
1311 :    
1312 :     return if defined($self->{_ss_bindings});
1313 :    
1314 :     my $fig = $self->{_fig};
1315 :     my $newG = $self->{_genome};
1316 :     my $newGdir = $self->{_orgdir};
1317 :    
1318 :     open(S, "<$newGdir/Subsystems/bindings") or die "Cannot open $newGdir/Subsystems/bindings: $!";
1319 :    
1320 :     my $peg_index;
1321 :     my $bindings;
1322 : olson 1.14 my $genome_index;
1323 :     my %seen;
1324 : olson 1.9 while (<S>)
1325 :     {
1326 :     chomp;
1327 :     my($sname, $role, $peg) = split(/\t/);
1328 : olson 1.14
1329 : olson 1.34 my($genome) = ($peg =~ /fig\|(\d+\.\d+)/);
1330 :    
1331 :     # my $genome = &FIG::genome_of($peg);
1332 : olson 1.9
1333 :     push(@{$bindings->{$sname}->{$role}}, $peg);
1334 :     push(@{$peg_index->{$peg}}, [$sname, $role]);
1335 : olson 1.14
1336 :     if (!$seen{$genome, $sname, $role})
1337 :     {
1338 : parrello 1.24 push(@{$genome_index->{$genome}}, [$sname, $role]);
1339 : olson 1.14 }
1340 :    
1341 : olson 1.9 }
1342 :     close(S);
1343 :    
1344 :     open(S, "<$newGdir/Subsystems/subsystems") or die "Cannot open $newGdir/Subsystems/subsystems: $!";
1345 :     my $variant;
1346 :     while (<S>)
1347 :     {
1348 :     chomp;
1349 :     my($sname, $v) = split(/\t/);
1350 :     $variant->{$sname} = $v;
1351 :     }
1352 :     close(S);
1353 :    
1354 :     $self->{_ss_bindings} = $bindings;
1355 :     $self->{_ss_variants} = $variant;
1356 :     $self->{_ss_peg_index} = $peg_index;
1357 : olson 1.14 $self->{_ss_genome_index} = $genome_index;
1358 : olson 1.9 }
1359 :    
1360 : overbeek 1.1 sub load_tbl {
1361 :     my($self) = @_;
1362 :    
1363 :     if ($self->{_tbl}) { return };
1364 :    
1365 :     my $newGdir = $self->{_orgdir};
1366 :     my $tbl = {};
1367 :     foreach my $x (`cat $newGdir/Features/*/tbl`)
1368 :     {
1369 :     if ($x =~ /^(\S+)\t(\S+)(\t(\S.*\S))?/)
1370 :     {
1371 :     my $fid = $1;
1372 :     my $loc = [split(/,/,$2)];
1373 : paczian 1.12 my $aliases = $4 ? [split(/\t/,$4)] : [];
1374 : overbeek 1.1 $tbl->{$fid} = [$loc,$aliases];
1375 :     }
1376 :     }
1377 :     $self->{_tbl} = $tbl;
1378 :     }
1379 :    
1380 :     sub load_functions {
1381 :     my($self) = @_;
1382 :    
1383 :     if ($self->{_functions}) { return };
1384 :    
1385 :     my $newG = $self->{_genome};
1386 :     my $newGdir = $self->{_orgdir};
1387 :    
1388 :     my $functions = {};
1389 :     foreach my $x (`cat $newGdir/*functions`)
1390 :     {
1391 :     if (($x =~ /^(fig\|(\d+\.\d+)\.\S+)\t(\S[^\t]*\S)/) && ($2 eq $newG))
1392 :     {
1393 :     $functions->{$1} = $3;
1394 :     }
1395 :     }
1396 :     $self->{_functions} = $functions;
1397 :     }
1398 :    
1399 : olson 1.16
1400 :     sub bbhs
1401 :     {
1402 :     my($self,$peg,$cutoff) = @_;
1403 :    
1404 :     my $fig = $self->{_fig};
1405 :     my $newG = $self->{_genome};
1406 :     my $newGdir = $self->{_orgdir};
1407 :    
1408 :     $self->load_bbh_index();
1409 :    
1410 :     $cutoff = 1.0e-10 unless defined($cutoff);
1411 :    
1412 :     my $tie = $self->{_bbh_tie};
1413 :    
1414 :     my @bbhs = $tie->get_dup($peg);
1415 :    
1416 :     # @bbhs now a list of comma-separated tuples (id2, score, bitscore)
1417 :     # print Dumper(\@bbhs);
1418 :    
1419 :     @bbhs = grep { $_->[1] < $cutoff } map { [split(/,/)] } @bbhs;
1420 :    
1421 :     if (not ($peg =~ /^fig\|(\d+\.\d+)/ && ($1 eq $newG)))
1422 :     {
1423 :     #
1424 :     # If this isn't one of our pegs, we retrieve the bbhs from the underlying
1425 :     # SEED and merge in the bbhs that we have here.
1426 :     #
1427 :    
1428 :     my @fbbhs = $fig->bbhs($peg, $cutoff);
1429 :     # print Dumper(\@fbbhs);
1430 :     push(@bbhs, @fbbhs);
1431 :     }
1432 :     return sort { $a->[1] <=> $b->[1] } @bbhs;
1433 :     }
1434 :    
1435 : olson 1.13 sub sims
1436 :     {
1437 : olson 1.27 my($self,$pegarg,$max,$maxP,$select, $max_expand, $filters) = @_;
1438 : olson 1.13
1439 :     my $fig = $self->{_fig};
1440 :     my $newG = $self->{_genome};
1441 :     my $newGdir = $self->{_orgdir};
1442 : parrello 1.24 $max = $max ? $max : 10000;
1443 :     $maxP = $maxP ? $maxP : 1.0e-5;
1444 : olson 1.13 $select = $select ? $select : "all";
1445 :    
1446 : olson 1.27 my $prefix = "_sims";
1447 :     my $flip_prefix = "_flips";
1448 :    
1449 :     if ($select eq 'raw')
1450 :     {
1451 :     $prefix = "_raw_sims";
1452 :     $flip_prefix = "_raw_flips";
1453 :     }
1454 :    
1455 :     #
1456 :     # Partition pegs into one set that is part of this organism
1457 :     # and another set that is not.
1458 :     #
1459 :     my @pegs = ref($pegarg) ? @$pegarg : ($pegarg);
1460 :     my(@part, @not_part, %part);
1461 :     my %sims;
1462 :    
1463 :     for my $peg (@pegs)
1464 :     {
1465 :     if ($peg =~ /^fig\|(\d+\.\d+)/ and $1 eq $newG)
1466 :     {
1467 :     push(@part, $peg);
1468 :     $part{$peg}++;
1469 :     }
1470 :     else
1471 :     {
1472 :     push(@not_part, $peg);
1473 :     $sims{$peg} = [];
1474 :     }
1475 :     }
1476 :    
1477 :     #
1478 :     # Retrieve a batch of the sims for the not-part pegs, and partition
1479 :     # into %sims.
1480 :     #
1481 :     for my $sim ($fig->sims(\@not_part, $max, $maxP, $select, $max_expand, $filters))
1482 :     {
1483 :     push(@{$sims{$sim->id1}}, $sim);
1484 :     }
1485 :    
1486 :     my @out;
1487 :     my $start_len;
1488 :    
1489 :     for my $peg (@pegs)
1490 : olson 1.13 {
1491 : olson 1.27 $start_len = @out;
1492 :     if (not $part{$peg})
1493 : overbeek 1.23 {
1494 : olson 1.27 my @flips = $self->retrieve_sims($peg, $flip_prefix, $maxP, $select);
1495 :    
1496 :     @flips = sort { $b->bsc <=> $a->bsc } @flips;
1497 :    
1498 :     # my @old = $fig->sims($peg,$max,$maxP,$select);
1499 :    
1500 :     my @old = sort { $b->bsc <=> $a->bsc } @{$sims{$peg}};
1501 :    
1502 :     # my @merged = ();
1503 :     my $i1 = 0;
1504 :     my $i2 = 0;
1505 :     while (($i1 < @flips) || ($i2 < @old))
1506 : overbeek 1.23 {
1507 : olson 1.27 if (($i1 == @flips) || (($i2 < @old) && ($flips[$i1]->[11] < $old[$i2]->[11])))
1508 :     {
1509 :     # push(@merged,$old[$i2]);
1510 :     push(@out,$old[$i2]);
1511 :     $i2++;
1512 :     }
1513 :     else
1514 :     {
1515 :     # push(@merged,$flips[$i1]);
1516 :     push(@out,$flips[$i1]);
1517 :     $i1++;
1518 :     }
1519 : overbeek 1.23 }
1520 :     }
1521 : olson 1.27 else
1522 :     {
1523 :     my @sims = $self->retrieve_sims($peg, $prefix, $maxP, $select);
1524 :     push(@out, @sims);
1525 :     }
1526 :    
1527 :     if (@out - $start_len > $max)
1528 :     {
1529 :     $#out = $start_len + $max - 1;
1530 :     }
1531 : olson 1.13 }
1532 : olson 1.27
1533 :     return @out;
1534 :     }
1535 : olson 1.13
1536 : olson 1.27 sub retrieve_sims
1537 :     {
1538 :     my($self, $peg, $prefix, $maxP, $select) = @_;
1539 :    
1540 : olson 1.13 $self->load_sims_index();
1541 :    
1542 : olson 1.27 my $info = $self->{"${prefix}_index"}->{$peg};
1543 : olson 1.13
1544 :     if ($info !~ /^(\d+),(\d+)$/)
1545 :     {
1546 :     return;
1547 :     }
1548 :     my($seek, $len) = ($1, $2);
1549 :    
1550 : olson 1.27 my $sims_txt = &FIG::read_block($self->{"${prefix}_fh"}, $seek, $len);
1551 : olson 1.13
1552 :     my @sims = map { bless $_, 'Sim' } grep { $_->[10] <= $maxP } map { [split(/\t/)] } @$sims_txt;
1553 :    
1554 :     if ($select =~ /^figx?$/)
1555 :     {
1556 :     @sims = grep { $_->[1] =~ /^fig/ } @sims;
1557 :     }
1558 :    
1559 :     return @sims;
1560 :     }
1561 :    
1562 :     sub sims_old {
1563 : overbeek 1.1 my($self,$peg,$max,$maxP,$select) = @_;
1564 :    
1565 :     my $fig = $self->{_fig};
1566 :     my $newG = $self->{_genome};
1567 :     my $newGdir = $self->{_orgdir};
1568 : parrello 1.24 $max = $max ? $max : 10000;
1569 :     $maxP = $maxP ? $maxP : 1.0e-5;
1570 : overbeek 1.1 $select = $select ? $select : "all";
1571 :    
1572 :     if (($peg =~ /^fig\|(\d+\.\d+)/) && ($1 eq $newG))
1573 :     {
1574 :     &load_sims($self);
1575 :     my @sims = ();
1576 :     my $raw_sims = $self->{_sims}->{$peg};
1577 :     if ($raw_sims)
1578 :     {
1579 :     foreach my $sim ( grep { $_->[10] <= $maxP } @$raw_sims )
1580 :     {
1581 :     my $id2 = $sim->id2;
1582 :     my $id1 = $sim->id1;
1583 :     my @relevant = ();
1584 :    
1585 :     my @maps_to = $fig->mapped_prot_ids( $id2 );
1586 :     my $ref_len = $maps_to[0]->[1];
1587 :    
1588 :     @maps_to = grep { $_->[0] !~ /^xxx\d+/ } @maps_to;
1589 :    
1590 :     if ( $select =~ /^figx?$/ ) # Only fig
1591 :     {
1592 :     @relevant = grep { $_->[0] =~ /^fig/ } @maps_to;
1593 :     }
1594 :     else # All
1595 :     {
1596 :     @relevant = @maps_to;
1597 :     }
1598 :    
1599 :     my $seen = {};
1600 :     foreach my $x ( @relevant )
1601 :     {
1602 :     my ( $x_id, $x_ln ) = @$x;
1603 :    
1604 :     next if $seen->{$x_id};
1605 :     $seen->{$x_id} = 1;
1606 :    
1607 :     my $delta2 = $ref_len - $x_ln; # Coordinate shift
1608 :     my $sim1 = [ @$sim ]; # Make a copy
1609 :     $sim1->[1] = $x_id;
1610 :     $sim1->[8] -= $delta2;
1611 :     $sim1->[9] -= $delta2;
1612 :     bless( $sim1, "Sim" );
1613 :     push( @sims, $sim1 );
1614 :     }
1615 :     }
1616 :     }
1617 :    
1618 :     if (@sims > $max) { $#sims = $max-1; }
1619 :     return @sims;
1620 :     }
1621 :     else
1622 :     {
1623 :     return $fig->sims($peg,$max,$maxP,$select);
1624 :     }
1625 :     }
1626 :    
1627 : olson 1.27
1628 :     sub coupled_to
1629 :     {
1630 :     my($self,$peg) = @_;
1631 :    
1632 :     my $fig = $self->{_fig};
1633 :     my $newG = $self->{_genome};
1634 :     my $newGdir = $self->{_orgdir};
1635 :    
1636 :     if ($peg =~ /^fig\|$newG\.peg/)
1637 :     {
1638 :    
1639 :     $self->load_coupling_index();
1640 :    
1641 :     my $tie = $self->{_pch_tie};
1642 : olson 1.30
1643 :     if (not defined($tie))
1644 :     {
1645 :     return;
1646 :     }
1647 : olson 1.27
1648 :     my @dat = $tie->get_dup($peg);
1649 :    
1650 :     return map { [split(/$;/, $_)] } @dat;
1651 :     }
1652 :     else
1653 :     {
1654 :     return $fig->coupled_to($peg);
1655 :     }
1656 :    
1657 :     }
1658 :    
1659 :     sub coupling_evidence
1660 :     {
1661 :     my($self,$peg1, $peg2) = @_;
1662 :    
1663 :     my $fig = $self->{_fig};
1664 :     my $newG = $self->{_genome};
1665 :     my $newGdir = $self->{_orgdir};
1666 :    
1667 :     if ($peg1 =~ /^fig\|$newG\.peg/)
1668 :     {
1669 :     $self->load_coupling_index();
1670 :    
1671 :     my $tie = $self->{_pch_ev_tie};
1672 :    
1673 : olson 1.30 if (not defined($tie))
1674 :     {
1675 :     return;
1676 :     }
1677 :    
1678 : olson 1.27 my @dat = $tie->get_dup("$peg1$;$peg2");
1679 :    
1680 :     my @a;
1681 :     return map { @a = split(/$;/, $_); [@a[0,1,4]] } @dat;
1682 :     }
1683 :     else
1684 :     {
1685 :     return $fig->coupling_evidence($peg1, $peg2);
1686 :     }
1687 :    
1688 :     }
1689 :    
1690 :     sub coupling_and_evidence
1691 :     {
1692 :     my($self,$peg1) = @_;
1693 :    
1694 :     my $fig = $self->{_fig};
1695 :     my $newG = $self->{_genome};
1696 :     my $newGdir = $self->{_orgdir};
1697 :    
1698 :     if ($peg1 =~ /^fig\|$newG\.peg/)
1699 :     {
1700 :     $self->load_coupling_index();
1701 :    
1702 :     my $tie = $self->{_pch_tie};
1703 :     my $evtie = $self->{_pch_ev_tie};
1704 :    
1705 : olson 1.30 if (not(defined($tie) and defined($evtie)))
1706 :     {
1707 :     return;
1708 :     }
1709 :    
1710 : olson 1.27 my @out;
1711 :     my @coupled_to = $tie->get_dup($peg1);
1712 :     for my $ent (@coupled_to)
1713 :     {
1714 :     my ($peg2, $score) = split(/$;/, $ent);
1715 :    
1716 :     my @ev = $evtie->get_dup("$peg1$;$peg2");
1717 :    
1718 :     my $l = [];
1719 :     for my $event (@ev)
1720 :     {
1721 :     my($peg3, $peg4, $iden3, $iden4, $rep) = split(/$;/, $event);
1722 :     push(@$l, [$peg3, $peg4]);
1723 :     }
1724 :     push(@out, [$score, $peg2, $l]);
1725 :     }
1726 :     return @out;
1727 :     }
1728 :     else
1729 :     {
1730 :     return $fig->coupling_and_evidence($peg1);
1731 :     }
1732 :    
1733 :     }
1734 :    
1735 : olson 1.45 sub in_pch_pin_with
1736 :     {
1737 :     my($self, $peg1, $diverse) = @_;
1738 :    
1739 :     my @all = $self->in_pch_pin_with_and_evidence($peg1);
1740 :    
1741 :     if ($diverse)
1742 :     {
1743 :     return map { $_->[0] } grep { $_->[1] == 1 } @all;
1744 :     }
1745 :     else
1746 :     {
1747 :     return map { $_->[0] } @all;
1748 :     }
1749 :     }
1750 :    
1751 :     sub in_pch_pin_with_and_evidence
1752 :     {
1753 :     my($self,$peg1) = @_;
1754 :    
1755 :     my $fig = $self->{_fig};
1756 :     my $newG = $self->{_genome};
1757 :     my $newGdir = $self->{_orgdir};
1758 :    
1759 :     if ($peg1 !~ /^fig\|$newG\.peg/)
1760 :     {
1761 :     return $fig->in_pch_pin_with_and_evidence($peg1);
1762 :     }
1763 :    
1764 :     $self->load_coupling_index();
1765 :    
1766 :     my $evtie = $self->{_pch_ev_tie};
1767 :    
1768 :     my($key, $value, $st);
1769 :    
1770 :     my %max;
1771 :    
1772 :     $key = "$peg1$;";
1773 :     my $qkey = quotemeta($key);
1774 :     for ($st = $evtie->seq($key, $value, R_CURSOR); $st == 0 and $key =~ /^$qkey/; $st = $evtie->seq($key, $value, R_NEXT))
1775 :     {
1776 :     # print "key=$key value=$value\n";
1777 :    
1778 :     my($peg3, $peg4, $iden3, $iden4, $rep) = split(/$;/, $value);
1779 :    
1780 :     if (exists($max{$peg3}))
1781 :     {
1782 :     $max{$peg3} = $rep if $rep > $max{$peg3};
1783 :     }
1784 :     else
1785 :     {
1786 :     $max{$peg3} = $rep;
1787 :     }
1788 :     }
1789 :    
1790 :     return map { [$_, $max{$_}] } keys %max;
1791 :     }
1792 :    
1793 :    
1794 : olson 1.27 sub load_coupling_index
1795 :     {
1796 :     my($self) = @_;
1797 :    
1798 :     return if defined($self->{_pch_index});
1799 :    
1800 :     my $fig = $self->{_fig};
1801 :     my $newG = $self->{_genome};
1802 :     my $newGdir = $self->{_orgdir};
1803 :    
1804 :     my $pch_btree = "$newGdir/pchs.btree";
1805 :     my $ev_btree = "$newGdir/pchs.evidence.btree";
1806 :    
1807 :     my $pch_index = {};
1808 :     my $ev_index = {};
1809 :    
1810 :     my $tied = tie %$pch_index, 'DB_File', $pch_btree, O_RDONLY, 0666, $DB_BTREE;
1811 :     my $evtied = tie %$ev_index, 'DB_File', $ev_btree, O_RDONLY, 0666, $DB_BTREE;
1812 :    
1813 :     #
1814 :     # Set these even if failed so we don't keep trying to open and failing.
1815 :     #
1816 :     $self->{_pch_index} = $pch_index;
1817 :     $self->{_pch_tie} = $tied;
1818 :     $self->{_pch_ev_index} = $ev_index;
1819 :     $self->{_pch_ev_tie} = $evtied;
1820 :    
1821 :     if (not $tied)
1822 :     {
1823 :     warn "Cannot tie pch index $pch_btree: $!\n";
1824 :     }
1825 :     }
1826 :    
1827 : olson 1.16 sub load_bbh_index
1828 :     {
1829 :     my($self) = @_;
1830 :    
1831 :     return if defined($self->{_bbh_index});
1832 :    
1833 :     my $fig = $self->{_fig};
1834 :     my $newG = $self->{_genome};
1835 :     my $newGdir = $self->{_orgdir};
1836 :    
1837 :     my $bbh_file = "$newGdir/bbhs";
1838 :     my $bbh_index_file = "$bbh_file.index";
1839 :     my $bbh_index = {};
1840 :    
1841 :     my $tied = tie %$bbh_index, 'DB_File', $bbh_index_file, O_RDONLY, 0666, $DB_BTREE;
1842 :    
1843 :     #
1844 :     # Set these even if failed so we don't keep trying to open and failing.
1845 :     #
1846 :     $self->{_bbh_index} = $bbh_index;
1847 :     $self->{_bbh_tie} = $tied;
1848 :    
1849 :     if (not $tied)
1850 :     {
1851 :     warn "Cannot tie bbh index $bbh_index_file: $!\n";
1852 :     }
1853 :     }
1854 :    
1855 : olson 1.36 sub load_contigs_index
1856 :     {
1857 :     my($self) = @_;
1858 :    
1859 :     return if defined($self->{_contigs_index});
1860 :    
1861 :     my $fig = $self->{_fig};
1862 :     my $newG = $self->{_genome};
1863 :     my $newGdir = $self->{_orgdir};
1864 :    
1865 :     my $contig_index_file = "$newGdir/contigs.btree";
1866 :     my $contig_index = {};
1867 : olson 1.37 my $contig_len_index_file = "$newGdir/contig_len.btree";
1868 :     my $contig_len_index = {};
1869 : olson 1.36
1870 :     my $tied = tie %$contig_index, 'DB_File', $contig_index_file, O_RDONLY, 0666, $DB_BTREE;
1871 : olson 1.37 if (not $tied)
1872 :     {
1873 : olson 1.47 # warn "Cannot tie contig index $contig_index_file: $!\n";
1874 : olson 1.37 }
1875 :    
1876 :     my $ltied = tie %$contig_len_index, 'DB_File', $contig_len_index_file, O_RDONLY, 0666, $DB_BTREE;
1877 :     if (not $ltied)
1878 :     {
1879 : olson 1.47 # warn "Cannot tie contig length index $contig_len_index_file: $!\n";
1880 : olson 1.37 }
1881 : olson 1.36
1882 :     #
1883 :     # Set these even if failed so we don't keep trying to open and failing.
1884 :     #
1885 : olson 1.37 $self->{_contigs_index} = $contig_index;
1886 :     $self->{_contigs_tie} = $tied;
1887 :     $self->{_contig_len_index} = $contig_len_index;
1888 :     $self->{_contig_len_tie} = $ltied;
1889 : olson 1.36
1890 :     }
1891 :    
1892 : olson 1.13 sub load_sims_index
1893 :     {
1894 :     my($self) = @_;
1895 :    
1896 :     return if defined($self->{_sims_index});
1897 :    
1898 :     my $fig = $self->{_fig};
1899 :     my $newG = $self->{_genome};
1900 :     my $newGdir = $self->{_orgdir};
1901 :    
1902 :     my $sims_file = "$newGdir/expanded_similarities";
1903 :     my $sims_index_file = "$sims_file.index";
1904 :     my $sims_index = {};
1905 :    
1906 : olson 1.27 my $flips_file = "$newGdir/expanded_similarities.flips";
1907 :     my $flips_index_file = "$flips_file.index";
1908 :     my $flips_index = {};
1909 :    
1910 :     my $raw_sims_file = "$newGdir/similarities";
1911 :     my $raw_sims_index_file = "$raw_sims_file.index";
1912 :     my $raw_sims_index = {};
1913 :    
1914 :     my $raw_flips_file = "$newGdir/similarities.flips";
1915 :     my $raw_flips_index_file = "$raw_flips_file.index";
1916 :     my $raw_flips_index = {};
1917 :    
1918 :     $self->open_sims($sims_index, $sims_file, $sims_index_file, "_sims");
1919 :     $self->open_sims($flips_index, $flips_file, $flips_index_file, "_flips");
1920 :     $self->open_sims($raw_sims_index, $raw_sims_file, $raw_sims_index_file, "_raw_sims");
1921 :     $self->open_sims($raw_flips_index, $raw_flips_file, $raw_flips_index_file, "_raw_flips");
1922 :     }
1923 :    
1924 :     #
1925 :     # Open a sims file, tie it to a hash, and store the info into the $self obj.
1926 :     #
1927 :     sub open_sims
1928 :     {
1929 :     my($self, $hash, $sims_file, $index_file, $prefix) = @_;
1930 :    
1931 :     my $tied = tie %$hash, 'DB_File', $index_file, O_RDONLY, 0666, $DB_BTREE;
1932 : olson 1.13
1933 :     #
1934 :     # Set these even if failed so we don't keep trying to open and failing.
1935 :     #
1936 : olson 1.27 $self->{"${prefix}_index"} = $hash;
1937 :     $self->{"${prefix}_tie"} = $tied;
1938 : olson 1.13
1939 :     if (not $tied)
1940 :     {
1941 : olson 1.27 warn "Cannot tie sims index $index_file: $!\n";
1942 : olson 1.13 }
1943 :    
1944 :     #
1945 :     # open the sims file as well.
1946 :     #
1947 :    
1948 : olson 1.27 $self->{"${prefix}_fh"} = new FileHandle("<$sims_file");
1949 : olson 1.13
1950 : olson 1.27 if (!$self->{"${prefix}_fh"})
1951 : olson 1.13 {
1952 :     warn "Cannot open sims file $sims_file: $!\n";
1953 :     }
1954 :     }
1955 :    
1956 : overbeek 1.1 sub load_sims {
1957 :     my($self) = @_;
1958 :    
1959 :     if ($self->{_sims}) { return };
1960 :    
1961 :     my $newGdir = $self->{_orgdir};
1962 :    
1963 :     my $sims = {};
1964 :     foreach my $x (`cat $newGdir/similarities`)
1965 :     {
1966 : olson 1.46 chomp $x;
1967 : overbeek 1.1 if ($x =~ /^(\S+)/)
1968 :     {
1969 :     push(@{$sims->{$1}},bless([split(/\t/,$x)],'Sim'));
1970 :     }
1971 :     }
1972 :     $self->{_sims} = $sims;
1973 :     }
1974 :    
1975 :     sub get_attributes {
1976 : olson 1.44 my($self, $id, $attr) = @_;
1977 : overbeek 1.1
1978 :     my $fig = $self->{_fig};
1979 :     my $newG = $self->{_genome};
1980 :     my $newGdir = $self->{_orgdir};
1981 :    
1982 : olson 1.44 if ((($id =~ /^fig\|(\d+\.\d+)\.peg\.\d+$/) and ($1 eq $newG)) or
1983 :     (($id =~ /^(\d+\.\d+)/ and $1 eq $newG)))
1984 : overbeek 1.1 {
1985 :     &load_attr($self);
1986 : olson 1.44
1987 :     my $t = $self->{_attr_id_tie};
1988 :     if (!$t)
1989 : overbeek 1.1 {
1990 : olson 1.44 #
1991 :     # Tied index not present, bail back to simple attributes.
1992 :     #
1993 :    
1994 :     my $l = $self->{_attr_id}->{$id};
1995 :     return $l ? @$l : ();
1996 : overbeek 1.1 }
1997 : olson 1.44
1998 :     my @v = $t->get_dup($id);
1999 :    
2000 :     my @out;
2001 :     for my $v (@v)
2002 : overbeek 1.1 {
2003 : olson 1.44 my @a = split(/$;/, $v);
2004 :    
2005 :     if (!defined($attr) or $attr eq $a[1])
2006 :     {
2007 :     push(@out, [@a]);
2008 :     }
2009 : overbeek 1.1 }
2010 : olson 1.44 return @out;
2011 : overbeek 1.1 }
2012 :     else
2013 :     {
2014 : olson 1.44 my @f = $fig->get_attributes($id, $attr);
2015 :     if (!defined($id) and defined(my $t = $self->{_attr_key_tie}))
2016 :     {
2017 :     #
2018 :     # lookup locally for $attr matches and merge with other output.
2019 :     #
2020 :    
2021 :     my @mine = $t->get_dup($attr);
2022 :     @mine = map { [split(/$;/, $_)] } @mine;
2023 :    
2024 :     return (@mine, @f);
2025 :     }
2026 :     else
2027 :     {
2028 :     return @f;
2029 :     }
2030 : overbeek 1.1 }
2031 :     }
2032 :    
2033 :     sub load_attr {
2034 :     my($self) = @_;
2035 :    
2036 : olson 1.44 if ($self->{_attr_id}) { return };
2037 : overbeek 1.1
2038 :     my $newGdir = $self->{_orgdir};
2039 : olson 1.44
2040 :     my $id = {};
2041 :     my $key = {};
2042 :    
2043 :     $self->{_attr_id_tie} = tie %$id, 'DB_File', "$newGdir/attr_id.btree", O_RDONLY, 0, $DB_BTREE;
2044 :     $self->{_attr_key_tie} = tie %$key, 'DB_File', "$newGdir/attr_key.btree", O_RDONLY, 0, $DB_BTREE;
2045 :    
2046 :     $self->{_attr_id} = $id;
2047 :     $self->{_attr_key} = $key;
2048 :    
2049 :     #
2050 :     # If the tie didn't work for ids, at least load up the evidence codes.
2051 :     #
2052 :    
2053 :     if (!$self->{_attr_id_tie})
2054 : overbeek 1.1 {
2055 : olson 1.44 foreach my $x (`cat $newGdir/evidence.codes`)
2056 : overbeek 1.1 {
2057 : olson 1.44 if ($x =~ /^(\S+)\t(\S+)/)
2058 :     {
2059 :     push(@{$id->{$1}},[$1,"evidence_code",$2,""]);
2060 :     }
2061 : overbeek 1.1 }
2062 :     }
2063 : olson 1.44
2064 : overbeek 1.1 }
2065 :    
2066 :     sub load_ann {
2067 :     my($self) = @_;
2068 :    
2069 :     if ($self->{_ann}) { return };
2070 :    
2071 :     my $newGdir = $self->{_orgdir};
2072 :     my $ann = {};
2073 :     if (open(ANN,"<$newGdir/annotations"))
2074 :     {
2075 :     $/ = "\n//\n";
2076 :     while (defined(my $x = <ANN>))
2077 :     {
2078 :     chomp $x;
2079 :     if ($x =~ /^(\S+)\n([^\n]+)\n([^\n]+)\n(.*)/s)
2080 :     {
2081 :     push(@{$ann->{$1}},[$1,$2,$3,"$4\n"]);
2082 :     }
2083 :     }
2084 :     $/ = "\n";
2085 :     close(ANN);
2086 :     }
2087 :     $self->{_ann} = $ann;
2088 :     }
2089 : overbeek 1.6
2090 :     sub taxonomy_of {
2091 :     my($self,$genome) = @_;
2092 :     my $fig = $self->{_fig};
2093 :     my $newG = $self->{_genome};
2094 :     my $newGdir = $self->{_orgdir};
2095 :    
2096 :     if ($newG ne $genome)
2097 :     {
2098 :     return $fig->taxonomy_of($genome);
2099 :     }
2100 :    
2101 :     my $tax;
2102 :     if (open(TAX,"<$newGdir/TAXONOMY") && ($tax = <TAX>))
2103 :     {
2104 :     chop $tax;
2105 :     return $tax;
2106 :     }
2107 :     else
2108 :     {
2109 :     return "unknown";
2110 :     }
2111 :     }
2112 :    
2113 :     sub build_tree_of_complete {
2114 :     my($self,$min_for_label) = @_;
2115 :     return $self->build_tree_of_all($min_for_label, "complete");
2116 :     }
2117 :    
2118 :     sub build_tree_of_all {
2119 :     my($self, $min_for_label, $complete)=@_;
2120 :     my(@last,@tax,$i,$prefix,$lev,$genome,$tax);
2121 :    
2122 :     $min_for_label = $min_for_label ? $min_for_label : 10;
2123 :     open(TMP,">/tmp/tree$$") || die "could not open /tmp/tree$$";
2124 :     print TMP "1. root\n";
2125 :    
2126 :     @last = ();
2127 :    
2128 :    
2129 :     foreach $genome (grep { ! $self->is_environmental($_) } $self->sort_genomes_by_taxonomy($self->genomes($complete)))
2130 :     {
2131 :     $tax = $self->taxonomy_of($genome);
2132 :     @tax = split(/\s*;\s*/,$tax);
2133 :     push(@tax,$genome);
2134 :     for ($i=0; ((@last > $i) && (@tax > $i) && ($last[$i] eq $tax[$i])); $i++) {}
2135 :     while ($i < @tax)
2136 :     {
2137 :     $lev = $i+2;
2138 :     $prefix = " " x (4 * ($lev-1));
2139 :     print TMP "$prefix$lev\. $tax[$i]\n";
2140 :     $i++;
2141 :     }
2142 :     @last = @tax;
2143 :     }
2144 :     close(TMP);
2145 :     my $tree = &tree_utilities::build_tree_from_outline("/tmp/tree$$");
2146 :     $tree->[0] = 'All';
2147 :     &FIG::limit_labels($tree,$min_for_label);
2148 :     unlink("/tmp/tree$$");
2149 :     return ($tree,&tree_utilities::tips_of_tree($tree));
2150 :     }
2151 :    
2152 :     sub sort_genomes_by_taxonomy {
2153 :     my($self,@genomes) = @_;
2154 :    
2155 :     return map { $_->[0] }
2156 :     sort { $a->[1] cmp $b->[1] }
2157 :     map { [$_,$self->taxonomy_of($_)] }
2158 :     @genomes;
2159 :     }
2160 : overbeek 1.1
2161 : overbeek 1.6 sub taxonomic_groups_of_complete {
2162 :     my($self,$min_for_labels) = @_;
2163 :    
2164 :     my($tree,undef) = $self->build_tree_of_complete($min_for_labels);
2165 :     return &FIG::taxonomic_groups($tree);
2166 :     }
2167 :    
2168 : olson 1.31 =head2 Search Database
2169 :    
2170 :     Searches the database for objects that match the query string in some way.
2171 :    
2172 :     Returns a list of results if the query is ambiguous or an unique identifier
2173 :     otherwise.
2174 :    
2175 :     =cut
2176 :    
2177 :     sub search_database {
2178 :     # get parameters
2179 :     my ($self, $query, $options) = @_;
2180 :    
2181 :     my $fig = $self->{_fig};
2182 :     my $newG = $self->{_genome};
2183 :     my $newGdir = $self->{_orgdir};
2184 :    
2185 :     #
2186 :     # Check for local ids and genome.
2187 :     #
2188 :    
2189 :    
2190 :     if ($query eq $newG or lc($query) eq lc($self->genus_species($newG)))
2191 :     {
2192 :     return { type => 'organism', result => $newG };
2193 :     }
2194 :    
2195 :     if ($query =~ /^fig\|(\d+\.\d+)\.peg/ and $1 eq $newG)
2196 :     {
2197 :     return { type => 'feature', result => $query };
2198 :     }
2199 :    
2200 :     #
2201 :     # Match on functions
2202 :     #
2203 :    
2204 :     &load_functions($self);
2205 :    
2206 :     my @protlist;
2207 :     my $fn = $self->{_functions};
2208 :     for my $id (keys %$fn)
2209 :     {
2210 :     if ($fn->{$id} =~ /$query/i)
2211 :     {
2212 :     push @protlist, [$id, $fn->{$id}, $newG];
2213 :     }
2214 :     }
2215 :    
2216 :     #
2217 :     # Pull FIG lookups.
2218 :     #
2219 :    
2220 :     my $res = $fig->search_database($query, $options);
2221 :    
2222 :     #
2223 :     # And integrate our protein list
2224 :     #
2225 :    
2226 :     my $res_prot;
2227 :     if (ref($res) eq 'ARRAY')
2228 :     {
2229 :     for my $ent (@$res)
2230 :     {
2231 :     if ($ent->{type} eq 'proteins')
2232 :     {
2233 :     $res_prot = $ent->{result};
2234 :     }
2235 :     }
2236 :     if (!$res_prot)
2237 :     {
2238 :     $res_prot = {type => 'proteins', result => \@protlist}
2239 :     }
2240 :     else
2241 :     {
2242 :     push(@$res_prot, @protlist);
2243 :     }
2244 :     }
2245 :    
2246 :     return $res;
2247 :    
2248 :     }
2249 :    
2250 :    
2251 : paczian 1.29 =head3 model_directory
2252 : paczian 1.26
2253 : paczian 1.29 C<< FIG->model_directory($organism); >>
2254 : paczian 1.26
2255 : paczian 1.29 Returns the model directory of an organism.
2256 : paczian 1.26
2257 :     =over 4
2258 :    
2259 :     =item $organism
2260 :    
2261 :     The seed-taxonomy id of the organism, e.g. 83333.1. If you pass the
2262 :     string 'All', you will get the directory for the global model.
2263 :    
2264 :     =back
2265 :    
2266 :     =cut
2267 :    
2268 : paczian 1.29 sub model_directory {
2269 : paczian 1.26 my ($self, $organism) = @_;
2270 :    
2271 : olson 1.27 my $directory;
2272 : paczian 1.26
2273 : olson 1.27 if ($organism eq $self->{_genome})
2274 :     {
2275 :     $directory = $self->{_orgdir} . "/Models";
2276 :     $directory .= "/$organism" if defined($organism);
2277 : formsma 1.38 }
2278 :     elsif(!defined $organism && defined $self->{_genome})
2279 :     {
2280 :     $directory = $self->{_orgdir} . "/Models";
2281 :     }
2282 :     else {
2283 : paczian 1.29 $directory = $self->{_fig}->model_directory($organism);
2284 : paczian 1.26 }
2285 :    
2286 :     return $directory;
2287 :     }
2288 :    
2289 : overbeek 1.1 1;

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3