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

Annotation of /FigKernelPackages/FIG.pm

Parent Directory Parent Directory | Revision Log Revision Log


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

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