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

Annotation of /FigKernelPackages/SeedV.pm

Parent Directory Parent Directory | Revision Log Revision Log


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

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

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3