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

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

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3