<h1>Arbitrary record access in constant time</h1>
<p>To get a random record in constant time, it is sufficient to get an arbitrary record in constant time.</p>
<p>I have two solutions here: One with <code>tabix</code> and one with <code>grabix</code>. I think the <code>grabix</code> solution is more elegant, but I am keeping the <code>tabix</code> solution below because <code>tabix</code> is a more mature tool than <code>grabix</code>.</p>
<p>Thanks to <a href="https://bioinformatics.stackexchange.com/users/37/user172818">user172818</a> for suggesting <code>grabix</code>.</p>
<h3>Update</h3>
<p>This answer previously stated that <code>tabix</code> and <code>grabix</code> perform lookups in <code>log(n)</code> time. After taking a closer look at the grabix source code and the tabix paper I am now convinced that lookups are independent of <code>n</code> in complexity. However, both tools use an index that scales in size proportionally to <code>n</code>. So, the loading of the index is order <code>n</code>. 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.</p>
<h2>Using <code>grabix</code></h2>
<ol>
<li>Compress with <code>bgzip</code>.</li>
<li>Index the file and perform lookups with <a href="https://github.com/arq5x/grabix" rel="nofollow noreferrer"><code>grabix</code></a></li>
</ol>
<p>In bash:</p>
<pre class="lang-bash prettyprint-override"><code>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)
</code></pre>
<h2>Using <code>tabix</code></h2>
<p>The tabix code is more complicated and relies on the iffy assumption that <code>\t</code> is an acceptable character for replacement of <code>\n</code> 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:</p>
<ol>
<li>Paste each record into a single line.</li>
<li>Add a dummy chromosome and line number as the first and second column.</li>
<li>Compress with <code>bgzip</code>.</li>
<li>Index the file and perform lookups with <a href="http://www.htslib.org/doc/tabix.html" rel="nofollow noreferrer"><code>tabix</code></a></li>
</ol>
<p>Note that we need to remove leading spaces introduced by <code>nl</code> and we need to introduce a dummy chromosome column to keep <code>tabix</code> happy:</p>
<pre class="lang-bash prettyprint-override"><code>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)
</code></pre>
<h1>Random record in constant time</h1>
<p>Now that we have a way of retrieving an arbitrary record in <code>log(n)</code> 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:</p>
<h2>Using grabix</h2>
<pre class="lang-python prettyprint-override"><code># 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))
</code></pre>
<h2>Using tabix</h2>
<pre class="lang-python prettyprint-override"><code># 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)
)
</code></pre>
<p>And this works for me:</p>
<pre class="lang-bash prettyprint-override"><code>python3.5 random_read.py
</code></pre>
<h1>Disclaimer</h1>
<p>Please note that <code>os.system</code> calls a system shell and is vulnerable to <a href="https://en.wikipedia.org/wiki/Code_injection#Shell_injection" rel="nofollow noreferrer">shell injection vulnerabilities</a>. If you are writing production code, then you probably want to <a href="https://docs.python.org/3/library/subprocess.html#security-considerations" rel="nofollow noreferrer">take extra precautions</a>.</p>
<p>Thanks to <a href="https://bioinformatics.stackexchange.com/users/104/chris-rands">Chris_Rands</a> for raising this issue.</p>

CC BY-SA 4.0
2017-06-19T14:06:33.253Z
2021-04-27T08:46:13.457Z
2021-04-27T08:46:13.457Z
winni2k
<p>You could shuffle the FASTQ once and then read sequences off the top of the file as you need them:</p>

<pre class="lang-bash prettyprint-override"><code>gzip -dc input.fastq.gz | paste - - - - | shuf | tr '\t' '\n'| gzip -c &gt; output.fastq.gz
</code></pre>

<p>I would recommend <a href="https://zlib.net/pigz/" rel="nofollow noreferrer">pigz</a> as a replacement for gzip in the compression step if you have it available.</p>

