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

Annotation of /FigKernelPackages/FigFam.pm

Parent Directory Parent Directory | Revision Log Revision Log


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

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3