[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.2 - (view) (download) (as text)

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

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3