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

Annotation of /FigKernelPackages/FIG.pm

Parent Directory Parent Directory | Revision Log Revision Log


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

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