15
$\begingroup$

We have a lot of Illumina sequenced exome data. Currently we are using spring for its great lossless compression, but we are looking if there is anything better (and most preferably opensource) which can let us compress our fastq files.

We also want to compress or reduce file size of BAM, but using CRAM in lossless mode didn't yield good results. The space savings were decent, roughly ~10% (IMO was expecting more), but the biggest hit came when using it with IGV. We have set up a http-server and stream bams to IGV (Helped us immensely as people wouldn't have to download and visualize bam files). With Cram, the loading times were significantly high, requiring users to commit more memory for IGV. Sometimes we would have to wait around a minute for reads to show up. With Bam it is near instantaneous.

Can't use lossy compression as these BAM files might be subject to future downstream analysis (Ex. CNV calling, variant calling or something else).

What are some good ways to compress fastq and bam files so we can maximize our storage?

$\endgroup$
0

5 Answers 5

12
$\begingroup$

EDIT: I am rewriting the answer in response to updates to the original question.

TL;DR: use CRAM

Background 1: quality binning and FASTQ compression

In the old days, base callers outputted base quality at full resolution – you could see quality from Q2 to Q40 in full range. As a result, quality strings were like semi-random strings and very difficult to compress. Later people gradually realized that keeping base quality in low resolution wouldn't affect downstream analysis. The Illumina basecaller started to output quality in 8 distinct values and later changed that to 4 bins. This change greatly simplified quality string and made them compressed better. For example, in old days, a 30X BAM would take ~100 GB. With quality binning, it would only take ~60 GB.

Background 2: GATK Base quality recalibration

In early 2010s, Illumina base quality was not calibrated well. GATK people introduced BQSR to correct that and observed noticeable improvement in SNP accuracy. Nonetheless, with improved Illumina base caller, their base quality became more accurate. Meanwhile, the world moved to 30X deep sequencing. The depth overwhelms slight inaccuracy in quality. I would say around 2015, BQSR was already unnecessary for data produced at the time.

Does it hurt to apply BQSR? Yes. First, BQSR introduces subtle biases towards the reference and towards known SNPs. Second, BQSR distorts the data. At least for some datasets, I observed that SNP accuracy dropped with variant quality after BQSR; I didn't observe this with raw quality. Third, BQSR is slow. Fourth, for new sequencers producing data at higher quality, BQSR is likely to decrease data quality. Last, related to the question, BQSR added another semi-random quality string and made compression even harder. Nowadays, running BQSR is a waste of resource for worse results. The official GATK best practice no longer uses BQSR according to their WDL file.

CRAM

These days a 30X human CRAM only takes ~15 GB (see this file). This is a huge contrast to ~100 GB BAM in early 2010s. OP only saw ~10% saving probably due to a) BQSR and/or b) old data with full quality resolution.

On encoding/decoding speed, CRAM was much slower than BAM. Not any more. The latest htslib implementation of CRAM is faster than BAM on encoding and only slightly slower on decoding. The poor performance of IGV on CRAM could be that the Java CRAM decoder is not as optimized.

It is true that CRAM is not as widely supported as BAM. However, all the other alternatives are much worse. Petagene said they had IGV-PG (PDF) for their format. That is not official IGV and I couldn't find more recent update beyond the 2019 press release. I don't see other viable options.

Note that the common practice is to keep all raw reads, mapped or not, in BAM or CRAM such that you can get raw reads back later. BAM/CRAM additionally keeps metadata like read group, sample name, run information etc and is actually more popular than FASTQ in large sequencing centers. Also note that you don't need to sort CRAM by coordinate. Unsorted CRAM is only a little larger than sorted CRAM.

CRAM and its competitors

The core CRAM developer, James Bonfield, is one of the most knowledgeable researchers on compression (and one of the best C programmers) in this field. He has done a lot of compression evaluation over the years. The conclusion is that on a fair benchmark, CRAM is comparable to the best tools so far in terms of compression ratio.

Petagene could compress better in the plot @terdon showed mostly because it has a special treatment of the OQ tag generated by GATK BQSR. It is a typical trick marketing people use to make their methods look better. With BQSR phased out, this plot is no longer relevant.

