Annotation of /FigKernelPackages/FIG.pm

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