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

Annotation of /FigKernelPackages/FigFam.pm

Parent Directory Parent Directory | Revision Log Revision Log


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

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3