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

Annotation of /FigKernelPackages/FIGV.pm

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1.7 - (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 : overbeek 1.1
30 :     sub new {
31 :     my($class,$org_dir,$low_level) = @_;
32 :    
33 :     my $fig;
34 :     if ($low_level && ($low_level =~ /sprout/i))
35 :     {
36 :     $fig = new SproutFIG($FIG_Config::sproutDB, $FIG_Config::sproutData);
37 :     }
38 :     else
39 :     {
40 :     $fig = new FIG;
41 :     }
42 :    
43 :     my $self = {};
44 :     $self->{_fig} = $fig;
45 :     $self->{_orgdir} = $org_dir;
46 :    
47 :     $org_dir =~ /(\d+\.\d+)$/;
48 :     $self->{_genome} = $1;
49 :     return bless $self, $class;
50 :     }
51 :    
52 : olson 1.4 #
53 :     # Redirect any method invocations that we don't implement out to the
54 :     # underlying FIG object.
55 :     #
56 :     sub AUTOLOAD
57 :     {
58 :     my($self, @args) = @_;
59 :    
60 : overbeek 1.6 if (ref($self) ne "FIGV") {
61 :     confess "BAD FIGV object passed to AUTOLOAD";
62 :     }
63 :    
64 : olson 1.4 my $meth = $AUTOLOAD;
65 :     $meth =~ s/.*:://;
66 :    
67 :     my $fig = $self->{_fig};
68 :     my $args = Dumper(\@args);
69 :     if (wantarray)
70 :     {
71 :     my @res = $fig->$meth(@args);
72 : olson 1.5 # warn "FIGV invoke $meth($args) returns\n", Dumper(\@res);
73 : olson 1.4 return @res;
74 :     }
75 :     else
76 :     {
77 :     my $res = $fig->$meth(@args);
78 : olson 1.5 # warn "FIGV invoke $meth($args) returns\n", Dumper($res);
79 : olson 1.4 return $res;
80 :     }
81 :     }
82 :    
83 : olson 1.5 #
84 :     # To retrieve a subsystem in FIGV, we create the subsystem as normal via $fig->get_subsystem,
85 :     # then insert the row for the virtual org dir we are processing.
86 :     #
87 :    
88 :     sub get_subsystem
89 :     {
90 :     my($self,$ssa) = @_;
91 :    
92 :     my $fig = $self->{_fig};
93 :     my $newG = $self->{_genome};
94 :     my $newGdir = $self->{_orgdir};
95 :    
96 :     my $ss = $fig->get_subsystem($ssa);
97 :     return undef unless $ss;
98 :    
99 :     return $ss;
100 :     }
101 : olson 1.4
102 : overbeek 1.1 sub genomes {
103 :     my($self,$complete) = @_;
104 :     my $fig = $self->{_fig};
105 :    
106 :     return ($self->{_genome},$fig->genomes($complete));
107 :     }
108 :    
109 :     sub genus_species {
110 :     my($self,$genome) = @_;
111 :    
112 :     my $fig = $self->{_fig};
113 :     my $newG = $self->{_genome};
114 :     my $newGdir = $self->{_orgdir};
115 :    
116 :     if (($genome eq $newG) && open(GENOME,"<$newGdir/GENOME"))
117 :     {
118 :     my $x = <GENOME>;
119 :     close(GENOME);
120 :     chop $x;
121 :     return $x;
122 :     }
123 :     else
124 :     {
125 :     return $fig->genus_species($genome);
126 :     }
127 :     }
128 :    
129 : overbeek 1.7 sub get_genome_assignment_data {
130 :     my($self,$genome) = @_;
131 :    
132 :     my $fig = $self->{_fig};
133 :     my $newG = $self->{_genome};
134 :     my $newGdir = $self->{_orgdir};
135 :    
136 :     if ($genome eq $newG)
137 :     {
138 :     my %assign = map { ( $_ =~ /^(fig\|\d+\.\d+\.peg\.\d+)\t(\S.*\S)/) ? ($1 => $2) : () }
139 :     `cat $newGdir/proposed*functions`;
140 :     return [map { [$_,$assign{$_}] } sort { &FIG::by_fig_id($a,$b) } keys(%assign)];
141 :     }
142 :     else
143 :     {
144 :     return $fig->get_genome_assignment_data($genome);
145 :     }
146 :     }
147 :    
148 :     sub org_of {
149 :     my($self,$peg) = @_;
150 :    
151 :     if ($peg =~ /^fig\|(\d+\.\d+)\.peg\.\d+/)
152 :     {
153 :     return $self->genus_species($1);
154 :     }
155 :     return "";
156 :     }
157 :    
158 :     sub get_genome_subsystem_data {
159 :     my($self,$genome) = @_;
160 :    
161 :     my $fig = $self->{_fig};
162 :     my $newG = $self->{_genome};
163 :     my $newGdir = $self->{_orgdir};
164 :    
165 :     if ($genome eq $newG)
166 :     {
167 :     return [map { ($_ =~ /^(\S[^\t]+\S)\t(\S[^\t]*\S)\t(\S+)/) ? [$1,$2,$3] : () }
168 :     `cat $newGdir/Subsystems/bindings`];
169 :     }
170 :     else
171 :     {
172 :     return $fig->get_genome_subsystem_data($genome);
173 :     }
174 :     }
175 :    
176 :     sub orgname_of_orgid {
177 :     my($self,$genome) = @_;
178 :    
179 :     return $self->genus_species($genome);
180 :     }
181 :    
182 :     sub orgid_of_orgname {
183 :     my($self,$genome_name) = @_;
184 :    
185 :     my @genomes = $self->genomes('complete');
186 :     my $i;
187 :     for ($i=0; ($i < @genomes) && ($genome_name ne $self->genus_species($genomes[$i])); $i++) {}
188 :     return ($i == @genomes) ? undef : $genomes[$i];
189 :     }
190 :    
191 :     sub genus_species_domain {
192 :     my($self,$genome) = @_;
193 :    
194 :     return [$self->genus_species($genome),$self->genome_domain($genome)];
195 :     }
196 :    
197 :     sub protein_subsystem_to_roles {
198 :     my ($self,$peg,$subsystem) = @_;
199 :    
200 :     my $fig = $self->{_fig};
201 :     my $newG = $self->{_genome};
202 :     my $newGdir = $self->{_orgdir};
203 :    
204 :     if (&FIG::genome_of($peg) ne $newG)
205 :     {
206 :     return $fig->protein_subsystem_to_roles($peg,$subsystem);
207 :     }
208 :     else
209 :     {
210 :     my @roles = map { (($_ =~ /^([^\t]+)\t([^\t]+)\t(\S+)$/) && ($1 eq $subsystem) && ($3 eq $peg)) ?
211 :     $2 : () } `cat $newGdir/Subsystems/bindings`;
212 :     my %roles = map { $_ => 1 } @roles;
213 :     return [sort keys(%roles)];
214 :     }
215 :     }
216 :    
217 :     sub contig_lengths {
218 :     my ($self, $genome) = @_;
219 :    
220 :     my $fig = $self->{_fig};
221 :     my $newG = $self->{_genome};
222 :     my $newGdir = $self->{_orgdir};
223 :    
224 :     if ($genome ne $newG)
225 :     {
226 :     return $fig->contig_lengths($genome);
227 :     }
228 :     else
229 :     {
230 :     my $contig_lengths = {};
231 :     if (open(CONTIGS,"<$newGdir/contigs"))
232 :     {
233 :     $/ = "\n>";
234 :     while (defined(my $x = <CONTIGS>))
235 :     {
236 :     if ($x =~ />?(\S+)[^\n]*\n(.*)/s)
237 :     {
238 :     my $id = $1;
239 :     my $seq = $2;
240 :     $seq =~ s/\s//gs;
241 :     $contig_lengths->{$id} = length($seq);
242 :     }
243 :     }
244 :     close(CONTIGS);
245 :     $/ = "\n";
246 :     }
247 :     return $contig_lengths;
248 :     }
249 :     }
250 :    
251 :     sub genome_szdna {
252 :     my ($self, $genome) = @_;
253 :    
254 :     my $fig = $self->{_fig};
255 :     my $newG = $self->{_genome};
256 :     my $newGdir = $self->{_orgdir};
257 :    
258 :     if ($genome ne $newG)
259 :     {
260 :     return $fig->genome_szdna($genome);
261 :     }
262 :     else
263 :     {
264 :     my $contig_lens = $self->contig_lengths($genome);
265 :     my $tot = 0;
266 :     while ( my($contig,$len) = each %$contig_lens)
267 :     {
268 :     $tot += $len;
269 :     }
270 :     return $tot;
271 :     }
272 :     }
273 :    
274 :     sub genome_version {
275 :     my ($self, $genome) = @_;
276 :    
277 :     my $fig = $self->{_fig};
278 :     my $newG = $self->{_genome};
279 :     my $newGdir = $self->{_orgdir};
280 :    
281 :     if ($genome ne $newG)
282 :     {
283 :     return $fig->genome_version($genome);
284 :     }
285 :     else
286 :     {
287 :     return "$genome.0";
288 :     }
289 :     }
290 :    
291 :     sub genome_pegs {
292 :     my ($self, $genome) = @_;
293 :    
294 :     my $fig = $self->{_fig};
295 :     my $newG = $self->{_genome};
296 :     my $newGdir = $self->{_orgdir};
297 :    
298 :     if ($genome ne $newG)
299 :     {
300 :     return $fig->genome_pegs($genome);
301 :     }
302 :     else
303 :     {
304 :     my @tmp = $self->all_features($genome,"peg");
305 :     my $n = @tmp;
306 :     return $n;
307 :     }
308 :     }
309 :    
310 :     sub genome_rnas {
311 :     my ($self, $genome) = @_;
312 :    
313 :     my $fig = $self->{_fig};
314 :     my $newG = $self->{_genome};
315 :     my $newGdir = $self->{_orgdir};
316 :    
317 :     if ($genome ne $newG)
318 :     {
319 :     return $fig->genome_rnas($genome);
320 :     }
321 :     else
322 :     {
323 :     my @tmp = $self->all_features($genome,"rna");
324 :     my $n = @tmp;
325 :     return $n;
326 :     }
327 :     }
328 :    
329 :     sub genome_domain {
330 :     my ($self, $genome) = @_;
331 :    
332 :     my $fig = $self->{_fig};
333 :     my $newG = $self->{_genome};
334 :     my $newGdir = $self->{_orgdir};
335 :    
336 :     if ($genome ne $newG)
337 :     {
338 :     return $fig->genome_domain($genome);
339 :     }
340 :     else
341 :     {
342 :     my $tax = $self->taxonomy_of($genome);
343 :     return ($tax =~ /^([^ \t;]+)/) ? $1 : "unknown";
344 :     }
345 :     }
346 : overbeek 1.2
347 :     sub genes_in_region {
348 :     my($self,$genome,$contig,$beg,$end) = @_;
349 :    
350 :     my $fig = $self->{_fig};
351 :     my $newG = $self->{_genome};
352 :     my $newGdir = $self->{_orgdir};
353 :    
354 :     if ($genome ne $newG)
355 :     {
356 :     return $fig->genes_in_region($genome,$contig,$beg,$end);
357 :     }
358 :     else
359 :     {
360 :     &load_tbl($self);
361 :     my $tblH = $self->{_tbl};
362 :     my $maxV = 0;
363 :     my $minV = 1000000000;
364 :     my $genes = [];
365 :     while ( my($fid,$tuple) = each %$tblH)
366 :     {
367 :     if (($tuple->[0]->[0] =~ /^(\S+)_(\d+)_\d+$/) && ($1 eq $contig))
368 :     {
369 :     my $beg1 = $2;
370 :     my $last = @{$tuple->[0]} - 1;
371 :     if (($tuple->[0]->[$last] =~ /^(\S+)_\d+_(\d+)$/) && ($1 eq $contig))
372 :     {
373 :     my $end1 = $2;
374 :     if (&overlaps($beg,$end,$beg1,$end1))
375 :     {
376 :     $minV = &FIG::min($minV,$beg1,$end1);
377 :     $maxV = &FIG::max($maxV,$beg1,$end1);
378 :     push(@$genes,$fid);
379 :     }
380 :     }
381 :     }
382 :     }
383 :     return ($genes,$minV,$maxV);
384 :     }
385 :     }
386 :    
387 :     sub overlaps {
388 :     my($b1,$e1,$b2,$e2) = @_;
389 :    
390 :     if ($b1 > $e1) { ($b1,$e1) = ($e1,$b1) }
391 :     if ($b2 > $e2) { ($b2,$e2) = ($e2,$b2) }
392 :     return &FIG::between($b1,$b2,$e1) || &FIG::between($b2,$b1,$e2);
393 :     }
394 :    
395 : overbeek 1.7 sub all_contigs {
396 :     my($self,$genome) = @_;
397 :    
398 :     my $fig = $self->{_fig};
399 :     my $newG = $self->{_genome};
400 :     my $newGdir = $self->{_orgdir};
401 :    
402 :     if ($genome ne $newG)
403 :     {
404 :     return $fig->all_contigs($genome);
405 :     }
406 :     else
407 :     {
408 :     &load_tbl($self);
409 :     my %contigs;
410 :     my $tblH = $self->{_tbl};
411 :     while ( my($fid,$tuple) = each %$tblH)
412 :     {
413 :     if ($tuple->[0]->[0] =~ /^(\S+)_\d+_\d+$/)
414 :     {
415 :     $contigs{$1} = 1;
416 :     }
417 :     }
418 :     return keys(%contigs);
419 :     }
420 :     }
421 :    
422 :     sub all_features {
423 :     my($self,$genome,$type) = @_;
424 :    
425 :     my $fig = $self->{_fig};
426 :     my $newG = $self->{_genome};
427 :     my $newGdir = $self->{_orgdir};
428 :    
429 :     if ($genome ne $newG)
430 :     {
431 :     return $fig->all_features($genome,$type);
432 :     }
433 :     else
434 :     {
435 :     &load_tbl($self);
436 :     my $tblH = $self->{_tbl};
437 :     return sort { &FIG::by_fig_id($a,$b) }
438 :     grep { ($_ =~ /^fig\|\d+\.\d+\.([^\.]+)/) && ($1 eq $type) }
439 :     keys(%$tblH);
440 :     }
441 :     }
442 :    
443 :     sub all_features_detailed_fast {
444 :     my($self,$genome) = @_;
445 :    
446 :     my $fig = $self->{_fig};
447 :     my $newG = $self->{_genome};
448 :     my $newGdir = $self->{_orgdir};
449 :    
450 :     if ($genome ne $newG)
451 :     {
452 :     return $fig->all_features_detailed_fast($genome);
453 :     }
454 :     else
455 :     {
456 :     &load_tbl($self);
457 :     my $tblH = $self->{_tbl};
458 :     my $feat_details = [];
459 :     while ( my($fid,$tuple) = each %$tblH)
460 :     {
461 :     if ($fid =~ /^fig\|\d+\.\d+\.(\S+)\.\d+/)
462 :     {
463 :     my $type = $1;
464 :     if ($tuple->[0]->[0] =~ /^\S+_(\d+)_(\d+)$/)
465 :     {
466 :     my($min,$max);
467 :     if ($1 < $2)
468 :     {
469 :     $min = $1;
470 :     $max = $2;
471 :     }
472 :     else
473 :     {
474 :     $min = $2;
475 :     $max = $1;
476 :     }
477 :     push(@$feat_details,[$fid,$tuple->[0]->[0],join(",",@{$tuple->[1]}),$type,$min,$max,$self->function_of($fid),'master','']);
478 :     }
479 :     }
480 :     }
481 :     return $feat_details;
482 :     }
483 :     }
484 :    
485 :     sub compute_clusters {
486 :     # Get the parameters.
487 :     my ($self, $pegList, $subsystem, $distance) = @_;
488 :     if (! defined $distance) {
489 :     $distance = 5000;
490 :     }
491 :    
492 :     my($peg,%by_contig);
493 :     foreach $peg (@$pegList)
494 :     {
495 :     my $loc;
496 :     if ($loc = $self->feature_location($peg))
497 :     {
498 :     my ($contig,$beg,$end) = &FIG::boundaries_of($loc);
499 :     my $genome = &FIG::genome_of($peg);
500 :     push(@{$by_contig{"$genome\t$contig"}},[($beg+$end)/2,$peg]);
501 :     }
502 :     }
503 :    
504 :     my @clusters = ();
505 :     foreach my $tuple (keys(%by_contig))
506 :     {
507 :     my $x = $by_contig{$tuple};
508 :     my @pegs = sort { $a->[0] <=> $b->[0] } @$x;
509 :     while ($x = shift @pegs)
510 :     {
511 :     my $clust = [$x->[1]];
512 :     while ((@pegs > 0) && (abs($pegs[0]->[0] - $x->[0]) <= $distance))
513 :     {
514 :     $x = shift @pegs;
515 :     push(@$clust,$x->[1]);
516 :     }
517 :    
518 :     if (@$clust > 1)
519 :     {
520 :     push(@clusters,$clust);
521 :     }
522 :     }
523 :     }
524 :     return sort { @$b <=> @$a } @clusters;
525 :     }
526 :    
527 :    
528 : overbeek 1.1 sub feature_location {
529 :     my($self,$fid) = @_;
530 :    
531 :     my $fig = $self->{_fig};
532 :     my $newG = $self->{_genome};
533 :     my $newGdir = $self->{_orgdir};
534 :    
535 :     if (($fid =~ /^fig\|(\d+\.\d+)/) && ($1 eq $newG))
536 :     {
537 :     &load_tbl($self);
538 :     if (my $x = $self->{_tbl}->{$fid})
539 :     {
540 :     return join(",",@{$x->[0]});
541 :     }
542 :     else
543 :     {
544 :     return undef;
545 :     }
546 :     }
547 :     else
548 :     {
549 :     return scalar $fig->feature_location($fid);
550 :     }
551 :     }
552 :    
553 :     sub function_of {
554 :     my($self,$fid) = @_;
555 :    
556 :     my $fig = $self->{_fig};
557 :     my $newG = $self->{_genome};
558 :     my $newGdir = $self->{_orgdir};
559 :    
560 :     if (($fid =~ /^fig\|(\d+\.\d+)/) && ($1 eq $newG))
561 :     {
562 :     &load_functions($self);
563 :     return $self->{_functions}->{$fid};
564 :     }
565 :     else
566 :     {
567 :     return scalar $fig->function_of($fid);
568 :     }
569 :     }
570 :    
571 :     sub feature_aliases {
572 :     my($self,$fid) = @_;
573 :    
574 :     my $fig = $self->{_fig};
575 :     my $newG = $self->{_genome};
576 :     my $newGdir = $self->{_orgdir};
577 :    
578 :     my @aliases;
579 :     if (($fid =~ /^fig\|(\d+\.\d+)/) && ($1 eq $newG))
580 :     {
581 :     &load_tbl($self);
582 :     if (my $x = $self->{_tbl}->{$fid})
583 :     {
584 :     @aliases = @{$x->[1]};
585 :     }
586 :     else
587 :     {
588 :     @aliases = ();
589 :     }
590 :     }
591 :     else
592 :     {
593 :     @aliases = $fig->feature_aliases($fid);
594 :     }
595 :     return wantarray() ? @aliases : join(",",@aliases);
596 :     }
597 :    
598 :     sub feature_annotations {
599 :     my($self,$fid,$rawtime) = @_;
600 :    
601 :     my $fig = $self->{_fig};
602 :     my $newG = $self->{_genome};
603 :     my $newGdir = $self->{_orgdir};
604 :    
605 :     my @annotations;
606 :     if (($fid =~ /^fig\|(\d+\.\d+)/) && ($1 eq $newG))
607 :     {
608 :     &load_ann($self);
609 :     if (my $x = $self->{_ann}->{$fid})
610 :     {
611 :     @annotations = @{$x};
612 :     }
613 :     else
614 :     {
615 :     @annotations = ();
616 :     }
617 :    
618 :     if ($rawtime)
619 :     {
620 :     return @annotations;
621 :     }
622 :     else
623 :     {
624 :     return map { $_->[1] = localtime($_->[1]); $_ } @annotations;
625 :     }
626 :     }
627 :     else
628 :     {
629 :     return $fig->feature_annotations($fid);
630 :     }
631 :     }
632 :    
633 :     sub get_translation {
634 :     my($self,$peg) = @_;
635 :    
636 :     my $fig = $self->{_fig};
637 :     my $newG = $self->{_genome};
638 :     my $newGdir = $self->{_orgdir};
639 :    
640 :     if (($peg =~ /^fig\|(\d+\.\d+)/) && ($1 eq $newG))
641 :     {
642 :     &load_pseq($self);
643 :     return $self->{_pseq}->{$peg};
644 :     }
645 :     else
646 :     {
647 :     return $fig->get_translation($peg);
648 :     }
649 :     }
650 :    
651 : olson 1.4 sub translation_length
652 :     {
653 :     my($self, $peg) = @_;
654 :    
655 :     my $fig = $self->{_fig};
656 :     my $newG = $self->{_genome};
657 :     my $newGdir = $self->{_orgdir};
658 :    
659 :     if (($peg =~ /^fig\|(\d+\.\d+)/) && ($1 eq $newG))
660 :     {
661 :     &load_pseq($self);
662 :     return length($self->{_pseq}->{$peg});
663 :     }
664 :     else
665 :     {
666 :     return $fig->translation_length($peg);
667 :     }
668 :     }
669 :    
670 :     sub translatable
671 :     {
672 :     my($self, $peg) = @_;
673 :    
674 :     my $fig = $self->{_fig};
675 :     my $newG = $self->{_genome};
676 :     my $newGdir = $self->{_orgdir};
677 :    
678 :     if (($peg =~ /^fig\|(\d+\.\d+)/) && ($1 eq $newG))
679 :     {
680 :     &load_pseq($self);
681 :     return length($self->{_pseq}->{$peg}) > 0 ? 1 : 0;;
682 :     }
683 :     else
684 :     {
685 :     return $fig->translatable($peg);
686 :     }
687 :     }
688 :    
689 :     sub is_real_feature
690 :     {
691 :     my($self, $fid) = @_;
692 :    
693 :     if ($self->is_virtual_feature($fid))
694 :     {
695 :     return 1;
696 :     }
697 :     else
698 :     {
699 :     return $self->{_fig}->is_real_feature($fid);
700 :     }
701 :    
702 :     }
703 :    
704 :     sub is_virtual_feature
705 :     {
706 :     my($self, $peg) = @_;
707 :    
708 :     my $fig = $self->{_fig};
709 :     my $newG = $self->{_genome};
710 :     my $newGdir = $self->{_orgdir};
711 :    
712 :     if (($peg =~ /^fig\|(\d+\.\d+)/) && ($1 eq $newG))
713 :     {
714 :     &load_pseq($self);
715 :     return 1;
716 :     }
717 :     else
718 :     {
719 :     return 0;
720 :     }
721 :     }
722 :    
723 : overbeek 1.1 sub load_pseq {
724 :     my($self) = @_;
725 :    
726 :     if ($self->{_pseq}) { return };
727 :    
728 :     my $newGdir = $self->{_orgdir};
729 :     my $pseq = {};
730 :     if (open(FASTA,"<$newGdir/Features/peg/fasta"))
731 :     {
732 :     $/ = "\n>";
733 :     my $x;
734 :     while (defined($x = <FASTA>))
735 :     {
736 :     if ($x =~ /^>?(\S+)[^\n]*\n(.*)/)
737 :     {
738 :     my $peg = $1;
739 :     my $seq = $2;
740 :     $seq =~ s/\s//gs;
741 :     $pseq->{$peg} = $seq;
742 :     }
743 :     }
744 :     close(FASTA);
745 :     $/ = "\n";
746 :     }
747 :     $self->{_pseq} = $pseq;
748 :     }
749 :    
750 :     sub load_tbl {
751 :     my($self) = @_;
752 :    
753 :     if ($self->{_tbl}) { return };
754 :    
755 :     my $newGdir = $self->{_orgdir};
756 :     my $tbl = {};
757 :     foreach my $x (`cat $newGdir/Features/*/tbl`)
758 :     {
759 :     if ($x =~ /^(\S+)\t(\S+)(\t(\S.*\S))?/)
760 :     {
761 :     my $fid = $1;
762 :     my $loc = [split(/,/,$2)];
763 :     my $aliases = $4 ? [split(/\t/,@$4)] : [];
764 :     $tbl->{$fid} = [$loc,$aliases];
765 :     }
766 :     }
767 :     $self->{_tbl} = $tbl;
768 :     }
769 :    
770 :     sub load_functions {
771 :     my($self) = @_;
772 :    
773 :     if ($self->{_functions}) { return };
774 :    
775 :     my $newG = $self->{_genome};
776 :     my $newGdir = $self->{_orgdir};
777 :    
778 :     my $functions = {};
779 :     foreach my $x (`cat $newGdir/*functions`)
780 :     {
781 :     if (($x =~ /^(fig\|(\d+\.\d+)\.\S+)\t(\S[^\t]*\S)/) && ($2 eq $newG))
782 :     {
783 :     $functions->{$1} = $3;
784 :     }
785 :     }
786 :     $self->{_functions} = $functions;
787 :     }
788 :    
789 :     sub sims {
790 :     my($self,$peg,$max,$maxP,$select) = @_;
791 :    
792 :     my $fig = $self->{_fig};
793 :     my $newG = $self->{_genome};
794 :     my $newGdir = $self->{_orgdir};
795 :     my $max = $max ? $max : 10000;
796 :     my $maxP = $maxP ? $maxP : 1.0e-5;
797 :     $select = $select ? $select : "all";
798 :    
799 :     if (($peg =~ /^fig\|(\d+\.\d+)/) && ($1 eq $newG))
800 :     {
801 :     &load_sims($self);
802 :     my @sims = ();
803 :     my $raw_sims = $self->{_sims}->{$peg};
804 :     if ($raw_sims)
805 :     {
806 :     foreach my $sim ( grep { $_->[10] <= $maxP } @$raw_sims )
807 :     {
808 :     my $id2 = $sim->id2;
809 :     my $id1 = $sim->id1;
810 :     my @relevant = ();
811 :    
812 :     my @maps_to = $fig->mapped_prot_ids( $id2 );
813 :     my $ref_len = $maps_to[0]->[1];
814 :    
815 :     @maps_to = grep { $_->[0] !~ /^xxx\d+/ } @maps_to;
816 :    
817 :     if ( $select =~ /^figx?$/ ) # Only fig
818 :     {
819 :     @relevant = grep { $_->[0] =~ /^fig/ } @maps_to;
820 :     }
821 :     else # All
822 :     {
823 :     @relevant = @maps_to;
824 :     }
825 :    
826 :     my $seen = {};
827 :     foreach my $x ( @relevant )
828 :     {
829 :     my ( $x_id, $x_ln ) = @$x;
830 :    
831 :     next if $seen->{$x_id};
832 :     $seen->{$x_id} = 1;
833 :    
834 :     my $delta2 = $ref_len - $x_ln; # Coordinate shift
835 :     my $sim1 = [ @$sim ]; # Make a copy
836 :     $sim1->[1] = $x_id;
837 :     $sim1->[8] -= $delta2;
838 :     $sim1->[9] -= $delta2;
839 :     bless( $sim1, "Sim" );
840 :     push( @sims, $sim1 );
841 :     }
842 :     }
843 :     }
844 :    
845 :     if (@sims > $max) { $#sims = $max-1; }
846 :     return @sims;
847 :     }
848 :     else
849 :     {
850 :     return $fig->sims($peg,$max,$maxP,$select);
851 :     }
852 :     }
853 :    
854 :     sub load_sims {
855 :     my($self) = @_;
856 :    
857 :     if ($self->{_sims}) { return };
858 :    
859 :     my $newGdir = $self->{_orgdir};
860 :    
861 :     my $sims = {};
862 :     foreach my $x (`cat $newGdir/similarities`)
863 :     {
864 :     chop; $x;
865 :     if ($x =~ /^(\S+)/)
866 :     {
867 :     push(@{$sims->{$1}},bless([split(/\t/,$x)],'Sim'));
868 :     }
869 :     }
870 :     $self->{_sims} = $sims;
871 :     }
872 :    
873 :     sub get_attributes {
874 :     my($self,$peg) = @_;
875 :    
876 :     my $fig = $self->{_fig};
877 :     my $newG = $self->{_genome};
878 :     my $newGdir = $self->{_orgdir};
879 :    
880 :     if (($peg =~ /^fig\|(\d+\.\d+)\.peg\.\d+$/) && ($1 eq $newG))
881 :     {
882 :     &load_attr($self);
883 :     if (my $x = $self->{_attr}->{$peg})
884 :     {
885 :     return @$x;
886 :     }
887 :     else
888 :     {
889 :     return ();
890 :     }
891 :     }
892 :     else
893 :     {
894 :     return $fig->get_attributes($peg);
895 :     }
896 :     }
897 :    
898 :     sub load_attr {
899 :     my($self) = @_;
900 :    
901 :     if ($self->{_attr}) { return };
902 :    
903 :     my $newGdir = $self->{_orgdir};
904 :     my $attr = {};
905 : overbeek 1.3 foreach my $x (`cat $newGdir/evidence.codes`)
906 : overbeek 1.1 {
907 :     if ($x =~ /^(\S+)\t(\S+)/)
908 :     {
909 :     push(@{$attr->{$1}},[$1,"evidence_code",$2,""]);
910 :     }
911 :     }
912 :     $self->{_attr} = $attr;
913 :     }
914 :    
915 :     sub load_ann {
916 :     my($self) = @_;
917 :    
918 :     if ($self->{_ann}) { return };
919 :    
920 :     my $newGdir = $self->{_orgdir};
921 :     my $ann = {};
922 :     if (open(ANN,"<$newGdir/annotations"))
923 :     {
924 :     $/ = "\n//\n";
925 :     while (defined(my $x = <ANN>))
926 :     {
927 :     chomp $x;
928 :     if ($x =~ /^(\S+)\n([^\n]+)\n([^\n]+)\n(.*)/s)
929 :     {
930 :     push(@{$ann->{$1}},[$1,$2,$3,"$4\n"]);
931 :     }
932 :     }
933 :     $/ = "\n";
934 :     close(ANN);
935 :     }
936 :     $self->{_ann} = $ann;
937 :     }
938 : overbeek 1.6
939 :     sub taxonomy_of {
940 :     my($self,$genome) = @_;
941 :     my $fig = $self->{_fig};
942 :     my $newG = $self->{_genome};
943 :     my $newGdir = $self->{_orgdir};
944 :    
945 :     if ($newG ne $genome)
946 :     {
947 :     return $fig->taxonomy_of($genome);
948 :     }
949 :    
950 :     my $tax;
951 :     if (open(TAX,"<$newGdir/TAXONOMY") && ($tax = <TAX>))
952 :     {
953 :     chop $tax;
954 :     return $tax;
955 :     }
956 :     else
957 :     {
958 :     return "unknown";
959 :     }
960 :     }
961 :    
962 :     sub build_tree_of_complete {
963 :     my($self,$min_for_label) = @_;
964 :     return $self->build_tree_of_all($min_for_label, "complete");
965 :     }
966 :    
967 :     sub build_tree_of_all {
968 :     my($self, $min_for_label, $complete)=@_;
969 :     my(@last,@tax,$i,$prefix,$lev,$genome,$tax);
970 :    
971 :     $min_for_label = $min_for_label ? $min_for_label : 10;
972 :     open(TMP,">/tmp/tree$$") || die "could not open /tmp/tree$$";
973 :     print TMP "1. root\n";
974 :    
975 :     @last = ();
976 :    
977 :    
978 :     foreach $genome (grep { ! $self->is_environmental($_) } $self->sort_genomes_by_taxonomy($self->genomes($complete)))
979 :     {
980 :     $tax = $self->taxonomy_of($genome);
981 :     @tax = split(/\s*;\s*/,$tax);
982 :     push(@tax,$genome);
983 :     for ($i=0; ((@last > $i) && (@tax > $i) && ($last[$i] eq $tax[$i])); $i++) {}
984 :     while ($i < @tax)
985 :     {
986 :     $lev = $i+2;
987 :     $prefix = " " x (4 * ($lev-1));
988 :     print TMP "$prefix$lev\. $tax[$i]\n";
989 :     $i++;
990 :     }
991 :     @last = @tax;
992 :     }
993 :     close(TMP);
994 :     my $tree = &tree_utilities::build_tree_from_outline("/tmp/tree$$");
995 :     $tree->[0] = 'All';
996 :     &FIG::limit_labels($tree,$min_for_label);
997 :     unlink("/tmp/tree$$");
998 :     return ($tree,&tree_utilities::tips_of_tree($tree));
999 :     }
1000 :    
1001 :     sub sort_genomes_by_taxonomy {
1002 :     my($self,@genomes) = @_;
1003 :    
1004 :     return map { $_->[0] }
1005 :     sort { $a->[1] cmp $b->[1] }
1006 :     map { [$_,$self->taxonomy_of($_)] }
1007 :     @genomes;
1008 :     }
1009 : overbeek 1.1
1010 : overbeek 1.6 sub taxonomic_groups_of_complete {
1011 :     my($self,$min_for_labels) = @_;
1012 :    
1013 :     my($tree,undef) = $self->build_tree_of_complete($min_for_labels);
1014 :     return &FIG::taxonomic_groups($tree);
1015 :     }
1016 :    
1017 : overbeek 1.1 1;

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3