On commercial software

In general, I welcome commercial tools and think they are invaluable to users. I also have huge respect to DRAGEN developers. However, on FASTQ storage, I would strongly recommend against closed-source compressors. If those tools go under, you may lose your data. Not worth it.

$\endgroup$
9
  • $\begingroup$ Thanks for your input, but mainly for the link. WDL seemed too complicated for me. We also noticed that BQSR was not recommended by nf-core people as well, but none provided a concrete reason beyond "its not useful". As I've also commented on @gringer post, CRAM didn't work well when using it with IGV, as well as some downstream tools don't work with CRAM. $\endgroup$ Apr 12, 2023 at 9:55
  • $\begingroup$ @KarthikNair It is interesting that you are criticizing CRAM for slowness in IGV, while you are compressing fastq with spring and you can't even load other alignment formats into IGV. What are you expecting? Also, you don't need to understand WDL. Just search/grep "BaseRecalibrator" or "BQSR" to see if it gets called. $\endgroup$
    – user172818
    Apr 12, 2023 at 15:07
  • $\begingroup$ @KarthikNair sadly, the GATK best practices, which were updated very recently, still recommend BQSR so it is really surprising that their recommended pipeline doesn't use it. That said, in my place of work, we analyze around 2k human samples a month, and we have found that BQSR often causes us to miss variants. I don't have statistics for this yet so it's anecdotal, but there is strong evidence that BQSR not only isn't particularly helpful today, but can be actively harmful. And not only with small panels, this also happens with whole exome sequencing samples. $\endgroup$
    – terdon
    Apr 12, 2023 at 15:34
  • $\begingroup$ @terdon That page is referencing their older pipelines, which do use BQSR. The WDL file in my link is two years old. It is the latest I can find. IMO, BQSR has been hurting for 10 years. Well, better late than never. $\endgroup$
    – user172818
    Apr 12, 2023 at 15:49
  • $\begingroup$ you focus on BAM but what would you suggest for FASTQ if not DRAGEN ORA reference based approach? Arguably long term of storage of FASTQ is more important since BAMs can be regenerated $\endgroup$ Apr 12, 2023 at 17:02
5
$\begingroup$

I'd recommend reconsidering CRAM. It's an open standard that is most likely to still be accessible many years in the future. Regarding your point about lossy compression...

Can't use lossy [CRAM] compression as these BAM files might be subject to future downstream analysis (Ex. CNV calling, variant calling or something else).

By default, CRAM compression is not lossy.

CRAM does, however, have an optional archive settings mode (samtools view --output-fmt cram --output-fmt-option archive ...) which is a lossy compression, doing things like removing read names, removing additional accessory fields, and additional compression of quality scores. In all cases, the base sequence of the reads is preserved.

https://www.htslib.org/doc/samtools.html

$\endgroup$
0
3
$\begingroup$

An alternative to the good option of PetaGene that @terdon mentioned and also not open source is the DRAGEN ORA compression. This is loseless reference based compression and can be ~5 x smaller over fastq.gz, but it's only for fastq files and only those generated by Illumina machines (and will be easiest within the DRAGEN ecosystem).

Another point to consider is where you store the data. Even with 'standard' compression you might get more benefits moving archival data to a system like AWS S3 Glacier, depending on your current infastructure setup of course

$\endgroup$
2
  • 1
    $\begingroup$ This is exactly what I do. I archive FASTQ.GZ (with added content-md5 checksum) to S3. I also archive processed outputs to a different location that is less fault tolerant than S3 as those files can be recreated. $\endgroup$
    – Ram RS
    Apr 11, 2023 at 15:18
  • $\begingroup$ We did check this out, but we do not sequence the NGS data ourselves. We ask external labs to perform the sequencing and they share the fastq.gz files with us. What is the approximate pricing for this? I can't clearly identify where those details are! $\endgroup$ Apr 12, 2023 at 9:41
3
$\begingroup$

The only free and open source tool I know that can help is zstd. Their github repository's README describes it as:

Zstandard, or zstd as short version, is a fast lossless compression algorithm, targeting real-time compression scenarios at zlib-level and better compression ratios. It's backed by a very fast entropy stage, provided by Huff0 and FSE library.

