Annotation of /FigKernelPackages/FIG.pm

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

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