[Bio] / Sprout / km_build_Data.pl Repository:
ViewVC logotype

Annotation of /Sprout/km_build_Data.pl

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1.2 - (view) (download) (as text)

1 : parrello 1.1 ########################################################################
2 :     ########################################################################
3 :     use strict;
4 :     use Data::Dumper;
5 :     use SeedUtils;
6 :     use Getopt::Long;
7 :    
8 :     my $usage = "usage: km_build_Data -d DataDir -k 8 -y reduction_count\n";
9 :     my $dataD;
10 :     my $k = 8;
11 :     my $aaLen = 0;
12 :     my $rast_dirs;
13 :     my $aaList = "LAGVEISRDTKPFNQYMHWC";
14 :    
15 :     my $use_pub_seed = 0; #### default is set to core_seed ###
16 :     my $rc = GetOptions('d=s' => \$dataD,
17 :     'p' => \$use_pub_seed,
18 :     'r=s' => \$rast_dirs,
19 :     'k=i' => \$k,
20 :     'y=i' => \$aaLen);
21 :    
22 :     # properties vs normal-search is determined by the presence of
23 :     # $dataD/genomes and $dataD/properties;
24 :    
25 :     if ((! $rc) || (! -d $dataD) || (! $k))
26 :     {
27 :     print STDERR $usage; exit ;
28 :     }
29 :     (((-s "$dataD/genomes") || (-s "$dataD/properties")) &&
30 :     (! ((-s "$dataD/genomes") && (-s "$dataD/properties"))))
31 :     || die "you need $dataD to contain exactly one of { genomes,properties }";
32 :    
33 :     unlink("$dataD/final.kmers");
34 :     unlink("$dataD/kmer.table.mem_map"); ### force rebuild of memory map
35 :     # print STDERR "Building kmers of size $k\n";
36 :    
37 :     my $distinct_signatures = 0; # number of distinct signature Kmers
38 :     my $distinct_functions = {}; # functions having a signature
39 :     my $seqs_with_func = {}; # hash giving number of seqs with a given function
40 :     my $seqs_with_a_signature = {}; # hash with keys of seqIDs - used to get number seqs with signatures
41 :     print "Building OTUs.\n";
42 :     my $g_to_o = &build_otu_index($dataD);
43 :     print "Building function index.\n";
44 :     &build_function_index($dataD,$use_pub_seed,$rast_dirs); # builds the function index
45 :     print "Building reduced kmers.\n";
46 :     &build_reduced_kmers($dataD,
47 :     $use_pub_seed,
48 :     $rast_dirs,
49 :     $g_to_o,
50 :     \$distinct_signatures,
51 :     $distinct_functions,
52 :     $seqs_with_func,
53 :     $seqs_with_a_signature);
54 :     print "Making final files.\n";
55 :     my $sz = &make_final($dataD,
56 :     $seqs_with_a_signature,
57 :     $distinct_signatures,
58 :     $distinct_functions,
59 :     $seqs_with_func);
60 :     open(SZ,">$dataD/size") || die "could not write size = $sz";
61 :     print SZ "$sz\n";
62 :     close(SZ);
63 :     print "ALL DONE.\n";
64 :    
65 :     sub build_otu_index {
66 :     my($dataD) = @_;
67 :     my $g_to_o = {};
68 :     my %counts;
69 :     if (-s "$dataD/genomes")
70 :     {
71 :     foreach $_ (`cat $dataD/genomes`)
72 :     {
73 :     chomp;
74 :     my($name,$genome) = split(/\t/,$_);
75 :     if ($name =~ /^(\S+)\s(\S+)/)
76 :     {
77 :     my $genus = $1;
78 :     my $species = $2;
79 :     if ($species !~ /^sp\.?$/i)
80 :     {
81 :     if ($genus !~ /^(other|unclassified)/)
82 :     {
83 :     $counts{"$genus $species"}++;
84 :     $g_to_o->{$genome} = "$genus $species";
85 :     }
86 :     }
87 :     }
88 :     }
89 :     }
90 :     else # if a property search
91 :     {
92 :     (-s "$dataD/properties") || die "you need a file $dataD/genomes or a file $dataD/properties";
93 :     foreach $_ (`cat $dataD/properties`)
94 :     {
95 :     chomp;
96 :     my($property,$genome) = split(/\t/,$_);
97 :     if (defined($genome) && defined($property))
98 :     {
99 :     $counts{$property}++;
100 :     $g_to_o->{$genome} = $property;
101 :     }
102 :     }
103 :     }
104 :     open(WTS,">$dataD/otu.occurrences") || die "could not open $dataD/otu.occurrences";
105 :     my $nxt = 0;
106 :     foreach my $gs_or_prop (sort keys(%counts))
107 :     {
108 :     print WTS join("\t",($nxt,$counts{$gs_or_prop},$gs_or_prop)),"\n";
109 :     $nxt++;
110 :     }
111 :     close(WTS);
112 :    
113 :     system "cut -f1,3 $dataD/otu.occurrences > $dataD/otu.index";
114 :     return $g_to_o;
115 :     }
116 :    
117 :     sub build_function_index {
118 :     my($dataD,$use_pub_seed,$rast_dirs) = @_;
119 :    
120 :     open(FI,">$dataD/function.index") || die "could not open $dataD/function.index";
121 :     if (-s "$dataD/genomes")
122 :     {
123 :     my $cmd;
124 :     if ($use_pub_seed)
125 :     {
126 :     $cmd = "cut -f2 $dataD/genomes | svr_all_features peg | svr_function_of";
127 :     }
128 :     elsif ($rast_dirs)
129 :     {
130 :     $cmd = "cat $rast_dirs/*/*functions";
131 :     }
132 :     else
133 :     {
134 :     $cmd = "cut -f2 $dataD/genomes | genomes_to_fids | grep peg | fids_to_functions | cut -f2,3";
135 :     }
136 :     open(ASSIGNMENTS,"$cmd |")
137 :     || die "cannot access assignments";
138 :    
139 :     my %funcs;
140 :     while (defined($_ = <ASSIGNMENTS>))
141 :     {
142 :     if ($_ =~ /^\S+\t(\S.*\S+)/)
143 :     {
144 :     my $stripped = &SeedUtils::strip_func_comment($1);
145 :     $funcs{$stripped}++;
146 :     }
147 :     }
148 :     close(ASSIGNMENTS);
149 :    
150 :     my $nxt = 0;
151 :     foreach my $f (sort keys(%funcs))
152 :     {
153 :     if ($funcs{$f} > 5) # && ((! &SeedUtils::hypo($f)) || ($f =~ /FIG/)))
154 :     {
155 :     print FI "$nxt\t$f\n";
156 :     $nxt++;
157 :     }
158 :     }
159 :     }
160 :     else # we are in a property Data directory
161 :     {
162 :     my $sets = {};
163 :     (-s "$dataD/properties") || die "something is wrong";
164 :     foreach $_ (`cat $dataD/properties`)
165 :     {
166 :     if ($_ =~ /^(\S+)\t(\S+)$/)
167 :     {
168 :     my($p,$g) = ($1,$2);
169 :     $sets->{$p}->{$g} = 1;
170 :     }
171 :     else
172 :     {
173 :     die "Bad entry in properties file: $_";
174 :     }
175 :     }
176 :     my $nxt = 0;
177 :     foreach my $p (sort keys(%$sets))
178 :     {
179 :     print FI join("\t",($nxt,$p)),"\n";
180 :     $nxt++;
181 :     }
182 :     }
183 :     close(FI);
184 :     }
185 :    
186 :     sub build_reduced_kmers {
187 :     my($dataD,
188 :     $use_pub_seed,
189 :     $rast_dirs,
190 :     $g_to_o,
191 :     $distinct_signatures,
192 :     $distinct_functions,
193 :     $seqs_with_func,
194 :     $seqs_with_a_signature) = @_;
195 :    
196 :     my %to_oI;
197 :     foreach $_ (`cat $dataD/otu.occurrences`)
198 :     {
199 :     if ($_ =~ /^(\S+)\t(\S+)\t(\S.*\S)/)
200 :     {
201 :     $to_oI{$3} = $1;
202 :     }
203 :     }
204 :    
205 :     my %to_fI;
206 :     foreach $_ (`cat $dataD/function.index`)
207 :     {
208 :     if ($_ =~ /^(\S+)\t(\S.*\S)/)
209 :     {
210 :     $to_fI{$2} = $1;
211 :     }
212 :     }
213 :    
214 :     my %id_to_fI;
215 :     &load_id_to_fI($dataD,$use_pub_seed,$rast_dirs,\%to_fI,\%id_to_fI);
216 :     my $properties = -s "$dataD/properties";
217 :     ###
218 :     # First, we build a file containing [kmer,fI,oI,off-set,sequence-ID]
219 :     # sorted by kmer
220 :     ###
221 :     open(RAW,"| sort -T . -S 1G > $dataD/sorted.kmers") || die "could not open $dataD/sorted.kmers";
222 :     my $seqID=0;
223 :     my @genomes = map { chomp; $_ } ($properties ? `cut -f2 $dataD/properties` : `cut -f2 $dataD/genomes`);
224 :     foreach my $g (sort @genomes)
225 :     {
226 :     my $cmd;
227 :     if ($use_pub_seed)
228 :     {
229 :     $cmd = "echo '$g' | svr_all_features peg | svr_translations_of";
230 :     }
231 :     elsif ($rast_dirs)
232 :     {
233 :     $cmd = "flatten_fasta < $rast_dirs/$g/Features/peg/fasta";
234 :     }
235 :     else
236 :     {
237 :     $cmd = "echo '$g' | genomes_to_fids | grep peg | fids_to_protein_sequences -fasta 0 | cut -f2,3";
238 :     }
239 :     # print "TRACE: command is $cmd\n";
240 :     foreach $_ (`$cmd`)
241 :     {
242 :     # print "TRACE: processing input line: $_\n";
243 :     if ($_ =~ /^(fig\|\d+\.\d+\.peg\.\d+)\t(\S.*\S)$/)
244 :     {
245 :     chomp;
246 :     my $seq = $2;
247 :     my $id = $1;
248 :     my $sub;
249 :     # The eval trick enables the command-line to specify how drastic the letter elimination should be.
250 :     # A value of 2 removes only L and A. A value of 4 removes L, A, G, and V. A value of 20 removes
251 :     # everything. (This would be bad.)
252 :     if ($aaLen > 0) {
253 :     $sub = substr($aaList, 0, $aaLen);
254 :     my $seq1 = $seq;
255 :     eval("\$seq =~ tr/$sub//d");
256 :     # print "($sub) $seq1 => $seq\n";
257 :     }
258 :     ($id =~ /^fig\|(\d+\.\d+)/) || die "bad peg $_";
259 :     my $oI = $to_oI{$g_to_o->{$g}};
260 :     my $fI = $id_to_fI{$id};
261 :     if (! defined($fI)) { $fI = $id_to_fI{"hypothetical protein"} } # default is hypothetical
262 :     if (defined($fI))
263 :     {
264 :     $seqs_with_func->{$fI}++; # NFj
265 :     for (my $i=0; ($i < (length($seq) - $k)); $i++)
266 :     {
267 :     my $kmer = uc substr($seq,$i,$k);
268 :     if ($kmer !~ /[^ACDEFGHIKLMNPQRSTVWY]/)
269 :     {
270 :     print RAW join("\t",($kmer,
271 :     $fI,
272 :     $oI,
273 :     (length($seq)-$i),
274 :     $seqID)),"\n";
275 :     }
276 :     }
277 :     $seqID++;
278 :     }
279 :     }
280 :     }
281 :     }
282 :     close(RAW);
283 :    
284 :     ### Now, we reduce sets of adjacent lines with the same kmer to a single line (or none)
285 :     ###
286 :     open(RAW,"<$dataD/sorted.kmers") || die "could not open sorted.kmers";
287 :     open(REDUCED,">$dataD/reduced.kmers") || die "could not open reduced kmers";
288 :     my $last = <RAW>;
289 :    
290 :     while ($last && ($last =~ /^(\S+)/))
291 :     {
292 :     my $curr = $1;
293 :     my @set;
294 :     while ($last && ($last =~ /^(\S+)\t(\S*)\t(\S*)\t(\S*)\t(\S+)$/) && ($1 eq $curr))
295 :     {
296 :     push(@set,[$2,$3,$4,$5]);
297 :     $last = <RAW>;
298 :     }
299 :     &process_set($curr,\@set,\*REDUCED,$seqs_with_a_signature,$distinct_signatures,$distinct_functions);
300 :     }
301 :     close(REDUCED);
302 :     unlink("$dataD/sorted.kmers");
303 :     }
304 :     sub process_set {
305 :     my($kmer,$set,$fh,$seqs_with_a_signature,$distinct_signatures,$distinct_functions) = @_;
306 :    
307 :     my %funcs;
308 :     my $tot = 0;
309 :     foreach my $tuple (@$set)
310 :     {
311 :     my($fI,$oI,$off,$seqID) = @$tuple;
312 :     $tot++;
313 :     $funcs{$fI}++;
314 :     }
315 :    
316 :     my $seqID;
317 :     my @tmp = sort { $funcs{$b} <=> $funcs{$a} } keys(%funcs);
318 :     if (defined($tmp[0]) && ($funcs{$tmp[0]} >= (0.8 * $tot)))
319 :     {
320 :     my $best_fI = $tmp[0];
321 :     # my $func_wt = $funcs{$tmp[0]};
322 :     my %otus;
323 :     my $otu = '';
324 :     my $seqs_containing_func = 0;
325 :     my $seqs_containing_sig = @$set;
326 :     foreach my $tuple (@$set)
327 :     {
328 :     my($fI,$oI,$off,$seqID) = @$tuple;
329 :     $seqs_with_a_signature->{$seqID} = 1;
330 :     if (($fI == $best_fI) && defined($oI))
331 :     {
332 :     $seqs_containing_func++;
333 :     $otus{$oI}++;
334 :     }
335 :     }
336 :     @tmp = sort { $otus{$b} <=> $otus{$a} } keys(%otus);
337 :     if (defined($tmp[0]) && ($otus{$tmp[0]} >= (0.8 * $tot)))
338 :     {
339 :     $otu = $tmp[0];
340 :     }
341 :     my @offsets = sort { $a <=> $b } map { $_->[2] } @$set;
342 :     my $median_off = $offsets[int(scalar @offsets / 2)];
343 :     print $fh join("\t",($kmer,
344 :     $median_off,
345 :     $best_fI,
346 :     $otu,
347 :     $seqs_containing_sig,
348 :     $seqs_containing_func)),"\n";
349 :     $$distinct_signatures++;
350 :     $distinct_functions->{$best_fI} = 1;
351 :     }
352 :     }
353 :    
354 :     sub load_id_to_fI {
355 :     my($dataD,$use_pub_seed,$rast_dirs,$to_fI,$id_to_fI) = @_;
356 :    
357 :     my $properties = -s "$dataD/properties";
358 :     my %genomes = map { ($_ =~ /^(\S[^\t]*\S)\t(\S+)/) ? ($2 => $1) : () }
359 :     ($properties ? `cat $dataD/properties` : `cat $dataD/genomes`);
360 :     foreach my $genome (keys(%genomes))
361 :     {
362 :     my $cmd;
363 :     if ($use_pub_seed)
364 :     {
365 :     $cmd = "echo '$genome' | svr_all_features peg | svr_function_of";
366 :     }
367 :     elsif ($rast_dirs)
368 :     {
369 :     $cmd = "cat $rast_dirs/$genome/*functions";
370 :     }
371 :     else
372 :     {
373 :     $cmd = "echo '$genome' | genomes_to_fids | grep peg | fids_to_functions | cut -f2,3";
374 :     }
375 :     foreach $_ (`$cmd`)
376 :     {
377 :     if ($_ =~ /^(\S+)\t(\S.*\S)/)
378 :     {
379 :     my $peg = $1;
380 :     my $func_or_prop = $properties ? $genomes{$genome} : $2;
381 :     $func_or_prop = &SeedUtils::strip_func_comment($func_or_prop);
382 :     if (defined($to_fI->{$func_or_prop}))
383 :     {
384 :     $id_to_fI->{$peg} = $to_fI->{$func_or_prop};
385 :     }
386 :     }
387 :     }
388 :     }
389 :     }
390 :    
391 :     sub make_final {
392 :     my($dataD,
393 :     $seqs_with_a_signature,
394 :     $distinct_signatures,
395 :     $distinct_functions,
396 :     $seqs_with_func) = @_;
397 :    
398 :     my $size = 0;
399 :     my %to_oI;
400 :     foreach $_ (`cat $dataD/otu.occurrences`)
401 :     {
402 :     if ($_ =~ /^(\d+)\t\d+\t(\S+)/)
403 :     {
404 :     $to_oI{$2} = $1;
405 :     }
406 :     }
407 :    
408 :     my $num_seqs_with_a_signature = keys(%$seqs_with_a_signature);
409 :     my $num_distinct_functions = keys(%$distinct_functions);
410 :     open(IN,"<$dataD/reduced.kmers") || die "could not open final.kmers";
411 :     open(OUT,">$dataD/final.kmers") || die "could not open $dataD/final.signatures";
412 :     while (defined($_ = <IN>))
413 :     {
414 :     chop;
415 :     my($kmer,$off,$fI,$otu,$seqs_containing_sig,$seqs_containing_func) = split(/\t/,$_);
416 :     my $fI_wt = &compute_weight_of_signature($num_seqs_with_a_signature,
417 :     $distinct_signatures,
418 :     $num_distinct_functions,
419 :     $seqs_containing_sig,
420 :     $seqs_with_func->{$fI},
421 :     $seqs_containing_func);
422 :     print OUT join("\t",($kmer,
423 :     $off,
424 :     $fI,
425 :     $fI_wt,
426 :     $otu)),"\n";
427 :     $size++;
428 :     }
429 :     close(IN);
430 :     unlink("$dataD/reduced.kmers");
431 :     close(OUT);
432 :     return $size;
433 :     }
434 :    
435 :    
436 :     # The meaning of the variables is now as follows:
437 :    
438 :     # $NSF := number of training sequences with at least one signature and an assigned function.
439 :     # $KS := number of signatures of function.
440 :     # $KF := number of distinct functions with at least one signature.
441 :     # $NSi := number of training sequences with at least one occurrence of signature Si.
442 :     # $NFj := number of training sequences having function Fj.
443 :     # $NSiFj := number of sequences with at least one occurrence of Si and having function Fj
444 :    
445 :     # Logical consistency requires that the following identities be satisfied:
446 :    
447 :     # $NSF = Sum(i=1..$KS, Sum(j=1..$KF, $NSiFj));
448 :     # $NSi = Sum(j=1..$KF, $NSiFj);
449 :     # $NFj = Sum(i=1..$KS, $NSiFj);
450 :    
451 :     sub compute_weight_of_signature {
452 :     my ($NSF, $KS, $KF, $NSi, $NFj, $NSiFj) = @_;
453 :    
454 :     return sprintf("%0.4f",(log( ($NSiFj + 1.0) / ($NSi - $NSiFj + 1.0) ) +
455 :     log( ($NSF - $NFj + $KS) / ($NFj + $KS) )));
456 :     }
457 :    
458 :     sub compute_weight_of_prior {
459 :     my ($NSF, $KF, $NFj) = @_;
460 :    
461 :     return (log( ($NFj + 1.0) / ($NSF + $KF - $NFj - 1.0) ) );
462 :     }

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3