Annotation of /FigKernelPackages/FIG.pm

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

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