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

Diff of /FigKernelScripts/FFB2_create_binary_kmers.pl

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 1.3, Fri Oct 29 17:10:34 2010 UTC revision 1.4, Mon Nov 22 17:41:30 2010 UTC
# Line 21  Line 21 
21  use KmersC;  use KmersC;
22  use strict;  use strict;
23  use Data::Dumper;  use Data::Dumper;
24    use PerlIO::gzip;
25    use SortMerge;
26  use Getopt::Long;  use Getopt::Long;
27    
28  #  #
# Line 40  Line 42 
42  # It is 1-based.  # It is 1-based.
43  #  #
44  my $multiline;  my $multiline;
45    my $merge = '';
46    
47  my $rc = GetOptions("size=s" => \$column_sizes,  my $rc = GetOptions("size=s" => \$column_sizes,
48                      "multiline" => \$multiline,                      "multiline" => \$multiline,
49                        "merge=s" => \$merge,
50                      "length=s" => \$motif_len);                      "length=s" => \$motif_len);
51    
52  if (!$rc || @ARGV < 2)  if (!$rc || @ARGV < 2)
# Line 60  Line 64 
64  my $out_file = shift;  my $out_file = shift;
65  my @cols = @ARGV;  my @cols = @ARGV;
66    
67    #
68    # for merge, we don't look at input fh.
69    #
70  my $input_fh;  my $input_fh;
71    if ($merge eq '')
72    {
73  if ($in_file eq '-')  if ($in_file eq '-')
74  {  {
75      if (@lens == 0 || !defined($motif_len))      if (@lens == 0 || !defined($motif_len))
# Line 100  Line 109 
109    
110      seek($input_fh, 0, 0);      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        }
119    }
120    
121  if (!@lens)  if (!@lens)
122  {  {
# Line 136  Line 153 
153    
154  my $count = 0;  my $count = 0;
155    
156  if ($multiline)  if ($merge)
157  {  {
158      my @blank_val = map { -1 } 1..$num_keys;      my @code;
159      my @val = @blank_val;      my @fh;
160    
161      my $cur;      if (!open(MF, "<", $merge))
     while (defined(my $l = <$input_fh>))  
162      {      {
163          if ($l =~ /^(\S+)\t(\d+)\t(\d+)/)          die "Cannot open merge file $merge: $!";
164        }
165    
166        my %by_col;
167    
168        while (<MF>)
169          {          {
170              my($motif, $col, $val) = ($1, $2, $3);          chomp;
171              if ($cur ne $motif)          my($file, $col) = split(/\t/);
172            my $fh;
173            if ($file =~ /\.gz$/)
174              {              {
175                  if (defined($cur))              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                  {                  {
180                      if (@val != $num_keys)              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                          die "Val became invalid: " . Dumper($val);          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                {
195                    close($fh);
196                    shift @$fhlist;
197                    if (@$fhlist == 0)
198                    {
199                        return ();
200                      }                      }
201                      if ($val[0] >= 0)                  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                          $cr->write_entry($cur, \@val);                  my($motif, $val) = ($1, $2);
210                          $count++;                  return ($motif, [$motif, $col, $val]);
211                }
212                else
213                {
214                    die "Badly formatted input at line $. file $file: $l\n";
215                      }                      }
216                      @val = @blank_val;          };
217                  }                  }
218                  $cur = $motif;  
219        my $writer = MultilineWriter->new($cr, $num_keys);
220    
221        Sort::Merge::sort_coderefs(\@code,
222                                   sub {
223                                       my($code, $key, $value) = @{$_[0]};
224                                       #print "@$value\n";
225                                       $writer->push_one_value(@$value);
226                                   });
227        $writer->finish();
228    
229        $count = $writer->{count};
230              }              }
231              $val[$col - 1] = 0 + $val;  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          }          }
241          else          else
242          {          {
243              warn "Badly formatted input at $.: $l\n";              warn "Badly formatted input at $.: $l\n";
244          }          }
245      }      }
246      if ($val[0] >= 0)      close($input_fh);
247      {      $writer->finish();
248          $cr->write_entry($cur, \@val);      $count = $writer->count();
         $count++;  
     }  
249  }  }
250  else  else
251  {  {
# Line 194  Line 264 
264  $cr->close_file();  $cr->close_file();
265    
266  print "Loaded $count oligos\n";  print "Loaded $count oligos\n";
267    
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        if ($self->{cur} ne $motif)
300        {
301            if ($self->{cur} ne '')
302            {
303                if (@{$self->{val}} != $self->{num_keys})
304                {
305                    die "Val became invalid: " . Dumper($self->{val});
306                }
307                if ($self->{val}->[0] >= 0)
308                {
309                    #print "write @{$self->{val}}\n";
310                    $self->{cr}->write_entry($self->{cur}, $self->{val});
311                    $self->{count}++;
312                }
313                $self->set_blank;
314            }
315            $self->{cur} = $motif;
316        }
317        $self->{val}->[$col - 1] = 0 + $val;
318    }
319    
320    sub finish
321    {
322        my($self) = @_;
323        if ($self->{val}->[0] >= 0)
324        {
325            $self->{cr}->write_entry($self->{cur}, $self->{val});
326            $self->{count}++;
327            $self->set_blank;
328        }
329    }

Legend:
Removed from v.1.3  
changed lines
  Added in v.1.4

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3