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

Annotation of /FigKernelPackages/FigFam.pm

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1.5 - (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 :    
145 : overbeek 1.1 if (! -s "$dir/bounds")
146 :     {
147 :     open(BOUNDS,">$dir/bounds") || die "could not open $dir/bounds";
148 :     foreach $peg (@$pegsL)
149 :     {
150 :     my @sims = $fig->sims($peg,1000,1.0e-20,"fig");
151 :    
152 :     my($i,$hits,$bad,$bad_sc,$last_good,$func2);
153 :     for ($i=0, $hits=0, $bad="", $bad_sc=0, $last_good=1; ($i < @sims) && (! $bad); $i++)
154 :     {
155 :     my $sim = $sims[$i];
156 :     my $id2 = $sim->id2;
157 :     if ($pegsH->{$id2})
158 :     {
159 :     $last_good = $sim->psc;
160 :     $hits++;
161 :     }
162 :     else
163 :     {
164 :     $func2 = $fig->function_of($id2);
165 :     if (! &ok_func($fig,$func,$func2,$id2))
166 :     {
167 :     $bad = $id2;
168 :     $bad_sc = $sim->psc;
169 :     }
170 :     }
171 :     }
172 :    
173 :     if ($hits > 0)
174 :     {
175 :     print BOUNDS join("\t",($peg,$hits,$last_good,$bad,$bad_sc)),"\n";
176 :     }
177 :     }
178 :     close(BOUNDS);
179 :     }
180 :     open(BOUNDS,"<$dir/bounds") || die "could not open $dir/bounds";
181 :    
182 :     my $bounds = {};
183 :     my $x;
184 :     while (defined($x = <BOUNDS>))
185 :     {
186 :     chomp $x;
187 :     my @flds = split(/\t/,$x);
188 :     my $peg =shift @flds;
189 :     $bounds->{$peg} = [@flds];
190 :     }
191 :     close(BOUNDS);
192 :     $fam->{bounds} = $bounds;
193 :    
194 : overbeek 1.3 my $blast_file = "$dir/blastout";
195 :     my $rep_file = "$dir/reps";
196 :     if (! -d "$dir/Split")
197 : overbeek 1.1 {
198 : overbeek 1.3 &FIG::run("split_and_trim_sequences $dir/Split blastout=$blast_file < $dir/PEGs.fasta");
199 :     if (-s $blast_file)
200 :     {
201 :     my @out = map { chop; [split(/\t/,$_) ] } `cat $blast_file`;
202 :     my @keep = ();
203 :     my $write = 0;
204 :     my $i;
205 :     for ($i=0; ($i < @out); $i++)
206 :     {
207 :     if ((($i == 0) ||
208 :     (($out[$i-1]->[0] ne $out[$i]->[0]) || ($out[$i-1]->[1] ne $out[$i]->[1]))) &&
209 :     ($out[$i-1]->[6] <= 1.0e-10))
210 :     {
211 :     push(@keep,$out[$i]);
212 :     }
213 :     else
214 :     {
215 :     $write = 1;
216 :     }
217 :     }
218 :     if ($write)
219 :     {
220 :     open(OUT,">$blast_file") || die $blast_file;
221 :     foreach my $x (@keep)
222 :     {
223 :     print OUT join("\t",@$x),"\n";
224 :     }
225 :     close(OUT);
226 :     }
227 :     }
228 : overbeek 1.1 }
229 : overbeek 1.3
230 :     if ((-s $blast_file) && (! -s $rep_file))
231 :     {
232 :     my $n = 1;
233 :     open(REP,">$rep_file") || die "could not open $rep_file";
234 :     open(FASTA,"<$dir/PEGs.fasta") || die "could not open $dir/PEGs.fasta";
235 :     my %seen;
236 :     my($seq,$peg,$sim);
237 :     $/ = "\n>";
238 :     while (defined($_ = <FASTA>))
239 :     {
240 :     chomp;
241 :     if ($_ =~ /^>?(\S+)[^\n]*\n(.*)/s)
242 :     {
243 :     $peg = $1;
244 :     if (! $seen{$peg})
245 :     {
246 :     $seq = $2;
247 :     $seq =~ s/\s//gs;
248 :     print REP ">$fam_id-$n\n$seq\n";
249 :     $n++;
250 :    
251 :     $/ = "\n";
252 :     $seen{$peg} = 1;
253 :    
254 :     open(BLAST,"<$blast_file") || die "could not open $blast_file";
255 :     while (defined($sim = <BLAST>))
256 :     {
257 :     if (($sim =~ /^(\S+)\t(\S+)(\t\d+){4}\t(\S+)/) && ($1 eq $peg) &&
258 :     ($4 < 1.0e-10))
259 :     {
260 :     $seen{$2} = 1;
261 :     }
262 :     }
263 :     close(BLAST);
264 :     $/ = "\n>";
265 :     }
266 :     }
267 :     }
268 :     $/ = "\n";
269 :     close(FASTA);
270 :     close(REP);
271 :     }
272 :     $fam->{reps} = [map { $_ =~ /^(\S+)/; $1 } `grep -v "^>" $rep_file`];
273 :    
274 : overbeek 1.1 bless $fam,$class;
275 :     return $fam;
276 :     }
277 :    
278 :     sub ok_func {
279 :     my($fig,$func,$func2,$id2) = @_;
280 :    
281 :     my $i;
282 :     if (&SameFunc::same_func($func,$func2)) { return 1 }
283 :     if (! &FIG::hypo($func2)) { return 0 }
284 :     my @sims = $fig->sims($id2,5,1.0e-30,"fig");
285 :     for ($i=0; ($i < @sims) && ($func eq $fig->function_of($sims[$i]->id2)); $i++) {}
286 :     return ($i == 5);
287 :     }
288 :    
289 : overbeek 1.3 sub representatives {
290 :     my($self) = @_;
291 :    
292 :     my $reps = $self->{reps};
293 :     return @$reps;
294 :     }
295 :    
296 :     sub should_be_member {
297 :     my($self,$seq) = @_;
298 : overbeek 1.4
299 : overbeek 1.3 my $dir = $self->{dir};
300 :     open(TMP,">$FIG_Config::temp/tmp$$.fasta")
301 :     || die "could not open $FIG_Config::temp/tmp$$.fasta";
302 :     print TMP ">query\n$seq\n";
303 :     close(TMP);
304 :    
305 :     my $query_ln = length($seq);
306 :     open(BLAST,"blastall -i $FIG_Config::temp/tmp$$.fasta -d $dir/PEGs.fasta -m 8 -FF -p blastp |")
307 :     || die "could not open blast";
308 :     my $sims = [];
309 :     my %seen;
310 :     my $should_be = 0;
311 :     while (defined($_ = <BLAST>))
312 :     {
313 :     chop;
314 :     my $sim = [split(/\t/,$_)];
315 :     my $peg = $sim->[1];
316 :     my $sc = $sim->[10];
317 :     if ((! $seen{$peg}) && ($sc <= 1.0e-10))
318 :     {
319 :     $seen{$peg} = 1;
320 : overbeek 1.4 push @$sim, $query_ln, $self->{peg_lengths}->{$peg};
321 :     bless $sim, 'Sim';
322 :     push @$sims, $sim;
323 : overbeek 1.3 my $bounds = $self->{bounds}->{$peg};
324 : overbeek 1.5 if ($bounds && ($sc <= $bounds->[1]) && ((! $bounds->[2]) || $bounds->[3] > $sc))
325 : overbeek 1.3 {
326 : overbeek 1.5 $should_be = $peg;
327 : overbeek 1.3 }
328 :     }
329 :     }
330 :     return ($should_be,$sims);
331 :     }
332 :    
333 : overbeek 1.5 sub family_function {
334 :     my($self) = @_;
335 :    
336 :     return $self->{function}
337 :     }
338 :    
339 : overbeek 1.1 1;

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3