[Bio] / FigKernelScripts / FFB2_create_binary_kmers.pl Repository:
ViewVC logotype

Annotation of /FigKernelScripts/FFB2_create_binary_kmers.pl

Parent Directory Parent Directory | Revision Log Revision Log


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

1 : overbeek 1.1 ########################################################################
2 :     #
3 :     # Copyright (c) 2003-2006 University of Chicago and Fellowship
4 :     # for Interpretations of Genomes. All Rights Reserved.
5 :     #
6 :     # This file is part of the SEED Toolkit.
7 :     #
8 :     # The SEED Toolkit is free software. You can redistribute
9 :     # it and/or modify it under the terms of the SEED Toolkit
10 :     # Public License.
11 :     #
12 :     # You should have received a copy of the SEED Toolkit Public License
13 :     # along with this program; if not write to the University of Chicago
14 :     # at info@ci.uchicago.edu or the Fellowship for Interpretation of
15 :     # Genomes at veronika@thefig.info or download a copy from
16 :     # http://www.theseed.org/LICENSE.TXT.
17 :     #
18 :     ########################################################################
19 :    
20 :    
21 :     use KmersC;
22 :     use strict;
23 :     use Data::Dumper;
24 : olson 1.4 use PerlIO::gzip;
25 :     use SortMerge;
26 : olson 1.2 use Getopt::Long;
27 : overbeek 1.1
28 :     #
29 :     # Scan input once to determine the max value of the value fields. Then create
30 :     # the output file and scan again, writing values.
31 :     #
32 :    
33 :     my @max;
34 :     my @lens;
35 : olson 1.2 my $motif_len;
36 :     my $column_sizes;
37 :    
38 : olson 1.3 #
39 :     # Multiline input is of the form
40 :     # motif column value
41 :     # where the column field defines the column into which the data should go.
42 :     # It is 1-based.
43 :     #
44 :     my $multiline;
45 : olson 1.4 my $merge = '';
46 : olson 1.3
47 : olson 1.2 my $rc = GetOptions("size=s" => \$column_sizes,
48 : olson 1.3 "multiline" => \$multiline,
49 : olson 1.4 "merge=s" => \$merge,
50 : olson 1.3 "length=s" => \$motif_len);
51 : overbeek 1.1
52 : olson 1.2 if (!$rc || @ARGV < 2)
53 : overbeek 1.1 {
54 : olson 1.2 die "Usage: $0 [-s column-sizes] [-l motif-len] input-data-file output-binary-file [cols]\n";
55 :     }
56 :    
57 :     if ($column_sizes)
58 :     {
59 :     @lens = split(/,/, $column_sizes);
60 : overbeek 1.1 print "Set column lengths to @lens\n";
61 :     }
62 :    
63 :     my $in_file = shift;
64 :     my $out_file = shift;
65 :     my @cols = @ARGV;
66 :    
67 : olson 1.4 #
68 :     # for merge, we don't look at input fh.
69 :     #
70 : olson 1.2 my $input_fh;
71 : olson 1.4 if ($merge eq '')
72 : overbeek 1.1 {
73 : olson 1.4 if ($in_file eq '-')
74 : overbeek 1.1 {
75 : olson 1.4 if (@lens == 0 || !defined($motif_len))
76 :     {
77 :     die "In order to read from stdin, both the -s and the -l arguments must be specified\n";
78 :     }
79 :     $input_fh = \*STDIN;
80 : overbeek 1.1 }
81 : olson 1.4 else
82 : overbeek 1.1 {
83 : olson 1.4 open($input_fh, "<", $in_file) or die "Cannot open $in_file: $!";
84 : olson 1.2
85 : olson 1.4 while (<$input_fh>)
86 : olson 1.2 {
87 : olson 1.4 chomp;
88 :     my($motif, @vals) = split(/\t/);
89 :     if (!defined($motif_len))
90 :     {
91 :     $motif_len = length($motif);
92 :     }
93 :    
94 :     if (@lens)
95 :     {
96 :     last;
97 :     }
98 :    
99 :     if (@cols)
100 :     {
101 :     @vals = @vals[@cols];
102 :     }
103 :    
104 :     for my $i (0 .. $#vals)
105 :     {
106 :     $max[$i] = $vals[$i] if $vals[$i] > $max[$i];
107 :     }
108 : olson 1.2 }
109 :    
110 : olson 1.4 seek($input_fh, 0, 0);
111 :     }
112 :     }
113 :     else
114 :     {
115 :     if (@lens == 0 || !defined($motif_len))
116 :     {
117 :     die "In order to merge, both the -s and the -l arguments must be specified\n";
118 : overbeek 1.1 }
119 :     }
120 :    
121 :     if (!@lens)
122 :     {
123 :     for my $m (@max)
124 :     {
125 :     my $l;
126 :     if ($m < 128)
127 :     {
128 :     $l = 1;
129 :     }
130 :     elsif ($m < 32768)
131 :     {
132 :     $l = 2;
133 :     }
134 :     elsif ($m < 2147483648)
135 :     {
136 :     $l = 4;
137 :     }
138 :     else
139 :     {
140 :     die "Max value $m is greater than a 4-byte signed int";
141 :     }
142 :     push(@lens, $l);
143 :     }
144 :     }
145 : olson 1.3
146 :     my $num_keys = @lens;
147 :    
148 : overbeek 1.1 print STDERR "Computed max lengths: @max and byte sizes: @lens\n";
149 :    
150 :     my $cr = new KmersFileCreator(0xfeedface, $motif_len, 0, \@lens);
151 :     $cr->open_file($out_file);
152 :     $cr->write_file_header();
153 :    
154 :     my $count = 0;
155 : olson 1.3
156 : olson 1.4 if ($merge)
157 : overbeek 1.1 {
158 : olson 1.4 my @code;
159 :     my @fh;
160 : olson 1.3
161 : olson 1.4 if (!open(MF, "<", $merge))
162 :     {
163 :     die "Cannot open merge file $merge: $!";
164 :     }
165 :    
166 :     my %by_col;
167 :    
168 :     while (<MF>)
169 : olson 1.3 {
170 : olson 1.4 chomp;
171 :     my($file, $col) = split(/\t/);
172 :     my $fh;
173 :     if ($file =~ /\.gz$/)
174 :     {
175 :     open($fh, "<:gzip", $file) or die "Cannot open $file with gzip: $!";
176 :     # open($fh, "-|", "gunzip", "-c", $file) or die "Cannot open $file with gzip: $!";
177 :     }
178 :     else
179 : olson 1.3 {
180 : olson 1.4 open($fh, "<", $file) or die "Cannot open $file: $!";
181 :     }
182 :     push(@{$by_col{$col}}, [$file, $fh]);
183 :     }
184 :     close(MF);
185 :     for my $col (sort keys %by_col)
186 :     {
187 :     my $fhlist = $by_col{$col};
188 :    
189 :     push @code, sub {
190 :     my $ent = $fhlist->[0];
191 :     my($file, $fh) = @$ent;
192 :     my $l = <$fh>;
193 :     while (!defined($l))
194 : olson 1.3 {
195 : olson 1.4 close($fh);
196 :     shift @$fhlist;
197 :     if (@$fhlist == 0)
198 : olson 1.3 {
199 : olson 1.4 return ();
200 : olson 1.3 }
201 : olson 1.4 my $ent = $fhlist->[0];
202 :     ($file, $fh) = @$ent;
203 :     print STDERR "Moving to $file for column $col\n";
204 :     $l = <$fh>;
205 :     }
206 :    
207 :     if ($l =~ /^(\S+)\t(\d+)/)
208 :     {
209 :     my($motif, $val) = ($1, $2);
210 :     return ($motif, [$motif, $col, $val]);
211 : olson 1.3 }
212 : olson 1.4 else
213 :     {
214 :     die "Badly formatted input at line $. file $file: $l\n";
215 :     }
216 :     };
217 :     }
218 :    
219 :     my $writer = MultilineWriter->new($cr, $num_keys);
220 :    
221 :     Sort::Merge::sort_coderefs(\@code,
222 :     sub {
223 :     my($code, $key, $value) = @{$_[0]};
224 : olson 1.5 # print "@$value\n";
225 : olson 1.4 $writer->push_one_value(@$value);
226 :     });
227 :     $writer->finish();
228 :    
229 :     $count = $writer->{count};
230 :     }
231 :     elsif ($multiline)
232 :     {
233 :     my $writer = MultilineWriter->new($cr, $num_keys);
234 :     while (defined(my $l = <$input_fh>))
235 :     {
236 :     if ($l =~ /^(\S+)\t(\d+)\t(\d+)/)
237 :     {
238 :     my($motif, $col, $val) = ($1, $2, $3);
239 :     $writer->push_one_value($motif, $col, $val);
240 : olson 1.3 }
241 :     else
242 :     {
243 :     warn "Badly formatted input at $.: $l\n";
244 :     }
245 :     }
246 : olson 1.4 close($input_fh);
247 :     $writer->finish();
248 :     $count = $writer->count();
249 : olson 1.3 }
250 :     else
251 :     {
252 :     while (<$input_fh>)
253 :     {
254 :     chomp;
255 :     my($motif, @vals) = split(/\t/);
256 :     if (@cols)
257 :     {
258 :     @vals = @vals[@cols];
259 :     }
260 :     $cr->write_entry($motif, \@vals);
261 :     $count++;
262 : overbeek 1.1 }
263 :     }
264 :     $cr->close_file();
265 :    
266 :     print "Loaded $count oligos\n";
267 : olson 1.4
268 :     package MultilineWriter;
269 :     use strict;
270 :     use base 'Class::Accessor';
271 :    
272 :     __PACKAGE__->mk_accessors(qw(cr cur val count));
273 :    
274 :     sub new
275 :     {
276 :     my($class, $cr, $num_keys) = @_;
277 :    
278 :     my $self = {
279 :     cr => $cr,
280 :     cur => '',
281 :     count => 0,
282 :     num_keys => $num_keys,
283 :     blank_val => [map { -1 } 1..$num_keys],
284 :     };
285 :     bless $self, $class;
286 :     $self->set_blank();
287 :     return $self;
288 :     }
289 :    
290 :     sub set_blank
291 :     {
292 :     my($self) = @_;
293 :     @{$self->{val}} = @{$self->{blank_val}};
294 :     }
295 :    
296 :     sub push_one_value
297 :     {
298 :     my($self, $motif, $col, $val) = @_;
299 : olson 1.5 # print "$self->{cur} $motif $col $val\n";
300 : olson 1.4 if ($self->{cur} ne $motif)
301 :     {
302 :     if ($self->{cur} ne '')
303 :     {
304 :     if (@{$self->{val}} != $self->{num_keys})
305 :     {
306 :     die "Val became invalid: " . Dumper($self->{val});
307 :     }
308 :     if ($self->{val}->[0] >= 0)
309 :     {
310 : olson 1.5 # print "write @{$self->{val}}\n";
311 : olson 1.4 $self->{cr}->write_entry($self->{cur}, $self->{val});
312 :     $self->{count}++;
313 :     }
314 :     $self->set_blank;
315 :     }
316 :     $self->{cur} = $motif;
317 :     }
318 :     $self->{val}->[$col - 1] = 0 + $val;
319 :     }
320 :    
321 :     sub finish
322 :     {
323 :     my($self) = @_;
324 :     if ($self->{val}->[0] >= 0)
325 :     {
326 :     $self->{cr}->write_entry($self->{cur}, $self->{val});
327 :     $self->{count}++;
328 :     $self->set_blank;
329 :     }
330 :     }

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3