Source-linked AI summary
Exploring Genome Characteristics and Sequence Quality Without a Reference
Jared T. Simpson
TL;DR
De novo assembly of large, complex genomes is challenging, and practical assistance for users has received comparatively little attention. The paper introduces reference-free read-quality and genome-characteristic analyses, including probabilistic classification of graph branches. The resulting software supports quality checks and assembly-difficulty assessment, while its branch classifications could help adapt assembler algorithms to genome structure or data quality.
Problem
De novo assembly of large, complex genomes is challenging, while comparatively little attention has been paid to practical assistance for assembly users.
Method
The software performs reference-free quality assessment and uses a probabilistic classifier to distinguish graph branches caused by sequencing errors, variants, or repeats.
Results
The software provides quality checks and assesses assembly difficulty by measuring de Bruijn graph branching structure; branch classifications could improve assembler performance by adapting algorithms to genome structure or data quality.
Takeaways & Limitations
The software helps users understand their data, select software and parameters, and troubleshoot assembly problems.
Takeaways & Limitations
The classifier has limited power to distinguish sequencing errors from variants when λk is small and outputs no classifications when λk < 10.
Abstract
from arXiv · showhide
The de novo assembly of large, complex genomes is a significant challenge with currently available DNA sequencing technology. While many de novo assembly software packages are available, comparatively little attention has been paid to assisting the user with the assembly. This paper addresses the practical aspects of de novo assembly by introducing new ways to perform quality assessment on a collection of DNA sequence reads. The software implementation calculates per-base error rates, paired-end fragment size histograms and coverage metrics in the absence of a reference genome. Additionally, the software will estimate characteristics of the sequenced genome, such as repeat content and heterozygosity, that are key determinants of assembly difficulty. The software described is freely available and open source under the GNU Public License.
1 Introduction
De novo genome assembly remains difficult and variable, while users receive comparatively little practical assistance. The paper introduces reference-free analyses of graph structure and read quality to better characterize assembly difficulty and input data.
- Motivation: De novo assembly remains a significant challenge, with assembly quality varying substantially across genomes and software submissions.Assemblathon2 produced 43 assemblies from 21 teams, with highly variable quality between submissions and across species.
- Assembly difficulty: Graph-based assemblers model genomes as walks through graphs whose branches arise from repeats, sequence variation, and sequencing errors.These branches obscure the true genomic walks and increase assembly difficulty.
- Approach: The software estimates branch contributions from variants, repeats, and sequencing errors while assessing coverage, per-base error rates, insert sizes, and coverage bias without a reference.These analyses provide a broad quality assessment of sequence reads and assembly-relevant genome properties.
- Implementation: The implementation uses an FM-index to analyze 170 Gbp of human genome data in under 24 hours with 56 GB of memory.It produces a PDF report for visually exploring genome properties and data quality.
2 Results
The analyses characterize genome variation, repeats, coverage, size, and paired-end libraries without prior genome knowledge, revealing factors that shape assembly difficulty. Across the test genomes, oyster is consistently difficult, while yeast is comparatively uncomplicated.
- 2.1 Input data: Six datasets span genomes with varying heterozygosity and repeat content, from straightforward yeast to highly heterozygous, repeat-rich oyster.The human and Assemblathon2 datasets represent large eukaryotic genomes with differing heterozygosity and repeat content.
- 2.2 Exploring heterozygosity: A probabilistic classifier separates de Bruijn-graph branches caused by sequencing errors, sequence variants, or genomic repeats.Variant-induced branches are quantified as a function of k; approximately 1 in 1000 human graph vertices has a variant-induced branch.
- 2.2 Exploring heterozygosity: Oyster has the highest variant-branch density, bird the highest heterozygosity among Assemblathon2 genomes, and yeast less than 10^-4 sequence-variation branch rate.In haploid yeast, the low branching likely reflects misclassified systematic errors or repeats rather than allelic variation.
- 2.3 Exploring genome repeat content: Repeat-branch rates decrease with k; oyster has a repeat rate comparable to human despite being approximately one-fifth its size, while yeast branches infrequently.Fish is more repetitive than expected for its relatively small genome size, whereas yeast appears straightforward to assemble from its low variation and repeat branching.
- 2.6 Simulating contig assembly: Simulated contig N50 generally peaks at an intermediate k, and oyster is again the most difficult dataset to assemble.The peak reflects a balance between resolving short repeats and maintaining graph connectivity; high-depth snake data supports very large k.
- 2.7 Paired-end data and 2.8 Model accuracy: Paired-end insert-size distributions are inferred by graph walks connecting read pairs, with three oyster modes reflecting three separate libraries.The model comparison found repeat branch rates for simulated and real data closely matched the diploid reference, while real-data variant rates were consistently higher.
3 Discussion
The program addresses the practical gap in helping users assess sequencing data and assembly difficulty. It classifies graph branches to distinguish sequencing errors from variants, potentially adapting assembly algorithms to genome structure and data quality.
- Practical assembly support: The program combines data-quality checks with assessment of assembly difficulty from de Bruijn graph branching.It is intended to make genome assembly easier and more consistent for users.
- Branch classification: Most assemblers use heuristics to decide whether graph branches arise from sequencing errors or variants.
- Branch classification: Branch classification models could replace these heuristics and may improve assembler performance by adapting algorithms to genome structure or data quality.
4 Methods
The methods use an FM-index to query read collections, sample reads, and implicitly represent de Bruijn graph structure. They estimate coverage and genome size from k-mer counts, with corrections for sequencing errors and safeguards for high-error or highly heterozygous data.
- FM-index framework: The FM-index counts occurrences of strings in a read collection and supports random extraction of complete reads.Sampling extracts a read after drawing an index uniformly from the n reads.
- De Bruijn graph representation: The de Bruijn graph uses k-mers as vertices, connects k-mers by k −1 overlaps, and treats reverse complements as the same vertex.Graph neighborhoods are recovered through k-mer count queries on the FM-index.
- De Bruijn graph representation: A vertex is a suffix or prefix branch when it has multiple corresponding neighbors.
- Coverage assessment: Coverage is assessed by sampling 100,000 reads, counting every k-mer, and plotting the frequency of each occurrence count.
- Genome-size estimation: Genome-size estimation relates n(l −k + 1) observed k-mers to the mean count λ_k of unique genomic k-mers.The relationship assumes equal read length and error-free reads before sequencing-error correction.
- Genome-size estimation: Sequencing errors inflate the observed k-mer total, so the basic genome-size calculation overestimates G unless the genomic k-mer proportion is estimated.
- Genome-size estimation: The implementation estimates λ_k from the k-mer-count mode, avoiding error modes and using 2m when a secondary peak at 2m is sufficiently high.For a second peak, the criterion is N_2m ≥ N_m/2.
- Genome-size estimation: The method uses k = 31 and samples 20,000 reads to estimate the 31-mer count histogram.
4.4 Branch classification
The classifier estimates whether de Bruijn graph branches arise from sequencing errors, sequence variants, or repeats. It combines k-mer coverage, branch structure, and probabilistic models to estimate branch rates while excluding low-information cases.
- Branch classification: The method classifies de Bruijn graph branches as caused by sequencing errors, sequence variants, or repeats under a diploid-genome assumption.It processes candidate branches after identifying k-mers expected to occur once and uses probabilistic classification.
- Probabilistic model: The classifier replaces total branch coverage with d, the number of reads containing a neighboring k-mer but not the current k-mer, reducing dependence on current-k-mer coverage.Under variant and error models, d is expected to be small; repeats provide more opportunities for neighboring k-mers to occur elsewhere.
- Coverage evidence: Coverage balance distinguishes variants from errors because variant-supported neighbors should have similar coverage, whereas errors should favor one higher-coverage neighbor.The two highest-coverage neighbors are denoted ka and kb, with ca ≥ cb.
- Probabilistic model: Posterior probabilities combine independent evidence from d and coverage counts using equal one-third priors for the three classifications.The outputs are expected branch counts Ne, Nv, and Nr, plus Nh, the expected number of checked homozygous k-mers.
- Output: Branch rates are reported as Nr/Nh and Nv/Nh across k = 21 to 71 in increments of 5, using 50,000 randomly sampled reads.Values are omitted when the expected number of branches for a classification is below 2, and classifications are not output when λk < 10.
4.5 Estimating per-base error rates
Per-base error rates are estimated without a reference by comparing reads that share exact 31-mer seeds, constructing local multiple alignments, and identifying base calls that disagree with supported consensus bases.
- Overlap detection: Read overlaps are seeded by shared 31-mers and retained only when they are at least 50 bp long with at least 95% identity.The retained overlaps are used to construct a multiple alignment around a sampled read.
- Error identification: A base is marked incorrect when it disagrees with the consensus, at least three reads support that consensus, and fewer than four support the original call.The rule uses the multiple alignment column to update the estimated error rate.
4.6 Calculating per-base quality scores
Average per-base quality scores are estimated by sampling reads from the input FASTQ file rather than evaluating every read.
- Quality-score sampling: The average quality score per base is computed from every 20th input FASTQ read, up to a maximum of 10,000,000 reads.This sampling procedure defines the data used for the quality-score estimate.
4.7 GC vs coverage plots
GC-versus-coverage plots summarize the relationship between read GC content and 31-mer count using sampled reads and two-dimensional histograms. Reads whose first 31-mer is likely erroneous are excluded.
- Histogram construction: The plots encode sampled reads as two-dimensional histograms of GC proportion and first-31-mer count.The count is obtained with countDNA, while GC proportion is calculated across the entire read.
- Read filtering: Reads with first-31-mer count one are rejected because that k-mer is likely to represent a sequencing error.The remaining GC-count pairs are passed to a two-dimensional histogram function.
4.9 Simulating de Bruijn assembly
The study simulates de Bruijn graph assembly by extending sampled reads through graph neighborhoods, while testing branch classification against a diploid reference and real-data pipeline reproducibility.
- Simulation procedure: Assembly walks begin from sampled reads whose first k-mer has at least a 0.50 probability of being homozygous.Reads below this threshold are discarded before contig construction.
- Simulation procedure: Contig extension follows the highest-coverage neighbor unless a repeat branch, missing neighbor, or termination condition stops the walk.The procedure extends in both suffix and prefix directions.
- Simulation procedure: 20,000 walks are performed for each k from 21 to 91 in increments of 5, with maximum walk length capped at 50,000.Repeated k-mers and previously used walks are excluded using a bloom filter.
- Validation: Accuracy is assessed using simulated reads from the diploid NA12878 reference processed through the same pipeline as real datasets.Branch rates are also calculated directly from the diploid reference graph using count-based classification.
- Reproducibility: The implementation writes genome-characteristic and quality-control results to JSON and generates PDF reports through a reproducible Makefile.The largest human dataset required 18 hours to build the FM-index and reached a 56GB memory highwater mark.
Figures
The figures collectively show genome-branch characteristics, k-mer and GC-dependent coverage, simulated contig continuity, sequencing quality, and paired-end fragment-size estimates.
- Genome graph characteristics: Figures 1 and 2 plot estimated variation and repeat branch rates for each genome as functions of k.Figure 2 notes that yeast estimates stop at k=51 when repeat branches fall below the emission threshold.
- Coverage diagnostics: Figures 3 and 4 visualize 51-mer frequencies and two-dimensional GC-count distributions for 31-mers across the datasets.The GC-count figure includes fish, yeast, and oyster data.
- Simulated assembly: Figure 5 reports N50 lengths of simulated contigs for k values from 21 to 91 in increments of 5.This figure examines how simulated contiguity changes with k.
- Read quality: Figures 6 and 7 show per-base error rates and mean quality scores for each dataset.The error-rate figure derives estimates from read-read overlaps.
- Library and validation analyses: Figure 8 estimates paired-end fragment sizes, while Figure 9 compares variation and repeat branch rates across real, simulated, and reference data.The oyster fragment-size estimate reflects a mixture of three libraries.
Tables
Table 1 presents genome-size estimates produced by the method alongside previously published estimates.
- Genome-size estimates: Table 1 compares genome size estimates from the method with previously published estimates.The supplied table description identifies the comparison but provides no numerical values.
- Genome-size estimates: The table is intended to assess agreement between the method's estimates and published estimates.No individual genome or estimate is specified in the supplied passage.
- Genome-size estimates: Table 1 reports method-derived estimates in the context of prior genome-size measurements.The supplied passage does not include the table rows or values.