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

Annotation of /FigKernelPackages/FIG.pm

Parent Directory Parent Directory | Revision Log Revision Log


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