Source-linked AI summary
Exploring single-sample SNP and INDEL calling with whole-genome de novo assembly
Heng Li
TL;DR
The paper asks whether de novo assembly can support variant calling without relying exclusively on read-to-reference mapping. It develops fermi to assemble short reads into information-preserving unitigs, then maps those unitigs to a reference. On human resequencing data, the approach achieves SNP accuracy close to standard mapping and particular strength in INDEL calling, supporting assembly-based calling as a complement to mapping-based pipelines.
Problem
Existing de novo assemblers had not matched the sensitivity of standard mapping-based calling for diploid mammalian genomes, while mapping can mishandle read relationships and difficult reference differences.
Method
Fermi assembles Illumina short reads into unitigs that preserve most input-read information, and variants are called by aligning the unitigs to a reference genome.
Results
Assembly-based calling achieves SNP accuracy close to standard mapping and particular strength in INDEL calling; fermi unitigs had N50 1,022bp and totaled 3.83Gb.
Takeaways & Limitations
Assembly-based variant calling is a practical and beneficial complement to mapping-based calling, with especially strong reported performance for INDELs.
Takeaways & Limitations
The approach depends on correctly constructing unitigs, while large-scale de novo assembly still faces challenges in error correction, read-pair use, memory control, and heuristic decisions.
Abstract
from arXiv · showhide
Motivation: Eugene Myers in his string graph paper (Myers, 2005) suggested that in a string graph or equivalently a unitig graph, any path spells a valid assembly. As a string/unitig graph also encodes every valid assembly of reads, such a graph, provided that it can be constructed correctly, is in fact a lossless representation of reads. In principle, every analysis based on whole-genome shotgun sequencing (WGS) data, such as SNP and insertion/deletion (INDEL) calling, can also be achieved with unitigs. Results: To explore the feasibility of using de novo assembly in the context of resequencing, we developed a de novo assembler, fermi, that assembles Illumina short reads into unitigs while preserving most of information of the input reads. SNPs and INDELs can be called by mapping the unitigs against a reference genome. By applying the method on 35-fold human resequencing data, we showed that in comparison to the standard pipeline, our approach yields similar accuracy for SNP calling and better results for INDEL calling. It has higher sensitivity than other de novo assembly based methods for variant calling. Our work suggests that variant calling with de novo assembly be a beneficial complement to the standard variant calling pipeline for whole-genome resequencing. In the methodological aspects, we proposed FMD-index for forward-backward extension of DNA sequences, a fast algorithm for finding all super-maximal exact matches and one-pass construction of unitigs from an FMD-index. Availability: http://github.com/lh3/fermi Contact: hengli@broadinstitute.org
1 INTRODUCTION
The paper examines de novo assembly as an alternative to mapping-based variant calling, motivated by limitations in how mapping handles read correlations and difficult reference differences. It introduces fermi and reports near-standard SNP accuracy with particular strength in INDEL calling.
- Motivation: Mapping-based variant calling ignores correlations between reads and may produce inconsistent outputs.The standard approach maps raw reads to a reference before detecting differences.
- Motivation: Mapping-dependent improvements remain vulnerable to mismapped or unmapped reads caused by divergence, long insertions, structural variation, copy-number changes, or reference misassemblies.These methods still rely on the initial mapping and therefore do not solve the problem from the root.
- Related approaches: De novo assembly avoids the mapping flaw by discovering variants through assembly-to-assembly alignment.The approach was historically less used for short reads but has re-emerged with longer paired-end reads and improved methods.
- Related approaches: Human-genome assembly is challenging because of genome size, repeats, and diploidy, while existing assemblers have limited heterozygote retention or sensitivity.The introduction notes that assembly-based variant discovery has attracted relatively little attention for human genomes.
- Study objective: The study shows assembly-based calling can approach standard mapping accuracy for SNPs and has particular strength in INDEL calling.The authors present this as the first such demonstration and also describe fermi as capable of human assembly.
2 METHODS
The methods review overlap-based and de Bruijn assembly, then motivate a unitig-based strategy that preserves read relationships while addressing short-read computational costs and practical assembly challenges.
- Overview: The methods section motivates fermi by treating unitigs as a lossless representation of reads and introducing a bidirectional FM-index for DNA sequences.It also presents algorithms for assembly using that index.
- Overlap-layout-consensus: Overlap-layout-consensus assembly can require O(N^2) overlap construction and at least O(E) transitive reduction, limiting its appeal when N approaches 10^9 reads.Later steps are roughly linear in practice, but the graph-building stages remain costly for short-read data.
- Overlap-layout-consensus: Simpson and Durbin’s O(N) algorithm finds all irreducible exact-overlap edges, replacing overlapping and transitive-reduction phases.This makes the overlap-based paradigm practically solvable in time roughly linear in total read length.
- De Bruijn graphs: De Bruijn graphs are attractive for many short reads because they can be constructed in linear time, but they lose long-range read information.Consequently, graph paths may be invalidated by the reads.
- De Bruijn graphs: Read threading can theoretically restore coherence to de Bruijn graphs, but complex repeats make it difficult and potentially inefficient; unitig graphs are simpler to construct coherently.The comparison highlights the trade-off between compact construction and preserving read relationships.
- Practical limitations: Practical de novo assembly still faces challenges in error correction, read-pair use, memory control, and heuristic final assembly.The field remains under active development, and many practical problems are not solved perfectly.
2.2 Rationale
A perfectly constructed unitig graph with per-unitig read counts preserves all read information and excludes paths invalidated by the reads.
- Rationale: A coherent unitig graph with per-unitig read counts encapsulates all information in the reads.The paper therefore treats unitig-based and read-based analyses as equivalent in principle.
2.3 Strings and FM-index
This section defines the string and suffix-index notation used to represent collections of DNA reads and support bidirectional sequence searches.
- Strings: The DNA alphabet contains $, A, C, G, T, and N in predefined lexicographic order, with N denoting ambiguity and $ marking string termination.A string terminated by $ is called a text.
- Strings with multiple sentinels: An ordered collection of reads is represented by concatenating each read with a sentinel, such as ACG$GTG$.The sentinel ordering distinguishes suffixes from different reads.
- Suffix arrays: The suffix array is a permutation listing text suffix positions in lexicographic order, and a string’s suffix-array interval identifies matching suffixes.The interval size is defined as I_s(P) = I_u(P) − I_l(P) + 1.
- FM-index: The BWT stores symbols preceding suffixes, while C(a) counts lexicographically smaller symbols and O(a, i) counts occurrences in a BWT prefix.The FM-index compresses the BWT, occurrence array, and suffix array.
- FM-index: Backward extension updates a pattern’s suffix-array interval using the cumulative count and occurrence arrays.The interval remains nonempty exactly when the extended pattern occurs as a substring of the text.
- Sequence retrieval: A collection can be queried to retrieve a sequence in linear time together with its lexicographic rank.The retrieval procedure is given as Algorithm 1.
2.4 FMD-index
The FMD-index represents forward and reverse-strand DNA sequences in one FM-index and supports bidirectional sequence extension through bi-intervals. Compared with two-index bidirectional BWT, it improves exact-matching speed in the reported example.
- A bi-interval records the interval information for a string and supports computing the bi-intervals of its one-base extensions.
- Forward and backward extensions are dual operations in the FMD-index, making it bidirectional.
- The FMD-index builds both forward and reverse-strand DNA sequences in one index.
- 80% speedup was obtained in BWA-SW after adopting the FMD-index as its data structure.
2.5 Unitig construction
Fermi labels reads and overlaps with integer pairs derived from bi-intervals, then constructs unitigs by unambiguously merging irreducibly overlapping reads. The overlap-finding algorithm uses bidirectional extension and assumes contained reads have been removed.
- 2.5.1 Labeling reads and overlaps: Fermi labels each read with an ordered integer pair [k, l], called its bi-interval, and represents the reverse read with [l, k].
- 2.5.1 Labeling reads and overlaps: Overlap types are encoded by unordered pairs such as ⟨l, k′⟩, with four types corresponding to bidirectional overlap-graph edges.
- 2.5.1 Labeling reads and overlaps: The extension algorithms transform a bi-interval [k, l, s] into the interval for either Wa or aW.
- 2.5.2 Finding irreducible overlaps: Irreducible-overlap detection uses backward extensions to find overlapping reads and forward extensions to exclude reducible overlaps.
- 2.5.2 Finding irreducible overlaps: The overlap algorithm requires that contained reads be absent, although fermi detects containment on the fly in its modified implementation.
- 2.5.3 Unitig construction: Unitig construction merges unambiguously overlapping reads and labels a merge from [k, l] and [k′, l′] with [k, l′].The labeling preserves graph topology without staging the graph in RAM and can be multi-threaded.
2.6 Finding the supermaximal exact matches
Fermi uses the FMD-index to find supermaximal exact matches between sequences and maps reads back to unitigs through these matches. Its algorithm extends exact matches in both directions while tracking interval-size changes.
- A maximal exact match cannot be extended in either direction, while a supermaximal exact match is not contained in another query-sequence MEM.
- Fermi uses SMEMs to map reads back to unitigs.
- The SMEM algorithm uses forward-backward extension to extend an exact match and detect the boundary of a maximal match.
- Fermi’s implementation finds full-length read matches and can optionally exclude matches identical to the query sequence.
2.7 Other implementation details
Fermi combines compressed-index construction, quality-aware error correction, bubble simplification, paired-end information, and reference alignment in its implementation. These procedures support read correction, assembly refinement, and assembly evaluation.
- 2.7.1 Error correction: The BWT is run-length encoded, and suffix arrays and merged BWTs are constructed using established algorithms adapted for multiple read subsets.
- 2.7.1 Error correction: Fermi collects frequent 23-mers, counts their following bases, and uses a quality-aware heuristic to correct sequencing errors in reads.
- 2.7.3 Bubble popping: Fermi simplifies detected bubbles by retaining the top two paths containing most reads and dropping vertices on other paths.
- 2.7.4 Using the paired-end information: At least five linking read pairs trigger local end assembly, while regions lacking bridging pairs can be used to break contigs.
- 2.7.5 Assembly evaluation: Contigs longer than 150bp are aligned to GRCh37 with BWA-SW to identify two types of alignment break points.
3 RESULTS
Fermi produced human short-read assemblies and used unitig-to-reference alignment for SNP and INDEL calling. Its SNP performance was comparable to standard approaches, while its INDEL calls showed broader confirmation and apparent higher overall accuracy.
- Assembly: 35-fold NA12878 reads were assembled in about 5 days using 12 CPUs and 96GB memory, with 92GB peak usage.
- Assembly: Fermi generated unitigs with N50 1,022bp totaling 3.83Gb, breaking mainly at heterozygotes.Collapsing heterozygotes and closing paired-end gaps produced longer contigs.
- Variant calling: Fermi called SNPs by mapping unitigs to the reference and using SAMtools, while INDELs were counted directly from pileup output.A standard INDEL caller was not used because short-read INDEL callers do not work well with long contig sequences.
- Assembly: Among fermi, SGA, and SOAPdenovo on the same dataset, fermi had longer N50 with similar alignment-break-point counts and comparable overall assembly quality.
- SNP calling: Fermi missed 3% of SNPs called in SS but found more additional SNPs, partly because assembled sequences increased BWA-SW alignment power.Some additional calls may reflect misassemblies, reference errors, or copy-number variations and could contain more errors.
- INDEL calling: Fermi’s INDEL call set was slightly smaller than BS but had larger overlap with other call sets and more calls confirmed by others, suggesting higher overall accuracy.Even combined short-read call sets missed as many as 14% of Mills et al. double-hit INDELs; manual checks indicated some of those calls may have high error rates.
4 DISCUSSIONS
The paper presents FMD-index and fermi as a basis for assembly-based variant calling, reporting SNP accuracy comparable to mapping-based calling and stronger INDEL performance. It also frames unitigs as lossless read representations in theory, while acknowledging practical construction challenges.
- FMD-index: FMD-index stores forward and reverse-complement DNA sequences in one FM-index, enabling faster forward-backward search than bidirectional BWT.The authors describe it as a more natural representation of DNA sequences.
- Fermi: Fermi uses FMD-index to assemble reads into unitigs and achieves similar quality to other mainstream assemblers.
- Variant calling: Assembly-based calling achieves similar SNP accuracy to standard mapping-based calling and reportedly improves INDEL sensitivity and precision.The approach aligns assembled unitigs to the reference genome to call SNPs and short INDELs.
- Unitig representation: With read counts retained, unitigs are theoretically a lossless but reduced representation of sequence reads.The representation preserves small variants, copy numbers, and structural changes when unitigs are constructed correctly.
- Limitations: Practical challenges for unitig-based storage and analysis include high computational cost and insufficient unitig accuracy.