Random access on a FASTQ file 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.
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 &quot;...a single limited transformation of the data to another file format...&quot;, 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

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

In bash:
gzip -dc input.fastq.gz | bgzip -c &gt; 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:

Paste each record into a single line.
Add a dummy chromosome and line number as the first and second column.
Compress with bgzip.
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 &gt; 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[&quot;N_LINES&quot;]) // 4
rand_record_start = random.randrange(0, n_records) * 4 + 1
rand_record_end = rand_record_start + 3
os.system(&quot;grabix grab output.fastq.gz {0} {1}&quot;.format(rand_record_start, rand_record_end))

Using tabix
# random_read.py
import os
import random

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