[Bio] / Kmers2 / build_data_directory.pl Repository:
ViewVC logotype

Annotation of /Kmers2/build_data_directory.pl

Parent Directory Parent Directory | Revision Log Revision Log


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

1 : overbeek 1.1 use strict;
2 :     use Data::Dumper;
3 :     use SeedEnv;
4 : overbeek 1.4 use FIG;
5 :     my $fig = new FIG;
6 : overbeek 1.1
7 : overbeek 1.2 my $usage = "usage: build_data_directory DataDir K";
8 : overbeek 1.1 my($genus,$dataD,$k);
9 :     (
10 :     ($dataD = shift @ARGV) &&
11 :     ($k = shift @ARGV) && ($k =~ /^\d{1,2}$/)
12 :     )
13 :     || die $usage;
14 :     print STDERR "building $dataD for Kmers of size $k\n";
15 :     mkdir($dataD,0777);
16 :    
17 : overbeek 1.4 my %seen;
18 :     my $existing = $ENV{'SAS_SERVER'};
19 :     $ENV{'SAS_SERVER'} = 'PUBSEED';
20 :     open(G,">$dataD/all.genomes") || die "could not open $dataD/all.genomes";
21 :     foreach $_ (`svr_all_genomes`)
22 :     {
23 :     if (($_ !~ /[pP]hage\b/) && ($_ !~ /[pP]lasmid\b/))
24 :     {
25 :     if (($_ =~ /^(\S+)\s(\S+).*\s(\d+\.\d+)$/) && ($2 ne "sp."))
26 :     {
27 :     my $genus = $1;
28 :     my $species = $2;
29 :     my $genome = $3;
30 :     if ($species !~ /^sp\.?$/i)
31 :     {
32 :     my $md5 = $fig->genome_md5sum($genome);
33 :     if ((! $seen{$md5}) && $fig->is_prokaryotic($genome))
34 :     {
35 :     $seen{$md5} = 1;
36 :     print G $_;
37 :     }
38 :     }
39 :     }
40 :     }
41 :     }
42 :     close(G);
43 :     $ENV{'SAS_SERVER'} = $existing;
44 :    
45 : overbeek 1.2 &build_function_index($dataD);
46 :     &build_otu_index($dataD);
47 : overbeek 1.1 &build_reduced_kmers($dataD,$k);
48 :     &extend_otu_index($dataD);
49 :     &update_kmers_with_extended_OTUs($dataD);
50 :    
51 :     sub build_function_index {
52 :     my($dataD) = @_;
53 :    
54 :     open(ASSIGNMENTS,"cut -f2 $dataD/all.genomes | svr_all_features peg | svr_function_of |")
55 :     || die "cannot access assignments";
56 :    
57 :     my %funcs;
58 :     while (defined($_ = <ASSIGNMENTS>))
59 :     {
60 :     if ($_ =~ /^\S+\t(\S[^\#]*\S+)/)
61 :     {
62 :     $funcs{$1}++;
63 :     }
64 :     }
65 :     close(ASSIGNMENTS);
66 :     open(FI,">$dataD/function.index") || die "could not open $dataD/function.index";
67 :    
68 :     my $nxt = 0;
69 :     foreach my $f (sort keys(%funcs))
70 :     {
71 :     if (($funcs{$f} > 5) && ((! &SeedUtils::hypo($f)) || ($f =~ /FIG/)))
72 :     {
73 :     print FI "$nxt\t$f\n";
74 :     $nxt++;
75 :     }
76 :     }
77 :     close(FI);
78 :     }
79 :    
80 :     sub build_otu_index {
81 :     my($dataD) = @_;
82 :    
83 :     my %otus;
84 :     foreach $_ (`cat $dataD/all.genomes`)
85 :     {
86 :     if ($_ =~ /(\S+\s\S+)/)
87 :     {
88 :     $otus{$1}++;
89 :     }
90 :     }
91 :    
92 :     open(OI,">$dataD/otu.index") || die "could not open $dataD/otu.index";
93 :     my $nxt = 0;
94 :     foreach my $otu (keys(%otus))
95 :     {
96 :     if ($otu !~ / sp\.$/)
97 :     {
98 :     print OI "$nxt\t$otu\n";
99 :     $nxt++;
100 :     }
101 :     }
102 :     close(OI);
103 :     }
104 :    
105 :     sub build_reduced_kmers {
106 :     my($dataD,$k) = @_;
107 :    
108 :     my %to_oI;
109 :     foreach $_ (`cat $dataD/otu.index`)
110 :     {
111 :     if ($_ =~ /^(\S+)\t(\S.*\S)/)
112 :     {
113 :     $to_oI{$2} = $1;
114 :     }
115 :     }
116 :    
117 :     my %g_to_oI;
118 :    
119 :     foreach $_ (`cat $dataD/all.genomes`)
120 :     {
121 :     if ($_ =~ /^(\S[^\t]+)\t(\S.*\S)/)
122 :     {
123 :     my $g = $2;
124 :     my $gs = $1;
125 :     if (($gs =~ /^(\S+ \S+)/) && defined($to_oI{$1}))
126 :     {
127 :     $g_to_oI{$g} = $to_oI{$1};
128 :     }
129 :     }
130 :     }
131 :    
132 :     my %to_fI;
133 :     foreach $_ (`cat $dataD/function.index`)
134 :     {
135 :     if ($_ =~ /^(\S+)\t(\S.*\S)/)
136 :     {
137 :     $to_fI{$2} = $1;
138 :     }
139 :     }
140 : overbeek 1.4 #
141 :     # we take assignments from both the pubSEED and the ASEED. ASEED functions
142 :     # override PSEED functions
143 :     #
144 : overbeek 1.1 my %peg_to_fI;
145 : overbeek 1.4 &load_peg_to_fI($dataD,\%to_fI,'PUBSEED',\%peg_to_fI);
146 :     &load_peg_to_fI($dataD,\%to_fI,'SEED',\%peg_to_fI);
147 : overbeek 1.1
148 :     open(RAW,"| sort -T . > $dataD/sorted.kmers") || die "could not open $dataD/sorted.kmers";
149 :     foreach my $g (`cut -f2 $dataD/all.genomes`)
150 :     {
151 : overbeek 1.2 # next if (! $g_to_oI{$g});
152 : overbeek 1.1 foreach $_ (`echo '$g' | svr_all_features peg | svr_translations_of`)
153 :     {
154 :     if ($_ =~ /^(fig\|\d+\.\d+\.peg\.\d+)\t(\S.*\S)$/)
155 :     {
156 :     chomp;
157 :     my $seq = $2;
158 :     my $id = $1;
159 :     ($id =~ /^fig\|(\d+\.\d+)/) || die "bad peg $_";
160 :     my $g = $1;
161 :     my $oI = $g_to_oI{$g};
162 :     my $fI = $peg_to_fI{$id};
163 :     for (my $i=0; ($i < (length($seq) - $k)); $i++)
164 :     {
165 :     my $kmer = uc substr($seq,$i,$k);
166 :     if ($kmer !~ /[^ACDEFGHIKLMNPQRSTVWY]/)
167 :     {
168 :     print RAW join("\t",($kmer,$fI,$oI,length($seq)-$i)),"\n";
169 :     }
170 :     }
171 :     }
172 :     }
173 :     }
174 :     close(RAW);
175 :    
176 :     open(RAW,"<$dataD/sorted.kmers") || die "could not open sorted.kmers";
177 :     open(REDUCED,">$dataD/reduced.kmers") || die "could not open reduced kmers";
178 :     my $last = <RAW>;
179 :     while ($last && ($last =~ /^(\S+)/))
180 :     {
181 :     my $curr = $1;
182 :     my @set;
183 :     while ($last && ($last =~ /^(\S+)\t(\S*)\t(\S*)\t(\S*)$/) && ($1 eq $curr))
184 :     {
185 :     push(@set,[$2,$3,$4]);
186 :     $last = <RAW>;
187 :     }
188 :     &process_set($curr,\@set,\*REDUCED);
189 :     }
190 :     close(REDUCED);
191 :     }
192 :    
193 : overbeek 1.4 sub load_peg_to_fI {
194 :     my($dataD,$to_fI,$which_seed,$peg_to_fI) = @_;
195 :     my $existing = $ENV{'SAS_SERVER'};
196 :     $ENV{'SAS_SERVER'} = $which_seed;
197 :    
198 :     foreach $_ (`cut -f2 $dataD/all.genomes | svr_all_features peg | svr_function_of`)
199 :     {
200 :     if ($_ =~ /^(\S+)\t(\S.*\S)/)
201 :     {
202 :     my $peg = $1;
203 :     my $func = $2;
204 :     if ($to_fI->{$func})
205 :     {
206 :     $peg_to_fI->{$peg} = $to_fI->{$func};
207 :     }
208 :     }
209 :     }
210 :     $ENV{'SAS_SERVER'} = $existing;
211 :     }
212 :    
213 :    
214 : overbeek 1.1 sub process_set {
215 :     my($kmer,$set,$fh) = @_;
216 :    
217 :     my %funcs;
218 :     foreach my $tuple (@$set)
219 :     {
220 :     my($fI,$oI,$off) = @$tuple;
221 :     $funcs{$fI}++;
222 :     }
223 :     my @tmp = sort { $funcs{$b} <=> $funcs{$a} } keys(%funcs);
224 : overbeek 1.3 if ($tmp[0] && ($funcs{$tmp[0]} > (0.5 * @$set)))
225 : overbeek 1.1 {
226 :     my $best_fI = $tmp[0];
227 :     my %otus;
228 :     my $otu = '';
229 :     foreach my $tuple (@$set)
230 :     {
231 :     my($fI,$oI,$off) = @$tuple;
232 :     if (($fI == $best_fI) && $oI)
233 :     {
234 :     $otus{$oI}++;
235 :     }
236 :     }
237 :     @tmp = sort { $otus{$b} <=> $otus{$a} } keys(%otus);
238 :     if ($tmp[0] && ($_ = &pickable_otu(\@tmp,\%otus,scalar @$set)))
239 :     {
240 :     $otu = $_;
241 :     }
242 :     my $sum_off = 0;
243 :     my $n = 0;
244 :     foreach my $tuple (@$set)
245 :     {
246 :     my($fI,$oI,$off) = @$tuple;
247 :     if ($fI == $best_fI)
248 :     {
249 :     $sum_off += $off;
250 :     $n++;
251 :     }
252 :     }
253 :     print $fh join("\t",($kmer,int($sum_off/$n),$best_fI,$otu)),"\n";
254 :     }
255 :     }
256 :    
257 :     sub pickable_otu {
258 :     my($poss_otus,$counts,$in_set) = @_;
259 :    
260 :     my @called;
261 :     if (@$poss_otus == 1) { return $poss_otus->[0] }
262 :     my $best;
263 :     while (($best = shift @$poss_otus) && ($counts->{$best} >= (0.2 * $in_set)))
264 :     {
265 :     push(@called,$best);
266 :     }
267 :     if (@called > 0)
268 :     {
269 : overbeek 1.2 return join(",",sort { $a <=> $b } @called);
270 : overbeek 1.1 }
271 :     return '';
272 :     }
273 :    
274 :     sub extend_otu_index {
275 :     my($DataD) = @_;
276 :    
277 :     my @otus = map { chop; [split(/\t/,$_)] } `cat $dataD/otu.index`;
278 :     my $nxt = $otus[-1]->[0] + 1;
279 :     open(EXT,">>$DataD/otu.index") || die "could not extend otu.index";
280 :     foreach $_ (`cut -f4 $DataD/reduced.kmers | grep ',' | sort -u`)
281 :     {
282 :     print EXT join("\t",($nxt++,$_));
283 :     }
284 :     close(EXT);
285 :     }
286 :    
287 :     sub update_kmers_with_extended_OTUs {
288 :     my($dataD) = @_;
289 :    
290 :     my %to_oI;
291 :     foreach $_ (`cat $dataD/otu.index`)
292 :     {
293 :     if ($_ =~ /^(\d+)\t(\S+)/)
294 :     {
295 :     $to_oI{$2} = $1;
296 :     }
297 :     }
298 :     open(IN,"<$dataD/reduced.kmers") || die "bad";
299 :     open(OUT,">$dataD/final.kmers") || die "could not open final.kmers";
300 :     while (defined($_ = <IN>))
301 :     {
302 :     if ($_ !~ /,/) { print OUT $_ }
303 :     else
304 :     {
305 :     chop;
306 :     my($kmer,$off,$fI,$otus) = split(/\t/,$_);
307 :     my $otu = $to_oI{$otus};
308 :     print OUT join("\t",($kmer,$off,$fI,$otu)),"\n";
309 :     }
310 :     }
311 :     close(IN);
312 :     close(OUT);
313 :     }
314 :    
315 :    
316 :    

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3