Take the 2-minute tour ×

I am stuck with a rather unique problem. I have 2 files which I am reading. A small version of those 2 files look like the following:


chr1    9873    12227   11873   2354    +   NR_046018   DDX11L1
chr1    760970  763155  762970  2185    +   NR_047520   LOC643837


chr1    9871    0   
chr1    9872    1
chr1    9873    1
chr1    9874    2
chr1    9875    1
chr1    9876    3
chr1    9877    3
chr1    760970  1
chr1    760971  1
chr1    760972  1
chr1    760973  2
chr1    760974  3
chr1    760975  3
chr1    760976  4
chr1    760977  5
chr1    760978  6
chr1    760979  7
chr1    760980  6
chr1    760981  7
chr1    760982  8
chr1    760983  9
chr1    760984  10
chr1    760985  11
chr1    760986  12
chr1    760987  10
chr1    760988  9
chr1    760989  6


  1. From 1st file, I have to pick up the 2nd element from each row and take it as $start. An ending position is determined by $end = $start + 10.

  2. Based on $start, I now have to take the 2nd file, and look at 2nd element of each row. Once $start is found, I need to sum the next 5 corresponding values of 3rd element in groups of 5, upto $end.

So as $end is $start + 10 and I am summing in groups of 5, 2 summation values would be obtained.

In case some values upto $end is not present in the 2nd element of 2nd file, the code should not stop, it should continue to perform summation and display sum as 0 (in case a continuous group of 5 elements is not present).

Taking the example of the files here, from File1, 2nd element = 9873, which is assigned to $start. Thus $end would be $start+10 ie 9883.

From File2, once $start is found in the 2nd element of the row, the 3rd element for the next 5 rows have to be summed as 1 group, and the next 5 values summed as 2nd group till $end.


Here as can be seen in File2, $end i.e 9883 is not present. Hence sum of values from 9879 to 9883 must be zero. It must not sum the values of 760970 onwards...

Desired Output

chr1    9873    12227   11873   2354    +   NR_046018   DDX11L1      10   0
chr1    760970  763155  762970  2185    +   NR_047520   LOC643837    8   25

Points to Note

  1. While dealing with actual files, $end = $start+10,000(instead of $end = $start+10)
  2. Also,in the same note, groups of 25 values will be summed(instead of 5), obtaining total 400 values while working with the actual files.
  3. In case there are a range of values which are not present in the 2nd element of $file2, summation should proceed as normal, if a continuous pair of 25 values are absent, 0 should be printed out.
  4. The files contain > 1 million rows each.


The code I've written so far manages to do the following :

  1. Read from files.
  2. Assign $start and $end from file1
  3. From file2 , push all 2nd elements into array @c_posn ; all 3rd elements into array @peak.
  4. Check if $start is present in @c_posn

I am not able to figure out how to do the summation part. I had thought of creating a hash, where all 2nd elements of 2nd file go into keys and 3rd elements into values. But the hash is coming unordered. So I created the 2 arrays namely @c_posn for 2nd elements, @peaks for 3rd elements. But now I don't know how to simultaneously compare the 2 arrays( to ensure values of 760970 don't get summed)

use 5.012;
use warnings;
use List::Util qw/first/;

my $file1 = 'chr1trialS.out';
my $file2 = 'b1.wig';

open my $fh1,'<',$file1 or die qw /Can't_open_file_$file1/;
open my $fh2,'<',$file2 or die qw /Can't_open_file_$file2/;

my($start, $end);
    my @val1 = split;
    $start = $val1[1]; #Assign start value
    $end = $start + 10; #Assign end value
    say $start,"->",$end; #Can be commented out

my @c_posn;
my @peak;

    my @val2 = split;   
    push @c_posn,$val2[1]; #Push all 2nd elements 
    push @peak, $val2[2];  #Push all 3rd elements        

if (first { $_ eq $start} @c_posn) { say "I found it! " } #To check if $start is present in @c_posn

