Are these files already just diffs from baseline/standard full DNA scans, or are they including all the data? Seems like most scans will only differ a little…
These are FASTQ files[0], which are the actual raw bases read by the sequencer, and a quality score for each base. This is upstream of comparison to a known genome any other processing.
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.
These snippets can be used to reconstructing the whole genome sequence but it's done statistically so they're not 100% confident. Some sections have more coverage than others. You want to keep the raw reads around just in case.
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.
Depends on the use case. If one is trying to broadly sequence a genome or exome, all of that is correct.
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.
The whole point behind Next Generation Sequencing is to get a lot of redundant and cheap raw data, but with some errors. This way you sequence the same base of the genome say 10 times or more. If 9 of them are the same as the reference, and 1 is another base (and this base is low quality), the Variant Caller can be confident that was a sequencing error, and the position is the same base as the reference. But if you have 5 reference and 5 of other base, it's probably heterozigous for that position. If you have more of other base, is homozigous for that other base.
They are FASTQ files, so they are the reads that came off the sequencer- plus metadata and error scores- represented as ASCII A, T, G, C.
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).
> 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.
BAM/CRAM reads are often mapped against a genome, but they needn't be. It's not uncommon in industry to convert FASTQ to unmapped BAM in order to only have a single file format, improve compression, etc.
Yep, I've been advocating zstd since I started working in the field. It compresses and decompresses so much faster than xz and is much much more compact than gzip.
All public SARS-CoV-2 consensus sequences are ~300GB uncompressed. If you compress using gzip you end up with ~30GB. With xz/zstd you can get it down to ~2GB. However xz takes ~40min to uncompress, whereas zstd can do it in ~8min.
Because large companies are willing to pay people market salaries to work on this stuff, and then give it away? Isn't that what we want, a model where free software development is properly compensated?
Free software gets made, to the degree that it benefits Facebook, in the ways which benefit Facebook. Software development in areas which don't benefit Facebook suffers. Facebook gets people to use their software, and therefore gains influence in places where they don't belong. Yes, the developer gets paid, but it's not a good solution that's healthy for an independent free software movement.
It should be no surprise that Facebook works on projects that benefit Facebook. By making it free software, it also happens to benefit the wider world at the same time. Win-win.
Is your problem with the idea that developers are getting paid at all, or just that this instance is at Facebook?
My main issue is that I would rather not use Facebook products, and I would rather we don't create standards which rely on Facebook products.
The FOSS funding issue I more of a side-point which I didn't even bring up, but I don't see "get hired by Facebook if you want to work on FOSS" to be a healthy solution to the FOSS funding problem, we need a way to ensure funding for FOSS developers who don't work for Facebook. But that's not strictly related to the zstd thing specifically.
How is it begging the question? Is it the word "product" you're unhappy with? I can call it a Facebook project instead if you prefer, it doesn't change the substance of anything I've written. The problem I have is with the "Facebook" part.
TBH, I worry a lot more about this kind of thing in the case of Google, where they can afford to give away things that other companies would have made an entire business out of, which devalues those markets, then they lose interest and sunset it all— and only in a few cases (like Wave) does any of it end up open source.
But FB isn't really competing with other companies here. They're giving away frameworks, libraries, tools— this is building block stuff, and the fact that it's perhaps particularly optimized for FB's workloads is just their privilege, having been the ones to fund it.
A lot of consumer surplus is derived from companies trying to undermine their competitors by commoditizing non-core competency. To not take advantage of that is to have to bad parts of capitalism without the good parts.
I didn't say Facebook shouldn't do what they need to reduce their own network traffic. They don't belong in our file systems and our gene sequencing standards and our non-Facebook web services.
It's completely open-source, general purpose software that is essentially static (lower chance of drama or ideological/management concerns ala golang).
It would be like avoiding C because it was developed by ATT.
If AT&T still controlled C compilers and the C standard and and contributions to the C compiler everyone uses was gated behind AT&T I might strongly have considered avoiding C.
If a good non-Facebook implementation of ZSTD shows up, I might strongly consider using (that implementation of) ZSTD.
I simply don't think it matters. This isn't golang/react. Features, bugs and upgrades don't move with nearly the velocity of something like a large, opinionated library or language. But I rest my case.
1. There are several implementations of Zstandard other than the "reference" libzstd that we maintain, tagged as "Ports" in the table here [0]. You are of course welcome to create or fork your own. We will happily take a PR adding yours to the table.
2. The format spec is an RFC [1]. It's frozen. We're not going to pull the rug out from under you somehow (and why would we??).
3. I think we do a pretty good job accepting contributions. But again, if we aren't, you can always fork it.
All in all, we've worked pretty hard to make zstd universal--portable, stable, flexible, robust, etc. I guess I don't really have a point... I guess I hope you reconsider your rejection of zstd? It's pretty cool.
This is going to sound like I'm moving the goal posts, because, well, I am. I apologize. But what I notice with all those ports are that they're all in languages other than C (or C++), which means they're pretty much confined to their language's ecosystem; there is no port that's prepared to take the place of the Facebook zstd library. If I want to write, say, Common Lisp bindings, I'm not going to write them against the Java port. The one which comes the closest is probably the Rust port, since Rust is establishing itself as a core low-level language which people write bindings to, but the Rust port is only a decompressor.
I'm also concerned with the term "port". I don't know exactly what you put in the term, but to me, it sounds like someone basically translated the C code directly to other languages. What I would be looking for is more like a "clean re-implementation" by an independent group of people who properly understand all the algorithms which are used and where all the code is written from scratch without the Facebook code as a "crutch". This would ensure that whatever Facebook does, there will be a group of people that is capable of fixing tricky bugs in the algorithms or implementing further optimizations. Basically, I would want the alternative implementation to be as unrelated to Facebook's libzstd as Clang is to GCC.
All of that said though, from a purely technical perspective, zstd seems like a truly great piece of engineering. Kudos. Though I will never understand what drives great developers to dedicate their life to producing value for Facebook.
I have some internal conflict about zstd. On the one hand, it's truly become the superior general purpose compression codec, and I hate to look a gift horse in the mouth there. Using it for everything from data-at-rest compression on filesystems and zswap to channel compression with rsync, it truly works wonders. And yet, I can't really deny how absolutely at odds I am with the business model and actions of the company that has sponsored its development. Personally, I've resolved that we should extract as much good from things I perceive as evil. I've stopped using Facebook platforms myself, so I'm doing as much as I can to reduce the amount of money I generate for the platform as I can. That's about the most "voting with my dollars" I can do. However, because zstd has been such a net good, I will continue to use it in projects going forward for compression needs.
Petty perhaps, but I'm really annoyed that he didn't follow the celebrated command line tool tradition of compress, gzip, xz, lzip, bzip2, etc. I always have to remember that zstd behaves differently and you have to take extra steps to achieve the same thing.
What differences do you have in mind? Zstd is largely compatible as a drop-in replacement for those commands. There are differences, but it's our goal to minimize them.
I'm a Facebook hater too, but I can't think of any even theoretical concerns around zstd. It's freely licensed, no patent concerns, it's straight up free software. You're not giving them money or supporting their ad model by using it. It's not like Chrome where Google has an interest in making web browsing suck.
What theoretical path to harm is there to adopting it more widely?
I think the best way I can describe my discomfort is this: If someone asked you: "What organization do you want to be the owner of, and own the copyright to, the compression library which underpins all our open standards and systems?", what would you have responded? I don't know for sure which organization I would've chosen, but Facebook would certainly be at the very bottom of the list. GNU and Mozilla might be closer to the top.
Zstd doesn't use a login, it doesn't track you across the internet, it doesn't make shadow profiles on you, and it doesn't have ads or telemetry. Zstd is even part of the Linux kernel at this point[0], so it's safe to say that it isn't a Facebook product anymore. They can't pull the rug on zstd, since it's BSD/GPL licensed. Unless I'm missing something, nobody is feeding the beast by using zstd.
There's lots of similar cool technologies made by other "evil" companies. I don't use ZFS, but its origins in Sun/Oracle isn't the reason I don't.
I don't use ZFS, and the Oracle origins is the reason why. I'm very worried that the best filesystem available is forever stuck with a hostile license, and that so many of the good filesystem engineers are contributing to Oracle's product rather than to a non-Oracle alternative. (EDIT: I know zstd is under a good license, so its situation isn't the same a ZFS's, I just commented on ZFS since you brought it up.)
> that so many of the good filesystem engineers are contributing to Oracle's product rather than to a non-Oracle alternative.
OpenZFS isn't Oracle's product; it certainly includes a great deal of Sun/Oracle code, but Oracle themselves can't use the OpenZFS commits without honoring the CDDL license, which to my knowledge they've not done.
ZFS was made by Sun Microsystems, which was a dear, awesome company that educated, fed, clothed and put bread on the table for an entire generation of computer professionals outside of Sun Microsystems, not to mention paid millions out of their own coffers to open source a reference implementation of a System V Release 4.0 UNIX. And to top it all off, Sun Microsystems made awesome hardware.
Sun kicked butt
had fun
didn't cheat
loved their customers
changed computing forever.
I hear people say that, yet Sun Microsystems is the company which created the CDDL to be intentionally GPL-incompatible[1] and then made a filesystem which they licensed under the CDDL. Then they sold themselves to one of the worst tech companies. Not exactly the behaviour I would expect for an awesome company.
It's a bit more complicated than that. Yes, the code that's out right now under the GPL will always be licensed under the GPL. But Facebook can at any time change the license of the zstd project, meaning all future releases will be under a different license. There will be an opportunity to fork the code base pre license change, but those forks are always risky and there's a very large change that the ecosystem won't end up with a de facto standard zstd fork.
They're big enough that it shouldn't help. Dictionaries are useful when you're compressing tiny files such that the overhead of the compressor learning the statistics from the input is significant. Basically you hope to save a fraction of the size of the dictionary from each file. From a 2GB file that's negligible. From a 1000 byte file it could be gigantic.
Using a custom dictionary works well if you want to compress single consensus sequences (30kbp), say before putting them into an sqlite db. Then you can get each sequence down to ~200 bytes.
If you compress a large dataset, there's no point in using a custom dictionary, the auto-embedded one works just fine.
Decompression speed is one of zstd's biggest advantages (in its class, lz4 is obviously faster in return for worse compression ratios). Zstd decompression is generally at least factor 3 faster than gzip.
The bitstream format for Gzip files (PKZIP's "deflate") was established in 1993. In comparison, the much newer Zstd doesn't represent a different tradeoff between compression time / decompression time / size; it's just better overall. There's a nice graph of the different time/compression tradeoffs both can achieve about 1/3rd of the way into https://engineering.fb.com/2016/08/31/core-data/smaller-and-...
To put that in perspective, consider that a large portion of the file is literally a nucleotide sequence, e.g. "GGGTGATGGCCGCTGCCGATGGCGTCAAATCCCACC" and that the rest of the file format is also fairly repetitive. If we just imagine only the nucleotide sequence, one way to compress that would be to take advantage of going from an 8 bit alphabet (ASCII encoding) to a 4 bit alphabet (CGAT). Naively, you'd expect ~50% compression from that alone, and be left with a bit stream that's still compressible.
Yeah, it would seem to me (as an uneducated moron looking from the outside) that if you have domain knowledge about what is being compressed, a specialized compression algorithm would be much better overall than a generic one.
Especially if major portions of the files are repeated.
Oh I wish it was easier to build these custom compression algorithms (or to help the existing ones). Whenever you know something about the file format, content, ranges of data, even sequence information...
Yes, two bits (assuming you only need to encode those four bases; unfortunately, biological reality requires FASTQ to encode other nucleotide values). There's even a format from UCSC called "2bit" (https://genome.ucsc.edu/goldenPath/help/twoBit.html)
Unfortunately, however, FASTQ files have more characters than ACGT and (N). They include a second line that encodes sequence quality. In practice, that line is even more compressible, but it complicates things.
(And of course the label on each entry is just standard 7-bit ascii), and also very compressible, but it contains higher complexity text.
> In practice, that line is even more compressible, but it complicates things.
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.
Wouldn't CGAT be a 2-bit alphabet, bringing file size to 25% of the original? Assuming you go for a binary format with length headers, so you don't need to reserve control characters.
Yeah, 2 bit vs 4 bit is my bad. It doesn't really change much though, except to strengthen the conclusion and just makes gzip's 10:1 compression ratio even less impressive.
I'm counting ASCII as 8 bits because it's still stored as 8 bits regardless of whether you only use 7 of them.
Apart from CGAT being a 2-bit alphabet, changing the alphabet does not change the information density. Expect that this kind of transformation has no impact on the compressed result, with most general purpose algorithms.
Same. If nothing else zstd universally supports compression and decompression with threading across any number of cores. It will even adjust threads dynamically based on the next limiting factor (I/O speed). Meanwhile it's also generally faster and provides significantly better compression ratios (as you describe).
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.
The 300GB is for "all public sequences", as explained.
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
For compressing similar genomes, specialized tools will do better. AGC [1] can compress 620k SARS-CoV-2 genomes (19 GB uncompressed) into tens of MB. More importantly, AGC allows you to extract individual sequences without decompressing the entire file. In addition, zstd/xz wouldn't have a big enough dictionary to compress multiple human genomes efficiently. AGC can compress 300 GB human assemblies to 1.5 GB, a 200-fold compression ratio.
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
Fair enough. For the relatively few tools that don’t support streaming, you may generate temporary plain files in your pipeline. Spades itself is a pipeline and writes temporary files. It is anyway not realistic to expect wide adoption of zstd in the near future. If you want to save storage cost now, you have to take some trade off.
I've never seen a 15x difference in density between gzip and zstd on genomic data in an apples-to-apples comparison; I'd be really curious to learn more about how you achieved that.
Glad to see zstd getting some love here. At a previous job I needed to capture and store large amounts of JSON about social media. I did an extensive compression bake-off taking into account both compression and storage costs, and for us the winner was zstd with the compression dial turned up a fair bit. (Sorry I don't have numbers here, but they're all in the hands of the previous employer. But I think we settled on -13 as cost-optimal if we were storing data on AWS for a year.)
I'm unfamiliar with the bioinformatics scene, but I can imagine there's a lot of "legacy" hardware and software lying around that can't easily be updated to support superior formats.
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.
Yeah that's exactly the issue. A lot of tools are written once and rarely updated, they're considered "done". And they only support gz out of the box because they were written 10 years ago.
I think the tools people actually use are rarely that old.
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.
But that's hardly an issue, since you can just add the new compression tool to your pipeline before and after the older tool. Many tools don't even support gzip because they assume the user will be piping in uncompressed data from zcat anyway.
What metadata do you mean? Since both zstd and gzip are lossless compression tools all you need if you have an old tool that only takes in gzipped fasta would be "cat foo.fasta.zstd | zstdcat | gzip | oldtool".
Compressing gzip is (relatively) slow if you want to actually compress the data. The biggest slowdown comes from searching for the best backreference at any given point. If you could encode some relatively compact data saying "here's where the good backrefs are", you could optimise your throughput/compression-ratio ratio (at the cost of slightly more data on-disk).
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.
Processes like alignment are alot slower than compression. This is a space where some tooling takes days or weeks to run, and you are running them on a cluster of xeons that can gzip or gunzip fast enough for your use cases anyhow.
The bioinformatics scene is generally legacy in the good way, like how bash tooling is legacy: consistent and highly modular. Sure, the tool you are using might not be able to support compressed data or some new compression algorithm, but then so what, you can just add whatever hotnewcompression tool before and after this tool in your pipeline, and throttle the work to maximize available storage and compute. Many pipelines are written simply in bash, or use something like snakemake or nextflow to ensure consistent environments.
A much bigger problem than the compression format is the fact that biologists over-sequence. You don't need 500x coverage, trust me, please learn to design your experiments.
That's a pretty expansive statement. Are you doing error-corrected sequencing to look for rare mosaic variants? Do you care about transcripts expressed at very low levels? Are you trying to tease apart subclonal variation in cancer? All of these could be good reasons to sequence deeply.
If anything, I see the opposite problem much more frequently: People who cheap out on the sequencing depth, neutering their statistical power!
> You don't need 500x coverage, trust me, please learn to design your experiments
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.
I would assume most datasets people are using aren't large enough to justify using a database versus parsing a flat file. You also have to realize this is a field of scientific code, where it matters more that you spend less time coding and more time interpreting results, versus spending time optimizing the pipeline to minimize compute time, and you might be working on a university cluster where your compute is powerful and quite cheap. For people who might work on clinical pipelines that will continue to be reran time again over a vast growing amount of patient data, they probably already put their data into databases. For your academic post doc working on 2000 samples from an experiment for one paper before they find another job doing something else entirely in two years, a flat file is fine.
And it’s much easier to teach a grad student who already knows some basic bash, R, or Python how to read a flat file into a data frame and make some plots or grep a few lines versus dealing with databases. While bioinformatics tooling is outdated in many ways, modern software engineering could do with more config.txt and fewer hidden SQLite databases holding settings only accessible by an electron GUI.
I was thinking for example in VCF files. A metadata header, a main table with eight clear columns and a ninth column that works as a "put here whatever you need", and then the related data for each sample in extra columns.
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.
It's just the letters A, G, T, and C over and over again, right? surely we can squeeze some better compression out of that, dictionary coding or something?
As someone who helped build a fastq data processing flow for a genetics company during the sequencing price drop era, iirc we had issues with data corruption with some of the other formats in tests. Funny to think, hey, I got my first taste of very large data management because of this!
There are compressors specialized for FASTQ that are faster and denser than zstd. FASTQ is the the most common format for storing DNA sequencing data-- a text file including metadata, sequence, confidence scores, etc. http://kirr.dyndns.org/sequence-compression-benchmark/
fastqz (not the best name!) supports reference compression too, for smaller files when they're reads compared to a reference genome.
Reference compression is interesting... technically, this type of compression should be possible without a reference. E.g. in a worst case scenario, you do a genome assembly from the reads, and then compress against that reference.
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.
That would be true for FASTA compression but FASTQ files contain other data. With a MiSeq, for example, 1/3 of the data is the DNA sequences, 1/3 is a quality score (0-37) and 1/3 is a header which includes things like flow cell ID but also a lot of info on the physical registration address of the cluster on the flow cell. In some applications this info is not required, but is from a scientific/data integrity perspective
Who would implement this suggestion? Is it an appeal to folks at large writing tools that interact with genomics data?
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).
For those toolmakers they would probably say something like "our tool supports uncompressed data specifically so you can add whatever compression you want in your pipeline for your use cases."
It addresses different bioinformatic use cases, but plink2 is an open-source tool that is widely used and makes extensive use of zstd. In particular, it accepts both gzip- and zstd-compressed text files as input, but output commands that can generate large text files often only have a zstd compression option.
Ubuntu, Fedora, and Arch use zstd for their packages, all switching in the last few years. FreeBSD seems to have it in base. OpenBSD doesn't have a man page for it. I think Mac OS doesn't include it.
I don't know anything about DNA files, but I wonder how gzip (deflate) would fare if one of the non-default compression strategies were used.
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.
I'm not sure why gzip still pops up for FASTQ data, as it is quite easy to bin the quality scores, align it against a reference genome and compress it as e.g. CRAM [1,2].
With 8 bins, the variant calling accuraccy seems to be preserved, while drastically reducing the file size.
You don't necessarily have a reference genome to align to. For example, I've recently been working with wastewater metagenomics where (a) the sample consists of a very large number of organisms and (b) we don't have reference genomes for most of these organisms anyway.
That can be a challenge, but you can also build an "artificial" reference genome. You just use it for compression, not for any real analyses. This would allow you to still use alignment-based compression.
But I agree with you: it really depends on the type of the data.
It would be nice also that the artificial reference represented global population structure- for example, the larger the genetic distance between an individual who is sequenced, and the identity of the person who makes up the reference (an amalgam of several individuals from a common US population), the less compression you get. Instead, it seems like you could create the "genome that is the shortest distance to all other genomes" (a centroid of cluster centroids) and then the standard deviation of your compressed sizes should be much smaller.
Well I think the issue with wastewater and other screening tech is that there is no global average reference genome. In that case they're sequencing everything from phages, viruses (human and plant), bacteria, fungi, plants/animals and human...its an everything soup.
oh. From what I can tell, the total world storage for non-human genome data is trivially small (a few petabytes and not growing rapidly). Human is huge- O(petabytes)/year for a single org is not out of the question.
Thats true, but we do tremendous amounts of human DNA sequencing for certain causes at scale(e.g. understanding/treating cancer) whereas environmental sequencing is usually done to monitor/search for things at a much lower sample rate(e.g. disease load in wastewater, biodiversity from environmental samples, and looking for natural products produced by the zillions of bacteria/archaea in the oceans). From e.g. a wastewater sample perspective the latter type is going to be the majority of data, we just filter out the stuff of interest and analyze it in situ - but theres no reason to store 1B E coli genomes whereas this is necessary if we want to understand cancer evolution.
If you want to use untargeted metagenomics to detect novel human viruses you're going to be generating petabytes all by yourself: https://arxiv.org/pdf/2108.02678.pdf
I can't see any reason why you would need to save petabytes. Remember- at that scale, people think really hard about whether to pay the long-term storage and associated costs (the value of having this system should exceed its costs). The case for this already exists in (for example) cancer and other pharma.
The storage is massively cheaper than the sequencing. At some point it could be worth going back and trying to figure out how much of the raw data you can safely discarded, but at least at first there are so many more other things that are more urgent.
(The paper I linked describes more or less what I'm currently working on)
I have some questions for folks who are working in this field. In particular, are you holding onto your FASTQs for a long time (> 1 year), and if so, why? Is max compression, or ease of analysis more important (IE, do you access the data during the retention period, do you ETL it to another format, etc). How much of your budget represents storage costs which could be affected by compression?
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.
Usually fastQ filed that aren’t publicly available are retained pretty much indefinitely. The reason is, in my experience, the ability to align the sequence to new and improved genome assemblies or using/benchmarking new tools altogether. I’ve had various reasons to go back to raw sequences that had been stored for 8+ years to reanalyze them alongside newly generated data.
Is the cost worth it? Many groups are dealing with literally petabytes of fastq.gz and they usually sit around un-re-analyzed (I keep my own personal genome BAMs around and occasionally un-map them and re-map them to new references).
It really depends on the population or potential future use-case. I’ve known many labs where cost wasn’t a concern because their institution was associated with (or had their own) cluster that was under-utilized. I think there’s also something to the fact that the cost of storing sequence pales in comparison to the sample collection, DNA extraction, library prep, and sequencing itself that if you can afford to generate that much sequence then storing it isn’t much of an issue. And if you really wanted to cut costs it’s easy to upload to NCBI rather than deleting it.
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.
One useful perspective to consider is the costs of computation compared to the costs of the experiment.
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 most cases* the raw data can be recreated from the aligned BAM/CRAM file, so it's often sensible to toss the raw FASTQs after alignment. But yeah, I go back to 10+ year old genome sequences more than some might expect.
*In the absence of read trimming or discarding unmapped reads
Genome storage and compression is interesting. Most long sequences have a lot of internal redundancy/repetition that can be effectively compressed by standard algorithms, but something tuned specifically to the task, can do better. Then there is the question of storing many sequences. Whole genome sequences of thousands of bacteria from an experiment, for example. Each is almost identical to the other. This benefits from an algorithm that can compress that efficiently in terms of space, but also allows any part to be recalled random-access reasonably efficiently in terms of time. It'd also be nice to be able to add more genes to the database without having to re-compress the entire thing.
The Zarr format is used in some genomics workflows (see https://github.com/zarr-developers/community/issues/19) and supports a wide range of modern compressors (e.g. Zstd, Zlib, BZ2, LZMA, ZFPY, Blosc, as well as many filters.)
I would zarr for dense matrices, mostly (I use them with microscope images). I see it also used for frequency/spatial observations in genomic imaging.
But I prefer parquet for most direct analysis of sequence, since it's the format best integrated with big data analytics. I care much less about total compression size than I do the ability to decompress the data I need quickly (say, to ETL it to a featurization pipeline).
Many thanks. Q: what's the story around versioning, provenance, reproducibility, (etc.ops) in your domain? I've seen various bolt x on/along git variants. Wondering if its worth the effort to make something to address that.
A few years ago I built some tools https://github.com/tf318/tamtools to store alignments against two different reference assemblies in an efficient way (taking advantage of the fact that the majority of each alignment to different assemblies would in fact be the same, just shifted in position).
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.
I'm slightly surprised that speed and file size would be the only considerations. For an archival format wouldn't you at least want a format with very robust error detection, even error correction? None of the proposed formats have this, basically just having a simple 4 byte CRC or XXHASH. Some kind of random access to the compressed file might be useful too (which xz has).
I wonder if making using of zstd's long mode (`--long`) and/or multithreading support (`-T0` to use all threads) would close the displayed performance gap with `pigz`. Doesn't seem like either were used, which is odd considering the comparison made and the file used.
One thing that's missing from the article when comparing to `pigz` is that you can use the `-T0` flag in `zstd` and it will parallelize according to the number of CPUs. On some limited benchmarks, I found it to be much faster and worth using. Some `zstd` installations come with `zstdmt`[0] as an alias to `zstd -T0`.
I'm not in the bioinformatic domain but maybe there are some legacy tools that depend on some specifics of gzip.
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.
Can confirm that the data is often, though not always, stored in a shared directory and that yes, many tools can/do read gzip files directly which probably hinders the adoption of better compression algorithms.
If they do then they are just probably running zcat for you for convenience sake. You can always incorporate better compression by adding it to your pipeline, and piping in uncompressed data to the legacy tool then compressing the resulting output with whatever tool you want.
The best storage format for genomics data is DNA. At some point in the future, it might actually be cheaper to just re-sequence than to store the fastq. Instead of optimizing the digital infrastructure, it might be better to just optimize the lab infrastructure for storing the physical DNA.
Compression of FASTQ files can be greatly improved by sorting/clustering the reads. I use clumpify from BBMap for that. The bad: clumpify does not support zstd at this point.
170 comments
[ 5.2 ms ] story [ 216 ms ] threadWhen I can get away with -16 for large file long term storage, I use it.
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.
[0]: https://en.wikipedia.org/wiki/FASTQ_format
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.
Is your problem with the idea that developers are getting paid at all, or just that this instance is at Facebook?
The FOSS funding issue I more of a side-point which I didn't even bring up, but I don't see "get hired by Facebook if you want to work on FOSS" to be a healthy solution to the FOSS funding problem, we need a way to ensure funding for FOSS developers who don't work for Facebook. But that's not strictly related to the zstd thing specifically.
But FB isn't really competing with other companies here. They're giving away frameworks, libraries, tools— this is building block stuff, and the fact that it's perhaps particularly optimized for FB's workloads is just their privilege, having been the ones to fund it.
Doing this kind of thing requires a shitload of time and money.
It would be like avoiding C because it was developed by ATT.
If a good non-Facebook implementation of ZSTD shows up, I might strongly consider using (that implementation of) ZSTD.
2. The format spec is an RFC [1]. It's frozen. We're not going to pull the rug out from under you somehow (and why would we??).
3. I think we do a pretty good job accepting contributions. But again, if we aren't, you can always fork it.
All in all, we've worked pretty hard to make zstd universal--portable, stable, flexible, robust, etc. I guess I don't really have a point... I guess I hope you reconsider your rejection of zstd? It's pretty cool.
[0] https://facebook.github.io/zstd/
[1] https://www.rfc-editor.org/rfc/rfc8878.txt
I'm also concerned with the term "port". I don't know exactly what you put in the term, but to me, it sounds like someone basically translated the C code directly to other languages. What I would be looking for is more like a "clean re-implementation" by an independent group of people who properly understand all the algorithms which are used and where all the code is written from scratch without the Facebook code as a "crutch". This would ensure that whatever Facebook does, there will be a group of people that is capable of fixing tricky bugs in the algorithms or implementing further optimizations. Basically, I would want the alternative implementation to be as unrelated to Facebook's libzstd as Clang is to GCC.
All of that said though, from a purely technical perspective, zstd seems like a truly great piece of engineering. Kudos. Though I will never understand what drives great developers to dedicate their life to producing value for Facebook.
What theoretical path to harm is there to adopting it more widely?
There's lots of similar cool technologies made by other "evil" companies. I don't use ZFS, but its origins in Sun/Oracle isn't the reason I don't.
[0] https://lore.kernel.org/lkml/20181109190304.8573-1-kilobyte@...
Anyways, I prefer to not use Facebook products.
OpenZFS isn't Oracle's product; it certainly includes a great deal of Sun/Oracle code, but Oracle themselves can't use the OpenZFS commits without honoring the CDDL license, which to my knowledge they've not done.
[1]: https://fedoraproject.org/wiki/Licensing/CDDL, https://en.wikipedia.org/wiki/Common_Development_and_Distrib...
It is a BSD license, not GPL.
https://raw.githubusercontent.com/facebook/zstd/dev/LICENSE
This is wrong? Once something is GPL, it's always GPL. https://softwareengineering.stackexchange.com/a/98777/347066
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
Those who prefer zstd can use it now.[1] https://www.biorxiv.org/content/10.1101/2022.04.07.487441v1
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.
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.
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>
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.
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).
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.
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)
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
Wikipedia has a page on the topic: https://en.wikipedia.org/wiki/Compression_of_genomic_sequenc...
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.
Pretrained dictionaries could also be useful too.
[0]: https://manpages.debian.org/bullseye/zstd/zstdmt.1.en.html
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.