Random access on a FASTQ file

Question

I would like to select a random record from a large set of n unaligned sequencing reads in log(n) time complexity (big O notation) or less. A record is defined as the equivalent of four lines in FASTQ format. The records do not fit in RAM and would need to be stored on disk. Ideally, I would like to store the reads in a compressed format.

I would prefer a solution that does not require any extra files such as for example a reference genome.

The title of this question mentions a FASTQ only because FASTQ is a common format for storing unaligned reads on disk. I am happy with answers that require a single limited transformation of the data to another file format in time complexity order n.

Update

A clarification: I want the random record to be selected with probability 1/n.

Answer 1

As wkretzsch suggested this was worthy of an actual answer, I feel the obvious solution is missing here; index the FASTQ.

Index it

As much as I typically hesitate to jump to a solution that requires a script or framework (as opposed to just unix command line tools), there is sadly no samtools fqidx (perhaps there should be), and existing answers suggest a lot of munging. Whilst they probably work, some appear cumbersome and have many steps in which you could make a mistake.

So, to keep things simple - a quick and dirty alternative approach might be to use biopython, seeing as it already has this functionality implemented to do this, and if installed, is trivial to use:

from Bio import SeqIO
fq = SeqIO.index("myfastq.fq", "fastq")

Once you've acquired an index, you'll get fast random access for any read:

# Random access by read name (I like owls)
record = fq["HOOT"]
record
#> SeqRecord(seq=Seq('ACGTACGT', SingleLetterAlphabet()), id='HOOT', name='HOOT', description='HOOT', dbxrefs=[])

# We can get the sequence
record.seq
#> Seq('ACGTACGT', SingleLetterAlphabet())

# and qualities
record.letter_annotations
#> {'phred_quality': [41, 41, 41, 41, 41, 41, 41, 41]}

If you want to select arbitrary random records, you could use something like randrange to select between 0 and the length of the references list.

from random import randint

# Coerce keys to a list, as `dictionary-keyiterator` has no
#   __getitem__ attribute that would allow integer indices
# Note also that this doesn't necessarily guarantee a sorted order
#  but I guess that doesn't matter if you just want random records
key_list = list(fq.keys())

# Select a random key
# Note we use len(key_list)-1 as randint endpoints are inclusive
random_readname = key_list[ randint(0, len(key_list)-1) ]

# Get your record
rand_record = fq.get( random_readname )

If you want multiple records, you'll probably want sample (to avoid replacement) instead:

from random import sample
N = 100    
random_indices = sample(xrange(len(key_list)), N)

for key_i in random_indices:
    random_readname = key_list[ key_i ]
    rand_record = fq.get( random_readname )
    # ...

For what it's worth, I think biopython holds this index in RAM, so if your file is absolutely massive, you might need to be more clever. If that's the case, you could iterate through the FASTQ once, and output the readname, file offset and length - akin to a FASTA fai.

Update See samtools fqidx

Answer 2

Arbitrary record access in constant time

To get a random record in constant time, it is sufficient to get an arbitrary record in constant time.

I have two solutions here: One with tabix and one with grabix. I think the grabix solution is more elegant, but I am keeping the tabix solution below because tabix is a more mature tool than grabix.

Thanks to user172818 for suggesting grabix.

Update

This answer previously stated that tabix and grabix perform lookups in log(n) time. After taking a closer look at the grabix source code and the tabix paper I am now convinced that lookups are independent of n in complexity. However, both tools use an index that scales in size proportionally to n. So, the loading of the index is order n. However, if we consider the loading of the index as "...a single limited transformation of the data to another file format...", then I think this answer is still a valid one. If more than one record is to be retrieved, then the index needs to be stored in memory, perhaps with a framework such as pysam or htslib.

Using grabix

1. Compress with bgzip.
2. Index the file and perform lookups with grabix

In bash:

gzip -dc input.fastq.gz | bgzip -c > output.fastq.gz

grabix index output.fastq.gz

# retrieve 5-th record (1-based) in log(n) time
# requires some math to convert indices (4*4 + 1, 4*4 + 4) = (17, 20)
grabix grab output.fastq.gz 17 20

# Count the number of records for part two of this question
export N_LINES=$(gzip -dc output.fastq.gz | wc -l)

Using tabix

The tabix code is more complicated and relies on the iffy assumption that \t is an acceptable character for replacement of \n in a FASTQ record. If you are happy with a file format that is close to but not exactly FASTQ, then you could do the following:

1. Paste each record into a single line.
2. Add a dummy chromosome and line number as the first and second column.
3. Compress with bgzip.
4. Index the file and perform lookups with tabix

Note that we need to remove leading spaces introduced by nl and we need to introduce a dummy chromosome column to keep tabix happy:

gzip -dc input.fastq.gz | paste - - - - | nl | sed 's/^ *//' | sed 's/^/dummy\t/' | bgzip -c > output.fastq.gz
tabix -s 1 -b 2 -e 2 output.fastq.gz 

