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

Annotation of /FigKernelPackages/SeedV.pm

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1.5 - (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 : olson 1.3 use FileHandle;
27 : overbeek 1.1 use Carp;
28 :     use strict;
29 :     use Data::Dumper;
30 :     use DB_File;
31 :     use File::Basename;
32 : olson 1.3 use SeedUtils;
33 :     use CorrTableEntry;
34 :     use SAPserver;
35 :    
36 :     my $pseed_url = 'http://servers.nmpdr.org/pseed/sapling/server.cgi';
37 : overbeek 1.1
38 :     sub new {
39 :     my ($class, $org_dir ) = @_;
40 :    
41 :     $org_dir || confess("$org_dir must be a valid SEED directory");
42 :     my $self = {};
43 :     $self->{_orgdir} = $org_dir;
44 :     ($org_dir =~ /(\d+\.\d+)$/) || confess("$org_dir must be a path ending in the genome ID");
45 :     $self->{_genome} = $1;
46 : olson 1.3
47 : overbeek 1.1 return bless $self, $class;
48 :     }
49 :    
50 :     sub genome_id
51 :     {
52 :     return $_[0]->{_genome};
53 :     }
54 :    
55 :     # return the path of the organism directory
56 :     sub organism_directory {
57 :     return $_[0]->{_orgdir};
58 :     }
59 :    
60 :     sub get_basic_statistics
61 :     {
62 :     my($self) = @_;
63 :    
64 :     my $newG = $self->{_genome};
65 :     my $newGdir = $self->{_orgdir};
66 :    
67 :     #
68 :     # Check cache.
69 :     #
70 :    
71 :     my $cache = "$newGdir/cache.basic_statistics";
72 :     my $fh = new FileHandle($cache);
73 :     if ($fh)
74 :     {
75 :     my $stats = {};
76 :     while (<$fh>)
77 :     {
78 :     chomp;
79 :     my($k, $v) = split(/\t/);
80 :     $stats->{$k} = $v;
81 :     }
82 :     close($fh);
83 :     return $stats;
84 :     }
85 :    
86 :    
87 :     my $subsystem_data = $self->get_genome_subsystem_data();
88 :    
89 :     my %sscount = map { $_->[0] => 1 } @$subsystem_data;
90 :     my $nss=scalar(keys(%sscount));
91 :    
92 :     my $statistics = {
93 :     num_subsystems => $nss,
94 :     num_contigs => scalar($self->all_contigs()),
95 :     num_basepairs => $self->genome_szdna(),
96 :     genome_name => $self->genus_species(),
97 :     genome_domain => $self->genome_domain(),
98 :     genome_pegs => $self->genome_pegs(),
99 :     genome_rnas => $self->genome_rnas(),
100 :     };
101 :    
102 :     $fh = new FileHandle(">$cache");
103 :     if ($fh)
104 :     {
105 :     while (my($k, $v) = each %$statistics)
106 :     {
107 :     print $fh join("\t", $k, $v), "\n";
108 :     }
109 :     close($fh);
110 :     }
111 :     return $statistics;
112 :     }
113 :    
114 :    
115 :     sub get_peg_statistics {
116 :     my ($self) = @_;
117 :    
118 :     my $newG = $self->{_genome};
119 :     my $newGdir = $self->{_orgdir};
120 :    
121 :     my $cache = "$newGdir/cache.peg_statistics";
122 :     my $fh = new FileHandle($cache);
123 :     if ($fh)
124 :     {
125 :     my $stats = {};
126 :     while (<$fh>)
127 :     {
128 :     chomp;
129 :     my($k, $v) = split(/\t/);
130 :     $stats->{$k} = $v;
131 :     }
132 :     close($fh);
133 :     return $stats;
134 :     }
135 :    
136 :    
137 :     my $subsystem_data = $self->get_genome_subsystem_data();
138 :     my $assignment_data = $self->get_genome_assignment_data();
139 :    
140 :     my $hypo_sub = 0;
141 :     my $hypo_nosub = 0;
142 :     my $nothypo_sub = 0;
143 :     my $nothypo_nosub = 0;
144 :     my %in = map { $_->[2] => 1 } @$subsystem_data;
145 :     my $in = keys(%in);
146 :    
147 :     my %sscount = map { $_->[0] => 1 } @$subsystem_data;
148 :    
149 :     foreach $_ (@$assignment_data)
150 :     {
151 :     my($peg,$func) = @$_;
152 :    
153 :     my $is_hypo = &SeedUtils::hypo($func);
154 :    
155 :     if ($is_hypo && $in{$peg}) { $hypo_sub++ }
156 :     elsif ($is_hypo && ! $in{$peg}) { $hypo_nosub++ }
157 :     elsif ((! $is_hypo) && (! $in{$peg})) { $nothypo_nosub++ }
158 :     elsif ((! $is_hypo) && $in{$peg}) { $nothypo_sub++ }
159 :     }
160 :     my $tot = $hypo_sub + $nothypo_sub + $hypo_nosub + $nothypo_nosub;
161 :    
162 :     my ($fracHS, $fracNHS, $fracHNS, $fracNHNS);
163 :    
164 :     if ($tot == 0) {
165 :     $fracHS = sprintf "%.2f", 0.0;
166 :     $fracNHS = sprintf "%.2f", 0.0;
167 :     $fracHNS = sprintf "%.2f", 0.0;
168 :     $fracNHNS = sprintf "%.2f", 0.0;
169 :     } else {
170 :     $fracHS = sprintf "%.2f", $hypo_sub / $tot * 100;
171 :     $fracNHS = sprintf "%.2f", $nothypo_sub / $tot * 100;
172 :     $fracHNS = sprintf "%.2f", $hypo_nosub / $tot * 100;
173 :     $fracNHNS = sprintf "%.2f", $nothypo_nosub / $tot * 100;
174 :     }
175 :    
176 :     my $statistics = {
177 :     hypothetical_in_subsystem => $hypo_sub,
178 :     hypothetical_not_in_subsystem => $hypo_nosub,
179 :     non_hypothetical_in_subsystem => $nothypo_sub,
180 :     non_hypothetical_not_in_subsystem => $nothypo_nosub,
181 :     hypothetical_in_subsystem_percent => $fracHS,
182 :     hypothetical_not_in_subsystem_percent => $fracHNS,
183 :     non_hypothetical_in_subsystem_percent => $fracNHS,
184 :     non_hypothetical_not_in_subsystem_percent => $fracNHNS
185 :     };
186 :    
187 :     $fh = new FileHandle(">$cache");
188 :     if ($fh)
189 :     {
190 :     while (my($k, $v) = each %$statistics)
191 :     {
192 :     print $fh join("\t", $k, $v), "\n";
193 :     }
194 :     close($fh);
195 :     }
196 :    
197 :     return $statistics;
198 :     }
199 :     sub get_variant_and_bindings
200 :     {
201 :     my($self, $ssa) = @_;
202 :    
203 :     my $newG = $self->{_genome};
204 :     my $newGdir = $self->{_orgdir};
205 :    
206 :     $self->load_ss_data();
207 :    
208 :     my $bindings = $self->{_ss_bindings}->{$ssa};
209 :     my $variant = $self->{_ss_variants}->{$ssa};
210 :    
211 :     unless ($bindings) {
212 :     $variant = '*-1';
213 :     $bindings = {};
214 :     }
215 :    
216 :     return ($variant, $bindings);
217 :     }
218 :    
219 :     sub active_subsystems
220 :     {
221 :     my($self, $all) = @_;
222 :    
223 :     my $newG = $self->{_genome};
224 :     my $newGdir = $self->{_orgdir};
225 :    
226 :     $self->load_ss_data();
227 :    
228 :     my $slist = {};
229 :    
230 :     if ($self->{_ss_variants})
231 :     {
232 :     %{$slist} = %{$self->{_ss_variants}};
233 :     }
234 :    
235 :     if (not $all)
236 :     {
237 :     for my $ss (keys %$slist)
238 :     {
239 :     my $var = $slist->{$ss};
240 :     delete $slist->{$ss} if $var eq 0 or $var eq -1;
241 :     }
242 :     }
243 :     return $slist;
244 :     }
245 :    
246 :     sub peg_to_subsystems
247 :     {
248 :     my($self, $peg) = @_;
249 :    
250 :     my $newG = $self->{_genome};
251 :     my $newGdir = $self->{_orgdir};
252 :    
253 :     $self->load_ss_data();
254 :    
255 :     my $variant = $self->{_ss_variants};
256 :     my %sub = map { $_ => 1 }
257 :     grep { $variant->{$_} !~ /^(-1)|0$/ }
258 :     map { $_->[0] }
259 :     $self->peg_to_roles_in_subsystems($peg);
260 :     return sort keys(%sub);
261 :     }
262 :    
263 :     sub peg_to_roles_in_subsystems
264 :     {
265 :     my($self,$peg) = @_;
266 :    
267 :     my $newG = $self->{_genome};
268 :     my $newGdir = $self->{_orgdir};
269 :    
270 :     $self->load_ss_data();
271 :    
272 :     my $ret = $self->{_ss_peg_index}->{$peg};
273 :    
274 :     return $ret ? @$ret : ();
275 :     }
276 :    
277 :     sub subsystems_for_peg
278 :     {
279 :     my($self,$peg) = @_;
280 :     return $self->peg_to_roles_in_subsystems($peg);
281 :     }
282 :    
283 :     sub subsystems_for_peg_complete {
284 :     my ($self, $peg) = @_;
285 :    
286 :     my $newG = $self->{_genome};
287 :    
288 :     $self->load_ss_data();
289 :    
290 :     my $ret = $self->{_ss_peg_index}->{$peg};
291 :     if ($ret) {
292 :     return map { [ $_->[0], $_->[1], $self->{_ss_variants}->{$_->[0] }, $peg ] } @$ret;
293 :     } else {
294 :     return ();
295 :     }
296 :     }
297 :    
298 :     sub genus_species {
299 :     my($self) = @_;
300 :    
301 :     my $newG = $self->{_genome};
302 :     my $newGdir = $self->{_orgdir};
303 :    
304 :     if (open(GENOME,"<$newGdir/GENOME"))
305 :     {
306 :     my $x = <GENOME>;
307 :     close(GENOME);
308 :     chop $x;
309 :     return $x;
310 :     }
311 :     else
312 :     {
313 :     return "";
314 :     }
315 :     }
316 :    
317 :     sub get_genome_assignment_data {
318 :     my($self) = @_;
319 :    
320 :     my $newG = $self->{_genome};
321 :     my $newGdir = $self->{_orgdir};
322 :    
323 :     my @fnfiles = <$newGdir/proposed*functions>;
324 :     if (@fnfiles == 0)
325 :     {
326 :     @fnfiles = <$newGdir/assigned*functions>;
327 :     }
328 :    
329 :     if (@fnfiles == 0)
330 :     {
331 :     return [];
332 :     }
333 :     my %assign;
334 : olson 1.4 for my $file (@fnfiles)
335 : overbeek 1.1 {
336 : olson 1.4 my $fh;
337 :     if (!open($fh, "<", $file))
338 : overbeek 1.1 {
339 : olson 1.4 warn "Cannot open $file: $!";
340 :     next;
341 :     }
342 :    
343 :     while (<$fh>)
344 :     {
345 :     if ( $_ =~ /^(fig\|\d+\.\d+\.peg\.\d+)\t(\S.*\S)/)
346 :     {
347 :     my($fid,$func) = ($1,$2);
348 :     $assign{$fid} = $func;
349 :     }
350 : overbeek 1.1 }
351 : olson 1.4 close($fh);
352 : overbeek 1.1 }
353 : olson 1.4 return [map { [$_,$assign{$_}] } sort { &SeedUtils::by_fig_id($a,$b) } keys(%assign)];
354 : overbeek 1.1 }
355 :    
356 :     sub get_genome_subsystem_data {
357 :     my($self,$genome) = @_;
358 :    
359 :     my $newG = $self->{_genome};
360 :     my $newGdir = $self->{_orgdir};
361 :    
362 : olson 1.4 my $fh;
363 :     open($fh, "<", "$newGdir/Subsystems/subsystems");
364 :     my %operational = map { (($_ =~ /^(\S.*\S)\t(\S+)/) && (($2 ne '-1') && ($2 ne '0'))) ? ($1 => 1) : () } <$fh>;
365 :     close($fh);
366 :    
367 :     open($fh, "<", "$newGdir/Subsystems/bindings");
368 :     my $rc = [grep { ! $self->is_deleted_fid($_->[2]) }
369 : overbeek 1.1 map { (($_ =~ /^(\S[^\t]+\S)\t(\S[^\t]*\S)\t(\S+)/) && $operational{$1} ) ? [$1,$2,$3] : () }
370 : olson 1.4 <$fh>];
371 :     close($fh);
372 :     return $rc;
373 : overbeek 1.1 }
374 :    
375 :     sub get_genome_subsystem_count
376 :     {
377 :     my($self,$genome) = @_;
378 :    
379 :     my $newG = $self->{_genome};
380 :     my $newGdir = $self->{_orgdir};
381 :    
382 :    
383 :     my $count = 0;
384 :     my ($entry, $vc);
385 :     open(SUBSYSTEMS, "<$newGdir/Subsystems/subsystems");
386 :     while (defined($entry = <SUBSYSTEMS>))
387 :     {
388 :     chomp $entry;
389 :     (undef, $vc) = split /\t/, $entry;
390 :     if ($vc != -1) { ++$count; }
391 :     }
392 :     close(SUBSYSTEMS);
393 :     return $count;
394 :     }
395 :    
396 :     sub orgname_of_orgid {
397 :     my($self,$genome) = @_;
398 :    
399 :     return $self->genus_species();
400 :     }
401 :    
402 :     sub genus_species_domain {
403 :     my($self,$genome) = @_;
404 :    
405 :     return [$self->genus_species(),$self->genome_domain()];
406 :     }
407 :    
408 :     sub protein_subsystem_to_roles {
409 :     my ($self,$peg,$subsystem) = @_;
410 :    
411 :     my $newG = $self->{_genome};
412 :     my $newGdir = $self->{_orgdir};
413 :    
414 : olson 1.4 open(my $fh, "<", "$newGdir/Subsystems/bindings");
415 :    
416 :     my @roles = map { (($_ =~ /^([^\t]+)\t([^\t]+)\t(\S+)$/) && ($1 eq $subsystem) && ($3 eq $peg)) ?
417 :     $2 : () } <$fh>;
418 :     my %roles = map { $_ => 1 } @roles;
419 :     return [sort keys(%roles)];
420 : overbeek 1.1 }
421 :    
422 :     sub contig_lengths {
423 :     my ($self) = @_;
424 :    
425 :     my $newG = $self->{_genome};
426 :     my $newGdir = $self->{_orgdir};
427 :    
428 : olson 1.3 my $contig_lengths = $self->{_contig_len_index};
429 : overbeek 1.1
430 : olson 1.3 # $self->load_contigs_index();
431 :     if (!tied(%$contig_lengths))
432 :     {
433 :     $self->load_contig_len_cache();
434 :     return $self->{_contig_len_cache};
435 :     }
436 :     return $contig_lengths;
437 : overbeek 1.1 }
438 :    
439 :     sub load_contig_len_cache
440 :     {
441 :     my($self) = @_;
442 :    
443 :     return if ref $self->{_contig_len_cache};
444 :    
445 :     my $newGdir = $self->{_orgdir};
446 :    
447 :     my $contig_lengths = {};
448 :     if (open(CONTIGS,"<$newGdir/contigs"))
449 :     {
450 :     local $/ = "\n>";
451 :     while (defined(my $x = <CONTIGS>))
452 :     {
453 :     chomp $x;
454 :     if ($x =~ />?(\S+)[^\n]*\n(.*)/s)
455 :     {
456 :     my $id = $1;
457 :     my $seq = $2;
458 :     $seq =~ s/\s//gs;
459 :     $contig_lengths->{$id} = length($seq);
460 :     }
461 :     }
462 :     close(CONTIGS);
463 :     }
464 :     $self->{_contig_len_cache} = $contig_lengths;
465 :     }
466 :    
467 :     sub contig_ln {
468 :     my ($self, $contig) = @_;
469 :    
470 :     my $newG = $self->{_genome};
471 :     my $newGdir = $self->{_orgdir};
472 :    
473 : olson 1.5 $self->load_contig_len_cache();
474 :     return $self->{_contig_len_cache}->{$contig};
475 : overbeek 1.1 }
476 :    
477 :     sub contigs_of
478 :     {
479 :     my ($self) = @_;
480 :    
481 :     my $newG = $self->{_genome};
482 :     my $newGdir = $self->{_orgdir};
483 :    
484 :     my @out;
485 :     $self->load_contigs_index();
486 :    
487 :     my $contigs = $self->{_contigs_index};
488 :     if (tied(%$contigs))
489 :     {
490 :     return keys %$contigs;
491 :     }
492 :    
493 :     $self->load_contig_len_cache();
494 :    
495 :     return keys %{$self->{_contig_len_cache}};
496 :     }
497 :    
498 :     =head3 dna_seq
499 :    
500 :     usage: $seq = dna_seq(@locations)
501 :    
502 :     Returns the concatenated subsequences described by the list of locations. Each location
503 :     must be of the form
504 :    
505 :     Contig_Beg_End
506 :    
507 :     where Contig must be the ID of a contig for genome $genome. If Beg > End the location
508 :     describes a stretch of the complementary strand.
509 :    
510 :     =cut
511 :     #: Return Type $;
512 :     sub dna_seq {
513 :     my($self,@locations) = @_;
514 :    
515 :     my $newG = $self->{_genome};
516 :     my $newGdir = $self->{_orgdir};
517 :    
518 :     my $contigs = $self->{_contigs_index};
519 :     if (!tied %$contigs)
520 :     {
521 :     $self->load_contig_seq();
522 :     $contigs = $self->{_contigs_seq_cache};
523 :     }
524 :    
525 :     my(@pieces,$loc,$contig,$beg,$end,$ln,$rdbH);
526 :    
527 :     @locations = map { split(/,/,$_) } @locations;
528 :     @pieces = ();
529 :     foreach $loc (@locations)
530 :     {
531 :     if ($loc =~ /^(\S+)_(\d+)_(\d+)$/)
532 :     {
533 :     ($contig,$beg,$end) = ($1,$2,$3);
534 :     my $seq = $contigs->{$contig};
535 :    
536 :     $ln = length($seq);
537 :    
538 :     if (! $ln) {
539 :     print STDERR "$contig: could not get length\n";
540 :     return "";
541 :     }
542 :    
543 :     if (&SeedUtils::between(1,$beg,$ln) && &SeedUtils::between(1,$end,$ln))
544 :     {
545 :     if ($beg < $end)
546 :     {
547 :     push(@pieces, substr($seq, $beg - 1, ($end - $beg) + 1));
548 :     }
549 :     else
550 :     {
551 :     push(@pieces, &SeedUtils::reverse_comp(substr($seq, $end - 1, ($beg - $end) + 1)));
552 :     }
553 :     }
554 :     }
555 :     }
556 :     return lc(join("",@pieces));
557 :     }
558 :    
559 :     sub load_contig_seq
560 :     {
561 :     my($self) = @_;
562 :    
563 :     return if ref($self->{_contigs_seq_cache});
564 :    
565 :     my $newGdir = $self->{_orgdir};
566 :    
567 :     my $contigs = {};
568 :    
569 :     if (open(CONTIGS,"<$newGdir/contigs"))
570 :     {
571 :     local $/ = "\n>";
572 :     while (defined(my $x = <CONTIGS>))
573 :     {
574 :     chomp $x;
575 :     if ($x =~ />?(\S+)[^\n]*\n(.*)/s)
576 :     {
577 :     my $id = $1;
578 :     my $seq = $2;
579 :     $seq =~ s/\s//gs;
580 :     $contigs->{$id} = $seq;
581 :     }
582 :     }
583 :     close(CONTIGS);
584 :     }
585 :     $self->{_contigs_seq_cache} = $contigs;
586 :     }
587 :    
588 :     sub genome_szdna {
589 :     my ($self) = @_;
590 :    
591 :     my $newG = $self->{_genome};
592 :     my $newGdir = $self->{_orgdir};
593 :    
594 :     my $contig_lens = $self->contig_lengths();
595 :     my $tot = 0;
596 :     while ( my($contig,$len) = each %$contig_lens)
597 :     {
598 :     $tot += $len;
599 :     }
600 :     return $tot;
601 :     }
602 :    
603 :     sub genome_pegs {
604 :     my ($self) = @_;
605 :    
606 :     my $newG = $self->{_genome};
607 :     my $newGdir = $self->{_orgdir};
608 :    
609 :     my @tmp = $self->all_features("peg");
610 :     my $n = @tmp;
611 :     return $n;
612 :     }
613 :    
614 :     sub genome_rnas {
615 :     my ($self) = @_;
616 :    
617 :     my $newG = $self->{_genome};
618 :     my $newGdir = $self->{_orgdir};
619 :    
620 :     my @tmp = $self->all_features("rna");
621 :     my $n = @tmp;
622 :     return $n;
623 :     }
624 :    
625 :     sub genome_domain {
626 :     my ($self) = @_;
627 :    
628 :     my $newG = $self->{_genome};
629 :     my $newGdir = $self->{_orgdir};
630 :    
631 :     my $tax = $self->taxonomy_of();
632 :     return ($tax =~ /^([^ \t;]+)/) ? $1 : "unknown";
633 :     }
634 :    
635 :     sub genes_in_region {
636 :     my($self,$contig,$beg,$end) = @_;
637 :    
638 :     my $newG = $self->{_genome};
639 :     my $newGdir = $self->{_orgdir};
640 :    
641 :     $self->load_feature_indexes();
642 :     #
643 :     # Use the recno index if exists.
644 :     #
645 :    
646 :     my $maxV = 0;
647 :     my $minV = 1000000000;
648 :     my $genes = [];
649 :    
650 :     my $recnos = $self->{_feat_recno};
651 :    
652 :     if (ref($recnos))
653 :     {
654 :     while (my($ftype, $list) = each (%$recnos))
655 :     {
656 :     #
657 :     # Look up start/end of this contig in the btree index.
658 :     #
659 :    
660 :     my $inf = $self->{_feat_btree}->{$ftype}->{$contig};
661 :     my($istart, $iend);
662 :     if ($inf)
663 :     {
664 :     ($istart, $iend) = split(/$;/, $inf);
665 :     }
666 :     else
667 :     {
668 :     $istart = 0;
669 :     $iend = $#$list;
670 :     }
671 :    
672 :     for (my $idx = $istart; $idx <= $iend; $idx++)
673 :     {
674 :     my($fid, $fcontig, $fbeg, $fend, $fstrand) = split(/$;/, $list->[$idx]);
675 :     if ($contig eq $fcontig and &overlaps($beg, $end, $fbeg, $fend))
676 :     {
677 :     $minV = &SeedUtils::min($minV,$fbeg,$fend);
678 :     $maxV = &SeedUtils::max($maxV,$fbeg,$fend);
679 :     push(@$genes,$fid);
680 :     }
681 :     }
682 :     }
683 :     }
684 :     else
685 :     {
686 :     &load_tbl($self);
687 :     my $tblH = $self->{_tbl};
688 :     while ( my($fid,$tuple) = each %$tblH)
689 :     {
690 :     if (($tuple->[0]->[0] =~ /^(\S+)_(\d+)_\d+$/) && ($1 eq $contig))
691 :     {
692 :     my $beg1 = $2;
693 :     my $last = @{$tuple->[0]} - 1;
694 :     if (($tuple->[0]->[$last] =~ /^(\S+)_\d+_(\d+)$/) && ($1 eq $contig))
695 :     {
696 :     my $end1 = $2;
697 :     if (&overlaps($beg,$end,$beg1,$end1))
698 :     {
699 :     $minV = &SeedUtils::min($minV,$beg1,$end1);
700 :     $maxV = &SeedUtils::max($maxV,$beg1,$end1);
701 :     push(@$genes,$fid);
702 :     }
703 :     }
704 :     }
705 :     }
706 :     }
707 :     return ($genes,$minV,$maxV);
708 :     }
709 :    
710 :     sub overlaps {
711 :     my($b1,$e1,$b2,$e2) = @_;
712 :    
713 :     if ($b1 > $e1) { ($b1,$e1) = ($e1,$b1) }
714 :     if ($b2 > $e2) { ($b2,$e2) = ($e2,$b2) }
715 :     return &SeedUtils::between($b1,$b2,$e1) || &SeedUtils::between($b2,$b1,$e2);
716 :     }
717 :    
718 :     sub all_contigs {
719 :     my($self) = @_;
720 :    
721 :     my $newG = $self->{_genome};
722 :     my $newGdir = $self->{_orgdir};
723 :    
724 :     $self->load_feature_indexes();
725 :    
726 :     my %contigs;
727 :     my $btrees = $self->{_feat_btree};
728 :    
729 :     if (ref($btrees))
730 :     {
731 :     while (my($ftype, $btree) = each (%$btrees))
732 :     {
733 :     map { $contigs{$_} = 1 } grep { !/^fig/ } keys %$btree ;
734 :     }
735 :     }
736 :     else
737 :     {
738 :     &load_tbl($self);
739 :     my $tblH = $self->{_tbl};
740 :     while ( my($fid,$tuple) = each %$tblH)
741 :     {
742 :     if ($tuple->[0]->[0] =~ /^(\S+)_\d+_\d+$/)
743 :     {
744 :     $contigs{$1} = 1;
745 :     }
746 :     }
747 :     }
748 :     return keys(%contigs);
749 :     }
750 :    
751 :     sub all_features {
752 :     my($self,$type) = @_;
753 :     if (not defined($type)) { $type = qq(); }
754 :    
755 :     my $newG = $self->{_genome};
756 :     my $newGdir = $self->{_orgdir};
757 :    
758 :     #warn "Loading feature indices";
759 :     $self->load_feature_indexes();
760 :    
761 :     my %contigs;
762 :     my $btrees = $self->{_feat_btree};
763 :    
764 :     if (ref($btrees)) {
765 :     warn "B-tree already loaded" if $ENV{FIG_VERBOSE};
766 :    
767 :     my $btree = $btrees->{$type};
768 :     return sort { &SeedUtils::by_fig_id($a, $b) }
769 :     grep { /^fig/ } keys %$btree;
770 :     }
771 :     else {
772 :     warn "Loading contig B-tree" if $ENV{FIG_VERBOSE};
773 :    
774 :     &load_tbl($self);
775 :     my $tblH = $self->{_tbl};
776 :    
777 :     return sort {
778 :     &SeedUtils::by_fig_id($a,$b)
779 :     } grep {
780 :     #...NOTE: Matches all feature types if $type is the null string
781 :     ($_ =~ /^fig\|\d+\.\d+\.([^\.]+)/) && ($1 =~ m/$type/)
782 :     } keys(%$tblH);
783 :     }
784 :     }
785 :    
786 :     sub all_features_detailed_fast {
787 :     my($self,$regmin, $regmax, $contig) = @_;
788 :    
789 :     my $newG = $self->{_genome};
790 :     my $newGdir = $self->{_orgdir};
791 :    
792 :     $self->load_feature_indexes();
793 :     my $feat_details = [];
794 :     my $recnos = $self->{_feat_recno};
795 :    
796 :     if (ref($recnos))
797 :     {
798 :     while (my($ftype, $list) = each (%$recnos))
799 :     {
800 :     #
801 :     # Look up start/end of this contig in the btree index.
802 :     #
803 :    
804 :     my $inf = $self->{_feat_btree}->{$ftype}->{$contig};
805 :     my($istart, $iend);
806 :     if ($inf)
807 :     {
808 :     ($istart, $iend) = split(/$;/, $inf);
809 :     }
810 :     else
811 :     {
812 :     $istart = 0;
813 :     $iend = $#$list;
814 :     }
815 :    
816 :     for (my $idx = $istart; $idx <= $iend; $idx++)
817 :     {
818 :     my($fid, $fcontig, $fbeg, $fend, $fstrand) = split(/$;/, $list->[$idx]);
819 :    
820 :     if (not defined($regmin) or not defined($regmax) or not defined($contig) or
821 :     (($contig eq $fcontig) and
822 :     ($fbeg < $regmin and $regmin < $fend) or ($fbeg < $regmax and $regmax < $fend) or ($fbeg > $regmin and $fend < $regmax)))
823 :     {
824 :     my($loc, $index, @aliases) = split(/$;/, $self->{_feat_btree}->{$ftype}->{$fid});
825 :    
826 :     my $function = $self->function_of($fid);
827 :     push(@$feat_details,[$fid, $loc, join(",", @aliases), $ftype, $fbeg, $fend, $function,'master','']);
828 :     }
829 :     }
830 :     }
831 :     }
832 :     else
833 :     {
834 :     &load_tbl($self);
835 :     my $tblH = $self->{_tbl};
836 :     while ( my($fid,$tuple) = each %$tblH)
837 :     {
838 :     if ($fid =~ /^fig\|\d+\.\d+\.(\S+)\.\d+/)
839 :     {
840 :     my $type = $1;
841 :     if ($tuple->[0]->[0] =~ /^(\S+)_(\d+)_(\d+)$/)
842 :     {
843 :     my($ctg, $beg, $end) = ($1, $2, $3);
844 :     next if (defined($contig) and $contig ne $ctg);
845 :    
846 :     my($min,$max);
847 :     if ($beg < $end)
848 :     {
849 :     $min = $beg;
850 :     $max = $end;
851 :     }
852 :     else
853 :     {
854 :     $min = $end;
855 :     $max = $beg;
856 :     }
857 :    
858 :     if (not defined($regmin) or not defined($regmax) or
859 :     ($min < $regmin and $regmin < $max) or ($min < $regmax and $regmax < $max) or ($min > $regmin and $max < $regmax))
860 :     {
861 :     my $function = $self->function_of($fid);
862 :     push(@$feat_details,[$fid,$tuple->[0]->[0],join(",",@{$tuple->[1]}),$type,$min,$max,$function,'master','']);
863 :     }
864 :     }
865 :     }
866 :     }
867 :     }
868 :     return $feat_details;
869 :     }
870 :    
871 :     sub feature_location {
872 :     my($self,$fid) = @_;
873 :    
874 :     my $newG = $self->{_genome};
875 :     my $newGdir = $self->{_orgdir};
876 :    
877 :     if (($fid =~ /^fig\|(\d+\.\d+)\.([^.]+)/) && ($1 eq $newG))
878 :     {
879 :     my $ftype = $2;
880 :    
881 :     $self->load_feature_indexes();
882 :    
883 : olson 1.5 my $btree = exists($self->{_feat_btree}) ? $self->{_feat_btree}->{$ftype} : undef;
884 : overbeek 1.1 if ($btree)
885 :     {
886 :     my($loc, $idx, @aliases) = split(/$;/, $btree->{$fid});
887 :     return wantarray ? split(/,/, $loc) : $loc;
888 :     }
889 :     else
890 :     {
891 :     &load_tbl($self);
892 :     if (my $x = $self->{_tbl}->{$fid})
893 :     {
894 :     if (wantarray)
895 :     {
896 :     return @{$x->[0]};
897 :     }
898 :     else
899 :     {
900 :     return join(",",@{$x->[0]});
901 :     }
902 :     }
903 :     else
904 :     {
905 :     return undef;
906 :     }
907 :     }
908 :     }
909 :     return undef;
910 :     }
911 :    
912 :     sub function_of {
913 :     my($self,$fid) = @_;
914 :    
915 :     my $newG = $self->{_genome};
916 :     my $newGdir = $self->{_orgdir};
917 :    
918 :     if (($fid =~ /^fig\|(\d+\.\d+)/) && ($1 eq $newG))
919 :     {
920 :     &load_functions($self);
921 :    
922 :     my $fn = $self->{_functions}->{$fid};
923 :     if (wantarray)
924 :     {
925 :     return ['master', $fn];
926 :     }
927 :     else
928 :     {
929 :     return $fn;
930 :     }
931 :     }
932 :     return "";
933 :     }
934 :    
935 : olson 1.4 sub assign_function
936 :     {
937 :     my($self, $fid, $user, $function, $confidence) = @_;
938 :    
939 :     $confidence = $confidence ? $confidence : "";
940 :    
941 :     my $newG = $self->{_genome};
942 :     my $newGdir = $self->{_orgdir};
943 :    
944 :     if (($fid =~ /^fig\|(\d+\.\d+)/) && ($1 ne $newG))
945 :     {
946 :     warn "assign_function on non-seedv fid\n";
947 :     return 0;
948 :     }
949 :    
950 :     $function =~ s/\s+/ /sg; # No multiple spaces
951 :     $function =~ s/^\s+//; # No space at begining
952 :     $function =~ s/\s+$//; # No space at end
953 :     $function =~ s/ ; /; /g; # No space before semicolon
954 :    
955 :     my $file = "$newGdir/proposed_user_functions";
956 :     my $status = 1;
957 :     if ( open( TMP, ">>$file" ) )
958 :     {
959 :     print TMP "$fid\t$function\t$confidence\n";
960 :     close(TMP);
961 :     }
962 :     else
963 :     {
964 :     print STDERR "FAILED ASSIGNMENT: $fid\t$function\t$confidence\n";
965 :     $status = 0;
966 :     }
967 :    
968 :     # mdj: force reload of functions to pick up new assignment
969 :     $self->load_functions(1);
970 :    
971 :     # We are not getting annotations logged. So, we will impose it here.
972 :     $self->add_annotation( $fid, $user, "Set master function to\n$function\n" );
973 :    
974 :     #
975 :     # Mark the genome directory as in need of having bindings recomputed.
976 :     #
977 :     # if (open(S, "<$newGdir/Subsystems/subsystems"))
978 :     # {
979 :     # while (<S>)
980 :     # {
981 :     # chomp;
982 :     # my($sname, $v) = split(/\t/);
983 :     # open(SFILE, ">$self->{_orgdir}/Subsystems/${sname}_bindings_need_recomputation");
984 :     # close(SFILE);
985 :     # }
986 :     # close(S);
987 :     # }
988 :     return $status;
989 :     }
990 :    
991 :     sub add_annotation {
992 :     my($self,$feature_id,$user,$annotation, $time_made) = @_;
993 :    
994 :     my $newG = $self->{_genome};
995 :     my $newGdir = $self->{_orgdir};
996 :    
997 :     if (($feature_id =~ /^fig\|(\d+\.\d+)/) && ($1 ne $newG))
998 :     {
999 :     warn "add_annotation on non-seedv fid\n";
1000 :     return 0;
1001 :     }
1002 :    
1003 :     $time_made = time unless $time_made =~ /^\d+$/;
1004 :    
1005 :     if ($self->is_deleted_fid($feature_id)) { return 0 }
1006 :    
1007 :     # print STDERR "add: fid=$feature_id user=$user annotation=$annotation\n";
1008 :    
1009 :     my $file = "$newGdir/annotations";
1010 :     my $ma = ($annotation =~ /^Set master function to/);
1011 :    
1012 :     if (open(TMP,">>$file"))
1013 :     {
1014 :     my $dataLine = "$feature_id\n$time_made\n$user\n$annotation" . ((substr($annotation,-1) eq "\n") ? "" : "\n");
1015 :     print TMP $dataLine . "//\n";
1016 :     close(TMP);
1017 :    
1018 :     #
1019 :     # Update local cache.
1020 :     #
1021 :     my $ann = $self->{_ann};
1022 :     push(@{$ann->{$feature_id}}, [$feature_id, $time_made, $user, $annotation . "\n"]);
1023 :     }
1024 :     return 0;
1025 :     }
1026 :    
1027 : overbeek 1.1
1028 :     =pod
1029 :    
1030 :     find_features_by_annotation
1031 :    
1032 :     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.
1033 :    
1034 :     If the case boolean is set the search is case insensitive. Otherwise case sensitive.
1035 :    
1036 :     =cut
1037 :    
1038 :     sub find_features_by_annotation {
1039 :     my($self,$anno_hash, $case)=@_;
1040 :     $self->load_functions;
1041 :    
1042 :     if ($case) {map {$anno_hash->{uc($_)}=1} keys %$anno_hash}
1043 :    
1044 :     my $res={};
1045 :     foreach my $id (keys %{$self->{_functions}})
1046 :     {
1047 :     my $fn = $self->{_functions}->{$id};
1048 :     $case ? $fn = uc($fn) : 1;
1049 :     if ($anno_hash->{$fn}) {push @{$res->{$fn}}, $id}
1050 :     }
1051 :    
1052 :     return $res;
1053 :     }
1054 :    
1055 :    
1056 :     =pod
1057 :    
1058 :     search_features_by_annotation
1059 :    
1060 : olson 1.3 Takes a string to find and an optional case boolean, and returns a
1061 :     hash with keys are the function and values are a reference to an array
1062 :     of the IDs that have that function.
1063 : overbeek 1.1
1064 :     If the case boolean is set the search is case insensitive. Otherwise case sensitive.
1065 :    
1066 : olson 1.3 Note that this was originally based on the find above, but this uses a
1067 :     regexp for matching. Will likely be a lot slower which is why I only
1068 :     search for a single term. There may be an SQL way of doing this, if so
1069 :     let Rob know how to do it and I'll replace this method.
1070 : overbeek 1.1
1071 :    
1072 :     =cut
1073 :    
1074 :     sub search_features_by_annotation {
1075 :     my($self,$term, $case)=@_;
1076 :     $self->load_functions;
1077 :    
1078 :     # to make case insensitive convert everything to uppercase
1079 :     # alternative is to use two regexps, one for case insens and one for not case insens
1080 :     # but the bad thing about that approach is that if you have a case insensitive search
1081 :     # you do two regexps for each failed match
1082 :    
1083 :     $case ? $term = uc($term) : 1;
1084 :    
1085 :     my $res={};
1086 :     foreach my $id (keys %{$self->{_functions}})
1087 :     {
1088 :     # we set two variables, one that has case changed for case insensitive searches
1089 :     my $fn = my $fnc = $self->{_functions}->{$id};
1090 :     $case ? $fn = uc($fn) : 1;
1091 :     if ($fn =~ m/$term/) {push @{$res->{$fnc}}, $id}
1092 :     }
1093 :    
1094 :     return $res;
1095 :     }
1096 :    
1097 :    
1098 :     sub feature_aliases {
1099 :     my($self,$fid) = @_;
1100 :    
1101 :     my $newG = $self->{_genome};
1102 :     my $newGdir = $self->{_orgdir};
1103 :    
1104 :     my @aliases;
1105 :    
1106 :     if (($fid =~ /^fig\|(\d+\.\d+)/) && ($1 eq $newG))
1107 :     {
1108 :     &load_tbl($self);
1109 :     if (my $x = $self->{_tbl}->{$fid})
1110 :     {
1111 :     @aliases = @{$x->[1]};
1112 :     }
1113 :     else
1114 :     {
1115 :     @aliases = ();
1116 :     }
1117 :     }
1118 :     else
1119 :     {
1120 :     @aliases = ();
1121 :     }
1122 :     return wantarray() ? @aliases : join(",",@aliases);
1123 :     }
1124 :    
1125 :     sub get_corresponding_ids {
1126 :     my($self, $id, $with_type_info) = @_;
1127 :    
1128 :     my $newG = $self->{_genome};
1129 :    
1130 :     if (($id =~ /^fig\|(\d+\.\d+)/) && ($1 eq $newG)) {
1131 :     my @aliases = $self->feature_aliases($id);
1132 :     my @corresponding_ids = ();
1133 :     foreach my $alias (@aliases) {
1134 :     if ($alias =~ /^gi\|/) {
1135 :     if ($with_type_info) {
1136 :     push(@corresponding_ids, [$alias, 'NCBI']);
1137 :     } else {
1138 :     push(@corresponding_ids, $alias);
1139 :     }
1140 :     last;
1141 :     }
1142 :     }
1143 :     return @corresponding_ids;
1144 :     }
1145 :     }
1146 :    
1147 :     sub feature_annotations {
1148 :     my($self,$fid,$rawtime) = @_;
1149 :    
1150 :     my $newG = $self->{_genome};
1151 :     my $newGdir = $self->{_orgdir};
1152 :    
1153 :     my @annotations;
1154 :     if (($fid =~ /^fig\|(\d+\.\d+)/) && ($1 eq $newG))
1155 :     {
1156 :     &load_ann($self);
1157 :     if (my $x = $self->{_ann}->{$fid})
1158 :     {
1159 :     @annotations = @{$x};
1160 :     }
1161 :     else
1162 :     {
1163 :     @annotations = ();
1164 :     }
1165 :    
1166 :     if ($rawtime)
1167 :     {
1168 :     return @annotations;
1169 :     }
1170 :     else
1171 :     {
1172 :     return map { my $r = [@$_]; $r->[1] = localtime($r->[1]); $r } @annotations;
1173 :     }
1174 :     }
1175 :     return ();
1176 :     }
1177 :    
1178 :     sub get_translation {
1179 :     my($self,$peg) = @_;
1180 :    
1181 :     my $newG = $self->{_genome};
1182 :     my $newGdir = $self->{_orgdir};
1183 :    
1184 :     if (($peg =~ /^fig\|(\d+\.\d+)/) && ($1 eq $newG))
1185 :     {
1186 :     &load_feature_indexes($self);
1187 :    
1188 :     my $out = $self->{_feat}->{peg}->{$peg};
1189 :     if (defined($out))
1190 :     {
1191 :     return $out;
1192 :     }
1193 :    
1194 :     #
1195 :     # If we have a blast-formatted fasta, use fastacmd to
1196 :     # do the lookup, and cache the output for later use.
1197 :     #
1198 :    
1199 :     if ($self->{_feat_fasta}->{peg})
1200 :     {
1201 :     my $id = "gnl|$peg";
1202 :     my $cmd = "$FIG_Config::ext_bin/fastacmd -d $self->{_feat_fasta}->{peg} -s '$id'";
1203 :     open(P, "$cmd|") or die "get_translation: cmd failed with $!: $cmd";
1204 :     $_ = <P>;
1205 :     my $out;
1206 :     while (<P>)
1207 :     {
1208 :     s/\s+//g;
1209 :     $out .= $_;
1210 :     }
1211 :     close(P);
1212 :     $self->{_feat}->{$peg} = $out;
1213 :     return $out;
1214 :     }
1215 :     else
1216 :     {
1217 :     return $self->{_feat}->{$peg};
1218 :     }
1219 :     }
1220 :     else
1221 :     {
1222 :     return undef;
1223 :     }
1224 :     }
1225 :    
1226 :     sub translation_length
1227 :     {
1228 :     my($self, $peg) = @_;
1229 :    
1230 :     my $newG = $self->{_genome};
1231 :     my $newGdir = $self->{_orgdir};
1232 :    
1233 :     if (($peg =~ /^fig\|(\d+\.\d+)/) && ($1 eq $newG))
1234 :     {
1235 :     my $t = $self->get_translation($peg);
1236 :     return length($t);
1237 :     }
1238 :     return undef;
1239 :     }
1240 :    
1241 :     sub load_feature_hash {
1242 :     my($self,$type) = @_;
1243 :    
1244 :     my $newGdir = $self->{_orgdir};
1245 :     my $typeH = {};
1246 :     if (open(FIDS,"<$newGdir/Features/$type/tbl"))
1247 :     {
1248 :     while ($_ = <FIDS>)
1249 :     {
1250 :     if ($_ =~ /^(\S+)/)
1251 :     {
1252 :     my $fid = $1;
1253 :     if (! $self->is_deleted_fid($fid))
1254 :     {
1255 :     $typeH->{$fid} = 1;
1256 :     }
1257 :     }
1258 :     }
1259 :     close(FIDS);
1260 :     }
1261 :     return $typeH;
1262 :     }
1263 :    
1264 :     sub load_feature_indexes
1265 :     {
1266 :     my($self) = @_;
1267 :    
1268 :     if ($self->{_feat}) { return };
1269 :    
1270 :     my $newGdir = $self->{_orgdir};
1271 :    
1272 :     for my $fdir (<$newGdir/Features/*>)
1273 :     {
1274 :     my $ftype = basename($fdir);
1275 :    
1276 :     #
1277 :     # If we have a tbl.btree, tie that for our use.
1278 :     #
1279 :    
1280 :     my $tbl_idx = {};
1281 :     my $tie = tie %$tbl_idx, 'DB_File', "$fdir/tbl.btree", O_RDONLY, 0666, $DB_BTREE;
1282 :     if ($tie)
1283 :     {
1284 :     $self->{_feat_tie}->{$ftype} = $tie;
1285 :     $self->{_feat_btree}->{$ftype} = $tbl_idx;
1286 :     }
1287 :    
1288 :     my $tbl_list = [];
1289 :     my $ltie = tie @$tbl_list, 'DB_File', "$fdir/tbl.recno", O_RDONLY, 0666, $DB_RECNO;
1290 :     if ($tie)
1291 :     {
1292 :     $self->{_feat_ltie}->{$ftype} = $ltie;
1293 :     $self->{_feat_recno}->{$ftype} = $tbl_list;
1294 :     }
1295 :    
1296 :     #
1297 :     # If we have fasta.norm.phr, set _pseq_fasta to the fasta file to use with fastacmd.
1298 :     #
1299 :    
1300 :     my $pseq = {};
1301 :    
1302 :     if (-f "$fdir/fasta.norm.phr")
1303 :     {
1304 :     $self->{_feat_fasta}->{$ftype} = "$fdir/fasta.norm";
1305 :    
1306 :     }
1307 :     else
1308 :     {
1309 :     #
1310 :     # Otherwise, we need to load the data.
1311 :     #
1312 :    
1313 :     if (open(FASTA,"<$newGdir/Features/peg/fasta"))
1314 :     {
1315 :     local $/ = "\n>";
1316 :     my $x;
1317 :     while (defined($x = <FASTA>))
1318 :     {
1319 :     chomp $x;
1320 :     if ($x =~ />?(\S+)[^\n]*\n(.*)/s)
1321 :     {
1322 :     my $peg = $1;
1323 :     my $seq = $2;
1324 :     $seq =~ s/\s//gs;
1325 :     if (! $self->is_deleted_fid($peg))
1326 :     {
1327 :     $pseq->{$peg} = $seq;
1328 :     }
1329 :     }
1330 :     }
1331 :     close(FASTA);
1332 :     }
1333 :     }
1334 :     $self->{_feat}->{$ftype} = $pseq;
1335 :     }
1336 :     }
1337 :    
1338 :     sub load_pseq {
1339 :     my($self) = @_;
1340 :    
1341 :     if ($self->{_pseq}) { return };
1342 :    
1343 :     my $newGdir = $self->{_orgdir};
1344 :     my $fdir = "$newGdir/Features/peg";
1345 :    
1346 :     #
1347 :     # If we have a tbl.btree, tie that for our use.
1348 :     #
1349 :    
1350 :     my $tbl_idx = {};
1351 :     my $tie = tie %$tbl_idx, 'DB_File', "$fdir/tbl.btree", O_RDONLY, 0666, $DB_BTREE;
1352 :     if ($tie)
1353 :     {
1354 :     $self->{_tbl_tie} = $tie;
1355 :     $self->{_tbl_btree} = $tbl_idx;
1356 :     }
1357 :    
1358 :     #
1359 :     # If we have fasta.norm.phr, set _pseq_fasta to the fasta file to use with fastacmd.
1360 :     #
1361 :    
1362 :     my $pseq = {};
1363 :    
1364 :     if (-f "$fdir/fasta.norm.phr")
1365 :     {
1366 :     $self->{_pseq_fasta} = "$fdir/fasta.norm";
1367 :     }
1368 :     else
1369 :     {
1370 :     #
1371 :     # Otherwise, we need to load the data.
1372 :     #
1373 :    
1374 :     if (open(FASTA,"<$newGdir/Features/peg/fasta"))
1375 :     {
1376 :     local $/ = "\n>";
1377 :     my $x;
1378 :     while (defined($x = <FASTA>))
1379 :     {
1380 :     chomp $x;
1381 :     if ($x =~ />?(\S+)[^\n]*\n(.*)/s)
1382 :     {
1383 :     my $peg = $1;
1384 :     my $seq = $2;
1385 :     $seq =~ s/\s//gs;
1386 :     if (! $self->is_deleted_fid($peg))
1387 :     {
1388 :     $pseq->{$peg} = $seq;
1389 :     }
1390 :     }
1391 :     }
1392 :     close(FASTA);
1393 :     }
1394 :     }
1395 :     $self->{_pseq} = $pseq;
1396 :     }
1397 :    
1398 :     sub load_ss_data
1399 :     {
1400 :     my($self, $force) = @_;
1401 :    
1402 :     return if defined($self->{_ss_bindings}) && (! $force);
1403 :    
1404 :     my $newG = $self->{_genome};
1405 :     my $newGdir = $self->{_orgdir};
1406 :    
1407 :     open(S, "<$newGdir/Subsystems/bindings") or die "Cannot open $newGdir/Subsystems/bindings: $!";
1408 :    
1409 :     my $peg_index;
1410 :     my $bindings;
1411 :    
1412 :     while (<S>)
1413 :     {
1414 :     chomp;
1415 :     my($sname, $role, $peg) = split(/\t/);
1416 :     if (! $self->is_deleted_fid($peg))
1417 :     {
1418 :     push(@{$bindings->{$sname}->{$role}}, $peg);
1419 :     push(@{$peg_index->{$peg}}, [$sname, $role]);
1420 :     }
1421 :    
1422 :     }
1423 :     close(S);
1424 :    
1425 :     open(S, "<$newGdir/Subsystems/subsystems") or die "Cannot open $newGdir/Subsystems/subsystems: $!";
1426 :     my $variant;
1427 :     while (<S>)
1428 :     {
1429 :     chomp;
1430 :     my($sname, $v) = split(/\t/);
1431 :     $variant->{$sname} = $v;
1432 :     }
1433 :     close(S);
1434 :    
1435 :     $self->{_ss_bindings} = $bindings;
1436 :     $self->{_ss_variants} = $variant;
1437 :     $self->{_ss_peg_index} = $peg_index;
1438 :     }
1439 :    
1440 :     sub load_tbl {
1441 :     my($self) = @_;
1442 :    
1443 :     if ($self->{_tbl}) { return };
1444 :    
1445 :     my $newGdir = $self->{_orgdir};
1446 :     my $tbl = {};
1447 :    
1448 : olson 1.4 for my $tbl_file (<$newGdir/Features/*/tbl>)
1449 : overbeek 1.1 {
1450 : olson 1.4 if (open(my $fh, "<", $tbl_file))
1451 : overbeek 1.1 {
1452 : olson 1.4 while (defined(my $x = <$fh>))
1453 : overbeek 1.1 {
1454 : olson 1.4 chomp $x;
1455 :     if ($x =~ /^(\S+)\t(\S+)(\t(\S.*\S))?/)
1456 :     {
1457 :     my $fid = $1;
1458 :     my $loc = [split(/,/,$2)];
1459 :     my $aliases = $4 ? [split(/\t/,$4)] : [];
1460 :     if (! $self->is_deleted_fid($fid))
1461 :     {
1462 :     $tbl->{$fid} = [$loc,$aliases];
1463 :     }
1464 :     }
1465 :     else {
1466 :     warn "Bad feature line in $newGdir:$x:\n";
1467 :     }
1468 : overbeek 1.1 }
1469 : olson 1.4 close($fh);
1470 : overbeek 1.1 }
1471 : olson 1.4 else
1472 :     {
1473 :     warn "Cannot open $tbl_file: $!";
1474 : overbeek 1.1 }
1475 :     }
1476 :     print STDERR ("Loaded ", (scalar keys %$tbl), " features from $newGdir\n") if $ENV{FIG_VERBOSE};
1477 :     $self->{_tbl} = $tbl;
1478 :     }
1479 :    
1480 :     sub load_functions {
1481 :     my($self, $force) = @_;
1482 :    
1483 :     if (!$force && $self->{_functions}) { return };
1484 :    
1485 :     my $newG = $self->{_genome};
1486 :     my $newGdir = $self->{_orgdir};
1487 :     my $functions = {};
1488 :     my $roles = {};
1489 :    
1490 :     # order of "cat" is important - proposed_user_functions must be last
1491 : olson 1.4 for my $fn_file (<$newGdir/*functions>)
1492 : overbeek 1.1 {
1493 : olson 1.4 if (open(my $fh, "<", $fn_file))
1494 : overbeek 1.1 {
1495 : olson 1.4 while (defined(my $x = <$fh>))
1496 : overbeek 1.1 {
1497 : olson 1.4 if (($x =~ /^(fig\|(\d+\.\d+)\.\S+)\t(\S[^\t]*\S)/) && ($2 eq $newG))
1498 :     {
1499 :     my $peg = $1;
1500 :     my $f = $3;
1501 :     if (! $self->is_deleted_fid($peg))
1502 : overbeek 1.1 {
1503 : olson 1.4 # user has overridden a function, so delete old roles
1504 :     if (exists($functions->{$peg})) {
1505 :     my @roles = &SeedUtils::roles_of_function($functions->{$peg});
1506 :     foreach $_ (@roles)
1507 :     {
1508 :     delete $roles->{$_};
1509 :     }
1510 :     }
1511 :    
1512 :     # add new roles
1513 :     my @roles = &SeedUtils::roles_of_function($f);
1514 :     foreach $_ (@roles)
1515 :     {
1516 :     push(@{$roles->{$_}},$peg);
1517 :     }
1518 :    
1519 :     # set function
1520 :     $functions->{$peg} = $f;
1521 : overbeek 1.1 }
1522 :     }
1523 :     }
1524 : olson 1.4 close($fh);
1525 :     }
1526 :     else
1527 :     {
1528 :     warn "Cannot open $fn_file: $!";
1529 : overbeek 1.1 }
1530 :     }
1531 :     $self->{_functions} = $functions;
1532 :     $self->{_roles} = $roles;
1533 :     }
1534 :    
1535 :     sub seqs_with_role {
1536 :     my($self,$role,$who,$genome) = @_;
1537 :    
1538 :     my $newG = $self->{_genome};
1539 :     if ($genome && $newG && ($genome eq $newG)) {
1540 :     &load_functions($self);
1541 :     my $pegsL = $self->{_roles}->{$role};
1542 :     return ($pegsL ? @{$pegsL} : ());
1543 :     }
1544 :     return ();
1545 :     }
1546 :    
1547 :     sub load_ann {
1548 :     my($self) = @_;
1549 :    
1550 :     if ($self->{_ann}) { return };
1551 :    
1552 :     my $newGdir = $self->{_orgdir};
1553 :     my $ann = {};
1554 :     if (open(ANN,"<$newGdir/annotations"))
1555 :     {
1556 :     $/ = "\n//\n";
1557 :     while (defined(my $x = <ANN>))
1558 :     {
1559 :     chomp $x;
1560 :     if ($x =~ /^(\S+)\n([^\n]+)\n([^\n]+)\n(.*)/s)
1561 :     {
1562 :     push(@{$ann->{$1}},[$1,$2,$3,"$4\n"]);
1563 :     }
1564 :     }
1565 :     $/ = "\n";
1566 :     close(ANN);
1567 :     }
1568 :     $self->{_ann} = $ann;
1569 :     }
1570 :    
1571 :     sub taxonomy_of {
1572 :     my($self) = @_;
1573 :    
1574 :     my $newG = $self->{_genome};
1575 :     my $newGdir = $self->{_orgdir};
1576 :    
1577 :     my $tax;
1578 :     if (open(TAX,"<$newGdir/TAXONOMY") && ($tax = <TAX>))
1579 :     {
1580 :     chop $tax;
1581 :     return $tax;
1582 :     }
1583 :     else
1584 :     {
1585 :     return "unknown";
1586 :     }
1587 :     }
1588 :    
1589 :     sub is_deleted_fid {
1590 :    
1591 :     return 0;
1592 :     }
1593 :    
1594 : olson 1.3 sub load_correspondences
1595 :     {
1596 :     my($self) = @_;
1597 :    
1598 : olson 1.4 #return if exists($self->{correspondence_index});
1599 : olson 1.3
1600 :     my $dir = $self->{_orgdir};
1601 :     my $index = {};
1602 :    
1603 :     for my $cfile (<$dir/CorrToReferenceGenomes/*>)
1604 :     {
1605 :     if ($cfile =~ m,/\d+\.\d+$,)
1606 :     {
1607 :     if (open(my $fh, "<", $cfile))
1608 :     {
1609 :     while (<$fh>)
1610 :     {
1611 :     my $ent = CorrTableEntry->new($_);
1612 :     push(@{$index->{$ent->id1}}, $ent);
1613 :     }
1614 :     close($fh);
1615 :     }
1616 :     else
1617 :     {
1618 :     warn "Could not open $cfile: $!\n";
1619 :     }
1620 :     }
1621 :     }
1622 :     $self->{correspondence_index} = $index;
1623 :     }
1624 :    
1625 :     sub get_correspondences
1626 :     {
1627 :     my($self, $id) = @_;
1628 :    
1629 :     $self->load_correspondences();
1630 :    
1631 :     return $self->{correspondence_index}->{$id};
1632 :     }
1633 :    
1634 :     sub get_best_correspondences
1635 :     {
1636 :     my($self, $id, $cutoff) = @_;
1637 :     my $corrs = $self->get_correspondences($id);
1638 :    
1639 :     my $lg;
1640 :     my $out = [];
1641 :    
1642 :     for my $hit (sort { SeedUtils::genome_of($a->id2) <=> SeedUtils::genome_of($b->id2) or
1643 :     $a->psc <=> $b->psc }
1644 :     grep { (!defined($cutoff) or $_->psc < $cutoff) and
1645 :     ($_->hitinfo eq '<=>' or $_->func1 eq $_->func2) }
1646 :     @$corrs)
1647 :     {
1648 :     if (defined($lg) && $lg ne SeedUtils::genome_of($hit->id2))
1649 :     {
1650 :     push(@$out, $hit);
1651 :     }
1652 :     $lg = SeedUtils::genome_of($hit->id2);
1653 :     }
1654 :     return $out;
1655 :    
1656 :     }
1657 :    
1658 :     sub get_pin
1659 :     {
1660 :     my($self, $peg, $n, $cutoff) = @_;
1661 :    
1662 :     my $hits = $self->get_best_correspondences($peg, $cutoff);
1663 :     # print join("\t", $_->id2, $_->func2) . "\n" for @$hits;
1664 :     my @pegs = map { $_->id2 }
1665 :     sort { $b->bsc <=> $a->bsc or
1666 :     abs($a->len1 - $a->len2) <=> abs($b->len1 - $b->len2) }
1667 :     @$hits;
1668 :     $#pegs = ($n - 1) if @pegs > $n;
1669 :     return \@pegs;
1670 :     }
1671 :    
1672 :     sub get_context
1673 :     {
1674 :     my($self, $peg, $pin, $width) = @_;
1675 :    
1676 :     #
1677 :     # Determine local context.
1678 :     #
1679 :    
1680 :     my $peg_loc = $self->feature_location($peg);
1681 :     my ($contig, $min, $max, $dir) = &SeedUtils::boundaries_of($peg_loc);
1682 :     my $center = ($min + $max) / 2;
1683 :    
1684 :     my ($my_genes, $minV, $maxV) = $self->genes_in_region($contig,
1685 :     int($center - $width / 2),
1686 :     int($center + $width / 2));
1687 :    
1688 :     my $left_extent = $center - $minV;
1689 :     my $right_extent = $maxV - $center;
1690 :    
1691 :     #
1692 :     # Determine other context. Get locations for the pin, then find the
1693 :     # regions.
1694 :     #
1695 :     my $sap = $self->{_sap};
1696 : olson 1.4 if (!defined($sap))
1697 :     {
1698 :     $sap = $self->{_sap} = SAPserver->new($pseed_url);
1699 :     }
1700 :    
1701 :     my $genome_names = $sap->genome_names(-ids => [map { &SeedUtils::genome_of($_) } @$pin]);
1702 :     my $my_genome = $self->genus_species;
1703 :     $my_genome = $self->genome_id if $my_genome eq '';
1704 :    
1705 :     $genome_names->{SeedUtils::genome_of($peg)} = $my_genome;
1706 : olson 1.3
1707 :     my $locs = $sap->fid_locations(-ids => $pin, -boundaries => 1);
1708 :    
1709 :     my @extended_locs;
1710 :     for my $lpeg (keys %$locs)
1711 :     {
1712 :     my $loc = $locs->{$lpeg};
1713 :     my($lcontig, $lmin, $lmax, $ldir) = &SeedUtils::boundaries_of($loc);
1714 :     my $center = ($lmin + $lmax) / 2;
1715 :     my $left_end = int($center - $left_extent);
1716 :     $left_end = 0 if $left_end < 0;
1717 :     my $right_end = int($center + $right_extent);
1718 :     push(@extended_locs, &SeedUtils::location_string($lcontig, $left_end, $right_end));
1719 :     }
1720 :    
1721 :     my $regions = $sap->genes_in_region(-locations => \@extended_locs, -includeLocation => 1);
1722 :    
1723 :     #
1724 :     # Determine functions.
1725 :     #
1726 :     my @ext_pegs = map { keys %$_ } values %$regions;
1727 :     my $funcs = $sap->ids_to_functions(-ids => \@ext_pegs);
1728 :    
1729 :     #
1730 :     # Overlay with local functions.
1731 :    
1732 :     $funcs->{$_} = $self->function_of($_) for @$my_genes;
1733 :    
1734 :     #
1735 :     # We have the whole shebang now. We can assemble and return.
1736 :     #
1737 :     # Each genome is a list of genes.
1738 :     # Each gene is a list [fid, contig, beg, end, dir, func].
1739 :     #
1740 :    
1741 :     my @context;
1742 :     my $cref = {};
1743 :     my $row = 0;
1744 :    
1745 :     #
1746 :     # Start with local context.
1747 :     #
1748 :    
1749 :     my @loc = map { my($contig, $min, $max, $dir) = &SeedUtils::boundaries_of($self->feature_location($_));
1750 :     [$_, $contig, $min, $max, $dir, $funcs->{$_}, $row] } @$my_genes;
1751 :     push @context, [ sort { $a->[2] <=> $b->[2] } @loc ];
1752 :     $cref->{$_->[0]} = $_ for @loc;
1753 :    
1754 :     $row++;
1755 :     #
1756 :     # And include the pinned region.
1757 :     #
1758 :     for my $loc (@extended_locs)
1759 :     {
1760 :     my $region = $regions->{$loc};
1761 :     my @row;
1762 :     for my $peg (keys %$region)
1763 :     {
1764 :     my($contig, $min, $max, $dir) = &SeedUtils::boundaries_of($region->{$peg});
1765 :     my $ent = [$peg, $contig, $min, $max, $dir, $funcs->{$peg}, $row];
1766 :     $cref->{$peg} = $ent;
1767 :     push(@row, $ent);
1768 :     }
1769 :     $row++;
1770 :    
1771 :     push @context, [ sort { $a->[2] <=> $b->[2] } @row ];
1772 :     }
1773 :    
1774 : olson 1.4 return \@context, $cref, $genome_names;
1775 : olson 1.3 }
1776 :    
1777 :    
1778 : overbeek 1.1 1;

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3