12
$\begingroup$

This has come up repeatedly recently: I have a very large text file (in the order of several GiB) and I need to perform line-based subsetting for around 10,000 lines. There exist solutions for specific scenarios (e.g. samtools view -s for randomly sampling BAM files) but sometimes my use-case doesn’t fit into these categories.

Unfortunately a naïve sed-based solution is extremely slow:

time sed -n -f <(awk -vOFS='' '{print $0, "p"}' line_numbers.txt) input_file > selected_lines.txt

Where line_numbers.txt is a file containing one line number per line.

Forget running this for 10,000 lines; it’s already grinding to a halt for a mere 1000.

How can I speed this up, ideally so that it scales only with the size of the input file, and has more or less constant runtime n the number of lines that I subset?

$\endgroup$
  • 2
    $\begingroup$ I think this question belongs more to stackoverflow than here, because it wants the general case to be addressed, not some specific bioinformatics file format. I tried to vote to close, but somehow, I don't have the choice to select stackoverflow as a more appropriate site. $\endgroup$ – bli Jun 5 '17 at 13:53
  • 3
    $\begingroup$ @bli I disagree, I think we should be open to any question that can arise in the course of a bioinformatician's work. Many, many bioinformatics questions can be boiled down to simple text parsing operations (consider translating between different sequence formats, for example) but I still feel those should be on topic. Also, you can't migrate to another site unless a specific migration path has been set up and we don't have any of those yet. $\endgroup$ – terdon Jun 5 '17 at 14:35
  • 2
    $\begingroup$ @bli I disagree as well, because Stack Overflow questions are pretty much always tied to a language (or the opposite, tied to none, and no code solution is required). By contrast, I'm interested in a solution but I don't care about the technology. $\endgroup$ – Konrad Rudolph Jun 5 '17 at 14:51
  • 1
    $\begingroup$ It would be helpful to include a bit more bioinformatics context / story into this question. You've given an example of what wouldn't be a use case, but not what would. I can't think off the top of my head about a situation where I've needed to subset based on the line number. It's more common for me to be subsetting based on a value in one of the columns. $\endgroup$ – gringer Jun 5 '17 at 20:05
  • 2
    $\begingroup$ With awk only: awk 'BEGIN{while((getline<"line_num.txt")>0)l[$1]=1}NR in l' input_file. $\endgroup$ – user172818 Jun 6 '17 at 2:29
7
$\begingroup$

Turns out, simply keeping track of the next candidate line (after sorting the sample line numbers) fixes the performance issue, and most of the remaining slowness seems to be due to the overhead of actually reading the file so there’s not very much to improve.

Since I don’t know how how to do this in sed, and it’s not trivial in awk either, here’s a Perl script:

#!/usr/bin/env perl

use strict;
use warnings;

my $file = $ARGV[0];
my $lines_file = $ARGV[1];

open my $lines_fh, '<', $lines_file or die "Cannot read file $lines_file";
chomp (my @lines = <$lines_fh>);
close $lines_fh;

@lines = sort {$a <=> $b} @lines;

open my $fh, '<', $file or die "Cannot read file $file";
my $line = 1;
my $next_line = 0;
while (<$fh>) {
    last if $next_line == scalar @lines;
    if ($line++ == $lines[$next_line]) {
        $next_line++;
        print;
    }
}
close $fh;

I’ve implemented a similar function in C++ for an R package, that's only slightly longer than the Perl script. It is ~3 times faster than the Perl script on my test file.