say "@c_posn"; #just to check all 2nd elements are obtained
say "@peak"; #just to check all 3rd elements are obtained   

Thank you for taking the time to go through my problem. If any clarifications are needed, please do ask me. I will be grateful for any and every comment/answer.

share|improve this question
Clarification needed: which items do you want to sum? Am I assuming correctly, that you want to sum the items in groups of five, $start inclusive and $end exclusive? However, this would cause the first sum for 760970 to be 1+1+1+2+3 = 8, and not 11. To get to 11, we would have to sum six values, which would be a bit weird: How should the groups be distributed then? –  amon Jan 2 '13 at 11:58
@amon You are right. The sum should be 8. I will correct the desired output. Summation needs to be performed in groups of 5. (For actual files, it will need to be done in groups of 25) –  Neal Jan 2 '13 at 12:13
add comment

3 Answers

up vote 2 down vote accepted

This is straightforward to do if b1.wig is small enough to be read into a hash in memory, taking the keys from column 2 and the values from column 3. Then all that must be done is to access each key in each sequence, using zero if a corresponding hash element is non-existent (and so accessing it returns undef).

You haven't said how you want to separate the new totals from the existing data from chr1trialS.out so I have used spaces. Of course this is easy to change if necessary.

use strict;
use warnings;

use constant SAMPLE_SIZE => 10;
use constant CHUNK_SIZE => 5;

my $file1 = 'chr1trialS.out';
my $file2 = 'b1.wig';

