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

Annotation of /FigKernelPackages/FigFam.pm

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1.38 - (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.32 my $kernel_file = "$dir/kernel_alignment.fasta";
330 : overbeek 1.36 my $ali_file = "$dir/alignment.fasta";
331 : overbeek 1.32
332 :     if ((! -s $kernel_file) && (! -s "$dir/built") && (-s "$dir/PEGs.fasta"))
333 :     {
334 : overbeek 1.36 my $seqs = &gjoseqlib::read_fasta("$dir/PEGs.fasta");
335 : overbeek 1.37 my $ali = &PHOB::trimmed_aligned_kernel(seqs => $seqs, retrim => 1 );
336 : overbeek 1.36 &gjoseqlib::print_alignment_as_fasta("$kernel_file",$ali);
337 :     }
338 : overbeek 1.32
339 : overbeek 1.36 if ((! -s $ali_file) && (-s $kernel_file) && (! -s "$dir/built") && (-s "$dir/PEGs.fasta"))
340 :     {
341 :     my $seqs = &gjoseqlib::read_fasta("$dir/PEGs.fasta");
342 :     my $ali;
343 :     if (-s "$ali_file.last")
344 : overbeek 1.32 {
345 : overbeek 1.36 $ali = &gjoseqlib::read_fasta("$ali_file.last");
346 : overbeek 1.32 }
347 : overbeek 1.36 else
348 :     {
349 :     $ali = &gjoseqlib::read_fasta($kernel_file);
350 :     }
351 :    
352 :     my($seq);
353 :     my %in_ali = map { $_->[0] => 1 } @$ali;
354 :     foreach $seq (grep { ! $in_ali{$_->[0]} } @$seqs)
355 :     {
356 : overbeek 1.38 # print STDERR "adding $seq->[0] to alignment\n";
357 : overbeek 1.36 $ali = &gjoalignment::add_to_alignment($seq,$ali,1);
358 :     &gjoseqlib::print_alignment_as_fasta("$ali_file.last",$ali);
359 :     }
360 : overbeek 1.38
361 :     if (! -s "$ali_file.last")
362 :     {
363 :     &gjoseqlib::print_alignment_as_fasta($ali_file,$ali);
364 :     }
365 :     else
366 :     {
367 :     rename("$ali_file.last",$ali_file);
368 :     }
369 : overbeek 1.32 }
370 :    
371 : overbeek 1.28 if (! -s "$dir/built") { system "echo 1 > $dir/built" }
372 : overbeek 1.1 bless $fam,$class;
373 :     return $fam;
374 :     }
375 :    
376 : overbeek 1.31
377 :     =head3 reset_function
378 :    
379 :     usage:
380 :     C<< $figfam_obj->reset_function(); >>
381 :    
382 :     Resets the function of a FIGfam to its "consensus function,"
383 :     where the weight of a vote toward the "consensus" is doubled
384 :     for PEGs connected to a subsytem.
385 :    
386 :     =cut
387 :    
388 : overbeek 1.28 sub reset_function {
389 :     my($self) = @_;
390 :    
391 :     my $fam_dir = $self->{dir};
392 :     my $fig = $self->{fig};
393 :    
394 :     my @pegs = grep { $fig->is_real_feature($_) } map { $_ =~ /^\S+\t(\S+)/; $1 } `cat "$fam_dir/PEGs"`;
395 :     my($peg,%funcs,$func,@subs);
396 :     foreach $peg (@pegs)
397 :     {
398 :     $func = $fig->function_of($peg);
399 :     @subs = grep { $fig->usable_subsystem($_) } $fig->peg_to_subsystems($peg);
400 :     my $incr = (@subs > 0) ? 2 : 1;
401 :     $funcs{$func} += $incr;
402 :     }
403 :     my @funcs = sort { $funcs{$b} <=> $funcs{$a} } keys(%funcs);
404 :     $func = (@funcs > 0) ? $funcs[0] : "hypothetical protein";
405 :     open(FUNC,">$fam_dir/function") || die "could not open $fam_dir/function";
406 : overbeek 1.29 print FUNC "$func\n";
407 : overbeek 1.28 close(FUNC);
408 :     }
409 :    
410 : overbeek 1.31
411 :     =head3 display
412 :    
413 :     usage:
414 :     C<< print $figfam_obj->display(); >>
415 :    
416 :     Returns a list-formatted table describing a family and its members.
417 :    
418 :     =cut
419 :    
420 : overbeek 1.11 sub display {
421 :     my($self) = @_;
422 :    
423 :     my $fam_id = $self->family_id;
424 :     my $fig = $self->{fig};
425 :     my $fam_func = $self->family_function;
426 :     my $fam_dir = $self->{dir};
427 :     my @pegs = map { chomp; $_ } `cut -f2 $fam_dir/PEGs`;
428 :    
429 :     my $peg;
430 :     my @disp = ( "ID: $fam_id",
431 :     "Function: $fam_func",
432 :     "Directory: $fam_dir",
433 :     "PEGs:",
434 :     map { $peg = $_; " " . join("\t",($peg,
435 :     $fig->genus_species(&FIG::genome_of($peg)),
436 :     scalar $fig->function_of($_)))
437 :     } @pegs
438 :     );
439 :    
440 :     return join("\n", @disp) . "\n";
441 :     }
442 :    
443 : overbeek 1.31
444 :     =head3 fasta_of
445 :    
446 :     usage:
447 :    
448 :     C<< $seqs_of_members = $figfam_obj->fasta_of(); >>
449 :    
450 :     or
451 :    
452 :     C<< $figfam_obj->fasta_of($filehandle_glob); >>
453 :    
454 :     Returns a ptr to a hash containing the sequences for all members of a FIGfam,
455 :     keyed by FIG feature ID (FID).
456 :     The filehandle glob is optional;
457 :     if included, the sequences will also print to the handle's file,
458 :     in fasta format.
459 :    
460 :     =cut
461 :    
462 : overbeek 1.30 sub fasta_of {
463 :     my($self, $fh) = @_;
464 :    
465 :     my $fam_id = $self->family_id;
466 :     my $fam_dir = $self->{dir};
467 :    
468 :     if (-s "$fam_dir/PEGs.fasta") {
469 :     open(FASTA, "<$fam_dir/PEGs.fasta")
470 :     || die "Could not read-open $fam_dir/PEGs.fasta";
471 :    
472 :     my $seq_of = {};
473 :     my ($id, $seqP);
474 :     while (($id, $seqP) = &FIG::read_fasta_record(\*FASTA)) {
475 :     $seq_of->{$id} = $$seqP;
476 :     if (defined($fh)) {
477 :     print $fh ">$id\n$$seqP\n";
478 :     }
479 :     }
480 :     close(FASTA) || die "Could not close $fam_dir/PEGs.fasta";
481 :     return $seq_of;
482 :     }
483 :     else {
484 :     return undef;
485 :     }
486 :     }
487 :    
488 : overbeek 1.31
489 :    
490 :     =head3 ok_func
491 :    
492 :     usage:
493 :     C<< if ( &FigFam::ok_func( $fig, $func, $func2, $id2, $fam) ) { #...do something... } >>
494 :    
495 :     C<$fig> is an object constriucted by C<FIG->new()>.
496 :    
497 :     C<$func> and C<$func2> are two possible assigned functions.
498 :    
499 :     C<$id2> is the FIG ID (FID) of the PEG having function C<$func2>.
500 :    
501 :     C<$fam> is the family that C<$id2> is presumed to be a member of.
502 :    
503 :     Returns C<TRUE> is C<$func> and C$func2> are loosely "the same,"
504 :     or if C<$func> is judged to be "ignorable."
505 :    
506 :     =cut
507 :    
508 : overbeek 1.1 sub ok_func {
509 : overbeek 1.8 my($fig,$func,$func2,$id2,$fam) = @_;
510 : overbeek 1.1
511 :     my $i;
512 : overbeek 1.8 if (&loose_same_func($func,$func2) || ($func eq $fig->function_of($id2))) { return 1 }
513 :     #
514 :     my $funcI;
515 :     if (defined($funcI = $fam->{ignorable_func}->{"$id2\t$func"}))
516 :     {
517 :     return $funcI;
518 :     }
519 : overbeek 1.7
520 : overbeek 1.1 my @sims = $fig->sims($id2,5,1.0e-30,"fig");
521 : overbeek 1.22 if ($ENV{VERBOSE} && ($ENV{VERBOSE} > 1)) { print STDERR &FIG::flatten_dumper(\@sims), "\n"; }
522 : overbeek 1.8 my($func3);
523 : overbeek 1.7 for ($i=0; ($i < @sims) && (defined($func3 = $fig->function_of($sims[$i]->id2))) &&
524 :     (&FIG::hypo($func3) || &SameFunc::same_func($func,$func3)); $i++) {}
525 : overbeek 1.17 if (($i == @sims) && ($ENV{'VERBOSE'})) { print STDERR "make assignment: $id2\t$func\n" }
526 : overbeek 1.8
527 :     $fam->{ignorable_func}->{"$id2\t$func"} = ($i == @sims);
528 :    
529 :     return $fam->{ignorable_func}->{"$id2\t$func"}
530 : overbeek 1.1 }
531 :    
532 : overbeek 1.31
533 :     =head3 member
534 :    
535 :     usage:
536 :     C<< if ( &FigFam::member($x, \@y)) { #...do something... } >>
537 :    
538 :     C<$x> is a scalar.
539 :    
540 :     C<\@y> is a ptr to a list.
541 :    
542 :     Returns C<TRUE> if C<$x> is in list C<@y>.
543 :    
544 :     =cut
545 :    
546 : overbeek 1.8 sub member {
547 :     my($x,$yL) = @_;
548 :     my $i;
549 :    
550 :     for ($i=0; ($i < @$yL) && ($yL->[$i] ne $x); $i++) {}
551 :     return ($i < @$yL);
552 :     }
553 :    
554 : overbeek 1.31
555 :    
556 :     =head3 list_members
557 :    
558 :     usage:
559 :     C<< @ids = $figfam_obj->list_members(); >>
560 :    
561 :     Returns a list of the PEG FIDs in a family.
562 :    
563 :     =cut
564 :    
565 : overbeek 1.26 sub list_members {
566 :     my ($self) = @_;
567 :    
568 :     my $fam_dir = $self->{dir};
569 :     my @pegs = map { chomp; $_ } `cut -f2 $fam_dir/PEGs`;
570 :    
571 :     return @pegs;
572 :     }
573 :    
574 : overbeek 1.31
575 :    
576 :     =head3 loose_same_func
577 :    
578 :     usage:
579 :    
580 :     C<< if ( &FigFam::loose_same_func($f1, $f2) ) { #...do something... } >>
581 :    
582 :     C<$f1> and C<$f2> are function strings.
583 :    
584 :     Returns C<TRUE> if C<$f1> and C<$f2> are "more or less the same."
585 :    
586 :     =cut
587 :    
588 : overbeek 1.8 sub loose_same_func {
589 :     my($f1,$f2) = @_;
590 :    
591 : overbeek 1.16 if (&SameFunc::same_func($f1,$f2,'strong')) { return 1 }
592 : overbeek 1.8 my @s1 = split(/\s*[;\/@]\s*/,$f1);
593 :     my @s2 = split(/\s*[;\/@]\s*/,$f2);
594 :     if ((@s1 == 1) && (@s2 > 1) && &member($s1[0],\@s2))
595 :     {
596 :     return 1;
597 :     }
598 :     elsif ((@s2 == 1) && (@s1 > 1) && &member($s2[0],\@s1))
599 :     {
600 :     return 1;
601 :     }
602 :     else
603 :     {
604 :     return 0;
605 :     }
606 :     }
607 :    
608 :    
609 :    
610 : overbeek 1.31 =head3 representatives
611 :    
612 :     usage:
613 :     C<< @rep_seqs = $figfam_obj->representatives(); >>
614 :    
615 :     Returns a list of the "representative sequences" characterizing a FIGfam.
616 :    
617 :     =cut
618 :    
619 : overbeek 1.3 sub representatives {
620 :     my($self) = @_;
621 :    
622 :     my $reps = $self->{reps};
623 :     return @$reps;
624 :     }
625 :    
626 : overbeek 1.31
627 :    
628 :     =head3 should_be_member
629 :    
630 :     usage:
631 :     C<< if ( $figfam_obj->should_be_member( $seq ) ) { #...do something... } >>
632 :    
633 :     Returns C<TRUE> if the protein sequence in C<$seq> is judged to be
634 :     "similar enough" to the members of a family to potentially be included.
635 :    
636 :     =cut
637 :    
638 : overbeek 1.3 sub should_be_member {
639 : overbeek 1.9 my($self,$seq,$debug) = @_;
640 : overbeek 1.19
641 : overbeek 1.9 if ($debug) { $ENV{'VERBOSE'} = $ENV{'DEBUG'} = 1 }
642 : overbeek 1.19
643 : overbeek 1.8 my $fig = $self->{fig};
644 : overbeek 1.3 my $dir = $self->{dir};
645 :     open(TMP,">$FIG_Config::temp/tmp$$.fasta")
646 :     || die "could not open $FIG_Config::temp/tmp$$.fasta";
647 :     print TMP ">query\n$seq\n";
648 :     close(TMP);
649 : overbeek 1.19
650 : overbeek 1.25 if ($ENV{'DEBUG'}) { print STDERR &Dumper([$self->family_id,$seq]) }
651 : overbeek 1.19
652 : overbeek 1.3 my $query_ln = length($seq);
653 : overbeek 1.19
654 :     my @tmp = `blastall -i $FIG_Config::temp/tmp$$.fasta -d $dir/PEGs.fasta -m 8 -FF -p blastp`;
655 :    
656 : overbeek 1.3 my %seen;
657 :     my $should_be = 0;
658 : overbeek 1.7 my $yes = 0;
659 :     my $no = 0;
660 : overbeek 1.19
661 : overbeek 1.8 my $ln1 = length($seq);
662 : overbeek 1.15 my $good = 0;
663 : overbeek 1.19
664 :     my $sims = [];
665 :     foreach $_ (@tmp)
666 : overbeek 1.3 {
667 : overbeek 1.7 if ($ENV{'DEBUG'}) { print STDERR $_ }
668 : overbeek 1.3 chop;
669 : overbeek 1.19
670 : overbeek 1.3 my $sim = [split(/\t/,$_)];
671 :     my $peg = $sim->[1];
672 : overbeek 1.8 next if (($_ = $ENV{'IGNORE'}) && (&FIG::genome_of($peg) eq $_)); # for debugging
673 :    
674 : overbeek 1.3 my $sc = $sim->[10];
675 : overbeek 1.6 my $bit_score = $sim->[11];
676 : overbeek 1.8
677 :     my $matched1 = abs($sim->[7] - $sim->[6]) + 1;
678 :     my $matched2 = abs($sim->[9] - $sim->[8]) + 1;
679 :     my $ln2 = $self->{peg_lengths}->{$peg};
680 :    
681 :     if ((! $seen{$peg}) && ($sc <= 1.0e-10) &&
682 : overbeek 1.24 # ($matched1 >= (0.7) * $ln1) &&
683 :     ($matched2 >= (0.7) * $ln2))
684 : overbeek 1.3 {
685 :     $seen{$peg} = 1;
686 : overbeek 1.4 push @$sim, $query_ln, $self->{peg_lengths}->{$peg};
687 :     bless $sim, 'Sim';
688 :     push @$sims, $sim;
689 : overbeek 1.19
690 : overbeek 1.3 my $bounds = $self->{bounds}->{$peg};
691 : overbeek 1.19 if ($bounds && (@$sims <= 10))
692 : overbeek 1.3 {
693 : overbeek 1.25 if (($bit_score >= $bounds->[1]) && ((! $bounds->[2]) || $bounds->[3] < $bit_score))
694 : overbeek 1.7 {
695 : overbeek 1.8 if ($ENV{'DEBUG'}) { print STDERR " yes - $peg\n" }
696 : overbeek 1.19 ++$yes;
697 : overbeek 1.15 if ($yes > ($no+1)) { $good = 1 }
698 : overbeek 1.7 }
699 :     else
700 :     {
701 : overbeek 1.8 my $bad_func = $bounds->[2] ? $fig->function_of($bounds->[2]) : "";
702 :     if ($ENV{'DEBUG'}) { print STDERR " no - $peg ",join(",",@$bounds)," $bad_func\n" }
703 : overbeek 1.19 ++$no;
704 : overbeek 1.7 }
705 : overbeek 1.3 }
706 :     }
707 :     }
708 : overbeek 1.25 if ($ENV{'DEBUG'}) { print STDERR " yes=$yes no=$no good=$good\n" }
709 : overbeek 1.15 return ($good,$sims);
710 : overbeek 1.3 }
711 :    
712 : overbeek 1.31
713 :    
714 :     =head3 family_function
715 :    
716 :     usage:
717 :     C<< $func = $figfam_obj->family_function(); >>
718 :    
719 :     Returns the "consensus function" assigned to a FIGfam object.
720 :    
721 :     =cut
722 :    
723 : overbeek 1.5 sub family_function {
724 :     my($self) = @_;
725 :    
726 :     return $self->{function}
727 :     }
728 :    
729 : overbeek 1.31
730 :    
731 :     =head3 family_id
732 :    
733 :     usage:
734 :     C<< $fam_id = $figfam_obj->family_id(); >>
735 :    
736 :     Returns the FIGfam ID of a FIGfam object.
737 :    
738 :     =cut
739 :    
740 : overbeek 1.7 sub family_id {
741 :     my($self) = @_;
742 :    
743 :     return $self->{id};
744 :     }
745 :    
746 : overbeek 1.31
747 :    
748 :     =head3 family_dir
749 :    
750 :     usage:
751 :    
752 :     C<< $dir = &FigFam::family_dir( $fam_data, $fam_id ); >>
753 :    
754 :     Returns the path to the subdirectory of C<$fam_data>
755 :     that the FIGfam data for a FIGfam with ID C<$fam_id>
756 :     would be stored in.
757 :    
758 :     =cut
759 :    
760 : overbeek 1.10 sub fam_dir {
761 :     my($fam_data,$fam_id) = @_;
762 :    
763 :     $fam_id =~ /(\d{3})$/;
764 :     return "$fam_data/FIGFAMS/$1/$fam_id";
765 :     }
766 :    
767 : overbeek 1.1 1;

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3