At BOSC 2012, we heard a
report from Richard Holland on the Pistoia Alliance Sequence Squeeze
competition. I'd run across this
a couple of times before -- most notably in the Quip paper
-- and was interested in hearing the results.
What was the problem being tackled? To quote,
The volume of next-generation sequencing data is a big
problem. Data volumes are growing rapidly as sequencing technology
improves. Individual runs are providing many more reads than before,
and decreasing run times mean that more data can today be generated by
a single machine in one day than a single machine could have produced
in the whole of 2005.
Storing millions of reads and their quality
scores uncompressed is impractical, yet current compression
technologies are becoming inadequate. There is a need for a new and
novel method of compressing sequence reads and their quality scores in
a way that preserves 100% of the information whilst achieving
much-improved linear (or, even better, non-linear) compression ratios.
In the end, though, I found the results
disappointing. The main conclusions?
"First, that the gzip algorithm that served as the competition
baseline is actually quite sufficient for run-of-the-mill
compression, and second, that it's extraordinarily difficult to
make huge improvements in all three of the judged dimensions," said
Nick Lynch, external liaison of the Pistoia Alliance and chair of
the judging panel.
In some sense, this is pretty reassuring -- if you're compressing
data, it's nice that gzip
is hard to beat! Even in situations like this, where you have
domain-specific data, gzip can compete reasonably
well. Better compression algorithms do exist, but there's a severe
In another sense, though, it's insane. We should be able to
compress the hell out of most short-read sequencing data! Short-read
sequencing data is massively redundant for two reasons. First,
because to achieve sufficient coverage with shotgun sequencing you need to
sample deeply on average. For example, if you have a 100x coverage
genomic sample, then more than half of your content is represented
more than 100 times in your data set. Second, because most of the
time what you're sequencing has already been sequenced before -- if
you're looking for genomic variants (in resequencing data), you expect
most of the data to be identical to what you already have. That is,
if you are trying to figure out how my genome is different from, say,
Craig Venter's, you can rely on the fact that, as two white guys,
we're probably about 99.99% identical at the genomic level. Even if
you're sequencing an Australian aboriginal
you can guess that our genomes would be 99.9% identical (this number
is made up -- corrections welcome. Is it 99.8 instead? ;) More
generally, the vast majority of sequencing is biomedical (because
that's where the money is) and with only a few exceptions, the vast
majority of that data is going to be from a few well-sequenced
species, and (from an information content perspective) essentially the
same as what everyone else is sequencing. Even if you want to keep
all the single-nucleotide variants and predicted insertion-deletion
variants, we're talking about a few megabytes per data set --
Some short-read data is more interesting - there are some very
interesting applications of sequencing, mostly in environmental,
ecological, and evolutionary studies, for which the theoretical
redundancy of the data sets is lower than that for straight genomic or
biomedical sequencing. For example, I would venture to say that our
current soil metagenomic data sets are only about 80-90% compressible;
this is compared to the common genomic and RNAseq data sets, which
should be (info-theoretically) between 99-99.9% compressible.
(Footnote: this is my mental calculation for the ratio between an
appropriately 2-bits-per-base uncompressed storage, and the maximal
information theoretic lossless compression. There's some fancier term
for this that I'm forgetting...)
But the main point is this -- most short-read sequencing data is crazy
redundant, WAY past the few 30x maximum compression ratio achieved by
even the best of the sequence squeeze participants.
So, why didn't the sequence squeeze competitors do better?
The answer is simple: FASTQ. Specifically, quality scores.
A digression: how Illumina sequencing works, kinda
So, you have some DNA. And you want to sequence it.
One very common way is to use an Illumina machine, like so:
- You fragment the DNA into molecules of length ~300-600 bases.
- You stick the DNA randomly to a 2-D surface.
- You amplify the DNA where it sits, so now you have lots of mostly identical
fragments of DNA in close proximity to each other.
- You strip the DNA down to single-stranded fragments.
- Then you repeatedly add individually fluorescently colored A, C, G,
and T bases to every single-stranded DNA molecule on the surface,
thus reading off the next base of sequence on the opposite strand.
- As you do #5, you image the 2-D surface with a microscope.
- You take the resulting movie and simply "read" the DNA sequence present
at each spot on the 2-D surface off from the twinkling little lights
at each spot.
- These DNA sequences are transcoded into strings of DNA, together
with a quality score that correlates with how confident the reading
algorithm in #6 is about the identify of each particular character of
DNA. This is all put into a FASTQ file, which is one "output" of the
You can view this as an incredibly complex but ridiculously effective
way to digitize DNA into bits, and it's how most sequencing is done
these days. (There are many other sequencing schemes, and no doubt
this will change. But as of October 2012, Illumina is by far the most
used sequencing approach in terms of volumes of data.)
Back to compression
The problem that the Sequence Squeeze participants ran into is,
basically, that the terms of the context specified lossless
compression -- whatever process you used had to retain all the bytes.
This is really troublesome because of one main thing -- the noise
inherent in the data.
First, the FASTQ format specifies a range for
quality scores of between 60-90 in length, and the information that
these quality scores encode is noisy. Noise is hard to encode
efficiently, and the more resolution you bring to bear on the noise,
the less compression you can apply. While the ultimate compression
you can do on the 60-character quality-score alphabet is going to be
very dependent on the frequency distribution of those 60 characters,
no matter how you cut it a 60-character alphabet is much less
encodable than a 4-character alphabet like DNA, especially if
there's randomness in the 60-character alphabet.
Second, the DNA sequence itself is noisy, too. The error rate of
Illumina is in the 0.1% to 1% range, which automatically limits the
compressibility of the data. This noise essentially limits the degree
to which you can use prior knowledge (like reference genomes) to
compress things, because the noise is random (i.e. uncorrelated with
the reference genome). But these limitations are relatively minor
compared to the noise in the FASTQ format.
There are (at least) two solutions to these problems. The first is to
reduce the alphabet size of the quality scores. The second (more
general) solution is to allow lossy compression.
The first is obvious, and is itself a kind of lossy compression, in the
sense that along with discarding the noise in the quality scores,
you're also going to discard some true information. But I think the
second is much more interesting.
Lossy vs lossless compression
The Wikipedia page on lossy compression is excellent.
The term essentially applies to discarding information in the
service of compression, where "information" is defined information
theoretically to mean pretty anything non-random. Lossless
compression means you can recover exactly the set of bits you had
to begin with; lossy compression means that you lose some of
the data. The trick with lossy compression is to pick which of the
data you discard in order to minimize the effect on downstream usage.
It turns out we all use lossy compression on a daily basis -- the
image and video encoding used for JPEGs and MPEGs are lossy! Most if
not all lossy compression schemes allow you to specify how much
information you lose, and a cool graphic demonstration of that is
available here, for JPEG images. You
can compress photographs a surprisingly large amount before they
With images, the JPEG compression algorithm is designed to keep the
image as recognizable as possible. Or, to put it another way, JPEG is
attempting to discard as much data as possible while retaining as much
mutual information as it can between the image and the downstream
application -- in this case, your visual cortex and the image
recognition capabilities of your brain. If you can throw away 95% of
the data and still have your lolcatz, then why not? Again, the key
point for lossy compression is downstream application -- what
information are you trying to keep, vs which information are you
trying to discard?
For sequencing applications, there's actually a pretty solid answer to
"what do we want to do with the data?" Scientists generally have one
of a few immediate downstream applications in mind -- sequencing
things for assembly, resequencing things to find novel variants, and
resequencing things to count them -- and it is relatively
straightforward to envision what lossy compression would be in these
cases. Assembly is an excellent example of lossy compression, in
the sense that you are trying to maximize the mutual information
between your sequencing data and the predicted genome output by the
assembler. For resequencing, the pattern of variants from some
reference is, basically, a compressed version of the sequencing data.
And, for counting applications, you can summarize counts on a
per-molecule or per-region basis.
Put another way, sequencing error is the information you don't
care about. You really only care about the bases that might be
real. With the high error rate of sequencing, the vast majority of
"information" in your data set is actually erroneous with respect to
the source biological data. And this puts the Pistoia Sequence
Squeeze competition into sharp relief: by insisting that software
perform lossless compression, they hamstrung themselves by effectively
requiring that all the error be kept around.
What do I think we should be doing?
A pipeline view
If you view sequence analysis as a pipeline -- produce data, subject
data to some analysis that leads to a prediction, and use the
predictions in further analyses -- then you're probably either
confused or concerned right now. Sequencing and bioinformatics is a
rapidly evolving field, and new sequencing technologies and analysis
packages are being developed up all the time; why would we want to
compress our data into a result by saving the output of, say, the
GATK toolkit? That would
prevent us from taking advantage of improvements to GATK, or
competitors to GATK! This thinking leads to the conclusion that I bet
motivated the Sequence Squeeze competition: we need to save all the
raw data, in case we come up with better ways to analyze it!
Fig 1. A somewhat abstract view of an analysis pipeline.
I'm sympathetic to the idea that we don't want to save all the output
from any given analysis package. All too often we are feeding our
data into a gigantic black box of an algorithm/software package/blob
and we have very little idea of how it will work specifically on our
data. Even when we do, we know that it's going to be making some
choices with respect to sensitivity and specificity, and perhaps
those will be the wrong choices.
However, given the contents of the FASTQ format, it is also pretty
crazy to insist that all of that stuff is worth keeping. Complex
proprietary algorithms are used to call bases and assign quality
scores, and we can be reasonably certain that there are some serious
heuristics (i.e. guessing) being applied; choosing to keep the exact
number spit out by the sequencer is pretty arbitrary. There is little
reason to believe that a 30-character quality score alphabet would be
significantly worse than a 60-character quality score alphabet. Heck,
in most cases, a variant call or a base call can be made with greater
than 99.99% certainty; I bet you could quantize the quality score
encoding down to "rock solid" and 7 shades of "less solid" without
significantly affecting your base calls.
On the scale of things, I don't think we should make any one part of
the pipeline sacred. We should be figuring out how to optimize the
data we keep in order to meet specific limits on analysis sensitivity;
and, in particular, we should be thinking about how to maximize the
re-analyzability of the data in light of the specific sensitivity we
need for each downstream application.
For example, if we are sequencing a single diploid individual, we can
pretty confidently throw away minor variant information; if you have
100-1 allele calling at a given location your answer is pretty much
guaranteed correct. If you're sequencing a population sample, e.g.
a tumor, you might have more trouble picking a high
threshold; heck, you might even want to keep all that data. But the
choice should be made for each application.
The interesting thing is that such an "arbitrary" line in the sand has
already been drawn! We don't keep the raw image files from the
sequencing movies, which actually contain the maximum amount of
digital information; we don't keep them because they're freakin' huge
(multiple terabytes per run), and it doesn't make economic sense.
It's easier and cheaper to resequence the 0.1% of samples that you
need to reanalyze from scratch than it is to keep around the 99.9% of
unnecessary data. It's just silly to argue that all the data is
sacred, and choosing the FASTQ files as your "perfect" stage is
not particularly meaningful.
So what should we keep? I bet we can compress the data down to a
certain point, discarding things that are obvious errors, while still
retaining all the signal and a certain amount of noise. This would
still make it possible to do any downstream analysis you wanted --
just on a vastly decreased set of data.
Our digital normalization approach approximates this for assembly: it
discards quite a bit of data and noise while keeping 99.999% of the
information. However, diginorm discards the abundance information.
Error correction, especially error correction that is tunable to
retain polymorphism, would be a more sensible long term approach: it
would enable lots of compression while still retaining abundance.
What we're working on in my lab is technology to do this efficiently
One of the key perspectives that I think our work (partitioning and
digital normalization) brings to the table is the idea of
prefiltering and eliminating known-bad data, rather than trying to
guess at what is going to be good. We are trying to keep all of the
information, in as generic a sense as possible, while discarding
noise and scaling the problem; whether or not that information is
actually used by any particular downstream analysis software is
irrelevant. Another way to look at it is that we're trying to discard
obvious noise and errors, rather than make specific calls about what
is going to be interesting -- see Wikipedia on the LHC Computing
Grid for an
example of this in particle physics. This perspective is good because
it gives us the ability to choose what kinds of tradeoffs we want
without being unduly tied to any one downstream analysis platform.
Fig 2. Inserting lossy compression into an analysis pipeline.
It's important to note that this is not a new idea by any means -- I
vaguely remember reading about cram
a while back. This is a package that Ewan Birney et al. are working
on that does reference based compression -- exactly what I'm talking
about above wrt human DNA. (Thankfully Nick Loman pointed it out to
me when I asked him to review a draft of this post and so I can refer
to it properly ;). Cram, like Quip, attempts to make use of assembly
techniques to generate a new reference for unmappable reads, while
taking the mapped reads back to the reference for reference-based
compression. Perhaps most interestingly to me, they are thinking
about quantizing quality scores for exactly the reasons I point out
above. You can read the full paper here.
What should the next competition be?
I'd love to see an open competition like the Pistoia Sequence Squeeze
contest, but focused on lossy compression, as a prefilter to one or
more specific downstream packages. That is, provide a way to
calculate sensitivity and specificity from the output of (say) GATK,
and also measure the time and space complexity of the entire pipeline
from raw data through to end analysis. You could of course do this
all on the cloud, to ensure a reasonably uniform and widely available
On the one hand, this would encourage people to shave data down to the
bare minimum such that GATK with the given parameters "got it right"
-- not necessarily a good outcome. But, on the other hand, it would
show what the true limits were, and it would probably drive a whole
bunch of research into what the various tradeoffs were and how to
design such algorithms, choose parameters, and evaluate their
Triage, compression, and Big Data
As one BOSC attendee told me, with the size of the data we're starting
to generate, we have to figure out how to compress the data, scale or
parallelize the algorithms, or discard some of it. Some data will be
gold, some will be copper, and some will be dross; it should at least
possible to efficiently separate much of the dross from the
My real interest is in enabling the appropriate economic decisions to
be made. Storage and processing time are considerations that should
to be discussed and explored for data-driven biology to become
less handicapped, and so far we haven't done much of that exploration.
I think this could end up having serious and dramatic consequences
for clinical data processing down the road, in addition to obvious
One caveat is that there are also legal issues surrounding data
retention: for example, clinical folk are apparently mandated to keep
all their "raw data" (which means FASTQ, not images, apparently). I'm
not prepared to take on the social and legal aspects of this proposal
:). But even so, you can dump your big data to very cold storage and
retain the filtered data for daily or monthly or yearly reanalysis.
I don't have any, really; I'll post updates from our work as they
occur. I'd love to hear about other people's attempts to do this
(it's not a new idea! see cram, above) and also about why you think
it's a horrible idea for NGS analysis...
Special thanks to Nick Loman, Jared Simpson, and Zamin Iqbal for sanity
checking and helpful references!
p.s. After writing all of this, I discovered Ewan Birney's May post
on the same topic
in which he basically reprises everything I said above. So, um, YEAH.
The question is, can our lab put all of that on a streaming basis?
Answer: yes :).
p.p.s. Jared Simpson also pointed out this very interesting paper on
using the Burrows-Wheeler Transform to compress large genomic databases:
Cox et al., Bioinformatics, 2012.
There are comments.