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

Annotation of /FigKernelPackages/FigFam.pm

Parent Directory Parent Directory | Revision Log Revision Log


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

1 : gdpusch 1.52 # -*- perl -*-
2 :     ########################################################################
3 : overbeek 1.1 # Copyright (c) 2003-2006 University of Chicago and Fellowship
4 :     # for Interpretations of Genomes. All Rights Reserved.
5 :     #
6 :     # This file is part of the SEED Toolkit.
7 :     #
8 :     # The SEED Toolkit is free software. You can redistribute
9 :     # it and/or modify it under the terms of the SEED Toolkit
10 :     # Public License.
11 :     #
12 :     # You should have received a copy of the SEED Toolkit Public License
13 :     # along with this program; if not write to the University of Chicago
14 :     # at info@ci.uchicago.edu or the Fellowship for Interpretation of
15 :     # Genomes at veronika@thefig.info or download a copy from
16 :     # http://www.theseed.org/LICENSE.TXT.
17 : gdpusch 1.52 ########################################################################
18 : overbeek 1.1
19 :     package FigFam;
20 :    
21 :     use strict;
22 :     use FIG;
23 :     use SameFunc;
24 : overbeek 1.32 use gjoseqlib;
25 :     use gjoalignment;
26 :     use PHOB;
27 : overbeek 1.45 use Tracer;
28 : overbeek 1.4
29 :     use Carp;
30 : overbeek 1.1 use Data::Dumper;
31 :    
32 : overbeek 1.31 =head1 Module to access FIGfams
33 :    
34 :     =head3 new
35 :    
36 :     usage:
37 :     C<< my $figfam_obj = FigFam->new($fig, $fam_id, $fam_dir); >>
38 :    
39 :     This is the constructor. Presumably, C<$class> is 'FigFam'.
40 :    
41 :     C<$fig> is an object reference created by a C<< FIG->new() >> constructor.
42 :    
43 :     C<$fam_id> is the ID of the family, of the form C<FIGnnnnnn> where C<n> is a digit;
44 :     it is required.
45 :    
46 :     C<$fam_data> is optional. If it is defined, it will be treated
47 :     as the directory that contains (or will contain) FigFam data.
48 :     If omitted, the FIGfam directory defaults to
49 :     C<$FIG_Config::data/FigfamsData>.
50 :    
51 :     A families file C<families.2c> must exist within C<$fam_data.>
52 :     It contains a 3-column tab-separated "tuple representation" of the FIG families,
53 :     of the form
54 :     C<<< FIGfam_IDE<lt>tabE<gt>FIG_Protein_IDE<lt>tabE<gt>Assigned_function. >>>
55 :    
56 :     If the C<$fam_dir/yyy/FIGxxxyyy> directory (which is where the subdirectory
57 :     where the specific family will go) does not exist, the specification of the PEGs
58 :     in C<$fam_file> will be used to initialize the family.
59 :    
60 :     =cut
61 : overbeek 1.1
62 : overbeek 1.50
63 :     #### To turn on tracing, at the command line type
64 :     ###
65 :     ### export TRACING=Ross
66 :     ### trace 3 FigFam
67 :     ###
68 :    
69 : overbeek 1.1 sub new {
70 : overbeek 1.3 my($class,$fig,$fam_id,$fam_data) = @_;
71 : overbeek 1.1
72 : overbeek 1.25 ($fam_id =~ /^FIG\d{6}$/) || confess "invalid family id: $fam_id";
73 :    
74 : overbeek 1.1 my $fam = {};
75 :     $fam->{id} = $fam_id;
76 :     $fam->{fig} = $fig;
77 :    
78 : overbeek 1.3 if (! defined($fam_data))
79 : overbeek 1.1 {
80 : overbeek 1.3 $fam_data = "$FIG_Config::data/FigfamsData";
81 : overbeek 1.1 }
82 : overbeek 1.3 $fam->{data} = $fam_data;
83 :    
84 :     my $fam_dir = "$fam_data/FIGFAMS";
85 :     my $fam_file = "$fam_data/families.2c";
86 : overbeek 1.10 my $dir = &fam_dir($fam_data,$fam_id);
87 : overbeek 1.3 $fam->{dir} = $dir;
88 :    
89 : overbeek 1.1 &FIG::verify_dir($dir);
90 : overbeek 1.50
91 : overbeek 1.45 if ((! -s "$dir/built") && (! -s "$dir/PEGs") && defined($fam_file))
92 : overbeek 1.1 {
93 : overbeek 1.4 if (-s $fam_file)
94 : overbeek 1.3 {
95 : overbeek 1.4 system("grep $fam_id \"$fam_file\" > \"$dir/PEGs\"");
96 :     if (! -s "$dir/PEGs")
97 :     {
98 :     system "rm -rf $dir";
99 :     confess "PEG file $dir/PEGs does not exist --- $dir deleted";
100 :     return undef;
101 :     }
102 : overbeek 1.50 &Tracer::Trace('created PEGs file') if (T(3));
103 : overbeek 1.4 }
104 :     else
105 :     {
106 : overbeek 1.13 confess "Family file $fam_file does not exist";
107 : overbeek 1.3 return undef;
108 :     }
109 : overbeek 1.18 }
110 :    
111 : overbeek 1.45 if ((! -s "$dir/function") && (! -s "$dir/built"))
112 : overbeek 1.18 {
113 : overbeek 1.28 &reset_function($fam);
114 : overbeek 1.50 &Tracer::Trace("reset function for $fam_id") if (T(3));
115 : overbeek 1.1 }
116 : overbeek 1.18
117 : overbeek 1.1 open(FUNC,"<$dir/function") || die "could not open $dir/function";
118 :     my $func = <FUNC>;
119 :     chomp $func;
120 :     close(FUNC);
121 :     $fam->{function} = $func;
122 :    
123 :     my($peg,$pegs);
124 : overbeek 1.32 my $pegsL = [
125 :     sort { &FIG::by_fig_id($a,$b) }
126 : overbeek 1.59 # grep { scalar $fig->function_of($_) eq $func }
127 : overbeek 1.32 grep { $fig->is_real_feature($_) }
128 :     map { $_ =~ /^\S+\t(\S+)/; $1 }
129 :     `cat \"$dir/PEGs\"`
130 :     ];
131 :    
132 : overbeek 1.46 if (@$pegsL < 2)
133 :     {
134 : overbeek 1.53 if (-w $dir) {
135 : gdpusch 1.52 open( TMP, ">$dir/built") || die "Could not write-open $dir/built";
136 :     print TMP "1\n";
137 :     close(TMP) || die "Could not close $dir/built";
138 :     }
139 :     else {
140 : olson 1.58 # warn "WARNING: $dir/built is not writable";
141 : gdpusch 1.52 }
142 : overbeek 1.47 return undef;
143 : overbeek 1.46 }
144 :    
145 : overbeek 1.1 $fam->{pegsL} = $pegsL;
146 :     my $pegsH = {};
147 :     foreach $peg (@$pegsL)
148 :     {
149 :     $pegsH->{$peg} = 1;
150 :     }
151 :     $fam->{pegsH} = $pegsH;
152 : overbeek 1.4
153 :    
154 :     if (-s "$dir/PEGs.fasta")
155 :     {
156 :     open(FASTA,"<$dir/PEGs.fasta") || die "could not read-open $dir/PEGs.fasta";
157 :     while (my ($peg, $seqP) = &FIG::read_fasta_record(\*FASTA))
158 :     {
159 :     $fam->{peg_lengths}->{$peg} = length($$seqP);
160 :     }
161 :     close(FASTA);
162 :     }
163 : overbeek 1.28 elsif (! -s "$dir/built")
164 : overbeek 1.1 {
165 : overbeek 1.4 open(FASTA,">$dir/PEGs.fasta") || die "could not write-open $dir/PEGs.fasta";
166 : overbeek 1.1 foreach $peg (@$pegsL)
167 :     {
168 :     my $seq = $fig->get_translation($peg);
169 :     if ($seq)
170 :     {
171 :     print FASTA ">$peg\n$seq\n";
172 : overbeek 1.3 $fam->{peg_lengths}->{$peg} = length($seq);
173 : overbeek 1.1 }
174 : overbeek 1.4 else
175 :     {
176 :     confess "Could not get seq for $peg";
177 :     }
178 : overbeek 1.1 }
179 :     close(FASTA);
180 : overbeek 1.4 }
181 : overbeek 1.45
182 : overbeek 1.28 if ((! -s "$dir/built") && ((! -s "$dir/PEGs.fasta.pin") || ((-M "$dir/PEGs.fasta.pin") > (-M "$dir/PEGs.fasta"))))
183 : overbeek 1.4 {
184 : overbeek 1.3 &FIG::run("$FIG_Config::ext_bin/formatdb -i \"$dir/PEGs.fasta\" -p");
185 : overbeek 1.1 }
186 : overbeek 1.4
187 : overbeek 1.28 if ((! -s "$dir/bounds.sims") && (! -s "$dir/built"))
188 : overbeek 1.1 {
189 : overbeek 1.6 open(SIMS,">$dir/bounds.sims")
190 :     || die "could not open $dir/bounds.sims";
191 : overbeek 1.53
192 :     my($i,$n,$j,@sims,$req);
193 :     for ($i=0; ($i < @$pegsL); $i += 50)
194 : overbeek 1.1 {
195 : overbeek 1.53 $n = ((@$pegsL - $i) > 50) ? $i+49 : (@$pegsL - 1);
196 :     $req = [@$pegsL[$i..$n]];
197 :     @sims = sort { ($a->id1 cmp $b->id1) } $fig->sims($req,50000,1.0e-10,"fig");
198 :     if (@sims == 0) { print STDERR "no similarities returned for ",join(",",@$req),"\n"; }
199 :    
200 :     if ($ENV{'DEBUG'})
201 :     {
202 :     print STDERR "extracting similarities\n";
203 :     print STDERR &Dumper($req);
204 :     foreach $_ (@sims) { print join("\t",($_->id1,$_->id2,$_->psc)),"\n";}
205 :     }
206 :    
207 :     for ($j=0; ($j < @sims); $j++)
208 : overbeek 1.6 {
209 : overbeek 1.53 my $sim = $sims[$j];
210 : overbeek 1.6 print SIMS join("\t",$sim->id1,
211 :     $sim->id2,
212 :     $sim->bit_score,
213 :     scalar $fig->function_of($sim->id2)),"\n";
214 : overbeek 1.53 if (($j == $#sims) || ($sim->id1 ne $sims[$j+1]->id1))
215 :     {
216 :     print SIMS "//\n";
217 :     }
218 : overbeek 1.6 }
219 :     }
220 :     close(SIMS);
221 :     }
222 : overbeek 1.32
223 : overbeek 1.45 if ((! -s "$dir/built") && (! -s "$dir/bounds"))
224 : overbeek 1.6 {
225 : overbeek 1.8 print STDERR " building bounds for $fam_id\n";
226 : overbeek 1.6 open(BOUNDS,">$dir/bounds") || die "could not open $dir/bounds";
227 :     open(SIMS,"<$dir/bounds.sims") || die "could not open $dir/bounds.sims";
228 : overbeek 1.1
229 : overbeek 1.6 my($sims);
230 :     $/ = "\n//\n";
231 :     while (defined($sims = <SIMS>))
232 :     {
233 :     chomp $sims;
234 :     my @sims = sort { $b->[2] <=> $a->[2] } map { [split(/\t/,$_)] } split(/\n/,$sims);
235 :     if ($peg = $sims[0]->[0])
236 :     {
237 : overbeek 1.53 # print STDERR "processing $peg\n";
238 : overbeek 1.6 my($i,$hits,$bad,$bad_sc,$last_good,$func2);
239 :     for ($i=0, $hits=0, $bad="", $bad_sc=0, $last_good=1; ($i < @sims) && (! $bad); $i++)
240 :     {
241 :     my $sim = $sims[$i];
242 :     my $id2 = $sim->[1];
243 :     if ($pegsH->{$id2})
244 :     {
245 :     $last_good = $sim->[2];
246 :     $hits++;
247 :     }
248 :     else
249 : overbeek 1.1 {
250 : overbeek 1.6 $func2 = $sim->[3];
251 : overbeek 1.8 if (! &ok_func($fig,$func,$func2,$id2,$fam))
252 : overbeek 1.6 {
253 :     $bad = $id2;
254 :     $bad_sc = $sim->[2];
255 :     }
256 : overbeek 1.1 }
257 :     }
258 :    
259 : overbeek 1.6 if ($hits > 0)
260 :     {
261 :     print BOUNDS join("\t",($peg,$hits,$last_good,$bad,$bad_sc)),"\n";
262 :     }
263 : overbeek 1.53 else
264 :     {
265 :     # print STDERR " failed\n";
266 :     }
267 : overbeek 1.1 }
268 :     }
269 :     close(BOUNDS);
270 : overbeek 1.6 close(SIMS);
271 :     $/ = "\n";
272 : overbeek 1.1 }
273 : overbeek 1.6
274 : overbeek 1.1 my $bounds = {};
275 : overbeek 1.45 if (open(BOUNDS,"<$dir/bounds"))
276 : overbeek 1.1 {
277 : overbeek 1.45 my $x;
278 :     while (defined($x = <BOUNDS>))
279 :     {
280 :     chomp $x;
281 :     my @flds = split(/\t/,$x);
282 : gdpusch 1.54 my $peg = shift @flds;
283 : overbeek 1.45 $bounds->{$peg} = [@flds];
284 :     }
285 :     close(BOUNDS);
286 :     $fam->{bounds} = $bounds;
287 : overbeek 1.1 }
288 : gdpusch 1.54
289 : overbeek 1.50 if ((! -d "$dir/SplitNew") && (! -s "$dir/built"))
290 : overbeek 1.1 {
291 : overbeek 1.50 mkdir("$dir/SplitNew",0777) || die "could not make $dir/SplitNew";
292 :     my $seqs = &gjoseqlib::read_fasta("$dir/PEGs.fasta");
293 :     my %seqs = map { $_->[0] => $_ } @$seqs;
294 :     my(undef,$sets) = &representative_sequences::rep_seq_2($seqs,{ max_sim => 0.3 });
295 :     my $n = 1;
296 :     foreach my $id (sort { @{$sets->{$b}} <=> @{$sets->{$a}} } keys(%$sets))
297 : overbeek 1.3 {
298 : overbeek 1.50 my @ext_seqs = map { $seqs{$_} } ($id,@{$sets->{$id}});
299 :     if (@ext_seqs > 1)
300 : overbeek 1.3 {
301 : overbeek 1.50 my $trimmed = &PHOB::kernel_of_trimmed_seqs(seqs => \@ext_seqs);
302 :     if (@$trimmed > 1)
303 : overbeek 1.3 {
304 : overbeek 1.50 my $ali = &PHOB::align_seqs(seqs => $trimmed);
305 :     &gjoseqlib::print_alignment_as_fasta("$dir/SplitNew/$n.ali",$ali);
306 :     $n++;
307 : overbeek 1.3 }
308 :     }
309 :     }
310 : overbeek 1.1 }
311 : overbeek 1.3
312 : overbeek 1.50
313 :     my $blast_file = "$dir/blastout";
314 :     my $rep_file = "$dir/reps";
315 :     if ((! -f $blast_file) && (! -s "$dir/built"))
316 :     {
317 :     &FIG::run("split_and_trim_sequences $dir/Split blastout=$blast_file < $dir/PEGs.fasta");
318 :     if (-s $blast_file)
319 :     {
320 :     my @out = map { chop; [split(/\t/,$_) ] } `cat $blast_file`;
321 :     my @keep = ();
322 :     my $write = 0;
323 :     my $i;
324 :     for ($i=0; ($i < @out); $i++)
325 :     {
326 :     if ((($i == 0) ||
327 :     (($out[$i-1]->[0] ne $out[$i]->[0]) || ($out[$i-1]->[1] ne $out[$i]->[1]))) &&
328 :     ($out[$i-1]->[6] <= 1.0e-10))
329 :     {
330 :     push(@keep,$out[$i]);
331 :     }
332 :     else
333 :     {
334 :     $write = 1;
335 :     }
336 :     }
337 :     if ($write)
338 :     {
339 :     open(OUT,">$blast_file") || die $blast_file;
340 :     foreach my $x (@keep)
341 :     {
342 :     print OUT join("\t",@$x),"\n";
343 :     }
344 :     close(OUT);
345 :     }
346 :     }
347 :     }
348 :    
349 : overbeek 1.51 # if ((-s "$dir/Split/set.sizes") && (! -s "$dir/built"))
350 :     # {
351 :     # my @to_align = map { (($_ =~ /^(\d+)\t(\d+)/) && ($2 > 1)) ? $1 : () } `cat $dir/Split/set.sizes`;
352 :     # foreach my $n (@to_align)
353 :     # {
354 :     # if (! -s "$dir/Split/$n.aln")
355 :     # {
356 :     # system "runclustalw $dir/Split/$n";
357 :     # }
358 :     # }
359 :     # }
360 : overbeek 1.45
361 : overbeek 1.28 if ((-s $blast_file) && (! -s $rep_file) && (! -s "$dir/built"))
362 : overbeek 1.3 {
363 : overbeek 1.50 my $n = 1;
364 :     open(REP,">$rep_file") || die "could not open $rep_file";
365 :     open(FASTA,"<$dir/PEGs.fasta") || die "could not open $dir/PEGs.fasta";
366 :     my %seen;
367 :     my($seq,$peg,$sim);
368 :     $/ = "\n>";
369 :     while (defined($_ = <FASTA>))
370 :     {
371 :     chomp;
372 :     if ($_ =~ /^>?(\S+)[^\n]*\n(.*)/s)
373 :     {
374 :     $peg = $1;
375 :     if (! $seen{$peg})
376 :     {
377 :     $seq = $2;
378 :     $seq =~ s/\s//gs;
379 :     print REP ">$fam_id-$n\n$seq\n";
380 :     $n++;
381 :    
382 :     $/ = "\n";
383 :     $seen{$peg} = 1;
384 :    
385 :     open(BLAST,"<$blast_file") || die "could not open $blast_file";
386 :     while (defined($sim = <BLAST>))
387 :     {
388 :     if (($sim =~ /^(\S+)\t(\S+)(\t\d+){4}\t(\S+)/) && ($1 eq $peg) &&
389 :     ($4 < 1.0e-10))
390 :     {
391 :     $seen{$2} = 1;
392 :     }
393 :     }
394 :     close(BLAST);
395 :     $/ = "\n>";
396 :     }
397 :     }
398 :     }
399 :     $/ = "\n";
400 :     close(FASTA);
401 :     close(REP);
402 :     }
403 : overbeek 1.3
404 : overbeek 1.28 $fam->{reps} = (-s $rep_file) ? [map { $_ =~ /^(\S+)/; $1 } `grep -v "^>" $rep_file`] : [];
405 : overbeek 1.3
406 : overbeek 1.45 # if (@$pegsL < 4)
407 :     # {
408 :     # if (! -s "$dir/built") { system "echo 1 > $dir/built" }
409 :     # return undef;
410 :     # }
411 :    
412 : gdpusch 1.52 if ((! -s "$dir/built") && $fam->{bounds}) {
413 : overbeek 1.53 if (-w $dir) {
414 : gdpusch 1.52 open( TMP, ">$dir/built") || die "Could not write-open $dir/built";
415 :     print TMP "1\n";
416 :     close(TMP) || die "Could not close $dir/built";
417 :     }
418 :     else {
419 :     warn "WARNING: $dir/built is not writable";
420 :     }
421 :     }
422 : overbeek 1.39 bless $fam,$class;
423 :     return $fam;
424 :     }
425 :    
426 :     sub alignment {
427 :     my($self) = @_;
428 :    
429 :     my $dir = $self->{dir};
430 :     my $fig = $self->{fig};
431 :    
432 : overbeek 1.32 my $kernel_file = "$dir/kernel_alignment.fasta";
433 : overbeek 1.36 my $ali_file = "$dir/alignment.fasta";
434 : overbeek 1.32
435 : overbeek 1.39 if ((! -s $kernel_file) && (! -s "$dir/built.ali") && (-s "$dir/PEGs.fasta"))
436 : overbeek 1.32 {
437 : overbeek 1.36 my $seqs = &gjoseqlib::read_fasta("$dir/PEGs.fasta");
438 : overbeek 1.42 my $ali = &PHOB::trimmed_aligned_kernel(seqs => $seqs, retrim => 1, tmp => $FIG_Config::temp );
439 : overbeek 1.36 &gjoseqlib::print_alignment_as_fasta("$kernel_file",$ali);
440 :     }
441 : overbeek 1.32
442 : overbeek 1.39 if ((! -s $ali_file) && (-s $kernel_file) && (! -s "$dir/built.ali") && (-s "$dir/PEGs.fasta"))
443 : overbeek 1.36 {
444 :     my $seqs = &gjoseqlib::read_fasta("$dir/PEGs.fasta");
445 :     my $ali;
446 :     if (-s "$ali_file.last")
447 : overbeek 1.32 {
448 : overbeek 1.36 $ali = &gjoseqlib::read_fasta("$ali_file.last");
449 : overbeek 1.32 }
450 : overbeek 1.36 else
451 :     {
452 :     $ali = &gjoseqlib::read_fasta($kernel_file);
453 :     }
454 :    
455 :     my($seq);
456 :     my %in_ali = map { $_->[0] => 1 } @$ali;
457 :     foreach $seq (grep { ! $in_ali{$_->[0]} } @$seqs)
458 :     {
459 : overbeek 1.39 print STDERR "adding $seq->[0] to alignment\n";
460 : overbeek 1.36 $ali = &gjoalignment::add_to_alignment($seq,$ali,1);
461 : overbeek 1.39 my @tmp = @$ali;
462 : overbeek 1.36 &gjoseqlib::print_alignment_as_fasta("$ali_file.last",$ali);
463 : overbeek 1.39 system "cat $ali_file.last >> $dir/log; echo '====' >> $dir/log";
464 : overbeek 1.36 }
465 : overbeek 1.38
466 :     if (! -s "$ali_file.last")
467 :     {
468 :     &gjoseqlib::print_alignment_as_fasta($ali_file,$ali);
469 :     }
470 :     else
471 :     {
472 :     rename("$ali_file.last",$ali_file);
473 :     }
474 : overbeek 1.39 system "echo 1 > $dir/built.ali";
475 :     }
476 :     if (-s $ali_file)
477 :     {
478 :     return &gjoseqlib::read_fasta($ali_file);
479 :     }
480 :     else
481 :     {
482 :     return undef;
483 : overbeek 1.32 }
484 : overbeek 1.1 }
485 :    
486 : overbeek 1.31
487 :     =head3 reset_function
488 :    
489 :     usage:
490 :     C<< $figfam_obj->reset_function(); >>
491 :    
492 :     Resets the function of a FIGfam to its "consensus function,"
493 :     where the weight of a vote toward the "consensus" is doubled
494 :     for PEGs connected to a subsytem.
495 :    
496 :     =cut
497 :    
498 : overbeek 1.28 sub reset_function {
499 :     my($self) = @_;
500 :    
501 :     my $fam_dir = $self->{dir};
502 :     my $fig = $self->{fig};
503 :    
504 :     my @pegs = grep { $fig->is_real_feature($_) } map { $_ =~ /^\S+\t(\S+)/; $1 } `cat "$fam_dir/PEGs"`;
505 :     my($peg,%funcs,$func,@subs);
506 :     foreach $peg (@pegs)
507 :     {
508 :     $func = $fig->function_of($peg);
509 :     @subs = grep { $fig->usable_subsystem($_) } $fig->peg_to_subsystems($peg);
510 :     my $incr = (@subs > 0) ? 2 : 1;
511 :     $funcs{$func} += $incr;
512 :     }
513 :     my @funcs = sort { $funcs{$b} <=> $funcs{$a} } keys(%funcs);
514 :     $func = (@funcs > 0) ? $funcs[0] : "hypothetical protein";
515 :     open(FUNC,">$fam_dir/function") || die "could not open $fam_dir/function";
516 : overbeek 1.29 print FUNC "$func\n";
517 : overbeek 1.28 close(FUNC);
518 :     }
519 :    
520 : overbeek 1.31
521 :     =head3 display
522 :    
523 :     usage:
524 :     C<< print $figfam_obj->display(); >>
525 :    
526 :     Returns a list-formatted table describing a family and its members.
527 :    
528 :     =cut
529 :    
530 : overbeek 1.11 sub display {
531 :     my($self) = @_;
532 :    
533 :     my $fam_id = $self->family_id;
534 :     my $fig = $self->{fig};
535 :     my $fam_func = $self->family_function;
536 :     my $fam_dir = $self->{dir};
537 :     my @pegs = map { chomp; $_ } `cut -f2 $fam_dir/PEGs`;
538 :    
539 :     my $peg;
540 :     my @disp = ( "ID: $fam_id",
541 :     "Function: $fam_func",
542 :     "Directory: $fam_dir",
543 :     "PEGs:",
544 :     map { $peg = $_; " " . join("\t",($peg,
545 :     $fig->genus_species(&FIG::genome_of($peg)),
546 :     scalar $fig->function_of($_)))
547 :     } @pegs
548 :     );
549 :    
550 :     return join("\n", @disp) . "\n";
551 :     }
552 :    
553 : mkubal 1.40 =head3 pegs_of
554 :    
555 :     usage:
556 :     C<< print $figfam_obj->pegs_of(); >>
557 :    
558 :     Returns a list of just pegs.
559 :    
560 :     =cut
561 :    
562 :     sub pegs_of {
563 :     my($self) = @_;
564 :     my $fam_dir = $self->{dir};
565 :     my @pegs = map { chomp; $_ } `cut -f2 $fam_dir/PEGs`;
566 :     return [@pegs];
567 :     }
568 :    
569 : overbeek 1.31
570 :     =head3 fasta_of
571 :    
572 :     usage:
573 :    
574 :     C<< $seqs_of_members = $figfam_obj->fasta_of(); >>
575 :    
576 :     or
577 :    
578 :     C<< $figfam_obj->fasta_of($filehandle_glob); >>
579 :    
580 :     Returns a ptr to a hash containing the sequences for all members of a FIGfam,
581 :     keyed by FIG feature ID (FID).
582 :     The filehandle glob is optional;
583 :     if included, the sequences will also print to the handle's file,
584 :     in fasta format.
585 :    
586 :     =cut
587 :    
588 : overbeek 1.30 sub fasta_of {
589 :     my($self, $fh) = @_;
590 :    
591 :     my $fam_id = $self->family_id;
592 :     my $fam_dir = $self->{dir};
593 :    
594 :     if (-s "$fam_dir/PEGs.fasta") {
595 :     open(FASTA, "<$fam_dir/PEGs.fasta")
596 :     || die "Could not read-open $fam_dir/PEGs.fasta";
597 :    
598 :     my $seq_of = {};
599 :     my ($id, $seqP);
600 :     while (($id, $seqP) = &FIG::read_fasta_record(\*FASTA)) {
601 :     $seq_of->{$id} = $$seqP;
602 :     if (defined($fh)) {
603 :     print $fh ">$id\n$$seqP\n";
604 :     }
605 :     }
606 :     close(FASTA) || die "Could not close $fam_dir/PEGs.fasta";
607 :     return $seq_of;
608 :     }
609 :     else {
610 :     return undef;
611 :     }
612 :     }
613 :    
614 : overbeek 1.31
615 :    
616 :     =head3 ok_func
617 :    
618 :     usage:
619 :     C<< if ( &FigFam::ok_func( $fig, $func, $func2, $id2, $fam) ) { #...do something... } >>
620 :    
621 :     C<$fig> is an object constriucted by C<FIG->new()>.
622 :    
623 :     C<$func> and C<$func2> are two possible assigned functions.
624 :    
625 :     C<$id2> is the FIG ID (FID) of the PEG having function C<$func2>.
626 :    
627 :     C<$fam> is the family that C<$id2> is presumed to be a member of.
628 :    
629 :     Returns C<TRUE> is C<$func> and C$func2> are loosely "the same,"
630 :     or if C<$func> is judged to be "ignorable."
631 :    
632 :     =cut
633 :    
634 : overbeek 1.1 sub ok_func {
635 : overbeek 1.8 my($fig,$func,$func2,$id2,$fam) = @_;
636 : overbeek 1.1
637 :     my $i;
638 : overbeek 1.51 if ($func eq $fig->function_of($id2)) { return 1 }
639 :     my @roles = split(/(\s*;\s+)|( [\@\/] )/,$func2);
640 :     for ($i=0; ($i < @roles) && (! &in_sub($fig,$roles[$i],$fam)); $i++) {}
641 :     if ($i < @roles) { return 0 }
642 :     if (&loose_same_func($func,$func2)) { return 1 }
643 : overbeek 1.8 #
644 :     my $funcI;
645 :     if (defined($funcI = $fam->{ignorable_func}->{"$id2\t$func"}))
646 :     {
647 :     return $funcI;
648 :     }
649 : overbeek 1.7
650 : overbeek 1.1 my @sims = $fig->sims($id2,5,1.0e-30,"fig");
651 : overbeek 1.22 if ($ENV{VERBOSE} && ($ENV{VERBOSE} > 1)) { print STDERR &FIG::flatten_dumper(\@sims), "\n"; }
652 : overbeek 1.8 my($func3);
653 : overbeek 1.7 for ($i=0; ($i < @sims) && (defined($func3 = $fig->function_of($sims[$i]->id2))) &&
654 :     (&FIG::hypo($func3) || &SameFunc::same_func($func,$func3)); $i++) {}
655 : gdpusch 1.52 if (($i == @sims) && ($ENV{VERBOSE})) { print STDERR "make assignment: $id2\t$func\n" }
656 : overbeek 1.8
657 :     $fam->{ignorable_func}->{"$id2\t$func"} = ($i == @sims);
658 :    
659 :     return $fam->{ignorable_func}->{"$id2\t$func"}
660 : overbeek 1.1 }
661 :    
662 : overbeek 1.51 sub in_sub {
663 :     my($fig,$role,$fam) = @_;
664 :    
665 :     my $x;
666 :     if (! defined($x = $fam->{in_sub}->{$role}))
667 :     {
668 :     my @sub = grep { $fig->usable_subsystem($_) } $fig->role_to_subsystems($role);
669 :     $x = $fam->{in_sub}->{$role} = (@sub > 0);
670 :     }
671 :     return $x;
672 :     }
673 : overbeek 1.31
674 :     =head3 member
675 :    
676 :     usage:
677 :     C<< if ( &FigFam::member($x, \@y)) { #...do something... } >>
678 :    
679 :     C<$x> is a scalar.
680 :    
681 :     C<\@y> is a ptr to a list.
682 :    
683 :     Returns C<TRUE> if C<$x> is in list C<@y>.
684 :    
685 :     =cut
686 :    
687 : overbeek 1.8 sub member {
688 :     my($x,$yL) = @_;
689 :     my $i;
690 :    
691 :     for ($i=0; ($i < @$yL) && ($yL->[$i] ne $x); $i++) {}
692 :     return ($i < @$yL);
693 :     }
694 :    
695 : overbeek 1.31
696 :    
697 :     =head3 list_members
698 :    
699 :     usage:
700 :     C<< @ids = $figfam_obj->list_members(); >>
701 :    
702 :     Returns a list of the PEG FIDs in a family.
703 :    
704 :     =cut
705 :    
706 : overbeek 1.26 sub list_members {
707 :     my ($self) = @_;
708 :    
709 : gdpusch 1.56 my $fig = $self->{fig};
710 : overbeek 1.26 my $fam_dir = $self->{dir};
711 : gdpusch 1.56 my @pegs = map { chomp;
712 :     (not $fig->is_deleted_fid($_)) ? ($_) : ()
713 :     } `cut -f2 $fam_dir/PEGs`;
714 : overbeek 1.26
715 :     return @pegs;
716 :     }
717 :    
718 : overbeek 1.31
719 :    
720 :     =head3 loose_same_func
721 :    
722 :     usage:
723 :    
724 :     C<< if ( &FigFam::loose_same_func($f1, $f2) ) { #...do something... } >>
725 :    
726 :     C<$f1> and C<$f2> are function strings.
727 :    
728 :     Returns C<TRUE> if C<$f1> and C<$f2> are "more or less the same."
729 :    
730 :     =cut
731 :    
732 : overbeek 1.8 sub loose_same_func {
733 :     my($f1,$f2) = @_;
734 :    
735 : overbeek 1.16 if (&SameFunc::same_func($f1,$f2,'strong')) { return 1 }
736 : overbeek 1.8 my @s1 = split(/\s*[;\/@]\s*/,$f1);
737 :     my @s2 = split(/\s*[;\/@]\s*/,$f2);
738 :     if ((@s1 == 1) && (@s2 > 1) && &member($s1[0],\@s2))
739 :     {
740 :     return 1;
741 :     }
742 :     elsif ((@s2 == 1) && (@s1 > 1) && &member($s2[0],\@s1))
743 :     {
744 :     return 1;
745 :     }
746 :     else
747 :     {
748 :     return 0;
749 :     }
750 :     }
751 :    
752 :    
753 :    
754 : overbeek 1.31 =head3 representatives
755 :    
756 :     usage:
757 :     C<< @rep_seqs = $figfam_obj->representatives(); >>
758 :    
759 :     Returns a list of the "representative sequences" characterizing a FIGfam.
760 :    
761 :     =cut
762 :    
763 : overbeek 1.3 sub representatives {
764 :     my($self) = @_;
765 :    
766 :     my $reps = $self->{reps};
767 :     return @$reps;
768 :     }
769 :    
770 : overbeek 1.31
771 :    
772 :     =head3 should_be_member
773 :    
774 :     usage:
775 :     C<< if ( $figfam_obj->should_be_member( $seq ) ) { #...do something... } >>
776 :    
777 :     Returns C<TRUE> if the protein sequence in C<$seq> is judged to be
778 :     "similar enough" to the members of a family to potentially be included.
779 :    
780 : overbeek 1.57 I have added the "loose" argument as an optional last argument. This means that
781 :    
782 :     C<< if ( $figfam_obj->should_be_member( $seq,0,1 ) ) { #...do something... } >>
783 :    
784 :     will return true, even if the input sequence is truncated (i.e., we do not force the
785 :     similarity to go across 80% of matched sequences in the family).
786 :    
787 : overbeek 1.31 =cut
788 :    
789 : overbeek 1.3 sub should_be_member {
790 : overbeek 1.57 my($self,$seq,$debug,$loose) = @_;
791 : overbeek 1.19
792 : gdpusch 1.52 if ($debug) {
793 :     $ENV{DEBUG} = 1;
794 :     $ENV{VERBOSE} = $ENV{VERBOSE} ? $ENV{VERBOSE} : 1;
795 :     print STDERR "\nTesting for membership:\n";
796 :     }
797 : overbeek 1.19
798 : overbeek 1.8 my $fig = $self->{fig};
799 : overbeek 1.3 my $dir = $self->{dir};
800 :     open(TMP,">$FIG_Config::temp/tmp$$.fasta")
801 :     || die "could not open $FIG_Config::temp/tmp$$.fasta";
802 :     print TMP ">query\n$seq\n";
803 :     close(TMP);
804 : overbeek 1.19
805 : gdpusch 1.52 if ($ENV{DEBUG}) { print STDERR &Dumper([$self->family_id,$seq]) }
806 : overbeek 1.19
807 : overbeek 1.3 my $query_ln = length($seq);
808 : overbeek 1.19
809 :     my @tmp = `blastall -i $FIG_Config::temp/tmp$$.fasta -d $dir/PEGs.fasta -m 8 -FF -p blastp`;
810 :    
811 : overbeek 1.3 my %seen;
812 :     my $should_be = 0;
813 : overbeek 1.7 my $yes = 0;
814 :     my $no = 0;
815 : overbeek 1.19
816 : overbeek 1.8 my $ln1 = length($seq);
817 : overbeek 1.15 my $good = 0;
818 : overbeek 1.42 my $bad = 0;
819 :    
820 : overbeek 1.19 my $sims = [];
821 : overbeek 1.55 for (my $simI=0; ($simI < @tmp) && (! $good) && (! $bad); $simI++)
822 : overbeek 1.3 {
823 : overbeek 1.55 $_ = $tmp[$simI];
824 : gdpusch 1.52 if ($ENV{DEBUG}) { print STDERR $_ }
825 : overbeek 1.3 chop;
826 : overbeek 1.19
827 : overbeek 1.3 my $sim = [split(/\t/,$_)];
828 :     my $peg = $sim->[1];
829 : gdpusch 1.52 next if (($_ = $ENV{IGNORE}) && (&FIG::genome_of($peg) eq $_)); # for debugging
830 : overbeek 1.8
831 : overbeek 1.3 my $sc = $sim->[10];
832 : overbeek 1.6 my $bit_score = $sim->[11];
833 : overbeek 1.8
834 :     my $matched1 = abs($sim->[7] - $sim->[6]) + 1;
835 :     my $matched2 = abs($sim->[9] - $sim->[8]) + 1;
836 :     my $ln2 = $self->{peg_lengths}->{$peg};
837 :    
838 : overbeek 1.57 if ((! $seen{$peg}) &&
839 :     ($loose ||
840 :     (($sc <= 1.0e-10) &&
841 :     ($matched1 >= (0.8 * $ln1)) &&
842 :     ($matched2 >= (0.8 * $ln2)))
843 :     )
844 :     )
845 : overbeek 1.3 {
846 :     $seen{$peg} = 1;
847 : overbeek 1.4 push @$sim, $query_ln, $self->{peg_lengths}->{$peg};
848 :     bless $sim, 'Sim';
849 :     push @$sims, $sim;
850 : overbeek 1.19
851 : overbeek 1.3 my $bounds = $self->{bounds}->{$peg};
852 : overbeek 1.19 if ($bounds && (@$sims <= 10))
853 : overbeek 1.3 {
854 : overbeek 1.57 if ((($bit_score >= $bounds->[1]) && ((! $bounds->[2]) || $bounds->[3] < $bit_score)) ||
855 :     ($loose &&
856 :     ((($bit_score/$ln1) >= ($bounds->[1] / $ln2)) &&
857 :     ((! $bounds->[2]) || (($bounds->[3]/$ln2) < ($bit_score/$ln1))))))
858 : overbeek 1.7 {
859 : gdpusch 1.52 if ($ENV{DEBUG}) { print STDERR " yes - $peg matched1=$matched1 ln1=$ln1 matched2=$matched2 ln2=$ln2\n" }
860 : overbeek 1.19 ++$yes;
861 : overbeek 1.15 if ($yes > ($no+1)) { $good = 1 }
862 : overbeek 1.7 }
863 :     else
864 :     {
865 : overbeek 1.8 my $bad_func = $bounds->[2] ? $fig->function_of($bounds->[2]) : "";
866 : gdpusch 1.52 if ($ENV{DEBUG}) { print STDERR " no - $peg ", join(",",@$bounds)," $bad_func\n" }
867 : overbeek 1.19 ++$no;
868 : overbeek 1.42 if ($no > $yes) { $bad = 1 }
869 : overbeek 1.7 }
870 : overbeek 1.3 }
871 : gdpusch 1.52 else {
872 :     if ($ENV{DEBUG}) {
873 :     print STDERR " bounds=", ($bounds ? qq(non-nil) : qq(nil))
874 :     , ", num_sims=", (scalar @$sims), "\n";
875 :     }
876 :     }
877 :     }
878 :     else {
879 :     if ($ENV{DEBUG}) {
880 :     print STDERR
881 :     " seen=", ($seen{$peg} ? $seen{$peg} : 0), " score=$sc,"
882 :     , " matched1=$matched1, ln1=$ln1,"
883 :     , " matched2=$matched2, ln2=$ln2,"
884 :     , "\n";
885 :     }
886 : overbeek 1.3 }
887 :     }
888 : gdpusch 1.52 if ($ENV{DEBUG}) { print STDERR " yes=$yes no=$no good=$good\n"; }
889 : overbeek 1.15 return ($good,$sims);
890 : overbeek 1.3 }
891 :    
892 : overbeek 1.31 =head3 family_function
893 :    
894 :     usage:
895 :     C<< $func = $figfam_obj->family_function(); >>
896 :    
897 :     Returns the "consensus function" assigned to a FIGfam object.
898 :    
899 :     =cut
900 :    
901 : overbeek 1.5 sub family_function {
902 :     my($self) = @_;
903 :    
904 :     return $self->{function}
905 :     }
906 :    
907 : overbeek 1.31
908 :    
909 :     =head3 family_id
910 :    
911 :     usage:
912 :     C<< $fam_id = $figfam_obj->family_id(); >>
913 :    
914 :     Returns the FIGfam ID of a FIGfam object.
915 :    
916 :     =cut
917 :    
918 : overbeek 1.7 sub family_id {
919 :     my($self) = @_;
920 :    
921 :     return $self->{id};
922 :     }
923 :    
924 : overbeek 1.31
925 :    
926 :     =head3 family_dir
927 :    
928 :     usage:
929 :    
930 :     C<< $dir = &FigFam::family_dir( $fam_data, $fam_id ); >>
931 :    
932 :     Returns the path to the subdirectory of C<$fam_data>
933 :     that the FIGfam data for a FIGfam with ID C<$fam_id>
934 :     would be stored in.
935 :    
936 :     =cut
937 :    
938 : overbeek 1.10 sub fam_dir {
939 :     my($fam_data,$fam_id) = @_;
940 :    
941 :     $fam_id =~ /(\d{3})$/;
942 :     return "$fam_data/FIGFAMS/$1/$fam_id";
943 :     }
944 :    
945 : overbeek 1.1 1;

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3