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

Annotation of /FigKernelPackages/SeedV.pm

Parent Directory Parent Directory | Revision Log Revision Log


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

1 : overbeek 1.1 # -*- perl -*-
2 : olson 1.2
3 :     #
4 :     # This is a SAS Component
5 :     #
6 :    
7 : overbeek 1.1 #########################################################################
8 :     # Copyright (c) 2003-2008 University of Chicago and Fellowship
9 :     # for Interpretations of Genomes. All Rights Reserved.
10 :     #
11 :     # This file is part of the SEED Toolkit.
12 :     #
13 :     # The SEED Toolkit is free software. You can redistribute
14 :     # it and/or modify it under the terms of the SEED Toolkit
15 :     # Public License.
16 :     #
17 :     # You should have received a copy of the SEED Toolkit Public License
18 :     # along with this program; if not write to the University of Chicago
19 :     # at info@ci.uchicago.edu or the Fellowship for Interpretation of
20 :     # Genomes at veronika@thefig.info or download a copy from
21 :     # http://www.theseed.org/LICENSE.TXT.
22 :     #########################################################################
23 :    
24 :     package SeedV;
25 :    
26 :     use Carp;
27 :     use strict;
28 :     use Data::Dumper;
29 :     use DB_File;
30 :     use File::Basename;
31 :    
32 :     sub new {
33 :     my ($class, $org_dir ) = @_;
34 :    
35 :     $org_dir || confess("$org_dir must be a valid SEED directory");
36 :     my $self = {};
37 :     $self->{_orgdir} = $org_dir;
38 :     ($org_dir =~ /(\d+\.\d+)$/) || confess("$org_dir must be a path ending in the genome ID");
39 :     $self->{_genome} = $1;
40 :     return bless $self, $class;
41 :     }
42 :    
43 :     sub genome_id
44 :     {
45 :     return $_[0]->{_genome};
46 :     }
47 :    
48 :     # return the path of the organism directory
49 :     sub organism_directory {
50 :     return $_[0]->{_orgdir};
51 :     }
52 :    
53 :     sub get_basic_statistics
54 :     {
55 :     my($self) = @_;
56 :    
57 :     my $newG = $self->{_genome};
58 :     my $newGdir = $self->{_orgdir};
59 :    
60 :     #
61 :     # Check cache.
62 :     #
63 :    
64 :     my $cache = "$newGdir/cache.basic_statistics";
65 :     my $fh = new FileHandle($cache);
66 :     if ($fh)
67 :     {
68 :     my $stats = {};
69 :     while (<$fh>)
70 :     {
71 :     chomp;
72 :     my($k, $v) = split(/\t/);
73 :     $stats->{$k} = $v;
74 :     }
75 :     close($fh);
76 :     return $stats;
77 :     }
78 :    
79 :    
80 :     my $subsystem_data = $self->get_genome_subsystem_data();
81 :    
82 :     my %sscount = map { $_->[0] => 1 } @$subsystem_data;
83 :     my $nss=scalar(keys(%sscount));
84 :    
85 :     my $statistics = {
86 :     num_subsystems => $nss,
87 :     num_contigs => scalar($self->all_contigs()),
88 :     num_basepairs => $self->genome_szdna(),
89 :     genome_name => $self->genus_species(),
90 :     genome_domain => $self->genome_domain(),
91 :     genome_pegs => $self->genome_pegs(),
92 :     genome_rnas => $self->genome_rnas(),
93 :     };
94 :    
95 :     $fh = new FileHandle(">$cache");
96 :     if ($fh)
97 :     {
98 :     while (my($k, $v) = each %$statistics)
99 :     {
100 :     print $fh join("\t", $k, $v), "\n";
101 :     }
102 :     close($fh);
103 :     }
104 :     return $statistics;
105 :     }
106 :    
107 :    
108 :     sub get_peg_statistics {
109 :     my ($self) = @_;
110 :    
111 :     my $newG = $self->{_genome};
112 :     my $newGdir = $self->{_orgdir};
113 :    
114 :     my $cache = "$newGdir/cache.peg_statistics";
115 :     my $fh = new FileHandle($cache);
116 :     if ($fh)
117 :     {
118 :     my $stats = {};
119 :     while (<$fh>)
120 :     {
121 :     chomp;
122 :     my($k, $v) = split(/\t/);
123 :     $stats->{$k} = $v;
124 :     }
125 :     close($fh);
126 :     return $stats;
127 :     }
128 :    
129 :    
130 :     my $subsystem_data = $self->get_genome_subsystem_data();
131 :     my $assignment_data = $self->get_genome_assignment_data();
132 :    
133 :     my $hypo_sub = 0;
134 :     my $hypo_nosub = 0;
135 :     my $nothypo_sub = 0;
136 :     my $nothypo_nosub = 0;
137 :     my %in = map { $_->[2] => 1 } @$subsystem_data;
138 :     my $in = keys(%in);
139 :    
140 :     my %sscount = map { $_->[0] => 1 } @$subsystem_data;
141 :    
142 :     foreach $_ (@$assignment_data)
143 :     {
144 :     my($peg,$func) = @$_;
145 :    
146 :     my $is_hypo = &SeedUtils::hypo($func);
147 :    
148 :     if ($is_hypo && $in{$peg}) { $hypo_sub++ }
149 :     elsif ($is_hypo && ! $in{$peg}) { $hypo_nosub++ }
150 :     elsif ((! $is_hypo) && (! $in{$peg})) { $nothypo_nosub++ }
151 :     elsif ((! $is_hypo) && $in{$peg}) { $nothypo_sub++ }
152 :     }
153 :     my $tot = $hypo_sub + $nothypo_sub + $hypo_nosub + $nothypo_nosub;
154 :    
155 :     my ($fracHS, $fracNHS, $fracHNS, $fracNHNS);
156 :    
157 :     if ($tot == 0) {
158 :     $fracHS = sprintf "%.2f", 0.0;
159 :     $fracNHS = sprintf "%.2f", 0.0;
160 :     $fracHNS = sprintf "%.2f", 0.0;
161 :     $fracNHNS = sprintf "%.2f", 0.0;
162 :     } else {
163 :     $fracHS = sprintf "%.2f", $hypo_sub / $tot * 100;
164 :     $fracNHS = sprintf "%.2f", $nothypo_sub / $tot * 100;
165 :     $fracHNS = sprintf "%.2f", $hypo_nosub / $tot * 100;
166 :     $fracNHNS = sprintf "%.2f", $nothypo_nosub / $tot * 100;
167 :     }
168 :    
169 :     my $statistics = {
170 :     hypothetical_in_subsystem => $hypo_sub,
171 :     hypothetical_not_in_subsystem => $hypo_nosub,
172 :     non_hypothetical_in_subsystem => $nothypo_sub,
173 :     non_hypothetical_not_in_subsystem => $nothypo_nosub,
174 :     hypothetical_in_subsystem_percent => $fracHS,
175 :     hypothetical_not_in_subsystem_percent => $fracHNS,
176 :     non_hypothetical_in_subsystem_percent => $fracNHS,
177 :     non_hypothetical_not_in_subsystem_percent => $fracNHNS
178 :     };
179 :    
180 :     $fh = new FileHandle(">$cache");
181 :     if ($fh)
182 :     {
183 :     while (my($k, $v) = each %$statistics)
184 :     {
185 :     print $fh join("\t", $k, $v), "\n";
186 :     }
187 :     close($fh);
188 :     }
189 :    
190 :     return $statistics;
191 :     }
192 :     sub get_variant_and_bindings
193 :     {
194 :     my($self, $ssa) = @_;
195 :    
196 :     my $newG = $self->{_genome};
197 :     my $newGdir = $self->{_orgdir};
198 :    
199 :     $self->load_ss_data();
200 :    
201 :     my $bindings = $self->{_ss_bindings}->{$ssa};
202 :     my $variant = $self->{_ss_variants}->{$ssa};
203 :    
204 :     unless ($bindings) {
205 :     $variant = '*-1';
206 :     $bindings = {};
207 :     }
208 :    
209 :     return ($variant, $bindings);
210 :     }
211 :    
212 :     sub active_subsystems
213 :     {
214 :     my($self, $all) = @_;
215 :    
216 :     my $newG = $self->{_genome};
217 :     my $newGdir = $self->{_orgdir};
218 :    
219 :     $self->load_ss_data();
220 :    
221 :     my $slist = {};
222 :    
223 :     if ($self->{_ss_variants})
224 :     {
225 :     %{$slist} = %{$self->{_ss_variants}};
226 :     }
227 :    
228 :     if (not $all)
229 :     {
230 :     for my $ss (keys %$slist)
231 :     {
232 :     my $var = $slist->{$ss};
233 :     delete $slist->{$ss} if $var eq 0 or $var eq -1;
234 :     }
235 :     }
236 :     return $slist;
237 :     }
238 :    
239 :     sub peg_to_subsystems
240 :     {
241 :     my($self, $peg) = @_;
242 :    
243 :     my $newG = $self->{_genome};
244 :     my $newGdir = $self->{_orgdir};
245 :    
246 :     $self->load_ss_data();
247 :    
248 :     my $variant = $self->{_ss_variants};
249 :     my %sub = map { $_ => 1 }
250 :     grep { $variant->{$_} !~ /^(-1)|0$/ }
251 :     map { $_->[0] }
252 :     $self->peg_to_roles_in_subsystems($peg);
253 :     return sort keys(%sub);
254 :     }
255 :    
256 :     sub peg_to_roles_in_subsystems
257 :     {
258 :     my($self,$peg) = @_;
259 :    
260 :     my $newG = $self->{_genome};
261 :     my $newGdir = $self->{_orgdir};
262 :    
263 :     $self->load_ss_data();
264 :    
265 :     my $ret = $self->{_ss_peg_index}->{$peg};
266 :    
267 :     return $ret ? @$ret : ();
268 :     }
269 :    
270 :     sub subsystems_for_peg
271 :     {
272 :     my($self,$peg) = @_;
273 :     return $self->peg_to_roles_in_subsystems($peg);
274 :     }
275 :    
276 :     sub subsystems_for_peg_complete {
277 :     my ($self, $peg) = @_;
278 :    
279 :     my $newG = $self->{_genome};
280 :    
281 :     $self->load_ss_data();
282 :    
283 :     my $ret = $self->{_ss_peg_index}->{$peg};
284 :     if ($ret) {
285 :     return map { [ $_->[0], $_->[1], $self->{_ss_variants}->{$_->[0] }, $peg ] } @$ret;
286 :     } else {
287 :     return ();
288 :     }
289 :     }
290 :    
291 :     sub genus_species {
292 :     my($self) = @_;
293 :    
294 :     my $newG = $self->{_genome};
295 :     my $newGdir = $self->{_orgdir};
296 :    
297 :     if (open(GENOME,"<$newGdir/GENOME"))
298 :     {
299 :     my $x = <GENOME>;
300 :     close(GENOME);
301 :     chop $x;
302 :     return $x;
303 :     }
304 :     else
305 :     {
306 :     return "";
307 :     }
308 :     }
309 :    
310 :     sub get_genome_assignment_data {
311 :     my($self) = @_;
312 :    
313 :     my $newG = $self->{_genome};
314 :     my $newGdir = $self->{_orgdir};
315 :    
316 :     my @fnfiles = <$newGdir/proposed*functions>;
317 :     if (@fnfiles == 0)
318 :     {
319 :     @fnfiles = <$newGdir/assigned*functions>;
320 :     }
321 :    
322 :     if (@fnfiles == 0)
323 :     {
324 :     return [];
325 :     }
326 :     my %assign;
327 :     foreach $_ (`cat @fnfiles`)
328 :     {
329 :     if ( $_ =~ /^(fig\|\d+\.\d+\.peg\.\d+)\t(\S.*\S)/)
330 :     {
331 :     my($fid,$func) = ($1,$2);
332 :     $assign{$fid} = $func;
333 :     }
334 :     return [map { [$_,$assign{$_}] } sort { &SeedUtils::by_fig_id($a,$b) } keys(%assign)];
335 :     }
336 :     }
337 :    
338 :     sub get_genome_subsystem_data {
339 :     my($self,$genome) = @_;
340 :    
341 :     my $newG = $self->{_genome};
342 :     my $newGdir = $self->{_orgdir};
343 :    
344 :     my %operational = map { (($_ =~ /^(\S.*\S)\t(\S+)/) && (($2 ne '-1') && ($2 ne '0'))) ? ($1 => 1) : () }
345 :     `cat $newGdir/Subsystems/subsystems`;
346 :    
347 :     return [grep { ! $self->is_deleted_fid($_->[2]) }
348 :     map { (($_ =~ /^(\S[^\t]+\S)\t(\S[^\t]*\S)\t(\S+)/) && $operational{$1} ) ? [$1,$2,$3] : () }
349 :     `cat $newGdir/Subsystems/bindings`];
350 :     }
351 :    
352 :     sub get_genome_subsystem_count
353 :     {
354 :     my($self,$genome) = @_;
355 :    
356 :     my $newG = $self->{_genome};
357 :     my $newGdir = $self->{_orgdir};
358 :    
359 :    
360 :     my $count = 0;
361 :     my ($entry, $vc);
362 :     open(SUBSYSTEMS, "<$newGdir/Subsystems/subsystems");
363 :     while (defined($entry = <SUBSYSTEMS>))
364 :     {
365 :     chomp $entry;
366 :     (undef, $vc) = split /\t/, $entry;
367 :     if ($vc != -1) { ++$count; }
368 :     }
369 :     close(SUBSYSTEMS);
370 :     return $count;
371 :     }
372 :    
373 :     sub orgname_of_orgid {
374 :     my($self,$genome) = @_;
375 :    
376 :     return $self->genus_species();
377 :     }
378 :    
379 :     sub genus_species_domain {
380 :     my($self,$genome) = @_;
381 :    
382 :     return [$self->genus_species(),$self->genome_domain()];
383 :     }
384 :    
385 :     sub protein_subsystem_to_roles {
386 :     my ($self,$peg,$subsystem) = @_;
387 :    
388 :     my $newG = $self->{_genome};
389 :     my $newGdir = $self->{_orgdir};
390 :    
391 :     my @roles = map { (($_ =~ /^([^\t]+)\t([^\t]+)\t(\S+)$/) && ($1 eq $subsystem) && ($3 eq $peg)) ?
392 :     $2 : () } `cat $newGdir/Subsystems/bindings`;
393 :     my %roles = map { $_ => 1 } @roles;
394 :     return [sort keys(%roles)];
395 :     }
396 :    
397 :     sub contig_lengths {
398 :     my ($self) = @_;
399 :    
400 :     my $newG = $self->{_genome};
401 :     my $newGdir = $self->{_orgdir};
402 :    
403 :     my $contig_lengths = $self->{_contig_len_index};
404 :    
405 :     $self->load_contigs_index();
406 :     if (!tied(%$contig_lengths))
407 :     {
408 :     $self->load_contig_len_cache();
409 :     return $self->{_contig_len_cache};
410 :     }
411 :     return $contig_lengths;
412 :     }
413 :    
414 :     sub load_contig_len_cache
415 :     {
416 :     my($self) = @_;
417 :    
418 :     return if ref $self->{_contig_len_cache};
419 :    
420 :     my $newGdir = $self->{_orgdir};
421 :    
422 :     my $contig_lengths = {};
423 :     if (open(CONTIGS,"<$newGdir/contigs"))
424 :     {
425 :     local $/ = "\n>";
426 :     while (defined(my $x = <CONTIGS>))
427 :     {
428 :     chomp $x;
429 :     if ($x =~ />?(\S+)[^\n]*\n(.*)/s)
430 :     {
431 :     my $id = $1;
432 :     my $seq = $2;
433 :     $seq =~ s/\s//gs;
434 :     $contig_lengths->{$id} = length($seq);
435 :     }
436 :     }
437 :     close(CONTIGS);
438 :     }
439 :     $self->{_contig_len_cache} = $contig_lengths;
440 :     }
441 :    
442 :     sub contig_ln {
443 :     my ($self, $contig) = @_;
444 :    
445 :     my $newG = $self->{_genome};
446 :     my $newGdir = $self->{_orgdir};
447 :    
448 :     my $contig_len = $self->{_contig_len_index};
449 :     $self->load_contigs_index();
450 :     if (tied(%$contig_len))
451 :     {
452 :     return $contig_len->{$contig};
453 :     }
454 :     $self->load_contig_len_cache();
455 :     return $self->{_contig_len_cache}->{$contig};
456 :     }
457 :    
458 :     sub contigs_of
459 :     {
460 :     my ($self) = @_;
461 :    
462 :     my $newG = $self->{_genome};
463 :     my $newGdir = $self->{_orgdir};
464 :    
465 :     my @out;
466 :     $self->load_contigs_index();
467 :    
468 :     my $contigs = $self->{_contigs_index};
469 :     if (tied(%$contigs))
470 :     {
471 :     return keys %$contigs;
472 :     }
473 :    
474 :     $self->load_contig_len_cache();
475 :    
476 :     return keys %{$self->{_contig_len_cache}};
477 :     }
478 :    
479 :     =head3 dna_seq
480 :    
481 :     usage: $seq = dna_seq(@locations)
482 :    
483 :     Returns the concatenated subsequences described by the list of locations. Each location
484 :     must be of the form
485 :    
486 :     Contig_Beg_End
487 :    
488 :     where Contig must be the ID of a contig for genome $genome. If Beg > End the location
489 :     describes a stretch of the complementary strand.
490 :    
491 :     =cut
492 :     #: Return Type $;
493 :     sub dna_seq {
494 :     my($self,@locations) = @_;
495 :    
496 :     my $newG = $self->{_genome};
497 :     my $newGdir = $self->{_orgdir};
498 :    
499 :     my $contigs = $self->{_contigs_index};
500 :     if (!tied %$contigs)
501 :     {
502 :     $self->load_contig_seq();
503 :     $contigs = $self->{_contigs_seq_cache};
504 :     }
505 :    
506 :     my(@pieces,$loc,$contig,$beg,$end,$ln,$rdbH);
507 :    
508 :     @locations = map { split(/,/,$_) } @locations;
509 :     @pieces = ();
510 :     foreach $loc (@locations)
511 :     {
512 :     if ($loc =~ /^(\S+)_(\d+)_(\d+)$/)
513 :     {
514 :     ($contig,$beg,$end) = ($1,$2,$3);
515 :     my $seq = $contigs->{$contig};
516 :    
517 :     $ln = length($seq);
518 :    
519 :     if (! $ln) {
520 :     print STDERR "$contig: could not get length\n";
521 :     return "";
522 :     }
523 :    
524 :     if (&SeedUtils::between(1,$beg,$ln) && &SeedUtils::between(1,$end,$ln))
525 :     {
526 :     if ($beg < $end)
527 :     {
528 :     push(@pieces, substr($seq, $beg - 1, ($end - $beg) + 1));
529 :     }
530 :     else
531 :     {
532 :     push(@pieces, &SeedUtils::reverse_comp(substr($seq, $end - 1, ($beg - $end) + 1)));
533 :     }
534 :     }
535 :     }
536 :     }
537 :     return lc(join("",@pieces));
538 :     }
539 :    
540 :     sub load_contig_seq
541 :     {
542 :     my($self) = @_;
543 :    
544 :     return if ref($self->{_contigs_seq_cache});
545 :    
546 :     my $newGdir = $self->{_orgdir};
547 :    
548 :     my $contigs = {};
549 :    
550 :     if (open(CONTIGS,"<$newGdir/contigs"))
551 :     {
552 :     local $/ = "\n>";
553 :     while (defined(my $x = <CONTIGS>))
554 :     {
555 :     chomp $x;
556 :     if ($x =~ />?(\S+)[^\n]*\n(.*)/s)
557 :     {
558 :     my $id = $1;
559 :     my $seq = $2;
560 :     $seq =~ s/\s//gs;
561 :     $contigs->{$id} = $seq;
562 :     }
563 :     }
564 :     close(CONTIGS);
565 :     }
566 :     $self->{_contigs_seq_cache} = $contigs;
567 :     }
568 :    
569 :     sub genome_szdna {
570 :     my ($self) = @_;
571 :    
572 :     my $newG = $self->{_genome};
573 :     my $newGdir = $self->{_orgdir};
574 :    
575 :     my $contig_lens = $self->contig_lengths();
576 :     my $tot = 0;
577 :     while ( my($contig,$len) = each %$contig_lens)
578 :     {
579 :     $tot += $len;
580 :     }
581 :     return $tot;
582 :     }
583 :    
584 :     sub genome_pegs {
585 :     my ($self) = @_;
586 :    
587 :     my $newG = $self->{_genome};
588 :     my $newGdir = $self->{_orgdir};
589 :    
590 :     my @tmp = $self->all_features("peg");
591 :     my $n = @tmp;
592 :     return $n;
593 :     }
594 :    
595 :     sub genome_rnas {
596 :     my ($self) = @_;
597 :    
598 :     my $newG = $self->{_genome};
599 :     my $newGdir = $self->{_orgdir};
600 :    
601 :     my @tmp = $self->all_features("rna");
602 :     my $n = @tmp;
603 :     return $n;
604 :     }
605 :    
606 :     sub genome_domain {
607 :     my ($self) = @_;
608 :    
609 :     my $newG = $self->{_genome};
610 :     my $newGdir = $self->{_orgdir};
611 :    
612 :     my $tax = $self->taxonomy_of();
613 :     return ($tax =~ /^([^ \t;]+)/) ? $1 : "unknown";
614 :     }
615 :    
616 :     sub genes_in_region {
617 :     my($self,$contig,$beg,$end) = @_;
618 :    
619 :     my $newG = $self->{_genome};
620 :     my $newGdir = $self->{_orgdir};
621 :    
622 :     $self->load_feature_indexes();
623 :     #
624 :     # Use the recno index if exists.
625 :     #
626 :    
627 :     my $maxV = 0;
628 :     my $minV = 1000000000;
629 :     my $genes = [];
630 :    
631 :     my $recnos = $self->{_feat_recno};
632 :    
633 :     if (ref($recnos))
634 :     {
635 :     while (my($ftype, $list) = each (%$recnos))
636 :     {
637 :     #
638 :     # Look up start/end of this contig in the btree index.
639 :     #
640 :    
641 :     my $inf = $self->{_feat_btree}->{$ftype}->{$contig};
642 :     my($istart, $iend);
643 :     if ($inf)
644 :     {
645 :     ($istart, $iend) = split(/$;/, $inf);
646 :     }
647 :     else
648 :     {
649 :     $istart = 0;
650 :     $iend = $#$list;
651 :     }
652 :    
653 :     for (my $idx = $istart; $idx <= $iend; $idx++)
654 :     {
655 :     my($fid, $fcontig, $fbeg, $fend, $fstrand) = split(/$;/, $list->[$idx]);
656 :     if ($contig eq $fcontig and &overlaps($beg, $end, $fbeg, $fend))
657 :     {
658 :     $minV = &SeedUtils::min($minV,$fbeg,$fend);
659 :     $maxV = &SeedUtils::max($maxV,$fbeg,$fend);
660 :     push(@$genes,$fid);
661 :     }
662 :     }
663 :     }
664 :     }
665 :     else
666 :     {
667 :     &load_tbl($self);
668 :     my $tblH = $self->{_tbl};
669 :     while ( my($fid,$tuple) = each %$tblH)
670 :     {
671 :     if (($tuple->[0]->[0] =~ /^(\S+)_(\d+)_\d+$/) && ($1 eq $contig))
672 :     {
673 :     my $beg1 = $2;
674 :     my $last = @{$tuple->[0]} - 1;
675 :     if (($tuple->[0]->[$last] =~ /^(\S+)_\d+_(\d+)$/) && ($1 eq $contig))
676 :     {
677 :     my $end1 = $2;
678 :     if (&overlaps($beg,$end,$beg1,$end1))
679 :     {
680 :     $minV = &SeedUtils::min($minV,$beg1,$end1);
681 :     $maxV = &SeedUtils::max($maxV,$beg1,$end1);
682 :     push(@$genes,$fid);
683 :     }
684 :     }
685 :     }
686 :     }
687 :     }
688 :     return ($genes,$minV,$maxV);
689 :     }
690 :    
691 :     sub overlaps {
692 :     my($b1,$e1,$b2,$e2) = @_;
693 :    
694 :     if ($b1 > $e1) { ($b1,$e1) = ($e1,$b1) }
695 :     if ($b2 > $e2) { ($b2,$e2) = ($e2,$b2) }
696 :     return &SeedUtils::between($b1,$b2,$e1) || &SeedUtils::between($b2,$b1,$e2);
697 :     }
698 :    
699 :     sub all_contigs {
700 :     my($self) = @_;
701 :    
702 :     my $newG = $self->{_genome};
703 :     my $newGdir = $self->{_orgdir};
704 :    
705 :     $self->load_feature_indexes();
706 :    
707 :     my %contigs;
708 :     my $btrees = $self->{_feat_btree};
709 :    
710 :     if (ref($btrees))
711 :     {
712 :     while (my($ftype, $btree) = each (%$btrees))
713 :     {
714 :     map { $contigs{$_} = 1 } grep { !/^fig/ } keys %$btree ;
715 :     }
716 :     }
717 :     else
718 :     {
719 :     &load_tbl($self);
720 :     my $tblH = $self->{_tbl};
721 :     while ( my($fid,$tuple) = each %$tblH)
722 :     {
723 :     if ($tuple->[0]->[0] =~ /^(\S+)_\d+_\d+$/)
724 :     {
725 :     $contigs{$1} = 1;
726 :     }
727 :     }
728 :     }
729 :     return keys(%contigs);
730 :     }
731 :    
732 :     sub all_features {
733 :     my($self,$type) = @_;
734 :     if (not defined($type)) { $type = qq(); }
735 :    
736 :     my $newG = $self->{_genome};
737 :     my $newGdir = $self->{_orgdir};
738 :    
739 :     #warn "Loading feature indices";
740 :     $self->load_feature_indexes();
741 :    
742 :     my %contigs;
743 :     my $btrees = $self->{_feat_btree};
744 :    
745 :     if (ref($btrees)) {
746 :     warn "B-tree already loaded" if $ENV{FIG_VERBOSE};
747 :    
748 :     my $btree = $btrees->{$type};
749 :     return sort { &SeedUtils::by_fig_id($a, $b) }
750 :     grep { /^fig/ } keys %$btree;
751 :     }
752 :     else {
753 :     warn "Loading contig B-tree" if $ENV{FIG_VERBOSE};
754 :    
755 :     &load_tbl($self);
756 :     my $tblH = $self->{_tbl};
757 :    
758 :     return sort {
759 :     &SeedUtils::by_fig_id($a,$b)
760 :     } grep {
761 :     #...NOTE: Matches all feature types if $type is the null string
762 :     ($_ =~ /^fig\|\d+\.\d+\.([^\.]+)/) && ($1 =~ m/$type/)
763 :     } keys(%$tblH);
764 :     }
765 :     }
766 :    
767 :     sub all_features_detailed_fast {
768 :     my($self,$regmin, $regmax, $contig) = @_;
769 :    
770 :     my $newG = $self->{_genome};
771 :     my $newGdir = $self->{_orgdir};
772 :    
773 :     $self->load_feature_indexes();
774 :     my $feat_details = [];
775 :     my $recnos = $self->{_feat_recno};
776 :    
777 :     if (ref($recnos))
778 :     {
779 :     while (my($ftype, $list) = each (%$recnos))
780 :     {
781 :     #
782 :     # Look up start/end of this contig in the btree index.
783 :     #
784 :    
785 :     my $inf = $self->{_feat_btree}->{$ftype}->{$contig};
786 :     my($istart, $iend);
787 :     if ($inf)
788 :     {
789 :     ($istart, $iend) = split(/$;/, $inf);
790 :     }
791 :     else
792 :     {
793 :     $istart = 0;
794 :     $iend = $#$list;
795 :     }
796 :    
797 :     for (my $idx = $istart; $idx <= $iend; $idx++)
798 :     {
799 :     my($fid, $fcontig, $fbeg, $fend, $fstrand) = split(/$;/, $list->[$idx]);
800 :    
801 :     if (not defined($regmin) or not defined($regmax) or not defined($contig) or
802 :     (($contig eq $fcontig) and
803 :     ($fbeg < $regmin and $regmin < $fend) or ($fbeg < $regmax and $regmax < $fend) or ($fbeg > $regmin and $fend < $regmax)))
804 :     {
805 :     my($loc, $index, @aliases) = split(/$;/, $self->{_feat_btree}->{$ftype}->{$fid});
806 :    
807 :     my $function = $self->function_of($fid);
808 :     push(@$feat_details,[$fid, $loc, join(",", @aliases), $ftype, $fbeg, $fend, $function,'master','']);
809 :     }
810 :     }
811 :     }
812 :     }
813 :     else
814 :     {
815 :     &load_tbl($self);
816 :     my $tblH = $self->{_tbl};
817 :     while ( my($fid,$tuple) = each %$tblH)
818 :     {
819 :     if ($fid =~ /^fig\|\d+\.\d+\.(\S+)\.\d+/)
820 :     {
821 :     my $type = $1;
822 :     if ($tuple->[0]->[0] =~ /^(\S+)_(\d+)_(\d+)$/)
823 :     {
824 :     my($ctg, $beg, $end) = ($1, $2, $3);
825 :     next if (defined($contig) and $contig ne $ctg);
826 :    
827 :     my($min,$max);
828 :     if ($beg < $end)
829 :     {
830 :     $min = $beg;
831 :     $max = $end;
832 :     }
833 :     else
834 :     {
835 :     $min = $end;
836 :     $max = $beg;
837 :     }
838 :    
839 :     if (not defined($regmin) or not defined($regmax) or
840 :     ($min < $regmin and $regmin < $max) or ($min < $regmax and $regmax < $max) or ($min > $regmin and $max < $regmax))
841 :     {
842 :     my $function = $self->function_of($fid);
843 :     push(@$feat_details,[$fid,$tuple->[0]->[0],join(",",@{$tuple->[1]}),$type,$min,$max,$function,'master','']);
844 :     }
845 :     }
846 :     }
847 :     }
848 :     }
849 :     return $feat_details;
850 :     }
851 :    
852 :     sub feature_location {
853 :     my($self,$fid) = @_;
854 :    
855 :     my $newG = $self->{_genome};
856 :     my $newGdir = $self->{_orgdir};
857 :    
858 :     if (($fid =~ /^fig\|(\d+\.\d+)\.([^.]+)/) && ($1 eq $newG))
859 :     {
860 :     my $ftype = $2;
861 :    
862 :     $self->load_feature_indexes();
863 :    
864 :     my $btree = $self->{_feat_btree}->{$ftype};
865 :     if ($btree)
866 :     {
867 :     my($loc, $idx, @aliases) = split(/$;/, $btree->{$fid});
868 :     return wantarray ? split(/,/, $loc) : $loc;
869 :     }
870 :     else
871 :     {
872 :     &load_tbl($self);
873 :     if (my $x = $self->{_tbl}->{$fid})
874 :     {
875 :     if (wantarray)
876 :     {
877 :     return @{$x->[0]};
878 :     }
879 :     else
880 :     {
881 :     return join(",",@{$x->[0]});
882 :     }
883 :     }
884 :     else
885 :     {
886 :     return undef;
887 :     }
888 :     }
889 :     }
890 :     return undef;
891 :     }
892 :    
893 :     sub function_of {
894 :     my($self,$fid) = @_;
895 :    
896 :     my $newG = $self->{_genome};
897 :     my $newGdir = $self->{_orgdir};
898 :    
899 :     if (($fid =~ /^fig\|(\d+\.\d+)/) && ($1 eq $newG))
900 :     {
901 :     &load_functions($self);
902 :    
903 :     my $fn = $self->{_functions}->{$fid};
904 :     if (wantarray)
905 :     {
906 :     return ['master', $fn];
907 :     }
908 :     else
909 :     {
910 :     return $fn;
911 :     }
912 :     }
913 :     return "";
914 :     }
915 :    
916 :    
917 :     =pod
918 :    
919 :     find_features_by_annotation
920 :    
921 :     Takes a reference to a hash of functions to find and an optional case boolean, and returns a hash with keys are the function and values are a reference to an array of the IDs that have that function.
922 :    
923 :     If the case boolean is set the search is case insensitive. Otherwise case sensitive.
924 :    
925 :     =cut
926 :    
927 :     sub find_features_by_annotation {
928 :     my($self,$anno_hash, $case)=@_;
929 :     $self->load_functions;
930 :    
931 :     if ($case) {map {$anno_hash->{uc($_)}=1} keys %$anno_hash}
932 :    
933 :     my $res={};
934 :     foreach my $id (keys %{$self->{_functions}})
935 :     {
936 :     my $fn = $self->{_functions}->{$id};
937 :     $case ? $fn = uc($fn) : 1;
938 :     if ($anno_hash->{$fn}) {push @{$res->{$fn}}, $id}
939 :     }
940 :    
941 :     return $res;
942 :     }
943 :    
944 :    
945 :     =pod
946 :    
947 :     search_features_by_annotation
948 :    
949 :     Takes a string to find and an optional case boolean, and returns a hash with keys are the function and values are a reference to an array of the IDs that have that function.
950 :    
951 :     If the case boolean is set the search is case insensitive. Otherwise case sensitive.
952 :    
953 :     Note that this was originally based on the find above, but this uses a regexp for matching. Will likely be a lot slower which is why I only search for a single term. There may be an SQL way of doing this, if so let Rob know how to do it and I'll replace this method.
954 :    
955 :    
956 :     =cut
957 :    
958 :     sub search_features_by_annotation {
959 :     my($self,$term, $case)=@_;
960 :     $self->load_functions;
961 :    
962 :     # to make case insensitive convert everything to uppercase
963 :     # alternative is to use two regexps, one for case insens and one for not case insens
964 :     # but the bad thing about that approach is that if you have a case insensitive search
965 :     # you do two regexps for each failed match
966 :    
967 :     $case ? $term = uc($term) : 1;
968 :    
969 :     my $res={};
970 :     foreach my $id (keys %{$self->{_functions}})
971 :     {
972 :     # we set two variables, one that has case changed for case insensitive searches
973 :     my $fn = my $fnc = $self->{_functions}->{$id};
974 :     $case ? $fn = uc($fn) : 1;
975 :     if ($fn =~ m/$term/) {push @{$res->{$fnc}}, $id}
976 :     }
977 :    
978 :     return $res;
979 :     }
980 :    
981 :    
982 :     sub feature_aliases {
983 :     my($self,$fid) = @_;
984 :    
985 :     my $newG = $self->{_genome};
986 :     my $newGdir = $self->{_orgdir};
987 :    
988 :     my @aliases;
989 :    
990 :     if (($fid =~ /^fig\|(\d+\.\d+)/) && ($1 eq $newG))
991 :     {
992 :     &load_tbl($self);
993 :     if (my $x = $self->{_tbl}->{$fid})
994 :     {
995 :     @aliases = @{$x->[1]};
996 :     }
997 :     else
998 :     {
999 :     @aliases = ();
1000 :     }
1001 :     }
1002 :     else
1003 :     {
1004 :     @aliases = ();
1005 :     }
1006 :     return wantarray() ? @aliases : join(",",@aliases);
1007 :     }
1008 :    
1009 :     sub get_corresponding_ids {
1010 :     my($self, $id, $with_type_info) = @_;
1011 :    
1012 :     my $newG = $self->{_genome};
1013 :    
1014 :     if (($id =~ /^fig\|(\d+\.\d+)/) && ($1 eq $newG)) {
1015 :     my @aliases = $self->feature_aliases($id);
1016 :     my @corresponding_ids = ();
1017 :     foreach my $alias (@aliases) {
1018 :     if ($alias =~ /^gi\|/) {
1019 :     if ($with_type_info) {
1020 :     push(@corresponding_ids, [$alias, 'NCBI']);
1021 :     } else {
1022 :     push(@corresponding_ids, $alias);
1023 :     }
1024 :     last;
1025 :     }
1026 :     }
1027 :     return @corresponding_ids;
1028 :     }
1029 :     }
1030 :    
1031 :     sub feature_annotations {
1032 :     my($self,$fid,$rawtime) = @_;
1033 :    
1034 :     my $newG = $self->{_genome};
1035 :     my $newGdir = $self->{_orgdir};
1036 :    
1037 :     my @annotations;
1038 :     if (($fid =~ /^fig\|(\d+\.\d+)/) && ($1 eq $newG))
1039 :     {
1040 :     &load_ann($self);
1041 :     if (my $x = $self->{_ann}->{$fid})
1042 :     {
1043 :     @annotations = @{$x};
1044 :     }
1045 :     else
1046 :     {
1047 :     @annotations = ();
1048 :     }
1049 :    
1050 :     if ($rawtime)
1051 :     {
1052 :     return @annotations;
1053 :     }
1054 :     else
1055 :     {
1056 :     return map { my $r = [@$_]; $r->[1] = localtime($r->[1]); $r } @annotations;
1057 :     }
1058 :     }
1059 :     return ();
1060 :     }
1061 :    
1062 :     sub get_translation {
1063 :     my($self,$peg) = @_;
1064 :    
1065 :     my $newG = $self->{_genome};
1066 :     my $newGdir = $self->{_orgdir};
1067 :    
1068 :     if (($peg =~ /^fig\|(\d+\.\d+)/) && ($1 eq $newG))
1069 :     {
1070 :     &load_feature_indexes($self);
1071 :    
1072 :     my $out = $self->{_feat}->{peg}->{$peg};
1073 :     if (defined($out))
1074 :     {
1075 :     return $out;
1076 :     }
1077 :    
1078 :     #
1079 :     # If we have a blast-formatted fasta, use fastacmd to
1080 :     # do the lookup, and cache the output for later use.
1081 :     #
1082 :    
1083 :     if ($self->{_feat_fasta}->{peg})
1084 :     {
1085 :     my $id = "gnl|$peg";
1086 :     my $cmd = "$FIG_Config::ext_bin/fastacmd -d $self->{_feat_fasta}->{peg} -s '$id'";
1087 :     open(P, "$cmd|") or die "get_translation: cmd failed with $!: $cmd";
1088 :     $_ = <P>;
1089 :     my $out;
1090 :     while (<P>)
1091 :     {
1092 :     s/\s+//g;
1093 :     $out .= $_;
1094 :     }
1095 :     close(P);
1096 :     $self->{_feat}->{$peg} = $out;
1097 :     return $out;
1098 :     }
1099 :     else
1100 :     {
1101 :     return $self->{_feat}->{$peg};
1102 :     }
1103 :     }
1104 :     else
1105 :     {
1106 :     return undef;
1107 :     }
1108 :     }
1109 :    
1110 :     sub translation_length
1111 :     {
1112 :     my($self, $peg) = @_;
1113 :    
1114 :     my $newG = $self->{_genome};
1115 :     my $newGdir = $self->{_orgdir};
1116 :    
1117 :     if (($peg =~ /^fig\|(\d+\.\d+)/) && ($1 eq $newG))
1118 :     {
1119 :     my $t = $self->get_translation($peg);
1120 :     return length($t);
1121 :     }
1122 :     return undef;
1123 :     }
1124 :    
1125 :     sub load_feature_hash {
1126 :     my($self,$type) = @_;
1127 :    
1128 :     my $newGdir = $self->{_orgdir};
1129 :     my $typeH = {};
1130 :     if (open(FIDS,"<$newGdir/Features/$type/tbl"))
1131 :     {
1132 :     while ($_ = <FIDS>)
1133 :     {
1134 :     if ($_ =~ /^(\S+)/)
1135 :     {
1136 :     my $fid = $1;
1137 :     if (! $self->is_deleted_fid($fid))
1138 :     {
1139 :     $typeH->{$fid} = 1;
1140 :     }
1141 :     }
1142 :     }
1143 :     close(FIDS);
1144 :     }
1145 :     return $typeH;
1146 :     }
1147 :    
1148 :     sub load_feature_indexes
1149 :     {
1150 :     my($self) = @_;
1151 :    
1152 :     if ($self->{_feat}) { return };
1153 :    
1154 :     my $newGdir = $self->{_orgdir};
1155 :    
1156 :     for my $fdir (<$newGdir/Features/*>)
1157 :     {
1158 :     my $ftype = basename($fdir);
1159 :    
1160 :     #
1161 :     # If we have a tbl.btree, tie that for our use.
1162 :     #
1163 :    
1164 :     my $tbl_idx = {};
1165 :     my $tie = tie %$tbl_idx, 'DB_File', "$fdir/tbl.btree", O_RDONLY, 0666, $DB_BTREE;
1166 :     if ($tie)
1167 :     {
1168 :     $self->{_feat_tie}->{$ftype} = $tie;
1169 :     $self->{_feat_btree}->{$ftype} = $tbl_idx;
1170 :     }
1171 :    
1172 :     my $tbl_list = [];
1173 :     my $ltie = tie @$tbl_list, 'DB_File', "$fdir/tbl.recno", O_RDONLY, 0666, $DB_RECNO;
1174 :     if ($tie)
1175 :     {
1176 :     $self->{_feat_ltie}->{$ftype} = $ltie;
1177 :     $self->{_feat_recno}->{$ftype} = $tbl_list;
1178 :     }
1179 :    
1180 :     #
1181 :     # If we have fasta.norm.phr, set _pseq_fasta to the fasta file to use with fastacmd.
1182 :     #
1183 :    
1184 :     my $pseq = {};
1185 :    
1186 :     if (-f "$fdir/fasta.norm.phr")
1187 :     {
1188 :     $self->{_feat_fasta}->{$ftype} = "$fdir/fasta.norm";
1189 :    
1190 :     }
1191 :     else
1192 :     {
1193 :     #
1194 :     # Otherwise, we need to load the data.
1195 :     #
1196 :    
1197 :     if (open(FASTA,"<$newGdir/Features/peg/fasta"))
1198 :     {
1199 :     local $/ = "\n>";
1200 :     my $x;
1201 :     while (defined($x = <FASTA>))
1202 :     {
1203 :     chomp $x;
1204 :     if ($x =~ />?(\S+)[^\n]*\n(.*)/s)
1205 :     {
1206 :     my $peg = $1;
1207 :     my $seq = $2;
1208 :     $seq =~ s/\s//gs;
1209 :     if (! $self->is_deleted_fid($peg))
1210 :     {
1211 :     $pseq->{$peg} = $seq;
1212 :     }
1213 :     }
1214 :     }
1215 :     close(FASTA);
1216 :     }
1217 :     }
1218 :     $self->{_feat}->{$ftype} = $pseq;
1219 :     }
1220 :     }
1221 :    
1222 :     sub load_pseq {
1223 :     my($self) = @_;
1224 :    
1225 :     if ($self->{_pseq}) { return };
1226 :    
1227 :     my $newGdir = $self->{_orgdir};
1228 :     my $fdir = "$newGdir/Features/peg";
1229 :    
1230 :     #
1231 :     # If we have a tbl.btree, tie that for our use.
1232 :     #
1233 :    
1234 :     my $tbl_idx = {};
1235 :     my $tie = tie %$tbl_idx, 'DB_File', "$fdir/tbl.btree", O_RDONLY, 0666, $DB_BTREE;
1236 :     if ($tie)
1237 :     {
1238 :     $self->{_tbl_tie} = $tie;
1239 :     $self->{_tbl_btree} = $tbl_idx;
1240 :     }
1241 :    
1242 :     #
1243 :     # If we have fasta.norm.phr, set _pseq_fasta to the fasta file to use with fastacmd.
1244 :     #
1245 :    
1246 :     my $pseq = {};
1247 :    
1248 :     if (-f "$fdir/fasta.norm.phr")
1249 :     {
1250 :     $self->{_pseq_fasta} = "$fdir/fasta.norm";
1251 :     }
1252 :     else
1253 :     {
1254 :     #
1255 :     # Otherwise, we need to load the data.
1256 :     #
1257 :    
1258 :     if (open(FASTA,"<$newGdir/Features/peg/fasta"))
1259 :     {
1260 :     local $/ = "\n>";
1261 :     my $x;
1262 :     while (defined($x = <FASTA>))
1263 :     {
1264 :     chomp $x;
1265 :     if ($x =~ />?(\S+)[^\n]*\n(.*)/s)
1266 :     {
1267 :     my $peg = $1;
1268 :     my $seq = $2;
1269 :     $seq =~ s/\s//gs;
1270 :     if (! $self->is_deleted_fid($peg))
1271 :     {
1272 :     $pseq->{$peg} = $seq;
1273 :     }
1274 :     }
1275 :     }
1276 :     close(FASTA);
1277 :     }
1278 :     }
1279 :     $self->{_pseq} = $pseq;
1280 :     }
1281 :    
1282 :     sub load_ss_data
1283 :     {
1284 :     my($self, $force) = @_;
1285 :    
1286 :     return if defined($self->{_ss_bindings}) && (! $force);
1287 :    
1288 :     my $newG = $self->{_genome};
1289 :     my $newGdir = $self->{_orgdir};
1290 :    
1291 :     open(S, "<$newGdir/Subsystems/bindings") or die "Cannot open $newGdir/Subsystems/bindings: $!";
1292 :    
1293 :     my $peg_index;
1294 :     my $bindings;
1295 :    
1296 :     while (<S>)
1297 :     {
1298 :     chomp;
1299 :     my($sname, $role, $peg) = split(/\t/);
1300 :     if (! $self->is_deleted_fid($peg))
1301 :     {
1302 :     push(@{$bindings->{$sname}->{$role}}, $peg);
1303 :     push(@{$peg_index->{$peg}}, [$sname, $role]);
1304 :     }
1305 :    
1306 :     }
1307 :     close(S);
1308 :    
1309 :     open(S, "<$newGdir/Subsystems/subsystems") or die "Cannot open $newGdir/Subsystems/subsystems: $!";
1310 :     my $variant;
1311 :     while (<S>)
1312 :     {
1313 :     chomp;
1314 :     my($sname, $v) = split(/\t/);
1315 :     $variant->{$sname} = $v;
1316 :     }
1317 :     close(S);
1318 :    
1319 :     $self->{_ss_bindings} = $bindings;
1320 :     $self->{_ss_variants} = $variant;
1321 :     $self->{_ss_peg_index} = $peg_index;
1322 :     }
1323 :    
1324 :     sub load_tbl {
1325 :     my($self) = @_;
1326 :    
1327 :     if ($self->{_tbl}) { return };
1328 :    
1329 :     my $newGdir = $self->{_orgdir};
1330 :     my $tbl = {};
1331 :    
1332 :     foreach my $x (`cat $newGdir/Features/*/tbl`)
1333 :     {
1334 :     chomp $x;
1335 :     if ($x =~ /^(\S+)\t(\S+)(\t(\S.*\S))?/)
1336 :     {
1337 :     my $fid = $1;
1338 :     my $loc = [split(/,/,$2)];
1339 :     my $aliases = $4 ? [split(/\t/,$4)] : [];
1340 :     if (! $self->is_deleted_fid($fid))
1341 :     {
1342 :     $tbl->{$fid} = [$loc,$aliases];
1343 :     }
1344 :     }
1345 :     else {
1346 :     warn "Bad feature line in $newGdir:$x:\n";
1347 :     }
1348 :     }
1349 :     print STDERR ("Loaded ", (scalar keys %$tbl), " features from $newGdir\n") if $ENV{FIG_VERBOSE};
1350 :     $self->{_tbl} = $tbl;
1351 :     }
1352 :    
1353 :     sub load_functions {
1354 :     my($self, $force) = @_;
1355 :    
1356 :     if (!$force && $self->{_functions}) { return };
1357 :    
1358 :     my $newG = $self->{_genome};
1359 :     my $newGdir = $self->{_orgdir};
1360 :     my $functions = {};
1361 :     my $roles = {};
1362 :    
1363 :     # order of "cat" is important - proposed_user_functions must be last
1364 :     foreach my $x (`cat $newGdir/*functions`)
1365 :     {
1366 :     if (($x =~ /^(fig\|(\d+\.\d+)\.\S+)\t(\S[^\t]*\S)/) && ($2 eq $newG))
1367 :     {
1368 :     my $peg = $1;
1369 :     my $f = $3;
1370 :     if (! $self->is_deleted_fid($peg))
1371 :     {
1372 :     # user has overridden a function, so delete old roles
1373 :     if (exists($functions->{$peg})) {
1374 :     my @roles = &SeedUtils::roles_of_function($functions->{$peg});
1375 :     foreach $_ (@roles)
1376 :     {
1377 :     delete $roles->{$_};
1378 :     }
1379 :     }
1380 :    
1381 :     # add new roles
1382 :     my @roles = &SeedUtils::roles_of_function($f);
1383 :     foreach $_ (@roles)
1384 :     {
1385 :     push(@{$roles->{$_}},$peg);
1386 :     }
1387 :    
1388 :     # set function
1389 :     $functions->{$peg} = $f;
1390 :     }
1391 :     }
1392 :     }
1393 :     $self->{_functions} = $functions;
1394 :     $self->{_roles} = $roles;
1395 :     }
1396 :    
1397 :     sub seqs_with_role {
1398 :     my($self,$role,$who,$genome) = @_;
1399 :    
1400 :     my $newG = $self->{_genome};
1401 :     if ($genome && $newG && ($genome eq $newG)) {
1402 :     &load_functions($self);
1403 :     my $pegsL = $self->{_roles}->{$role};
1404 :     return ($pegsL ? @{$pegsL} : ());
1405 :     }
1406 :     return ();
1407 :     }
1408 :    
1409 :     sub load_ann {
1410 :     my($self) = @_;
1411 :    
1412 :     if ($self->{_ann}) { return };
1413 :    
1414 :     my $newGdir = $self->{_orgdir};
1415 :     my $ann = {};
1416 :     if (open(ANN,"<$newGdir/annotations"))
1417 :     {
1418 :     $/ = "\n//\n";
1419 :     while (defined(my $x = <ANN>))
1420 :     {
1421 :     chomp $x;
1422 :     if ($x =~ /^(\S+)\n([^\n]+)\n([^\n]+)\n(.*)/s)
1423 :     {
1424 :     push(@{$ann->{$1}},[$1,$2,$3,"$4\n"]);
1425 :     }
1426 :     }
1427 :     $/ = "\n";
1428 :     close(ANN);
1429 :     }
1430 :     $self->{_ann} = $ann;
1431 :     }
1432 :    
1433 :     sub taxonomy_of {
1434 :     my($self) = @_;
1435 :    
1436 :     my $newG = $self->{_genome};
1437 :     my $newGdir = $self->{_orgdir};
1438 :    
1439 :     my $tax;
1440 :     if (open(TAX,"<$newGdir/TAXONOMY") && ($tax = <TAX>))
1441 :     {
1442 :     chop $tax;
1443 :     return $tax;
1444 :     }
1445 :     else
1446 :     {
1447 :     return "unknown";
1448 :     }
1449 :     }
1450 :    
1451 :     sub is_deleted_fid {
1452 :    
1453 :     return 0;
1454 :     }
1455 :    
1456 :     1;

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3