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

Annotation of /FigKernelPackages/FIGV.pm

Parent Directory Parent Directory | Revision Log Revision Log


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

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3