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

Annotation of /FigKernelPackages/FIG.pm

Parent Directory Parent Directory | Revision Log Revision Log


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

1 : efrank 1.1 package FIG;
2 :    
3 : olson 1.111 use strict;
4 :    
5 : overbeek 1.135 use Fcntl qw/:flock/; # import LOCK_* constants
6 :    
7 : olson 1.116 use POSIX;
8 : olson 1.158 use IPC::Open2;
9 : olson 1.116
10 : efrank 1.1 use DBrtns;
11 :     use Sim;
12 :     use Blast;
13 :     use FIG_Config;
14 : overbeek 1.36 use tree_utilities;
15 : olson 1.93 use Subsystem;
16 : olson 1.162 use SeedDas;
17 : olson 1.183 use Construct;
18 : parrello 1.200 use FIGRules;
19 : parrello 1.210 use Tracer;
20 : olson 1.260
21 : olson 1.245 eval { require FigGFF; };
22 : olson 1.260 if ($@ and $ENV{USER} eq "olson")
23 :     {
24 :     warn $@;
25 :     }
26 : olson 1.79
27 :     #
28 :     # Conditionally evaluate this in case its prerequisites are not available.
29 :     #
30 :    
31 :     our $ClearinghouseOK = eval {
32 :     require Clearinghouse;
33 :     };
34 : efrank 1.1
35 : olson 1.10 use IO::Socket;
36 :    
37 : efrank 1.1 use FileHandle;
38 :    
39 :     use Carp;
40 :     use Data::Dumper;
41 : overbeek 1.25 use Time::Local;
42 : olson 1.93 use File::Spec;
43 : olson 1.123 use File::Copy;
44 : olson 1.112 #
45 :     # Try to load the RPC stuff; it might fail on older versions of the software.
46 :     #
47 :     eval {
48 :     require FIGrpc;
49 :     };
50 :    
51 :     my $xmlrpc_available = 1;
52 :     if ($@ ne "")
53 :     {
54 :     $xmlrpc_available = 0;
55 :     }
56 :    
57 : efrank 1.1
58 : olson 1.111 use FIGAttributes;
59 :     use base 'FIGAttributes';
60 :    
61 :     use vars qw(%_FunctionAttributes);
62 :    
63 :     use Data::Dumper;
64 :    
65 : olson 1.124 #
66 :     # Force all new files to be all-writable.
67 :     #
68 :    
69 :     umask 0;
70 :    
71 : parrello 1.210 =head1 FIG Genome Annotation System
72 :    
73 :     =head2 Introduction
74 :    
75 :     This is the main object for access to the SEED data store. The data store
76 :     itself is a combination of flat files and a database. The flat files can
77 :     be moved easily between systems and the database rebuilt as needed.
78 :    
79 :     A reduced set of this object's functions are available via the B<SFXlate>
80 :     object. The SFXlate object uses a single database to represent all its
81 :     genomic information. It provides a much smaller capability for updating
82 :     the data, and eliminates all similarities except for bidirectional best
83 :     hits.
84 :    
85 :     The key to making the FIG system work is proper configuration of the
86 :     C<FIG_Config.pm> file. This file contains names and URLs for the key
87 :     directories as well as the type and login information for the database.
88 :    
89 :     =cut
90 :    
91 :     #: Constructor FIG->new();
92 :    
93 :     =head2 Public Methods
94 :    
95 :     =head3 new
96 :    
97 :     C<< my $fig = FIG->new(); >>
98 :    
99 :     This is the constructor for a FIG object. It uses no parameters.
100 :    
101 :     =cut
102 :    
103 : efrank 1.1 sub new {
104 :     my($class) = @_;
105 :    
106 : olson 1.102 #
107 :     # Check to see if we have a FIG_URL environment variable set.
108 :     # If we do, don't actually create a FIG object, but rather
109 :     # create a FIGrpc and return that as the return from this constructor.
110 :     #
111 :    
112 : parrello 1.210 if ($ENV{FIG_URL} ne "" && $xmlrpc_available) {
113 :     Trace("Creating figrpc for '$ENV{FIG_URL}'") if T(0);
114 :     my $figrpc = new FIGrpc($ENV{FIG_URL});
115 :     return $figrpc;
116 : olson 1.102 }
117 :    
118 : efrank 1.1 my $rdbH = new DBrtns;
119 :     bless {
120 : parrello 1.210 _dbf => $rdbH,
121 :     }, $class;
122 : efrank 1.1 }
123 :    
124 : parrello 1.210
125 :     =head3 get_system_name
126 :    
127 :     C<< my $name = $fig->get_system_name; >>
128 :    
129 :     Returns C<seed>, indicating that this is object is using the SEED
130 :     database. The same method on an SFXlate object will return C<sprout>.
131 :    
132 :     =cut
133 :     #: Return Type $;
134 :     sub get_system_name {
135 : olson 1.207 return "seed";
136 : olson 1.205 }
137 : parrello 1.210
138 :     # Destructor: releases the database handle.
139 : olson 1.205
140 : efrank 1.1 sub DESTROY {
141 :     my($self) = @_;
142 :     my($rdbH);
143 :    
144 : parrello 1.210 if ($rdbH = $self->db_handle) {
145 :     $rdbH->DESTROY;
146 : efrank 1.1 }
147 :     }
148 :    
149 : parrello 1.210 =head3 delete_genomes
150 :    
151 :     C<< $fig->delete_genomes(\@genomes); >>
152 :    
153 :     Delete the specified genomes from the data store. This requires making
154 :     system calls to move and delete files.
155 :    
156 :     =cut
157 :     #: Return Type ;
158 : overbeek 1.7 sub delete_genomes {
159 :     my($self,$genomes) = @_;
160 :     my $tmpD = "$FIG_Config::temp/tmp.deleted.$$";
161 :     my $tmp_Data = "$FIG_Config::temp/Data.$$";
162 :    
163 :     my %to_del = map { $_ => 1 } @$genomes;
164 :     open(TMP,">$tmpD") || die "could not open $tmpD";
165 :    
166 :     my $genome;
167 :     foreach $genome ($self->genomes)
168 :     {
169 :     if (! $to_del{$genome})
170 :     {
171 :     print TMP "$genome\n";
172 :     }
173 :     }
174 :     close(TMP);
175 :    
176 :     &run("extract_genomes $tmpD $FIG_Config::data $tmp_Data");
177 : parrello 1.200
178 : overbeek 1.47 # &run("mv $FIG_Config::data $FIG_Config::data.deleted; mv $tmp_Data $FIG_Config::data; fig load_all; rm -rf $FIG_Config::data.deleted");
179 : parrello 1.200
180 :     &run("mv $FIG_Config::data $FIG_Config::data.deleted");
181 : overbeek 1.47 &run("mv $tmp_Data $FIG_Config::data");
182 :     &run("fig load_all");
183 :     &run("rm -rf $FIG_Config::data.deleted");
184 : overbeek 1.7 }
185 : parrello 1.200
186 : parrello 1.210 =head3 add_genome
187 :    
188 :     C<< my $ok = $fig->add_genome($genomeF); >>
189 :    
190 :     Add a new genome to the data store. A genome's data is kept in a directory
191 :     by itself, underneath the main organism directory.
192 :    
193 :     =over 4
194 :    
195 :     =item genomeF
196 :    
197 :     Name of the directory containing the genome files.
198 :    
199 :     =item RETURN
200 :    
201 :     Returns TRUE if successful, else FALSE.
202 :    
203 :     =back
204 :    
205 :     =cut
206 :     #: Return Type $;
207 : efrank 1.1 sub add_genome {
208 :     my($self,$genomeF) = @_;
209 :    
210 :     my $rc = 0;
211 : olson 1.93
212 :     my(undef, $path, $genome) = File::Spec->splitpath($genomeF);
213 :    
214 :     if ($genome !~ /^\d+\.\d+$/)
215 :     {
216 :     warn "Invalid genome filename $genomeF\n";
217 :     return $rc;
218 :     }
219 :    
220 :     if (-d $FIG_Config::organisms/$genome)
221 :     {
222 :     warn "Organism already exists for $genome\n";
223 :     return $rc;
224 :     }
225 : parrello 1.200
226 : olson 1.93
227 :     #
228 :     # We're okay, it doesn't exist.
229 :     #
230 :    
231 :     my @errors = `$FIG_Config::bin/verify_genome_directory $genomeF`;
232 :    
233 :     if (@errors)
234 : efrank 1.1 {
235 : olson 1.93 warn "Errors found while verifying genome directory $genomeF:\n";
236 :     print join("", @errors);
237 :     return $rc;
238 :     }
239 : parrello 1.200
240 : olson 1.93 &run("cp -r $genomeF $FIG_Config::organisms");
241 :     &run("chmod -R 777 $FIG_Config::organisms/$genome");
242 :    
243 :     &run("index_contigs $genome");
244 :     &run("compute_genome_counts $genome");
245 :     &run("load_features $genome");
246 :    
247 :     $rc = 1;
248 :     if (-s "$FIG_Config::organisms/$genome/Features/peg/fasta")
249 :     {
250 :     &run("index_translations $genome");
251 :     my @tmp = `cut -f1 $FIG_Config::organisms/$genome/Features/peg/tbl`;
252 :     chomp @tmp;
253 :     &run("cat $FIG_Config::organisms/$genome/Features/peg/fasta >> $FIG_Config::data/Global/nr");
254 :     &enqueue_similarities(\@tmp);
255 :     }
256 :     if ((-s "$FIG_Config::organisms/$genome/assigned_functions") ||
257 :     (-d "$FIG_Config::organisms/$genome/UserModels"))
258 :     {
259 :     &run("add_assertions_of_function $genome");
260 : efrank 1.1 }
261 : parrello 1.200
262 : efrank 1.1 return $rc;
263 :     }
264 :    
265 : parrello 1.210 =head3 enqueue_similarities
266 : olson 1.93
267 :     usage: enqueue_similarities(\@sims)
268 :    
269 :     Queue the passed fids (a reference to a list) for similarity
270 :     computation.
271 :    
272 :     =cut
273 : parrello 1.210 #: Return Type ;
274 : olson 1.93 sub enqueue_similarities {
275 : efrank 1.1 my($fids) = @_;
276 :     my $fid;
277 :    
278 : olson 1.93 my $sim_q = "$FIG_Config::global/queued_similarities";
279 :    
280 :     open(TMP,">>$sim_q")
281 :     || die "could not open $sim_q";
282 :    
283 :     #
284 :     # We need to lock here so that if a computation is creating a snapshot of the
285 :     # queue, we block until it's done.
286 :     #
287 :    
288 :     flock(TMP, LOCK_EX) or die "Cannot lock $sim_q\n";
289 :    
290 : efrank 1.1 foreach $fid (@$fids)
291 :     {
292 :     print TMP "$fid\n";
293 :     }
294 :     close(TMP);
295 : olson 1.10 }
296 :    
297 : olson 1.281 =head3 export_similarity_request
298 :    
299 :     Creates a similarity computation request from the queued similarities and
300 :     the current NR.
301 :    
302 :     We keep track of the exported requests in case one gets lost.
303 :    
304 :     =cut
305 :    
306 :     sub export_similarity_request
307 :     {
308 :     my($self, $nr_file, $fasta_file) = @_;
309 :    
310 :     my $req_dir = "$FIG_Config::fig/var/sim_requests";
311 :     &verify_dir("$FIG_Config::fig/var");
312 :     &verify_dir($req_dir);
313 :    
314 :     $req_dir = "$req_dir/" . time;
315 :     &verify_dir($req_dir);
316 :    
317 :     #
318 :     # Open all of our output files before zeroing out the sim queue, in case
319 :     # there is a problem.
320 :     #
321 :    
322 :     open(my $user_fasta_fh, ">$fasta_file") or confess "Cannot open $fasta_file for writing: $!";
323 :     open(my $fasta_fh, ">$req_dir/fasta.in");
324 :    
325 :     open(my $user_nr_fh, ">$nr_file") or confess "Cannot open $nr_file for writing: $!";
326 :     open(my $nr_fh, ">$req_dir/nr") or confess "Cannot open $req_dir/nr for writing: $!";
327 :    
328 :     open(my $nr_read_fh, "<$FIG_Config::data/Global/nr") or die "Cannot open $FIG_Config::data/Global/nr for reading: $!";
329 :    
330 :     my $sim_q = "$FIG_Config::global/queued_similarities";
331 :    
332 :     #
333 :     # We need to lock here so that if a computation is creating a snapshot of the
334 :     # queue, we block until it's done.
335 :     #
336 :    
337 :     open(my $sim_q_lock, ">>$sim_q") or confess "could not open $sim_q";
338 :     flock($sim_q_lock, LOCK_EX) or confess "Cannot lock $sim_q\n";
339 :    
340 :     #
341 :     # Everything open & locked, start copying.
342 :     #
343 :    
344 :     copy("$sim_q", "$req_dir/q") or confess "Copy $sim_q $req_dir/q failed: $!";
345 :    
346 :     my($buf);
347 :     while (1)
348 :     {
349 :     my $n = read($nr_read_fh, $buf, 4096);
350 :     defined($n) or confess "Error reading nr: $!";
351 :     last unless $n;
352 :     syswrite($user_nr_fh, $buf) or confess "Error writing $nr_file: $!";
353 :     syswrite($nr_fh, $buf) or confess "Error writing $req_dir/nr: $!";
354 :     }
355 :    
356 :     close($nr_read_fh);
357 :     close($nr_fh);
358 :     close($user_nr_fh);
359 :    
360 :     #
361 :     # We can zero out the queue and unlock now.
362 :     #
363 :    
364 :     open(F, ">$sim_q") or die "Cannot open $sim_q to truncate it: $!\n";
365 :     close(F);
366 :    
367 :     close($sim_q_lock);
368 :    
369 :     #
370 :     # Generate the fasta input from the queued ids.
371 :     #
372 :    
373 :     open(my $q_fh, "<$req_dir/q");
374 :     while (my $id = <$q_fh>)
375 :     {
376 :     chomp $id;
377 :    
378 :     my $seq = $self->get_translation($id);
379 :    
380 :     display_id_and_seq($id, \$seq, $user_fasta_fh);
381 :     display_id_and_seq($id, \$seq, $fasta_fh);
382 :     }
383 :     close($q_fh);
384 :    
385 :     close($user_fasta_fh);
386 :     close($fasta_fh);
387 :     }
388 :    
389 : parrello 1.210 =head3 create_sim_askfor_pool
390 : olson 1.93
391 :     usage: create_sim_askfor_pool()
392 :    
393 : olson 1.123 Creates an askfor pool, a snapshot of the current NR and similarity
394 :     queue. Zeros out the old queue.
395 :    
396 :     The askfor pool needs to keep track of which sequences need to be
397 :     calculated, which have been handed out, etc. To simplify this task we
398 : olson 1.279 chunk the sequences into fairly small numbers (20k characters) and
399 : olson 1.123 allocate work on a per-chunk basis. We make use of the relational
400 :     database to keep track of chunk status as well as the seek locations
401 :     into the file of sequence data. The initial creation of the pool
402 :     involves indexing the sequence data with seek offsets and lengths and
403 :     populating the sim_askfor_index table with this information and with
404 :     initial status information.
405 : olson 1.93
406 : parrello 1.200 =cut
407 : parrello 1.210 #: Return Type $;
408 : olson 1.93 sub create_sim_askfor_pool
409 :     {
410 : olson 1.123 my($self, $chunk_size) = @_;
411 :    
412 : olson 1.279 $chunk_size = 20000 unless $chunk_size =~ /^\d+$/;
413 : olson 1.93
414 : olson 1.279 my $pool_dir = "$FIG_Config::fig/var/sim_pools";
415 : olson 1.93 &verify_dir($pool_dir);
416 :    
417 :     #
418 :     # Lock the pool directory.
419 :     #
420 :     open(my $lock, ">$pool_dir/lockfile");
421 :    
422 :     flock($lock, LOCK_EX);
423 :    
424 :     my $num = 0;
425 :     if (open(my $toc, "<$pool_dir/TOC"))
426 :     {
427 :     while (<$toc>)
428 :     {
429 :     chomp;
430 : olson 1.123 # print STDERR "Have toc entry $_\n";
431 : olson 1.93 my ($idx, $time, $str) = split(/\s+/, $_, 3);
432 :    
433 :     $num = max($num, $idx);
434 :     }
435 :     close($toc);
436 :     }
437 :     $num++;
438 :     open(my $toc, ">>$pool_dir/TOC") or die "Cannot write $pool_dir/TOC: $!\n";
439 :    
440 :     print $toc "$num ", time(), " New toc entry\n";
441 :     close($toc);
442 :    
443 : olson 1.123 my $cpool_id = sprintf "%04d", $num;
444 :     my $cpool_dir = "$pool_dir/$cpool_id";
445 : olson 1.93
446 :     #
447 :     # All set, create the directory for this pool.
448 :     #
449 :    
450 :     &verify_dir($cpool_dir);
451 :    
452 :     #
453 :     # Now we can copy the nr and sim queue here.
454 :     # Do this stuff inside an eval so we can clean up
455 :     # the lockfile.
456 :     #
457 :    
458 :     eval {
459 :     my $sim_q = "$FIG_Config::global/queued_similarities";
460 :    
461 : olson 1.123 copy("$sim_q", "$cpool_dir/q");
462 :     copy("$FIG_Config::data/Global/nr", "$cpool_dir/nr");
463 : olson 1.93
464 :     open(F, ">$sim_q") or die "Cannot open $sim_q to truncate it: $!\n";
465 :     close(F);
466 :     };
467 : parrello 1.200
468 : olson 1.93 unlink("$pool_dir/lockfile");
469 :     close($lock);
470 : olson 1.123
471 :     #
472 :     # We've created our pool; we can now run the formatdb and
473 :     # extract the sequences for the blast run.
474 :     #
475 :     my $child_pid = $self->run_in_background(sub {
476 :     #
477 :     # Need to close db or there's all sorts of trouble.
478 :     #
479 :    
480 :     my $cmd = "$FIG_Config::ext_bin/formatdb -i $cpool_dir/nr -p T -l $cpool_dir/formatdb.log";
481 :     print "Will run '$cmd'\n";
482 :     &run($cmd);
483 :     print "finished. Logfile:\n";
484 :     print &FIG::file_read("$cpool_dir/formatdb.log");
485 :     unlink("$cpool_dir/formatdb.pid");
486 : olson 1.279 });
487 :     warn "Running formatdb in background job $child_pid\n";
488 : olson 1.123 open(FPID, ">$cpool_dir/formatdb.pid");
489 :     print FPID "$child_pid\n";
490 :     close(FPID);
491 :    
492 :     my $db = $self->db_handle();
493 :     if (!$db->table_exists("sim_queue"))
494 :     {
495 :     $db->create_table(tbl => "sim_queue",
496 :     flds => "qid varchar(32), chunk_id INTEGER, seek INTEGER, len INTEGER, " .
497 :     "assigned BOOL, finished BOOL, output_file varchar(255), " .
498 :     "assignment_expires INTEGER, worker_info varchar(255)"
499 :     );
500 :     }
501 :    
502 :     #
503 :     # Write the fasta input file. Keep track of how many have been written,
504 :     # and write seek info into the database as appropriate.
505 :     #
506 :    
507 :     open(my $seq_fh, ">$cpool_dir/fasta.in");
508 :    
509 :     my($chunk_idx, $chunk_begin, $seq_idx);
510 :    
511 : olson 1.279 my $cur_size = 0;
512 :    
513 : olson 1.123 $chunk_idx = 0;
514 :     $chunk_begin = 0;
515 :     $seq_idx = 0;
516 :    
517 :     my(@seeks);
518 :    
519 : olson 1.279 my $tmpfile = "$FIG_Config::temp/simseek.$$";
520 :     open(my $tmpfh, ">$tmpfile") or confess "Cannot open tmpfile $tmpfile: $!";
521 :    
522 : olson 1.123 open(my $q_fh, "<$cpool_dir/q");
523 :     while (my $id = <$q_fh>)
524 :     {
525 :     chomp $id;
526 :    
527 :     my $seq = $self->get_translation($id);
528 :    
529 :     #
530 :     # check if we're at the beginning of a chunk
531 :     #
532 :    
533 :     print $seq_fh ">$id\n$seq\n";
534 :    
535 :     #
536 :     # Check if we're at the end of a chunk
537 :     #
538 :    
539 : olson 1.279 $cur_size += length($seq);
540 :     if ($cur_size >= $chunk_size)
541 : olson 1.123 {
542 :     my $chunk_end = tell($seq_fh);
543 :     my $chunk_len = $chunk_end - $chunk_begin;
544 :    
545 :     push(@seeks, [$cpool_id, $chunk_idx, $chunk_begin, $chunk_len]);
546 : olson 1.279 print $tmpfh join("\t", $cpool_id, $chunk_idx, $chunk_begin, $chunk_len, 'FALSE', 'FALSE'), "\n";
547 : olson 1.123 $chunk_idx++;
548 :     $chunk_begin = $chunk_end;
549 : olson 1.279 $cur_size = 0;
550 : olson 1.123 }
551 :     $seq_idx++;
552 :     }
553 :    
554 : olson 1.279 if ($cur_size > 0)
555 : olson 1.123 {
556 :     my $chunk_end = tell($seq_fh);
557 :     my $chunk_len = $chunk_end - $chunk_begin;
558 :    
559 : olson 1.279 print $tmpfh join("\t", $cpool_id, $chunk_idx, $chunk_begin, $chunk_len, 'FALSE', 'FALSE'), "\n";
560 : olson 1.123 push(@seeks, [$cpool_id, $chunk_idx, $chunk_begin, $chunk_len]);
561 :     }
562 :    
563 :     close($q_fh);
564 :     close($seq_fh);
565 : olson 1.279 close($tmpfh);
566 : olson 1.123
567 : olson 1.279 warn "Write seqs from $tmpfile\n";
568 : olson 1.123
569 : olson 1.279 $self->db_handle->load_table(tbl => 'sim_queue',
570 :     file => $tmpfile);
571 : parrello 1.200
572 : olson 1.279 unlink($tmpfile);
573 :    
574 :     # for my $seek (@seeks)
575 :     # {
576 :     # my($cpool_id, $chunk_idx, $chunk_begin, $chunk_len) = @$seek;
577 :    
578 :     # $db->SQL("insert into sim_queue (qid, chunk_id, seek, len, assigned, finished) " .
579 :     # "values('$cpool_id', $chunk_idx, $chunk_begin, $chunk_len, FALSE, FALSE)");
580 :     # }
581 : parrello 1.200
582 : olson 1.123 return $cpool_id;
583 :     }
584 :    
585 : parrello 1.210 #=head3 get_sim_queue
586 :     #
587 :     #usage: get_sim_queue($pool_id, $all_sims)
588 :     #
589 :     #Returns the sims in the given pool. If $all_sims is true, return the entire queue. Otherwise,
590 :     #just return the sims awaiting processing.
591 :     #
592 :     #=cut
593 : olson 1.123
594 :     sub get_sim_queue
595 :     {
596 :     my($self, $pool_id, $all_sims) = @_;
597 : olson 1.279
598 :    
599 :     }
600 :    
601 :    
602 :    
603 :     =head3 get_sim_work
604 :    
605 :     Get the next piece of sim computation work to be performed. Returned are
606 :     the path to the NR and a string containing the fasta data.
607 :    
608 :     =cut
609 :    
610 :     sub get_sim_work
611 :     {
612 :     my($self) = @_;
613 :    
614 :     #
615 :     # For now, just don't care about order of data that we get back.
616 :     #
617 :    
618 :     my $db = $self->db_handle();
619 :     my $lock = FIG::SimLock->new;
620 :    
621 :     my $work = $db->SQL(qq(SELECT qid, chunk_id, seek, len
622 :     FROM sim_queue
623 :     WHERE not finished
624 :     LIMIT 1));
625 :     print "Got work ", Dumper($work), "\n";
626 :    
627 :     if (not $work or @$work == 0)
628 :     {
629 :     return undef;
630 :     }
631 :    
632 :     my($cpool_id, $chunk_id, $seek, $len) = @{$work->[0]};
633 :    
634 :     my $pool_dir = "$FIG_Config::fig/var/sim_pools";
635 :     my $cpool_dir = "$pool_dir/$cpool_id";
636 :    
637 :     my $nr = "$cpool_dir/nr";
638 :     open(my $fh, "<$cpool_dir/fasta.in");
639 :     seek($fh, $seek, 0);
640 :     my $fasta;
641 :     read($fh, $fasta, $len);
642 :    
643 :     return($cpool_id, $chunk_id, $nr, $fasta, "$cpool_dir/out.$chunk_id");
644 :     }
645 :    
646 :     =head3 sim_work_done
647 :    
648 :     Declare that the work in pool_id/chunk_id has been completed, and output written
649 :     to the pool directory (get_sim_work gave it the path).
650 :    
651 :     =cut
652 :    
653 :     sub sim_work_done
654 :     {
655 :     my($self, $pool_id, $chunk_id, $out_file) = @_;
656 :    
657 :     if (! -f $out_file)
658 :     {
659 :     confess "sim_work_done: output file $out_file does not exist";
660 :     }
661 :    
662 :     my $db = $self->db_handle();
663 :     my $lock = FIG::SimLock->new;
664 :    
665 :     my $dbh = $db->{_dbh};
666 :    
667 :     my $rows = $dbh->do(qq(UPDATE sim_queue
668 :     SET finished = TRUE, output_file = ?
669 :     WHERE qid = ? and chunk_id = ?), undef, $out_file, $pool_id, $chunk_id);
670 :     if ($rows != 1)
671 :     {
672 :     if ($dbh->errstr)
673 :     {
674 :     confess "Update not able to set finished=TRUE: ", $dbh->errstr;
675 :     }
676 :     else
677 :     {
678 :     confess "Update not able to set finished=TRUE";
679 :     }
680 :     }
681 :    
682 :     #
683 :     # Determine if this was the last piece of work for this pool. If so, we can
684 :     # schedule the postprocessing work.
685 :     #
686 :     # Note we're still holding the lock.
687 :     #
688 :    
689 :     my $out = $db->SQL(qq(SELECT chunk_id
690 :     FROM sim_queue
691 :     WHERE qid = ? AND not finished), undef, $pool_id);
692 :     if (@$out == 0)
693 :     {
694 :     #
695 :     # Pool is done.
696 :     #
697 :     $self->schedule_sim_pool_postprocessing($pool_id);
698 :     }
699 : olson 1.123 }
700 :    
701 : olson 1.279 =head3 schedule_sim_pool_postprocessing
702 :    
703 :     Schedule a job to do the similarity postprocessing for pool $pool_id.
704 :    
705 :    
706 :     =cut
707 :    
708 :     sub schedule_sim_pool_postprocessing
709 :     {
710 :     my($self, $pool_id) = @_;
711 :    
712 :     my $pool_dir = "$FIG_Config::fig/var/sim_pools";
713 :     my $cpool_dir = "$pool_dir/$pool_id";
714 :    
715 :     my $js = JobScheduler->new();
716 :     my $job = $js->job_create();
717 :    
718 :     my $spath = $job->get_script_path();
719 :     open(my $sfh, ">$spath");
720 :     print $sfh <<END;
721 :     #!/bin/sh
722 :     . $FIG_Config::fig_disk/config/fig-user-env.sh
723 :     $FIG_Config::bin/postprocess_computed_sims $pool_id
724 :     END
725 :    
726 :     close($sfh);
727 :     chmod(0775, $spath);
728 :    
729 :     #
730 :     # Write the job ID to the subsystem queue dir.
731 :     #
732 :    
733 :     open(J, ">$cpool_dir/postprocess_jobid");
734 :     print J $job->get_id(), "\n";
735 :     close(J);
736 :    
737 :     $job->enqueue();
738 :     }
739 :    
740 :     =head3 postprocess_computed_sims
741 :    
742 :     We build a pipe to this pipeline:
743 :    
744 :     reduce_sims peg.synonyms 300 | reformat_sims nr | split_sims dest prefix
745 :    
746 :     We put the new sims in the pool directory, and then copy to NewSims.
747 :    
748 :     =cut
749 :    
750 :     sub postprocess_computed_sims
751 :     {
752 :     my($self, $pool_id) = @_;
753 :    
754 :     #
755 :     # We don't lock here because the job is already done, and we
756 :     # shouldn't (ha, ha) ever postprocess twice.
757 :     #
758 :    
759 :     my $pool_dir = "$FIG_Config::fig/var/sim_pools";
760 :     my $cpool_dir = "$pool_dir/$pool_id";
761 :    
762 :     my $sim_dir = "$cpool_dir/NewSims";
763 :     &verify_dir($sim_dir);
764 :    
765 :     #
766 :     # Open the processing pipeline.
767 :     #
768 :    
769 :     my $reduce = "$FIG_Config::bin/reduce_sims $FIG_Config::global/peg.synonyms 300";
770 :     my $reformat = "$FIG_Config::bin/reformat_sims $cpool_dir/nr";
771 :     my $split = "$FIG_Config::bin/split_sims $sim_dir sims.$pool_id";
772 :     open(my $process, "| $reduce | $reformat | $split");
773 :    
774 :     #
775 :     # Iterate over all the sims files, taken from the database.
776 :     #
777 :    
778 :     my $dbh = $self->db_handle()->{_dbh};
779 :     my $files = $dbh->selectcol_arrayref(qq(SELECT output_file
780 :     FROM sim_queue
781 :     WHERE qid = ? and output_file IS NOT NULL
782 :     ORDER BY chunk_id), undef, $pool_id);
783 :     for my $file (@$files)
784 :     {
785 :     my $buf;
786 :     open(my $fh, "<$file") or confess "Cannot sim input file $file: $!";
787 :     while (read($fh, $buf, 4096))
788 :     {
789 :     print $process $buf;
790 :     }
791 :     close($fh);
792 :     }
793 :     my $res = close($process);
794 :     if (!$res)
795 :     {
796 :     if ($!)
797 :     {
798 :     confess "Error closing process pipeline: $!";
799 :     }
800 :     else
801 :     {
802 :     confess "Process pipeline exited with status $?";
803 :     }
804 :     }
805 :    
806 :     #
807 :     # If we got here, it worked. Copy the new sims files over to NewSims.
808 :     #
809 :    
810 :     opendir(my $simdh, $sim_dir) or confess "Cannot open $sim_dir: $!";
811 :     my @new_sims = grep { $_ !~ /^\./ } readdir($simdh);
812 :     closedir($simdh);
813 :    
814 :     &verify_dir("$FIG_Config::data/NewSims");
815 :    
816 :     for my $sim_file (@new_sims)
817 :     {
818 :     my $target = "$FIG_Config::data/NewSims/$sim_file";
819 :     if (-s $target)
820 :     {
821 :     confess "$target already exists";
822 :     }
823 :     print "copying sim file $sim_file\n";
824 :     &FIG::run("cp $sim_dir/$sim_file $target");
825 :     &FIG::run("$FIG_Config::bin/index_sims $target");
826 :     }
827 :     }
828 :    
829 :    
830 : parrello 1.210 =head3 get_active_sim_pools
831 : olson 1.123
832 :     usage: get_active_sim_pools()
833 :    
834 :     Return a list of the pool id's for the sim processing queues that have entries awaiting
835 :     computation.
836 :    
837 :     =cut
838 : parrello 1.210 #: Return Type @;
839 : olson 1.123 sub get_active_sim_pools
840 :     {
841 :     my($self) = @_;
842 :    
843 :     my $dbh = $self->db_handle();
844 :    
845 :     my $res = $dbh->SQL("select distinct qid from sim_queue where not finished");
846 :     return undef unless $res;
847 :    
848 :     return map { $_->[0] } @$res;
849 :     }
850 :    
851 : parrello 1.210 =head3 get_sim_pool_info
852 : olson 1.123
853 :     usage: get_sim_pool_info($pool_id)
854 :    
855 :     Return information about the given sim pool. Return value
856 :     is a list ($total_entries, $n_finished, $n_assigned, $n_unassigned)
857 :    
858 :     =cut
859 : parrello 1.210 #: Return Type @;
860 : olson 1.123 sub get_sim_pool_info
861 :     {
862 :     my($self, $pool_id) = @_;
863 :     my($dbh, $res, $total_entries, $n_finished, $n_assigned, $n_unassigned);
864 :    
865 :     $dbh = $self->db_handle();
866 :    
867 :     $res = $dbh->SQL("select count(chunk_id) from sim_queue where qid = '$pool_id'");
868 : parrello 1.200 $total_entries = $res->[0]->[0];
869 : olson 1.123
870 :     $res = $dbh->SQL("select count(chunk_id) from sim_queue where qid = '$pool_id' and finished");
871 :     $n_finished = $res->[0]->[0];
872 :    
873 :     $res = $dbh->SQL("select count(chunk_id) from sim_queue where qid = '$pool_id' and assigned and not finished");
874 :     $n_assigned = $res->[0]->[0];
875 :    
876 :     $res = $dbh->SQL("select count(chunk_id) from sim_queue where qid = '$pool_id' and not finished and not assigned");
877 :     $n_unassigned = $res->[0]->[0];
878 :    
879 :     return ($total_entries, $n_finished, $n_assigned, $n_unassigned);
880 : olson 1.93 }
881 :    
882 : parrello 1.210 #=head3 get_sim_chunk
883 :     #
884 :     #usage: get_sim_chunk($n_seqs, $worker_id)
885 :     #
886 :     #Returns a chunk of $n_seqs of work.
887 :     #
888 :     #From Ross, about how sims are processed:
889 :     #
890 :     #Here is how I process them:
891 :     #
892 :     #
893 :     # bash$ cd /Volumes/seed/olson/Sims/June22.out
894 :     # bash$ for i in really*
895 :     # > do
896 :     # > cat < $i >> /Volumes/laptop/new.sims
897 :     # > done
898 :     #
899 :     #
900 :     #Then, I need to "reformat" them by adding to columns to each one
901 :     # and split the result into files of about 3M each This I do using
902 :     #
903 :     #reduce_sims /Volumes/laptop/NR/NewNR/peg.synonyms.june21 300 < /Volumes/laptop/new.sims |
904 :     # reformat_sims /Volumes/laptop/NR/NewNR/checked.nr.june21 > /Volumes/laptop/reformated.sims
905 :     #rm /Volumes/laptop/new.sims
906 :     #split_sims /Volumes/laptop/NewSims sims.june24 reformated.sims
907 :     #rm reformatted.sims
908 :     #
909 :     #=cut
910 : olson 1.93
911 : parrello 1.210 sub get_sim_chunk
912 :     {
913 :     my($self, $n_seqs, $worker_id) = @_;
914 : olson 1.93
915 :    
916 : parrello 1.210 }
917 : olson 1.123
918 : parrello 1.210 =head3 get_local_hostname
919 : parrello 1.200
920 : parrello 1.213 usage: my $result = $fig->get_local_hostname();
921 : olson 1.123
922 : olson 1.93 =cut
923 : parrello 1.213 #: Return Type $;
924 : olson 1.10 sub get_local_hostname {
925 : olson 1.52
926 :     #
927 :     # See if there is a FIGdisk/config/hostname file. If there
928 :     # is, force the hostname to be that.
929 :     #
930 :    
931 :     my $hostfile = "$FIG_Config::fig_disk/config/hostname";
932 :     if (-f $hostfile)
933 :     {
934 :     my $fh;
935 :     if (open($fh, $hostfile))
936 :     {
937 :     my $hostname = <$fh>;
938 :     chomp($hostname);
939 :     return $hostname;
940 :     }
941 :     }
942 : parrello 1.200
943 : olson 1.10 #
944 :     # First check to see if we our hostname is correct.
945 :     #
946 :     # Map it to an IP address, and try to bind to that ip.
947 :     #
948 :    
949 :     my $tcp = getprotobyname('tcp');
950 : parrello 1.200
951 : olson 1.10 my $hostname = `hostname`;
952 : golsen 1.44 chomp($hostname);
953 : olson 1.10
954 :     my @hostent = gethostbyname($hostname);
955 :    
956 :     if (@hostent > 0)
957 :     {
958 :     my $sock;
959 :     my $ip = $hostent[4];
960 : parrello 1.200
961 : olson 1.10 socket($sock, PF_INET, SOCK_STREAM, $tcp);
962 :     if (bind($sock, sockaddr_in(0, $ip)))
963 :     {
964 :     #
965 :     # It worked. Reverse-map back to a hopefully fqdn.
966 :     #
967 :    
968 :     my @rev = gethostbyaddr($ip, AF_INET);
969 :     if (@rev > 0)
970 :     {
971 : olson 1.28 my $host = $rev[0];
972 :     #
973 :     # Check to see if we have a FQDN.
974 :     #
975 :    
976 :     if ($host =~ /\./)
977 :     {
978 :     #
979 :     # Good.
980 :     #
981 :     return $host;
982 :     }
983 :     else
984 :     {
985 :     #
986 :     # We didn't get a fqdn; bail and return the IP address.
987 :     #
988 :     return get_hostname_by_adapter()
989 :     }
990 : olson 1.10 }
991 :     else
992 :     {
993 :     return inet_ntoa($ip);
994 :     }
995 :     }
996 :     else
997 :     {
998 :     #
999 :     # Our hostname must be wrong; we can't bind to the IP
1000 :     # address it maps to.
1001 :     # Return the name associated with the adapter.
1002 :     #
1003 :     return get_hostname_by_adapter()
1004 :     }
1005 :     }
1006 :     else
1007 :     {
1008 :     #
1009 :     # Our hostname isn't known to DNS. This isn't good.
1010 :     # Return the name associated with the adapter.
1011 :     #
1012 :     return get_hostname_by_adapter()
1013 :     }
1014 :     }
1015 :    
1016 : parrello 1.213 =head3 get_hostname_by_adapter
1017 :    
1018 :     usage: my $name = $fig->get_hostname_by_adapter();
1019 :    
1020 :     =cut
1021 :     #: Return Type $;
1022 : olson 1.10 sub get_hostname_by_adapter {
1023 : olson 1.111 shift if UNIVERSAL::isa($_[0],__PACKAGE__);
1024 : olson 1.10 #
1025 :     # Attempt to determine our local hostname based on the
1026 :     # network environment.
1027 :     #
1028 :     # This implementation reads the routing table for the default route.
1029 :     # We then look at the interface config for the interface that holds the default.
1030 :     #
1031 :     #
1032 :     # Linux routing table:
1033 :     # [olson@yips 0.0.0]$ netstat -rn
1034 :     # Kernel IP routing table
1035 :     # Destination Gateway Genmask Flags MSS Window irtt Iface
1036 :     # 140.221.34.32 0.0.0.0 255.255.255.224 U 0 0 0 eth0
1037 :     # 169.254.0.0 0.0.0.0 255.255.0.0 U 0 0 0 eth0
1038 :     # 127.0.0.0 0.0.0.0 255.0.0.0 U 0 0 0 lo
1039 :     # 0.0.0.0 140.221.34.61 0.0.0.0 UG 0 0 0 eth0
1040 : parrello 1.200 #
1041 : olson 1.10 # Mac routing table:
1042 : parrello 1.200 #
1043 : olson 1.10 # bash-2.05a$ netstat -rn
1044 :     # Routing tables
1045 : parrello 1.200 #
1046 : olson 1.10 # Internet:
1047 :     # Destination Gateway Flags Refs Use Netif Expire
1048 :     # default 140.221.11.253 UGSc 12 120 en0
1049 :     # 127.0.0.1 127.0.0.1 UH 16 8415486 lo0
1050 :     # 140.221.8/22 link#4 UCS 12 0 en0
1051 :     # 140.221.8.78 0:6:5b:f:51:c4 UHLW 0 183 en0 408
1052 :     # 140.221.8.191 0:3:93:84:ab:e8 UHLW 0 92 en0 622
1053 :     # 140.221.8.198 0:e0:98:8e:36:e2 UHLW 0 5 en0 691
1054 :     # 140.221.9.6 0:6:5b:f:51:d6 UHLW 1 63 en0 1197
1055 :     # 140.221.10.135 0:d0:59:34:26:34 UHLW 2 2134 en0 1199
1056 :     # 140.221.10.152 0:30:1b:b0:ec:dd UHLW 1 137 en0 1122
1057 :     # 140.221.10.153 127.0.0.1 UHS 0 0 lo0
1058 :     # 140.221.11.37 0:9:6b:53:4e:4b UHLW 1 624 en0 1136
1059 :     # 140.221.11.103 0:30:48:22:59:e6 UHLW 3 973 en0 1016
1060 :     # 140.221.11.224 0:a:95:6f:7:10 UHLW 1 1 en0 605
1061 :     # 140.221.11.237 0:1:30:b8:80:c0 UHLW 0 0 en0 1158
1062 :     # 140.221.11.250 0:1:30:3:1:0 UHLW 0 0 en0 1141
1063 :     # 140.221.11.253 0:d0:3:e:70:a UHLW 13 0 en0 1199
1064 :     # 169.254 link#4 UCS 0 0 en0
1065 : parrello 1.200 #
1066 : olson 1.10 # Internet6:
1067 :     # Destination Gateway Flags Netif Expire
1068 :     # UH lo0
1069 :     # fe80::%lo0/64 Uc lo0
1070 :     # link#1 UHL lo0
1071 :     # fe80::%en0/64 link#4 UC en0
1072 :     # 0:a:95:a8:26:68 UHL lo0
1073 :     # ff01::/32 U lo0
1074 :     # ff02::%lo0/32 UC lo0
1075 :     # ff02::%en0/32 link#4 UC en0
1076 :    
1077 :     my($fh);
1078 :    
1079 :     if (!open($fh, "netstat -rn |"))
1080 :     {
1081 :     warn "Cannot run netstat to determine local IP address\n";
1082 :     return "localhost";
1083 :     }
1084 :    
1085 :     my $interface_name;
1086 : parrello 1.200
1087 : olson 1.10 while (<$fh>)
1088 :     {
1089 :     my @cols = split();
1090 :    
1091 :     if ($cols[0] eq "default" || $cols[0] eq "0.0.0.0")
1092 :     {
1093 :     $interface_name = $cols[$#cols];
1094 :     }
1095 :     }
1096 :     close($fh);
1097 : parrello 1.200
1098 : olson 1.11 # print "Default route on $interface_name\n";
1099 : olson 1.10
1100 :     #
1101 :     # Find ifconfig.
1102 :     #
1103 :    
1104 :     my $ifconfig;
1105 :    
1106 :     for my $dir ((split(":", $ENV{PATH}), "/sbin", "/usr/sbin"))
1107 :     {
1108 :     if (-x "$dir/ifconfig")
1109 :     {
1110 :     $ifconfig = "$dir/ifconfig";
1111 :     last;
1112 :     }
1113 :     }
1114 :    
1115 :     if ($ifconfig eq "")
1116 :     {
1117 :     warn "Ifconfig not found\n";
1118 :     return "localhost";
1119 :     }
1120 : olson 1.11 # print "Foudn $ifconfig\n";
1121 : olson 1.10
1122 :     if (!open($fh, "$ifconfig $interface_name |"))
1123 :     {
1124 :     warn "Could not run $ifconfig: $!\n";
1125 :     return "localhost";
1126 :     }
1127 :    
1128 :     my $ip;
1129 :     while (<$fh>)
1130 :     {
1131 :     #
1132 :     # Mac:
1133 :     # inet 140.221.10.153 netmask 0xfffffc00 broadcast 140.221.11.255
1134 :     # Linux:
1135 :     # inet addr:140.221.34.37 Bcast:140.221.34.63 Mask:255.255.255.224
1136 :     #
1137 :    
1138 :     chomp;
1139 :     s/^\s*//;
1140 :    
1141 : olson 1.11 # print "Have '$_'\n";
1142 : olson 1.10 if (/inet\s+addr:(\d+\.\d+\.\d+\.\d+)\s+/)
1143 :     {
1144 :     #
1145 :     # Linux hit.
1146 :     #
1147 :     $ip = $1;
1148 : olson 1.11 # print "Got linux $ip\n";
1149 : olson 1.10 last;
1150 :     }
1151 :     elsif (/inet\s+(\d+\.\d+\.\d+\.\d+)\s+/)
1152 :     {
1153 :     #
1154 :     # Mac hit.
1155 :     #
1156 :     $ip = $1;
1157 : olson 1.11 # print "Got mac $ip\n";
1158 : olson 1.10 last;
1159 :     }
1160 :     }
1161 :     close($fh);
1162 :    
1163 :     if ($ip eq "")
1164 :     {
1165 :     warn "Didn't find an IP\n";
1166 :     return "localhost";
1167 :     }
1168 :    
1169 :     return $ip;
1170 : efrank 1.1 }
1171 :    
1172 : parrello 1.213 =head3 get_seed_id
1173 :    
1174 :     usage: my $id = $fig->get_seed_id();
1175 :    
1176 :     =cut
1177 :     #: Return type $;
1178 : olson 1.38 sub get_seed_id {
1179 :     #
1180 :     # Retrieve the seed identifer from FIGdisk/config/seed_id.
1181 :     #
1182 :     # If it's not there, create one, and make it readonly.
1183 :     #
1184 :    
1185 : olson 1.111 shift if UNIVERSAL::isa($_[0],__PACKAGE__);
1186 : olson 1.38 my $id;
1187 :     my $id_file = "$FIG_Config::fig_disk/config/seed_id";
1188 :     if (! -f $id_file)
1189 :     {
1190 :     my $newid = `uuidgen`;
1191 :     if (!$newid)
1192 :     {
1193 :     die "Cannot run uuidgen: $!";
1194 :     }
1195 :    
1196 :     chomp($newid);
1197 :     my $fh = new FileHandle(">$id_file");
1198 :     if (!$fh)
1199 :     {
1200 :     die "error creating $id_file: $!";
1201 :     }
1202 :     print $fh "$newid\n";
1203 :     $fh->close();
1204 :     chmod(0444, $id_file);
1205 :     }
1206 :     my $fh = new FileHandle("<$id_file");
1207 :     $id = <$fh>;
1208 :     chomp($id);
1209 :     return $id;
1210 :     }
1211 :    
1212 : olson 1.155 =pod
1213 :    
1214 :     =head1 get_release_info
1215 :    
1216 : olson 1.195 Return the current data release information. It is returned as the list
1217 : parrello 1.200 ($name, $id, $inst, $email, $parent_id, $description).
1218 : olson 1.195
1219 :     The release info comes from the file FIG/Data/RELEASE. It is formatted as:
1220 :    
1221 :     <release-name>
1222 :     <unique id>
1223 :     <institution>
1224 :     <contact email>
1225 :     <unique id of data release this release derived from>
1226 :     <description>
1227 :    
1228 :     For instance:
1229 :     -----
1230 :     SEED Data Release, 09/15/2004.
1231 :     4148208C-1DF2-11D9-8417-000A95D52EF6
1232 :     ANL/FIG
1233 :     olson@mcs.anl.gov
1234 :    
1235 :     Test release.
1236 :     -----
1237 :    
1238 :     If no RELEASE file exists, this routine will create one with a new unique ID. This
1239 :     lets a peer optimize the data transfer by being able to cache ID translations
1240 :     from this instance.
1241 : olson 1.155
1242 :     =cut
1243 : parrello 1.213 #: Return Type @;
1244 : olson 1.155 sub get_release_info
1245 :     {
1246 : olson 1.196 my($fig, $no_create) = @_;
1247 : olson 1.195
1248 :     my $rel_file = "$FIG_Config::data/RELEASE";
1249 :    
1250 : olson 1.196 if (! -f $rel_file and !$no_create)
1251 : olson 1.195 {
1252 :     #
1253 :     # Create a new one.
1254 :     #
1255 :    
1256 :     my $newid = `uuidgen`;
1257 :     if (!$newid)
1258 :     {
1259 :     die "Cannot run uuidgen: $!";
1260 :     }
1261 :    
1262 :     chomp($newid);
1263 :    
1264 :     my $relinfo = "Automatically generated release info " . localtime();
1265 :     my $inst = "Unknown";
1266 :     my $contact = "Unknown";
1267 :     my $parent = "";
1268 :     my( $a, $b, $e, $v, $env ) = $fig->genome_counts;
1269 :     my $description = "Automatically generated release info\n";
1270 :     $description .= "Contains $a archaeal, $b bacterial, $e eukaryal, $v viral and $env environmental genomes.\n";
1271 :    
1272 :     my $fh = new FileHandle(">$rel_file");
1273 :     if (!$fh)
1274 :     {
1275 :     warn "error creating $rel_file: $!";
1276 :     return undef;
1277 :     }
1278 :     print $fh "$relinfo\n";
1279 :     print $fh "$newid\n";
1280 :     print $fh "$inst\n";
1281 :     print $fh "$contact\n";
1282 :     print $fh "$parent\n";
1283 :     print $fh $description;
1284 :     $fh->close();
1285 :     chmod(0444, $rel_file);
1286 :     }
1287 :    
1288 :     if (open(my $fh, $rel_file))
1289 :     {
1290 :     my(@lines) = <$fh>;
1291 :     close($fh);
1292 : parrello 1.200
1293 : olson 1.195 chomp(@lines);
1294 : parrello 1.200
1295 : olson 1.195 my($info, $id, $inst, $contact, $parent, @desc) = @lines;
1296 :    
1297 :     return ($info, $id, $inst, $contact, $parent, join("\n", @desc));
1298 :     }
1299 : olson 1.155
1300 :     return undef;
1301 :     }
1302 :    
1303 :     =pod
1304 :    
1305 :     =head1 get_peer_last_update
1306 :    
1307 : parrello 1.213 usage: my $date = $fig->get_peer_last_update($peer_id);
1308 :    
1309 : olson 1.155 Return the timestamp from the last successful peer-to-peer update with
1310 : parrello 1.200 the given peer.
1311 : olson 1.155
1312 :     We store this information in FIG/Data/Global/Peers/<peer-id>.
1313 :    
1314 :     =cut
1315 : parrello 1.213 #: Return Type $;
1316 : olson 1.155 sub get_peer_last_update
1317 :     {
1318 :     my($self, $peer_id) = @_;
1319 :    
1320 :     my $dir = "$FIG_Config::data/Global/Peers";
1321 :     &verify_dir($dir);
1322 :     $dir .= "/$peer_id";
1323 :     &verify_dir($dir);
1324 :    
1325 :     my $update_file = "$dir/last_update";
1326 :     if (-f $update_file)
1327 :     {
1328 :     my $time = file_head($update_file, 1);
1329 :     chomp $time;
1330 :     return $time;
1331 :     }
1332 :     else
1333 :     {
1334 :     return undef;
1335 :     }
1336 :     }
1337 :    
1338 : parrello 1.213 =pod
1339 :    
1340 :     =head1 set_peer_last_update
1341 :    
1342 :     usage: $fig->set_peer_last_update($peer_id, $time);
1343 :    
1344 :     =cut
1345 :     #: Return Type ;
1346 :    
1347 : olson 1.155 sub set_peer_last_update
1348 :     {
1349 :     my($self, $peer_id, $time) = @_;
1350 :    
1351 :     my $dir = "$FIG_Config::data/Global/Peers";
1352 :     &verify_dir($dir);
1353 :     $dir .= "/$peer_id";
1354 :     &verify_dir($dir);
1355 :    
1356 :     my $update_file = "$dir/last_update";
1357 :     open(F, ">$update_file");
1358 :     print F "$time\n";
1359 :     close(F);
1360 :     }
1361 :    
1362 : parrello 1.213 =head3 cgi_url
1363 :    
1364 :     usage: my $url = $fig->cgi_url();
1365 :    
1366 :     =cut
1367 :     #: Return Type $;
1368 : efrank 1.1 sub cgi_url {
1369 : olson 1.111 shift if UNIVERSAL::isa($_[0],__PACKAGE__);
1370 : efrank 1.1 return &plug_url($FIG_Config::cgi_url);
1371 :     }
1372 : parrello 1.200
1373 : parrello 1.213 =head3 temp_url
1374 :    
1375 :     usage: my $url = $fig->temp_url();
1376 :    
1377 :     =cut
1378 :     #: Return Type $;
1379 : efrank 1.1 sub temp_url {
1380 : olson 1.111 shift if UNIVERSAL::isa($_[0],__PACKAGE__);
1381 : efrank 1.1 return &plug_url($FIG_Config::temp_url);
1382 :     }
1383 : parrello 1.200
1384 : parrello 1.213 =head3 plug_url
1385 :    
1386 :     usage: my $url2 = $fig->plug_url($url);
1387 :    
1388 :     =cut
1389 :     #: Return Type $;
1390 : efrank 1.1 sub plug_url {
1391 : olson 1.111 shift if UNIVERSAL::isa($_[0],__PACKAGE__);
1392 : efrank 1.1 my($url) = @_;
1393 :    
1394 : golsen 1.44 my $name;
1395 :    
1396 :     # Revised by GJO
1397 :     # First try to get url from the current http request
1398 :    
1399 :     if ( defined( $ENV{ 'HTTP_HOST' } ) # This is where $cgi->url gets its value
1400 :     && ( $name = $ENV{ 'HTTP_HOST' } )
1401 :     && ( $url =~ s~^http://[^/]*~http://$name~ ) # ~ is delimiter
1402 :     ) {}
1403 :    
1404 :     # Otherwise resort to alternative sources
1405 :    
1406 :     elsif ( ( $name = &get_local_hostname )
1407 :     && ( $url =~ s~^http://[^/]*~http://$name~ ) # ~ is delimiter
1408 :     ) {}
1409 :    
1410 : efrank 1.1 return $url;
1411 :     }
1412 :    
1413 : parrello 1.213 =head3 file_read
1414 :    
1415 :     usage: my $text = $fig->file_read($fileName);
1416 :    
1417 :     usage: my @lines = $fig->file_read($fileName);
1418 :    
1419 :     =cut
1420 :     #: Return Type $;
1421 :     #: Return Type @;
1422 : olson 1.90 sub file_read
1423 :     {
1424 : olson 1.111 shift if UNIVERSAL::isa($_[0],__PACKAGE__);
1425 : olson 1.90 my($file) = @_;
1426 :    
1427 :     if (open(my $fh, "<$file"))
1428 :     {
1429 : parrello 1.200 if (wantarray)
1430 : olson 1.90 {
1431 :     my @ret = <$fh>;
1432 :     return @ret;
1433 :     }
1434 :     else
1435 :     {
1436 :     local $/;
1437 :     my $text = <$fh>;
1438 :     close($fh);
1439 :     return $text;
1440 :     }
1441 :     }
1442 :     }
1443 :    
1444 :    
1445 : parrello 1.213 =head3 file_head
1446 :    
1447 :     usage: my $text = $fig->file_read($fileName, $count);
1448 :    
1449 :     usage: my @lines = $fig->file_read($fileName, $count);
1450 :    
1451 :     =cut
1452 :     #: Return Type $;
1453 :     #: Return Type @;
1454 : olson 1.90 sub file_head
1455 :     {
1456 : olson 1.111 shift if UNIVERSAL::isa($_[0],__PACKAGE__);
1457 : olson 1.90 my($file, $n) = @_;
1458 :    
1459 :     if (!$n)
1460 :     {
1461 :     $n = 1;
1462 :     }
1463 :    
1464 :     if (open(my $fh, "<$file"))
1465 :     {
1466 :     my(@ret, $i);
1467 :    
1468 :     $i = 0;
1469 :     while (<$fh>)
1470 :     {
1471 :     push(@ret, $_);
1472 :     $i++;
1473 :     last if $i >= $n;
1474 :     }
1475 : olson 1.93 close($fh);
1476 : olson 1.155
1477 :     if (wantarray)
1478 :     {
1479 :     return @ret;
1480 :     }
1481 :     else
1482 :     {
1483 :     return join("", @ret);
1484 :     }
1485 : olson 1.90 }
1486 :     }
1487 :    
1488 :    
1489 : efrank 1.1 =pod
1490 :    
1491 :     =head1 hiding/caching in a FIG object
1492 :    
1493 :     We save the DB handle, cache taxonomies, and put a few other odds and ends in the
1494 :     FIG object. We expect users to invoke these services using the object $fig constructed
1495 :     using:
1496 :    
1497 :     use FIG;
1498 :     my $fig = new FIG;
1499 :    
1500 :     $fig is then used as the basic mechanism for accessing FIG services. It is, of course,
1501 :     just a hash that is used to retain/cache data. The most commonly accessed item is the
1502 :     DB filehandle, which is accessed via $self->db_handle.
1503 :    
1504 :     We cache genus/species expansions, taxonomies, distances (very crudely estimated) estimated
1505 : parrello 1.200 between genomes, and a variety of other things. I am not sure that using cached/2 was a
1506 : efrank 1.1 good idea, but I did it.
1507 :    
1508 :     =cut
1509 :    
1510 :     sub db_handle {
1511 :     my($self) = @_;
1512 :    
1513 :     return $self->{_dbf};
1514 :     }
1515 :    
1516 :     sub cached {
1517 :     my($self,$what) = @_;
1518 :    
1519 :     my $x = $self->{$what};
1520 :     if (! $x)
1521 :     {
1522 :     $x = $self->{$what} = {};
1523 :     }
1524 :     return $x;
1525 :     }
1526 :    
1527 :     ################ Basic Routines [ existed since WIT ] ##########################
1528 :    
1529 :    
1530 :     =pod
1531 :    
1532 :     =head1 min
1533 :    
1534 :     usage: $n = &FIG::min(@x)
1535 :    
1536 :     Assumes @x contains numeric values. Returns the minimum of the values.
1537 :    
1538 :     =cut
1539 : parrello 1.213 #: Return Type $;
1540 : efrank 1.1 sub min {
1541 : olson 1.111 shift if UNIVERSAL::isa($_[0],__PACKAGE__);
1542 : efrank 1.1 my(@x) = @_;
1543 :     my($min,$i);
1544 :    
1545 :     (@x > 0) || return undef;
1546 :     $min = $x[0];
1547 :     for ($i=1; ($i < @x); $i++)
1548 :     {
1549 :     $min = ($min > $x[$i]) ? $x[$i] : $min;
1550 :     }
1551 :     return $min;
1552 :     }
1553 :    
1554 :     =pod
1555 :    
1556 :     =head1 max
1557 :    
1558 :     usage: $n = &FIG::max(@x)
1559 :    
1560 :     Assumes @x contains numeric values. Returns the maximum of the values.
1561 :    
1562 :     =cut
1563 : parrello 1.213 #: Return Type $;
1564 : efrank 1.1 sub max {
1565 : olson 1.111 shift if UNIVERSAL::isa($_[0],__PACKAGE__);
1566 : efrank 1.1 my(@x) = @_;
1567 :     my($max,$i);
1568 :    
1569 :     (@x > 0) || return undef;
1570 :     $max = $x[0];
1571 :     for ($i=1; ($i < @x); $i++)
1572 :     {
1573 :     $max = ($max < $x[$i]) ? $x[$i] : $max;
1574 :     }
1575 :     return $max;
1576 :     }
1577 :    
1578 :     =pod
1579 :    
1580 :     =head1 between
1581 :    
1582 :     usage: &FIG::between($x,$y,$z)
1583 :    
1584 :     Returns true iff $y is between $x and $z.
1585 :    
1586 :     =cut
1587 : parrello 1.213 #: Return Type $;
1588 : efrank 1.1 sub between {
1589 : olson 1.111 shift if UNIVERSAL::isa($_[0],__PACKAGE__);
1590 : efrank 1.1 my($x,$y,$z) = @_;
1591 :    
1592 :     if ($x < $z)
1593 :     {
1594 :     return (($x <= $y) && ($y <= $z));
1595 :     }
1596 :     else
1597 :     {
1598 :     return (($x >= $y) && ($y >= $z));
1599 :     }
1600 :     }
1601 :    
1602 :     =pod
1603 :    
1604 :     =head1 standard_genetic_code
1605 :    
1606 :     usage: $code = &FIG::standard_genetic_code()
1607 :    
1608 :     Routines like "translate" can take a "genetic code" as an argument. I implemented such
1609 :     codes using hashes that assumed uppercase DNA triplets as keys.
1610 :    
1611 :     =cut
1612 : parrello 1.213 #: Return Type $;
1613 : efrank 1.1 sub standard_genetic_code {
1614 : olson 1.111 shift if UNIVERSAL::isa($_[0],__PACKAGE__);
1615 : parrello 1.200
1616 : efrank 1.1 my $code = {};
1617 :    
1618 :     $code->{"AAA"} = "K";
1619 :     $code->{"AAC"} = "N";
1620 :     $code->{"AAG"} = "K";
1621 :     $code->{"AAT"} = "N";
1622 :     $code->{"ACA"} = "T";
1623 :     $code->{"ACC"} = "T";
1624 :     $code->{"ACG"} = "T";
1625 :     $code->{"ACT"} = "T";
1626 :     $code->{"AGA"} = "R";
1627 :     $code->{"AGC"} = "S";
1628 :     $code->{"AGG"} = "R";
1629 :     $code->{"AGT"} = "S";
1630 :     $code->{"ATA"} = "I";
1631 :     $code->{"ATC"} = "I";
1632 :     $code->{"ATG"} = "M";
1633 :     $code->{"ATT"} = "I";
1634 :     $code->{"CAA"} = "Q";
1635 :     $code->{"CAC"} = "H";
1636 :     $code->{"CAG"} = "Q";
1637 :     $code->{"CAT"} = "H";
1638 :     $code->{"CCA"} = "P";
1639 :     $code->{"CCC"} = "P";
1640 :     $code->{"CCG"} = "P";
1641 :     $code->{"CCT"} = "P";
1642 :     $code->{"CGA"} = "R";
1643 :     $code->{"CGC"} = "R";
1644 :     $code->{"CGG"} = "R";
1645 :     $code->{"CGT"} = "R";
1646 :     $code->{"CTA"} = "L";
1647 :     $code->{"CTC"} = "L";
1648 :     $code->{"CTG"} = "L";
1649 :     $code->{"CTT"} = "L";
1650 :     $code->{"GAA"} = "E";
1651 :     $code->{"GAC"} = "D";
1652 :     $code->{"GAG"} = "E";
1653 :     $code->{"GAT"} = "D";
1654 :     $code->{"GCA"} = "A";
1655 :     $code->{"GCC"} = "A";
1656 :     $code->{"GCG"} = "A";
1657 :     $code->{"GCT"} = "A";
1658 :     $code->{"GGA"} = "G";
1659 :     $code->{"GGC"} = "G";
1660 :     $code->{"GGG"} = "G";
1661 :     $code->{"GGT"} = "G";
1662 :     $code->{"GTA"} = "V";
1663 :     $code->{"GTC"} = "V";
1664 :     $code->{"GTG"} = "V";
1665 :     $code->{"GTT"} = "V";
1666 :     $code->{"TAA"} = "*";
1667 :     $code->{"TAC"} = "Y";
1668 :     $code->{"TAG"} = "*";
1669 :     $code->{"TAT"} = "Y";
1670 :     $code->{"TCA"} = "S";
1671 :     $code->{"TCC"} = "S";
1672 :     $code->{"TCG"} = "S";
1673 :     $code->{"TCT"} = "S";
1674 :     $code->{"TGA"} = "*";
1675 :     $code->{"TGC"} = "C";
1676 :     $code->{"TGG"} = "W";
1677 :     $code->{"TGT"} = "C";
1678 :     $code->{"TTA"} = "L";
1679 :     $code->{"TTC"} = "F";
1680 :     $code->{"TTG"} = "L";
1681 :     $code->{"TTT"} = "F";
1682 : parrello 1.200
1683 : efrank 1.1 return $code;
1684 :     }
1685 :    
1686 :     =pod
1687 :    
1688 :     =head1 translate
1689 :    
1690 :     usage: $aa_seq = &FIG::translate($dna_seq,$code,$fix_start);
1691 :    
1692 :     If $code is undefined, I use the standard genetic code. If $fix_start is true, I
1693 :     will translate initial TTG or GTG to 'M'.
1694 :    
1695 :     =cut
1696 : parrello 1.213 #: Return Type $;
1697 : efrank 1.1 sub translate {
1698 : olson 1.111 shift if UNIVERSAL::isa($_[0],__PACKAGE__);
1699 : efrank 1.1 my( $dna,$code,$start) = @_;
1700 :     my( $i,$j,$ln );
1701 :     my( $x,$y );
1702 :     my( $prot );
1703 :    
1704 :     if (! defined($code))
1705 :     {
1706 :     $code = &FIG::standard_genetic_code;
1707 :     }
1708 :     $ln = length($dna);
1709 :     $prot = "X" x ($ln/3);
1710 :     $dna =~ tr/a-z/A-Z/;
1711 :    
1712 :     for ($i=0,$j=0; ($i < ($ln-2)); $i += 3,$j++)
1713 :     {
1714 :     $x = substr($dna,$i,3);
1715 :     if ($y = $code->{$x})
1716 :     {
1717 :     substr($prot,$j,1) = $y;
1718 :     }
1719 :     }
1720 : parrello 1.200
1721 : efrank 1.1 if (($start) && ($ln >= 3) && (substr($dna,0,3) =~ /^[GT]TG$/))
1722 :     {
1723 :     substr($prot,0,1) = 'M';
1724 :     }
1725 :     return $prot;
1726 :     }
1727 :    
1728 :     =pod
1729 :    
1730 : parrello 1.213 =head1 reverse_comp
1731 : efrank 1.1
1732 :     usage: $dnaR = &FIG::reverse_comp($dna) or
1733 :     $dnaRP = &FIG::rev_comp($seqP)
1734 :    
1735 :     In WIT, we implemented reverse complement passing a pointer to a sequence and returning
1736 :     a pointer to a sequence. In most cases the pointers are a pain (although in a few they
1737 :     are just what is needed). Hence, I kept both versions of the function to allow you
1738 :     to use whichever you like. Use rev_comp only for long strings where passing pointers is a
1739 :     reasonable effeciency issue.
1740 :    
1741 :     =cut
1742 : parrello 1.213 #: Return Type $;
1743 : efrank 1.1 sub reverse_comp {
1744 : olson 1.111 shift if UNIVERSAL::isa($_[0],__PACKAGE__);
1745 : efrank 1.1 my($seq) = @_;
1746 :    
1747 :     return ${&rev_comp(\$seq)};
1748 :     }
1749 :    
1750 : parrello 1.213 =head1 rev_comp
1751 :    
1752 :     =cut
1753 :     #: Return Type $;
1754 : efrank 1.1 sub rev_comp {
1755 : olson 1.111 shift if UNIVERSAL::isa($_[0],__PACKAGE__);
1756 : efrank 1.1 my( $seqP ) = @_;
1757 :     my( $rev );
1758 :    
1759 :     $rev = reverse( $$seqP );
1760 :     $rev =~ tr/a-z/A-Z/;
1761 :     $rev =~ tr/ACGTUMRWSYKBDHV/TGCAAKYWSRMVHDB/;
1762 :     return \$rev;
1763 :     }
1764 :    
1765 :     =pod
1766 :    
1767 :     =head1 verify_dir
1768 :    
1769 :     usage: &FIG::verify_dir($dir)
1770 :    
1771 :     Makes sure that $dir exists. If it has to create it, it sets permissions to 0777.
1772 :    
1773 :     =cut
1774 :    
1775 :     sub verify_dir {
1776 : olson 1.111 shift if UNIVERSAL::isa($_[0],__PACKAGE__);
1777 : efrank 1.1 my($dir) = @_;
1778 :    
1779 :     if (-d $dir) { return }
1780 :     if ($dir =~ /^(.*)\/[^\/]+$/)
1781 :     {
1782 :     &verify_dir($1);
1783 :     }
1784 : olson 1.153 mkdir($dir,0777) || die "could not make $dir: $!";
1785 : olson 1.184 # chmod 02777,$dir;
1786 : efrank 1.1 }
1787 :    
1788 :     =pod
1789 :    
1790 : overbeek 1.283 =head1 table_exists
1791 :    
1792 :     usage: if ($fig->table_exists($table)) { }
1793 :    
1794 :     Test for existence of table in the relational DB. This returns 0 if the table
1795 :     does not exist; 1 if it does
1796 :    
1797 :     =cut
1798 :    
1799 :    
1800 :     =pod
1801 :    
1802 : efrank 1.1 =head1 run
1803 :    
1804 :     usage: &FIG::run($cmd)
1805 :    
1806 :     Runs $cmd and fails (with trace) if the command fails.
1807 :    
1808 :     =cut
1809 :    
1810 :     sub run {
1811 : olson 1.111 shift if UNIVERSAL::isa($_[0],__PACKAGE__);
1812 : efrank 1.1 my($cmd) = @_;
1813 :    
1814 : overbeek 1.278 if ($ENV{VERBOSE})
1815 :     {
1816 :     my @tmp = `date`; chomp @tmp; print STDERR "$tmp[0]: running $cmd\n";
1817 :     }
1818 : efrank 1.1 (system($cmd) == 0) || confess "FAILED: $cmd";
1819 :     }
1820 :    
1821 : gdpusch 1.45
1822 :    
1823 :     =pod
1824 :    
1825 :     =head1 read_fasta_record(\*FILEHANDLE)
1826 :    
1827 : gdpusch 1.109 Usage: ( $seq_id, $seq_pointer, $comment ) = &read_fasta_record(\*FILEHANDLE);
1828 : gdpusch 1.45
1829 :     Function: Reads a FASTA-formatted sequence file one record at a time.
1830 :     The input filehandle defaults to STDIN if not specified.
1831 : parrello 1.200 Returns a sequence ID, a pointer to the sequence, and an optional
1832 :     record comment (NOTE: Record comments are deprecated, as some tools
1833 :     such as BLAST do not handle them gracefully). Returns an empty list
1834 :     if attempting to read a record results in an undefined value
1835 : gdpusch 1.45 (e.g., due to reaching the EOF).
1836 :    
1837 :     Author: Gordon D. Pusch
1838 :    
1839 :     Date: 2004-Feb-18
1840 :    
1841 :     =cut
1842 : parrello 1.213 #: Return Type @;
1843 : gdpusch 1.45 sub read_fasta_record
1844 :     {
1845 : olson 1.111 shift if UNIVERSAL::isa($_[0],__PACKAGE__);
1846 : gdpusch 1.45 my ($file_handle) = @_;
1847 : gdpusch 1.46 my ( $old_end_of_record, $fasta_record, @lines, $head, $sequence, $seq_id, $comment, @parsed_fasta_record );
1848 : parrello 1.200
1849 : gdpusch 1.45 if (not defined($file_handle)) { $file_handle = \*STDIN; }
1850 : parrello 1.200
1851 : gdpusch 1.45 $old_end_of_record = $/;
1852 :     $/ = "\n>";
1853 : parrello 1.200
1854 : gdpusch 1.45 if (defined($fasta_record = <$file_handle>))
1855 :     {
1856 :     chomp $fasta_record;
1857 :     @lines = split( /\n/, $fasta_record );
1858 :     $head = shift @lines;
1859 :     $head =~ s/^>?//;
1860 :     $head =~ m/^(\S+)/;
1861 :     $seq_id = $1;
1862 : parrello 1.200
1863 : gdpusch 1.45 if ($head =~ m/^\S+\s+(.*)$/) { $comment = $1; } else { $comment = ""; }
1864 : parrello 1.200
1865 : gdpusch 1.45 $sequence = join( "", @lines );
1866 : parrello 1.200
1867 : gdpusch 1.45 @parsed_fasta_record = ( $seq_id, \$sequence, $comment );
1868 :     }
1869 :     else
1870 :     {
1871 :     @parsed_fasta_record = ();
1872 :     }
1873 : parrello 1.200
1874 : gdpusch 1.45 $/ = $old_end_of_record;
1875 : parrello 1.200
1876 : gdpusch 1.45 return @parsed_fasta_record;
1877 :     }
1878 :    
1879 :    
1880 : efrank 1.1 =pod
1881 :    
1882 :     =head1 display_id_and_seq
1883 :    
1884 :     usage: &FIG::display_id_and_seq($id_and_comment,$seqP,$fh)
1885 :    
1886 :     This command has always been used to put out fasta sequences. Note that it
1887 :     takes a pointer to the sequence. $fh is optional and defalts to STDOUT.
1888 :    
1889 :     =cut
1890 :    
1891 : mkubal 1.53
1892 : efrank 1.1 sub display_id_and_seq {
1893 : olson 1.111 shift if UNIVERSAL::isa($_[0],__PACKAGE__);
1894 : efrank 1.1 my( $id, $seq, $fh ) = @_;
1895 : parrello 1.200
1896 : efrank 1.1 if (! defined($fh) ) { $fh = \*STDOUT; }
1897 : parrello 1.200
1898 : efrank 1.1 print $fh ">$id\n";
1899 :     &display_seq($seq, $fh);
1900 :     }
1901 :    
1902 :     sub display_seq {
1903 : olson 1.111 shift if UNIVERSAL::isa($_[0],__PACKAGE__);
1904 : efrank 1.1 my ( $seq, $fh ) = @_;
1905 :     my ( $i, $n, $ln );
1906 : parrello 1.200
1907 : efrank 1.1 if (! defined($fh) ) { $fh = \*STDOUT; }
1908 :    
1909 :     $n = length($$seq);
1910 :     # confess "zero-length sequence ???" if ( (! defined($n)) || ($n == 0) );
1911 :     for ($i=0; ($i < $n); $i += 60)
1912 :     {
1913 :     if (($i + 60) <= $n)
1914 :     {
1915 :     $ln = substr($$seq,$i,60);
1916 :     }
1917 :     else
1918 :     {
1919 :     $ln = substr($$seq,$i,($n-$i));
1920 :     }
1921 :     print $fh "$ln\n";
1922 :     }
1923 :     }
1924 :    
1925 :     ########## I commented the pods on the following routines out, since they should not
1926 :     ########## be part of the SOAP/WSTL interface
1927 :     #=pod
1928 :     #
1929 :     #=head1 file2N
1930 :     #
1931 :     #usage: $n = $fig->file2N($file)
1932 :     #
1933 :     #In some of the databases I need to store filenames, which can waste a lot of
1934 :     #space. Hence, I maintain a database for converting filenames to/from integers.
1935 :     #
1936 :     #=cut
1937 :     #
1938 : olson 1.111 sub file2N :scalar {
1939 : efrank 1.1 my($self,$file) = @_;
1940 :     my($relational_db_response);
1941 :    
1942 :     my $rdbH = $self->db_handle;
1943 :    
1944 :     if (($relational_db_response = $rdbH->SQL("SELECT fileno FROM file_table WHERE ( file = \'$file\')")) &&
1945 :     (@$relational_db_response == 1))
1946 :     {
1947 :     return $relational_db_response->[0]->[0];
1948 :     }
1949 :     elsif (($relational_db_response = $rdbH->SQL("SELECT MAX(fileno) FROM file_table ")) && (@$relational_db_response == 1) && ($relational_db_response->[0]->[0]))
1950 :     {
1951 :     my $fileno = $relational_db_response->[0]->[0] + 1;
1952 :     if ($rdbH->SQL("INSERT INTO file_table ( file, fileno ) VALUES ( \'$file\', $fileno )"))
1953 :     {
1954 :     return $fileno;
1955 :     }
1956 :     }
1957 :     elsif ($rdbH->SQL("INSERT INTO file_table ( file, fileno ) VALUES ( \'$file\', 1 )"))
1958 :     {
1959 :     return 1;
1960 :     }
1961 :     return undef;
1962 :     }
1963 :    
1964 :     #=pod
1965 :     #
1966 :     #=head1 N2file
1967 :     #
1968 :     #usage: $filename = $fig->N2file($n)
1969 :     #
1970 :     #In some of the databases I need to store filenames, which can waste a lot of
1971 :     #space. Hence, I maintain a database for converting filenames to/from integers.
1972 :     #
1973 :     #=cut
1974 :     #
1975 : olson 1.111 sub N2file :scalar {
1976 : efrank 1.1 my($self,$fileno) = @_;
1977 :     my($relational_db_response);
1978 :    
1979 :     my $rdbH = $self->db_handle;
1980 :    
1981 :     if (($relational_db_response = $rdbH->SQL("SELECT file FROM file_table WHERE ( fileno = $fileno )")) &&
1982 :     (@$relational_db_response == 1))
1983 :     {
1984 :     return $relational_db_response->[0]->[0];
1985 :     }
1986 :     return undef;
1987 :     }
1988 :    
1989 :    
1990 :     #=pod
1991 :     #
1992 :     #=head1 openF
1993 :     #
1994 :     #usage: $fig->openF($filename)
1995 :     #
1996 :     #Parts of the system rely on accessing numerous different files. The most obvious case is
1997 :     #the situation with similarities. It is important that the system be able to run in cases in
1998 :     #which an arbitrary number of files cannot be open simultaneously. This routine (with closeF) is
1999 :     #a hack to handle this. I should probably just pitch them and insist that the OS handle several
2000 :     #hundred open filehandles.
2001 :     #
2002 :     #=cut
2003 :     #
2004 :     sub openF {
2005 :     my($self,$file) = @_;
2006 :     my($fxs,$x,@fxs,$fh);
2007 :    
2008 :     $fxs = $self->cached('_openF');
2009 :     if ($x = $fxs->{$file})
2010 :     {
2011 :     $x->[1] = time();
2012 :     return $x->[0];
2013 :     }
2014 : parrello 1.200
2015 : efrank 1.1 @fxs = keys(%$fxs);
2016 :     if (defined($fh = new FileHandle "<$file"))
2017 :     {
2018 : overbeek 1.98 if (@fxs >= 50)
2019 : efrank 1.1 {
2020 :     @fxs = sort { $fxs->{$a}->[1] <=> $fxs->{$b}->[1] } @fxs;
2021 :     $x = $fxs->{$fxs[0]};
2022 :     undef $x->[0];
2023 :     delete $fxs->{$fxs[0]};
2024 :     }
2025 :     $fxs->{$file} = [$fh,time()];
2026 :     return $fh;
2027 :     }
2028 :     return undef;
2029 :     }
2030 :    
2031 :     #=pod
2032 :     #
2033 :     #=head1 closeF
2034 :     #
2035 :     #usage: $fig->closeF($filename)
2036 :     #
2037 :     #Parts of the system rely on accessing numerous different files. The most obvious case is
2038 :     #the situation with similarities. It is important that the system be able to run in cases in
2039 :     #which an arbitrary number of files cannot be open simultaneously. This routine (with openF) is
2040 :     #a hack to handle this. I should probably just pitch them and insist that the OS handle several
2041 :     #hundred open filehandles.
2042 :     #
2043 :     #=cut
2044 :     #
2045 :     sub closeF {
2046 :     my($self,$file) = @_;
2047 :     my($fxs,$x);
2048 :    
2049 : parrello 1.200 if (($fxs = $self->{_openF}) &&
2050 : efrank 1.1 ($x = $fxs->{$file}))
2051 :     {
2052 :     undef $x->[0];
2053 :     delete $fxs->{$file};
2054 :     }
2055 :     }
2056 :    
2057 :     =pod
2058 :    
2059 :     =head1 ec_name
2060 :    
2061 :     usage: $enzymatic_function = $fig->ec_name($ec)
2062 :    
2063 :     Returns enzymatic name for EC.
2064 :    
2065 :     =cut
2066 :    
2067 :     sub ec_name {
2068 :     my($self,$ec) = @_;
2069 :    
2070 :     ($ec =~ /^\d+\.\d+\.\d+\.\d+$/) || return "";
2071 :     my $rdbH = $self->db_handle;
2072 :     my $relational_db_response = $rdbH->SQL("SELECT name FROM ec_names WHERE ( ec = \'$ec\' )");
2073 :    
2074 :     return (@$relational_db_response == 1) ? $relational_db_response->[0]->[0] : "";
2075 :     return "";
2076 :     }
2077 :    
2078 :     =pod
2079 :    
2080 :     =head1 all_roles
2081 :    
2082 :     usage: @roles = $fig->all_roles
2083 :    
2084 : mkubal 1.54 Supposed to return all known roles. For now, we get all ECs with "names".
2085 : efrank 1.1
2086 :     =cut
2087 :    
2088 :     sub all_roles {
2089 :     my($self) = @_;
2090 :    
2091 :     my $rdbH = $self->db_handle;
2092 :     my $relational_db_response = $rdbH->SQL("SELECT ec,name FROM ec_names");
2093 :    
2094 :     return @$relational_db_response;
2095 :     }
2096 :    
2097 :     =pod
2098 :    
2099 :     =head1 expand_ec
2100 :    
2101 :     usage: $expanded_ec = $fig->expand_ec($ec)
2102 :    
2103 :     Expands "1.1.1.1" to "1.1.1.1 - alcohol dehydrogenase" or something like that.
2104 :    
2105 :     =cut
2106 :    
2107 :     sub expand_ec {
2108 :     my($self,$ec) = @_;
2109 :     my($name);
2110 :    
2111 :     return ($name = $self->ec_name($ec)) ? "$ec - $name" : $ec;
2112 :     }
2113 :    
2114 :    
2115 :     =pod
2116 :    
2117 :     =head1 clean_tmp
2118 :    
2119 :     usage: &FIG::clean_tmp
2120 :    
2121 :     We store temporary files in $FIG_Config::temp. There are specific classes of files
2122 :     that are created and should be saved for at least a few days. This routine can be
2123 :     invoked to clean out those that are over two days old.
2124 :    
2125 :     =cut
2126 :    
2127 :     sub clean_tmp {
2128 :    
2129 :     my($file);
2130 :     if (opendir(TMP,"$FIG_Config::temp"))
2131 :     {
2132 :     # change the pattern to pick up other files that need to be cleaned up
2133 :     my @temp = grep { $_ =~ /^(Geno|tmp)/ } readdir(TMP);
2134 :     foreach $file (@temp)
2135 :     {
2136 :     if (-M "$FIG_Config::temp/$file" > 2)
2137 :     {
2138 :     unlink("$FIG_Config::temp/$file");
2139 :     }
2140 :     }
2141 :     }
2142 :     }
2143 :    
2144 :     ################ Routines to process genomes and genome IDs ##########################
2145 :    
2146 :    
2147 :     =pod
2148 :    
2149 :     =head1 genomes
2150 :    
2151 : golsen 1.150 usage: @genome_ids = $fig->genomes( $complete, $restrictions, $domain );
2152 : efrank 1.1
2153 :     Genomes are assigned ids of the form X.Y where X is the taxonomic id maintained by
2154 :     NCBI for the species (not the specific strain), and Y is a sequence digit assigned to
2155 :     this particular genome (as one of a set with the same genus/species). Genomes also
2156 :     have versions, but that is a separate issue.
2157 :    
2158 :     =cut
2159 :    
2160 : olson 1.111 sub genomes :remote :list {
2161 : golsen 1.150 my( $self, $complete, $restrictions, $domain ) = @_;
2162 : overbeek 1.13
2163 :     my $rdbH = $self->db_handle;
2164 :    
2165 :     my @where = ();
2166 :     if ($complete)
2167 :     {
2168 :     push(@where,"( complete = \'1\' )")
2169 :     }
2170 :    
2171 :     if ($restrictions)
2172 :     {
2173 :     push(@where,"( restrictions = \'1\' )")
2174 :     }
2175 : golsen 1.150
2176 :     if ($domain)
2177 :     {
2178 :     push( @where, "( maindomain = '$domain' )" )
2179 :     }
2180 :    
2181 : overbeek 1.13 my $relational_db_response;
2182 :     if (@where > 0)
2183 :     {
2184 :     my $where = join(" AND ",@where);
2185 :     $relational_db_response = $rdbH->SQL("SELECT genome FROM genome where $where");
2186 :     }
2187 :     else
2188 :     {
2189 :     $relational_db_response = $rdbH->SQL("SELECT genome FROM genome");
2190 :     }
2191 :     my @genomes = sort { $a <=> $b } map { $_->[0] } @$relational_db_response;
2192 : efrank 1.1 return @genomes;
2193 :     }
2194 :    
2195 : overbeek 1.180 sub is_complete {
2196 :     my($self,$genome) = @_;
2197 :    
2198 :     my $rdbH = $self->db_handle;
2199 :     my $relational_db_response = $rdbH->SQL("SELECT genome FROM genome where (genome = '$genome') AND (complete = '1')");
2200 :     return (@$relational_db_response == 1)
2201 :     }
2202 : golsen 1.150
2203 : efrank 1.2 sub genome_counts {
2204 : overbeek 1.13 my($self,$complete) = @_;
2205 :     my($x,$relational_db_response);
2206 : efrank 1.2
2207 : overbeek 1.13 my $rdbH = $self->db_handle;
2208 :    
2209 : parrello 1.200 if ($complete)
2210 : overbeek 1.13 {
2211 : gdpusch 1.107 $relational_db_response = $rdbH->SQL("SELECT genome,maindomain FROM genome where complete = '1'");
2212 : overbeek 1.13 }
2213 :     else
2214 :     {
2215 : gdpusch 1.107 $relational_db_response = $rdbH->SQL("SELECT genome,maindomain FROM genome");
2216 : overbeek 1.13 }
2217 :    
2218 : gdpusch 1.107 my ($arch, $bact, $euk, $vir, $env, $unk) = (0, 0, 0, 0, 0, 0);
2219 : overbeek 1.13 if (@$relational_db_response > 0)
2220 : efrank 1.2 {
2221 : overbeek 1.13 foreach $x (@$relational_db_response)
2222 : efrank 1.2 {
2223 : gdpusch 1.107 if ($x->[1] =~ /^archaea/i) { ++$arch }
2224 :     elsif ($x->[1] =~ /^bacter/i) { ++$bact }
2225 :     elsif ($x->[1] =~ /^eukar/i) { ++$euk }
2226 :     elsif ($x->[1] =~ /^vir/i) { ++$vir }
2227 :     elsif ($x->[1] =~ /^env/i) { ++$env }
2228 :     else { ++$unk }
2229 : efrank 1.2 }
2230 :     }
2231 : parrello 1.200
2232 : gdpusch 1.107 return ($arch, $bact, $euk, $vir, $env, $unk);
2233 :     }
2234 :    
2235 :    
2236 :     =pod
2237 :    
2238 :     =head1 genome_domain
2239 :    
2240 :     usage: $domain = $fig->genome_domain($genome_id);
2241 :    
2242 :     Returns the domain of a genome ID, and 'undef' if it is not in the database.
2243 :    
2244 :     =cut
2245 :    
2246 :     sub genome_domain {
2247 :     my($self,$genome) = @_;
2248 :     my $relational_db_response;
2249 :     my $rdbH = $self->db_handle;
2250 : parrello 1.200
2251 :     if ($genome)
2252 : gdpusch 1.107 {
2253 :     if (($relational_db_response = $rdbH->SQL("SELECT genome,maindomain FROM genome WHERE ( genome = \'$genome\' )"))
2254 :     && (@$relational_db_response == 1))
2255 :     {
2256 :     # die Dumper($relational_db_response);
2257 :     return $relational_db_response->[0]->[1];
2258 :     }
2259 :     }
2260 :     return undef;
2261 : efrank 1.2 }
2262 :    
2263 : gdpusch 1.92
2264 :     =pod
2265 :    
2266 :     =head1 genome_pegs
2267 :    
2268 : gdpusch 1.107 usage: $num_pegs = $fig->genome_pegs($genome_id);
2269 : gdpusch 1.92
2270 : gdpusch 1.107 Returns the number of protein-encoding genes (PEGs) in $genome_id if
2271 :     "$genome_id" is indexed in the "genome" database, and 'undef' otherwise.
2272 : gdpusch 1.92
2273 :     =cut
2274 :    
2275 :     sub genome_pegs {
2276 :     my($self,$genome) = @_;
2277 :     my $relational_db_response;
2278 :     my $rdbH = $self->db_handle;
2279 : parrello 1.200
2280 :     if ($genome)
2281 : gdpusch 1.92 {
2282 :     if (($relational_db_response = $rdbH->SQL("SELECT pegs FROM genome WHERE ( genome = \'$genome\' )"))
2283 :     && (@$relational_db_response == 1))
2284 :     {
2285 :     return $relational_db_response->[0]->[0];
2286 :     }
2287 :     }
2288 :     return undef;
2289 :     }
2290 :    
2291 :    
2292 : efrank 1.1 =pod
2293 :    
2294 : gdpusch 1.92 =head1 genome_rnas
2295 :    
2296 : gdpusch 1.107 usage: $num_rnas = $fig->genome_rnas($genome_id);
2297 : gdpusch 1.92
2298 : gdpusch 1.107 Returns the number of RNA-encoding genes (RNAs) in $genome_id if
2299 :     "$genome_id" is indexed in the "genome" database, and 'undef' otherwise.
2300 : gdpusch 1.92
2301 :     =cut
2302 :    
2303 :     sub genome_rnas {
2304 :     my($self,$genome) = @_;
2305 :     my $relational_db_response;
2306 :     my $rdbH = $self->db_handle;
2307 : parrello 1.200
2308 :     if ($genome)
2309 : gdpusch 1.92 {
2310 :     if (($relational_db_response = $rdbH->SQL("SELECT rnas FROM genome WHERE ( genome = \'$genome\' )"))
2311 :     && (@$relational_db_response == 1))
2312 :     {
2313 :     return $relational_db_response->[0]->[0];
2314 :     }
2315 :     }
2316 :     return undef;
2317 :     }
2318 :    
2319 :    
2320 :     =pod
2321 :    
2322 :     =head1 genome_szdna
2323 : efrank 1.1
2324 : gdpusch 1.92 usage: $szdna = $fig->genome_szdna($genome_id);
2325 : gdpusch 1.91
2326 : gdpusch 1.107 Returns the number of DNA base-pairs in the genome contigs file(s) of $genome_id
2327 :     "$genome_id" is indexed in the "genome" database, and 'undef' otherwise.
2328 : gdpusch 1.91
2329 :     =cut
2330 :    
2331 : gdpusch 1.92 sub genome_szdna {
2332 : gdpusch 1.91 my($self,$genome) = @_;
2333 :     my $relational_db_response;
2334 :     my $rdbH = $self->db_handle;
2335 : parrello 1.200
2336 :     if ($genome)
2337 : gdpusch 1.91 {
2338 :     if (($relational_db_response = $rdbH->SQL("SELECT szdna FROM genome WHERE ( genome = \'$genome\' )"))
2339 :     && (@$relational_db_response == 1))
2340 :     {
2341 :     return $relational_db_response->[0]->[0];
2342 :     }
2343 :     }
2344 :     return undef;
2345 :     }
2346 :    
2347 :    
2348 :     =pod
2349 :    
2350 :     =head1 genome_version
2351 :    
2352 : efrank 1.1 usage: $version = $fig->genome_version($genome_id);
2353 :    
2354 :     Versions are incremented for major updates. They are put in as major
2355 :     updates of the form 1.0, 2.0, ...
2356 :    
2357 :     Users may do local "editing" of the DNA for a genome, but when they do,
2358 :     they increment the digits to the right of the decimal. Two genomes remain
2359 : parrello 1.200 comparable only if the versions match identically. Hence, minor updating should be
2360 : efrank 1.1 committed only by the person/group responsible for updating that genome.
2361 :    
2362 :     We can, of course, identify which genes are identical between any two genomes (by matching
2363 :     the DNA or amino acid sequences). However, the basic intent of the system is to
2364 :     support editing by the main group issuing periodic major updates.
2365 :    
2366 :     =cut
2367 :    
2368 : olson 1.113 sub genome_version :scalar {
2369 : efrank 1.1 my($self,$genome) = @_;
2370 :    
2371 :     my(@tmp);
2372 :     if ((-s "$FIG_Config::organisms/$genome/VERSION") &&
2373 :     (@tmp = `cat $FIG_Config::organisms/$genome/VERSION`) &&
2374 : overbeek 1.84 ($tmp[0] =~ /^(\S+)$/))
2375 : efrank 1.1 {
2376 :     return $1;
2377 :     }
2378 :     return undef;
2379 :     }
2380 :    
2381 : olson 1.236 =head1 genome_md5sum
2382 :    
2383 :     usage: $md5sum = $fig->genome_md5sum($genome_id);
2384 :    
2385 :     Returns the MD5 checksum of genome $genome_id.
2386 :    
2387 :     The checksum of a genome is defined as the checksum of its signature file. The signature
2388 :     file consists of tab-separated lines, one for each contig, ordered by the contig id.
2389 :     Each line contains the contig ID, the length of the contig in nucleotides, and the
2390 :     MD5 checksum of the nucleotide data, with uppercase letters forced to lower case.
2391 :    
2392 :    
2393 :     =cut
2394 :    
2395 : olson 1.237 sub genome_md5sum :scalar {
2396 : olson 1.236 my($self,$genome) = @_;
2397 :     my $relational_db_response;
2398 :     my $rdbH = $self->db_handle;
2399 :    
2400 :     if ($genome)
2401 :     {
2402 :     if (($relational_db_response = $rdbH->SQL("SELECT md5sum FROM genome_md5sum WHERE ( genome = \'$genome\' )"))
2403 :     && (@$relational_db_response == 1))
2404 :     {
2405 :     return $relational_db_response->[0]->[0];
2406 :     }
2407 :     }
2408 :     return undef;
2409 :     }
2410 :    
2411 : olson 1.260 =head1 genome_with_md5sum
2412 :    
2413 :     usage: $genome = $fig->genome_with_md5sum($cksum);
2414 :    
2415 :     Find a genome with checksum $cksum.
2416 :    
2417 :     =cut
2418 :    
2419 :     sub genome_with_md5sum :scalar {
2420 :     my($self,$cksum) = @_;
2421 :     my $relational_db_response;
2422 :     my $rdbH = $self->db_handle;
2423 :    
2424 :     if (($relational_db_response = $rdbH->SQL("SELECT genome FROM genome_md5sum WHERE ( md5sum = \'$cksum\' )"))
2425 :     && (@$relational_db_response == 1))
2426 :     {
2427 :     return $relational_db_response->[0]->[0];
2428 :     }
2429 :    
2430 :     return undef;
2431 :     }
2432 :    
2433 : olson 1.237 sub contig_md5sum :scalar {
2434 :     my($self, $genome, $contig) = @_;
2435 :     my $relational_db_response;
2436 :     my $rdbH = $self->db_handle;
2437 :    
2438 :     if ($genome)
2439 :     {
2440 :     if (($relational_db_response = $rdbH->SQL(qq(SELECT md5 FROM contig_md5sums
2441 :     WHERE (genome = ? AND
2442 :     contig = ?)), undef, $genome, $contig))
2443 :     && (@$relational_db_response == 1))
2444 :     {
2445 :     return $relational_db_response->[0]->[0];
2446 :     }
2447 :     }
2448 :     return undef;
2449 :     }
2450 :    
2451 : efrank 1.1 =pod
2452 :    
2453 :     =head1 genus_species
2454 :    
2455 :     usage: $gs = $fig->genus_species($genome_id)
2456 :    
2457 : parrello 1.200 Returns the genus and species (and strain if that has been properly recorded)
2458 : efrank 1.1 in a printable form.
2459 :    
2460 :     =cut
2461 :    
2462 : olson 1.111 sub genus_species :scalar {
2463 : efrank 1.1 my ($self,$genome) = @_;
2464 : overbeek 1.13 my $ans;
2465 : efrank 1.1
2466 :     my $genus_species = $self->cached('_genus_species');
2467 :     if (! ($ans = $genus_species->{$genome}))
2468 :     {
2469 : overbeek 1.13 my $rdbH = $self->db_handle;
2470 :     my $relational_db_response = $rdbH->SQL("SELECT genome,gname FROM genome");
2471 :     my $pair;
2472 :     foreach $pair (@$relational_db_response)
2473 : efrank 1.1 {
2474 : overbeek 1.13 $genus_species->{$pair->[0]} = $pair->[1];
2475 : efrank 1.1 }
2476 : overbeek 1.13 $ans = $genus_species->{$genome};
2477 : efrank 1.1 }
2478 :     return $ans;
2479 :     }
2480 :    
2481 :     =pod
2482 :    
2483 :     =head1 org_of
2484 :    
2485 :     usage: $org = $fig->org_of($prot_id)
2486 :    
2487 :     In the case of external proteins, we can usually determine an organism, but not
2488 : parrello 1.200 anything more precise than genus/species (and often not that). This routine takes
2489 : efrank 1.2 a protein ID (which may be a feature ID) and returns "the organism".
2490 : efrank 1.1
2491 :     =cut
2492 :    
2493 :     sub org_of {
2494 :     my($self,$prot_id) = @_;
2495 :     my $relational_db_response;
2496 :     my $rdbH = $self->db_handle;
2497 :    
2498 :     if ($prot_id =~ /^fig\|/)
2499 :     {
2500 : golsen 1.138 #
2501 :     # Trying to guess what Ross wanted (there was a servere bug):
2502 :     #
2503 :     # deleted -> undefined
2504 :     # failed lookup -> ""
2505 :     #
2506 :     return $self->is_deleted_fid( $prot_id) ? undef
2507 :     : $self->genus_species( $self->genome_of( $prot_id ) ) || "";
2508 : efrank 1.1 }
2509 :    
2510 :     if (($relational_db_response = $rdbH->SQL("SELECT org FROM external_orgs WHERE ( prot = \'$prot_id\' )")) &&
2511 :     (@$relational_db_response >= 1))
2512 :     {
2513 : overbeek 1.248 $relational_db_response->[0]->[0] =~ s/^\d+://;
2514 : efrank 1.1 return $relational_db_response->[0]->[0];
2515 :     }
2516 :     return "";
2517 :     }
2518 :    
2519 : golsen 1.130 #
2520 :     # Support for colorizing organisms by domain
2521 :     # -- GJO
2522 :     #
2523 :     =pod
2524 :    
2525 :     =head1 genus_species_domain
2526 :    
2527 :     usage: ($gs, $domain) = $fig->genus_species_domain($genome_id)
2528 :    
2529 : parrello 1.200 Returns the genus and species (and strain if that has been properly recorded)
2530 : golsen 1.130 in a printable form, and domain.
2531 :    
2532 :     =cut
2533 :    
2534 :     sub genus_species_domain {
2535 :     my ($self, $genome) = @_;
2536 :    
2537 :     my $genus_species_domain = $self->cached('_genus_species_domain');
2538 :     if ( ! $genus_species_domain->{ $genome } )
2539 :     {
2540 :     my $rdbH = $self->db_handle;
2541 :     my $relational_db_response = $rdbH->SQL("SELECT genome,gname,maindomain FROM genome");
2542 :     my $triple;
2543 :     foreach $triple ( @$relational_db_response )
2544 :     {
2545 :     $genus_species_domain->{ $triple->[0] } = [ $triple->[1], $triple->[2] ];
2546 :     }
2547 :     }
2548 :     my $gsdref = $genus_species_domain->{ $genome };
2549 :     return $gsdref ? @$gsdref : ( "", "" );
2550 :     }
2551 :    
2552 :    
2553 :     my %domain_color = ( AR => "#DDFFFF", BA => "#FFDDFF", EU => "#FFFFDD",
2554 :     VI => "#DDDDDD", EN => "#BBBBBB" );
2555 :    
2556 :     sub domain_color {
2557 :     my ($domain) = @_;
2558 :     defined $domain || return "#FFFFFF";
2559 :     return $domain_color{ uc substr($domain, 0, 2) } || "#FFFFFF";
2560 :     }
2561 :    
2562 :    
2563 :     =pod
2564 :    
2565 :     =head1 org_and_color_of
2566 :    
2567 :     usage: ($org, $color) = $fig->org_and_domain_of($prot_id)
2568 :    
2569 :     Return the best guess organism and domain html color string of an organism.
2570 :     In the case of external proteins, we can usually determine an organism, but not
2571 : parrello 1.200 anything more precise than genus/species (and often not that). This routine takes
2572 : golsen 1.130 a protein ID (which may be a feature ID) and returns "the organism".
2573 :    
2574 :     =cut
2575 :    
2576 :     sub org_and_color_of {
2577 :     my($self,$prot_id) = @_;
2578 :     my $relational_db_response;
2579 :     my $rdbH = $self->db_handle;
2580 :    
2581 :     if ($prot_id =~ /^fig\|/)
2582 :     {
2583 :     my( $gs, $domain ) = $self->genus_species_domain($self->genome_of($prot_id));
2584 :     return ( $gs, domain_color( $domain ) );
2585 :     }
2586 :    
2587 :     if (($relational_db_response = $rdbH->SQL("SELECT org FROM external_orgs WHERE ( prot = \'$prot_id\' )")) &&
2588 :     (@$relational_db_response >= 1))
2589 :     {
2590 :     return ($relational_db_response->[0]->[0], "#FFFFFF");
2591 :     }
2592 :     return ("", "#FFFFFF");
2593 :     }
2594 :    
2595 :     #
2596 :     # End of support for colorizing organisms by domain
2597 :     # -- GJO
2598 :     #
2599 :    
2600 : efrank 1.1 =pod
2601 :    
2602 :     =head1 abbrev
2603 :    
2604 :     usage: $abbreviated_name = $fig->abbrev($genome_name)
2605 :    
2606 :     For alignments and such, it is very useful to be able to produce an abbreviation of genus/species.
2607 :     That's what this does. Note that multiple genus/species might reduce to the same abbreviation, so
2608 :     be careful (disambiguate them, if you must).
2609 :    
2610 :     =cut
2611 :    
2612 : olson 1.111 sub abbrev :scalar {
2613 :     shift if UNIVERSAL::isa($_[0],__PACKAGE__);
2614 : efrank 1.1 my($genome_name) = @_;
2615 :    
2616 :     $genome_name =~ s/^(\S{3})\S+/$1./;
2617 : overbeek 1.198 $genome_name =~ s/^(\S+)\s+(\S{3})\S+/$1$2./;
2618 : overbeek 1.257 $genome_name =~ s/ //g;
2619 :     if (length($genome_name) > 10)
2620 : efrank 1.1 {
2621 : overbeek 1.257 $genome_name = substr($genome_name,0,10);
2622 : efrank 1.1 }
2623 :     return $genome_name;
2624 :     }
2625 :    
2626 :     ################ Routines to process Features and Feature IDs ##########################
2627 :    
2628 :     =pod
2629 :    
2630 :     =head1 ftype
2631 :    
2632 :     usage: $type = &FIG::ftype($fid)
2633 :    
2634 :     Returns the type of a feature, given the feature ID. This just amounts
2635 :     to lifting it out of the feature ID, since features have IDs of tghe form
2636 :    
2637 :     fig|x.y.f.n
2638 :    
2639 :     where
2640 :     x.y is the genome ID
2641 :     f is the type pf feature
2642 :     n is an integer that is unique within the genome/type
2643 :    
2644 :     =cut
2645 :    
2646 :     sub ftype {
2647 : olson 1.111 shift if UNIVERSAL::isa($_[0],__PACKAGE__);
2648 : efrank 1.1 my($feature_id) = @_;
2649 :    
2650 :     if ($feature_id =~ /^fig\|\d+\.\d+\.([^\.]+)/)
2651 :     {
2652 :     return $1;
2653 :     }
2654 :     return undef;
2655 :     }
2656 :    
2657 :     =pod
2658 :    
2659 :     =head1 genome_of
2660 :    
2661 :     usage: $genome_id = $fig->genome_of($fid)
2662 :    
2663 :     This just extracts the genome ID from a feature ID.
2664 :    
2665 :     =cut
2666 :    
2667 :    
2668 : olson 1.113 sub genome_of :scalar {
2669 : olson 1.111 shift if UNIVERSAL::isa($_[0],__PACKAGE__);
2670 : parrello 1.200 my $prot_id = (@_ == 1) ? $_[0] : $_[1];
2671 : efrank 1.1
2672 :     if ($prot_id =~ /^fig\|(\d+\.\d+)/) { return $1; }
2673 :     return undef;
2674 :     }
2675 :    
2676 : olson 1.96 =head1 genome_and_peg_of
2677 :    
2678 :     usage: ($genome_id, $peg_number = $fig->genome_and_peg_of($fid)
2679 :    
2680 :     This just extracts the genome ID and peg number from a feature ID.
2681 :    
2682 :     =cut
2683 :    
2684 :    
2685 :     sub genome_and_peg_of {
2686 : olson 1.111 shift if UNIVERSAL::isa($_[0],__PACKAGE__);
2687 : parrello 1.200 my $prot_id = (@_ == 1) ? $_[0] : $_[1];
2688 : olson 1.96
2689 :     if ($prot_id =~ /^fig\|(\d+\.\d+)\.peg\.(\d+)/)
2690 :     {
2691 :     return ($1, $2);
2692 :     }
2693 :     return undef;
2694 :     }
2695 :    
2696 : efrank 1.1 =pod
2697 :    
2698 :     =head1 by_fig_id
2699 :    
2700 :     usage: @sorted_by_fig_id = sort { &FIG::by_fig_id($a,$b) } @fig_ids
2701 :    
2702 :     This is a bit of a clutzy way to sort a list of FIG feature IDs, but it works.
2703 :    
2704 :     =cut
2705 :    
2706 :     sub by_fig_id {
2707 :     my($a,$b) = @_;
2708 :     my($g1,$g2,$t1,$t2,$n1,$n2);
2709 :     if (($a =~ /^fig\|(\d+\.\d+).([^\.]+)\.(\d+)$/) && (($g1,$t1,$n1) = ($1,$2,$3)) &&
2710 :     ($b =~ /^fig\|(\d+\.\d+).([^\.]+)\.(\d+)$/) && (($g2,$t2,$n2) = ($1,$2,$3)))
2711 :     {
2712 :     ($g1 <=> $g2) or ($t1 cmp $t2) or ($n1 <=> $n2);
2713 :     }
2714 :     else
2715 :     {
2716 :     $a cmp $b;
2717 :     }
2718 :     }
2719 :    
2720 :     =pod
2721 :    
2722 :     =head1 genes_in_region
2723 :    
2724 :     usage: ($features_in_region,$beg1,$end1) = $fig->genes_in_region($genome,$contig,$beg,$end)
2725 :    
2726 :     It is often important to be able to find the genes that occur in a specific region on
2727 :     a chromosome. This routine is designed to provide this information. It returns all genes
2728 :     that overlap the region ($genome,$contig,$beg,$end). $beg1 is set to the minimum coordinate of
2729 :     the returned genes (which may be before the given region), and $end1 the maximum coordinate.
2730 :    
2731 :     The routine assumes that genes are not more than 10000 bases long, which is certainly not true
2732 :     in eukaryotes. Hence, in euks you may well miss genes that overlap the boundaries of the specified
2733 :     region (sorry).
2734 :    
2735 :     =cut
2736 : parrello 1.213 #: Return Type @;
2737 : efrank 1.1 sub genes_in_region {
2738 :     my($self,$genome,$contig,$beg,$end) = @_;
2739 :     my($x,$relational_db_response,$feature_id,$b1,$e1,@feat,@tmp,$l,$u);
2740 :    
2741 :     my $pad = 10000;
2742 :     my $rdbH = $self->db_handle;
2743 :    
2744 : parrello 1.200 my $minV = $beg - $pad;
2745 : efrank 1.1 my $maxV = $end + $pad;
2746 : parrello 1.200 if (($relational_db_response = $rdbH->SQL("SELECT id FROM features
2747 : golsen 1.141 WHERE ( minloc > $minV ) AND ( minloc < $maxV ) AND ( maxloc < $maxV) AND
2748 : efrank 1.1 ( genome = \'$genome\' ) AND ( contig = \'$contig\' );")) &&
2749 :     (@$relational_db_response >= 1))
2750 :     {
2751 : parrello 1.200 @tmp = sort { ($a->[1] cmp $b->[1]) or
2752 : overbeek 1.129 (($a->[2]+$a->[3]) <=> ($b->[2]+$b->[3]))
2753 : efrank 1.1 }
2754 : parrello 1.200 map { $feature_id = $_->[0];
2755 :     $x = $self->feature_location($feature_id);
2756 :     $x ? [$feature_id,&boundaries_of($x)] : ()
2757 : efrank 1.1 } @$relational_db_response;
2758 :    
2759 :    
2760 :     ($l,$u) = (10000000000,0);
2761 :     foreach $x (@tmp)
2762 :     {
2763 :     ($feature_id,undef,$b1,$e1) = @$x;
2764 :     if (&between($beg,&min($b1,$e1),$end) || &between(&min($b1,$e1),$beg,&max($b1,$e1)))
2765 :     {
2766 : overbeek 1.136 if (! $self->is_deleted_fid($feature_id))
2767 :     {
2768 :     push(@feat,$feature_id);
2769 :     $l = &min($l,&min($b1,$e1));
2770 :     $u = &max($u,&max($b1,$e1));
2771 :     }
2772 : efrank 1.1 }
2773 :     }
2774 :     (@feat <= 0) || return ([@feat],$l,$u);
2775 :     }
2776 :     return ([],$l,$u);
2777 :     }
2778 :    
2779 : golsen 1.141
2780 :     #=============================================================================
2781 :     # Using the following version is better, but it brings out a very annoying
2782 :     # issue with some genomes. It already exists in the current code (above)
2783 :     # for some genes in some genomes. For example, visit fig|70601.1.peg.1.
2784 :     # This is true for any genome that has a feature that crosses the origin.
2785 :     # The root of the problem lies in boundaries_of. I am working on a fix that
2786 :     # replaces boundaries_of with a more sophisticated function. When it is
2787 :     # all done, genes_in_retion should behave as desired. -- GJO, Aug. 22, 2004
2788 :     #=============================================================================
2789 : parrello 1.200 #
2790 : golsen 1.141 # =pod
2791 : parrello 1.200 #
2792 : golsen 1.141 # =head1 genes_in_region
2793 : parrello 1.200 #
2794 : golsen 1.141 # usage: ( $features_in_region, $min_coord, $max_coord )
2795 :     # = $fig->genes_in_region( $genome, $contig, $beg, $end )
2796 : parrello 1.200 #
2797 : golsen 1.141 # It is often important to be able to find the genes that occur in a specific
2798 :     # region on a chromosome. This routine is designed to provide this information.
2799 :     # It returns all genes that overlap the region ( $genome, $contig, $beg, $end ).
2800 :     # $min_coord is set to the minimum coordinate of the returned genes (which may
2801 :     # preceed the given region), and $max_coord the maximum coordinate. Because
2802 :     # the database is indexed by the leftmost and rightmost coordinates of each
2803 :     # feature, the function makes no assumption about the length of the feature, but
2804 :     # it can (and probably will) miss features spanning multiple contigs.
2805 : parrello 1.200 #
2806 : golsen 1.141 # =cut
2807 : parrello 1.200 #
2808 :     #
2809 : golsen 1.141 # sub genes_in_region {
2810 :     # my ( $self, $genome, $contig, $beg, $end ) = @_;
2811 :     # my ( $x, $db_response, $feature_id, $b1, $e1, @tmp, @bounds );
2812 :     # my ( $min_coord, $max_coord );
2813 : parrello 1.200 #
2814 : golsen 1.141 # my @features = ();
2815 :     # my $rdbH = $self->db_handle;
2816 : parrello 1.200 #
2817 : golsen 1.141 # if ( ( $db_response = $rdbH->SQL( "SELECT id
2818 :     # FROM features
2819 :     # WHERE ( contig = '$contig' )
2820 :     # AND ( genome = '$genome' )
2821 : parrello 1.200 # AND ( minloc <= $end )
2822 : golsen 1.141 # AND ( maxloc >= $beg );"
2823 :     # )
2824 :     # )
2825 :     # && ( @$db_response > 0 )
2826 :     # )
2827 :     # {
2828 :     # # The sort is unnecessary, but provides a consistent ordering
2829 : parrello 1.200 #
2830 : golsen 1.141 # @tmp = sort { ( $a->[1] cmp $b->[1] ) # contig
2831 :     # || ( ($a->[2] + $a->[3] ) <=> ( $b->[2] + $b->[3] ) ) # midpoint
2832 :     # }
2833 :     # map { $feature_id = $_->[0];
2834 :     # ( ( ! $self->is_deleted_fid( $feature_id ) ) # not deleted
2835 :     # && ( $x = $self->feature_location( $feature_id ) ) # and has location
2836 :     # && ( ( @bounds = boundaries_of( $x ) ) == 3 ) # and has bounds
2837 : parrello 1.200 # ) ? [ $feature_id, @bounds ] : ()
2838 : golsen 1.141 # } @$db_response;
2839 : parrello 1.200 #
2840 : golsen 1.141 # ( $min_coord, $max_coord ) = ( 10000000000, 0 );
2841 : parrello 1.200 #
2842 : golsen 1.141 # foreach $x ( @tmp )
2843 :     # {
2844 :     # ( $feature_id, undef, $b1, $e1 ) = @$x;
2845 :     # push @features, $feature_id;
2846 :     # my ( $min, $max ) = ( $b1 <= $e1 ) ? ( $b1, $e1 ) : ( $e1, $b1 );
2847 :     # ( $min_coord <= $min ) || ( $min_coord = $min );
2848 :     # ( $max_coord >= $max ) || ( $max_coord = $max );
2849 :     # }
2850 :     # }
2851 : parrello 1.200 #
2852 : golsen 1.141 # return ( @features ) ? ( [ @features ], $min_coord, $max_coord )
2853 :     # : ( [], undef, undef );
2854 :     # }
2855 :    
2856 :     # These will be part of the fix to genes_in_region. -- GJO
2857 :    
2858 :     =pod
2859 :    
2860 :     =head1 regions_spanned
2861 :    
2862 :     usage: ( [ $contig, $beg, $end ], ... ) = $fig->regions_spanned( $loc )
2863 :    
2864 :     The location of a feature in a scalar context is
2865 :    
2866 :     contig_b1_e1,contig_b2_e2,... [one contig_b_e for each segment]
2867 :    
2868 :     This routine takes as input a fig location and reduces it to one or more
2869 :     regions spanned by the gene. Unlike boundaries_of, regions_spanned handles
2870 :     wrapping through the orgin, features split over contigs and exons that are
2871 :     not ordered nicely along the chromosome (ugly but true).
2872 :    
2873 :     =cut
2874 :    
2875 :     sub regions_spanned {
2876 :     shift if UNIVERSAL::isa( $_[0], __PACKAGE__ );
2877 :     my( $location ) = ( @_ == 1 ) ? $_[0] : $_[1];
2878 :     defined( $location ) || return undef;
2879 :    
2880 :     my @regions = ();
2881 :    
2882 :     my ( $cur_contig, $cur_beg, $cur_end, $cur_dir );
2883 :     my ( $contig, $beg, $end, $dir );
2884 :     my @segs = split( /\s*,\s*/, $location ); # should not have space, but ...
2885 :     @segs || return undef;
2886 :    
2887 :     # Process the first segment
2888 :    
2889 :     my $seg = shift @segs;
2890 :     ( ( $cur_contig, $cur_beg, $cur_end ) = ( $seg =~ /^(\S+)_(\d+)_\d+$/ ) )
2891 :     || return undef;
2892 :     $cur_dir = ( $cur_end >= $cur_beg ) ? 1 : -1;
2893 :    
2894 :     foreach $seg ( @segs ) {
2895 :     ( ( $contig, $beg, $end ) = ( $seg =~ /^(\S+)_(\d+)_\d+$/ ) ) || next;
2896 :     $dir = ( $end >= $beg ) ? 1 : -1;
2897 :    
2898 :     # Is this a continuation? Update end
2899 :    
2900 :     if ( ( $contig eq $cur_contig )
2901 :     && ( $dir == $cur_dir )
2902 :     && ( ( ( $dir > 0 ) && ( $end > $cur_end ) )
2903 :     || ( ( $dir < 0 ) && ( $end < $cur_end ) ) )
2904 :     )
2905 :     {
2906 :     $cur_end = $end;
2907 :     }
2908 :    
2909 :     # Not a continuation. Report previous and update current.
2910 :    
2911 :     else
2912 :     {
2913 :     push @regions, [ $cur_contig, $cur_beg, $cur_end ];
2914 :     ( $cur_contig, $cur_beg, $cur_end, $cur_dir )
2915 :     = ( $contig, $beg, $end, $dir );
2916 :     }
2917 :     }
2918 :    
2919 :     # There should alwasy be a valid, unreported region.
2920 :    
2921 :     push @regions, [ $cur_contig, $cur_beg, $cur_end ];
2922 :    
2923 :     return wantarray ? @regions : \@regions;
2924 :     }
2925 :    
2926 :    
2927 :     =pod
2928 :    
2929 :     =head1 filter_regions
2930 :    
2931 :     usage: @regions = filter_regions( $contig, $min, $max, @regions )
2932 :     \@regions = filter_regions( $contig, $min, $max, @regions )
2933 :     @regions = filter_regions( $contig, $min, $max, \@regions )
2934 :     \@regions = filter_regions( $contig, $min, $max, \@regions )
2935 :    
2936 :     This function provides a simple filter for extracting a list of genome regions
2937 :     for those that overlap a particular interval. Region definitions correspond
2938 :     to those produced by regions_spanned. That is, [ contig, beg, end ]. In the
2939 :     function call, either $contig or $min and $max can be undefined (permitting
2940 :     anything).
2941 :    
2942 :     =cut
2943 :    
2944 :     sub filter_regions {
2945 :     my ( $contig, $min, $max, @regions ) = @_;
2946 :    
2947 :     @regions || return ();
2948 :     ( ref( $regions[0] ) eq "ARRAY" ) || return undef;
2949 :    
2950 :     # Is it a region list, or a reference to a region list?
2951 :    
2952 :     if ( ref( $regions[0]->[0] ) eq "ARRAY" ) { @regions = @{ $regions[0] } }
2953 :    
2954 :     if ( ! defined( $contig ) )
2955 :     {
2956 :     ( defined( $min ) && defined( $max ) ) || return undef;
2957 :     }
2958 :     else # with a defined contig name, allow undefined range
2959 :     {
2960 :     defined( $min ) || ( $min = 1 );
2961 :     defined( $max ) || ( $max = 1000000000 );
2962 :     }
2963 :     ( $min <= $max ) || return ();
2964 :    
2965 :     my ( $c, $b, $e );
2966 :     my @filtered = grep { ( @$_ >= 3 ) # Allow extra fields?
2967 :     && ( ( $c, $b, $e ) = @$_ )
2968 :     && ( ( ! defined( $contig ) ) || ( $c eq $contig ) )
2969 :     && ( ( $e >= $b ) || ( ( $b, $e ) = ( $e, $b ) ) )
2970 :     && ( ( $b <= $max ) && ( $e >= $min ) )
2971 :     } @regions;
2972 :    
2973 :     return wantarray ? @filtered : \@filtered;
2974 :     }
2975 :    
2976 :    
2977 : efrank 1.1 sub close_genes {
2978 :     my($self,$fid,$dist) = @_;
2979 : parrello 1.200
2980 : mkubal 1.147 # warn "In close_genes, self=$self, fid=$fid";
2981 : parrello 1.200
2982 : efrank 1.1 my $loc = $self->feature_location($fid);
2983 :     if ($loc)
2984 :     {
2985 :     my($contig,$beg,$end) = &FIG::boundaries_of($loc);
2986 :     if ($contig && $beg && $end)
2987 :     {
2988 :     my $min = &min($beg,$end) - $dist;
2989 :     my $max = &max($beg,$end) + $dist;
2990 :     my $feat;
2991 :     ($feat,undef,undef) = $self->genes_in_region(&FIG::genome_of($fid),$contig,$min,$max);
2992 :     return @$feat;
2993 :     }
2994 :     }
2995 :     return ();
2996 :     }
2997 :    
2998 : mkubal 1.147 sub adjacent_genes
2999 :     {
3000 :     my ($self, $fid, $dist) = @_;
3001 :     my (@close, $strand, $i);
3002 : parrello 1.200
3003 : mkubal 1.147 # warn "In adjacent_genes, self=$self, fid=$fid";
3004 : parrello 1.200
3005 :    
3006 : mkubal 1.147 $strand = $self->strand_of($fid);
3007 : parrello 1.200
3008 : mkubal 1.147 $dist = $dist || 2000;
3009 :     @close = $self->close_genes($fid, $dist);
3010 :     for ($i=0; $i < @close; ++$i) { last if ($close[$i] eq $fid); }
3011 : parrello 1.200
3012 : redwards 1.157 # RAE note that if $self->strand_of($close[$i-1]) ne $strand then left/right neighbors
3013 :     # were never set! oops!
3014 : parrello 1.200
3015 : redwards 1.157 # I think the concept of Left and right is confused here. In my mind, left and right
3016 :     # are independent of strand ?? E.g. take a look at PEG fig|196600.1.peg.1806
3017 :     # this is something like
3018 :     #
3019 :     # ---> <--1805--- --1806--> <--1807-- <----
3020 :     #
3021 :     # 1805 is always the left neighbor, no?
3022 :    
3023 :     my ($left_neighbor, $right_neighbor) = ($close[$i-1], $close[$i+1]);
3024 : parrello 1.200
3025 : redwards 1.157 if (0) # this was if ($i > 0) I just skip this whole section!
3026 : mkubal 1.147 {
3027 :     if ($self->strand_of($close[$i-1]) eq $strand) { $left_neighbor = $close[$i-1]; }
3028 : redwards 1.157 # else {$left_neighbor=$close[$i+1]} # RAE: this is the alternative that is needed if you do it by strand
3029 : mkubal 1.147 }
3030 : parrello 1.200
3031 : mkubal 1.147 if ($i < $#close)
3032 :     {
3033 :     if ($self->strand_of($close[$i+1]) eq $strand) { $right_neighbor = $close[$i+1]; }
3034 : redwards 1.157 # else {$right_neighbor = $close[$i-1]} # RAE: this is the alternative that is needed if you do it by strand
3035 : mkubal 1.147 }
3036 : parrello 1.200
3037 : mkubal 1.147 # ...return genes in transcription order...
3038 : parrello 1.200 if ($strand eq '-')
3039 : mkubal 1.147 {
3040 :     ($left_neighbor, $right_neighbor) = ($right_neighbor, $left_neighbor);
3041 :     }
3042 : parrello 1.200
3043 : mkubal 1.147 return ($left_neighbor, $right_neighbor) ;
3044 :     }
3045 :    
3046 : efrank 1.1
3047 :     =pod
3048 :    
3049 :     =head1 feature_location
3050 :    
3051 :     usage: $loc = $fig->feature_location($fid) OR
3052 :     @loc = $fig->feature_location($fid)
3053 :    
3054 :     The location of a feature in a scalar context is
3055 :    
3056 :     contig_b1_e1,contig_b2_e2,... [one contig_b_e for each exon]
3057 :    
3058 :     In a list context it is
3059 :    
3060 :     (contig_b1_e1,contig_b2_e2,...)
3061 :    
3062 :     =cut
3063 :    
3064 : olson 1.111 sub feature_location :scalar :list {
3065 : efrank 1.1 my($self,$feature_id) = @_;
3066 :     my($relational_db_response,$locations,$location);
3067 : parrello 1.200
3068 : mkubal 1.147 # warn "In feature_location, self=$self, fid=$feature_id";
3069 : parrello 1.200
3070 : overbeek 1.136 if ($self->is_deleted_fid($feature_id)) { return undef }
3071 :    
3072 : efrank 1.1 $locations = $self->cached('_location');
3073 :     if (! ($location = $locations->{$feature_id}))
3074 :     {
3075 :     my $rdbH = $self->db_handle;
3076 :     if (($relational_db_response = $rdbH->SQL("SELECT location FROM features WHERE ( id = \'$feature_id\' )")) &&
3077 :     (@$relational_db_response == 1))
3078 :     {
3079 :     $locations->{$feature_id} = $location = $relational_db_response->[0]->[0];
3080 :     }
3081 :     }
3082 :    
3083 :     if ($location)
3084 :     {
3085 :     return wantarray() ? split(/,/,$location) : $location;
3086 :     }
3087 :     return undef;
3088 :     }
3089 :    
3090 : mkubal 1.147 sub contig_of
3091 :     {
3092 :     my ($self, $locus) = @_;
3093 : parrello 1.200
3094 : olson 1.159 $locus =~ m/^([^,]+)_\d+_\d+/;
3095 : parrello 1.200
3096 : mkubal 1.147 return $1;
3097 :     }
3098 :    
3099 :     sub beg_of
3100 :     {
3101 :     my ($self, $locus) = @_;
3102 : parrello 1.200
3103 : olson 1.159 $locus =~ m/^[^,]+_(\d+)_\d+/;
3104 : parrello 1.200
3105 : mkubal 1.147 return $1;
3106 :     }
3107 :    
3108 :     sub end_of
3109 :     {
3110 :     my ($self, $locus) = @_;
3111 : parrello 1.200
3112 : mkubal 1.147 $locus =~ m/\S+_\d+_(\d+)$/;
3113 : parrello 1.200
3114 : mkubal 1.147 return $1;
3115 :     }
3116 :    
3117 : parrello 1.200 sub strand_of
3118 : mkubal 1.147 {
3119 :     my ($self, $fid) = @_;
3120 :     my ($beg, $end);
3121 : parrello 1.200
3122 : mkubal 1.147 # warn "In strand_of, self=$self, fid=$fid";
3123 : parrello 1.200
3124 : mkubal 1.147 (undef, $beg, $end) = $self->boundaries_of($self->feature_location($fid));
3125 : parrello 1.200
3126 : mkubal 1.147 if ($beg < $end) { return '+'; } else { return '-'; }
3127 :     }
3128 :    
3129 : olson 1.158 =pod
3130 :    
3131 :     =head1 find_contig_with_checksum
3132 :    
3133 :     Find a contig in the given genome with the given checksum.
3134 :    
3135 :    
3136 :     =cut
3137 :    
3138 :     sub find_contig_with_checksum
3139 :     {
3140 :     my($self, $genome, $checksum) = @_;
3141 : parrello 1.200
3142 : olson 1.158 #
3143 :     # This implementation scans all the contig files for the organism; when
3144 :     # we have contig checksums in the database this will simplify
3145 :     # significantly.
3146 :     #
3147 :     # For some efficiency, we cache the checksums we compute along the way since
3148 :     # it's probably likely we'll poke at other contigs for this organism.
3149 :     #
3150 :    
3151 :     my $gdir = "$FIG_Config::organisms/$genome";
3152 :    
3153 :     my $cached_cksum = $self->cached('_contig_checksum');
3154 : parrello 1.200
3155 : olson 1.158 if (opendir(my $dh, $gdir))
3156 :     {
3157 :     for my $file (map { "$gdir/$_" } grep { $_ =~ /^contigs\d*$/ } readdir($dh))
3158 :     {
3159 :     local $/ = "\n>";
3160 : parrello 1.200
3161 : olson 1.158 if (open(my $fh, "<$file"))
3162 :     {
3163 :     while (<$fh>)
3164 :     {
3165 :     chomp;
3166 : olson 1.160
3167 :     #
3168 :     # Pull the contig identifier from the first line.
3169 :     # We need the >? to handle the first line in the file;
3170 :     # the others are removed by the chomp above because
3171 :     # $/ is set to "\n>";
3172 :     #
3173 : parrello 1.200
3174 : olson 1.160 if (s/^>?\s*(\S+)([^\n]*)\n//)
3175 : olson 1.158 {
3176 :     my $ident = $1;
3177 :     my $contig_txt = $_;
3178 : parrello 1.200
3179 : olson 1.158 $contig_txt =~ s/\s//sg;
3180 :     $contig_txt = uc($contig_txt);
3181 :    
3182 :     #
3183 :     # See if we have a cached value.
3184 :     #
3185 :    
3186 :     my $this_checksum;
3187 :     $this_checksum = $cached_cksum->{$genome, $ident};
3188 :     if (!$this_checksum)
3189 :     {
3190 : parrello 1.200
3191 : olson 1.158 my($rd, $wr, $pid);
3192 : parrello 1.200
3193 : olson 1.158 if (!($pid = open2($rd, $wr, "cksum")))
3194 :     {
3195 :     die "Cannot run open2 cksum: $!";
3196 :     }
3197 : parrello 1.200
3198 : olson 1.158 $wr->write($contig_txt, length($contig_txt));
3199 : parrello 1.200
3200 : olson 1.158 close($wr);
3201 : parrello 1.200
3202 : olson 1.158 $_ = <$rd>;
3203 :     close($rd);
3204 :     waitpid $pid, 0;
3205 :     chomp;
3206 : parrello 1.200
3207 : olson 1.158 my @vals = split(/\s+/, $_);
3208 :     $this_checksum = $vals[0];
3209 :     $cached_cksum->{$genome, $ident} = $this_checksum;
3210 :     }
3211 :     if ($this_checksum == $checksum)
3212 :     {
3213 :     return $ident;
3214 :     }
3215 :     }
3216 :     }
3217 :     }
3218 :     }
3219 :     }
3220 :     }
3221 :    
3222 :     sub contig_checksum
3223 :     {
3224 :     my($self, $genome, $contig) = @_;
3225 :    
3226 :     my $contig_txt = $self->read_contig($genome, $contig);
3227 :    
3228 :     $contig_txt =~ s/\s//sg;
3229 :     $contig_txt = uc($contig_txt);
3230 :    
3231 : olson 1.159 my($rd, $wr, $pid);
3232 : olson 1.158
3233 : olson 1.159 if (!($pid = open2($rd, $wr, "cksum")))
3234 : olson 1.158 {
3235 :     die "Cannot run open2 cksum: $!";
3236 :     }
3237 :    
3238 :     $wr->write($contig_txt, length($contig_txt));
3239 : parrello 1.200
3240 : olson 1.158 close($wr);
3241 :    
3242 :     $_ = <$rd>;
3243 : olson 1.159 close($rd);
3244 :     waitpid $pid, 0;
3245 :    
3246 : olson 1.158 chomp;
3247 :     my @vals = split(/\s+/, $_);
3248 :     if (wantarray)
3249 :     {
3250 :     return @vals;
3251 :     }
3252 :     else
3253 :     {
3254 :     return $vals[0];
3255 :     }
3256 :     }
3257 :    
3258 :     =pod
3259 :    
3260 :     =head1 read_contig
3261 :    
3262 :     Read a single contig from the contigs file.
3263 :    
3264 :     =cut
3265 :     sub read_contig
3266 :     {
3267 :     my($self, $genome, $contig) = @_;
3268 :    
3269 :     #
3270 :     # Read the contig. The database has it in a set of chunks, but we
3271 :     # just use the seek to the starting point, and read up to the next "\n>".
3272 :     #
3273 :    
3274 :     my $ret = $self->db_handle->SQL(qq!SELECT fileno, seek FROM contig_seeks
3275 :     WHERE genome = '$genome' and contig = '$contig' and
3276 :     startn = 0!);
3277 :     if (!$ret or @$ret != 1)
3278 :     {
3279 :     return undef;
3280 :     }
3281 :    
3282 :     my($fileno, $seek) = @{$ret->[0]};
3283 :     my $file = $self->N2file($fileno);
3284 :    
3285 :     my($fh, $contig_txt);
3286 : parrello 1.200
3287 : olson 1.158 if (!open($fh, "<$file"))
3288 :     {
3289 :     warn "contig_checksum: could not open $file: $!\n";
3290 :     return undef;
3291 :     }
3292 :    
3293 :     seek($fh, $seek, 0);
3294 :    
3295 :     {
3296 :     local $/ = "\n>";
3297 :    
3298 :     $contig_txt = <$fh>;
3299 :     chomp($contig_txt);
3300 :     }
3301 :    
3302 :     return $contig_txt;
3303 :     }
3304 : mkubal 1.147
3305 : efrank 1.1 =pod
3306 :    
3307 :     =head1 boundaries_of
3308 :    
3309 :     usage: ($contig,$beg,$end) = $fig->boundaries_of($loc)
3310 :    
3311 :     The location of a feature in a scalar context is
3312 :    
3313 :     contig_b1_e1,contig_b2_e2,... [one contig_b_e for each exon]
3314 :    
3315 :     This routine takes as input such a location and reduces it to a single
3316 :     description of the entire region containing the gene.
3317 :    
3318 :     =cut
3319 :    
3320 :     sub boundaries_of {
3321 : olson 1.111 shift if UNIVERSAL::isa($_[0],__PACKAGE__);
3322 : efrank 1.1 my($location) = (@_ == 1) ? $_[0] : $_[1];
3323 :     my($contigQ);
3324 :    
3325 :     if (defined($location))
3326 :     {
3327 :     my @exons = split(/,/,$location);
3328 :     my($contig,$beg,$end);
3329 : olson 1.203
3330 : parrello 1.200 if (($exons[0] =~ /^(\S+)_(\d+)_\d+$/) &&
3331 : efrank 1.1 (($contig,$beg) = ($1,$2)) && ($contigQ = quotemeta $contig) &&
3332 : parrello 1.200 ($exons[$#exons] =~ /^$contigQ\_\d+_(\d+)$/) &&
3333 : efrank 1.1 ($end = $1))
3334 :     {
3335 :     return ($contig,$beg,$end);
3336 :     }
3337 :     }
3338 :     return undef;
3339 :     }
3340 :    
3341 : heiko 1.175 =pod
3342 :    
3343 :     =head1 all_features_detailed
3344 :    
3345 :     usage: $fig->all_features_detailed($genome)
3346 :    
3347 :     Returns a list of all feature IDs, location, aliases, type in the designated genome.
3348 :     Used in gendb import to speed up the process
3349 :     =cut
3350 :    
3351 :     sub all_features_detailed {
3352 :     my($self,$genome) = @_;
3353 :    
3354 :     my $rdbH = $self->db_handle;
3355 :     my $relational_db_response = $rdbH->SQL("SELECT id, location, aliases, type FROM features WHERE (genome = \'$genome\')");
3356 :     my @features;
3357 :     foreach my $tuple (@$relational_db_response)
3358 :     {
3359 :     push @features, $tuple unless ($self->is_deleted_fid($tuple->[0]));
3360 :     }
3361 :     return \@features;
3362 :     }
3363 :    
3364 : efrank 1.1
3365 :     =pod
3366 :    
3367 :     =head1 all_features
3368 :    
3369 :     usage: $fig->all_features($genome,$type)
3370 :    
3371 :     Returns a list of all feature IDs of a specified type in the designated genome. You would
3372 : parrello 1.200 usually use just
3373 : efrank 1.1
3374 : parrello 1.200 $fig->pegs_of($genome) or
3375 : efrank 1.1 $fig->rnas_of($genome)
3376 :    
3377 :     which simply invoke this routine.
3378 :    
3379 :     =cut
3380 :    
3381 :     sub all_features {
3382 :     my($self,$genome,$type) = @_;
3383 :    
3384 :     my $rdbH = $self->db_handle;
3385 :     my $relational_db_response = $rdbH->SQL("SELECT id FROM features WHERE (genome = \'$genome\' AND (type = \'$type\'))");
3386 :    
3387 :     if (@$relational_db_response > 0)
3388 :     {
3389 : overbeek 1.136 return grep { ! $self->is_deleted_fid($_) } map { $_->[0] } @$relational_db_response;
3390 : efrank 1.1 }
3391 :     return ();
3392 :     }
3393 :    
3394 :    
3395 :     =pod
3396 :    
3397 : overbeek 1.152 =head1 pegs_of
3398 : efrank 1.1
3399 : overbeek 1.152 usage: $fig->pegs_of($genome)
3400 : efrank 1.1
3401 :     Returns a list of all PEGs in the specified genome. Note that order is not
3402 :     specified.
3403 :    
3404 :     =cut
3405 :    
3406 :     sub pegs_of {
3407 :     my($self,$genome) = @_;
3408 : parrello 1.200
3409 : efrank 1.1 return $self->all_features($genome,"peg");
3410 :     }
3411 :    
3412 :    
3413 :     =pod
3414 :    
3415 : overbeek 1.152 =head1 rnas_of
3416 : efrank 1.1
3417 : overbeek 1.152 usage: $fig->rnas_of($genome)
3418 : efrank 1.1
3419 :     Returns a list of all RNAs for the given genome.
3420 :    
3421 :     =cut
3422 :    
3423 :     sub rnas_of {
3424 :     my($self,$genome) = @_;
3425 : parrello 1.200
3426 : efrank 1.1 return $self->all_features($genome,"rna");
3427 :     }
3428 :    
3429 :     =pod
3430 :    
3431 :     =head1 feature_aliases
3432 :    
3433 :     usage: @aliases = $fig->feature_aliases($fid) OR
3434 : parrello 1.200 $aliases = $fig->feature_aliases($fid)
3435 : efrank 1.1
3436 :     Returns a list of aliases (gene IDs, arbitrary numbers assigned by authors, etc.) for the feature.
3437 :     These must come from the tbl files, so add them there if you want to see them here.
3438 :    
3439 :     In a scalar context, the aliases come back with commas separating them.
3440 :    
3441 :     =cut
3442 :    
3443 :     sub feature_aliases {
3444 :     my($self,$feature_id) = @_;
3445 : overbeek 1.87 my($rdbH,$relational_db_response,@aliases,$aliases,%aliases,$x);
3446 : efrank 1.1
3447 : overbeek 1.136 if ($self->is_deleted_fid($feature_id)) { return undef }
3448 :    
3449 : efrank 1.1 $rdbH = $self->db_handle;
3450 : overbeek 1.87 @aliases = ();
3451 : efrank 1.1 if (($relational_db_response = $rdbH->SQL("SELECT aliases FROM features WHERE ( id = \'$feature_id\' )")) &&
3452 :     (@$relational_db_response == 1))
3453 :     {
3454 :     $aliases = $relational_db_response->[0]->[0];
3455 : overbeek 1.87 %aliases = map { $_ => 1 } split(/,/,$aliases);
3456 : overbeek 1.173 }
3457 :    
3458 :     if (($relational_db_response = $rdbH->SQL("SELECT alias FROM ext_alias WHERE ( id = \'$feature_id\' )")) &&
3459 :     (@$relational_db_response > 0))
3460 :     {
3461 :     foreach $x (@$relational_db_response)
3462 : overbeek 1.87 {
3463 : overbeek 1.173 $aliases{$x->[0]} = 1;
3464 : overbeek 1.87 }
3465 : efrank 1.1 }
3466 : parrello 1.200
3467 : overbeek 1.173 @aliases = sort keys(%aliases);
3468 : overbeek 1.87
3469 : overbeek 1.131 return wantarray() ? @aliases : join(",",@aliases);
3470 : efrank 1.1 }
3471 :    
3472 :     =pod
3473 :    
3474 : overbeek 1.34 =head1 by_alias
3475 :    
3476 : parrello 1.200 usage: $peg = $fig->by_alias($alias)
3477 : overbeek 1.34
3478 :     Returns a FIG id if the alias can be converted. Right now we convert aliases
3479 : overbeek 1.211 of the form NP_* (RefSeq IDs), gi|* (GenBank IDs), sp|* (Swiss Prot), uni|* (UniProt),
3480 :     kegg|* (KEGG) and maybe a few more
3481 : overbeek 1.34
3482 :     =cut
3483 :    
3484 :     sub by_alias {
3485 : overbeek 1.148 my($self,$alias,$genome) = @_;
3486 : overbeek 1.34 my($rdbH,$relational_db_response,$peg);
3487 : parrello 1.200 my ($peg, $flag) = FIGRules::NormalizeAlias($alias);
3488 :     if ($flag) {
3489 :     return $peg;
3490 :     } else {
3491 : parrello 1.201 my $genomeQ = $genome ? quotemeta $genome : "";
3492 :     $rdbH = $self->db_handle;
3493 : parrello 1.200
3494 : parrello 1.201 if (($relational_db_response = $rdbH->SQL("SELECT id FROM ext_alias WHERE ( alias = ? )", undef, $peg)) &&
3495 :     (@$relational_db_response > 0)) {
3496 :    
3497 : overbeek 1.209 if (@$relational_db_response == 1) {
3498 :     $peg = $relational_db_response->[0]->[0];
3499 :     return wantarray() ? ($peg) : $peg;
3500 :     } elsif (wantarray()) {
3501 :     return map { $_->[0] } @$relational_db_response;
3502 :     }
3503 :     }
3504 :     return wantarray() ? () : "";
3505 : overbeek 1.148 }
3506 : overbeek 1.34 }
3507 :    
3508 : overbeek 1.170 sub to_alias {
3509 :     my($self,$fid,$type) = @_;
3510 : parrello 1.200
3511 : overbeek 1.170 my @aliases = grep { $_ =~ /^$type\|/ } $self->feature_aliases($fid);
3512 :    
3513 : overbeek 1.171 if (wantarray())
3514 :     {
3515 :     return @aliases;
3516 :     }
3517 :     elsif (@aliases > 0)
3518 : overbeek 1.170 {
3519 :     return $aliases[0];
3520 :     }
3521 :     else
3522 :     {
3523 :     return "";
3524 :     }
3525 :     }
3526 :    
3527 : overbeek 1.34 =pod
3528 :    
3529 : efrank 1.1 =head1 possibly_truncated
3530 :    
3531 :     usage: $fig->possibly_truncated($fid)
3532 :    
3533 :     Returns true iff the feature occurs near the end of a contig.
3534 :    
3535 :     =cut
3536 :    
3537 :     sub possibly_truncated {
3538 :     my($self,$feature_id) = @_;
3539 :     my($loc);
3540 :    
3541 :     if ($loc = $self->feature_location($feature_id))
3542 :     {
3543 : parrello 1.200 my $genome = $self->genome_of($feature_id);
3544 :     my ($contig,$beg,$end) = $self->boundaries_of($loc);
3545 : efrank 1.1 if ((! $self->near_end($genome,$contig,$beg)) && (! $self->near_end($genome,$contig,$end)))
3546 :     {
3547 :     return 0;
3548 :     }
3549 :     }
3550 :     return 1;
3551 :     }
3552 :    
3553 :     sub near_end {
3554 :     my($self,$genome,$contig,$x) = @_;
3555 :    
3556 :     return (($x < 300) || ($x > ($self->contig_ln($genome,$contig) - 300)));
3557 :     }
3558 :    
3559 : overbeek 1.27 sub is_real_feature {
3560 :     my($self,$fid) = @_;
3561 :     my($relational_db_response);
3562 :    
3563 : overbeek 1.136 if ($self->is_deleted_fid($fid)) { return 0 }
3564 :    
3565 : overbeek 1.27 my $rdbH = $self->db_handle;
3566 :     return (($relational_db_response = $rdbH->SQL("SELECT id FROM features WHERE ( id = \'$fid\' )")) &&
3567 : mkubal 1.53 (@$relational_db_response == 1)) ? 1 : 0;
3568 : overbeek 1.27 }
3569 :    
3570 : efrank 1.1 ################ Routines to process functional coupling for PEGs ##########################
3571 :    
3572 :     =pod
3573 :    
3574 : overbeek 1.284 =head1 coupled_to
3575 :    
3576 :     usage: @coupled_to = $fig->coupled_to($peg)
3577 :    
3578 :     The new form of coupling and evidence computation is based on precomputed data.
3579 :     The old form took minutes to dynamically compute things when needed. The old form
3580 :     still works, ikf the directory Data/CouplingData is not present. If it is present,
3581 :     it is assumed to contain comprehensive coupling data in the form of precomputed scores
3582 :     and PCHs.
3583 :    
3584 :     If Data/CouplingData is present, this routine returns a list of 2-tuples [Peg,Sc]. It
3585 :     returns the empty list if the peg is not coupled. It returns undef, if Data/CouplingData
3586 :     is not there.
3587 :    
3588 :     =cut
3589 :    
3590 :     sub coupled_to {
3591 :     my($self,$peg) = @_;
3592 :    
3593 :     my $rdbH = $self->db_handle;
3594 :     if (! $rdbH->table_exists('fc_pegs')) { return undef }
3595 :    
3596 :     my $relational_db_response = $rdbH->SQL("SELECT peg2,score FROM fc_pegs
3597 :     WHERE peg1 = \'$peg\' ");
3598 :     return @$relational_db_response;
3599 :     }
3600 :    
3601 :     =pod
3602 :    
3603 :     =head1 coupling_evidence
3604 :    
3605 :     usage: @evidence = $fig->coupling_evidence($peg1,$peg2)
3606 :    
3607 :     The new form of coupling and evidence computation is based on precomputed data.
3608 :     The old form took minutes to dynamically compute things when needed. The old form
3609 :     still works, ikf the directory Data/CouplingData is not present. If it is present,
3610 :     it is assumed to contain comprehensive coupling data in the form of precomputed scores
3611 :     and PCHs.
3612 :    
3613 :     If Data/CouplingData is present, this routine returns a list of 3-tuples [Peg3,Peg4,Rep].
3614 :     Here, Peg3 is similar to Peg1, Peg4 is similar to Peg2, and Rep == 1 iff this is a
3615 :     "representative pair". That is, we take all pairs and create a representative set
3616 :     in which each pair is not "too close" to any other pair in the representative set.
3617 :     Think of "too close" as being roughly 95% identical at the DNA level. This keeps (usually)
3618 :     a single pair from a set of different genomes from the same species.
3619 :    
3620 :     It returns the empty list if the peg is not coupled. It returns undef, if Data/CouplingData
3621 :     is not there.
3622 :    
3623 :     =cut
3624 :    
3625 :     sub coupling_evidence {
3626 :     my($self,$peg1,$peg2) = @_;
3627 :    
3628 :     my $rdbH = $self->db_handle;
3629 :     if (! $rdbH->table_exists('pchs')) { return undef }
3630 :    
3631 :     my $relational_db_response = $rdbH->SQL("SELECT peg3,peg4,rep FROM pchs
3632 :     WHERE peg1 = \'$peg1\' AND
3633 :     peg2 = \'$peg2\'");
3634 :     return @$relational_db_response;
3635 :     }
3636 :    
3637 :     =pod
3638 :    
3639 : efrank 1.1 =head1 coupling_and_evidence
3640 :    
3641 :     usage: @coupling_data = $fig->coupling_and_evidence($fid,$bound,$sim_cutoff,$coupling_cutoff,$keep_record)
3642 :    
3643 :     A computation of couplings and evidence starts with a given peg and produces a list of
3644 :     3-tuples. Each 3-tuple is of the form
3645 :    
3646 :     [Score,CoupledToFID,Evidence]
3647 :    
3648 :     Evidence is a list of 2-tuples of FIDs that are close in other genomes (producing
3649 :     a "pair of close homologs" of [$peg,CoupledToFID]). The maximum score for a single
3650 :     PCH is 1, but "Score" is the sum of the scores for the entire set of PCHs.
3651 :    
3652 : overbeek 1.284 NOTE: once the new version of precomputed coupling is installed (i.e., when Data/CouplingData
3653 :     is filled with the precomputed relations), the parameters on computing evidence are ignored.
3654 :    
3655 : parrello 1.200 If $keep_record is true, the system records the information, asserting coupling for each
3656 : efrank 1.1 of the pairs in the set of evidence, and asserting a pin from the given $fd through all
3657 :     of the PCH entries used in forming the score.
3658 :    
3659 :     =cut
3660 :    
3661 :     sub coupling_and_evidence {
3662 : overbeek 1.284 my($self,$peg1,$bound,$sim_cutoff,$coupling_cutoff,$keep_record) = @_;
3663 :    
3664 :     my $rdbH = $self->db_handle;
3665 :     if ($rdbH->table_exists('fc_pegs') && $self->is_complete(&FIG::genome_of($peg1)))
3666 :     {
3667 :     my @ans = ();
3668 :     my $tuple;
3669 :     foreach $tuple ($self->coupled_to($peg1))
3670 :     {
3671 :     my($peg2,$sc) = @$tuple;
3672 :     my $evidence = [map { [$_->[0],$_->[1]] } $self->coupling_evidence($peg1,$peg2)];
3673 :     push(@ans,[$sc,$peg2,$evidence]);
3674 :     }
3675 :     return @ans;
3676 :     }
3677 :     else
3678 :     {
3679 :     return &coupling_and_evidence1($self,$peg1,$bound,$sim_cutoff,$coupling_cutoff,$keep_record);
3680 :     }
3681 :     }
3682 :    
3683 :     sub coupling_and_evidence1 {
3684 : efrank 1.1 my($self,$feature_id,$bound,$sim_cutoff,$coupling_cutoff,$keep_record) = @_;
3685 :     my($neighbors,$neigh,$similar1,$similar2,@hits,$sc,$ev,$genome1);
3686 :    
3687 : overbeek 1.136 if ($self->is_deleted_fid($feature_id)) { return undef }
3688 :    
3689 : efrank 1.1 if ($feature_id =~ /^fig\|(\d+\.\d+)/)
3690 :     {
3691 :     $genome1 = $1;
3692 :     }
3693 : overbeek 1.136 else
3694 :     {
3695 :     return undef;
3696 :     }
3697 : parrello 1.200 my $locations = $self->feature_location($feature_id);
3698 :     my($contig,$beg,$end) = $self->boundaries_of($locations);
3699 : efrank 1.1 if (! $contig) { return () }
3700 :    
3701 : parrello 1.200 ($neighbors,undef,undef) = $self->genes_in_region($self->genome_of($feature_id),
3702 : efrank 1.1 $contig,
3703 : parrello 1.200 &min($beg,$end) - $bound,
3704 : efrank 1.1 &max($beg,$end) + $bound);
3705 :     if (@$neighbors == 0) { return () }
3706 :     $similar1 = $self->acceptably_close($feature_id,$sim_cutoff);
3707 :     @hits = ();
3708 :    
3709 :     foreach $neigh (grep { $_ =~ /peg/ } @$neighbors)
3710 :     {
3711 :     next if ($neigh eq $feature_id);
3712 :     $similar2 = $self->acceptably_close($neigh,$sim_cutoff);
3713 :     ($sc,$ev) = $self->coupling_ev($genome1,$similar1,$similar2,$bound);
3714 :     if ($sc >= $coupling_cutoff)
3715 :     {
3716 :     push(@hits,[$sc,$neigh,$ev]);
3717 :     }
3718 :     }
3719 :     if ($keep_record)
3720 :     {
3721 :     $self->add_chr_clusters_and_pins($feature_id,\@hits);
3722 :     }
3723 :     return sort { $b->[0] <=> $a->[0] } @hits;
3724 :     }
3725 :    
3726 : overbeek 1.35 sub fast_coupling {
3727 :     my($self,$peg,$bound,$coupling_cutoff) = @_;
3728 :     my($genome,$genome1,$genome2,$peg1,$peg2,$peg3,%maps,$loc,$loc1,$loc2,$loc3);
3729 :     my($pairs,$sc,%ev);
3730 :    
3731 : overbeek 1.136 if ($self->is_deleted_fid($peg)) { return undef }
3732 :    
3733 : overbeek 1.35 my @ans = ();
3734 :    
3735 :     $genome = &genome_of($peg);
3736 :     foreach $peg1 ($self->in_pch_pin_with($peg))
3737 :     {
3738 :     $peg1 =~ s/,.*$//;
3739 :     if ($peg ne $peg1)
3740 :     {
3741 :     $genome1 = &genome_of($peg1);
3742 :     $maps{$peg}->{$genome1} = $peg1;
3743 :     }
3744 :     }
3745 :    
3746 :     $loc = [&boundaries_of(scalar $self->feature_location($peg))];
3747 :     foreach $peg1 ($self->in_cluster_with($peg))
3748 :     {
3749 :     if ($peg ne $peg1)
3750 :     {
3751 :     # print STDERR "peg1=$peg1\n";
3752 :     $loc1 = [&boundaries_of(scalar $self->feature_location($peg1))];
3753 :     if (&close_enough($loc,$loc1,$bound))
3754 :     {
3755 :     foreach $peg2 ($self->in_pch_pin_with($peg1))
3756 :     {
3757 :     $genome2 = &genome_of($peg2);
3758 :     if (($peg3 = $maps{$peg}->{$genome2}) && ($peg2 ne $peg3))
3759 :     {
3760 :     $loc2 = [&boundaries_of(scalar $self->feature_location($peg2))];
3761 :     $loc3 = [&boundaries_of(scalar $self->feature_location($peg3))];
3762 :     if (&close_enough($loc2,$loc3,$bound))
3763 :     {
3764 :     push(@{$ev{$peg1}},[$peg3,$peg2]);
3765 :     }
3766 :     }
3767 :     }
3768 :     }
3769 :     }
3770 :     }
3771 :     foreach $peg1 (keys(%ev))
3772 :     {
3773 :     $pairs = $ev{$peg1};
3774 : overbeek 1.43 $sc = $self->score([$peg,map { $_->[0] } @$pairs]);
3775 : overbeek 1.35 if ($sc >= $coupling_cutoff)
3776 :     {
3777 :     push(@ans,[$sc,$peg1]);
3778 :     }
3779 :     }
3780 :     return sort { $b->[0] <=> $a->[0] } @ans;
3781 :     }
3782 :    
3783 :    
3784 :     sub score {
3785 : overbeek 1.43 my($self,$pegs) = @_;
3786 : overbeek 1.51 my(@ids);
3787 : overbeek 1.35
3788 : overbeek 1.51 if ($self->{_no9s_scoring})
3789 :     {
3790 :     @ids = map { $self->maps_to_id($_) } grep { $_ !~ /^fig\|999999/ } @$pegs;
3791 :     }
3792 :     else
3793 :     {
3794 :     @ids = map { $self->maps_to_id($_) } @$pegs;
3795 :     }
3796 : overbeek 1.43 return &score1($self,\@ids) - 1;
3797 :     }
3798 :    
3799 :     sub score1 {
3800 :     my($self,$pegs) = @_;
3801 :     my($sim);
3802 : overbeek 1.204
3803 :     my($iden_cutoff) = 97;
3804 :     my($iden_cutoff_gap) = 100 - $iden_cutoff;
3805 :    
3806 : overbeek 1.43 my($first,@rest) = @$pegs;
3807 :     my $count = 1;
3808 :     my %hits = map { $_ => 1 } @rest;
3809 :     my @ordered = sort { $b->[0] <=> $a->[0] }
3810 :     map { $sim = $_; [$sim->iden,$sim->id2] }
3811 :     grep { $hits{$_->id2} }
3812 :     $self->sims($first,1000,1,"raw");
3813 : overbeek 1.76 my %ordered = map { $_->[1] => 1 } @ordered;
3814 :     foreach $_ (@rest)
3815 :     {
3816 :     if (! $ordered{$_})
3817 :     {
3818 :     push(@ordered,[0,$_]);
3819 :     }
3820 :     }
3821 :    
3822 : overbeek 1.204 while ((@ordered > 0) && ($ordered[0]->[0] >= $iden_cutoff))
3823 : overbeek 1.35 {
3824 : overbeek 1.43 shift @ordered ;
3825 :     }
3826 :     while (@ordered > 0)
3827 :     {
3828 :     my $start = $ordered[0]->[0];
3829 :     $_ = shift @ordered;
3830 :     my @sub = ( $_->[1] );
3831 : overbeek 1.204 while ((@ordered > 0) && ($ordered[0]->[0] > ($start-$iden_cutoff_gap)))
3832 : overbeek 1.35 {
3833 : overbeek 1.43 $_ = shift @ordered;
3834 :     push(@sub, $_->[1]);
3835 : overbeek 1.35 }
3836 :    
3837 : overbeek 1.43 if (@sub == 1)
3838 :     {
3839 :     $count++;
3840 :     }
3841 :     else
3842 :     {
3843 :     $count += &score1($self,\@sub);
3844 :     }
3845 : overbeek 1.35 }
3846 : overbeek 1.43 return $count;
3847 : overbeek 1.35 }
3848 :    
3849 : efrank 1.1 =pod
3850 :    
3851 :     =head1 add_chr_clusters_and_pins
3852 :    
3853 :     usage: $fig->add_chr_clusters_and_pins($peg,$hits)
3854 :    
3855 :     The system supports retaining data relating to functional coupling. If a user
3856 :     computes evidence once and then saves it with this routine, data relating to
3857 :     both "the pin" and the "clusters" (in all of the organisms supporting the
3858 :     functional coupling) will be saved.
3859 :    
3860 :     $hits must be a pointer to a list of 3-tuples of the sort returned by
3861 :     $fig->coupling_and_evidence.
3862 :    
3863 :     =cut
3864 :    
3865 :     sub add_chr_clusters_and_pins {
3866 :     my($self,$peg,$hits) = @_;
3867 :     my(@clusters,@pins,$x,$sc,$neigh,$pairs,$y,@corr,@orgs,%projection);
3868 :     my($genome,$cluster,$pin,$peg2);
3869 :    
3870 :     if (@$hits > 0)
3871 :     {
3872 :     @clusters = ();
3873 :     @pins = ();
3874 :     push(@clusters,[$peg,map { $_->[1] } @$hits]);
3875 :     foreach $x (@$hits)
3876 :     {
3877 :     ($sc,$neigh,$pairs) = @$x;
3878 :     push(@pins,[$neigh,map { $_->[1] } @$pairs]);
3879 :     foreach $y (@$pairs)
3880 :     {
3881 :     $peg2 = $y->[0];
3882 :     if ($peg2 =~ /^fig\|(\d+\.\d+)/)
3883 :     {
3884 :     $projection{$1}->{$peg2} = 1;
3885 :     }
3886 :     }
3887 :     }
3888 :     @corr = ();
3889 :     @orgs = keys(%projection);
3890 :     if (@orgs > 0)
3891 :     {
3892 :     foreach $genome (sort { $a <=> $b } @orgs)
3893 :     {
3894 :     push(@corr,sort { &FIG::by_fig_id($a,$b) } keys(%{$projection{$genome}}));
3895 :     }
3896 :     push(@pins,[$peg,@corr]);
3897 :     }
3898 :    
3899 :     foreach $cluster (@clusters)
3900 :     {
3901 :     $self->add_chromosomal_cluster($cluster);
3902 :     }
3903 :    
3904 :     foreach $pin (@pins)
3905 :     {
3906 :     $self->add_pch_pin($pin);
3907 :     }
3908 :     }
3909 :     }
3910 :    
3911 :     sub coupling_ev {
3912 :     my($self,$genome1,$sim1,$sim2,$bound) = @_;
3913 :     my($ev,$sc,$i,$j);
3914 :    
3915 :     $ev = [];
3916 :    
3917 :     $i = 0;
3918 :     $j = 0;
3919 :     while (($i < @$sim1) && ($j < @$sim2))
3920 :     {
3921 :     if ($sim1->[$i]->[0] < $sim2->[$j]->[0])
3922 :     {
3923 :     $i++;
3924 :     }
3925 :     elsif ($sim1->[$i]->[0] > $sim2->[$j]->[0])
3926 :     {
3927 :     $j++;
3928 :     }
3929 :     else
3930 :     {
3931 : overbeek 1.193 $self->accumulate_ev($genome1,$sim1->[$i]->[1],$sim2->[$j]->[1],$bound,$ev);
3932 : efrank 1.1 $i++;
3933 :     $j++;
3934 :     }
3935 :     }
3936 : overbeek 1.43 return ($self->score([map { $_->[0] } @$ev]),$ev);
3937 : efrank 1.1 }
3938 :    
3939 :     sub accumulate_ev {
3940 :     my($self,$genome1,$feature_ids1,$feature_ids2,$bound,$ev) = @_;
3941 : overbeek 1.43 my($genome2,@locs1,@locs2,$i,$j,$x);
3942 : efrank 1.1
3943 :     if ((@$feature_ids1 == 0) || (@$feature_ids2 == 0)) { return 0 }
3944 :    
3945 :     $feature_ids1->[0] =~ /^fig\|(\d+\.\d+)/;
3946 :     $genome2 = $1;
3947 :     @locs1 = map { $x = $self->feature_location($_); $x ? [&boundaries_of($x)] : () } @$feature_ids1;
3948 :     @locs2 = map { $x = $self->feature_location($_); $x ? [&boundaries_of($x)] : () } @$feature_ids2;
3949 :    
3950 :     for ($i=0; ($i < @$feature_ids1); $i++)
3951 :     {
3952 :     for ($j=0; ($j < @$feature_ids2); $j++)
3953 :     {
3954 : parrello 1.200 if (($feature_ids1->[$i] ne $feature_ids2->[$j]) &&
3955 : efrank 1.1 &close_enough($locs1[$i],$locs2[$j],$bound))
3956 :     {
3957 :     push(@$ev,[$feature_ids1->[$i],$feature_ids2->[$j]]);
3958 :     }
3959 :     }
3960 :     }
3961 :     }
3962 :    
3963 :     sub close_enough {
3964 : olson 1.111 shift if UNIVERSAL::isa($_[0],__PACKAGE__);
3965 : efrank 1.1 my($locs1,$locs2,$bound) = @_;
3966 :    
3967 :     # print STDERR &Dumper(["close enough",$locs1,$locs2]);
3968 :     return (($locs1->[0] eq $locs2->[0]) && (abs((($locs1->[1]+$locs1->[2])/2) - (($locs2->[1]+$locs2->[2])/2)) <= $bound));
3969 :     }
3970 :    
3971 :     sub acceptably_close {
3972 :     my($self,$feature_id,$sim_cutoff) = @_;
3973 :     my(%by_org,$id2,$genome,$sim);
3974 :    
3975 :     my($ans) = [];
3976 :    
3977 : overbeek 1.31 foreach $sim ($self->sims($feature_id,1000,$sim_cutoff,"fig"))
3978 : efrank 1.1 {
3979 :     $id2 = $sim->id2;
3980 :     if ($id2 =~ /^fig\|(\d+\.\d+)/)
3981 :     {
3982 :     my $genome = $1;
3983 : overbeek 1.51 if (! $self->is_eukaryotic($genome))
3984 : efrank 1.1 {
3985 :     push(@{$by_org{$genome}},$id2);
3986 :     }
3987 :     }
3988 :     }
3989 :     foreach $genome (sort { $a <=> $b } keys(%by_org))
3990 :     {
3991 :     push(@$ans,[$genome,$by_org{$genome}]);
3992 :     }
3993 :     return $ans;
3994 :     }
3995 :    
3996 :     ################ Translations of PEGsand External Protein Sequences ##########################
3997 :    
3998 :    
3999 :     =pod
4000 :    
4001 :     =head1 translatable
4002 :    
4003 :     usage: $fig->translatable($prot_id)
4004 :    
4005 :     The system takes any number of sources of protein sequences as input (and builds an nr
4006 :     for the purpose of computing similarities). For each of these input fasta files, it saves
4007 :     (in the DB) a filename, seek address and length so that it can go get the translation if
4008 :     needed. This routine simply returns true iff info on the translation exists.
4009 :    
4010 :     =cut
4011 :    
4012 :     sub translatable {
4013 :     my($self,$prot) = @_;
4014 :    
4015 :     return &translation_length($self,$prot) ? 1 : 0;
4016 :     }
4017 : parrello 1.200
4018 : efrank 1.1
4019 :     =pod
4020 :    
4021 :     =head1 translation_length
4022 :    
4023 :     usage: $len = $fig->translation_length($prot_id)
4024 :    
4025 :     The system takes any number of sources of protein sequences as input (and builds an nr
4026 :     for the purpose of computing similarities). For each of these input fasta files, it saves
4027 :     (in the DB) a filename, seek address and length so that it can go get the translation if
4028 :     needed. This routine returns the length of a translation. This does not require actually
4029 :     retrieving the translation.
4030 :    
4031 :     =cut
4032 :    
4033 :     sub translation_length {
4034 :     my($self,$prot) = @_;
4035 :    
4036 : overbeek 1.136 if ($self->is_deleted_fid($prot)) { return undef }
4037 :    
4038 : efrank 1.1 $prot =~ s/^([^\|]+\|[^\|]+)\|.*$/$1/;
4039 :     my $rdbH = $self->db_handle;
4040 : parrello 1.200 my $relational_db_response = $rdbH->SQL("SELECT slen,seek FROM protein_sequence_seeks
4041 : efrank 1.1 WHERE id = \'$prot\' ");
4042 :    
4043 : overbeek 1.145 my @vals = sort { $b->[1] <=> $a->[1] } @$relational_db_response;
4044 :     return (@vals > 0) ? $vals[0]->[0] : undef;
4045 : efrank 1.1 }
4046 :    
4047 :    
4048 :     =pod
4049 :    
4050 :     =head1 get_translation
4051 :    
4052 :     usage: $translation = $fig->get_translation($prot_id)
4053 :    
4054 :     The system takes any number of sources of protein sequences as input (and builds an nr
4055 :     for the purpose of computing similarities). For each of these input fasta files, it saves
4056 :     (in the DB) a filename, seek address and length so that it can go get the translation if
4057 :     needed. This routine returns a protein sequence.
4058 :    
4059 :     =cut
4060 :    
4061 :     sub get_translation {
4062 :     my($self,$id) = @_;
4063 :     my($rdbH,$relational_db_response,$fileN,$file,$fh,$seek,$ln,$tran);
4064 :    
4065 : overbeek 1.136 if ($self->is_deleted_fid($id)) { return '' }
4066 :    
4067 : efrank 1.1 $rdbH = $self->db_handle;
4068 :     $id =~ s/^([^\|]+\|[^\|]+)\|.*$/$1/;
4069 :    
4070 :     $relational_db_response = $rdbH->SQL("SELECT fileno, seek, len FROM protein_sequence_seeks WHERE id = \'$id\' ");
4071 :    
4072 : overbeek 1.256 if ((! ($relational_db_response && @$relational_db_response > 0)) &&
4073 :     ($id !~ /^fig\|/) &&
4074 :     ($id = $self->by_alias($id)))
4075 :     {
4076 :     $relational_db_response = $rdbH->SQL("SELECT fileno, seek, len FROM protein_sequence_seeks WHERE id = \'$id\' ");
4077 :     }
4078 :    
4079 : overbeek 1.145 if ($relational_db_response && @$relational_db_response > 0)
4080 : efrank 1.1 {
4081 : overbeek 1.145 my @vals = sort { $b->[1] <=> $a->[1] } @$relational_db_response;
4082 :     ($fileN,$seek,$ln) = @{$vals[0]};
4083 : efrank 1.1 if (($fh = $self->openF($self->N2file($fileN))) &&
4084 :     ($ln > 10))
4085 :     {
4086 :     seek($fh,$seek,0);
4087 :     read($fh,$tran,$ln-1);
4088 :     $tran =~ s/\s//g;
4089 :     return $tran;
4090 :     }
4091 :     }
4092 :     return '';
4093 :     }
4094 :    
4095 :     =pod
4096 :    
4097 :     =head1 mapped_prot_ids
4098 :    
4099 :     usage: @mapped = $fig->mapped_prot_ids($prot)
4100 :    
4101 :     This routine is at the heart of maintaining synonyms for protein sequences. The system
4102 :     determines which protein sequences are "essentially the same". These may differ in length
4103 :     (presumably due to miscalled starts), but the tails are identical (and the heads are not "too" extended).
4104 :     Anyway, the set of synonyms is returned as a list of 2-tuples [Id,length] sorted
4105 : parrello 1.200 by length.
4106 : efrank 1.1
4107 :     =cut
4108 :    
4109 :     sub mapped_prot_ids {
4110 :     my($self,$id) = @_;
4111 : parrello 1.200
4112 : overbeek 1.136 if ($self->is_deleted_fid($id)) { return () }
4113 :    
4114 : efrank 1.1 my $rdbH = $self->db_handle;
4115 :     my $relational_db_response = $rdbH->SQL("SELECT maps_to FROM peg_synonyms WHERE syn_id = \'$id\' ");
4116 :     if ($relational_db_response && (@$relational_db_response == 1))
4117 :     {
4118 :     $id = $relational_db_response->[0]->[0];
4119 :     }
4120 :    
4121 :     $relational_db_response = $rdbH->SQL("SELECT syn_id,syn_ln,maps_to_ln FROM peg_synonyms WHERE maps_to = \'$id\' ");
4122 :     if ($relational_db_response && (@$relational_db_response > 0))
4123 :     {
4124 :     return ([$id,$relational_db_response->[0]->[2]],map { [$_->[0],$_->[1]] } @$relational_db_response);
4125 :     }
4126 :     else
4127 :     {
4128 : overbeek 1.241 my $len = $self->translation_length($id);
4129 :     if ($len)
4130 :     {
4131 :     return ([$id,$len]);
4132 :     }
4133 :     else
4134 :     {
4135 :     return ();
4136 :     }
4137 : efrank 1.1 }
4138 : overbeek 1.14 }
4139 :    
4140 :     sub maps_to_id {
4141 :     my($self,$id) = @_;
4142 : parrello 1.200
4143 : overbeek 1.14 my $rdbH = $self->db_handle;
4144 :     my $relational_db_response = $rdbH->SQL("SELECT maps_to FROM peg_synonyms WHERE syn_id = \'$id\' ");
4145 :     return ($relational_db_response && (@$relational_db_response == 1)) ? $relational_db_response->[0]->[0] : $id;
4146 : efrank 1.1 }
4147 :    
4148 : olson 1.244 ################ GFF3 utilities ##########################
4149 :    
4150 :     sub get_gff_writer
4151 :     {
4152 :     my($self, %options) = @_;
4153 :    
4154 : olson 1.246 my $w = GFFWriter->new($self, %options);
4155 : olson 1.244
4156 :     return $w;
4157 :     }
4158 :    
4159 :    
4160 : efrank 1.1 ################ Assignments of Function to PEGs ##########################
4161 :    
4162 : overbeek 1.146 # set to undef to unset user
4163 :     #
4164 :     sub set_user {
4165 :     my($self,$user) = @_;
4166 :    
4167 :     $self->{_user} = $user;
4168 :     }
4169 :    
4170 :     sub get_user {
4171 :     my($self) = @_;
4172 :    
4173 :     return $self->{_user};
4174 :     }
4175 :    
4176 : efrank 1.1 =pod
4177 :    
4178 :     =head1 function_of
4179 :    
4180 :     usage: @functions = $fig->function_of($peg) OR
4181 :     $function = $fig->function_of($peg,$user)
4182 : parrello 1.200
4183 : efrank 1.1 In a list context, you get back a list of 2-tuples. Each 2-tuple is of the
4184 :     form [MadeBy,Function].
4185 :    
4186 :     In a scalar context,
4187 :    
4188 :     1. user is "master" if not specified
4189 :     2. function returned is the user's, if one exists; otherwise, master's, if one exists
4190 :    
4191 :     In a scalar context, you get just the function.
4192 :    
4193 :     =cut
4194 :    
4195 :     # Note that we do not return confidence. I propose a separate function to get both
4196 :     # function and confidence
4197 :     #
4198 :     sub function_of {
4199 :     my($self,$id,$user) = @_;
4200 :     my($relational_db_response,@tmp,$entry,$i);
4201 :     my $wantarray = wantarray();
4202 :     my $rdbH = $self->db_handle;
4203 :    
4204 : overbeek 1.136 if ($self->is_deleted_fid($id)) { return $wantarray ? () : "" }
4205 :    
4206 : efrank 1.1 if (($id =~ /^fig\|(\d+\.\d+\.peg\.\d+)/) && ($wantarray || $user))
4207 :     {
4208 :     if (($relational_db_response = $rdbH->SQL("SELECT made_by,assigned_function FROM assigned_functions WHERE ( prot = \'$id\' )")) &&
4209 :     (@$relational_db_response >= 1))
4210 :     {
4211 : overbeek 1.191 @tmp = sort { $a->[0] cmp $b->[0] } map { $_->[1] =~ s/^\s//; $_->[1] =~ s/(\t\S)?\s*$//; [$_->[0],$_->[1]] } @$relational_db_response;
4212 : efrank 1.1 for ($i=0; ($i < @tmp) && ($tmp[$i]->[0] ne "master"); $i++) {}
4213 :     if ($i < @tmp)
4214 :     {
4215 :     $entry = splice(@tmp,$i,1);
4216 :     unshift @tmp, ($entry);
4217 :     }
4218 :    
4219 :     my $val;
4220 :     if ($wantarray) { return @tmp }
4221 :     elsif ($user && ($val = &extract_by_who(\@tmp,$user))) { return $val }
4222 :     elsif ($user && ($val = &extract_by_who(\@tmp,"master"))) { return $val }
4223 :     else { return "" }
4224 :     }
4225 :     }
4226 :     else
4227 :     {
4228 :     if (($relational_db_response = $rdbH->SQL("SELECT assigned_function FROM assigned_functions WHERE ( prot = \'$id\' AND made_by = \'master\' )")) &&
4229 :     (@$relational_db_response >= 1))
4230 :     {
4231 : parrello 1.200 $relational_db_response->[0]->[0] =~ s/^\s//; $relational_db_response->[0]->[0] =~ s/(\t\S)?\s*$//;
4232 : efrank 1.1 return $wantarray ? (["master",$relational_db_response->[0]->[0]]) : $relational_db_response->[0]->[0];
4233 :     }
4234 :     }
4235 : parrello 1.200
4236 : efrank 1.1 return $wantarray ? () : "";
4237 :     }
4238 :    
4239 :     =pod
4240 :    
4241 :     =head1 translated_function_of
4242 :    
4243 :     usage: $function = $fig->translated_function_of($peg,$user)
4244 : parrello 1.200
4245 : efrank 1.1 You get just the translated function.
4246 :    
4247 :     =cut
4248 :    
4249 :     sub translated_function_of {
4250 :     my($self,$id,$user) = @_;
4251 :    
4252 :     my $func = $self->function_of($id,$user);
4253 :     if ($func)
4254 :     {
4255 :     $func = $self->translate_function($func);
4256 :     }
4257 :     return $func;
4258 :     }
4259 : parrello 1.200
4260 : efrank 1.1
4261 :     sub extract_by_who {
4262 : olson 1.111 shift if UNIVERSAL::isa($_[0],__PACKAGE__);
4263 : efrank 1.1 my($xL,$who) = @_;
4264 :     my($i);
4265 :    
4266 :     for ($i=0; ($i < @$xL) && ($xL->[$i]->[0] ne $who); $i++) {}
4267 :     return ($i < @$xL) ? $xL->[$i]->[1] : "";
4268 :     }
4269 :    
4270 :    
4271 :     =pod
4272 :    
4273 :     =head1 translate_function
4274 :    
4275 :     usage: $translated_func = $fig->translate_function($func)
4276 :    
4277 : parrello 1.200 Translates a function based on the function.synonyms table.
4278 : efrank 1.1
4279 :     =cut
4280 :    
4281 :     sub translate_function {
4282 :     my($self,$function) = @_;
4283 :    
4284 :     my ($tran,$from,$to,$line);
4285 :     if (! ($tran = $self->{_function_translation}))
4286 :     {
4287 :     $tran = {};
4288 :     if (open(TMP,"<$FIG_Config::global/function.synonyms"))
4289 :     {
4290 :     while (defined($line = <TMP>))
4291 :     {
4292 : golsen 1.44 chomp $line;
4293 : efrank 1.1 ($from,$to) = split(/\t/,$line);
4294 :     $tran->{$from} = $to;
4295 :     }
4296 :     close(TMP);
4297 :     }
4298 : overbeek 1.22 foreach $from (keys(%$tran))
4299 :     {
4300 :     $to = $tran->{$from};
4301 :     if ($tran->{$to})
4302 :     {
4303 :     delete $tran->{$from};
4304 :     }
4305 :     }
4306 : efrank 1.1 $self->{_function_translation} = $tran;
4307 :     }
4308 : overbeek 1.4
4309 :     while ($to = $tran->{$function})
4310 :     {
4311 :     $function = $to;
4312 :     }
4313 :     return $function;
4314 : efrank 1.1 }
4315 : parrello 1.200
4316 : efrank 1.1 =pod
4317 :    
4318 :     =head1 assign_function
4319 :    
4320 :     usage: $fig->assign_function($peg,$user,$function,$confidence)
4321 :    
4322 :     Assigns a function. Note that confidence can (and should be if unusual) included.
4323 :     Note that no annotation is written. This should normally be done in a separate
4324 :     call of the form
4325 :    
4326 : parrello 1.200
4327 : efrank 1.1
4328 :     =cut
4329 :    
4330 :     sub assign_function {
4331 :     my($self,$peg,$user,$function,$confidence) = @_;
4332 : overbeek 1.229 my($role,$roleQ,$kvs,$kv,$k,$v);
4333 : efrank 1.1
4334 : overbeek 1.197 if (! $self->is_real_feature($peg)) { return 0 }
4335 : overbeek 1.136
4336 : overbeek 1.229 my $genome = $self->genome_of($peg);
4337 :    
4338 : overbeek 1.228 $function =~ s/\s+/ /sg;
4339 :     $function =~ s/^\s+//;
4340 :     $function =~ s/\s+$//;
4341 : overbeek 1.229 if ($function =~ /^(.*?)\!(.*)/)
4342 :     {
4343 :     ($function,$kvs) = ($1,$2);
4344 :     if ($kvs)
4345 :     {
4346 :     $kvs =~ s/^\s+//;
4347 :     $kvs =~ s/\s+$//;
4348 :     foreach $kv (split(/\s+\!\s+/,$kvs))
4349 :     {
4350 :     if ($kv =~ /^([A-Za-z0-9._\-\+\%]+)\s+\^\s+(.*)$/)
4351 :     {
4352 :     ($k,$v) = ($1,$2);
4353 :     if ($v !~ /\S/)
4354 :     {
4355 :     &replace_peg_key_value($self,$peg,$k,"");
4356 :     }
4357 :     else
4358 :     {
4359 :     &replace_peg_key_value($self,$peg,$k,$v);
4360 :     }
4361 :     }
4362 :     elsif ($kv =~ /^([A-Za-z0-9._\-\+\%]+)$/)
4363 :     {
4364 :     &replace_peg_key_value($self,$peg,$1,1);
4365 :     }
4366 :     }
4367 :     }
4368 :     }
4369 : overbeek 1.228
4370 : efrank 1.1 my $rdbH = $self->db_handle;
4371 :     $confidence = $confidence ? $confidence : "";
4372 :    
4373 :     $rdbH->SQL("DELETE FROM assigned_functions WHERE ( prot = \'$peg\' AND made_by = \'$user\' )");
4374 :    
4375 :     my $funcQ = quotemeta $function;
4376 :     $rdbH->SQL("INSERT INTO assigned_functions ( prot, made_by, assigned_function, quality, org ) VALUES ( \'$peg\', \'$user\', \'$funcQ\', \'$confidence\', \'$genome\' )");
4377 :     $rdbH->SQL("DELETE FROM roles WHERE ( prot = \'$peg\' AND made_by = \'$user\' )");
4378 :    
4379 :     foreach $role (&roles_of_function($function))
4380 :     {
4381 :     $roleQ = quotemeta $role;
4382 :     $rdbH->SQL("INSERT INTO roles ( prot, role, made_by, org ) VALUES ( \'$peg\', '$roleQ\', \'$user\', \'$genome\' )");
4383 :     }
4384 :    
4385 :     &verify_dir("$FIG_Config::organisms/$genome/UserModels");
4386 : parrello 1.200 if ($user ne "master")
4387 : efrank 1.1 {
4388 :     &verify_dir("$FIG_Config::organisms/$genome/UserModels/$user");
4389 :     }
4390 :    
4391 : overbeek 1.66 my $file;
4392 :     if ((($user eq "master") && ($file = "$FIG_Config::organisms/$genome/assigned_functions") && open(TMP,">>$file")) ||
4393 :     (($user ne "master") && ($file = "$FIG_Config::organisms/$genome/UserModels/$user/assigned_functions") && open(TMP,">>$file")))
4394 : efrank 1.1 {
4395 :     flock(TMP,LOCK_EX) || confess "cannot lock assigned_functions";
4396 :     seek(TMP,0,2) || confess "failed to seek to the end of the file";
4397 :     print TMP "$peg\t$function\t$confidence\n";
4398 :     close(TMP);
4399 : overbeek 1.66 chmod(0777,$file);
4400 : efrank 1.1 return 1;
4401 :     }
4402 : overbeek 1.125 else
4403 :     {
4404 :     print STDERR "FAILED ASSIGNMENT: $peg\t$function\t$confidence\n";
4405 :     }
4406 : efrank 1.1 return 0;
4407 :     }
4408 :    
4409 :     sub hypo {
4410 : olson 1.111 shift if UNIVERSAL::isa($_[0],__PACKAGE__);
4411 : efrank 1.1 my $x = (@_ == 1) ? $_[0] : $_[1];
4412 :    
4413 : overbeek 1.23 if (! $x) { return 1 }
4414 :     if ($x =~ /hypoth/i) { return 1 }
4415 :     if ($x =~ /conserved protein/i) { return 1 }
4416 : overbeek 1.63 if ($x =~ /gene product/i) { return 1 }
4417 :     if ($x =~ /interpro/i) { return 1 }
4418 :     if ($x =~ /B[sl][lr]\d/i) { return 1 }
4419 :     if ($x =~ /^U\d/) { return 1 }
4420 :     if ($x =~ /^orf/i) { return 1 }
4421 :     if ($x =~ /uncharacterized/i) { return 1 }
4422 : overbeek 1.285 if ($x =~ /pseudogene/i) { return 1 }
4423 : overbeek 1.63 if ($x =~ /^predicted/i) { return 1 }
4424 :     if ($x =~ /AGR_/) { return 1 }
4425 : overbeek 1.51 if ($x =~ /similar to/i) { return 1 }
4426 : overbeek 1.63 if ($x =~ /similarity/i) { return 1 }
4427 :     if ($x =~ /glimmer/i) { return 1 }
4428 : overbeek 1.23 if ($x =~ /unknown/i) { return 1 }
4429 : overbeek 1.204 if (($x =~ /domain/i) ||
4430 :     ($x =~ /^y[a-z]{2,4}\b/i) ||
4431 :     ($x =~ /complete/i) ||
4432 : overbeek 1.250 ($x =~ /ensang/i) ||
4433 :     ($x =~ /unnamed/i) ||
4434 :     ($x =~ /EG:/) ||
4435 : overbeek 1.261 ($x =~ /orf\d+/i) ||
4436 :     ($x =~ /RIKEN/) ||
4437 :     ($x =~ /Expressed/i) ||
4438 : overbeek 1.250 ($x =~ /[a-zA-Z]{2,3}\|/) ||
4439 : overbeek 1.204 ($x =~ /predicted by Psort/) ||
4440 : overbeek 1.285 ($x =~ /predicted coding region/) ||
4441 : overbeek 1.204 ($x =~ /^bh\d+/i) ||
4442 :     ($x =~ /cds_/i) ||
4443 :     ($x =~ /^[a-z]{2,3}\d+/i) ||
4444 :     ($x =~ /similar to/i) ||
4445 :     ($x =~ / identi/i) ||
4446 : overbeek 1.275 ($x =~ /ortholog of/i) ||
4447 : overbeek 1.204 ($x =~ /structural feature/i)) { return 1 }
4448 : overbeek 1.23 return 0;
4449 : efrank 1.1 }
4450 :    
4451 :     ############################ Similarities ###############################
4452 :    
4453 :     =pod
4454 :    
4455 :     =head1 sims
4456 :    
4457 :     usage: @sims = $fig->sims($peg,$maxN,$maxP,$select)
4458 :    
4459 :     Returns a list of similarities for $peg such that
4460 :    
4461 :     there will be at most $maxN similarities,
4462 :    
4463 :     each similarity will have a P-score <= $maxP, and
4464 :    
4465 :     $select gives processing instructions:
4466 :    
4467 :     "raw" means that the similarities will not be expanded (by far fastest option)
4468 :     "fig" means return only similarities to fig genes
4469 :     "all" means that you want all the expanded similarities.
4470 :    
4471 :     By "expanded", we refer to taking a "raw similarity" against an entry in the non-redundant
4472 : parrello 1.200 protein collection, and converting it to a set of similarities (one for each of the
4473 : efrank 1.1 proteins that are essentially identical to the representative in the nr).
4474 :    
4475 : redwards 1.188 Each entry in @sims is a refence to an array. These are the values in each array position:
4476 :    
4477 : golsen 1.238 0. The query peg
4478 :     1. The similar peg
4479 :