$\endgroup$
  • $\begingroup$ Wouldn't this be far faster if you used a hash instead? $line{$_}++ for @lines and then, in the while loop: print if $line{$.}. $\endgroup$ – terdon Jun 5 '17 at 13:07
  • 1
    $\begingroup$ @terdon No. Hash lookup time is constant on average, but slower than a direct index lookup in a contiguous array (plus the occasional index increment). So much for theory. I also tested it, and it holds in practice. (And, just as expected, the actual difference is almost negligible.) $\endgroup$ – Konrad Rudolph Jun 5 '17 at 13:08
  • $\begingroup$ Well, TIL. Thanks, I had always assumed hashes were faster by definition (that's what you get when you let biologists write code :). $\endgroup$ – terdon Jun 5 '17 at 13:50
5
$\begingroup$

Perl should be fairly fast with this when using a hash set to store the list of lines. A structure like this also works for subsetting based on a field value, where the comparison would be with the field rather than "$.":

#!/usr/bin/perl

use strict;
use warnings;

my $lines_file = $ARGV[0];
my %include_lines = ();

open my $lines_fh, '<', $lines_file or die "Cannot read file $lines_file";
while(<$lines_fh>){
  chomp;
  $include_lines{$_} = 1;
}
close $lines_fh;

while(<>){
  if($include_lines{$.}){ # "$." -- line number of current file
    print;
  }
}

Note that according to this SO answer, the "$." operator is not strictly the current line number, and can be influenced by different file operations or other settings.

Edit: just saw your comment about speed in your answer, comparing hash sets to a sorted list. The $lines[$next_line] bit feels a bit odd to me. Have you tried out using shift or pop on a sorted list to fetch the next line:

#!/usr/bin/perl

use strict;
use warnings;

my $lines_file = $ARGV[0];

open my $lines_fh, '<', $lines_file or die "Cannot read file $lines_file";
chomp (my @lines = <$lines_fh>);
close $lines_fh;

@lines = sort {$a <=> $b} @lines;
my $next_line = shift(@lines);

while (<>) {
    if ($. == $next_line) {
        $next_line = shift(@lines);
        print;
        last if (!@lines);
    }
}

I changed the shift to a pop (reversing the sort order to match), and got times of 4.8s, 5.8s and 4.1s for Konrad's original code, my hash code, and my pop code respectively, fetching 10,000 lines from /usr/share/dict/british-english-insane 25 times over (after copying input files to /tmp). They're all the same from my perspective: quick enough that it would take me longer to type out the command than it would to run it. Using shift instead of pop doesn't seem to change the time noticeably.

$\endgroup$
  • $\begingroup$ shift/pop modify the array which may be quite slow if memory is moved as a consequence (but I admit I don’t know whether Perl does that). It should definitely never be faster than indexed access; though to be honest it’s less code. $\endgroup$ – Konrad Rudolph Jun 5 '17 at 21:38
  • $\begingroup$ Sorry about that; I wrote these up quickly without testing. I've fixed those bugs (changed a % to a $ in line 12), and put a "my" at the front of line 18. $\endgroup$ – gringer Jun 6 '17 at 8:07
4
$\begingroup$

Some related questions appear in other sites, with potentially interesting solutions, which I report here:

To sample approximately 1% of the non-empty lines:

awk 'BEGIN {srand()} !/^$/ { if (rand() <= .01) print $0}' input_file

(from https://stackoverflow.com/a/692321/1878788)

To select 1000 random lines:

shuf -n 1000 input_file

(from https://stackoverflow.com/a/15065490/1878788, and https://unix.stackexchange.com/a/108604/55127)

Edit: Python solutions using a list of lines

Using a set of line indices and selecting lines by testing set membership:

#!/usr/bin/env python3

import sys

with open(sys.argv[2], "r") as line_numbers_file:
    line_indices = set(int(line) - 1 for line in line_numbers_file)

with open(sys.argv[1], "r") as input_file:
    print(*(line.strip() for (idx, line) in enumerate(input_file)
            if idx in line_indices), sep="\n")

Using a numpy boolean array together with itertools.compress:

#!/usr/bin/env python3

import sys
from itertools import compress
from numpy import zeros

with open(sys.argv[2], "r") as line_numbers_file:
    line_indices = [int(line) - 1 for line in line_numbers_file]

selector = zeros(max(line_indices) + 1, dtype=bool)
selector[line_indices] = 1

with open(sys.argv[1], "r") as input_file:
    print(*(line.strip() for line in compress(input_file, selector)), sep="\n")

I did some tests on a file containing 15774756 sam records and a list of 10000 pre-generated random line numbers.

The perl script proposed by Konrad Rudolph (https://bioinformatics.stackexchange.com/a/454/292) runs in about 5.3 seconds.

The set membership testing python solution runs in about 4.45 seconds.

The compress based solution runs in about 3.4 seconds. I suspect this may vary a lot depending on the highest line number we want, since the number of iterations will depend on the length of the boolean array. Here the highest line number was 15773768, so pretty high compared with the total number of lines.

I tried with python 3.6. I suspect that python 2.7 could be slightly faster, but haven't tested.

$\endgroup$
  • 4
    $\begingroup$ I think the question is about selecting a specific set of lines from the file, not just a random sample. $\endgroup$ – terdon Jun 5 '17 at 14:33
  • $\begingroup$ @terdon is right, the user cases I've come across recently required specific, nonrandom lines. Nevertheless, this is a good complement. $\endgroup$ – Konrad Rudolph Jun 5 '17 at 14:54
  • $\begingroup$ @KonradRudolph Since you mentioned samtools view -s, I thought that you were thinking of random selection for the rest of the question and that your list of line numbers was randomly generated. $\endgroup$ – bli Jun 5 '17 at 15:06
  • $\begingroup$ @bli My intention was to show this as a counter-example how I can't solve the problem. 😉 $\endgroup$ – Konrad Rudolph Jun 5 '17 at 15:10
  • $\begingroup$ int(line.strip()) can be simplified to int(line); int ignores leading/trailing whitespace anyway $\endgroup$ – Chris_Rands Jun 5 '17 at 19:01
1
$\begingroup$

I wrote a command-line (C++14) tool called subset which is up on Github: https://github.com/alexpreynolds/subset

This should be reasonably memory efficient and fast. The subset tool does not store input lines in a table, but instead streams through the file once, storing a 4 or 8k buffer chunk of the input file (depending on OS).

It stores line numbers in an array, but eight bytes per integer * 100k is 800kB, for that use case — not very much memory there.

There's an O(nlogn) sort penalty on the line number array, but again this list will be much smaller than the query file, and integer sorting is fairly optimized, so the hit should be small.

If your line number list is already sorted, I could add an option to skip sorting; let me know if that would be useful.

The filtering step walks through the line number array and input file linearly, printing lines where there are index matches, and skipping over the rest.

Indeed, subset will quit early in parsing the input file if there are no more line numbers to query. So this feature is especially useful for speeding up filtering of very large query files. (If your query file has 1M rows, say, and your last line number of interest is 12345, there's no reason to read through the rest of the file.)

You can grab, build and install it like so:

$ git clone https://github.com/alexpreynolds/subset.git
$ cd subset
$ make
$ cp subset /usr/local/bin

Once the binary is in your path, there are a couple ways to use it.

For example, you can specify a start index and length value. The following grabs seven lines starting with the 33rd line (32 as a 0-indexed value):

$ subset --prefix-with-indices -s 32 -n 7 -i query.txt > answer.txt

Or you can specify a text file containing line numbers, each on a separate line. The following reads in a file called line-numbers.txt and uses that to filter query.txt:

$ subset --prefix-with-indices -l line-numbers.txt -i query.txt > answer.txt

The indices in line-numbers.txt should be positive, 0-indexed integers. The list of numbers does not need to be sorted, as subset will sort the list of numbers for you. This is so that an efficient single pass through the input/query file can be done.

You can leave out --prefix-with-indices to leave out the debug prefix. This is there so that you can do a sanity check on the result.

The test/makefile tests demonstrate options and usages for the two types of filtering.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.