Motivation: The bulk of space taken up by NGS sequencing CRAM files consists of per-base quality values. Most of these are unnecessary for variant calling, offering an opportunity for space saving. Results: On the Syndip test set, a 17 fold reduction in the quality storage portion of a CRAM file can be achieved while maintaining variant calling accuracy. The size reduction of an entire CRAM file varied from 2.2 to 7.4 fold, depending on the non-quality content of the original file (see Supplementary Material S6 for details). Availability and implementation: Crumble is OpenSource and can be obtained from https://github.com/jkbonfield/crumble. Supplementary information: Supplementary data are available at Bioinformatics online.This work was funded by the Wellcome Trust [WT098051]

Bonfield, James K

Durbin, Richard

McCarthy, Shane A

Apollo (Cambridge)

Sequence analysisCrumble: reference free lossy compression ofsequence quality valuesJames K. Bonfield1,*, Shane A. McCarthy1,2 and Richard Durbin1,21DNA Pipelines, Wellcome Sanger Institute, Wellcome Genome Campus, Hinxton CB10 1SA, UK and 2Departmentof Genetics, University of Cambridge, Cambridge CB2 3EH, UK*To whom correspondence should be addressed.Associate Editor: Inanc BirolReceived on January 5, 2018; revised on June 18, 2018; editorial decision on July 4, 2018; accepted on July 9, 2018AbstractMotivation: The bulk of space taken up by NGS sequencing CRAM files consists of per-base qualityvalues. Most of these are unnecessary for variant calling, offering an opportunity for space saving.Results: On the Syndip test set, a 17 fold reduction in the quality storage portion of a CRAM file canbe achieved while maintaining variant calling accuracy. The size reduction of an entire CRAM filevaried from 2.2 to 7.4 fold, depending on the non-quality content of the original file (seeSupplementary Material S6 for details).Availability and implementation: Crumble is OpenSource and can be obtained from https://github.com/jkbonfield/crumble.Contact: jkb@sanger.ac.ukSupplementary information: Supplementary data are available at Bioinformatics online.1 IntroductionThe rapid reduction of costs for genome sequencing (Wetterstrand,2016) has led to a corresponding growth in storage costs, far out-stripping Moore’s Law for CPU and Kryder’s Law for storage. Thishas led to considerable research into DNA sequence data compres-sion (Numanagic et al., 2016).The most significant component in data storage cost is the per-nucleotide confidence values, which carry information about thelikelihood of each base call being in error. The original CRAM pro-posal (Fritz et al., 2011) introduced the term ‘quality budget’ forlossy compression. Given a fixed amount of storage we can decidehow to spend this budget, either by uniform degradation of all qual-ities or more targeted fidelity in important regions only. How to tar-get this has been the focus of lossy compression research, with twomain strategies: ‘horizontal’ and ‘vertical’.‘Horizontal’ compression smooths qualities along each sequencein turn, as implemented in libCSAM (Ca´novas et al., 2014), QVZ(Malysa et al., 2015) and FaStore (Roguski et al., 2018) or viaquantization (Illumina, 2014). This type of compression can beapplied before alignment and is entirely reference free.‘Vertical’ compression takes a slice through an aligned dataset inthe SAM format (Li et al., 2009) to determine which qualities tokeep and which to discard, as used in CALQ (Voges et al., 2018), orvia hashing techniques on unaligned data in Leon (Benoit et al.,2015) and GeneCodeq (Greenfield et al., 2016). Traditional lossmeasures, such as mean squared error, will appear very high, butthese tools focus on minimizing the changes in post-processed data(variant calling).We present Crumble as a mixture of both horizontal andvertical compression. It operates on coordinate sorted alignedSequence Alignment/Map (SAM), Binary Alignment/Map (BAM)or CRAM files. Although this approach does not explicitly use areference, the sequence aligner does, which may result in some ref-erence bias.2 Materials and methodsA variant caller evaluates the sequence base calls overlapping eachgenome locus along with their associated qualities to determinewhether that site represents a variant. Irrespective of whether thecall is a variant, if the same call is made with comparable confi-dence both with and without sequence quality values present then itcan be concluded that the qualities are not necessary in thatcolumn.VC The Author(s) 2018. Published by Oxford University Press. 337This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0/), which permitsunrestricted reuse, distribution, and reproduction in any medium, provided the original work is properly cited.Bioinformatics, 35(2), 2019, 337–339doi: 10.1093/bioinformatics/bty608Advance Access Publication Date: 10 July 2018Applications NoteThis requires running the variant caller twice to assess the change,but if limited to sites with high confidence calls the need for a secondtest can be avoided. We implemented a fast, but naı¨ve, caller derivedfrom Gap5’s consensus algorithm (Bonfield and Whitwham, 2010).This is a pure pileup-column oriented approach that treats the lack ofa base (deletion) as a fifth base type (‘*’) and then identifies the mostlikely homozygous or heterozygous combination of bases that matchthe observed base calls, confidence values and mapping qualities.Thus it assumes a single individual diploid sample and is not tuned towork with somatic variants. This is further modified by reducing theconfidence for the consensus by the bases which do not match the hy-pothesis, thus producing a deliberately pessimistic caller. The aim isnot to have a built-in high-quality caller, but to preserve quality valuesif any downstream variant caller may be uncertain while retaining in-dependence from any standard tool.Even when deemed unnecessary, qualities cannot be entirely dis-carded as tools expect them to exist. By replacing the qualities forbases that agree with a confident consensus call with constant highvalues, the entropy of the quality signal is reduced. Quality valuesfor bases that disagree with a confident consensus call may optional-ly be set to a constant low value, heavily quantized or left intact.There are sites where any variant caller may incorrectly give thewrong call with high confidence. Furthermore the reference itselfmay be incorrect and a subsequent realignment to an updated refer-ence may change read locations and alignment strings. We do notwish to replace qualities in such regions. We therefore have a set ofheuristics to try to find potentially unreliable calls and retain verba-tim the confidence values for these locations and surrounding basesdepending on sequence context. Similarly there may be places wherean entire read needs to have qualities retained as there is evidencefor it being misplaced or being part of a large structuralrearrangement.The heuristics used in Crumble to identify where confidence val-ues should be retained vary by compression level requested, butinclude:• Concordant soft clipping: many reads having soft clipped basesat the same site often indicates a large insertion (absent in the ref-erence) or contamination.• Excessive depth: possible contamination or collapsed repeat.Variant calls often appear unusually good in such data, evenwhen wrong.• Low mapping quality: possibly caused through poor reference.We optionally can also store quality values for the reads withhigh mapping quality that colocate with many low mappingquality reads.• Unexpected number of variants: we assume data from a singlediploid sample with at most two alleles at each locus. More thantwo alleles imply misaligned data, duplication or contamination.• Low quality variant calls: typically a single base locus where theconsensus is unclear.• Proximity to short tandem repeats: alignments are often poor insuch regions, especially if indels are present, leading to basesoccurring in the wrong pileup column.Finally for the quality values that we deem necessary to keep, weoptionally provide horizontal compression via the P-block algorithmfrom CSAM. This is most useful on older Illumina datasets withover 40 distinct levels of quality values.The nature of the Crumble algorithm makes it amenable tostreaming and it does not require large amounts of memory tooperate.3 ResultsAnalysis of how quality compression affects variant calling wasperformed on Syndip (Li et al., 2018), an Illumina sequenced li-brary artificially constructed from the haploid cell lines CHM1 andCHM13, with an associated high quality truth set based on twoPacBio assemblies (Schneider et al., 2017). When compared withthe Genome in a Bottle (GIAB) or Platinum Genomes (PlatGen)datasets this has a considerably larger set of tricky indels in thetruth set, giving SNP false positive rates 5–10 times higher (Li et al.,2018) than on GIAB or PlatGen truth sets. Although Syndip stillrequires a list of regions to exclude, the total number of excludednon-N reference bases is 40% fewer than GIAB 3.3.2. By restrictinganalysis to solely the regions within Syndip and not within GIABwe observe 65% of Chromosome 1 false positives occur within thisregion, but crumble still shows good performance (seeSupplementary Material).The input BAM file (ERR1341796) had previously been createdwith GATK best practices including IndelRealigner and BaseQuality Score Recalibration steps. To test the impact on raw variantcalling, we ran GATK HaplotypeCaller (Poplin et al., 2017),Bcftools (Li, 2011) and Freebayes (Garrison and Marth, 2012), fil-tering to calls of quality 30 or above, without use of GATK VariantQuality Score Recalibration. As a baseline we compare Crumble tothe original lossless results and against a single fixed quality value.This latter test demonstrates that quality values are important, butwe only need a small quality budget to achieve comparable resultsto lossless compression. Indeed, we observe that vertical qualityscore compression can marginally improve variant calling by stand-ard callers, as has been noted previously in the QVZ (Malysa et al.,2015) and Leon (Benoit et al., 2015) papers.Table 1 shows the GATK lossless results on the Syndip along withthe changes caused by lossy compression using a variety of Crumbleoptions on both the full Syndip data and a low coverage subset. Wechose the minimal compression level, an expected maximum compres-sion level and a set of manually tuned parameters optimised for thisdataset. The manual tuning traded false positives and false negatives inan attempt to get a call set comparable or better than the original in allregards. It is unknown if the tuned parameters are appropriate for alldatasets. More complete comparisons including against other tools areavailable in the online Supplementary Material.On the original BAM file with 50 coverage we observed a 17fold reduction in the size of CRAM compressed quality values, whileachieving a 6% drop in filtered SNP false positive rate (higher preci-sion) and 2% drop in false negative rates (higher recall). Indels alsosee a 1% improvement in both measures. At a sub-sampled 15coverage we see a 1% drop in filtered SNP false positive rates and a10% reduction in SNP false negatives. Indel calls were more com-parable, with 1% higher false positives and 3% lower falsenegatives.It is likely these gains to both SNP precision and recall only applyto data coming from a single individual, but they demonstrate the ef-ficacy of lossy quality compression.4 ConclusionWe have demonstrated that Crumble, when combined with CRAM,can greatly reduce file storage costs while having a minimal, if notbeneficial, impact on variant calling accuracy of individual samples.For maximum compression Crumble also permits discarding readidentifiers and some auxiliary tags, typically yielding files in the size338 J.K.Bonfield et al.of 5–10 Gb for a 30 whole genome processed with Crumble -9p8.Using this across a variety of BAM and CRAM files Crumble gavean overall file size reduction from 3- to 7.8-fold (details inSupplementary Material).Crumble is designed to operate on a single sample file. For mul-tiple samples, it is best to apply Crumble to each sample independ-ently, produce gVCF, and then jointly call from those. NoteCrumble is explicitly designed to operate on diploid data, so it is notappropriate for use on sequence from cancer or other samples withsubclonal genetic structure.AcknowledgementsWe would like to thank Yasin Memari for help testing and evaluating an ear-lier version of the program.FundingThis work was funded by the Wellcome Trust [WT098051].Conflict of Interest: none declared.ReferencesBenoit,G. et al. (2015) Reference-free compression of high throughputsequencing data with a probabilistic de Bruijn graph. BMC Bioinformatics,16, 288.Bonfield,J.K. and Whitwham,A. (2010) Gap5–editing the billion fragment se-quence assembly. Bioinformatics, 26, 1699–1703.Ca´novas,R. et al. (2014) Lossy compression of quality scores in genomic data.Bioinformatics, 30, 2130–2136.Fritz,M.H.-Y. et al. (2011) Efficient storage of high throughput dna sequenc-ing data using reference-based compression. Genome Res., 21, 734–740.Garrison,E. and Marth,G. (2012) Haplotype-based variant detection fromshort-read sequencing. arXiv Preprint arXiv, 1207, 3907.Greenfield,D.L. et al. (2016) GeneCodeq: quality score compression andimproved genotyping using a Bayesian framework. Bioinformatics, 32,3124–3132.Illumina (2014). Reducing whole-genome data storage footprint. Technicalreport. www.illumina.com/documents/products/whitepapers/whitepaper_datacompression.pdf (6 October 2017, date last accessed).Li,H. (2011) A statistical framework for snp calling, mutation discovery, asso-ciation mapping and population genetical parameter estimation fromsequencing data. Bioinformatics, 27, 2987–2993.Li,H. et al. (2009) The sequence alignment/map format and SAMtools.Bioinformatics, 25, 2078–2079.Li,H. et al. (2018) A synthetic-diploid benchmark for accurate variant callingevaluation. Nat. Methods. doi: 10.1038/s41592-018-0054-7.Malysa,G. et al. (2015) QVZ: lossy compression of quality values.Bioinformatics, 31, 3122–3129.Numanagic,I. et al. (2016) Comparison of high-throughput sequencing datacompression tools. Nat. Methods, 13, 1005–1008.Poplin,R. et al. (2017) Scaling accurate genetic variant discovery to tens ofthousands of samples. bioRxiv, doi: 10.1101/201178.Roguski,Ł. et al. (2018) FaStore: a space-saving solution for raw sequencingdata. Bioinformatics, 34, 2748–2756.Schneider,V.A. et al. (2017) Evaluation of GRCh38 and de novo haploid gen-ome assemblies demonstrates the enduring quality of the reference assembly.Genome Res., 27, 849–864.Voges,J. et al. (2018) Calq: compression of quality values of aligned sequenc-ing data. Bioinformatics, 34, 1650–1658.Wetterstrand,K.A. (2016). DNA sequencing costs: data from the NHGRI gen-ome sequencing program (GSP). www.genome.gov/sequencingcostsdata.(6 October 2017, date last accessed).Table 1. Effect of lossy quality compression on 50 and 15 Syndip data using GATK HaplotypeCallerCategory Original Original F Crumble-1 Crumble-1 F Crumble-9p8 Crumble-9p8 F Crumble* Crumble* F50 Qual size (MB) 4107 — 614 — 235 — 229 —50 SNP False Positive 6226 2968 –359 –79 –251 –67 –526 –18150 SNP False Negative 4648 7625 0 –53 –25 –184 þ41 –12350 Indel False Positive 3965 3649 –7 –41 þ19 þ9 –35 –3250 Indel False Negative 7881 7972 þ7 þ11 –103 –82 –93 –7215 Qual size (MB) 1211 – 260 — 77 — 72 —15 SNP False Positive 4798 2517 –10 þ63 þ347 þ225 –359 –2915 SNP False Negative 14985 27761 –205 –297 –3027 –4608 –1866 –286515 Indel False Positive 2781 2521 þ2 –14 þ109 þ60 þ53 þ2615 Indel False Negative 13136 13925 –8 þ5 –484 –427 –444 –410Note: Comparison of unfiltered and filtered (marked with ‘F’) calls on the Syndip truth set. GATK filtering rules are listed in the Supplementary Material.Crumble* refers to parameters optimized for this dataset: ‘crumble -9p8 -u30 -Q60 -D100’. The false positive/negative values of the GATK calls on the crumbleddataset are shown relative to their respective GATK called lossless dataset. The truth set for Chromosome 1 has 269 655 SNPs and 46 036 indels, counting multi-allelic sites once per allele. The quality sizes are absolute for all files.Crumble 339

Crumble: reference free lossy compression of sequence quality values

https://www.repository.cam.ac.uk/bitstream/handle/1810/299052/Crumble%20reference%20free%20lossy%20compression%20of%20sequence%20quality%20values.pdf?sequence=1&isAllowed=y

Crumble: reference free lossy compression of sequence quality values

Abstract

Similar works

Full text

Available Versions

Apollo (Cambridge)