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

Annotation of /FigKernelPackages/FigFam.pm

Parent Directory Parent Directory | Revision Log Revision Log


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

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3