<p>The downside of this approach is that you only get <code>n</code> reads before you need to run the shuffle again, and <a href="https://stackoverflow.com/a/24492814/528691">apparently</a> <code>shuf</code> 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.</p>

<p>Using <code>sort -R</code> is complexity <code>n log(n)</code> and uses temporary files, so it <a href="https://vkundeti.blogspot.co.uk/2008/03/tech-algorithmic-details-of-unix-sort.html" rel="nofollow noreferrer">should not</a> run out of memory:</p>

<pre class="lang-bash prettyprint-override"><code>gzip -dc input.fastq.gz | paste - - - - | nl | sort -R | perl -pe 's/\s*\d+\t//' | tr '\t' '\n'| gzip -c &gt; output.fastq.gz
</code></pre>

<p>The <code>nl</code> and <code>perl</code> commands are necessary to make sure that identical records are not sorted next to each other.</p>

CC BY-SA 3.0
2017-06-19T13:15:22.910Z
2017-06-22T11:06:17.097Z
2017-06-22T11:06:17.097Z
winni2k
<p>One possibility is to:</p>

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

<p>This would get you the log(n) lookup time you want.</p>

<p>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.</p>

CC BY-SA 3.0
2017-06-19T13:44:42.180Z
2017-06-19T13:44:42.180Z
kvg
<p>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. <a href="http://simpsonlab.github.io/2015/05/19/io-performance/" rel="nofollow noreferrer">http://simpsonlab.github.io/2015/05/19/io-performance/</a></p>

<p>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. </p>

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

CC BY-SA 3.0
2017-06-20T17:23:42.670Z
Caching can dramatically improve query times. It would be interesting if they test out a `mmap`-based approach, over `fseek`. Especially with SSD hardware.
Alex Reynolds
CC BY-SA 3.0
2017-06-20T17:20:24.900Z
2017-06-20T17:20:24.900Z
Daniel Standage
<p>I wrote <a href="https://github.com/alexpreynolds/sample" rel="nofollow noreferrer">a tool called <code>sample</code></a> that you can use to do random sampling without reading the entire file into memory. </p>

<p>It can be used where GNU <code>shuf</code> fails for lack of sufficient memory. </p>

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

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

<p>You might use it on a FASTQ file like so:</p>

<pre class="lang-bash prettyprint-override"><code>$ sample -k 1234 -l 4 in.fq &gt; out.fq
</code></pre>

<p>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.</p>

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

CC BY-SA 3.0
2017-06-19T20:51:54.910Z
Your tool `sample` uses reservoir sampling, which is order `n` if I am not mistaken.  However, the question asks for an order `log(n)` algorithm.  This may be an XY problem, since the tool does only need one pass to retrieve `k` records (1234 in your case). However, In my particular case I do not know `k` ahead of time, which is the reason for the `log(n)` complexity requirement.
winni2k
CC BY-SA 3.0
2017-06-19T21:08:22.290Z
You can omit `-k` to shuffle the entire file, or specify it as a shell variable if derived at runtime (i.e. `-k ${K}`).
Alex Reynolds
CC BY-SA 3.0
2017-06-19T20:19:14.163Z
2017-06-22T11:07:39.083Z
2017-06-22T11:07:39.083Z
Alex Reynolds
<p>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.</p>

<ol>
<li>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.</li>
<li>If the number of reads you need is large, read the entire file and decide for each read whether or not to sample it.</li>
<li>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.</li>
</ol>

