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

Annotation of /FigKernelPackages/SeedV.pm

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1.4 - (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 :     my $contig_len = $self->{_contig_len_index};
474 :     $self->load_contigs_index();
475 :     if (tied(%$contig_len))
476 :     {
477 :     return $contig_len->{$contig};
478 :     }
479 :     $self->load_contig_len_cache();
480 :     return $self->{_contig_len_cache}->{$contig};
481 :     }
482 :    
483 :     sub contigs_of
484 :     {
485 :     my ($self) = @_;
486 :    
487 :     my $newG = $self->{_genome};
488 :     my $newGdir = $self->{_orgdir};
489 :    
490 :     my @out;
491 :     $self->load_contigs_index();
492 :    
493 :     my $contigs = $self->{_contigs_index};
494 :     if (tied(%$contigs))
495 :     {
496 :     return keys %$contigs;
497 :     }
498 :    
499 :     $self->load_contig_len_cache();
500 :    
501 :     return keys %{$self->{_contig_len_cache}};
502 :     }
503 :    
504 :     =head3 dna_seq
505 :    
506 :     usage: $seq = dna_seq(@locations)
507 :    
508 :     Returns the concatenated subsequences described by the list of locations. Each location
509 :     must be of the form
510 :    
511 :     Contig_Beg_End
512 :    
513 :     where Contig must be the ID of a contig for genome $genome. If Beg > End the location
514 :     describes a stretch of the complementary strand.
515 :    
516 :     =cut
517 :     #: Return Type $;
518 :     sub dna_seq {
519 :     my($self,@locations) = @_;
520 :    
521 :     my $newG = $self->{_genome};
522 :     my $newGdir = $self->{_orgdir};
523 :    
524 :     my $contigs = $self->{_contigs_index};
525 :     if (!tied %$contigs)
526 :     {
527 :     $self->load_contig_seq();
528 :     $contigs = $self->{_contigs_seq_cache};
529 :     }
530 :    
531 :     my(@pieces,$loc,$contig,$beg,$end,$ln,$rdbH);
532 :    
533 :     @locations = map { split(/,/,$_) } @locations;
534 :     @pieces = ();
535 :     foreach $loc (@locations)
536 :     {
537 :     if ($loc =~ /^(\S+)_(\d+)_(\d+)$/)
538 :     {
539 :     ($contig,$beg,$end) = ($1,$2,$3);
540 :     my $seq = $contigs->{$contig};
541 :    
542 :     $ln = length($seq);
543 :    
544 :     if (! $ln) {
545 :     print STDERR "$contig: could not get length\n";
546 :     return "";
547 :     }
548 :    
549 :     if (&SeedUtils::between(1,$beg,$ln) && &SeedUtils::between(1,$end,$ln))
550 :     {
551 :     if ($beg < $end)
552 :     {
553 :     push(@pieces, substr($seq, $beg - 1, ($end - $beg) + 1));
554 :     }
555 :     else
556 :     {
557 :     push(@pieces, &SeedUtils::reverse_comp(substr($seq, $end - 1, ($beg - $end) + 1)));
558 :     }
559 :     }
560 :     }
561 :     }
562 :     return lc(join("",@pieces));
563 :     }
564 :    
565 :     sub load_contig_seq
566 :     {
567 :     my($self) = @_;
568 :    
569 :     return if ref($self->{_contigs_seq_cache});
570 :    
571 :     my $newGdir = $self->{_orgdir};
572 :    
573 :     my $contigs = {};
574 :    
575 :     if (open(CONTIGS,"<$newGdir/contigs"))
576 :     {
577 :     local $/ = "\n>";
578 :     while (defined(my $x = <CONTIGS>))
579 :     {
580 :     chomp $x;
581 :     if ($x =~ />?(\S+)[^\n]*\n(.*)/s)
582 :     {
583 :     my $id = $1;
584 :     my $seq = $2;
585 :     $seq =~ s/\s//gs;
586 :     $contigs->{$id} = $seq;
587 :     }
588 :     }
589 :     close(CONTIGS);
590 :     }
591 :     $self->{_contigs_seq_cache} = $contigs;
592 :     }
593 :    
594 :     sub genome_szdna {
595 :     my ($self) = @_;
596 :    
597 :     my $newG = $self->{_genome};
598 :     my $newGdir = $self->{_orgdir};
599 :    
600 :     my $contig_lens = $self->contig_lengths();
601 :     my $tot = 0;
602 :     while ( my($contig,$len) = each %$contig_lens)
603 :     {
604 :     $tot += $len;
605 :     }
606 :     return $tot;
607 :     }
608 :    
609 :     sub genome_pegs {
610 :     my ($self) = @_;
611 :    
612 :     my $newG = $self->{_genome};
613 :     my $newGdir = $self->{_orgdir};
614 :    
615 :     my @tmp = $self->all_features("peg");
616 :     my $n = @tmp;
617 :     return $n;
618 :     }
619 :    
620 :     sub genome_rnas {
621 :     my ($self) = @_;
622 :    
623 :     my $newG = $self->{_genome};
624 :     my $newGdir = $self->{_orgdir};
625 :    
626 :     my @tmp = $self->all_features("rna");
627 :     my $n = @tmp;
628 :     return $n;
629 :     }
630 :    
631 :     sub genome_domain {
632 :     my ($self) = @_;
633 :    
634 :     my $newG = $self->{_genome};
635 :     my $newGdir = $self->{_orgdir};
636 :    
637 :     my $tax = $self->taxonomy_of();
638 :     return ($tax =~ /^([^ \t;]+)/) ? $1 : "unknown";
639 :     }
640 :    
641 :     sub genes_in_region {
642 :     my($self,$contig,$beg,$end) = @_;
643 :    
644 :     my $newG = $self->{_genome};
645 :     my $newGdir = $self->{_orgdir};
646 :    
647 :     $self->load_feature_indexes();
648 :     #
649 :     # Use the recno index if exists.
650 :     #
651 :    
652 :     my $maxV = 0;
653 :     my $minV = 1000000000;
654 :     my $genes = [];
655 :    
656 :     my $recnos = $self->{_feat_recno};
657 :    
658 :     if (ref($recnos))
659 :     {
660 :     while (my($ftype, $list) = each (%$recnos))
661 :     {
662 :     #
663 :     # Look up start/end of this contig in the btree index.
664 :     #
665 :    
666 :     my $inf = $self->{_feat_btree}->{$ftype}->{$contig};
667 :     my($istart, $iend);
668 :     if ($inf)
669 :     {
670 :     ($istart, $iend) = split(/$;/, $inf);
671 :     }
672 :     else
673 :     {
674 :     $istart = 0;
675 :     $iend = $#$list;
676 :     }
677 :    
678 :     for (my $idx = $istart; $idx <= $iend; $idx++)
679 :     {
680 :     my($fid, $fcontig, $fbeg, $fend, $fstrand) = split(/$;/, $list->[$idx]);
681 :     if ($contig eq $fcontig and &overlaps($beg, $end, $fbeg, $fend))
682 :     {
683 :     $minV = &SeedUtils::min($minV,$fbeg,$fend);
684 :     $maxV = &SeedUtils::max($maxV,$fbeg,$fend);
685 :     push(@$genes,$fid);
686 :     }
687 :     }
688 :     }
689 :     }
690 :     else
691 :     {
692 :     &load_tbl($self);
693 :     my $tblH = $self->{_tbl};
694 :     while ( my($fid,$tuple) = each %$tblH)
695 :     {
696 :     if (($tuple->[0]->[0] =~ /^(\S+)_(\d+)_\d+$/) && ($1 eq $contig))
697 :     {
698 :     my $beg1 = $2;
699 :     my $last = @{$tuple->[0]} - 1;
700 :     if (($tuple->[0]->[$last] =~ /^(\S+)_\d+_(\d+)$/) && ($1 eq $contig))
701 :     {
702 :     my $end1 = $2;
703 :     if (&overlaps($beg,$end,$beg1,$end1))
704 :     {
705 :     $minV = &SeedUtils::min($minV,$beg1,$end1);
706 :     $maxV = &SeedUtils::max($maxV,$beg1,$end1);
707 :     push(@$genes,$fid);
708 :     }
709 :     }
710 :     }
711 :     }
712 :     }
713 :     return ($genes,$minV,$maxV);
714 :     }
715 :    
716 :     sub overlaps {
717 :     my($b1,$e1,$b2,$e2) = @_;
718 :    
719 :     if ($b1 > $e1) { ($b1,$e1) = ($e1,$b1) }
720 :     if ($b2 > $e2) { ($b2,$e2) = ($e2,$b2) }
721 :     return &SeedUtils::between($b1,$b2,$e1) || &SeedUtils::between($b2,$b1,$e2);
722 :     }
723 :    
724 :     sub all_contigs {
725 :     my($self) = @_;
726 :    
727 :     my $newG = $self->{_genome};
728 :     my $newGdir = $self->{_orgdir};
729 :    
730 :     $self->load_feature_indexes();
731 :    
732 :     my %contigs;
733 :     my $btrees = $self->{_feat_btree};
734 :    
735 :     if (ref($btrees))
736 :     {
737 :     while (my($ftype, $btree) = each (%$btrees))
738 :     {
739 :     map { $contigs{$_} = 1 } grep { !/^fig/ } keys %$btree ;
740 :     }
741 :     }
742 :     else
743 :     {
744 :     &load_tbl($self);
745 :     my $tblH = $self->{_tbl};
746 :     while ( my($fid,$tuple) = each %$tblH)
747 :     {
748 :     if ($tuple->[0]->[0] =~ /^(\S+)_\d+_\d+$/)
749 :     {
750 :     $contigs{$1} = 1;
751 :     }
752 :     }
753 :     }
754 :     return keys(%contigs);
755 :     }
756 :    
757 :     sub all_features {
758 :     my($self,$type) = @_;
759 :     if (not defined($type)) { $type = qq(); }
760 :    
761 :     my $newG = $self->{_genome};
762 :     my $newGdir = $self->{_orgdir};
763 :    
764 :     #warn "Loading feature indices";
765 :     $self->load_feature_indexes();
766 :    
767 :     my %contigs;
768 :     my $btrees = $self->{_feat_btree};
769 :    
770 :     if (ref($btrees)) {
771 :     warn "B-tree already loaded" if $ENV{FIG_VERBOSE};
772 :    
773 :     my $btree = $btrees->{$type};
774 :     return sort { &SeedUtils::by_fig_id($a, $b) }
775 :     grep { /^fig/ } keys %$btree;
776 :     }
777 :     else {
778 :     warn "Loading contig B-tree" if $ENV{FIG_VERBOSE};
779 :    
780 :     &load_tbl($self);
781 :     my $tblH = $self->{_tbl};
782 :    
783 :     return sort {
784 :     &SeedUtils::by_fig_id($a,$b)
785 :     } grep {
786 :     #...NOTE: Matches all feature types if $type is the null string
787 :     ($_ =~ /^fig\|\d+\.\d+\.([^\.]+)/) && ($1 =~ m/$type/)
788 :     } keys(%$tblH);
789 :     }
790 :     }
791 :    
792 :     sub all_features_detailed_fast {
793 :     my($self,$regmin, $regmax, $contig) = @_;
794 :    
795 :     my $newG = $self->{_genome};
796 :     my $newGdir = $self->{_orgdir};
797 :    
798 :     $self->load_feature_indexes();
799 :     my $feat_details = [];
800 :     my $recnos = $self->{_feat_recno};
801 :    
802 :     if (ref($recnos))
803 :     {
804 :     while (my($ftype, $list) = each (%$recnos))
805 :     {
806 :     #
807 :     # Look up start/end of this contig in the btree index.
808 :     #
809 :    
810 :     my $inf = $self->{_feat_btree}->{$ftype}->{$contig};
811 :     my($istart, $iend);
812 :     if ($inf)
813 :     {
814 :     ($istart, $iend) = split(/$;/, $inf);
815 :     }
816 :     else
817 :     {
818 :     $istart = 0;
819 :     $iend = $#$list;
820 :     }
821 :    
822 :     for (my $idx = $istart; $idx <= $iend; $idx++)
823 :     {
824 :     my($fid, $fcontig, $fbeg, $fend, $fstrand) = split(/$;/, $list->[$idx]);
825 :    
826 :     if (not defined($regmin) or not defined($regmax) or not defined($contig) or
827 :     (($contig eq $fcontig) and
828 :     ($fbeg < $regmin and $regmin < $fend) or ($fbeg < $regmax and $regmax < $fend) or ($fbeg > $regmin and $fend < $regmax)))
829 :     {
830 :     my($loc, $index, @aliases) = split(/$;/, $self->{_feat_btree}->{$ftype}->{$fid});
831 :    
832 :     my $function = $self->function_of($fid);
833 :     push(@$feat_details,[$fid, $loc, join(",", @aliases), $ftype, $fbeg, $fend, $function,'master','']);
834 :     }
835 :     }
836 :     }
837 :     }
838 :     else
839 :     {
840 :     &load_tbl($self);
841 :     my $tblH = $self->{_tbl};
842 :     while ( my($fid,$tuple) = each %$tblH)
843 :     {
844 :     if ($fid =~ /^fig\|\d+\.\d+\.(\S+)\.\d+/)
845 :     {
846 :     my $type = $1;
847 :     if ($tuple->[0]->[0] =~ /^(\S+)_(\d+)_(\d+)$/)
848 :     {
849 :     my($ctg, $beg, $end) = ($1, $2, $3);
850 :     next if (defined($contig) and $contig ne $ctg);
851 :    
852 :     my($min,$max);
853 :     if ($beg < $end)
854 :     {
855 :     $min = $beg;
856 :     $max = $end;
857 :     }
858 :     else
859 :     {
860 :     $min = $end;
861 :     $max = $beg;
862 :     }
863 :    
864 :     if (not defined($regmin) or not defined($regmax) or
865 :     ($min < $regmin and $regmin < $max) or ($min < $regmax and $regmax < $max) or ($min > $regmin and $max < $regmax))
866 :     {
867 :     my $function = $self->function_of($fid);
868 :     push(@$feat_details,[$fid,$tuple->[0]->[0],join(",",@{$tuple->[1]}),$type,$min,$max,$function,'master','']);
869 :     }
870 :     }
871 :     }
872 :     }
873 :     }
874 :     return $feat_details;
875 :     }
876 :    
877 :     sub feature_location {
878 :     my($self,$fid) = @_;
879 :    
880 :     my $newG = $self->{_genome};
881 :     my $newGdir = $self->{_orgdir};
882 :    
883 :     if (($fid =~ /^fig\|(\d+\.\d+)\.([^.]+)/) && ($1 eq $newG))
884 :     {
885 :     my $ftype = $2;
886 :    
887 :     $self->load_feature_indexes();
888 :    
889 :     my $btree = $self->{_feat_btree}->{$ftype};
890 :     if ($btree)
891 :     {
892 :     my($loc, $idx, @aliases) = split(/$;/, $btree->{$fid});
893 :     return wantarray ? split(/,/, $loc) : $loc;
894 :     }
895 :     else
896 :     {
897 :     &load_tbl($self);
898 :     if (my $x = $self->{_tbl}->{$fid})
899 :     {
900 :     if (wantarray)
901 :     {
902 :     return @{$x->[0]};
903 :     }
904 :     else
905 :     {
906 :     return join(",",@{$x->[0]});
907 :     }
908 :     }
909 :     else
910 :     {
911 :     return undef;
912 :     }
913 :     }
914 :     }
915 :     return undef;
916 :     }
917 :    
918 :     sub function_of {
919 :     my($self,$fid) = @_;
920 :    
921 :     my $newG = $self->{_genome};
922 :     my $newGdir = $self->{_orgdir};
923 :    
924 :     if (($fid =~ /^fig\|(\d+\.\d+)/) && ($1 eq $newG))
925 :     {
926 :     &load_functions($self);
927 :    
928 :     my $fn = $self->{_functions}->{$fid};
929 :     if (wantarray)
930 :     {
931 :     return ['master', $fn];
932 :     }
933 :     else
934 :     {
935 :     return $fn;
936 :     }
937 :     }
938 :     return "";
939 :     }
940 :    
941 : olson 1.4 sub assign_function
942 :     {
943 :     my($self, $fid, $user, $function, $confidence) = @_;
944 :    
945 :     $confidence = $confidence ? $confidence : "";
946 :    
947 :     my $newG = $self->{_genome};
948 :     my $newGdir = $self->{_orgdir};
949 :    
950 :     if (($fid =~ /^fig\|(\d+\.\d+)/) && ($1 ne $newG))
951 :     {
952 :     warn "assign_function on non-seedv fid\n";
953 :     return 0;
954 :     }
955 :    
956 :     $function =~ s/\s+/ /sg; # No multiple spaces
957 :     $function =~ s/^\s+//; # No space at begining
958 :     $function =~ s/\s+$//; # No space at end
959 :     $function =~ s/ ; /; /g; # No space before semicolon
960 :    
961 :     my $file = "$newGdir/proposed_user_functions";
962 :     my $status = 1;
963 :     if ( open( TMP, ">>$file" ) )
964 :     {
965 :     print TMP "$fid\t$function\t$confidence\n";
966 :     close(TMP);
967 :     }
968 :     else
969 :     {
970 :     print STDERR "FAILED ASSIGNMENT: $fid\t$function\t$confidence\n";
971 :     $status = 0;
972 :     }
973 :    
974 :     # mdj: force reload of functions to pick up new assignment
975 :     $self->load_functions(1);
976 :    
977 :     # We are not getting annotations logged. So, we will impose it here.
978 :     $self->add_annotation( $fid, $user, "Set master function to\n$function\n" );
979 :    
980 :     #
981 :     # Mark the genome directory as in need of having bindings recomputed.
982 :     #
983 :     # if (open(S, "<$newGdir/Subsystems/subsystems"))
984 :     # {
985 :     # while (<S>)
986 :     # {
987 :     # chomp;
988 :     # my($sname, $v) = split(/\t/);
989 :     # open(SFILE, ">$self->{_orgdir}/Subsystems/${sname}_bindings_need_recomputation");
990 :     # close(SFILE);
991 :     # }
992 :     # close(S);
993 :     # }
994 :     return $status;
995 :     }
996 :    
997 :     sub add_annotation {
998 :     my($self,$feature_id,$user,$annotation, $time_made) = @_;
999 :    
1000 :     my $newG = $self->{_genome};
1001 :     my $newGdir = $self->{_orgdir};
1002 :    
1003 :     if (($feature_id =~ /^fig\|(\d+\.\d+)/) && ($1 ne $newG))
1004 :     {
1005 :     warn "add_annotation on non-seedv fid\n";
1006 :     return 0;
1007 :     }
1008 :    
1009 :     $time_made = time unless $time_made =~ /^\d+$/;
1010 :    
1011 :     if ($self->is_deleted_fid($feature_id)) { return 0 }
1012 :    
1013 :     # print STDERR "add: fid=$feature_id user=$user annotation=$annotation\n";
1014 :    
1015 :     my $file = "$newGdir/annotations";
1016 :     my $ma = ($annotation =~ /^Set master function to/);
1017 :    
1018 :     if (open(TMP,">>$file"))
1019 :     {
1020 :     my $dataLine = "$feature_id\n$time_made\n$user\n$annotation" . ((substr($annotation,-1) eq "\n") ? "" : "\n");
1021 :     print TMP $dataLine . "//\n";
1022 :     close(TMP);
1023 :    
1024 :     #
1025 :     # Update local cache.
1026 :     #
1027 :     my $ann = $self->{_ann};
1028 :     push(@{$ann->{$feature_id}}, [$feature_id, $time_made, $user, $annotation . "\n"]);
1029 :     }
1030 :     return 0;
1031 :     }
1032 :    
1033 : overbeek 1.1
1034 :     =pod
1035 :    
1036 :     find_features_by_annotation
1037 :    
1038 :     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.
1039 :    
1040 :     If the case boolean is set the search is case insensitive. Otherwise case sensitive.
1041 :    
1042 :     =cut
1043 :    
1044 :     sub find_features_by_annotation {
1045 :     my($self,$anno_hash, $case)=@_;
1046 :     $self->load_functions;
1047 :    
1048 :     if ($case) {map {$anno_hash->{uc($_)}=1} keys %$anno_hash}
1049 :    
1050 :     my $res={};
1051 :     foreach my $id (keys %{$self->{_functions}})
1052 :     {
1053 :     my $fn = $self->{_functions}->{$id};
1054 :     $case ? $fn = uc($fn) : 1;
1055 :     if ($anno_hash->{$fn}) {push @{$res->{$fn}}, $id}
1056 :     }
1057 :    
1058 :     return $res;
1059 :     }
1060 :    
1061 :    
1062 :     =pod
1063 :    
1064 :     search_features_by_annotation
1065 :    
1066 : olson 1.3 Takes a string to find and an optional case boolean, and returns a
1067 :     hash with keys are the function and values are a reference to an array
1068 :     of the IDs that have that function.
1069 : overbeek 1.1
1070 :     If the case boolean is set the search is case insensitive. Otherwise case sensitive.
1071 :    
1072 : olson 1.3 Note that this was originally based on the find above, but this uses a
1073 :     regexp for matching. Will likely be a lot slower which is why I only
1074 :     search for a single term. There may be an SQL way of doing this, if so
1075 :     let Rob know how to do it and I'll replace this method.
1076 : overbeek 1.1
1077 :    
1078 :     =cut
1079 :    
1080 :     sub search_features_by_annotation {
1081 :     my($self,$term, $case)=@_;
1082 :     $self->load_functions;
1083 :    
1084 :     # to make case insensitive convert everything to uppercase
1085 :     # alternative is to use two regexps, one for case insens and one for not case insens
1086 :     # but the bad thing about that approach is that if you have a case insensitive search
1087 :     # you do two regexps for each failed match
1088 :    
1089 :     $case ? $term = uc($term) : 1;
1090 :    
1091 :     my $res={};
1092 :     foreach my $id (keys %{$self->{_functions}})
1093 :     {
1094 :     # we set two variables, one that has case changed for case insensitive searches
1095 :     my $fn = my $fnc = $self->{_functions}->{$id};
1096 :     $case ? $fn = uc($fn) : 1;
1097 :     if ($fn =~ m/$term/) {push @{$res->{$fnc}}, $id}
1098 :     }
1099 :    
1100 :     return $res;
1101 :     }
1102 :    
1103 :    
1104 :     sub feature_aliases {
1105 :     my($self,$fid) = @_;
1106 :    
1107 :     my $newG = $self->{_genome};
1108 :     my $newGdir = $self->{_orgdir};
1109 :    
1110 :     my @aliases;
1111 :    
1112 :     if (($fid =~ /^fig\|(\d+\.\d+)/) && ($1 eq $newG))
1113 :     {
1114 :     &load_tbl($self);
1115 :     if (my $x = $self->{_tbl}->{$fid})
1116 :     {
1117 :     @aliases = @{$x->[1]};
1118 :     }
1119 :     else
1120 :     {
1121 :     @aliases = ();
1122 :     }
1123 :     }
1124 :     else
1125 :     {
1126 :     @aliases = ();
1127 :     }
1128 :     return wantarray() ? @aliases : join(",",@aliases);
1129 :     }
1130 :    
1131 :     sub get_corresponding_ids {
1132 :     my($self, $id, $with_type_info) = @_;
1133 :    
1134 :     my $newG = $self->{_genome};
1135 :    
1136 :     if (($id =~ /^fig\|(\d+\.\d+)/) && ($1 eq $newG)) {
1137 :     my @aliases = $self->feature_aliases($id);
1138 :     my @corresponding_ids = ();
1139 :     foreach my $alias (@aliases) {
1140 :     if ($alias =~ /^gi\|/) {
1141 :     if ($with_type_info) {
1142 :     push(@corresponding_ids, [$alias, 'NCBI']);
1143 :     } else {
1144 :     push(@corresponding_ids, $alias);
1145 :     }
1146 :     last;
1147 :     }
1148 :     }
1149 :     return @corresponding_ids;
1150 :     }
1151 :     }
1152 :    
1153 :     sub feature_annotations {
1154 :     my($self,$fid,$rawtime) = @_;
1155 :    
1156 :     my $newG = $self->{_genome};
1157 :     my $newGdir = $self->{_orgdir};
1158 :    
1159 :     my @annotations;
1160 :     if (($fid =~ /^fig\|(\d+\.\d+)/) && ($1 eq $newG))
1161 :     {
1162 :     &load_ann($self);
1163 :     if (my $x = $self->{_ann}->{$fid})
1164 :     {
1165 :     @annotations = @{$x};
1166 :     }
1167 :     else
1168 :     {
1169 :     @annotations = ();
1170 :     }
1171 :    
1172 :     if ($rawtime)
1173 :     {
1174 :     return @annotations;
1175 :     }
1176 :     else
1177 :     {
1178 :     return map { my $r = [@$_]; $r->[1] = localtime($r->[1]); $r } @annotations;
1179 :     }
1180 :     }
1181 :     return ();
1182 :     }
1183 :    
1184 :     sub get_translation {
1185 :     my($self,$peg) = @_;
1186 :    
1187 :     my $newG = $self->{_genome};
1188 :     my $newGdir = $self->{_orgdir};
1189 :    
1190 :     if (($peg =~ /^fig\|(\d+\.\d+)/) && ($1 eq $newG))
1191 :     {
1192 :     &load_feature_indexes($self);
1193 :    
1194 :     my $out = $self->{_feat}->{peg}->{$peg};
1195 :     if (defined($out))
1196 :     {
1197 :     return $out;
1198 :     }
1199 :    
1200 :     #
1201 :     # If we have a blast-formatted fasta, use fastacmd to
1202 :     # do the lookup, and cache the output for later use.
1203 :     #
1204 :    
1205 :     if ($self->{_feat_fasta}->{peg})
1206 :     {
1207 :     my $id = "gnl|$peg";
1208 :     my $cmd = "$FIG_Config::ext_bin/fastacmd -d $self->{_feat_fasta}->{peg} -s '$id'";
1209 :     open(P, "$cmd|") or die "get_translation: cmd failed with $!: $cmd";
1210 :     $_ = <P>;
1211 :     my $out;
1212 :     while (<P>)
1213 :     {
1214 :     s/\s+//g;
1215 :     $out .= $_;
1216 :     }
1217 :     close(P);
1218 :     $self->{_feat}->{$peg} = $out;
1219 :     return $out;
1220 :     }
1221 :     else
1222 :     {
1223 :     return $self->{_feat}->{$peg};
1224 :     }
1225 :     }
1226 :     else
1227 :     {
1228 :     return undef;
1229 :     }
1230 :     }
1231 :    
1232 :     sub translation_length
1233 :     {
1234 :     my($self, $peg) = @_;
1235 :    
1236 :     my $newG = $self->{_genome};
1237 :     my $newGdir = $self->{_orgdir};
1238 :    
1239 :     if (($peg =~ /^fig\|(\d+\.\d+)/) && ($1 eq $newG))
1240 :     {
1241 :     my $t = $self->get_translation($peg);
1242 :     return length($t);
1243 :     }
1244 :     return undef;
1245 :     }
1246 :    
1247 :     sub load_feature_hash {
1248 :     my($self,$type) = @_;
1249 :    
1250 :     my $newGdir = $self->{_orgdir};
1251 :     my $typeH = {};
1252 :     if (open(FIDS,"<$newGdir/Features/$type/tbl"))
1253 :     {
1254 :     while ($_ = <FIDS>)
1255 :     {
1256 :     if ($_ =~ /^(\S+)/)
1257 :     {
1258 :     my $fid = $1;
1259 :     if (! $self->is_deleted_fid($fid))
1260 :     {
1261 :     $typeH->{$fid} = 1;
1262 :     }
1263 :     }
1264 :     }
1265 :     close(FIDS);
1266 :     }
1267 :     return $typeH;
1268 :     }
1269 :    
1270 :     sub load_feature_indexes
1271 :     {
1272 :     my($self) = @_;
1273 :    
1274 :     if ($self->{_feat}) { return };
1275 :    
1276 :     my $newGdir = $self->{_orgdir};
1277 :    
1278 :     for my $fdir (<$newGdir/Features/*>)
1279 :     {
1280 :     my $ftype = basename($fdir);
1281 :    
1282 :     #
1283 :     # If we have a tbl.btree, tie that for our use.
1284 :     #
1285 :    
1286 :     my $tbl_idx = {};
1287 :     my $tie = tie %$tbl_idx, 'DB_File', "$fdir/tbl.btree", O_RDONLY, 0666, $DB_BTREE;
1288 :     if ($tie)
1289 :     {
1290 :     $self->{_feat_tie}->{$ftype} = $tie;
1291 :     $self->{_feat_btree}->{$ftype} = $tbl_idx;
1292 :     }
1293 :    
1294 :     my $tbl_list = [];
1295 :     my $ltie = tie @$tbl_list, 'DB_File', "$fdir/tbl.recno", O_RDONLY, 0666, $DB_RECNO;
1296 :     if ($tie)
1297 :     {
1298 :     $self->{_feat_ltie}->{$ftype} = $ltie;
1299 :     $self->{_feat_recno}->{$ftype} = $tbl_list;
1300 :     }
1301 :    
1302 :     #
1303 :     # If we have fasta.norm.phr, set _pseq_fasta to the fasta file to use with fastacmd.
1304 :     #
1305 :    
1306 :     my $pseq = {};
1307 :    
1308 :     if (-f "$fdir/fasta.norm.phr")
1309 :     {
1310 :     $self->{_feat_fasta}->{$ftype} = "$fdir/fasta.norm";
1311 :    
1312 :     }
1313 :     else
1314 :     {
1315 :     #
1316 :     # Otherwise, we need to load the data.
1317 :     #
1318 :    
1319 :     if (open(FASTA,"<$newGdir/Features/peg/fasta"))
1320 :     {
1321 :     local $/ = "\n>";
1322 :     my $x;
1323 :     while (defined($x = <FASTA>))
1324 :     {
1325 :     chomp $x;
1326 :     if ($x =~ />?(\S+)[^\n]*\n(.*)/s)
1327 :     {
1328 :     my $peg = $1;
1329 :     my $seq = $2;
1330 :     $seq =~ s/\s//gs;
1331 :     if (! $self->is_deleted_fid($peg))
1332 :     {
1333 :     $pseq->{$peg} = $seq;
1334 :     }
1335 :     }
1336 :     }
1337 :     close(FASTA);
1338 :     }
1339 :     }
1340 :     $self->{_feat}->{$ftype} = $pseq;
1341 :     }
1342 :     }
1343 :    
1344 :     sub load_pseq {
1345 :     my($self) = @_;
1346 :    
1347 :     if ($self->{_pseq}) { return };
1348 :    
1349 :     my $newGdir = $self->{_orgdir};
1350 :     my $fdir = "$newGdir/Features/peg";
1351 :    
1352 :     #
1353 :     # If we have a tbl.btree, tie that for our use.
1354 :     #
1355 :    
1356 :     my $tbl_idx = {};
1357 :     my $tie = tie %$tbl_idx, 'DB_File', "$fdir/tbl.btree", O_RDONLY, 0666, $DB_BTREE;
1358 :     if ($tie)
1359 :     {
1360 :     $self->{_tbl_tie} = $tie;
1361 :     $self->{_tbl_btree} = $tbl_idx;
1362 :     }
1363 :    
1364 :     #
1365 :     # If we have fasta.norm.phr, set _pseq_fasta to the fasta file to use with fastacmd.
1366 :     #
1367 :    
1368 :     my $pseq = {};
1369 :    
1370 :     if (-f "$fdir/fasta.norm.phr")
1371 :     {
1372 :     $self->{_pseq_fasta} = "$fdir/fasta.norm";
1373 :     }
1374 :     else
1375 :     {
1376 :     #
1377 :     # Otherwise, we need to load the data.
1378 :     #
1379 :    
1380 :     if (open(FASTA,"<$newGdir/Features/peg/fasta"))
1381 :     {
1382 :     local $/ = "\n>";
1383 :     my $x;
1384 :     while (defined($x = <FASTA>))
1385 :     {
1386 :     chomp $x;
1387 :     if ($x =~ />?(\S+)[^\n]*\n(.*)/s)
1388 :     {
1389 :     my $peg = $1;
1390 :     my $seq = $2;
1391 :     $seq =~ s/\s//gs;
1392 :     if (! $self->is_deleted_fid($peg))
1393 :     {
1394 :     $pseq->{$peg} = $seq;
1395 :     }
1396 :     }
1397 :     }
1398 :     close(FASTA);
1399 :     }
1400 :     }
1401 :     $self->{_pseq} = $pseq;
1402 :     }
1403 :    
1404 :     sub load_ss_data
1405 :     {
1406 :     my($self, $force) = @_;
1407 :    
1408 :     return if defined($self->{_ss_bindings}) && (! $force);
1409 :    
1410 :     my $newG = $self->{_genome};
1411 :     my $newGdir = $self->{_orgdir};
1412 :    
1413 :     open(S, "<$newGdir/Subsystems/bindings") or die "Cannot open $newGdir/Subsystems/bindings: $!";
1414 :    
1415 :     my $peg_index;
1416 :     my $bindings;
1417 :    
1418 :     while (<S>)
1419 :     {
1420 :     chomp;
1421 :     my($sname, $role, $peg) = split(/\t/);
1422 :     if (! $self->is_deleted_fid($peg))
1423 :     {
1424 :     push(@{$bindings->{$sname}->{$role}}, $peg);
1425 :     push(@{$peg_index->{$peg}}, [$sname, $role]);
1426 :     }
1427 :    
1428 :     }
1429 :     close(S);
1430 :    
1431 :     open(S, "<$newGdir/Subsystems/subsystems") or die "Cannot open $newGdir/Subsystems/subsystems: $!";
1432 :     my $variant;
1433 :     while (<S>)
1434 :     {
1435 :     chomp;
1436 :     my($sname, $v) = split(/\t/);
1437 :     $variant->{$sname} = $v;
1438 :     }
1439 :     close(S);
1440 :    
1441 :     $self->{_ss_bindings} = $bindings;
1442 :     $self->{_ss_variants} = $variant;
1443 :     $self->{_ss_peg_index} = $peg_index;
1444 :     }
1445 :    
1446 :     sub load_tbl {
1447 :     my($self) = @_;
1448 :    
1449 :     if ($self->{_tbl}) { return };
1450 :    
1451 :     my $newGdir = $self->{_orgdir};
1452 :     my $tbl = {};
1453 :    
1454 : olson 1.4 for my $tbl_file (<$newGdir/Features/*/tbl>)
1455 : overbeek 1.1 {
1456 : olson 1.4 if (open(my $fh, "<", $tbl_file))
1457 : overbeek 1.1 {
1458 : olson 1.4 while (defined(my $x = <$fh>))
1459 : overbeek 1.1 {
1460 : olson 1.4 chomp $x;
1461 :     if ($x =~ /^(\S+)\t(\S+)(\t(\S.*\S))?/)
1462 :     {
1463 :     my $fid = $1;
1464 :     my $loc = [split(/,/,$2)];
1465 :     my $aliases = $4 ? [split(/\t/,$4)] : [];
1466 :     if (! $self->is_deleted_fid($fid))
1467 :     {
1468 :     $tbl->{$fid} = [$loc,$aliases];
1469 :     }
1470 :     }
1471 :     else {
1472 :     warn "Bad feature line in $newGdir:$x:\n";
1473 :     }
1474 : overbeek 1.1 }
1475 : olson 1.4 close($fh);
1476 : overbeek 1.1 }
1477 : olson 1.4 else
1478 :     {
1479 :     warn "Cannot open $tbl_file: $!";
1480 : overbeek 1.1 }
1481 :     }
1482 :     print STDERR ("Loaded ", (scalar keys %$tbl), " features from $newGdir\n") if $ENV{FIG_VERBOSE};
1483 :     $self->{_tbl} = $tbl;
1484 :     }
1485 :    
1486 :     sub load_functions {
1487 :     my($self, $force) = @_;
1488 :    
1489 :     if (!$force && $self->{_functions}) { return };
1490 :    
1491 :     my $newG = $self->{_genome};
1492 :     my $newGdir = $self->{_orgdir};
1493 :     my $functions = {};
1494 :     my $roles = {};
1495 :    
1496 :     # order of "cat" is important - proposed_user_functions must be last
1497 : olson 1.4 for my $fn_file (<$newGdir/*functions>)
1498 : overbeek 1.1 {
1499 : olson 1.4 if (open(my $fh, "<", $fn_file))
1500 : overbeek 1.1 {
1501 : olson 1.4 while (defined(my $x = <$fh>))
1502 : overbeek 1.1 {
1503 : olson 1.4 if (($x =~ /^(fig\|(\d+\.\d+)\.\S+)\t(\S[^\t]*\S)/) && ($2 eq $newG))
1504 :     {
1505 :     my $peg = $1;
1506 :     my $f = $3;
1507 :     if (! $self->is_deleted_fid($peg))
1508 : overbeek 1.1 {
1509 : olson 1.4 # user has overridden a function, so delete old roles
1510 :     if (exists($functions->{$peg})) {
1511 :     my @roles = &SeedUtils::roles_of_function($functions->{$peg});
1512 :     foreach $_ (@roles)
1513 :     {
1514 :     delete $roles->{$_};
1515 :     }
1516 :     }
1517 :    
1518 :     # add new roles
1519 :     my @roles = &SeedUtils::roles_of_function($f);
1520 :     foreach $_ (@roles)
1521 :     {
1522 :     push(@{$roles->{$_}},$peg);
1523 :     }
1524 :    
1525 :     # set function
1526 :     $functions->{$peg} = $f;
1527 : overbeek 1.1 }
1528 :     }
1529 :     }
1530 : olson 1.4 close($fh);
1531 :     }
1532 :     else
1533 :     {
1534 :     warn "Cannot open $fn_file: $!";
1535 : overbeek 1.1 }
1536 :     }
1537 :     $self->{_functions} = $functions;
1538 :     $self->{_roles} = $roles;
1539 :     }
1540 :    
1541 :     sub seqs_with_role {
1542 :     my($self,$role,$who,$genome) = @_;
1543 :    
1544 :     my $newG = $self->{_genome};
1545 :     if ($genome && $newG && ($genome eq $newG)) {
1546 :     &load_functions($self);
1547 :     my $pegsL = $self->{_roles}->{$role};
1548 :     return ($pegsL ? @{$pegsL} : ());
1549 :     }
1550 :     return ();
1551 :     }
1552 :    
1553 :     sub load_ann {
1554 :     my($self) = @_;
1555 :    
1556 :     if ($self->{_ann}) { return };
1557 :    
1558 :     my $newGdir = $self->{_orgdir};
1559 :     my $ann = {};
1560 :     if (open(ANN,"<$newGdir/annotations"))
1561 :     {
1562 :     $/ = "\n//\n";
1563 :     while (defined(my $x = <ANN>))
1564 :     {
1565 :     chomp $x;
1566 :     if ($x =~ /^(\S+)\n([^\n]+)\n([^\n]+)\n(.*)/s)
1567 :     {
1568 :     push(@{$ann->{$1}},[$1,$2,$3,"$4\n"]);
1569 :     }
1570 :     }
1571 :     $/ = "\n";
1572 :     close(ANN);
1573 :     }
1574 :     $self->{_ann} = $ann;
1575 :     }
1576 :    
1577 :     sub taxonomy_of {
1578 :     my($self) = @_;
1579 :    
1580 :     my $newG = $self->{_genome};
1581 :     my $newGdir = $self->{_orgdir};
1582 :    
1583 :     my $tax;
1584 :     if (open(TAX,"<$newGdir/TAXONOMY") && ($tax = <TAX>))
1585 :     {
1586 :     chop $tax;
1587 :     return $tax;
1588 :     }
1589 :     else
1590 :     {
1591 :     return "unknown";
1592 :     }
1593 :     }
1594 :    
1595 :     sub is_deleted_fid {
1596 :    
1597 :     return 0;
1598 :     }
1599 :    
1600 : olson 1.3 sub load_correspondences
1601 :     {
1602 :     my($self) = @_;
1603 :    
1604 : olson 1.4 #return if exists($self->{correspondence_index});
1605 : olson 1.3
1606 :     my $dir = $self->{_orgdir};
1607 :     my $index = {};
1608 :    
1609 :     for my $cfile (<$dir/CorrToReferenceGenomes/*>)
1610 :     {
1611 :     if ($cfile =~ m,/\d+\.\d+$,)
1612 :     {
1613 :     if (open(my $fh, "<", $cfile))
1614 :     {
1615 :     while (<$fh>)
1616 :     {
1617 :     my $ent = CorrTableEntry->new($_);
1618 :     push(@{$index->{$ent->id1}}, $ent);
1619 :     }
1620 :     close($fh);
1621 :     }
1622 :     else
1623 :     {
1624 :     warn "Could not open $cfile: $!\n";
1625 :     }
1626 :     }
1627 :     }
1628 :     $self->{correspondence_index} = $index;
1629 :     }
1630 :    
1631 :     sub get_correspondences
1632 :     {
1633 :     my($self, $id) = @_;
1634 :    
1635 :     $self->load_correspondences();
1636 :    
1637 :     return $self->{correspondence_index}->{$id};
1638 :     }
1639 :    
1640 :     sub get_best_correspondences
1641 :     {
1642 :     my($self, $id, $cutoff) = @_;
1643 :     my $corrs = $self->get_correspondences($id);
1644 :    
1645 :     my $lg;
1646 :     my $out = [];
1647 :    
1648 :     for my $hit (sort { SeedUtils::genome_of($a->id2) <=> SeedUtils::genome_of($b->id2) or
1649 :     $a->psc <=> $b->psc }
1650 :     grep { (!defined($cutoff) or $_->psc < $cutoff) and
1651 :     ($_->hitinfo eq '<=>' or $_->func1 eq $_->func2) }
1652 :     @$corrs)
1653 :     {
1654 :     if (defined($lg) && $lg ne SeedUtils::genome_of($hit->id2))
1655 :     {
1656 :     push(@$out, $hit);
1657 :     }
1658 :     $lg = SeedUtils::genome_of($hit->id2);
1659 :     }
1660 :     return $out;
1661 :    
1662 :     }
1663 :    
1664 :     sub get_pin
1665 :     {
1666 :     my($self, $peg, $n, $cutoff) = @_;
1667 :    
1668 :     my $hits = $self->get_best_correspondences($peg, $cutoff);
1669 :     # print join("\t", $_->id2, $_->func2) . "\n" for @$hits;
1670 :     my @pegs = map { $_->id2 }
1671 :     sort { $b->bsc <=> $a->bsc or
1672 :     abs($a->len1 - $a->len2) <=> abs($b->len1 - $b->len2) }
1673 :     @$hits;
1674 :     $#pegs = ($n - 1) if @pegs > $n;
1675 :     return \@pegs;
1676 :     }
1677 :    
1678 :     sub get_context
1679 :     {
1680 :     my($self, $peg, $pin, $width) = @_;
1681 :    
1682 :     #
1683 :     # Determine local context.
1684 :     #
1685 :    
1686 :     my $peg_loc = $self->feature_location($peg);
1687 :     my ($contig, $min, $max, $dir) = &SeedUtils::boundaries_of($peg_loc);
1688 :     my $center = ($min + $max) / 2;
1689 :    
1690 :     my ($my_genes, $minV, $maxV) = $self->genes_in_region($contig,
1691 :     int($center - $width / 2),
1692 :     int($center + $width / 2));
1693 :    
1694 :     my $left_extent = $center - $minV;
1695 :     my $right_extent = $maxV - $center;
1696 :    
1697 :     #
1698 :     # Determine other context. Get locations for the pin, then find the
1699 :     # regions.
1700 :     #
1701 :     my $sap = $self->{_sap};
1702 : olson 1.4 if (!defined($sap))
1703 :     {
1704 :     $sap = $self->{_sap} = SAPserver->new($pseed_url);
1705 :     }
1706 :    
1707 :     my $genome_names = $sap->genome_names(-ids => [map { &SeedUtils::genome_of($_) } @$pin]);
1708 :     my $my_genome = $self->genus_species;
1709 :     $my_genome = $self->genome_id if $my_genome eq '';
1710 :    
1711 :     $genome_names->{SeedUtils::genome_of($peg)} = $my_genome;
1712 : olson 1.3
1713 :     my $locs = $sap->fid_locations(-ids => $pin, -boundaries => 1);
1714 :    
1715 :     my @extended_locs;
1716 :     for my $lpeg (keys %$locs)
1717 :     {
1718 :     my $loc = $locs->{$lpeg};
1719 :     my($lcontig, $lmin, $lmax, $ldir) = &SeedUtils::boundaries_of($loc);
1720 :     my $center = ($lmin + $lmax) / 2;
1721 :     my $left_end = int($center - $left_extent);
1722 :     $left_end = 0 if $left_end < 0;
1723 :     my $right_end = int($center + $right_extent);
1724 :     push(@extended_locs, &SeedUtils::location_string($lcontig, $left_end, $right_end));
1725 :     }
1726 :    
1727 :     my $regions = $sap->genes_in_region(-locations => \@extended_locs, -includeLocation => 1);
1728 :    
1729 :     #
1730 :     # Determine functions.
1731 :     #
1732 :     my @ext_pegs = map { keys %$_ } values %$regions;
1733 :     my $funcs = $sap->ids_to_functions(-ids => \@ext_pegs);
1734 :    
1735 :     #
1736 :     # Overlay with local functions.
1737 :    
1738 :     $funcs->{$_} = $self->function_of($_) for @$my_genes;
1739 :    
1740 :     #
1741 :     # We have the whole shebang now. We can assemble and return.
1742 :     #
1743 :     # Each genome is a list of genes.
1744 :     # Each gene is a list [fid, contig, beg, end, dir, func].
1745 :     #
1746 :    
1747 :     my @context;
1748 :     my $cref = {};
1749 :     my $row = 0;
1750 :    
1751 :     #
1752 :     # Start with local context.
1753 :     #
1754 :    
1755 :     my @loc = map { my($contig, $min, $max, $dir) = &SeedUtils::boundaries_of($self->feature_location($_));
1756 :     [$_, $contig, $min, $max, $dir, $funcs->{$_}, $row] } @$my_genes;
1757 :     push @context, [ sort { $a->[2] <=> $b->[2] } @loc ];
1758 :     $cref->{$_->[0]} = $_ for @loc;
1759 :    
1760 :     $row++;
1761 :     #
1762 :     # And include the pinned region.
1763 :     #
1764 :     for my $loc (@extended_locs)
1765 :     {
1766 :     my $region = $regions->{$loc};
1767 :     my @row;
1768 :     for my $peg (keys %$region)
1769 :     {
1770 :     my($contig, $min, $max, $dir) = &SeedUtils::boundaries_of($region->{$peg});
1771 :     my $ent = [$peg, $contig, $min, $max, $dir, $funcs->{$peg}, $row];
1772 :     $cref->{$peg} = $ent;
1773 :     push(@row, $ent);
1774 :     }
1775 :     $row++;
1776 :    
1777 :     push @context, [ sort { $a->[2] <=> $b->[2] } @row ];
1778 :     }
1779 :    
1780 : olson 1.4 return \@context, $cref, $genome_names;
1781 : olson 1.3 }
1782 :    
1783 :    
1784 : overbeek 1.1 1;

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3