Annotation of /FigKernelPackages/FIG.pm

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

MCS Webmaster	ViewVC Help
Powered by ViewVC 1.0.3