How do I efficiently subset a very large line-based file?

Question

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?

Answer 1

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.

Answer 2

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.

Answer 3

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.

Answer 4

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.