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
$ 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
$ 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
$ 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
$ 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
$ 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
$ 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.