# now retrieve the 5th record (1-based) in log(n) time
tabix output.fastq.gz dummy:5-5 

# This command will retrieve the 5th record and convert it record back into FASTQ format
tabix output.fastq.gz dummy:5-5 | perl -pe 's/^dummy\t\d+\t//' | tr '\t' '\n'

# Count the number of records for part two of this question
export N_RECORDS=$(gzip -dc output.fastq.gz | wc -l)

Random record in constant time

Now that we have a way of retrieving an arbitrary record in log(n) time, retrieving a random record is simply a matter of getting a good random number generator and sampling. Here is some example code to do this in python:

Using grabix

# random_read.py
import os
import random

n_records = int(os.environ["N_LINES"]) // 4
rand_record_start = random.randrange(0, n_records) * 4 + 1
rand_record_end = rand_record_start + 3
os.system("grabix grab output.fastq.gz {0} {1}".format(rand_record_start, rand_record_end))

Using tabix

# random_read.py
import os
import random

n_records = int(os.environ["N_RECORDS"])
rand_record_index = random.randrange(0, n_records) + 1
# super ugly, but works...
os.system(
    "tabix output.fastq.gz dummy:{0}-{0} | perl -pe 's/^dummy\t\d+\t//' | tr '\t' '\n'".format(
        rand_record_index)
)

And this works for me:

python3.5 random_read.py

Disclaimer

Please note that os.system calls a system shell and is vulnerable to shell injection vulnerabilities. If you are writing production code, then you probably want to take extra precautions.

Thanks to Chris_Rands for raising this issue.

Answer 3

You could shuffle the FASTQ once and then read sequences off the top of the file as you need them:

gzip -dc input.fastq.gz | paste - - - - | shuf | tr '\t' '\n'| gzip -c > output.fastq.gz

I would recommend pigz as a replacement for gzip in the compression step if you have it available.

The downside of this approach is that you only get n reads before you need to run the shuffle again, and apparently shuf holds all data in RAM, so it would die with an out of memory error if the FASTQ file does not fit into RAM as is specified in the question.

Using sort -R is complexity n log(n) and uses temporary files, so it should not run out of memory:

gzip -dc input.fastq.gz | paste - - - - | nl | sort -R | perl -pe 's/\s*\d+\t//' | tr '\t' '\n'| gzip -c > output.fastq.gz

The nl and perl commands are necessary to make sure that identical records are not sorted next to each other.

Answer 4

One possibility is to:

1. reformat the data such that each record is a single line containing the read description, bases, and quality scores
2. pad out each record to a maximum length in each field such that every record in the file is the same number of bytes
3. the total number of records can now be calculated as file size / record size
4. choose a random record number between 0 and the total number of records
5. binary search over the reformatted file until you obtain your read

This would get you the log(n) lookup time you want.

Of course, the data wouldn't be compressed. You could 2-bit encode the bases and quantize the quals to save some space, but that'd be lossy and is perhaps not what you're looking for. Alternatively, you could block-gzip the reformatted data and keep a record of how many blocks are in the file and how many reads are in each block (since the filesize will no longer reflect the number of records in the file). Then to obtain a specific read, you'd calculate the block number it'll appear in, decompress the block, and return the appropriate read.

Answer 5

One of the most thorough treatments of this question (or a similar question: grabbing a random subset of reads) was given by Jared Simpson in a blog post a few years ago. http://simpsonlab.github.io/2015/05/19/io-performance/

If you just want to grab a single random read, Jared's benchmarks suggest that seeking to a random position in the file and then retrieving the next complete read should be the most performant option.

However, more generally, if you want to take a random subset of reads, there are many factors that affect performance.

Answer 6

I wrote a tool called sample that you can use to do random sampling without reading the entire file into memory.

It can be used where GNU shuf fails for lack of sufficient memory.

It requires two passes through the file to do a random sample, but the second pass is generally fast(er) as it uses mmap routines to do cached reads.

If you do repeated samples, the repeated samples are also mmap-ed (cached) and will run quickly.

You might use it on a FASTQ file like so:

$ sample -k 1234 -l 4 in.fq > out.fq

It parses the input file into records by every four newline characters (such as the format of a FASTQ file), reading line offset positions into memory. So the memory overhead is relatively very low.

It then applies reservoir sampling on those line offsets to write out a random sample (say, 1234 records in this example) to standard output.

Answer 7

I wrote a tool called strandex that matches entire FASTQ records starting at a random offset inside a decompressed file (since non-block compressed gzip files are not seek-able). It was a response to this comment about the screed FASTQ parser in khmer, and was initially just a proof of concept but it's been useful for me and a few other people. The general idea is to:

1. seek to a random offset in the file
2. read a chunk of the file
3. test if the regex pattern @.+[\n\r]+.+[\n\r]+\+.*?[\n\r].+[\n\r] matches
4. if so, extract the matching FASTQ record using regex capture groups
5. if no match, goto step 2

