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

Annotation of /FigKernelPackages/FigFam.pm

Parent Directory Parent Directory | Revision Log Revision Log


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

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3