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

Annotation of /FigKernelPackages/SimFC.pm

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1.3 - (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 SimFC;
19 :    
20 :     use DB_File;
21 :    
22 :     use FIG;
23 :     use Data::Dumper;
24 :    
25 :     sub generate_partitions {
26 :     my($fig,$max_sz) = @_;
27 :     my $partitionsF = "$FIG_Config::temp/partitions.$$";
28 :     open(PART,">$partitionsF") || die "could not open $partitionsF";
29 :     my %prot_hash;
30 :     my %seen_hash;
31 :     &build_hashes($fig,\%prot_hash,\%seen_hash,$max_sz);
32 :    
33 :     my @all = sort keys(%prot_hash);
34 :     my $n = @all;
35 :     print STDERR "$n proteins to process\n";
36 :    
37 :     my $i;
38 :     for ($i=0; ($i < @all); $i++)
39 :     {
40 :     if (! $seen_hash{$all[$i]})
41 :     {
42 :     print STDERR "processing $all[$i], $i of $#all\n";
43 :     &process_one_peg($fig,$all[$i],\%seen_hash,\%prot_hash,$max_sz,\$n,\*PART);
44 :     }
45 :     else
46 :     {
47 :     print STDERR "seen $all[$i]\n";
48 :     }
49 :     }
50 :     close(PART);
51 :     return $partitionsF;
52 :     }
53 :    
54 :     sub process_one_peg {
55 :     my($fig,$prot,$seen,$prot_hash,$max_sz,$nP,$fh) = @_;
56 :    
57 :     my %in;
58 :     $in{$prot} = 1;
59 :     my %closest;
60 :     $closest{$prot} = 0;
61 :    
62 :     my $sz = 1;
63 :     my $peg = $prot;
64 :     while (($sz < $max_sz) && defined($peg))
65 :     {
66 :     $$nP--;
67 :     delete $closest{$peg};
68 :     my @sims = $fig->sims($peg,10000,1,'raw');
69 :     foreach my $sim (@sims)
70 :     {
71 :     my $id2 = $sim->id2;
72 :     if ($prot_hash->{$id2})
73 :     {
74 :     if (! $in{$id2})
75 :     {
76 :     $in{$id2} = 1;
77 :     $sz++;
78 :     $closest{$id2} = $sim->bsc;
79 :     }
80 :     elsif (defined($closest{$id2}) && ($closest{$id2} < $sim->bsc))
81 :     {
82 :     $closest{$id2} = $sim->bsc;
83 :     }
84 :     }
85 :     }
86 :     $seen->{$peg} = 1;
87 :     print $fh "IN\t$prot\t$peg\n";
88 :     $peg = &the_closest(\%closest,$seen);
89 :     }
90 :     print STDERR "$$nP left\n";
91 :     foreach $_ (sort keys(%in))
92 :     {
93 :     print $fh join("\t",('SET',$prot,$_)),"\n";
94 :     }
95 :     }
96 :    
97 :     sub the_closest {
98 :     my($closest,$seen) = @_;
99 :    
100 :     my($peg,$x,$peg1,$x1);
101 :     while (($peg1,$x1) = each %$closest)
102 :     {
103 :     if ((! $seen->{$peg1}) && ((! defined($peg1)) || ($x1 > $x)))
104 :     {
105 :     $peg = $peg1;
106 :     $x = $x1;
107 :     }
108 :     }
109 :     return $peg;
110 :     }
111 :    
112 :     sub build_hashes {
113 :     my($fig,$prot_hash,$seen_hash,$sz) = @_;
114 :    
115 :     my($prot_hash_tie,$seen_hash_tie);
116 :     # $prot_hash_tie = tie %$prot_hash, 'DB_File',"partition_prot_hash_$sz.db",O_RDWR,0666,$DB_BTREE;
117 :     # $seen_hash_tie = tie %$seen_hash, 'DB_File',"seen_hash_$sz.db",O_RDWR,0666,$DB_BTREE;
118 :     # return;
119 :    
120 :     $prot_hash_tie = tie %$prot_hash, 'DB_File',"partition_prot_hash_$sz.db",O_CREAT,0666,$DB_BTREE;
121 :     $prot_hash_tie || die "tie failed";
122 :    
123 :     $seen_hash_tie = tie %$seen_hash, 'DB_File',"seen_hash_$sz.db",O_CREAT,0666,$DB_BTREE;
124 :     $seen_hash_tie || die "tie failed";
125 :    
126 :     open(SYMS,"<$FIG_Config::global/peg.synonyms")
127 :     || die "could not open peg.synonyms";
128 :     while (defined($_ = <SYMS>))
129 :     {
130 :     chop;
131 :     my($head,$rest) = split(/\t/,$_);
132 :     if ($rest =~ /fig\|/)
133 :     {
134 :     my @fig = map { ($_ =~ /(fig\|\d+\.\d+\.peg\.\d+)/) ? $1 : () } split(/;/,$rest);
135 :     @fig = grep { $fig->is_complete(&FIG::genome_of($_)) } @fig;
136 :     if (@fig > 0)
137 :     {
138 :     $head =~ s/,\d+$//;
139 :     $prot_hash->{$head} = 1;
140 :     foreach my $peg (@fig)
141 :     {
142 :     $seen_hash->{$peg} = 1;
143 :     }
144 :     }
145 :     }
146 :     }
147 :     foreach my $genome ($fig->genomes('complete'))
148 :     {
149 :     foreach my $peg ($fig->all_features($genome,'peg'))
150 :     {
151 :     if (! $seen_hash->{$peg})
152 :     {
153 :     $prot_hash->{$peg} = 1;
154 :     }
155 :     }
156 :     }
157 :     }
158 :    
159 :     sub build_partition_directories {
160 :     my($fig,$partitionsF,$dir) = @_;
161 :    
162 :     if (-d $dir) { die "$dir already exists" }
163 :    
164 :     mkdir($dir,0777) || die "could not make $dir";
165 :    
166 :     my $partD = "$dir/Partitions";
167 :     my $simF = "$dir/SimilarityPartitions";
168 :     my $relF = "$dir/ToSimilarityPartition";
169 :    
170 :     my %prot_hash;
171 :     &build_expand_hash($fig,\%prot_hash);
172 :    
173 :     open(SIMPART,">$simF")
174 :     || die "could not open $simF";
175 :     open(REL,">$relF")
176 :     || die "could not open $relF";
177 :    
178 :     mkdir("$partD",0777) || die "could not make $partD";
179 :     my $n = 1;
180 :     open(PART,"<$partitionsF") || die "could not open $partitionsF";
181 :     my $x = <PART>;
182 :     while (defined($x) && ($x =~ /^IN\t(\S+)\t(\S+)/))
183 :     {
184 :     my $currset = $1;
185 :     my @setI = ();
186 :     while (defined($x) && ($x =~ /^IN\t(\S+)\t(\S+)/) && ($1 eq $currset))
187 :     {
188 :     push(@setI,$2);
189 :     $x = <PART>;
190 :     }
191 :    
192 :     my $bug;
193 :     my @setS = ();
194 :     while (defined($x) && ($x =~ /^SET\t(\S+)\t(\S+)/) && ($1 eq $currset))
195 :     {
196 :     # if ($1 eq 'fig|314291.3.peg.2784') { $bug = 1 }
197 :     push(@setS,$2);
198 :     $x = <PART>;
199 :     }
200 :     &process_set($fig,\*SIMPART,\*REL,$partD,\@setI,\@setS,\$n,\%prot_hash,$bug);
201 :     }
202 :     close(SIMPART);
203 :     close(REL);
204 :     close(PART);
205 :    
206 :     print STDERR "starting to build representative sets for each partition: using 10 processors\n";
207 : overbeek 1.2 system "make_partition_reps $dir 10 2> $dir/make_partition_reps.stderr > $dir/make_partition_reps.stdout";
208 : overbeek 1.1 }
209 :    
210 :     sub process_set {
211 :     my($fig,$entity_fh,$rel_fh,$partD,$setI,$setS,$nP,$prot_hash,$bug) = @_;
212 :     # if ($bug) { print STDERR &Dumper($setI,$setS) }
213 :    
214 :     my @reduced = &reduce($fig,$setS,$prot_hash);
215 :     # if ($bug) { print STDERR &Dumper('reduced1',\@reduced); }
216 :     if (@reduced > 1)
217 :     {
218 :     my $md5;
219 :     foreach $md5 (@reduced)
220 :     {
221 :     print $entity_fh join("\t",($$nP,$md5)),"\n";
222 :     }
223 :     my @reduced_in = &reduce($fig,$setI,$prot_hash);
224 :     # if ($bug) { print STDERR &Dumper('reduced_in',\@reduced_in); }
225 :     foreach $md5 (@reduced_in)
226 :     {
227 :     print $rel_fh join("\t",($$nP,$md5)),"\n";
228 :     }
229 :     my $dir1 = $$nP % 1000;
230 :     &FIG::verify_dir("$partD/$dir1");
231 :     my $fasta = "$partD/$dir1/$$nP";
232 :     open(FASTA,">$fasta") || die "could not make $partD/$dir1/$$nP";
233 :     open(SUMMARY,">$fasta.summary") || die "could not open $fasta.summary";
234 :     my %locs;
235 : overbeek 1.3 my $count = 0;
236 : overbeek 1.1
237 :     foreach $md5 (@reduced)
238 :     {
239 : overbeek 1.3 my @pegs = grep { $fig->is_real_feature($_) } $fig->pegs_with_md5($md5);
240 : overbeek 1.1 # if ($bug) { print STDERR &Dumper([$md5,\@pegs]); }
241 :     if (@pegs == 0)
242 :     {
243 :     print STDERR "could not handle hash=$md5\n";
244 :     }
245 :     else
246 :     {
247 :     my($i,$x);
248 :     for ($i=0; ($i < @pegs) && (! ($x = $fig->get_translation($pegs[$i]))); $i++) {}
249 :     if ($i < @pegs)
250 :     {
251 : overbeek 1.3 if (length($x) > 10)
252 :     {
253 :     $count++;
254 :     print FASTA ">$md5\n$x\n";
255 :     }
256 : overbeek 1.1 }
257 :     my $hits = join(";",map { "$_," . $fig->feature_location($_) } @pegs);
258 :     print SUMMARY join("\t",($md5,length($x),$hits)),"\n";
259 :     }
260 :     }
261 :     close(FASTA);
262 :     close(SUMMARY);
263 : overbeek 1.3 if ($count)
264 :     {
265 :     system "formatdb -i $fasta -p T";
266 :     }
267 : overbeek 1.1 $$nP++;
268 :     }
269 :     # if ($bug) { die "aborted" }
270 :     }
271 :    
272 :     sub reduce {
273 :     my($fig,$set,$prot_hash) = @_;
274 :    
275 :     my @setE = ();
276 :     my $x;
277 :     foreach $x (@$set)
278 :     {
279 :     my $toL;
280 :     if ($toL = $prot_hash->{$x})
281 :     {
282 :     push(@setE,split(/,/,$toL));
283 :     }
284 :     elsif ($x =~ /^fig\|/)
285 :     {
286 :     push(@setE,$x);
287 :     }
288 :     }
289 :     # print STDERR &Dumper(['expanded',\@setE]);
290 :     my %reduced = map { $_ => 1 } @setE;
291 :     my @reduced_set = sort { &FIG::by_fig_id($a,$b) } keys(%reduced);
292 :     # print STDERR &Dumper(\@reduced_set);
293 :    
294 :     my %md5s;
295 :     foreach my $peg (@reduced_set)
296 :     {
297 :     # print STDERR &Dumper($peg);
298 :     my $md5 = $fig->md5_of_peg($peg);
299 :     if ($md5)
300 :     {
301 :     push(@{$md5s{$md5}},$peg);
302 :     }
303 :     }
304 :     my @md5_set = sort keys(%md5s);
305 :     return @md5_set;
306 :     }
307 :    
308 :     sub build_expand_hash {
309 :     my($fig,$prot_hash) = @_;
310 :    
311 :     # my $prot_hash_tie = tie %$prot_hash, 'DB_File','make_sim_part_prot_hash.db',O_RDWR,0666,$DB_BTREE; return;
312 :     my $prot_hash_tie = tie %$prot_hash, 'DB_File','make_sim_part_prot_hash.db',O_CREAT,0666,$DB_BTREE;
313 :     $prot_hash_tie || die "tie failed";
314 :    
315 :     open(SYMS,"<$FIG_Config::global/peg.synonyms")
316 :     || die "could not open peg.synonyms";
317 :     while (defined($_ = <SYMS>))
318 :     {
319 :     chop;
320 :     my($head,$rest) = split(/\t/,$_);
321 :     if ($rest =~ /fig\|/)
322 :     {
323 :     my @fig = map { ($_ =~ /(fig\|\d+\.\d+\.peg\.\d+)/) ? $1 : () } split(/;/,$rest);
324 :     @fig = grep { $fig->is_complete(&FIG::genome_of($_)) } @fig;
325 :     if (@fig > 0)
326 :     {
327 :     $head =~ s/,\d+$//;
328 :     if ($head =~ /^fig\|/) { push(@fig,$head) }
329 :     $prot_hash->{$head} = join(",",@fig);
330 :     }
331 :     }
332 :     }
333 :     close(SYMS);
334 :     # print STDERR "prot_hash has been constructed\n";
335 :     }
336 :    
337 :     sub make_fc_data {
338 :     my($fig,$dir,$min_occur,$max_dist) = @_;
339 :    
340 :     my %to_part_hash;
341 :     &load_to_part($dir,\%to_part_hash);
342 :    
343 :     my %pairs_hash;
344 :     my $pairs_hash_tie = tie %pairs_hash, 'DB_File',"pairs_hash.$$.db",O_CREAT,0666,$DB_BTREE;
345 :     ($pairs_hash_tie) || die "tieing hash failed";
346 :    
347 :    
348 :     &count_co_occurrences($fig,\%to_part_hash,\%pairs_hash,$max_dist);
349 :     open(FC,"| sort -k 1,2 -n >$dir/PartitionFC") || die "could not open $dir/PartitionFC";
350 :     while (my($pair,$n) = each(%pairs_hash))
351 :     {
352 :     if ($n >= $min_occur)
353 :     {
354 :     my($k1,$k2) = split(/\t/,$pair);
355 :     print FC join("\t",($k1,$k2,$n)),"\n";
356 :     print FC join("\t",($k2,$k1,$n)),"\n";
357 :     }
358 :     }
359 :     close(FC);
360 :     untie %pairs_hash;
361 :     untie %to_part_hash;
362 :     unlink("to_part_hash.$$.db","pairs_hash.$$.db");
363 :     }
364 :    
365 :     sub load_to_part {
366 :     my($dir,$to_part_hash) = @_;
367 :    
368 :     $to_part_hash_tie = tie %$to_part_hash, 'DB_File',"to_part_hash.28042.db",O_RDWR,0666,$DB_BTREE;
369 :     return;
370 :    
371 :     my $to_part_hash_tie = tie %$to_part_hash, 'DB_File',"to_part_hash.$$.db",O_CREAT,0666,$DB_BTREE;
372 :     ($to_part_hash_tie) || die "tieing hash failed";
373 :    
374 :     open(PARTS,"<$dir/ToSimilarityPartition")
375 :     || die "could not open $dir/ToSimilarityPartition";
376 :     while (defined($_ = <PARTS>))
377 :     {
378 :     if ($_ =~ /^(\d+)\t(\S+)/)
379 :     {
380 :     $to_part_hash->{$2} = $1;
381 :     }
382 :     }
383 :     close(PARTS);
384 :     }
385 :    
386 :    
387 :    
388 :     sub count_co_occurrences {
389 :     my($fig,$to_part_hash,$pairs_hash,$max_dist) = @_;
390 :    
391 :     my $genome;
392 :     foreach $genome (grep { $fig->is_prokaryotic($_) } $fig->genomes('complete'))
393 :     {
394 :     my @loc_md5_pairs;
395 :     foreach my $peg ($fig->all_features($genome,"peg"))
396 :     {
397 :     my $loc = $fig->feature_location($peg);
398 :     my($contig,$beg,$end) = &FIG::boundaries_of($loc);
399 :     my $md5 = $fig->md5_of_peg($peg);
400 :     my $p1;
401 :     if ($md5 && $contig && ($p1 = $to_part_hash->{$md5}))
402 :     {
403 :     push(@loc_md5_pairs,[$contig,&FIG::min($beg,$end),&FIG::max($beg,$end),$p1,$peg]);
404 :     }
405 :     }
406 :     @loc_md5_pairs = sort { ($a->[0] cmp $b->[0]) or ($a->[1] <=> $b->[1]) or ($a->[2] <=> $b->[2]) }
407 :     @loc_md5_pairs;
408 :    
409 :     my($i,$j);
410 :     for ($i=0; ($i < (@loc_md5_pairs - 1)); $i++)
411 :     {
412 :     for ($j = $i+1; ($j < @loc_md5_pairs) &&
413 :     ($loc_md5_pairs[$j]->[0] eq $loc_md5_pairs[$i]->[0]) &&
414 :     (($loc_md5_pairs[$j]->[1] - $loc_md5_pairs[$i]->[2]) <= $max_dist);
415 :     $j++)
416 :     {
417 :     if ($loc_md5_pairs[$i]->[3] ne $loc_md5_pairs[$j]->[3])
418 :     {
419 :     my $pair = ($loc_md5_pairs[$i]->[3] lt $loc_md5_pairs[$j]->[3]) ?
420 :     "$loc_md5_pairs[$i]->[3]\t$loc_md5_pairs[$j]->[3]" :
421 :     "$loc_md5_pairs[$j]->[3]\t$loc_md5_pairs[$i]->[3]";
422 :     $pairs_hash->{$pair}++;
423 :     }
424 :     }
425 :     }
426 :     }
427 :     }
428 :    
429 :     1;

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3