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

Annotation of /FigKernelPackages/SeedV.pm

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1.11 - (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 : olson 1.10 my($ctg, $min, $max, $strand) = &SeedUtils::boundaries_of($tuple->[0]);
852 :    
853 :     if (not defined($regmin) or not defined($regmax) or
854 :     ($min < $regmin and $regmin < $max) or ($min < $regmax and $regmax < $max) or ($min > $regmin and $max < $regmax))
855 : overbeek 1.1 {
856 : olson 1.10 my $function = $self->function_of($fid);
857 :     push(@$feat_details,[$fid,join(",", @{$tuple->[0]}),join(",",@{$tuple->[1]}),$type,$min,$max,$function,'master','']);
858 : overbeek 1.1 }
859 :     }
860 :     }
861 :     }
862 :     return $feat_details;
863 :     }
864 :    
865 :     sub feature_location {
866 :     my($self,$fid) = @_;
867 :    
868 :     my $newG = $self->{_genome};
869 :     my $newGdir = $self->{_orgdir};
870 :    
871 :     if (($fid =~ /^fig\|(\d+\.\d+)\.([^.]+)/) && ($1 eq $newG))
872 :     {
873 :     my $ftype = $2;
874 :    
875 :     $self->load_feature_indexes();
876 :    
877 : olson 1.5 my $btree = exists($self->{_feat_btree}) ? $self->{_feat_btree}->{$ftype} : undef;
878 : overbeek 1.1 if ($btree)
879 :     {
880 :     my($loc, $idx, @aliases) = split(/$;/, $btree->{$fid});
881 :     return wantarray ? split(/,/, $loc) : $loc;
882 :     }
883 :     else
884 :     {
885 :     &load_tbl($self);
886 :     if (my $x = $self->{_tbl}->{$fid})
887 :     {
888 :     if (wantarray)
889 :     {
890 :     return @{$x->[0]};
891 :     }
892 :     else
893 :     {
894 :     return join(",",@{$x->[0]});
895 :     }
896 :     }
897 :     else
898 :     {
899 :     return undef;
900 :     }
901 :     }
902 :     }
903 :     return undef;
904 :     }
905 :    
906 :     sub function_of {
907 :     my($self,$fid) = @_;
908 :    
909 :     my $newG = $self->{_genome};
910 :     my $newGdir = $self->{_orgdir};
911 :    
912 :     if (($fid =~ /^fig\|(\d+\.\d+)/) && ($1 eq $newG))
913 :     {
914 :     &load_functions($self);
915 :    
916 :     my $fn = $self->{_functions}->{$fid};
917 :     if (wantarray)
918 :     {
919 :     return ['master', $fn];
920 :     }
921 :     else
922 :     {
923 :     return $fn;
924 :     }
925 :     }
926 :     return "";
927 :     }
928 :    
929 : olson 1.4 sub assign_function
930 :     {
931 :     my($self, $fid, $user, $function, $confidence) = @_;
932 :    
933 :     $confidence = $confidence ? $confidence : "";
934 :    
935 :     my $newG = $self->{_genome};
936 :     my $newGdir = $self->{_orgdir};
937 :    
938 :     if (($fid =~ /^fig\|(\d+\.\d+)/) && ($1 ne $newG))
939 :     {
940 :     warn "assign_function on non-seedv fid\n";
941 :     return 0;
942 :     }
943 :    
944 :     $function =~ s/\s+/ /sg; # No multiple spaces
945 :     $function =~ s/^\s+//; # No space at begining
946 :     $function =~ s/\s+$//; # No space at end
947 :     $function =~ s/ ; /; /g; # No space before semicolon
948 :    
949 :     my $file = "$newGdir/proposed_user_functions";
950 :     my $status = 1;
951 :     if ( open( TMP, ">>$file" ) )
952 :     {
953 :     print TMP "$fid\t$function\t$confidence\n";
954 :     close(TMP);
955 :     }
956 :     else
957 :     {
958 :     print STDERR "FAILED ASSIGNMENT: $fid\t$function\t$confidence\n";
959 :     $status = 0;
960 :     }
961 :    
962 :     # mdj: force reload of functions to pick up new assignment
963 :     $self->load_functions(1);
964 :    
965 :     # We are not getting annotations logged. So, we will impose it here.
966 :     $self->add_annotation( $fid, $user, "Set master function to\n$function\n" );
967 :    
968 :     #
969 :     # Mark the genome directory as in need of having bindings recomputed.
970 :     #
971 :     # if (open(S, "<$newGdir/Subsystems/subsystems"))
972 :     # {
973 :     # while (<S>)
974 :     # {
975 :     # chomp;
976 :     # my($sname, $v) = split(/\t/);
977 :     # open(SFILE, ">$self->{_orgdir}/Subsystems/${sname}_bindings_need_recomputation");
978 :     # close(SFILE);
979 :     # }
980 :     # close(S);
981 :     # }
982 :     return $status;
983 :     }
984 :    
985 :     sub add_annotation {
986 :     my($self,$feature_id,$user,$annotation, $time_made) = @_;
987 :    
988 :     my $newG = $self->{_genome};
989 :     my $newGdir = $self->{_orgdir};
990 :    
991 :     if (($feature_id =~ /^fig\|(\d+\.\d+)/) && ($1 ne $newG))
992 :     {
993 : olson 1.6 warn "add_annotation on non-seedv fid '$feature_id'\n";
994 : olson 1.4 return 0;
995 :     }
996 :    
997 :     $time_made = time unless $time_made =~ /^\d+$/;
998 :    
999 :     if ($self->is_deleted_fid($feature_id)) { return 0 }
1000 :    
1001 :     # print STDERR "add: fid=$feature_id user=$user annotation=$annotation\n";
1002 :    
1003 :     my $file = "$newGdir/annotations";
1004 :     my $ma = ($annotation =~ /^Set master function to/);
1005 :    
1006 :     if (open(TMP,">>$file"))
1007 :     {
1008 :     my $dataLine = "$feature_id\n$time_made\n$user\n$annotation" . ((substr($annotation,-1) eq "\n") ? "" : "\n");
1009 :     print TMP $dataLine . "//\n";
1010 :     close(TMP);
1011 :    
1012 :     #
1013 :     # Update local cache.
1014 :     #
1015 :     my $ann = $self->{_ann};
1016 :     push(@{$ann->{$feature_id}}, [$feature_id, $time_made, $user, $annotation . "\n"]);
1017 :     }
1018 : olson 1.6 else
1019 :     {
1020 :     warn "Cannot write $file: $!";
1021 :     }
1022 : olson 1.4 return 0;
1023 :     }
1024 :    
1025 : overbeek 1.1
1026 :     =pod
1027 :    
1028 :     find_features_by_annotation
1029 :    
1030 :     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.
1031 :    
1032 :     If the case boolean is set the search is case insensitive. Otherwise case sensitive.
1033 :    
1034 :     =cut
1035 :    
1036 :     sub find_features_by_annotation {
1037 :     my($self,$anno_hash, $case)=@_;
1038 :     $self->load_functions;
1039 :    
1040 :     if ($case) {map {$anno_hash->{uc($_)}=1} keys %$anno_hash}
1041 :    
1042 :     my $res={};
1043 :     foreach my $id (keys %{$self->{_functions}})
1044 :     {
1045 :     my $fn = $self->{_functions}->{$id};
1046 :     $case ? $fn = uc($fn) : 1;
1047 :     if ($anno_hash->{$fn}) {push @{$res->{$fn}}, $id}
1048 :     }
1049 :    
1050 :     return $res;
1051 :     }
1052 :    
1053 :    
1054 :     =pod
1055 :    
1056 :     search_features_by_annotation
1057 :    
1058 : olson 1.3 Takes a string to find and an optional case boolean, and returns a
1059 :     hash with keys are the function and values are a reference to an array
1060 :     of the IDs that have that function.
1061 : overbeek 1.1
1062 :     If the case boolean is set the search is case insensitive. Otherwise case sensitive.
1063 :    
1064 : olson 1.3 Note that this was originally based on the find above, but this uses a
1065 :     regexp for matching. Will likely be a lot slower which is why I only
1066 :     search for a single term. There may be an SQL way of doing this, if so
1067 :     let Rob know how to do it and I'll replace this method.
1068 : overbeek 1.1
1069 :    
1070 :     =cut
1071 :    
1072 :     sub search_features_by_annotation {
1073 :     my($self,$term, $case)=@_;
1074 :     $self->load_functions;
1075 :    
1076 :     # to make case insensitive convert everything to uppercase
1077 :     # alternative is to use two regexps, one for case insens and one for not case insens
1078 :     # but the bad thing about that approach is that if you have a case insensitive search
1079 :     # you do two regexps for each failed match
1080 :    
1081 :     $case ? $term = uc($term) : 1;
1082 :    
1083 :     my $res={};
1084 :     foreach my $id (keys %{$self->{_functions}})
1085 :     {
1086 :     # we set two variables, one that has case changed for case insensitive searches
1087 :     my $fn = my $fnc = $self->{_functions}->{$id};
1088 :     $case ? $fn = uc($fn) : 1;
1089 :     if ($fn =~ m/$term/) {push @{$res->{$fnc}}, $id}
1090 :     }
1091 :    
1092 :     return $res;
1093 :     }
1094 :    
1095 :    
1096 :     sub feature_aliases {
1097 :     my($self,$fid) = @_;
1098 :    
1099 :     my $newG = $self->{_genome};
1100 :     my $newGdir = $self->{_orgdir};
1101 :    
1102 :     my @aliases;
1103 :    
1104 :     if (($fid =~ /^fig\|(\d+\.\d+)/) && ($1 eq $newG))
1105 :     {
1106 :     &load_tbl($self);
1107 :     if (my $x = $self->{_tbl}->{$fid})
1108 :     {
1109 :     @aliases = @{$x->[1]};
1110 :     }
1111 :     else
1112 :     {
1113 :     @aliases = ();
1114 :     }
1115 :     }
1116 :     else
1117 :     {
1118 :     @aliases = ();
1119 :     }
1120 :     return wantarray() ? @aliases : join(",",@aliases);
1121 :     }
1122 :    
1123 :     sub get_corresponding_ids {
1124 :     my($self, $id, $with_type_info) = @_;
1125 :    
1126 :     my $newG = $self->{_genome};
1127 :    
1128 :     if (($id =~ /^fig\|(\d+\.\d+)/) && ($1 eq $newG)) {
1129 :     my @aliases = $self->feature_aliases($id);
1130 :     my @corresponding_ids = ();
1131 :     foreach my $alias (@aliases) {
1132 :     if ($alias =~ /^gi\|/) {
1133 :     if ($with_type_info) {
1134 :     push(@corresponding_ids, [$alias, 'NCBI']);
1135 :     } else {
1136 :     push(@corresponding_ids, $alias);
1137 :     }
1138 :     last;
1139 :     }
1140 :     }
1141 :     return @corresponding_ids;
1142 :     }
1143 :     }
1144 :    
1145 :     sub feature_annotations {
1146 :     my($self,$fid,$rawtime) = @_;
1147 :    
1148 :     my $newG = $self->{_genome};
1149 :     my $newGdir = $self->{_orgdir};
1150 :    
1151 :     my @annotations;
1152 :     if (($fid =~ /^fig\|(\d+\.\d+)/) && ($1 eq $newG))
1153 :     {
1154 :     &load_ann($self);
1155 :     if (my $x = $self->{_ann}->{$fid})
1156 :     {
1157 :     @annotations = @{$x};
1158 :     }
1159 :     else
1160 :     {
1161 :     @annotations = ();
1162 :     }
1163 :    
1164 :     if ($rawtime)
1165 :     {
1166 :     return @annotations;
1167 :     }
1168 :     else
1169 :     {
1170 :     return map { my $r = [@$_]; $r->[1] = localtime($r->[1]); $r } @annotations;
1171 :     }
1172 :     }
1173 :     return ();
1174 :     }
1175 :    
1176 :     sub get_translation {
1177 :     my($self,$peg) = @_;
1178 :    
1179 :     my $newG = $self->{_genome};
1180 :     my $newGdir = $self->{_orgdir};
1181 :    
1182 :     if (($peg =~ /^fig\|(\d+\.\d+)/) && ($1 eq $newG))
1183 :     {
1184 :     &load_feature_indexes($self);
1185 :    
1186 :     my $out = $self->{_feat}->{peg}->{$peg};
1187 :     if (defined($out))
1188 :     {
1189 :     return $out;
1190 :     }
1191 :    
1192 :     #
1193 :     # If we have a blast-formatted fasta, use fastacmd to
1194 :     # do the lookup, and cache the output for later use.
1195 :     #
1196 :    
1197 :     if ($self->{_feat_fasta}->{peg})
1198 :     {
1199 :     my $id = "gnl|$peg";
1200 :     my $cmd = "$FIG_Config::ext_bin/fastacmd -d $self->{_feat_fasta}->{peg} -s '$id'";
1201 :     open(P, "$cmd|") or die "get_translation: cmd failed with $!: $cmd";
1202 :     $_ = <P>;
1203 :     my $out;
1204 :     while (<P>)
1205 :     {
1206 :     s/\s+//g;
1207 :     $out .= $_;
1208 :     }
1209 :     close(P);
1210 :     $self->{_feat}->{$peg} = $out;
1211 :     return $out;
1212 :     }
1213 :     else
1214 :     {
1215 :     return $self->{_feat}->{$peg};
1216 :     }
1217 :     }
1218 :     else
1219 :     {
1220 :     return undef;
1221 :     }
1222 :     }
1223 :    
1224 :     sub translation_length
1225 :     {
1226 :     my($self, $peg) = @_;
1227 :    
1228 :     my $newG = $self->{_genome};
1229 :     my $newGdir = $self->{_orgdir};
1230 :    
1231 :     if (($peg =~ /^fig\|(\d+\.\d+)/) && ($1 eq $newG))
1232 :     {
1233 :     my $t = $self->get_translation($peg);
1234 :     return length($t);
1235 :     }
1236 :     return undef;
1237 :     }
1238 :    
1239 :     sub load_feature_hash {
1240 :     my($self,$type) = @_;
1241 :    
1242 :     my $newGdir = $self->{_orgdir};
1243 :     my $typeH = {};
1244 :     if (open(FIDS,"<$newGdir/Features/$type/tbl"))
1245 :     {
1246 :     while ($_ = <FIDS>)
1247 :     {
1248 :     if ($_ =~ /^(\S+)/)
1249 :     {
1250 :     my $fid = $1;
1251 :     if (! $self->is_deleted_fid($fid))
1252 :     {
1253 :     $typeH->{$fid} = 1;
1254 :     }
1255 :     }
1256 :     }
1257 :     close(FIDS);
1258 :     }
1259 :     return $typeH;
1260 :     }
1261 :    
1262 :     sub load_feature_indexes
1263 :     {
1264 :     my($self) = @_;
1265 :    
1266 :     if ($self->{_feat}) { return };
1267 :    
1268 :     my $newGdir = $self->{_orgdir};
1269 :    
1270 : olson 1.8 for my $ftype ($self->feature_types())
1271 : overbeek 1.1 {
1272 : olson 1.8 my $fdir = "$newGdir/Features/$ftype";
1273 : overbeek 1.1 my $ftype = basename($fdir);
1274 :    
1275 :     #
1276 :     # If we have a tbl.btree, tie that for our use.
1277 :     #
1278 :    
1279 :     my $tbl_idx = {};
1280 :     my $tie = tie %$tbl_idx, 'DB_File', "$fdir/tbl.btree", O_RDONLY, 0666, $DB_BTREE;
1281 :     if ($tie)
1282 :     {
1283 :     $self->{_feat_tie}->{$ftype} = $tie;
1284 :     $self->{_feat_btree}->{$ftype} = $tbl_idx;
1285 :     }
1286 :    
1287 :     my $tbl_list = [];
1288 :     my $ltie = tie @$tbl_list, 'DB_File', "$fdir/tbl.recno", O_RDONLY, 0666, $DB_RECNO;
1289 :     if ($tie)
1290 :     {
1291 :     $self->{_feat_ltie}->{$ftype} = $ltie;
1292 :     $self->{_feat_recno}->{$ftype} = $tbl_list;
1293 :     }
1294 :    
1295 :     #
1296 :     # If we have fasta.norm.phr, set _pseq_fasta to the fasta file to use with fastacmd.
1297 :     #
1298 :    
1299 :     my $pseq = {};
1300 :    
1301 :     if (-f "$fdir/fasta.norm.phr")
1302 :     {
1303 :     $self->{_feat_fasta}->{$ftype} = "$fdir/fasta.norm";
1304 :    
1305 :     }
1306 :     else
1307 :     {
1308 :     #
1309 :     # Otherwise, we need to load the data.
1310 :     #
1311 :    
1312 :     if (open(FASTA,"<$newGdir/Features/peg/fasta"))
1313 :     {
1314 :     local $/ = "\n>";
1315 :     my $x;
1316 :     while (defined($x = <FASTA>))
1317 :     {
1318 :     chomp $x;
1319 :     if ($x =~ />?(\S+)[^\n]*\n(.*)/s)
1320 :     {
1321 :     my $peg = $1;
1322 :     my $seq = $2;
1323 :     $seq =~ s/\s//gs;
1324 :     if (! $self->is_deleted_fid($peg))
1325 :     {
1326 :     $pseq->{$peg} = $seq;
1327 :     }
1328 :     }
1329 :     }
1330 :     close(FASTA);
1331 :     }
1332 :     }
1333 :     $self->{_feat}->{$ftype} = $pseq;
1334 :     }
1335 :     }
1336 :    
1337 :     sub load_pseq {
1338 :     my($self) = @_;
1339 :    
1340 :     if ($self->{_pseq}) { return };
1341 :    
1342 :     my $newGdir = $self->{_orgdir};
1343 :     my $fdir = "$newGdir/Features/peg";
1344 :    
1345 :     #
1346 :     # If we have a tbl.btree, tie that for our use.
1347 :     #
1348 :    
1349 :     my $tbl_idx = {};
1350 :     my $tie = tie %$tbl_idx, 'DB_File', "$fdir/tbl.btree", O_RDONLY, 0666, $DB_BTREE;
1351 :     if ($tie)
1352 :     {
1353 :     $self->{_tbl_tie} = $tie;
1354 :     $self->{_tbl_btree} = $tbl_idx;
1355 :     }
1356 :    
1357 :     #
1358 :     # If we have fasta.norm.phr, set _pseq_fasta to the fasta file to use with fastacmd.
1359 :     #
1360 :    
1361 :     my $pseq = {};
1362 :    
1363 :     if (-f "$fdir/fasta.norm.phr")
1364 :     {
1365 :     $self->{_pseq_fasta} = "$fdir/fasta.norm";
1366 :     }
1367 :     else
1368 :     {
1369 :     #
1370 :     # Otherwise, we need to load the data.
1371 :     #
1372 :    
1373 :     if (open(FASTA,"<$newGdir/Features/peg/fasta"))
1374 :     {
1375 :     local $/ = "\n>";
1376 :     my $x;
1377 :     while (defined($x = <FASTA>))
1378 :     {
1379 :     chomp $x;
1380 :     if ($x =~ />?(\S+)[^\n]*\n(.*)/s)
1381 :     {
1382 :     my $peg = $1;
1383 :     my $seq = $2;
1384 :     $seq =~ s/\s//gs;
1385 :     if (! $self->is_deleted_fid($peg))
1386 :     {
1387 :     $pseq->{$peg} = $seq;
1388 :     }
1389 :     }
1390 :     }
1391 :     close(FASTA);
1392 :     }
1393 :     }
1394 :     $self->{_pseq} = $pseq;
1395 :     }
1396 :    
1397 :     sub load_ss_data
1398 :     {
1399 :     my($self, $force) = @_;
1400 :    
1401 :     return if defined($self->{_ss_bindings}) && (! $force);
1402 :    
1403 :     my $newG = $self->{_genome};
1404 :     my $newGdir = $self->{_orgdir};
1405 :    
1406 :     open(S, "<$newGdir/Subsystems/bindings") or die "Cannot open $newGdir/Subsystems/bindings: $!";
1407 :    
1408 :     my $peg_index;
1409 :     my $bindings;
1410 :    
1411 :     while (<S>)
1412 :     {
1413 :     chomp;
1414 :     my($sname, $role, $peg) = split(/\t/);
1415 :     if (! $self->is_deleted_fid($peg))
1416 :     {
1417 :     push(@{$bindings->{$sname}->{$role}}, $peg);
1418 :     push(@{$peg_index->{$peg}}, [$sname, $role]);
1419 :     }
1420 :    
1421 :     }
1422 :     close(S);
1423 :    
1424 :     open(S, "<$newGdir/Subsystems/subsystems") or die "Cannot open $newGdir/Subsystems/subsystems: $!";
1425 :     my $variant;
1426 :     while (<S>)
1427 :     {
1428 :     chomp;
1429 :     my($sname, $v) = split(/\t/);
1430 :     $variant->{$sname} = $v;
1431 :     }
1432 :     close(S);
1433 :    
1434 :     $self->{_ss_bindings} = $bindings;
1435 :     $self->{_ss_variants} = $variant;
1436 :     $self->{_ss_peg_index} = $peg_index;
1437 :     }
1438 :    
1439 :     sub load_tbl {
1440 :     my($self) = @_;
1441 :    
1442 :     if ($self->{_tbl}) { return };
1443 :    
1444 :     my $newGdir = $self->{_orgdir};
1445 :     my $tbl = {};
1446 :    
1447 : olson 1.6 my $contig_entries = [];
1448 :     my $cur_contig;
1449 :    
1450 : olson 1.8 for my $type ($self->feature_types())
1451 : overbeek 1.1 {
1452 : olson 1.8 my $tbl_file = "$newGdir/Features/$type/tbl";
1453 : olson 1.6 my($type) = $tbl_file =~ m,/([^/]+)/tbl$,;
1454 : gdpusch 1.9 print STDERR "Load $type\n" if $ENV{VERBOSE};
1455 : olson 1.4 if (open(my $fh, "<", $tbl_file))
1456 : overbeek 1.1 {
1457 : olson 1.4 while (defined(my $x = <$fh>))
1458 : overbeek 1.1 {
1459 : olson 1.4 chomp $x;
1460 :     if ($x =~ /^(\S+)\t(\S+)(\t(\S.*\S))?/)
1461 :     {
1462 :     my $fid = $1;
1463 :     my $loc = [split(/,/,$2)];
1464 : olson 1.6 my($contig) = $loc->[0] =~ /^(.*)_\d+_\d+$/;
1465 : olson 1.4 my $aliases = $4 ? [split(/\t/,$4)] : [];
1466 :     if (! $self->is_deleted_fid($fid))
1467 :     {
1468 :     $tbl->{$fid} = [$loc,$aliases];
1469 : olson 1.6 if ($type eq 'peg' && $contig ne $cur_contig)
1470 :     {
1471 :     push(@$contig_entries, [$contig, $fid]);
1472 :     $cur_contig = $contig;
1473 :     }
1474 : olson 1.4 }
1475 :     }
1476 :     else {
1477 :     warn "Bad feature line in $newGdir:$x:\n";
1478 :     }
1479 : overbeek 1.1 }
1480 : olson 1.4 close($fh);
1481 : overbeek 1.1 }
1482 : olson 1.4 else
1483 :     {
1484 :     warn "Cannot open $tbl_file: $!";
1485 : overbeek 1.1 }
1486 :     }
1487 :     print STDERR ("Loaded ", (scalar keys %$tbl), " features from $newGdir\n") if $ENV{FIG_VERBOSE};
1488 : olson 1.6
1489 :     $self->{_contig_entrypoints} = $contig_entries;
1490 :    
1491 : overbeek 1.1 $self->{_tbl} = $tbl;
1492 :     }
1493 :    
1494 :     sub load_functions {
1495 :     my($self, $force) = @_;
1496 :    
1497 :     if (!$force && $self->{_functions}) { return };
1498 :    
1499 :     my $newG = $self->{_genome};
1500 :     my $newGdir = $self->{_orgdir};
1501 :     my $functions = {};
1502 :     my $roles = {};
1503 :    
1504 :     # order of "cat" is important - proposed_user_functions must be last
1505 : olson 1.8
1506 :     opendir(my $dh, $newGdir);
1507 :     my @fn_files = map { "$newGdir/$_" } sort grep { /functions$/ } readdir($dh);
1508 :     closedir($dh);
1509 :     for my $fn_file (@fn_files)
1510 :     # for my $fn_file (<$newGdir/*functions>)
1511 : overbeek 1.1 {
1512 : olson 1.4 if (open(my $fh, "<", $fn_file))
1513 : overbeek 1.1 {
1514 : olson 1.4 while (defined(my $x = <$fh>))
1515 : overbeek 1.1 {
1516 : olson 1.4 if (($x =~ /^(fig\|(\d+\.\d+)\.\S+)\t(\S[^\t]*\S)/) && ($2 eq $newG))
1517 :     {
1518 :     my $peg = $1;
1519 :     my $f = $3;
1520 :     if (! $self->is_deleted_fid($peg))
1521 : overbeek 1.1 {
1522 : olson 1.4 # user has overridden a function, so delete old roles
1523 :     if (exists($functions->{$peg})) {
1524 :     my @roles = &SeedUtils::roles_of_function($functions->{$peg});
1525 :     foreach $_ (@roles)
1526 :     {
1527 :     delete $roles->{$_};
1528 :     }
1529 :     }
1530 :    
1531 :     # add new roles
1532 :     my @roles = &SeedUtils::roles_of_function($f);
1533 :     foreach $_ (@roles)
1534 :     {
1535 :     push(@{$roles->{$_}},$peg);
1536 :     }
1537 :    
1538 :     # set function
1539 :     $functions->{$peg} = $f;
1540 : overbeek 1.1 }
1541 :     }
1542 :     }
1543 : olson 1.4 close($fh);
1544 :     }
1545 :     else
1546 :     {
1547 :     warn "Cannot open $fn_file: $!";
1548 : overbeek 1.1 }
1549 :     }
1550 :     $self->{_functions} = $functions;
1551 :     $self->{_roles} = $roles;
1552 :     }
1553 :    
1554 :     sub seqs_with_role {
1555 :     my($self,$role,$who,$genome) = @_;
1556 :    
1557 :     my $newG = $self->{_genome};
1558 :     if ($genome && $newG && ($genome eq $newG)) {
1559 :     &load_functions($self);
1560 :     my $pegsL = $self->{_roles}->{$role};
1561 :     return ($pegsL ? @{$pegsL} : ());
1562 :     }
1563 :     return ();
1564 :     }
1565 :    
1566 :     sub load_ann {
1567 :     my($self) = @_;
1568 :    
1569 :     if ($self->{_ann}) { return };
1570 :    
1571 :     my $newGdir = $self->{_orgdir};
1572 :     my $ann = {};
1573 :     if (open(ANN,"<$newGdir/annotations"))
1574 :     {
1575 :     $/ = "\n//\n";
1576 :     while (defined(my $x = <ANN>))
1577 :     {
1578 :     chomp $x;
1579 :     if ($x =~ /^(\S+)\n([^\n]+)\n([^\n]+)\n(.*)/s)
1580 :     {
1581 :     push(@{$ann->{$1}},[$1,$2,$3,"$4\n"]);
1582 :     }
1583 :     }
1584 :     $/ = "\n";
1585 :     close(ANN);
1586 :     }
1587 :     $self->{_ann} = $ann;
1588 :     }
1589 :    
1590 :     sub taxonomy_of {
1591 :     my($self) = @_;
1592 :    
1593 :     my $newG = $self->{_genome};
1594 :     my $newGdir = $self->{_orgdir};
1595 :    
1596 :     my $tax;
1597 :     if (open(TAX,"<$newGdir/TAXONOMY") && ($tax = <TAX>))
1598 :     {
1599 :     chop $tax;
1600 :     return $tax;
1601 :     }
1602 :     else
1603 :     {
1604 :     return "unknown";
1605 :     }
1606 :     }
1607 :    
1608 : olson 1.6 sub is_deleted_fid
1609 :     {
1610 :     my($self, $fid) = @_;
1611 :     my $newG = $self->{_genome};
1612 : overbeek 1.1
1613 : olson 1.6 if (($fid =~ /^fig\|(\d+\.\d+)/) && ($1 eq $newG))
1614 :     {
1615 :     my $delH;
1616 :     $delH = $self->{_deleted_features};
1617 :     if (! $delH)
1618 :     {
1619 :     $delH = $self->load_deleted_features_hash();
1620 :     }
1621 :     return $delH->{$fid} ? 1 : 0;
1622 :     }
1623 : overbeek 1.1 return 0;
1624 :     }
1625 :    
1626 : olson 1.3 sub load_correspondences
1627 :     {
1628 :     my($self) = @_;
1629 :    
1630 : olson 1.4 #return if exists($self->{correspondence_index});
1631 : olson 1.3
1632 :     my $dir = $self->{_orgdir};
1633 :     my $index = {};
1634 : olson 1.6 my $rev_index = {};
1635 : olson 1.3
1636 : olson 1.8 opendir(my $dh, "$dir/CorrToReferenceGenomes");
1637 :     my @corr_files = map { "$dir/CorrToReferenceGenomes/$_" } sort grep { /^\d+\.\d+$/ } readdir($dh);
1638 :     closedir($dh);
1639 :     for my $cfile (@corr_files)
1640 :     # for my $cfile (<$dir/CorrToReferenceGenomes/*>)
1641 : olson 1.3 {
1642 :     if ($cfile =~ m,/\d+\.\d+$,)
1643 :     {
1644 :     if (open(my $fh, "<", $cfile))
1645 :     {
1646 :     while (<$fh>)
1647 :     {
1648 :     my $ent = CorrTableEntry->new($_);
1649 :     push(@{$index->{$ent->id1}}, $ent);
1650 : olson 1.6 push(@{$rev_index->{$ent->id2}}, $ent);
1651 : olson 1.3 }
1652 :     close($fh);
1653 :     }
1654 :     else
1655 :     {
1656 :     warn "Could not open $cfile: $!\n";
1657 :     }
1658 :     }
1659 :     }
1660 :     $self->{correspondence_index} = $index;
1661 : olson 1.6 $self->{correspondence_index_rev} = $rev_index;
1662 : olson 1.3 }
1663 :    
1664 :     sub get_correspondences
1665 :     {
1666 :     my($self, $id) = @_;
1667 :    
1668 :     $self->load_correspondences();
1669 :    
1670 :     return $self->{correspondence_index}->{$id};
1671 :     }
1672 :    
1673 : olson 1.6 sub get_correspondences_rev
1674 :     {
1675 :     my($self, $id) = @_;
1676 :    
1677 :     $self->load_correspondences();
1678 :    
1679 :     return $self->{correspondence_index_rev}->{$id};
1680 :     }
1681 :    
1682 : olson 1.3 sub get_best_correspondences
1683 :     {
1684 :     my($self, $id, $cutoff) = @_;
1685 :     my $corrs = $self->get_correspondences($id);
1686 :    
1687 :     my $lg;
1688 :     my $out = [];
1689 :    
1690 :     for my $hit (sort { SeedUtils::genome_of($a->id2) <=> SeedUtils::genome_of($b->id2) or
1691 :     $a->psc <=> $b->psc }
1692 :     grep { (!defined($cutoff) or $_->psc < $cutoff) and
1693 : olson 1.6 ($_->hitinfo eq '<=>' or SeedUtils::strip_func_comment($_->func1) eq SeedUtils::strip_func_comment($_->func2)) }
1694 : olson 1.3 @$corrs)
1695 :     {
1696 :     if (defined($lg) && $lg ne SeedUtils::genome_of($hit->id2))
1697 :     {
1698 :     push(@$out, $hit);
1699 :     }
1700 :     $lg = SeedUtils::genome_of($hit->id2);
1701 :     }
1702 :     return $out;
1703 :    
1704 :     }
1705 :    
1706 :     sub get_pin
1707 :     {
1708 :     my($self, $peg, $n, $cutoff) = @_;
1709 :    
1710 :     my $hits = $self->get_best_correspondences($peg, $cutoff);
1711 :     # print join("\t", $_->id2, $_->func2) . "\n" for @$hits;
1712 :     my @pegs = map { $_->id2 }
1713 :     sort { $b->bsc <=> $a->bsc or
1714 :     abs($a->len1 - $a->len2) <=> abs($b->len1 - $b->len2) }
1715 :     @$hits;
1716 :     $#pegs = ($n - 1) if @pegs > $n;
1717 :     return \@pegs;
1718 :     }
1719 :    
1720 :     sub get_context
1721 :     {
1722 :     my($self, $peg, $pin, $width) = @_;
1723 :    
1724 :     #
1725 :     # Determine local context.
1726 :     #
1727 :    
1728 :     my $peg_loc = $self->feature_location($peg);
1729 :     my ($contig, $min, $max, $dir) = &SeedUtils::boundaries_of($peg_loc);
1730 :     my $center = ($min + $max) / 2;
1731 :    
1732 :     my ($my_genes, $minV, $maxV) = $self->genes_in_region($contig,
1733 :     int($center - $width / 2),
1734 :     int($center + $width / 2));
1735 :    
1736 : olson 1.6 # my $left_extent = $center - $minV;
1737 :     # my $right_extent = $maxV - $center;
1738 :     my $left_extent = int($width / 2);
1739 :     my $right_extent = int($width / 2);
1740 : olson 1.3
1741 :     #
1742 :     # Determine other context. Get locations for the pin, then find the
1743 :     # regions.
1744 :     #
1745 :     my $sap = $self->{_sap};
1746 : olson 1.4 if (!defined($sap))
1747 :     {
1748 :     $sap = $self->{_sap} = SAPserver->new($pseed_url);
1749 :     }
1750 :    
1751 :     my $genome_names = $sap->genome_names(-ids => [map { &SeedUtils::genome_of($_) } @$pin]);
1752 :     my $my_genome = $self->genus_species;
1753 :     $my_genome = $self->genome_id if $my_genome eq '';
1754 :    
1755 :     $genome_names->{SeedUtils::genome_of($peg)} = $my_genome;
1756 : olson 1.3
1757 :     my $locs = $sap->fid_locations(-ids => $pin, -boundaries => 1);
1758 :    
1759 :     my @extended_locs;
1760 :     for my $lpeg (keys %$locs)
1761 :     {
1762 :     my $loc = $locs->{$lpeg};
1763 :     my($lcontig, $lmin, $lmax, $ldir) = &SeedUtils::boundaries_of($loc);
1764 :     my $center = ($lmin + $lmax) / 2;
1765 :     my $left_end = int($center - $left_extent);
1766 :     $left_end = 0 if $left_end < 0;
1767 :     my $right_end = int($center + $right_extent);
1768 :     push(@extended_locs, &SeedUtils::location_string($lcontig, $left_end, $right_end));
1769 :     }
1770 :    
1771 :     my $regions = $sap->genes_in_region(-locations => \@extended_locs, -includeLocation => 1);
1772 :    
1773 :     #
1774 :     # Determine functions.
1775 :     #
1776 :     my @ext_pegs = map { keys %$_ } values %$regions;
1777 :     my $funcs = $sap->ids_to_functions(-ids => \@ext_pegs);
1778 :    
1779 :     #
1780 :     # Overlay with local functions.
1781 :    
1782 :     $funcs->{$_} = $self->function_of($_) for @$my_genes;
1783 :    
1784 :     #
1785 :     # We have the whole shebang now. We can assemble and return.
1786 :     #
1787 :     # Each genome is a list of genes.
1788 :     # Each gene is a list [fid, contig, beg, end, dir, func].
1789 :     #
1790 :    
1791 :     my @context;
1792 :     my $cref = {};
1793 :     my $row = 0;
1794 :    
1795 :     #
1796 :     # Start with local context.
1797 :     #
1798 :    
1799 :     my @loc = map { my($contig, $min, $max, $dir) = &SeedUtils::boundaries_of($self->feature_location($_));
1800 :     [$_, $contig, $min, $max, $dir, $funcs->{$_}, $row] } @$my_genes;
1801 :     push @context, [ sort { $a->[2] <=> $b->[2] } @loc ];
1802 :     $cref->{$_->[0]} = $_ for @loc;
1803 :    
1804 :     $row++;
1805 :     #
1806 :     # And include the pinned region.
1807 :     #
1808 :     for my $loc (@extended_locs)
1809 :     {
1810 :     my $region = $regions->{$loc};
1811 :     my @row;
1812 :     for my $peg (keys %$region)
1813 :     {
1814 : olson 1.6 my $func = defined($funcs->{$peg}) ? $funcs->{$peg} : "";
1815 : olson 1.3 my($contig, $min, $max, $dir) = &SeedUtils::boundaries_of($region->{$peg});
1816 : olson 1.6 my $ent = [$peg, $contig, $min, $max, $dir, $func, $row];
1817 : olson 1.3 $cref->{$peg} = $ent;
1818 :     push(@row, $ent);
1819 :     }
1820 :     $row++;
1821 :    
1822 :     push @context, [ sort { $a->[2] <=> $b->[2] } @row ];
1823 :     }
1824 :    
1825 : olson 1.4 return \@context, $cref, $genome_names;
1826 : olson 1.3 }
1827 :    
1828 : olson 1.6 sub add_feature {
1829 :     my( $self, $user, $genome, $type, $location, $aliases, $sequence) = @_;
1830 :    
1831 :     my $newG = $self->{_genome};
1832 :     my $newGdir = $self->{_orgdir};
1833 :    
1834 :     if ($genome ne $newG) {
1835 :     return undef;
1836 :     }
1837 :    
1838 :     # perform sanity checks
1839 :     unless ($genome && $type && $location && $sequence) {
1840 :     print STDERR "SEED error: add_feature failed due to missing parameter\n";
1841 :     return undef;
1842 :     }
1843 :    
1844 :     if ( $type !~ /^[0-9A-Za-z_]+$/ )
1845 :     {
1846 :     print STDERR "SEED error: add_feature failed due to bad type: $type\n";
1847 :     return undef;
1848 :     }
1849 :    
1850 :     if ( length ( $location ) > 5000 )
1851 :     {
1852 :     print STDERR "SEED error: add_feature failed because location is over 5000 char:\n";
1853 :     print STDERR "$location\n";
1854 :     return undef;
1855 :     }
1856 :    
1857 :     my @loc = split( /,/, $location );
1858 : gdpusch 1.11 my @loc2 = grep { defined($_->[0]) && $_->[1] && $_->[2] }
1859 : olson 1.6 map { [ $_ =~ m/^(.+)_(\d+)_(\d+)$/ ] }
1860 :     @loc;
1861 :    
1862 :     if ( (@loc2 == 0) || ( @loc != @loc2 ) )
1863 :     {
1864 :     print STDERR "SEED error: add_feature failed because location is missing or malformed:\n";
1865 :     print STDERR "$location\n";
1866 :     return undef;
1867 :     }
1868 :    
1869 :     if ( my @bad_names = grep { length( $_->[0] ) > 96 } @loc2 )
1870 :     {
1871 :     print STDERR "SEED error: add_feature failed because location contains a contig name of over 96 char:\n";
1872 :     print STDERR join( ", ", @bad_names ) . "\n";
1873 :     return undef;
1874 :     }
1875 :    
1876 :     # create type directory if it does not exist
1877 :     unless (-d "$newGdir/Features/$type") {
1878 :     &SeedUtils::verify_dir("$newGdir/Features/$type");
1879 :     (-d "$newGdir/Features/$type")
1880 :     || die qq(Feature directory \'$newGdir/Features/$type\' does not exist, and could not be created);
1881 :    
1882 :     open(TMP, ">", "$newGdir/Features/$type/tbl")
1883 :     || die "Could not create empty $newGdir/Features/$type/tbl: $!";
1884 :     close(TMP);
1885 :    
1886 :     open(TMP, ">", "$newGdir/Features/$type/fasta")
1887 :     || die "Could not create empty $newGdir/Features/$type/fasta: $!";
1888 :     close(TMP);
1889 :     }
1890 :    
1891 :     # create an id
1892 :     my $id = "fig|$genome.$type.";
1893 :     my $feature_dir = "$newGdir/Features/$type";
1894 :     my $file = "$feature_dir/tbl";
1895 :     if (-f $file) {
1896 :     unless (open(FILE, "<$file")) {
1897 :     print STDERR "SEED error: could not open tbl file: $@\n";
1898 :     return undef;
1899 :     }
1900 :     my $entry;
1901 :     my $max = 0;
1902 :     while (defined($entry = <FILE>)) {
1903 :     chomp $entry;
1904 :     if ($entry =~ /^fig\|$genome\.$type\.(\d+)/) {
1905 :     my $curr_id = $1;
1906 :     if ($curr_id > $max) {
1907 :     $max = $curr_id;
1908 :     }
1909 :     }
1910 :     else {
1911 :     confess qq(Could not parse $type tbl entry: $entry);
1912 :     }
1913 :     }
1914 :     close FILE;
1915 :     ++$max;
1916 :     $id .= $max;
1917 :     } else {
1918 :     $id .= "1";
1919 :     }
1920 :    
1921 :     # append to tbl file
1922 :     unless (open(FILE, ">>$file")) {
1923 :     print STDERR "SEED error: could not open tbl file: $@\n";
1924 :     return undef;
1925 :     }
1926 :    
1927 :     $aliases =~ s/,/\t/g;
1928 :     print FILE "$id\t$location\t$aliases\n";
1929 :     close FILE;
1930 :     chmod(0777,$file);
1931 :    
1932 :     my $typeH = $self->{_features}->{$type};
1933 :     if ($typeH)
1934 :     {
1935 :     $typeH->{$id} = 1;
1936 :     }
1937 :    
1938 :     my $tbl = $self->{_tbl};
1939 :     if ($tbl)
1940 :     {
1941 :     $tbl->{$id} = [[@loc],[split(/\t/,$aliases)]];
1942 :     }
1943 :    
1944 :     # append to fasta file
1945 :     $sequence =~ s/\s//g;
1946 :    
1947 :     $file = "$feature_dir/fasta";
1948 :     unless (open(FILE, ">>$file")) {
1949 :     print STDERR "SEED error: could not open fasta file: $@\n";
1950 :     return undef;
1951 :     }
1952 :     print FILE ">$id\n$sequence\n";
1953 :     close FILE;
1954 :     chmod(0777,$file);
1955 :    
1956 :     # append to called_by file
1957 :     $file = "$newGdir/called_by";
1958 :     unless (open(FILE, ">>$file")) {
1959 :     print STDERR "SEED error: could not open called_by file: $@\n";
1960 :     return undef;
1961 :     }
1962 :     print FILE "$id\tmanual: $user\n";
1963 :     close FILE;
1964 :     chmod(0777,$file);
1965 :    
1966 :     #
1967 :     # If a btree was created for this, clear the ref to it and delete the files since
1968 :     # they are now invalid.
1969 :     #
1970 :    
1971 :     my $tie = $self->{_feat_tie}->{$type};
1972 :     my $btree = $self->{_feat_btree}->{$type};
1973 :    
1974 :     if ($tie)
1975 :     {
1976 :     untie $tie;
1977 :     delete $self->{$_}->{$type} for qw(_feat_tie _feat_btree _feat_ltie _feat_recno _feat_fasta);
1978 :    
1979 :     unlink("$feature_dir/tbl.btree");
1980 :     unlink("$feature_dir/tbl.recno");
1981 :     }
1982 :    
1983 :     if (-f "$feature_dir/fasta.norm.phr")
1984 :     {
1985 : olson 1.8 # todo fix dangerous
1986 :     #unlink(<$feature_dir/fasta.norm.*>);
1987 : olson 1.6 }
1988 :    
1989 :    
1990 :     # declare success
1991 :     return $id;
1992 :     }
1993 :    
1994 :     sub delete_feature {
1995 :     my($self,$user,$fid) = @_;
1996 :    
1997 :     my $newG = $self->{_genome};
1998 :     my $newGdir = $self->{_orgdir};
1999 :    
2000 :     if (($fid =~ /^fig\|(\d+\.\d+)/) && ($1 ne $newG))
2001 :     {
2002 :     warn "delete_feature on non-seedv fid\n";
2003 :     return 0;
2004 :     }
2005 :    
2006 :     if (open(DEL,">>$newGdir/deleted.fids"))
2007 :     {
2008 :     print DEL "$fid\t$user\n";
2009 :     close(DEL);
2010 :     $self->load_deleted_features_hash();
2011 :     }
2012 :     else
2013 :     {
2014 :     carp "could not open $newGdir/deleted.fids: failed to delete $fid (user=$user)\n";
2015 :     }
2016 :     }
2017 :    
2018 :     sub load_deleted_features_hash {
2019 :     my($self) = @_;
2020 :    
2021 :     my $newGdir = $self->{_orgdir};
2022 :     my $deletedH = {};
2023 :     if (open(DEL,"<$newGdir/deleted.fids"))
2024 :     {
2025 :     local $/ = "\n";
2026 :     while (<DEL>)
2027 :     {
2028 :     if ($_ =~ /^(\S+)/)
2029 :     {
2030 :     $deletedH->{$1} = 1;
2031 :     }
2032 :     }
2033 :     close(DEL);
2034 :     }
2035 :     $self->{_deleted_features} = $deletedH;
2036 :    
2037 :     return $deletedH;
2038 :     }
2039 :    
2040 :     #
2041 :     # List of feature types that exist in this genome.
2042 :     #
2043 :     sub feature_types
2044 :     {
2045 :     my($self) = @_;
2046 :    
2047 :     my $newGdir = $self->{_orgdir};
2048 :    
2049 :     my %sort_prefs = (peg => 2,
2050 :     rna => 1);
2051 :    
2052 :     if (opendir(my $dh, "$newGdir/Features"))
2053 :     {
2054 :     my @types = grep { ! /^\./ && -d "$newGdir/Features/$_" } readdir($dh);
2055 :     closedir($dh);
2056 :     return sort { $sort_prefs{$b} <=> $sort_prefs{$a} or $a cmp $b } @types;
2057 :     }
2058 :     return ();
2059 :     }
2060 : olson 1.3
2061 : olson 1.7 sub write_features_for_comparison
2062 :     {
2063 :     my($self, $fh) = @_;
2064 :    
2065 :     for my $ent (@{$self->all_features_detailed_fast()})
2066 :     {
2067 :     my($fid, $loc, undef, $type, $min, $max, $fn) = @$ent;
2068 :     next unless $type eq 'peg';
2069 :     print $fh join("\t", $fid, $loc, $fn), "\n";
2070 :     }
2071 :     }
2072 :    
2073 :    
2074 :    
2075 : overbeek 1.1 1;

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3