We read this blog post on using it to compress fastq files and were intrigued so we ran some tests. We couldn't reproduce the same level of compression that the article claims but that might be due to different input files or the way an index was/wasn't created.

Even in our tests though, we did see that it provides better compression than gzip. Here are the specific values for one of the tests, run using an R1 file of a human WGS run:

Uncompressed fastq bgzip-compressed zstd-compressed
179G 58G 47G

This was done using -z 19, to set the compression level to 19 (we found this to be the best tradeoff between size and time to compress). We also saw that using a human genome fasta file to create a compression index, and then using that when compressing human fastq files, actually resulted in worse performance (larger compressed file) than not using the index. But it is entirely possible that we just did it wrong, we didn't spend a lot of time on it.

So you might want to investigate zstd as well, maybe you can get better performance if you play around with the indexing. On the plus side, it is free and open source. On the other hand, it won't do much with bam or cram files.

$\endgroup$
1
  • $\begingroup$ Thanks for the input. I'm pretty sure that with Spring, the output size might be around 20-30 GB, if not slightly lesser. Granted, you would have to compress / decompress the data, but the final file size would be pretty small. We have compressed WGS samples with spring, it goes from total 50GB of gz files to around 15-20 GB in final spring file. $\endgroup$ Apr 12, 2023 at 9:37
3
$\begingroup$

You might want to consider PetaGene. This is a paid tool, and not open source, but it is something we use at work and it is really good! It can do lossless compression of gzip-compressed fastq files, bam files, and even cram files, and achieves really good size reduction (image from the PetaGene webpage):

petagene compression performance bar graph

In my own testing, using files of different sizes, I get an average size reduction of 62.1% (2.6 x smaller) for fastq.gz files (compared to the gzipped fastq file, not compared to the original, uncompressed fastq file) and an average reduction of 82.1% (5.6 x smaller) for bam files. I didn't test cram files.

If you have the budget for a paid service, these guys do good work! The greatest advantage, apart from the very significant size reduction, is that they also provide a system that allows you to use the compressed files natively. So you can feed the petagene-compressed bam file directly to bwa for example. They create special symlinks (foo.fastq.gz -> foo.fastq.petacompressed) and you can then use these symlinks and petagene handles the decompression on the fly.


The Wikipedia page on CRAM has a nice table with a few alternatives to CRAM which I reproduce below. I haven't tried any except PetaGene though so cannot comment on their performance:

Name Key advantage over CRAM Repo
Deez None. An early approach to SAM compression, introduced 2014. https://github.com/sfu-compbio/deez
Genozip Significantly better compression, but slower random access https://github.com/divonlan/genozip
GenomSys Utilizes ISO-standard MPEG-G, but compression inferior to CRAM 3.1 Source code not publicly available
PetaSuite Virtual filesystem interface Source code not publicly available
$\endgroup$
4
  • $\begingroup$ We considered this software, but sadly we don't have the budget for this. And yes, on-the-fly decompression is a great feature to have, makes storing all data in this format extremely viable. Right now, we get decent compression with spring, but it needs lots of de/compression time, not to mention the CPU resources it guzzles during the process lol $\endgroup$ Apr 12, 2023 at 9:43
  • $\begingroup$ @terdon genozip is quality work and is very impressive. Its full source code is also on github. If genozip goes under some day (I hope not!), you may create a fork and continue to use their format. $\endgroup$
    – user172818
    Apr 12, 2023 at 15:12
  • $\begingroup$ @user172818 ~terdon - can you share any further info on how PetaGene achieves compression and the critique mentioned of "special treatment of the OQ tag generated by GATK BQSR"? I know it's closed source but trying to understand conceptually how their compression works $\endgroup$ Sep 29, 2023 at 8:17
  • $\begingroup$ Not I, @Chris_Rands, sorry. I'm sure user172818 will know far more, but I have no idea and wasn't even aware there was any special treatment before reading their answer. I know at least one person who used to be an engineer working for them is active on the site. If there's something that can be shared, maybe they will answer if you post a question. $\endgroup$
    – terdon
    Sep 29, 2023 at 8:22

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.