The CLI script and python library can be installed with pip install strandex, and by default the script samples -n reads using file offsets starting at a random number (reproducable by setting -s seed), then moves a number of bytes to uniformly sample the entire file, determined by the file size. In this way the process is not truly random, but random enough, and avoids sampling any part of the file too much.

If you specify -n reads less than the total number in the file you get downsampling, and if -n is greater than the total number of reads you get upsampling.

The advantage of this method over others is - it's a one-pass method, and since it's using the python C internals regex engine it's quite fast.

Answer 8

To emphasise the issue (as outlined in the 1/n update), consider an input file with one record that is 5 million bases long, and one records that is 100 bases long. You want an equal probability of selecting either of these two records. Any random seek methods will overwhelmingly pick out the long record.

I expect that indexing the locations of record starts is really the only workable option here, particularly if multi-line records are possible. Create an index file containing the start locations of each record (as identically sized integers, e.g. 64-bit), then sample from the index file (which is an identical length for each record) to fetch the start location. I'd envisage that this file would only contain the start locations; any additional metadata (including sequence name) would require seeking in the original file.

Once the indexing is done, the file can be compressed with bgzip, with specific offsets retrieved using the -b and -s options. However, I expect that compression would not be particularly efficient if multiple random records were desired.

Answer 9

Here's another approach that doesn't require any indexing, using BEDOPS bedextract to do a log(n) sample on a sorted BED file. Your sample contains random records with equal probability 1/n.

This approach requires a single O(n) pass through the file to transform it to a BED file:

$ cat records.fastq | paste - - - - | awk '{ print "chrZ\t"s"\t"(s+1)"$0 }' > records.bed

Store the intervals in a separate file:

$ cut -f1-3 records.bed > intervals.bed

To do a random sample of k elements, shuffle the intervals file and preserve the order of shuffled elements.

You can do this with the sample tool I outlined earlier:

$ sample -k ${K} -s intervals.bed > intervals-sample.bed

Or you can shuf and sort-bed to do the same thing:

$ shuf -n ${K} intervals.bed | sort-bed - > intervals-sample.bed

There's an O(klog(k)) cost here, but if k <<< n, i.e., you're working with whole-genome scale input, this cost is amortized over the log(n) search performance.

Next, use bedextract to do a binary search on the records, and delinearize to get back to FASTQ:

$ bedextract records.bed intervals-sample.bed | cut -f4 | tr '\t' '\n' > sample.fq

With Unix I/O streams, this can be done in one pass:

$ sample -k ${K} -s intervals.bed | bedextract records.bed - | cut -f4 | tr '\t' '\n' > sample.fq

By baking the sort order into records.bed, you're guaranteed the ability to do a binary search, which is log(n).

Note: Further, by linearizing the FASTQ input to a BED file and querying on BED intervals, you have equal probability of picking any one interval (interval == FQ record). You can draw an unbiased sample without the hassle of creating and storing a separate index.

Answer 10

I believe some of these answers misinterpret the original question. The OP (I believe) doesn't wish to find particular sequences, so indexing is unnecessary. I believe the OP just wants a random sampling of the dataset.

1. If the number of reads you need in small, sure you can seek on the file and return the next complete record. Remember the read IDs to sample without replacement. Jumping around on disk is slow, however, and you'll be reading 4k or whatever your disk buffer size is every time.
2. If the number of reads you need is large, read the entire file and decide for each read whether or not to sample it.
3. The way I do it (for Illumina): estimate the number of reads using the file size and read length and header length, decide how many reads I want, skip the first bunch of reads (due to edge effect, e.g. 20-100k; you can use seek or just read from start and discard), and read the next batch sequentially. I suppose it depends on what you're using the reads for; i.e. if you want random reads from your sample for biological analysis, or you want random reads from your run for QC analysis. This method is more appropriate for the former, but works well enough for the latter too.

Answer 11

I suppose I can provide code for this if needed, but keep in mind that it'd probably have to be done in C.

One could make this easier by using a bgzf compressed fastq file. Yes, BAM uses that already, but since Illumina's software now defaults to producing that there should be a reasonable amount of it already available.

1. Count the number of blocks in the file.
2. Randomly select one of these blocks.
3. Perform reservoir selection on the entries in that block.

That won't be strictly O(log N), but it's lower on memory and IO requirements than parsing through the whole file and doing shuffling or otherwise completely munging the data into an otherwise not-terribly-useful format. This also has the benefit of allowing selecting more than one random entry a bit faster, since you can memoize things.

Answer 12

A quick and dirty solution could be to convert the FASTQ file to two FASTAs, one of the them storing the bases and the other one storing the qualities, and then use standard methods for random access for FASTA.

The only complication is that base qualities can contain >. However, this is a problem only when it is the first character of a line and we can fix it by prepending, e.g., I to each quality sequence.

If also other people think that this could be a sufficient solution, I will expand the answer and add all the commands.