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

Annotation of /FigKernelPackages/FIG.pm

Parent Directory Parent Directory | Revision Log Revision Log


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

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

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3