text: user492: 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.
  
  user73: 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.
  
  user776: 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.
  
  user267: 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:
  
  seek to a random offset in the file
  read a chunk of the file
  test if the regex pattern @.+[\n\r]+.+[\n\r]+\+.*?[\n\r].+[\n\r] matches
  if so, extract the matching FASTQ record using regex capture groups
  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.
  
  user2436: 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.
  
  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.
  If the number of reads you need is large, read the entire file and decide for each read whether or not to sample it.
  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.
  
  
  user776: 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.
  
  user96: 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.
  
  user822: One possibility is to:
  
  reformat the data such that each record is a single line containing the read description, bases, and quality scores
  pad out each record to a maximum length in each field such that every record in the file is the same number of bytes
  the total number of records can now be calculated as file size / record size
  choose a random record number between 0 and the total number of records
  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.
  
  user492: 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.
  
  user215: 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
  
  user492: 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.
  
original_text: 492: <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>
   73: <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>
   776: <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>
   267: <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>
   2436: <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>
   776: <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>
   96: <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>
   822: <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>
   492: <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>
   215: <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="nofollow 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="nofollow 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="nofollow 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>
   492: <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>
  