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

Annotation of /FigKernelPackages/FIG.pm

Parent Directory Parent Directory | Revision Log Revision Log


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

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

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3