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

Annotation of /FigKernelPackages/SeedV.pm

Parent Directory Parent Directory | Revision Log Revision Log


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

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3