my %data2;
  open my $fh, '<', $file2 or die $!;

  while (<$fh>) {
    my ($key, $val) = (split)[1,2];
    $data2{$key} = $val;

open my $fh, '<', $file1 or die $!;

while (<$fh>) {
  my $key = (split)[1];
  my @totals;
  my $n = 0;
  while ($n < SAMPLE_SIZE) {
    push @totals, 0 if $n++ % CHUNK_SIZE == 0;
    $totals[-1] += $data2{$key++} // 0;
  print "$_ @totals\n";


chr1    9873    12227   11873   2354    +   NR_046018   DDX11L1 10 0
chr1    760970  763155  762970  2185    +   NR_047520   LOC643837 8 25
share|improve this answer
Hello again sir! We ‘meet’ again. I have always admired the fantastic way in which you can convert seemingly complicated problems into simple, straightforward code. :) That being said, I do need some clarifications, and would be grateful if you could spare some more of your time. 1. Why are we using the use constant pragma? What benefit it provides and can we not simply use a variable to store sample_size & chunk_size? –  Neal Jan 3 '13 at 11:31
2. Please elaborate a little on this nifty piece of code sir. I have understood some part of it, but not the whole :( while ($n < SAMPLE_SIZE) { push @totals, 0 if $n++ % CHUNK_SIZE == 0; $totals[-1] += $data2{$key++} // 0; } –  Neal Jan 3 '13 at 11:31
@Neal On 1.: constants are good for named values that will not change during execution of the script. Also, they are inlined and possibly constant-folded, which can (generally, but not quite in this case) lead to a speed increase. It is just a good habit to use constants, even here. –  amon Jan 3 '13 at 12:23
@Neal On 2.: The inner while-loop can be seen as a for-loop like for my $n (0 .. SAMPLE_SIZE - 1){...}. Whenever $n is a multiple of CHUNK_SIZE, a new field in @totals is initialized with zero. $totals[-1] always accesses the last field. If any data is available for that specific key, it is added to this last field, else a zero is added. Then, $key is incremented. –  amon Jan 3 '13 at 12:23
@Neal On 3.: split takes 0–3 arguments. If no arguments are given, it splits on any whitespace, in the default variable $_, and into as may pieces as are possible. If a first argument is given, this is a regex, determining what a seperator is. The string containing of one whitespace " " is special, and is equivalent to /\s+/. If a second argument is given, this is the string you want to split. (The third arg can limit the number of pieces you produce). So split means split " ", $_, but split $_ means split /$_/, $_, which will treat $_ as a regex, and then match on itself –  amon Jan 3 '13 at 12:29
show 3 more comments

You had the right idea with the hash. Whether it's ordered or not isn't particularly relevant because, if I understand correctly, you're looking for 11 specific values (9873, 9874, 9875... 9883), not the start value and the next 10 in the file (9873,... 9877, 760970,... 760975).

Here's how I would go about it, based on your description:

#!/usr/bin/env perl

use strict;
use warnings;

my $sum_interval = 5;   # number of lines to group into each sum
my $sum_count = 2;      # number of sums to generate
my @sums;               # final results of the operation

my %lookup;
open my $fh2, '<', 'file2.txt' or die "Can't open file 2: $!";
while (<$fh2>) { 
  my @data = split;
  $lookup{$data[1]} = $data[2];
close $fh2;

open my $fh1, '<', 'file1.txt' or die "Can't open file 1: $!";
while (my $line = <$fh1>) { 
  my @line_sums;
  my $start = (split /\s+/, $line)[1];
  for my $interval_num (0 .. $sum_count - 1) {
    my $cur_sum = 0;
    my $interval_start = $start + ($sum_interval * $interval_num);
    for (0 .. $sum_interval - 1) {
      # use || instead of // for Perl older than 5.10
      $cur_sum += $lookup{$interval_start + $_} // 0;
    push @line_sums, $cur_sum;
  push @sums, \@line_sums;
use Data::Dumper; print Dumper(\@sums);

Variable names can probably be improved on, but you can just change $sum_interval and $sum_count to 25 and 400 and it should work identically in your real application.

If the sample data you provided is placed into file1.txt and file2.txt, this produces the output:

$VAR1 = [

This output matches the results I come up with if I do the sums by hand.

Note that I deviated slightly from your spec in that it sums from $start to $start + 9 rather than $start + 10 because you said it should sum for two groups of five and $start to $start + 10 is 11 items.

Edit: Revised initial pseudocode into a complete, runnable program.

share|improve this answer
Hello Dave! Many thanks for your answer. You are right that I am looking for specific values (9873..9883). I have also put the data of file2 into a hash. However, I am unable to get the pseudocode to work. I copied the exact same code, changing $file1_handle to $fh1 for sake of brevity. When I run the code, I get error warnings like Use of unintialized value $_ in split at program.pl line 47,<$fh1> line 1 ;Use of unintialized value $start in addition at program.pl line 50,<$fh1> line 1;Use of unintialized value $start in addition at program.pl line 50,<$fh1> line 1. –  Neal Jan 2 '13 at 11:57
These warnings are repeated twice, with line 1 changing to line 2 in the 2nd set of warnings. –  Neal Jan 2 '13 at 11:59
Like I said, it was pseudocode, not actual, tested, working code, so I'm not surprised that there are warnings. And even an actual error... the split line should set $start to (split /\s+/, $line)[1]; - I forgot that you have to explicitly include the pattern to split on when splitting a value other than $_. (Answer edited to correct this.) –  Dave Sherohman Jan 2 '13 at 12:43
Thank you Dave for taking the time to answer my query. It is quite interesting to note that both you and Borodin follow largely the same idea :) –  Neal Jan 3 '13 at 12:25
add comment

Here is my current solution:


use 5.012; use warnings;

my $file1 = Reader->open("<", "filename1");
my $file2 = Reader->open("<", "filename2");

my $groupsize = 5;
my $step = 10;
my $sum_number = int($step / $groupsize) + ($step % $groupsize ? 1 : 0); # ceil($step/$groupsize)

use constant DEBUG_FLAG => 0;
sub DEBUG (@)   { say STDERR "DEBUG: ", @_ if DEBUG_FLAG }

while (my $line1 = $file1->readline) {
    my (undef, $start) = split ' ', $line1, 3;
    my $end = $start + $step;
    my @sums = (0) x $sum_number; # initialize all fields to zero
    my $i = 0;
    my $last;
    while (my $line2 = $file2->readline) {
        my (undef, $key, $val) = split ' ', $line2, 4;
        if ($start > $key) { # throw away all keys that are too small
            DEBUG "key $key too small for start $start";
        } elsif ($key >= $end) { # termination condition
            DEBUG "key $key too large for end $end";
            last LINE2;
        } else {
            $last = $key unless defined $last;
            $i += $key - $last; # get interval. This may be set to "1" as an optimization
            DEBUG "counting ($i): $sums[$i/$groupsize] + $val at $key";
            $sums[$i/$groupsize] += $val;
            $last = $key;
    DEBUG "inner loop broken";
    say join "\t", $line1, @sums; # assuming tab-seperated output

    package Reader;
    # There is probably a CPAN module for this ... :/
    use Carp;
    use constant DEBUG_FLAG => 0;
    sub open :method {
        my ($class, $mode, $filename) = @_;
        open my $fh, $mode, $filename or die qq(Can't open "$filename": $!);
        bless [$fh, []] => $class;
    sub readline :method {
        my $self = shift;
        return shift @{ $self->[1] } if @{ $self->[1] };
        my $line = scalar readline $self->[0];
        chomp $line if defined $line;
        carp "readline: " . ($line // "undef") if DEBUG_FLAG;
        return $line;
    sub pushback {
        my ($self, $line) = @_;
        carp "pushback: " . ($line // "undef") if DEBUG_FLAG;
        unshift @{ $self->[1] }, $line;
        return $self;
    sub eof :method {
        my $self = shift;
        eof $self->[0];


chr1    9873    12227   11873   2354    +   NR_046018   DDX11L1         10      0
chr1    760970  763155  762970  2185    +   NR_047520   LOC643837       8       25

This solution assumes that both input files are sorted by the second field, in ascending order, and that no overlapping sequences will be requested. If these condition can be met, it performs in constant memory, and linear time. If not, it will produce garbage, and you might be interested in using the other answer (linear memory, linear time, no restrictions). In fact, the answer by Dave Sherohman is less fragile in general, and will probably perform faster on most input.

Depending on your system, you may get a speed increase if you throw out all the object orientation, and inline the code for buffering lines (or rather, a line).

About the $i = $key - $last: The code keeps working if keys are skipped, and still adds the numbers into the right bucket. If you can assert that no keys will be skipped, or that the correct sum is irrelevant (the first five lines with IDs less than $end, not the next five IDs should be added), then removing the $last variable and simply incrementing $i by one is ok.

share|improve this answer
Regardless of whether you use an array or a hash, you should definitely read file2 into memory unless it's too big to fit. Scrubbing back and forth through the file on disk is going to be pretty slow. –  Dave Sherohman Jan 2 '13 at 12:50
I believe you also have an assumption that the ranges summed for each line in file1 are non-overlapping (i.e., you wouldn't do both 9871 and 9874). If that assumption and the stated assumptions about both files being sorted on the second field hold, then I withdraw my earlier comment, since the file would only be read in once (modulo pushing one line back at the end of each grouping), so the performance impact relative to reading the whole thing up front would be minimal-to-zero. –  Dave Sherohman Jan 2 '13 at 13:10
@DaveSherohman You have a very good point about overlapping sequences. My code does look a bit fragile in comparision. –  amon Jan 2 '13 at 13:37
@amon Thank you for your well thought out answer. But alas! I am still not yet familiar with the object oriented style of Perl :( and hence, quite a few parts of the code are a bit difficult for me , at this stage, to understand. I do hope you will forgive me for this lapse as I am still pretty much at a novice level only, got a lot of ground to cover... –  Neal Jan 3 '13 at 12:29
add comment

Your Answer


