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

Annotation of /FigKernelPackages/FigFam.pm

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1.26 - (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.1 my @first = `head -n 1 \"$dir/PEGs\"`;
84 : overbeek 1.18 my $func = ($first[0] && ($first[0] =~ /^\S+\t\S+\t(\S.*\S)/)) ? $1 : "hypothetical protein";
85 :     open(FUNC,">$dir/function") || die "could not open $dir/function";
86 :     print FUNC "$1\n";
87 :     close(FUNC);
88 : overbeek 1.1 }
89 : overbeek 1.18
90 : overbeek 1.1 open(FUNC,"<$dir/function") || die "could not open $dir/function";
91 :     my $func = <FUNC>;
92 :     chomp $func;
93 :     close(FUNC);
94 :     $fam->{function} = $func;
95 :    
96 :     my($peg,$pegs);
97 : 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\"`];
98 : overbeek 1.1
99 :     $fam->{pegsL} = $pegsL;
100 :     my $pegsH = {};
101 :     foreach $peg (@$pegsL)
102 :     {
103 :     $pegsH->{$peg} = 1;
104 :     }
105 :     $fam->{pegsH} = $pegsH;
106 : overbeek 1.4
107 :    
108 :     if (-s "$dir/PEGs.fasta")
109 :     {
110 :     open(FASTA,"<$dir/PEGs.fasta") || die "could not read-open $dir/PEGs.fasta";
111 :     while (my ($peg, $seqP) = &FIG::read_fasta_record(\*FASTA))
112 :     {
113 :     $fam->{peg_lengths}->{$peg} = length($$seqP);
114 :     }
115 :     close(FASTA);
116 :     }
117 :     else
118 : overbeek 1.1 {
119 : overbeek 1.4 open(FASTA,">$dir/PEGs.fasta") || die "could not write-open $dir/PEGs.fasta";
120 : overbeek 1.1 foreach $peg (@$pegsL)
121 :     {
122 :     my $seq = $fig->get_translation($peg);
123 :     if ($seq)
124 :     {
125 :     print FASTA ">$peg\n$seq\n";
126 : overbeek 1.3 $fam->{peg_lengths}->{$peg} = length($seq);
127 : overbeek 1.1 }
128 : overbeek 1.4 else
129 :     {
130 :     confess "Could not get seq for $peg";
131 :     }
132 : overbeek 1.1 }
133 :     close(FASTA);
134 : overbeek 1.4 }
135 :    
136 :     if ((! -s "$dir/PEGs.fasta.pin") || ((-M "$dir/PEGs.fasta.pin") > (-M "$dir/PEGs.fasta")))
137 :     {
138 : overbeek 1.3 &FIG::run("$FIG_Config::ext_bin/formatdb -i \"$dir/PEGs.fasta\" -p");
139 : overbeek 1.1 }
140 : overbeek 1.4
141 : overbeek 1.6 if (! -s "$dir/bounds.sims")
142 : overbeek 1.1 {
143 : overbeek 1.6 open(SIMS,">$dir/bounds.sims")
144 :     || die "could not open $dir/bounds.sims";
145 :    
146 : overbeek 1.1 foreach $peg (@$pegsL)
147 :     {
148 :     my @sims = $fig->sims($peg,1000,1.0e-20,"fig");
149 : overbeek 1.6 foreach my $sim (@sims)
150 :     {
151 :     print SIMS join("\t",$sim->id1,
152 :     $sim->id2,
153 :     $sim->bit_score,
154 :     scalar $fig->function_of($sim->id2)),"\n";
155 :     }
156 :     print SIMS "//\n";
157 :     }
158 :     close(SIMS);
159 :     }
160 :    
161 :     if (! -s "$dir/bounds")
162 :     {
163 : overbeek 1.8 print STDERR " building bounds for $fam_id\n";
164 : overbeek 1.6 open(BOUNDS,">$dir/bounds") || die "could not open $dir/bounds";
165 :     open(SIMS,"<$dir/bounds.sims") || die "could not open $dir/bounds.sims";
166 : overbeek 1.1
167 : overbeek 1.6 my($sims);
168 :     $/ = "\n//\n";
169 :     while (defined($sims = <SIMS>))
170 :     {
171 :     chomp $sims;
172 :     my @sims = sort { $b->[2] <=> $a->[2] } map { [split(/\t/,$_)] } split(/\n/,$sims);
173 :     if ($peg = $sims[0]->[0])
174 :     {
175 :     my($i,$hits,$bad,$bad_sc,$last_good,$func2);
176 :     for ($i=0, $hits=0, $bad="", $bad_sc=0, $last_good=1; ($i < @sims) && (! $bad); $i++)
177 :     {
178 :     my $sim = $sims[$i];
179 :     my $id2 = $sim->[1];
180 :     if ($pegsH->{$id2})
181 :     {
182 :     $last_good = $sim->[2];
183 :     $hits++;
184 :     }
185 :     else
186 : overbeek 1.1 {
187 : overbeek 1.6 $func2 = $sim->[3];
188 : overbeek 1.8 if (! &ok_func($fig,$func,$func2,$id2,$fam))
189 : overbeek 1.6 {
190 :     $bad = $id2;
191 :     $bad_sc = $sim->[2];
192 :     }
193 : overbeek 1.16 elsif ($peg eq "fig|262719.3.peg.75")
194 :     {
195 :     print STDERR "ok: $func\n$func2\n$id2\n";
196 :     }
197 : overbeek 1.1 }
198 :     }
199 :    
200 : overbeek 1.6 if ($hits > 0)
201 :     {
202 :     print BOUNDS join("\t",($peg,$hits,$last_good,$bad,$bad_sc)),"\n";
203 :     }
204 : overbeek 1.1 }
205 :     }
206 :     close(BOUNDS);
207 : overbeek 1.6 close(SIMS);
208 :     $/ = "\n";
209 : overbeek 1.1 }
210 : overbeek 1.6
211 : overbeek 1.1 open(BOUNDS,"<$dir/bounds") || die "could not open $dir/bounds";
212 :    
213 :     my $bounds = {};
214 :     my $x;
215 :     while (defined($x = <BOUNDS>))
216 :     {
217 :     chomp $x;
218 :     my @flds = split(/\t/,$x);
219 :     my $peg =shift @flds;
220 :     $bounds->{$peg} = [@flds];
221 :     }
222 :     close(BOUNDS);
223 :     $fam->{bounds} = $bounds;
224 :    
225 : overbeek 1.3 my $blast_file = "$dir/blastout";
226 :     my $rep_file = "$dir/reps";
227 :     if (! -d "$dir/Split")
228 : overbeek 1.1 {
229 : overbeek 1.3 &FIG::run("split_and_trim_sequences $dir/Split blastout=$blast_file < $dir/PEGs.fasta");
230 :     if (-s $blast_file)
231 :     {
232 :     my @out = map { chop; [split(/\t/,$_) ] } `cat $blast_file`;
233 :     my @keep = ();
234 :     my $write = 0;
235 :     my $i;
236 :     for ($i=0; ($i < @out); $i++)
237 :     {
238 :     if ((($i == 0) ||
239 :     (($out[$i-1]->[0] ne $out[$i]->[0]) || ($out[$i-1]->[1] ne $out[$i]->[1]))) &&
240 :     ($out[$i-1]->[6] <= 1.0e-10))
241 :     {
242 :     push(@keep,$out[$i]);
243 :     }
244 :     else
245 :     {
246 :     $write = 1;
247 :     }
248 :     }
249 :     if ($write)
250 :     {
251 :     open(OUT,">$blast_file") || die $blast_file;
252 :     foreach my $x (@keep)
253 :     {
254 :     print OUT join("\t",@$x),"\n";
255 :     }
256 :     close(OUT);
257 :     }
258 :     }
259 : overbeek 1.1 }
260 : overbeek 1.3
261 :     if ((-s $blast_file) && (! -s $rep_file))
262 :     {
263 :     my $n = 1;
264 :     open(REP,">$rep_file") || die "could not open $rep_file";
265 :     open(FASTA,"<$dir/PEGs.fasta") || die "could not open $dir/PEGs.fasta";
266 :     my %seen;
267 :     my($seq,$peg,$sim);
268 :     $/ = "\n>";
269 :     while (defined($_ = <FASTA>))
270 :     {
271 :     chomp;
272 :     if ($_ =~ /^>?(\S+)[^\n]*\n(.*)/s)
273 :     {
274 :     $peg = $1;
275 :     if (! $seen{$peg})
276 :     {
277 :     $seq = $2;
278 :     $seq =~ s/\s//gs;
279 :     print REP ">$fam_id-$n\n$seq\n";
280 :     $n++;
281 :    
282 :     $/ = "\n";
283 :     $seen{$peg} = 1;
284 :    
285 :     open(BLAST,"<$blast_file") || die "could not open $blast_file";
286 :     while (defined($sim = <BLAST>))
287 :     {
288 :     if (($sim =~ /^(\S+)\t(\S+)(\t\d+){4}\t(\S+)/) && ($1 eq $peg) &&
289 :     ($4 < 1.0e-10))
290 :     {
291 :     $seen{$2} = 1;
292 :     }
293 :     }
294 :     close(BLAST);
295 :     $/ = "\n>";
296 :     }
297 :     }
298 :     }
299 :     $/ = "\n";
300 :     close(FASTA);
301 :     close(REP);
302 :     }
303 :     $fam->{reps} = [map { $_ =~ /^(\S+)/; $1 } `grep -v "^>" $rep_file`];
304 :    
305 : overbeek 1.1 bless $fam,$class;
306 :     return $fam;
307 :     }
308 :    
309 : overbeek 1.11 sub display {
310 :     my($self) = @_;
311 :    
312 :     my $fam_id = $self->family_id;
313 :     my $fig = $self->{fig};
314 :     my $fam_func = $self->family_function;
315 :     my $fam_dir = $self->{dir};
316 :     my @pegs = map { chomp; $_ } `cut -f2 $fam_dir/PEGs`;
317 :    
318 :     my $peg;
319 :     my @disp = ( "ID: $fam_id",
320 :     "Function: $fam_func",
321 :     "Directory: $fam_dir",
322 :     "PEGs:",
323 :     map { $peg = $_; " " . join("\t",($peg,
324 :     $fig->genus_species(&FIG::genome_of($peg)),
325 :     scalar $fig->function_of($_)))
326 :     } @pegs
327 :     );
328 :    
329 :     return join("\n", @disp) . "\n";
330 :     }
331 :    
332 : overbeek 1.1 sub ok_func {
333 : overbeek 1.8 my($fig,$func,$func2,$id2,$fam) = @_;
334 : overbeek 1.1
335 :     my $i;
336 : overbeek 1.8 if (&loose_same_func($func,$func2) || ($func eq $fig->function_of($id2))) { return 1 }
337 :     #
338 :     my $funcI;
339 :     if (defined($funcI = $fam->{ignorable_func}->{"$id2\t$func"}))
340 :     {
341 :     return $funcI;
342 :     }
343 : overbeek 1.7
344 : overbeek 1.1 my @sims = $fig->sims($id2,5,1.0e-30,"fig");
345 : overbeek 1.22 if ($ENV{VERBOSE} && ($ENV{VERBOSE} > 1)) { print STDERR &FIG::flatten_dumper(\@sims), "\n"; }
346 : overbeek 1.8 my($func3);
347 : overbeek 1.7 for ($i=0; ($i < @sims) && (defined($func3 = $fig->function_of($sims[$i]->id2))) &&
348 :     (&FIG::hypo($func3) || &SameFunc::same_func($func,$func3)); $i++) {}
349 : overbeek 1.17 if (($i == @sims) && ($ENV{'VERBOSE'})) { print STDERR "make assignment: $id2\t$func\n" }
350 : overbeek 1.8
351 :     $fam->{ignorable_func}->{"$id2\t$func"} = ($i == @sims);
352 :    
353 :     return $fam->{ignorable_func}->{"$id2\t$func"}
354 : overbeek 1.1 }
355 :    
356 : overbeek 1.8 sub member {
357 :     my($x,$yL) = @_;
358 :     my $i;
359 :    
360 :     for ($i=0; ($i < @$yL) && ($yL->[$i] ne $x); $i++) {}
361 :     return ($i < @$yL);
362 :     }
363 :    
364 : overbeek 1.26 sub list_members {
365 :     my ($self) = @_;
366 :    
367 :     my $fam_dir = $self->{dir};
368 :     my @pegs = map { chomp; $_ } `cut -f2 $fam_dir/PEGs`;
369 :    
370 :     return @pegs;
371 :     }
372 :    
373 : overbeek 1.8 sub loose_same_func {
374 :     my($f1,$f2) = @_;
375 :    
376 : overbeek 1.16 if (&SameFunc::same_func($f1,$f2,'strong')) { return 1 }
377 : overbeek 1.8 my @s1 = split(/\s*[;\/@]\s*/,$f1);
378 :     my @s2 = split(/\s*[;\/@]\s*/,$f2);
379 :     if ((@s1 == 1) && (@s2 > 1) && &member($s1[0],\@s2))
380 :     {
381 :     return 1;
382 :     }
383 :     elsif ((@s2 == 1) && (@s1 > 1) && &member($s2[0],\@s1))
384 :     {
385 :     return 1;
386 :     }
387 :     else
388 :     {
389 :     return 0;
390 :     }
391 :     }
392 :    
393 :    
394 :    
395 : overbeek 1.3 sub representatives {
396 :     my($self) = @_;
397 :    
398 :     my $reps = $self->{reps};
399 :     return @$reps;
400 :     }
401 :    
402 :     sub should_be_member {
403 : overbeek 1.9 my($self,$seq,$debug) = @_;
404 : overbeek 1.19
405 : overbeek 1.9 if ($debug) { $ENV{'VERBOSE'} = $ENV{'DEBUG'} = 1 }
406 : overbeek 1.19
407 : overbeek 1.8 my $fig = $self->{fig};
408 : overbeek 1.3 my $dir = $self->{dir};
409 :     open(TMP,">$FIG_Config::temp/tmp$$.fasta")
410 :     || die "could not open $FIG_Config::temp/tmp$$.fasta";
411 :     print TMP ">query\n$seq\n";
412 :     close(TMP);
413 : overbeek 1.19
414 : overbeek 1.25 if ($ENV{'DEBUG'}) { print STDERR &Dumper([$self->family_id,$seq]) }
415 : overbeek 1.19
416 : overbeek 1.3 my $query_ln = length($seq);
417 : overbeek 1.19
418 :     my @tmp = `blastall -i $FIG_Config::temp/tmp$$.fasta -d $dir/PEGs.fasta -m 8 -FF -p blastp`;
419 :    
420 : overbeek 1.3 my %seen;
421 :     my $should_be = 0;
422 : overbeek 1.7 my $yes = 0;
423 :     my $no = 0;
424 : overbeek 1.19
425 : overbeek 1.8 my $ln1 = length($seq);
426 : overbeek 1.15 my $good = 0;
427 : overbeek 1.19
428 :     my $sims = [];
429 :     foreach $_ (@tmp)
430 : overbeek 1.3 {
431 : overbeek 1.7 if ($ENV{'DEBUG'}) { print STDERR $_ }
432 : overbeek 1.3 chop;
433 : overbeek 1.19
434 : overbeek 1.3 my $sim = [split(/\t/,$_)];
435 :     my $peg = $sim->[1];
436 : overbeek 1.8 next if (($_ = $ENV{'IGNORE'}) && (&FIG::genome_of($peg) eq $_)); # for debugging
437 :    
438 : overbeek 1.3 my $sc = $sim->[10];
439 : overbeek 1.6 my $bit_score = $sim->[11];
440 : overbeek 1.8
441 :     my $matched1 = abs($sim->[7] - $sim->[6]) + 1;
442 :     my $matched2 = abs($sim->[9] - $sim->[8]) + 1;
443 :     my $ln2 = $self->{peg_lengths}->{$peg};
444 :    
445 :     if ((! $seen{$peg}) && ($sc <= 1.0e-10) &&
446 : overbeek 1.24 # ($matched1 >= (0.7) * $ln1) &&
447 :     ($matched2 >= (0.7) * $ln2))
448 : overbeek 1.3 {
449 :     $seen{$peg} = 1;
450 : overbeek 1.4 push @$sim, $query_ln, $self->{peg_lengths}->{$peg};
451 :     bless $sim, 'Sim';
452 :     push @$sims, $sim;
453 : overbeek 1.19
454 : overbeek 1.3 my $bounds = $self->{bounds}->{$peg};
455 : overbeek 1.19 if ($bounds && (@$sims <= 10))
456 : overbeek 1.3 {
457 : overbeek 1.25 if (($bit_score >= $bounds->[1]) && ((! $bounds->[2]) || $bounds->[3] < $bit_score))
458 : overbeek 1.7 {
459 : overbeek 1.8 if ($ENV{'DEBUG'}) { print STDERR " yes - $peg\n" }
460 : overbeek 1.19 ++$yes;
461 : overbeek 1.15 if ($yes > ($no+1)) { $good = 1 }
462 : overbeek 1.7 }
463 :     else
464 :     {
465 : overbeek 1.8 my $bad_func = $bounds->[2] ? $fig->function_of($bounds->[2]) : "";
466 :     if ($ENV{'DEBUG'}) { print STDERR " no - $peg ",join(",",@$bounds)," $bad_func\n" }
467 : overbeek 1.19 ++$no;
468 : overbeek 1.7 }
469 : overbeek 1.3 }
470 :     }
471 :     }
472 : overbeek 1.25 if ($ENV{'DEBUG'}) { print STDERR " yes=$yes no=$no good=$good\n" }
473 : overbeek 1.15 return ($good,$sims);
474 : overbeek 1.3 }
475 :    
476 : overbeek 1.5 sub family_function {
477 :     my($self) = @_;
478 :    
479 :     return $self->{function}
480 :     }
481 :    
482 : overbeek 1.7 sub family_id {
483 :     my($self) = @_;
484 :    
485 :     return $self->{id};
486 :     }
487 :    
488 : overbeek 1.10 sub fam_dir {
489 :     my($fam_data,$fam_id) = @_;
490 :    
491 :     $fam_id =~ /(\d{3})$/;
492 :     return "$fam_data/FIGFAMS/$1/$fam_id";
493 :     }
494 :    
495 : overbeek 1.1 1;

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3