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

Annotation of /FigKernelPackages/FigFam.pm

Parent Directory Parent Directory | Revision Log Revision Log


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

1 : overbeek 1.1 #
2 :     # Copyright (c) 2003-2006 University of Chicago and Fellowship
3 :     # for Interpretations of Genomes. All Rights Reserved.
4 :     #
5 :     # This file is part of the SEED Toolkit.
6 :     #
7 :     # The SEED Toolkit is free software. You can redistribute
8 :     # it and/or modify it under the terms of the SEED Toolkit
9 :     # Public License.
10 :     #
11 :     # You should have received a copy of the SEED Toolkit Public License
12 :     # along with this program; if not write to the University of Chicago
13 :     # at info@ci.uchicago.edu or the Fellowship for Interpretation of
14 :     # Genomes at veronika@thefig.info or download a copy from
15 :     # http://www.theseed.org/LICENSE.TXT.
16 :     #
17 :    
18 :     package FigFam;
19 :    
20 :     use strict;
21 :     use FIG;
22 :     use SameFunc;
23 : overbeek 1.4
24 :     use Carp;
25 : overbeek 1.1 use Data::Dumper;
26 :    
27 :     # This is the constructor. Presumably, $class is 'FigFam'.
28 :     #
29 :     # $fam_id is the ID of the family (of the form FIGnnnnnnn). It is required.
30 :     #
31 : overbeek 1.3 # $fam_data is optional. If it exists,
32 :     # it gives the directory that is (or will) contain the FigFam data.
33 : overbeek 1.5 # $FIG_Config::data/FigfamsData will be used as a default.
34 : overbeek 1.1 #
35 : overbeek 1.3 # $fam_file must exist within $fam_data. It contains a 2-column version of the FIG families.
36 :     # If the $fam_dir/yyy/FIGxxxxyyy directory (which is where the directory for the specific family
37 : overbeek 1.1 # will go) does not exist, the specification of the PEGs in $fam_file will be used to initialize
38 :     # the family.
39 :     #
40 :    
41 :     sub new {
42 : overbeek 1.3 my($class,$fig,$fam_id,$fam_data) = @_;
43 : overbeek 1.1
44 : overbeek 1.25 ($fam_id =~ /^FIG\d{6}$/) || confess "invalid family id: $fam_id";
45 :    
46 : overbeek 1.1 my $fam = {};
47 :     $fam->{id} = $fam_id;
48 :     $fam->{fig} = $fig;
49 :    
50 : overbeek 1.3 if (! defined($fam_data))
51 : overbeek 1.1 {
52 : overbeek 1.3 $fam_data = "$FIG_Config::data/FigfamsData";
53 : overbeek 1.1 }
54 : overbeek 1.3 $fam->{data} = $fam_data;
55 :    
56 :     my $fam_dir = "$fam_data/FIGFAMS";
57 :     my $fam_file = "$fam_data/families.2c";
58 : overbeek 1.10 my $dir = &fam_dir($fam_data,$fam_id);
59 : overbeek 1.3 $fam->{dir} = $dir;
60 :    
61 : overbeek 1.1 &FIG::verify_dir($dir);
62 :     if ((! -s "$dir/PEGs") && defined($fam_file))
63 :     {
64 : overbeek 1.4 if (-s $fam_file)
65 : overbeek 1.3 {
66 : overbeek 1.4 system("grep $fam_id \"$fam_file\" > \"$dir/PEGs\"");
67 :     if (! -s "$dir/PEGs")
68 :     {
69 :     system "rm -rf $dir";
70 :     confess "PEG file $dir/PEGs does not exist --- $dir deleted";
71 :     return undef;
72 :     }
73 :     }
74 :     else
75 :     {
76 : overbeek 1.13 confess "Family file $fam_file does not exist";
77 : overbeek 1.3 return undef;
78 :     }
79 : overbeek 1.18 }
80 :    
81 :     if (! -s "$dir/function")
82 :     {
83 : overbeek 1.28 &reset_function($fam);
84 : overbeek 1.1 }
85 : overbeek 1.18
86 : overbeek 1.1 open(FUNC,"<$dir/function") || die "could not open $dir/function";
87 :     my $func = <FUNC>;
88 :     chomp $func;
89 :     close(FUNC);
90 :     $fam->{function} = $func;
91 :    
92 :     my($peg,$pegs);
93 : overbeek 1.10 my $pegsL = [sort { &FIG::by_fig_id($a,$b) } grep { $fig->is_real_feature($_) } map { $_ =~ /^\S+\t(\S+)/; $1 } `cat \"$dir/PEGs\"`];
94 : overbeek 1.1
95 :     $fam->{pegsL} = $pegsL;
96 :     my $pegsH = {};
97 :     foreach $peg (@$pegsL)
98 :     {
99 :     $pegsH->{$peg} = 1;
100 :     }
101 :     $fam->{pegsH} = $pegsH;
102 : overbeek 1.4
103 :    
104 :     if (-s "$dir/PEGs.fasta")
105 :     {
106 :     open(FASTA,"<$dir/PEGs.fasta") || die "could not read-open $dir/PEGs.fasta";
107 :     while (my ($peg, $seqP) = &FIG::read_fasta_record(\*FASTA))
108 :     {
109 :     $fam->{peg_lengths}->{$peg} = length($$seqP);
110 :     }
111 :     close(FASTA);
112 :     }
113 : overbeek 1.28 elsif (! -s "$dir/built")
114 : overbeek 1.1 {
115 : overbeek 1.4 open(FASTA,">$dir/PEGs.fasta") || die "could not write-open $dir/PEGs.fasta";
116 : overbeek 1.1 foreach $peg (@$pegsL)
117 :     {
118 :     my $seq = $fig->get_translation($peg);
119 :     if ($seq)
120 :     {
121 :     print FASTA ">$peg\n$seq\n";
122 : overbeek 1.3 $fam->{peg_lengths}->{$peg} = length($seq);
123 : overbeek 1.1 }
124 : overbeek 1.4 else
125 :     {
126 :     confess "Could not get seq for $peg";
127 :     }
128 : overbeek 1.1 }
129 :     close(FASTA);
130 : overbeek 1.4 }
131 :    
132 : overbeek 1.28 if ((! -s "$dir/built") && ((! -s "$dir/PEGs.fasta.pin") || ((-M "$dir/PEGs.fasta.pin") > (-M "$dir/PEGs.fasta"))))
133 : overbeek 1.4 {
134 : overbeek 1.3 &FIG::run("$FIG_Config::ext_bin/formatdb -i \"$dir/PEGs.fasta\" -p");
135 : overbeek 1.1 }
136 : overbeek 1.4
137 : overbeek 1.28 if ((! -s "$dir/bounds.sims") && (! -s "$dir/built"))
138 : overbeek 1.1 {
139 : overbeek 1.6 open(SIMS,">$dir/bounds.sims")
140 :     || die "could not open $dir/bounds.sims";
141 :    
142 : overbeek 1.1 foreach $peg (@$pegsL)
143 :     {
144 :     my @sims = $fig->sims($peg,1000,1.0e-20,"fig");
145 : overbeek 1.6 foreach my $sim (@sims)
146 :     {
147 :     print SIMS join("\t",$sim->id1,
148 :     $sim->id2,
149 :     $sim->bit_score,
150 :     scalar $fig->function_of($sim->id2)),"\n";
151 :     }
152 :     print SIMS "//\n";
153 :     }
154 :     close(SIMS);
155 :     }
156 :    
157 : overbeek 1.28 if ((! "$dir/built") && (! -f "$dir/bounds"))
158 : overbeek 1.6 {
159 : overbeek 1.8 print STDERR " building bounds for $fam_id\n";
160 : overbeek 1.6 open(BOUNDS,">$dir/bounds") || die "could not open $dir/bounds";
161 :     open(SIMS,"<$dir/bounds.sims") || die "could not open $dir/bounds.sims";
162 : overbeek 1.1
163 : overbeek 1.6 my($sims);
164 :     $/ = "\n//\n";
165 :     while (defined($sims = <SIMS>))
166 :     {
167 :     chomp $sims;
168 :     my @sims = sort { $b->[2] <=> $a->[2] } map { [split(/\t/,$_)] } split(/\n/,$sims);
169 :     if ($peg = $sims[0]->[0])
170 :     {
171 :     my($i,$hits,$bad,$bad_sc,$last_good,$func2);
172 :     for ($i=0, $hits=0, $bad="", $bad_sc=0, $last_good=1; ($i < @sims) && (! $bad); $i++)
173 :     {
174 :     my $sim = $sims[$i];
175 :     my $id2 = $sim->[1];
176 :     if ($pegsH->{$id2})
177 :     {
178 :     $last_good = $sim->[2];
179 :     $hits++;
180 :     }
181 :     else
182 : overbeek 1.1 {
183 : overbeek 1.6 $func2 = $sim->[3];
184 : overbeek 1.8 if (! &ok_func($fig,$func,$func2,$id2,$fam))
185 : overbeek 1.6 {
186 :     $bad = $id2;
187 :     $bad_sc = $sim->[2];
188 :     }
189 : overbeek 1.1 }
190 :     }
191 :    
192 : overbeek 1.6 if ($hits > 0)
193 :     {
194 :     print BOUNDS join("\t",($peg,$hits,$last_good,$bad,$bad_sc)),"\n";
195 :     }
196 : overbeek 1.1 }
197 :     }
198 :     close(BOUNDS);
199 : overbeek 1.6 close(SIMS);
200 :     $/ = "\n";
201 : overbeek 1.1 }
202 : overbeek 1.6
203 : overbeek 1.1 open(BOUNDS,"<$dir/bounds") || die "could not open $dir/bounds";
204 :    
205 :     my $bounds = {};
206 :     my $x;
207 :     while (defined($x = <BOUNDS>))
208 :     {
209 :     chomp $x;
210 :     my @flds = split(/\t/,$x);
211 :     my $peg =shift @flds;
212 :     $bounds->{$peg} = [@flds];
213 :     }
214 :     close(BOUNDS);
215 :     $fam->{bounds} = $bounds;
216 :    
217 : overbeek 1.3 my $blast_file = "$dir/blastout";
218 :     my $rep_file = "$dir/reps";
219 : overbeek 1.28 if ((! -f $blast_file) && (! -s "$dir/built"))
220 : overbeek 1.1 {
221 : overbeek 1.3 &FIG::run("split_and_trim_sequences $dir/Split blastout=$blast_file < $dir/PEGs.fasta");
222 :     if (-s $blast_file)
223 :     {
224 :     my @out = map { chop; [split(/\t/,$_) ] } `cat $blast_file`;
225 :     my @keep = ();
226 :     my $write = 0;
227 :     my $i;
228 :     for ($i=0; ($i < @out); $i++)
229 :     {
230 :     if ((($i == 0) ||
231 :     (($out[$i-1]->[0] ne $out[$i]->[0]) || ($out[$i-1]->[1] ne $out[$i]->[1]))) &&
232 :     ($out[$i-1]->[6] <= 1.0e-10))
233 :     {
234 :     push(@keep,$out[$i]);
235 :     }
236 :     else
237 :     {
238 :     $write = 1;
239 :     }
240 :     }
241 :     if ($write)
242 :     {
243 :     open(OUT,">$blast_file") || die $blast_file;
244 :     foreach my $x (@keep)
245 :     {
246 :     print OUT join("\t",@$x),"\n";
247 :     }
248 :     close(OUT);
249 :     }
250 :     }
251 : overbeek 1.1 }
252 : overbeek 1.3
253 : overbeek 1.28 if ((-s $blast_file) && (! -s $rep_file) && (! -s "$dir/built"))
254 : overbeek 1.3 {
255 :     my $n = 1;
256 :     open(REP,">$rep_file") || die "could not open $rep_file";
257 :     open(FASTA,"<$dir/PEGs.fasta") || die "could not open $dir/PEGs.fasta";
258 :     my %seen;
259 :     my($seq,$peg,$sim);
260 :     $/ = "\n>";
261 :     while (defined($_ = <FASTA>))
262 :     {
263 :     chomp;
264 :     if ($_ =~ /^>?(\S+)[^\n]*\n(.*)/s)
265 :     {
266 :     $peg = $1;
267 :     if (! $seen{$peg})
268 :     {
269 :     $seq = $2;
270 :     $seq =~ s/\s//gs;
271 :     print REP ">$fam_id-$n\n$seq\n";
272 :     $n++;
273 :    
274 :     $/ = "\n";
275 :     $seen{$peg} = 1;
276 :    
277 :     open(BLAST,"<$blast_file") || die "could not open $blast_file";
278 :     while (defined($sim = <BLAST>))
279 :     {
280 :     if (($sim =~ /^(\S+)\t(\S+)(\t\d+){4}\t(\S+)/) && ($1 eq $peg) &&
281 :     ($4 < 1.0e-10))
282 :     {
283 :     $seen{$2} = 1;
284 :     }
285 :     }
286 :     close(BLAST);
287 :     $/ = "\n>";
288 :     }
289 :     }
290 :     }
291 :     $/ = "\n";
292 :     close(FASTA);
293 :     close(REP);
294 :     }
295 : overbeek 1.28
296 :     $fam->{reps} = (-s $rep_file) ? [map { $_ =~ /^(\S+)/; $1 } `grep -v "^>" $rep_file`] : [];
297 : overbeek 1.3
298 : overbeek 1.28 if (! -s "$dir/built") { system "echo 1 > $dir/built" }
299 : overbeek 1.1 bless $fam,$class;
300 :     return $fam;
301 :     }
302 :    
303 : overbeek 1.28 sub reset_function {
304 :     my($self) = @_;
305 :    
306 :     my $fam_dir = $self->{dir};
307 :     my $fig = $self->{fig};
308 :    
309 :     my @pegs = grep { $fig->is_real_feature($_) } map { $_ =~ /^\S+\t(\S+)/; $1 } `cat "$fam_dir/PEGs"`;
310 :     my($peg,%funcs,$func,@subs);
311 :     foreach $peg (@pegs)
312 :     {
313 :     $func = $fig->function_of($peg);
314 :     @subs = grep { $fig->usable_subsystem($_) } $fig->peg_to_subsystems($peg);
315 :     my $incr = (@subs > 0) ? 2 : 1;
316 :     $funcs{$func} += $incr;
317 :     }
318 :     my @funcs = sort { $funcs{$b} <=> $funcs{$a} } keys(%funcs);
319 :     $func = (@funcs > 0) ? $funcs[0] : "hypothetical protein";
320 :     open(FUNC,">$fam_dir/function") || die "could not open $fam_dir/function";
321 : overbeek 1.29 print FUNC "$func\n";
322 : overbeek 1.28 close(FUNC);
323 :     }
324 :    
325 : overbeek 1.11 sub display {
326 :     my($self) = @_;
327 :    
328 :     my $fam_id = $self->family_id;
329 :     my $fig = $self->{fig};
330 :     my $fam_func = $self->family_function;
331 :     my $fam_dir = $self->{dir};
332 :     my @pegs = map { chomp; $_ } `cut -f2 $fam_dir/PEGs`;
333 :    
334 :     my $peg;
335 :     my @disp = ( "ID: $fam_id",
336 :     "Function: $fam_func",
337 :     "Directory: $fam_dir",
338 :     "PEGs:",
339 :     map { $peg = $_; " " . join("\t",($peg,
340 :     $fig->genus_species(&FIG::genome_of($peg)),
341 :     scalar $fig->function_of($_)))
342 :     } @pegs
343 :     );
344 :    
345 :     return join("\n", @disp) . "\n";
346 :     }
347 :    
348 : overbeek 1.1 sub ok_func {
349 : overbeek 1.8 my($fig,$func,$func2,$id2,$fam) = @_;
350 : overbeek 1.1
351 :     my $i;
352 : overbeek 1.8 if (&loose_same_func($func,$func2) || ($func eq $fig->function_of($id2))) { return 1 }
353 :     #
354 :     my $funcI;
355 :     if (defined($funcI = $fam->{ignorable_func}->{"$id2\t$func"}))
356 :     {
357 :     return $funcI;
358 :     }
359 : overbeek 1.7
360 : overbeek 1.1 my @sims = $fig->sims($id2,5,1.0e-30,"fig");
361 : overbeek 1.22 if ($ENV{VERBOSE} && ($ENV{VERBOSE} > 1)) { print STDERR &FIG::flatten_dumper(\@sims), "\n"; }
362 : overbeek 1.8 my($func3);
363 : overbeek 1.7 for ($i=0; ($i < @sims) && (defined($func3 = $fig->function_of($sims[$i]->id2))) &&
364 :     (&FIG::hypo($func3) || &SameFunc::same_func($func,$func3)); $i++) {}
365 : overbeek 1.17 if (($i == @sims) && ($ENV{'VERBOSE'})) { print STDERR "make assignment: $id2\t$func\n" }
366 : overbeek 1.8
367 :     $fam->{ignorable_func}->{"$id2\t$func"} = ($i == @sims);
368 :    
369 :     return $fam->{ignorable_func}->{"$id2\t$func"}
370 : overbeek 1.1 }
371 :    
372 : overbeek 1.8 sub member {
373 :     my($x,$yL) = @_;
374 :     my $i;
375 :    
376 :     for ($i=0; ($i < @$yL) && ($yL->[$i] ne $x); $i++) {}
377 :     return ($i < @$yL);
378 :     }
379 :    
380 : overbeek 1.26 sub list_members {
381 :     my ($self) = @_;
382 :    
383 :     my $fam_dir = $self->{dir};
384 :     my @pegs = map { chomp; $_ } `cut -f2 $fam_dir/PEGs`;
385 :    
386 :     return @pegs;
387 :     }
388 :    
389 : overbeek 1.8 sub loose_same_func {
390 :     my($f1,$f2) = @_;
391 :    
392 : overbeek 1.16 if (&SameFunc::same_func($f1,$f2,'strong')) { return 1 }
393 : overbeek 1.8 my @s1 = split(/\s*[;\/@]\s*/,$f1);
394 :     my @s2 = split(/\s*[;\/@]\s*/,$f2);
395 :     if ((@s1 == 1) && (@s2 > 1) && &member($s1[0],\@s2))
396 :     {
397 :     return 1;
398 :     }
399 :     elsif ((@s2 == 1) && (@s1 > 1) && &member($s2[0],\@s1))
400 :     {
401 :     return 1;
402 :     }
403 :     else
404 :     {
405 :     return 0;
406 :     }
407 :     }
408 :    
409 :    
410 :    
411 : overbeek 1.3 sub representatives {
412 :     my($self) = @_;
413 :    
414 :     my $reps = $self->{reps};
415 :     return @$reps;
416 :     }
417 :    
418 :     sub should_be_member {
419 : overbeek 1.9 my($self,$seq,$debug) = @_;
420 : overbeek 1.19
421 : overbeek 1.9 if ($debug) { $ENV{'VERBOSE'} = $ENV{'DEBUG'} = 1 }
422 : overbeek 1.19
423 : overbeek 1.8 my $fig = $self->{fig};
424 : overbeek 1.3 my $dir = $self->{dir};
425 :     open(TMP,">$FIG_Config::temp/tmp$$.fasta")
426 :     || die "could not open $FIG_Config::temp/tmp$$.fasta";
427 :     print TMP ">query\n$seq\n";
428 :     close(TMP);
429 : overbeek 1.19
430 : overbeek 1.25 if ($ENV{'DEBUG'}) { print STDERR &Dumper([$self->family_id,$seq]) }
431 : overbeek 1.19
432 : overbeek 1.3 my $query_ln = length($seq);
433 : overbeek 1.19
434 :     my @tmp = `blastall -i $FIG_Config::temp/tmp$$.fasta -d $dir/PEGs.fasta -m 8 -FF -p blastp`;
435 :    
436 : overbeek 1.3 my %seen;
437 :     my $should_be = 0;
438 : overbeek 1.7 my $yes = 0;
439 :     my $no = 0;
440 : overbeek 1.19
441 : overbeek 1.8 my $ln1 = length($seq);
442 : overbeek 1.15 my $good = 0;
443 : overbeek 1.19
444 :     my $sims = [];
445 :     foreach $_ (@tmp)
446 : overbeek 1.3 {
447 : overbeek 1.7 if ($ENV{'DEBUG'}) { print STDERR $_ }
448 : overbeek 1.3 chop;
449 : overbeek 1.19
450 : overbeek 1.3 my $sim = [split(/\t/,$_)];
451 :     my $peg = $sim->[1];
452 : overbeek 1.8 next if (($_ = $ENV{'IGNORE'}) && (&FIG::genome_of($peg) eq $_)); # for debugging
453 :    
454 : overbeek 1.3 my $sc = $sim->[10];
455 : overbeek 1.6 my $bit_score = $sim->[11];
456 : overbeek 1.8
457 :     my $matched1 = abs($sim->[7] - $sim->[6]) + 1;
458 :     my $matched2 = abs($sim->[9] - $sim->[8]) + 1;
459 :     my $ln2 = $self->{peg_lengths}->{$peg};
460 :    
461 :     if ((! $seen{$peg}) && ($sc <= 1.0e-10) &&
462 : overbeek 1.24 # ($matched1 >= (0.7) * $ln1) &&
463 :     ($matched2 >= (0.7) * $ln2))
464 : overbeek 1.3 {
465 :     $seen{$peg} = 1;
466 : overbeek 1.4 push @$sim, $query_ln, $self->{peg_lengths}->{$peg};
467 :     bless $sim, 'Sim';
468 :     push @$sims, $sim;
469 : overbeek 1.19
470 : overbeek 1.3 my $bounds = $self->{bounds}->{$peg};
471 : overbeek 1.19 if ($bounds && (@$sims <= 10))
472 : overbeek 1.3 {
473 : overbeek 1.25 if (($bit_score >= $bounds->[1]) && ((! $bounds->[2]) || $bounds->[3] < $bit_score))
474 : overbeek 1.7 {
475 : overbeek 1.8 if ($ENV{'DEBUG'}) { print STDERR " yes - $peg\n" }
476 : overbeek 1.19 ++$yes;
477 : overbeek 1.15 if ($yes > ($no+1)) { $good = 1 }
478 : overbeek 1.7 }
479 :     else
480 :     {
481 : overbeek 1.8 my $bad_func = $bounds->[2] ? $fig->function_of($bounds->[2]) : "";
482 :     if ($ENV{'DEBUG'}) { print STDERR " no - $peg ",join(",",@$bounds)," $bad_func\n" }
483 : overbeek 1.19 ++$no;
484 : overbeek 1.7 }
485 : overbeek 1.3 }
486 :     }
487 :     }
488 : overbeek 1.25 if ($ENV{'DEBUG'}) { print STDERR " yes=$yes no=$no good=$good\n" }
489 : overbeek 1.15 return ($good,$sims);
490 : overbeek 1.3 }
491 :    
492 : overbeek 1.5 sub family_function {
493 :     my($self) = @_;
494 :    
495 :     return $self->{function}
496 :     }
497 :    
498 : overbeek 1.7 sub family_id {
499 :     my($self) = @_;
500 :    
501 :     return $self->{id};
502 :     }
503 :    
504 : overbeek 1.10 sub fam_dir {
505 :     my($fam_data,$fam_id) = @_;
506 :    
507 :     $fam_id =~ /(\d{3})$/;
508 :     return "$fam_data/FIGFAMS/$1/$fam_id";
509 :     }
510 :    
511 : overbeek 1.1 1;

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3