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

Annotation of /FigKernelPackages/FigFam.pm

Parent Directory Parent Directory | Revision Log Revision Log


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

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3