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

Annotation of /FigKernelPackages/FigFam.pm

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1.7 - (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.4
24 :     use Carp;
25 : overbeek 1.1 use Data::Dumper;
26 :    
27 :     # This is the constructor. Presumably, $class is 'FigFam'.
28 :     #
29 :     # $fam_id is the ID of the family (of the form FIGnnnnnnn). It is required.
30 :     #
31 : overbeek 1.3 # $fam_data is optional. If it exists,
32 :     # it gives the directory that is (or will) contain the FigFam data.
33 : overbeek 1.5 # $FIG_Config::data/FigfamsData will be used as a default.
34 : overbeek 1.1 #
35 : overbeek 1.3 # $fam_file must exist within $fam_data. It contains a 2-column version of the FIG families.
36 :     # If the $fam_dir/yyy/FIGxxxxyyy directory (which is where the directory for the specific family
37 : overbeek 1.1 # will go) does not exist, the specification of the PEGs in $fam_file will be used to initialize
38 :     # the family.
39 :     #
40 :    
41 :     sub new {
42 : overbeek 1.3 my($class,$fig,$fam_id,$fam_data) = @_;
43 : overbeek 1.1
44 :     ($fam_id =~ /^FIG\d{3}(\d{3})$/) || die "$fam_id is not a valid FIG family ID";
45 :    
46 :     my $yyy = $1;
47 :     my $fam = {};
48 :     $fam->{id} = $fam_id;
49 :     $fam->{fig} = $fig;
50 :    
51 : overbeek 1.3 if (! defined($fam_data))
52 : overbeek 1.1 {
53 : overbeek 1.3 $fam_data = "$FIG_Config::data/FigfamsData";
54 : overbeek 1.1 }
55 : overbeek 1.3 $fam->{data} = $fam_data;
56 :    
57 :     my $fam_dir = "$fam_data/FIGFAMS";
58 :     my $fam_file = "$fam_data/families.2c";
59 : overbeek 1.1 my $dir = "$fam_dir/$yyy/$fam_id";
60 : overbeek 1.3 $fam->{dir} = $dir;
61 :    
62 : overbeek 1.1 &FIG::verify_dir($dir);
63 :     if ((! -s "$dir/PEGs") && defined($fam_file))
64 :     {
65 : overbeek 1.4 if (-s $fam_file)
66 : overbeek 1.3 {
67 : overbeek 1.4 system("grep $fam_id \"$fam_file\" > \"$dir/PEGs\"");
68 :     if (! -s "$dir/PEGs")
69 :     {
70 :     system "rm -rf $dir";
71 :     confess "PEG file $dir/PEGs does not exist --- $dir deleted";
72 :     return undef;
73 :     }
74 :     }
75 :     else
76 :     {
77 :     confess "Family file $fam_file doe not exist";
78 : overbeek 1.3 return undef;
79 :     }
80 : overbeek 1.4
81 : overbeek 1.1 my @first = `head -n 1 \"$dir/PEGs\"`;
82 :     if ($first[0] =~ /^\S+\t\S+\t(\S.*\S)/)
83 :     {
84 :     open(FUNC,">$dir/function") || die "could not open $dir/function";
85 :     print FUNC "$1\n";
86 :     close(FUNC);
87 :     }
88 :     else
89 :     {
90 :     die "$dir has no valid function file";
91 :     }
92 :     }
93 :     open(FUNC,"<$dir/function") || die "could not open $dir/function";
94 :     my $func = <FUNC>;
95 :     chomp $func;
96 :     close(FUNC);
97 :     $fam->{function} = $func;
98 :    
99 :     my($peg,$pegs);
100 :     my $pegsL = [sort { &FIG::by_fig_id($a,$b) } map { $_ =~ /^\S+\t(\S+)/; $1 } `cat \"$dir/PEGs\"`];
101 :    
102 :     $fam->{pegsL} = $pegsL;
103 :     my $pegsH = {};
104 :     foreach $peg (@$pegsL)
105 :     {
106 :     $pegsH->{$peg} = 1;
107 :     }
108 :     $fam->{pegsH} = $pegsH;
109 : overbeek 1.4
110 :    
111 :     if (-s "$dir/PEGs.fasta")
112 :     {
113 :     open(FASTA,"<$dir/PEGs.fasta") || die "could not read-open $dir/PEGs.fasta";
114 :     while (my ($peg, $seqP) = &FIG::read_fasta_record(\*FASTA))
115 :     {
116 :     $fam->{peg_lengths}->{$peg} = length($$seqP);
117 :     }
118 :     close(FASTA);
119 :     }
120 :     else
121 : overbeek 1.1 {
122 : overbeek 1.4 open(FASTA,">$dir/PEGs.fasta") || die "could not write-open $dir/PEGs.fasta";
123 : overbeek 1.1 foreach $peg (@$pegsL)
124 :     {
125 :     my $seq = $fig->get_translation($peg);
126 :     if ($seq)
127 :     {
128 :     print FASTA ">$peg\n$seq\n";
129 : overbeek 1.3 $fam->{peg_lengths}->{$peg} = length($seq);
130 : overbeek 1.1 }
131 : overbeek 1.4 else
132 :     {
133 :     confess "Could not get seq for $peg";
134 :     }
135 : overbeek 1.1 }
136 :     close(FASTA);
137 : overbeek 1.4 }
138 :    
139 :     if ((! -s "$dir/PEGs.fasta.pin") || ((-M "$dir/PEGs.fasta.pin") > (-M "$dir/PEGs.fasta")))
140 :     {
141 : overbeek 1.3 &FIG::run("$FIG_Config::ext_bin/formatdb -i \"$dir/PEGs.fasta\" -p");
142 : overbeek 1.1 }
143 : overbeek 1.4
144 : overbeek 1.6 if (! -s "$dir/bounds.sims")
145 : overbeek 1.1 {
146 : overbeek 1.6 open(SIMS,">$dir/bounds.sims")
147 :     || die "could not open $dir/bounds.sims";
148 :    
149 : overbeek 1.1 foreach $peg (@$pegsL)
150 :     {
151 :     my @sims = $fig->sims($peg,1000,1.0e-20,"fig");
152 : overbeek 1.6 foreach my $sim (@sims)
153 :     {
154 :     print SIMS join("\t",$sim->id1,
155 :     $sim->id2,
156 :     $sim->bit_score,
157 :     scalar $fig->function_of($sim->id2)),"\n";
158 :     }
159 :     print SIMS "//\n";
160 :     }
161 :     close(SIMS);
162 :     }
163 :    
164 :     if (! -s "$dir/bounds")
165 :     {
166 :     open(BOUNDS,">$dir/bounds") || die "could not open $dir/bounds";
167 :     open(SIMS,"<$dir/bounds.sims") || die "could not open $dir/bounds.sims";
168 : overbeek 1.1
169 : overbeek 1.6 my($sims);
170 :     $/ = "\n//\n";
171 :     while (defined($sims = <SIMS>))
172 :     {
173 :     chomp $sims;
174 :     my @sims = sort { $b->[2] <=> $a->[2] } map { [split(/\t/,$_)] } split(/\n/,$sims);
175 :     if ($peg = $sims[0]->[0])
176 :     {
177 :     my($i,$hits,$bad,$bad_sc,$last_good,$func2);
178 :     for ($i=0, $hits=0, $bad="", $bad_sc=0, $last_good=1; ($i < @sims) && (! $bad); $i++)
179 :     {
180 :     my $sim = $sims[$i];
181 :     my $id2 = $sim->[1];
182 :     if ($pegsH->{$id2})
183 :     {
184 :     $last_good = $sim->[2];
185 :     $hits++;
186 :     }
187 :     else
188 : overbeek 1.1 {
189 : overbeek 1.6 $func2 = $sim->[3];
190 :     if (! &ok_func($fig,$func,$func2,$id2))
191 :     {
192 :     $bad = $id2;
193 :     $bad_sc = $sim->[2];
194 :     }
195 : overbeek 1.1 }
196 :     }
197 :    
198 : overbeek 1.6 if ($hits > 0)
199 :     {
200 :     print BOUNDS join("\t",($peg,$hits,$last_good,$bad,$bad_sc)),"\n";
201 :     }
202 : overbeek 1.1 }
203 :     }
204 :     close(BOUNDS);
205 : overbeek 1.6 close(SIMS);
206 :     $/ = "\n";
207 : overbeek 1.1 }
208 : overbeek 1.6
209 : overbeek 1.1 open(BOUNDS,"<$dir/bounds") || die "could not open $dir/bounds";
210 :    
211 :     my $bounds = {};
212 :     my $x;
213 :     while (defined($x = <BOUNDS>))
214 :     {
215 :     chomp $x;
216 :     my @flds = split(/\t/,$x);
217 :     my $peg =shift @flds;
218 :     $bounds->{$peg} = [@flds];
219 :     }
220 :     close(BOUNDS);
221 :     $fam->{bounds} = $bounds;
222 :    
223 : overbeek 1.3 my $blast_file = "$dir/blastout";
224 :     my $rep_file = "$dir/reps";
225 :     if (! -d "$dir/Split")
226 : overbeek 1.1 {
227 : overbeek 1.3 &FIG::run("split_and_trim_sequences $dir/Split blastout=$blast_file < $dir/PEGs.fasta");
228 :     if (-s $blast_file)
229 :     {
230 :     my @out = map { chop; [split(/\t/,$_) ] } `cat $blast_file`;
231 :     my @keep = ();
232 :     my $write = 0;
233 :     my $i;
234 :     for ($i=0; ($i < @out); $i++)
235 :     {
236 :     if ((($i == 0) ||
237 :     (($out[$i-1]->[0] ne $out[$i]->[0]) || ($out[$i-1]->[1] ne $out[$i]->[1]))) &&
238 :     ($out[$i-1]->[6] <= 1.0e-10))
239 :     {
240 :     push(@keep,$out[$i]);
241 :     }
242 :     else
243 :     {
244 :     $write = 1;
245 :     }
246 :     }
247 :     if ($write)
248 :     {
249 :     open(OUT,">$blast_file") || die $blast_file;
250 :     foreach my $x (@keep)
251 :     {
252 :     print OUT join("\t",@$x),"\n";
253 :     }
254 :     close(OUT);
255 :     }
256 :     }
257 : overbeek 1.1 }
258 : overbeek 1.3
259 :     if ((-s $blast_file) && (! -s $rep_file))
260 :     {
261 :     my $n = 1;
262 :     open(REP,">$rep_file") || die "could not open $rep_file";
263 :     open(FASTA,"<$dir/PEGs.fasta") || die "could not open $dir/PEGs.fasta";
264 :     my %seen;
265 :     my($seq,$peg,$sim);
266 :     $/ = "\n>";
267 :     while (defined($_ = <FASTA>))
268 :     {
269 :     chomp;
270 :     if ($_ =~ /^>?(\S+)[^\n]*\n(.*)/s)
271 :     {
272 :     $peg = $1;
273 :     if (! $seen{$peg})
274 :     {
275 :     $seq = $2;
276 :     $seq =~ s/\s//gs;
277 :     print REP ">$fam_id-$n\n$seq\n";
278 :     $n++;
279 :    
280 :     $/ = "\n";
281 :     $seen{$peg} = 1;
282 :    
283 :     open(BLAST,"<$blast_file") || die "could not open $blast_file";
284 :     while (defined($sim = <BLAST>))
285 :     {
286 :     if (($sim =~ /^(\S+)\t(\S+)(\t\d+){4}\t(\S+)/) && ($1 eq $peg) &&
287 :     ($4 < 1.0e-10))
288 :     {
289 :     $seen{$2} = 1;
290 :     }
291 :     }
292 :     close(BLAST);
293 :     $/ = "\n>";
294 :     }
295 :     }
296 :     }
297 :     $/ = "\n";
298 :     close(FASTA);
299 :     close(REP);
300 :     }
301 :     $fam->{reps} = [map { $_ =~ /^(\S+)/; $1 } `grep -v "^>" $rep_file`];
302 :    
303 : overbeek 1.1 bless $fam,$class;
304 :     return $fam;
305 :     }
306 :    
307 :     sub ok_func {
308 :     my($fig,$func,$func2,$id2) = @_;
309 :    
310 :     my $i;
311 :     if (&SameFunc::same_func($func,$func2)) { return 1 }
312 :     if (! &FIG::hypo($func2)) { return 0 }
313 : overbeek 1.7 return 1;
314 :    
315 : overbeek 1.1 my @sims = $fig->sims($id2,5,1.0e-30,"fig");
316 : overbeek 1.7 my $func3;
317 :     for ($i=0; ($i < @sims) && (defined($func3 = $fig->function_of($sims[$i]->id2))) &&
318 :     (&FIG::hypo($func3) || &SameFunc::same_func($func,$func3)); $i++) {}
319 :     return ($i == @sims);
320 : overbeek 1.1 }
321 :    
322 : overbeek 1.3 sub representatives {
323 :     my($self) = @_;
324 :    
325 :     my $reps = $self->{reps};
326 :     return @$reps;
327 :     }
328 :    
329 :     sub should_be_member {
330 :     my($self,$seq) = @_;
331 : overbeek 1.4
332 : overbeek 1.3 my $dir = $self->{dir};
333 :     open(TMP,">$FIG_Config::temp/tmp$$.fasta")
334 :     || die "could not open $FIG_Config::temp/tmp$$.fasta";
335 :     print TMP ">query\n$seq\n";
336 :     close(TMP);
337 :    
338 : overbeek 1.7 if ($ENV{'DEBUG'}) { print STDERR &Dumper([$self->family_id,$seq]) }
339 :    
340 : overbeek 1.3 my $query_ln = length($seq);
341 : overbeek 1.7 open(BLAST,"blastall -i $FIG_Config::temp/tmp$$.fasta -d $dir/PEGs.fasta -m 8 -FF -p blastp | head -n 10 |")
342 : overbeek 1.3 || die "could not open blast";
343 :     my $sims = [];
344 :     my %seen;
345 :     my $should_be = 0;
346 : overbeek 1.7 my $yes = 0;
347 :     my $no = 0;
348 :    
349 : overbeek 1.3 while (defined($_ = <BLAST>))
350 :     {
351 : overbeek 1.7 if ($ENV{'DEBUG'}) { print STDERR $_ }
352 : overbeek 1.3 chop;
353 :     my $sim = [split(/\t/,$_)];
354 :     my $peg = $sim->[1];
355 :     my $sc = $sim->[10];
356 : overbeek 1.6 my $bit_score = $sim->[11];
357 : overbeek 1.3 if ((! $seen{$peg}) && ($sc <= 1.0e-10))
358 :     {
359 :     $seen{$peg} = 1;
360 : overbeek 1.4 push @$sim, $query_ln, $self->{peg_lengths}->{$peg};
361 :     bless $sim, 'Sim';
362 :     push @$sims, $sim;
363 : overbeek 1.3 my $bounds = $self->{bounds}->{$peg};
364 : overbeek 1.7 if ($bounds)
365 : overbeek 1.3 {
366 : overbeek 1.7 if (($bit_score > $bounds->[1]) && ((! $bounds->[2]) || $bounds->[3] < $bit_score))
367 :     {
368 :     if ($ENV{'DEBUG'}) { print STDERR " yes\n" }
369 :     $yes++;
370 :     }
371 :     else
372 :     {
373 :     if ($ENV{'DEBUG'}) { print STDERR " no\n" }
374 :     $no++;
375 :     }
376 : overbeek 1.3 }
377 :     }
378 :     }
379 : overbeek 1.7 $should_be = ($yes > $no);
380 :     if ($ENV{'DEBUG'}) { print STDERR " should_be = $should_be\n" }
381 : overbeek 1.3 return ($should_be,$sims);
382 :     }
383 :    
384 : overbeek 1.5 sub family_function {
385 :     my($self) = @_;
386 :    
387 :     return $self->{function}
388 :     }
389 :    
390 : overbeek 1.7 sub family_id {
391 :     my($self) = @_;
392 :    
393 :     return $self->{id};
394 :     }
395 :    
396 : overbeek 1.1 1;

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3