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

Annotation of /FigKernelPackages/FIG.pm

Parent Directory Parent Directory | Revision Log Revision Log


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

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