We're wasting money by only supporting gzip for raw DNA files
bioinformaticszen.com
bioinformaticszen.com
fastqz (not the best name!) supports reference compression too, for smaller files when they're reads compared to a reference genome.
Here's an 2013 paper considering the problem: https://journals.plos.org/plosone/article?id=10.1371/journal...
Still, zstd is widely available and a simple drop-in replacement for gzip.
Fascinating to think that in principle, genome assembly and fastq file compression have many parallels, and are essentially the same problem. For example, once you assemble the genome, you can just track the coordinates and 'diff' for each read and the original FASTQ file can be regenerated exactly.
Although I'm not familiar with all the details of the method, my understanding is that the reference used by compression can be extremely low quality for biological purposes while still providing good compression -- making it possible to make a quick-and-dirty reference with much less computational time.
Especially if major portions of the files are repeated.
(And of course the label on each entry is just standard 7-bit ascii), and also very compressible, but it contains higher complexity text.
It depends: for some platforms (ex: NextSeq) the quality score in binned and highly compressible, but for others (ex: Nanopore) it's much higher entropy than the sequence.
ASCII in the strict sense is 7-bit, CGAT is 2-bit.
I'm counting ASCII as 8 bits because it's still stored as 8 bits regardless of whether you only use 7 of them.
In my mind "always threading, everywhere" is reason enough to always use zstd. Sure there are implementations for other compression algorithms that support threading but the reference zstd implementation from Facebook always does, everywhere.
That seems clear enough.
But to answer your question about "one virus", the initial isolate of SARS-CoV-2 from Wuhan, China has a genome containing 29,903 nucleotides. You can see the full raw sequence by clicking on "FASTA" on this page: https://www.ncbi.nlm.nih.gov/nuccore/1798174254
I agree for sequence reads, zstd is better than gzip in almost every aspect. The only problem is that zstd is not as widely available. Many Linux distros don't ship zstd by default but on managed clusters, common users don't have the root permission to install zstd system-wide.
It is worth mentioning that many mainstream bioinformatics tools can read FASTQ over pipes. You can use something like
mytool <(zstd -dc in.fastq.zst) > output.txt # bash; also works on multiple zstd input files
zstd -dc in.fastq.zst | mytool - > output.txt # general; for multiple input files, use mkfifo
[1] https://www.biorxiv.org/content/10.1101/2022.04.07.487441v1
Wikipedia has a page on the topic: https://en.wikipedia.org/wiki/Compression_of_genomic_sequenc...
I'm curious because in the past I've seen people are very cost sensitive but also want to keep lots of data that seems to go mostly unused.
I also think there’s a certain degree of fear in not being able to generate the data again if you needed to, due to lack of funding or the organism itself not being available.
Considering sequencing costs alone (ignoring costs of obtaining & processing a biological sample), Illumina's current cheapest sequencing kits - Novaseq S4 300 cycles - cost ~$15k for 3000 Gbases of data (whole-genome sequencing for about 24 people). As gziped fastqs, that data will occupy about 1.3TB of space.
Extremely high-volume purchasers may be able knock that sequencing price down some, but the storage challenge starts to seem less intimidating when you realize spending 5% of your sequencing budget on storage would give you a budget of $750/TB of data.
*In the absence of read trimming or discarding unmapped reads
With 8 bins, the variant calling accuraccy seems to be preserved, while drastically reducing the file size.
[1]: https://en.wikipedia.org/wiki/CRAM_%28file_format%29
[2]: https://lh3.github.io/2020/05/25/format-quality-binning-and-...
But I agree with you: it really depends on the type of the data.
(The paper I linked describes more or less what I'm currently working on)
Perhaps an interim solution would be something like "zstd+metadata", where the metadata is sufficient to transparently and efficiently reconstruct a gzip-compressed file on-demand. (similarly to how JPEG-XL allows oldschool JPEGs to be recompressed without losing any data)
This could have a bit of compute overhead, but since gzip decompression is so single-threaded I think you could do the conversion in parallel without actually slowing down the hot path. So, the performance would be approximately the same as using gzip (assuming decompression is compute-bound, not io-bound), but with all the storage benefits of using zstd.
Almost all new tools support gzip by default, because almost all public data is gzip-compressed. Developers often don't bother adding zstd support, because there is little zstd-compressed data around. Because the labs developing new tools are rarely the ones that have to store petabytes of data, they don't benefit from zstd directly. And once the active development is over, most tools receive only minimal support, because nobody wants to pay for maintaining the long tail of average tools.
Also, zstd is the latest reason why gzip is still widely used. Every time someone releases a new general-purpose data compressor that becomes popular, gzip gains a few more years. Because there have been too many gzip replacements over the years, none of them has managed to become popular enough to actually replace gzip before the attention moves to the next state-of-the art compressor. If people agreed this time that zstd is what everyone will be using for the next 20 years, bioinformatics would probably move to it in a few years.
Furthermore, if you want access to be truly transparent (e.g. something like a FUSE fs that transparently presents a set of gzip files, each backed by a corresponding zstd file on disk), then you (may) need to support seeking to an arbitrary offset. As far as I know, zstd doesn't support seeking into a stream by default, so you'd have to add some metadata to facilitate that.
It is possible that you don't actually care about this, in which case yes, you don't need any extra metadata.
Next thing you have is a set of tools to recreate a small subset of SQL, to index the file, to add in bulk, to edit the metadata...
The typical VCF has data enough to be a SQLite, and nobody parses the VCF directly but with tools.
This ends in a sad number of bio-scientists that cannot do the simplest SQL query, but know perfectly vcftools, samtools, bedtools and others (or have them hardcoded in shell scripts). Those formats start so simple you can "parse" them with grep, cut, wc and paste, but soon they need special tooling and get feature creep.
</joke>
[0]: https://manpages.debian.org/bullseye/zstd/zstdmt.1.en.html
Due to higher decompression cost, the opportunity here seems localized to long-term storage. It feels like it would make more sense as a project (or product!) that implemented an efficient long-term archive (perhaps with a less compressed LRU cache in front).
When I can get away with -16 for large file long term storage, I use it.
I was thinking something like the possibility to seek through a compressed file (after having built an index). Your DNA file is probably stored on a shared folder somewhere (if it's big you probably don't want to copy it to every workstation) and you point your software to the gz file directly, it create a seek-index on the first use, and then when you need view only a small section of the file you don't have to extract the whole file.
Some more advanced compression like zstd might use some form of dictionary to do the compression/decompression, which may be bigger (up to 32MB? max dictionary size for zstd whereas lz77 used by gzip use a dictionary based on a sliding window of the last 32K token) which you may have to transfer.
The intent was to enhance this to store alignments against multiple references as new references are published, and probably to rewrite in Rust or C rather than the initial Python version.
In retrospect I would be interested to know whether this domain-specific compression effort, with zstd to the resulting "hybrid" alignment, would be more efficient than just letting zstd do its own thing with a full set of individual alignments against the different references.
zstd has delta / patch mode, which creates a file that stores the "patch" to create a new file from an old (reference) file. See https://github.com/facebook/zstd/wiki/Zstandard-as-a-patchin...
See the man page: https://github.com/facebook/zstd/blob/dev/programs/zstd.1.md
More general (pre-alignment) sequence compression methods never really took of (e.g. https://github.com/BEETL/BEETL). Probably because it helps so much to have common format that most workflows can start with. Here, the replacement of gzip with zstd would be a lower hanging fruit.
`zstd` has a `--patch-from` mode, which does essentially the same thing, computing a new file using another file (typically an older version) as a reference. When similarities are larges, it leads to a huge reduction in compressed content.
I wonder what would be the performance of `zstd --patch-from=` when employing the same reference as these specialized DNA compressors.
e.g. the RLE strategy is good for data with many compressible patterns, while the filtered or Huffman-only strategies are good for data with few compressible patterns.
The default strategy is a balance between the two extremes, which isn't tuned for a specific kind of input.
For smaller players where maximum tool interoperability is important, gzip seems good enough.
Pretrained dictionaries could also be useful too.
You're correct that there's a lot of redundancy in these files, not least because each area of the genome is generally covered several times over. Removing this redundancy (or using it for error correction) is a big part of the downstream processing of this file.
Longer higher quality reads have been developed and for them less coverage will be needed and this will produce smaller files with higher degrees of confidence.
For something like a liquid biopsy test that's looking for circulating tumor DNA, a given mutation may only exist in 0.5% of the reads that cover a particular locus. Then you have to deduplicate to tell which are PCR duplicates of the same original strand, and which have different progenitor fragments. Low-coverage long reads are bad for this.
This isn't a diff from a reference genome, it's a ton of reads (files often 30+GB). It's really only useful as an archive format, and even then, it's a pretty bad one for many reasons.
An interesting thing about the file format is that it's three lines in a row each with different coding properties, so if you break the data into three streams and encode each one on its own, the overall compression is somewhat better, because the dictionary and lookback don't have to accomodate the metadata, DNA, and error lines together, giving better dictionary and lookback hit rates.
These formats are terrible for downstream analysis: they are typically not indexed or sorted in a useful way, so it's common to convert them to other formats (BAM, CRAM), which are optimized for compressing DNA data (and likely exceed zstd in some circumstances).
FASTQ are raw data reads. BAM and CRAM are reads mapped against a genome.
For those not aware one can, in general, use a predictive model to compress data using arithmetic coding.
No, it’s 100% longer, twice as long, and the inverse is 50%
If anything, I see the opposite problem much more frequently: People who cheap out on the sequencing depth, neutering their statistical power!
Maybe you don't need 500x coverage, and if you don't, then I'm glad you're not sequencing to that depth. But please don't tell other people they're doing something wrong when you don't understand their use cases.
Making the blanket statement that 500x coverage is generally not needed for any experiment is patently false.
I've built clinical liquid biopsy tests where it was a QC failure if particular sites of interest were under 5000x coverage. Sequencing such sites to only 500x coverage would be a waste of everybody's time, because it would be impossible to gain meaningful results from so little data.