CC BY-SA 3.0
2018-03-01T05:41:23.087Z
2018-03-01T05:41:23.087Z
Edward Kirton
<p>To emphasise the issue (as outlined in the <code>1/n</code> 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.</p>

<p>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 <em>only</em> contain the start locations; any additional metadata (including sequence name) would require seeking in the original file.</p>

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

CC BY-SA 3.0
2017-06-19T15:29:47.287Z
In reference to your third paragraph: I am probably missing something, but would indexing not be order `n`?  And I believe that's acceptable as a pre-processing step in the question?
winni2k
CC BY-SA 3.0
2017-06-19T20:26:42.700Z
`n` is defined as the number of "unaligned sequencing reads" right in the first sentence.  I'm not sure how to be more clear. But please do suggest an edit now that you know what is meant.
winni2k
CC BY-SA 3.0
2017-06-19T20:30:02.527Z
I am conflating "record" with "read" in my question.  I shall make an edit to be more specific. I was adopting the wording of biopython, which talks of a sequencing record as the combination of a sequence, a list of qualities and a description.
winni2k
CC BY-SA 3.0
2017-06-19T20:33:54.723Z
"Anything that needs to read through all the lines will need at least n, so anything that does more than that has complexity greater than n" - That is not true of complexity in big O notation, which is what I am using here.  I shall make that clearer as well.
winni2k
CC BY-SA 3.0
2017-06-20T07:38:01.663Z
I really like the random seek method with block compression.  The problem I have with it is that it does not guarantee uniform sampling of records if the records are not all of identical size.
winni2k
CC BY-SA 3.0
2017-06-19T13:57:44.610Z
2017-06-20T10:13:08.820Z
2017-06-20T10:13:08.820Z
gringer
<p>Here's another approach that doesn't require any indexing, using <a href="http://bedops.readthedocs.io/en/latest/content/reference/set-operations/bedextract.html" rel="nofollow noreferrer">BEDOPS <code>bedextract</code></a> to do a <code>log(n)</code> sample on a sorted BED file. Your sample contains random records with equal probability <code>1/n</code>.</p>

<p>This approach requires a single <code>O(n)</code> pass through the file to transform it to a BED file:</p>

<pre class="lang-bash prettyprint-override"><code>$ cat records.fastq | paste - - - - | awk '{ print "chrZ\t"s"\t"(s+1)"$0 }' &gt; records.bed
</code></pre>

<p>Store the intervals in a separate file:</p>

<pre class="lang-bash prettyprint-override"><code>$ cut -f1-3 records.bed &gt; intervals.bed
</code></pre>

<p>To do a random sample of <code>k</code> elements, shuffle the intervals file and preserve the order of shuffled elements. </p>

<p>You can do this with the <a href="https://github.com/alexpreynolds/sample" rel="nofollow noreferrer"><code>sample</code> tool</a> I outlined earlier:</p>

<pre class="lang-bash prettyprint-override"><code>$ sample -k ${K} -s intervals.bed &gt; intervals-sample.bed
</code></pre>

<p>Or you can <code>shuf</code> and <code>sort-bed</code> to do the same thing:</p>

<pre class="lang-bash prettyprint-override"><code>$ shuf -n ${K} intervals.bed | sort-bed - &gt; intervals-sample.bed
</code></pre>

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

<p>Next, use <code>bedextract</code> to do a binary search on the records, and delinearize to get back to FASTQ:</p>

<pre class="lang-bash prettyprint-override"><code>$ bedextract records.bed intervals-sample.bed | cut -f4 | tr '\t' '\n' &gt; sample.fq
</code></pre>

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

<pre class="lang-bash prettyprint-override"><code>$ sample -k ${K} -s intervals.bed | bedextract records.bed - | cut -f4 | tr '\t' '\n' &gt; sample.fq
</code></pre>

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

<p><em>Note:</em> 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.</p>

CC BY-SA 3.0
2017-06-19T21:30:20.647Z
2017-06-22T11:07:13.503Z
2017-06-22T11:07:13.503Z
Alex Reynolds
<p>I wrote a tool called <a href="https://github.com/mdshw5/strandex" rel="nofollow noreferrer">strandex</a> 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 <a href="https://github.com/dib-lab/khmer/issues/1002#issuecomment-102649433" rel="nofollow noreferrer">this comment</a> about the <code>screed</code> FASTQ parser in <a href="https://github.com/dib-lab/khmer" rel="nofollow noreferrer">khmer</a>, 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:</p>

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

<p>The CLI script and python library can be installed with <code>pip install strandex</code>, and by default the script samples <code>-n</code> reads using file offsets starting at a random number (reproducable by setting <code>-s</code> 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. </p>

<p>If you specify <code>-n</code> reads less than the total number in the file you get <em>down</em>sampling, and if <code>-n</code> is greater than the total number of reads you get <em>up</em>sampling. </p>

<p>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.</p>

CC BY-SA 3.0
2017-06-20T07:46:01.800Z
Cool! Thanks for this answer. Does this tool have a way of guaranteeing uniform sampling of records with uneven sizes?
winni2k
CC BY-SA 3.0
2017-06-20T08:02:58.347Z
I suppose you could perform rejection sampling of the record with probability `1/L` that you keep the record, where `L` is the number of chars in a record...
winni2k
CC BY-SA 3.0
2017-06-22T01:46:23.913Z
@wkretzsch unfortunately it does not. The tool even fails if you have paired FASTQ files that have been trimmed to different lengths.
Matt Shirley
CC BY-SA 3.0
2017-06-20T00:43:57.250Z
2017-06-20T00:43:57.250Z
Matt Shirley
<p>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.</p>

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

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

CC BY-SA 3.0
2017-06-19T15:36:46.087Z
I think this would be an acceptable solution if it is possible to reconstruct the entire randomly selected FASTQ record (including description) from this scheme?
winni2k
CC BY-SA 3.0
2017-06-19T18:53:22.423Z
Yes. Each record will be identical in length (assuming the 'I' hack included a corresponding '-' hack for sequence), so if it is known where in file A the start of a sequence record is, then the same location will be the start of a sequence record in file B.
gringer
CC BY-SA 3.0
2017-06-19T20:56:04.393Z
@gringer, What about the read descriptions?
winni2k
CC BY-SA 3.0
2017-06-19T21:35:34.153Z
@wkretzsch Do you have any specific read in mind? If yes, could you create a gist for it? Basically all FASTQ I have ever seen had the + and @ lines identical.
Karel Břinda
CC BY-SA 3.0
2017-06-20T02:22:30.993Z
The headers for the quality and sequence would be identical in such a paired situation.
gringer
CC BY-SA 3.0
2017-06-20T07:56:48.023Z
[This read](https://gist.github.com/wkretzsch/97fd6f182ea40a4dde937316fa026905) is from the Genome in a Bottle Ashkenazi Jew trio. It has an @ line with several characters, but only a single character on the + line.  Furthermore, the assumption that all @ lines are of identical length seems [evil](http://www.cs.technion.ac.il/users/yechiel/c++-faq/defn-evil.html) to me.
winni2k
CC BY-SA 3.0
2017-06-20T09:42:38.863Z
It's not that the header lines would be the same length; it's that the header line for each read would be the same as its quality counterpart. Regardless, this also doesn't deal with the updated request for an equal probability for each sequence.
gringer
CC BY-SA 3.0
2017-06-19T15:02:15.580Z
2017-06-19T15:02:15.580Z
Karel Břinda
<p>I suppose I can provide code for this if needed, but keep in mind that it'd probably have to be done in C.</p>

<p>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.</p>

<ol>
<li>Count the number of blocks in the file.</li>
<li>Randomly select one of these blocks.</li>
<li>Perform reservoir selection on the entries in that block.</li>
</ol>

<p>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.</p>

CC BY-SA 3.0
2017-06-19T13:51:48.503Z
2017-06-19T13:51:48.503Z
Devon Ryan
<p>As <a href="https://bioinformatics.stackexchange.com/users/492/wkretzsch">wkretzsch</a> suggested this was worthy of an actual answer, I feel the obvious solution is missing here; index the <code>FASTQ</code>.</p>

<h2>Index it</h2>

<p>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 <code>samtools fqidx</code> (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.</p>

<p>So, to keep things simple - a quick and dirty alternative approach might be to use <code>biopython</code>, seeing as it <a href="http://biopython.org/DIST/docs/tutorial/Tutorial.html#htoc291" rel="noreferrer">already has this functionality implemented</a> to do this, and if installed, is trivial to use:</p>

<pre class="lang-python prettyprint-override"><code>from Bio import SeqIO
fq = SeqIO.index("myfastq.fq", "fastq")
</code></pre>

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

<pre class="lang-python prettyprint-override"><code># Random access by read name (I like owls)
record = fq["HOOT"]
record
#&gt; SeqRecord(seq=Seq('ACGTACGT', SingleLetterAlphabet()), id='HOOT', name='HOOT', description='HOOT', dbxrefs=[])

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

# and qualities
record.letter_annotations
#&gt; {'phred_quality': [41, 41, 41, 41, 41, 41, 41, 41]}
</code></pre>

<p>If you want to select arbitrary random records, you could use something like <a href="https://docs.python.org/3.1/library/random.html#random.randrange" rel="noreferrer"><code>randrange</code></a> to select between 0 and the length of the references list.</p>

<pre class="lang-python prettyprint-override"><code>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 )
</code></pre>

<p>If you want multiple records, you'll probably want <a href="https://docs.python.org/2/library/random.html#random.sample" rel="noreferrer"><code>sample</code></a> (to avoid replacement) instead:</p>

<pre class="lang-python prettyprint-override"><code>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 )
    # ...
</code></pre>

<p>For what it's worth, I think <code>biopython</code> 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 <code>FASTQ</code> once, and output the readname, file offset and length - akin to a <code>FASTA</code> <code>fai</code>.</p>

<p><strong>Update</strong>
See <code>samtools fqidx</code></p>

CC BY-SA 3.0
2017-06-21T14:27:11.310Z
``Bio.SeqIO.index`` does hold the index in RAM as a Python dict, but ``Bio.SeqIO.index_db`` puts the index in a reusable SQLite3 database file.
Peter Cock
CC BY-SA 3.0
2017-06-21T16:49:24.093Z
@peterjc Aha! That's good to know.
Sam Nicholls
CC BY-SA 4.0
2017-06-21T12:36:08.063Z
2018-08-30T11:10:07.147Z
2018-08-30T11:10:07.147Z
Sam Nicholls
<p>I would like to select a random record from a large set of <code>n</code> unaligned sequencing reads in <code>log(n)</code> time complexity (<a href="https://en.wikipedia.org/wiki/Big_O_notation" rel="noreferrer">big O notation</a>) 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.</p>

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

<p>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 <code>n</code>.</p>

<h2>Update</h2>

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

CC BY-SA 3.0
2017-06-20T07:26:24.903Z
Related question: https://bioinformatics.stackexchange.com/q/453/292
bli
<sam><fastq><reads><benchmarking>
Random access on a FASTQ file
# Random access on a FASTQ file

(Asked by: username_0 on 2017-06-19)

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](https://en.wikipedia.org/wiki/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`.

**Tags:** sam, fastq, reads, benchmarking
> **Comments:**
> **username_1**: Related question: https://bioinformatics.stackexchange.com/q/453/292  

---

## 12 answers

### Answer by username_10 on 2017-06-21. Score: 8 | Accepted answer

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

> **Comments:**
> **username_11**: ``Bio.SeqIO.index`` does hold the index in RAM as a Python dict, but ``Bio.SeqIO.index_db`` puts the index in a reusable SQLite3 database file.  
> **username_10**: @peterjc Aha! That's good to know.  

---

### Answer by username_0 on 2017-06-19. Score: 11

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

Compress with bgzip.
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:

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 > 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 by username_0 on 2017-06-19. Score: 5

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 by username_2 on 2017-06-19. Score: 3

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 by username_3 on 2017-06-20. Score: 2

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.

> **Comments:**
> **username_4**: Caching can dramatically improve query times. It would be interesting if they test out a `mmap`-based approach, over `fseek`. Especially with SSD hardware.  

---

### Answer by username_4 on 2017-06-19. Score: 2

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.

> **Comments:**
> **username_0**: Your tool `sample` uses reservoir sampling, which is order `n` if I am not mistaken.  However, the question asks for an order `log(n)` algorithm.  This may be an XY problem, since the tool does only need one pass to retrieve `k` records (1234 in your case). However, In my particular case I do not know `k` ahead of time, which is the reason for the `log(n)` complexity requirement.  
> **username_4**: You can omit `-k` to shuffle the entire file, or specify it as a shell variable if derived at runtime (i.e. `-k ${K}`).  

---

### Answer by username_5 on 2018-03-01. Score: 1

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 by username_6 on 2017-06-19. Score: 1

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.

> **Comments:**
> **username_0**: In reference to your third paragraph: I am probably missing something, but would indexing not be order `n`?  And I believe that's acceptable as a pre-processing step in the question?  
> **username_0**: `n` is defined as the number of "unaligned sequencing reads" right in the first sentence.  I'm not sure how to be more clear. But please do suggest an edit now that you know what is meant.  
> **username_0**: I am conflating "record" with "read" in my question.  I shall make an edit to be more specific. I was adopting the wording of biopython, which talks of a sequencing record as the combination of a sequence, a list of qualities and a description.  
> **username_0**: "Anything that needs to read through all the lines will need at least n, so anything that does more than that has complexity greater than n" - That is not true of complexity in big O notation, which is what I am using here.  I shall make that clearer as well.  
> **username_0**: I really like the random seek method with block compression.  The problem I have with it is that it does not guarantee uniform sampling of records if the records are not all of identical size.  

---

### Answer by username_4 on 2017-06-19. Score: 1

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 by username_7 on 2017-06-20. Score: 1

I wrote a tool called [strandex](https://github.com/mdshw5/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](https://github.com/dib-lab/khmer/issues/1002#issuecomment-102649433) about the `screed` FASTQ parser in [khmer](https://github.com/dib-lab/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 *down*sampling, and if `-n` is greater than the total number of reads you get *up*sampling.

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.

> **Comments:**
> **username_0**: Cool! Thanks for this answer. Does this tool have a way of guaranteeing uniform sampling of records with uneven sizes?  
> **username_0**: I suppose you could perform rejection sampling of the record with probability `1/L` that you keep the record, where `L` is the number of chars in a record...  
> **username_7**: @wkretzsch unfortunately it does not. The tool even fails if you have paired FASTQ files that have been trimmed to different lengths.  

---

### Answer by username_8 on 2017-06-19. Score: 0

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.

> **Comments:**
> **username_0**: I think this would be an acceptable solution if it is possible to reconstruct the entire randomly selected FASTQ record (including description) from this scheme?  
> **username_6**: Yes. Each record will be identical in length (assuming the 'I' hack included a corresponding '-' hack for sequence), so if it is known where in file A the start of a sequence record is, then the same location will be the start of a sequence record in file B.  
> **username_0**: @gringer, What about the read descriptions?  
> **username_8**: @wkretzsch Do you have any specific read in mind? If yes, could you create a gist for it? Basically all FASTQ I have ever seen had the + and @ lines identical.  
> **username_6**: The headers for the quality and sequence would be identical in such a paired situation.  
> **username_0**: [This read](https://gist.github.com/wkretzsch/97fd6f182ea40a4dde937316fa026905) is from the Genome in a Bottle Ashkenazi Jew trio. It has an @ line with several characters, but only a single character on the + line.  Furthermore, the assumption that all @ lines are of identical length seems [evil](http://www.cs.technion.ac.il/users/yechiel/c++-faq/defn-evil.html) to me.  
> **username_6**: It's not that the header lines would be the same length; it's that the header line for each read would be the same as its quality counterpart. Regardless, this also doesn't deal with the updated request for an equal probability for each sequence.  

---

### Answer by username_9 on 2017-06-19. Score: 0

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.


---

