Source-linked AI summary
Estimation of genomic characteristics by analyzing k-mer frequency in de novo genome projects
Binghang Liu, Yujian Shi, Jianying Yuan, Xuesong Hu, Hao Zhang, Nan Li, Zhenyu Li, Yanxiang Chen, Desheng Mu, Wei Fan
TL;DR
Fragmentary assemblies make genome size and repeat content difficult to estimate, especially for repetitive or heterozygous genomes. This paper models k-mer frequencies in raw sequencing data to estimate genomic characteristics and reports improved accuracy over existing approaches, with k-mer analysis offering an assembly-independent framework.
Problem
Fragmentary assemblies make accurate genome-size and repeat-content estimation difficult, particularly for highly repetitive, polyploid, or heterozygous genomes.
Method
The paper introduces a framework that models k-mer-frequency distributions from raw sequencing data to estimate genome size, repeat structure, and heterozygous rate.
Results
K-mer frequency analysis showed a significant accuracy advantage over traditional experimental methods, while new techniques could increase estimation accuracy by one or two magnitudes.
Takeaways & Limitations
K-mer frequency analysis provides a general, assembly-independent approach for estimating genomic characteristics from sequencing data.
Takeaways & Limitations
Applying the method to real sequencing data remains a major challenge.
Abstract
from arXiv · showhide
Background: With the fast development of next generation sequencing technologies, increasing numbers of genomes are being de novo sequenced and assembled. However, most are in fragmental and incomplete draft status, and thus it is often difficult to know the accurate genome size and repeat content. Furthermore, many genomes are highly repetitive or heterozygous, posing problems to current assemblers utilizing short reads. Therefore, it is necessary to develop efficient assembly-independent methods for accurate estimation of these genomic characteristics. Results: Here we present a framework for modeling the distribution of k-mer frequency from sequencing data and estimating the genomic characteristics such as genome size, repeat structure and heterozygous rate. By introducing novel techniques of k-mer individuals, float precision estimation, and proper treatment of sequencing error and coverage bias, the estimation accuracy of our method is significantly improved over existing methods. We also studied how the various genomic and sequencing characteristics affect the estimation accuracy using simulated sequencing data, and discussed the limitations on applying our method to real sequencing data. Conclusion: Based on this research, we show that the k-mer frequency analysis can be used as a general and assembly-independent method for estimating genomic characteristics, which can improve our understanding of a species genome, help design the sequencing strategy of genome projects, and guide the development of assembly algorithms. The programs developed in this research are written using C/C++, and freely accessible at Github URL (https://github.com/fanagislab/GCE) or BGI ftp ( ftp://ftp.genomics.org.cn/pub/gce).
Background · Methods
The paper develops an assembly-independent k-mer frequency framework to estimate genome size, repeat structure, and heterozygous rate despite incomplete assemblies and complex sequencing data. Its methods introduce k-mer individuals, float precision estimation, and treatments for sequencing error and coverage bias, with applications to sequencing-strategy and assembler selection.
- Background: Highly repetitive, polyploid, and heterozygous genomes make short-read assembly difficult, often producing incomplete and fragmental sequences.These characteristics also make genomic characteristics difficult to estimate properly for non-bioinformaticians.
- Background: The framework improves k-mer-based genomic estimation by introducing k-mer individuals and float precision estimation while accounting for sequencing error and coverage bias.It extends estimation beyond genome size and repeat structure to heterozygous rate.
- Counting k-mer frequency: K-mer frequencies are counted with available tools or the in-house Kmerfreq software, using a k-mer size small enough to limit memory use but large enough for mostly unique genomic k-mers.The maximum number of k-mers is fixed as 4^K once K is selected.
- Estimating sequencing depth and genome size: Under ideal random, error-free, unbiased sequencing assumptions, sequencing depth and genome size are derived from total k-mer count and expected k-mer coverage depth.The model infers expected coverage depth from the Poisson-distributed k-mer frequency curve.
- Estimating sequencing depth and genome size: The k-mer individuals curve is introduced as an equivalent coverage-depth representation whose points weight k-mer species by their depth values.It has the same shape as a Poisson distribution but shifts rightward by one unit, and can be combined with the species curve for estimation.
- Estimating sequencing depth and genome size: Float precision estimation uses neighboring points around the coverage-depth peak because the real expected depth is often non-integer, unlike integer precision estimation.The method can average estimates from 5 to 10 adjacent point pairs for greater accuracy.
- Exploring repetitive genomes: For repetitive genomes, repeat-derived k-mers generate multiple peaks in species and individuals curves, with peak heights related to repeat-class ratios ai and bi.The species and individuals distributions are modeled as compound discrete Poisson and varied-Poisson distributions, respectively.
- Exploring the heterozygous genomes: The framework estimates heterozygous rate from SNP-restricted heterozygous and homozygous k-mer classes and addresses sequencing error and coverage bias to improve genomic estimation.Heterozygosity is simulated by introducing random SNP sites into an ideal haploid reference to create an ideal diploid genome.
Results and discussion
Simulations showed that float-precision k-mer methods improve genome-size estimation, while repeats, heterozygosity, sequencing errors, and k-mer size affect accuracy. Tests on simulated and real sequencing data demonstrated accurate genomic-characteristic estimates with practical value for assembly guidance.
- Coverage and precision: Float precision reduced ΔG/G to 0.01%-0.1%, versus 0.1%-10% for integer precision, with roughly two orders of magnitude advantage below ck-mer 20.The advantage remained roughly one order of magnitude when ck-mer was above 20 but below 85.
- Simulated-data accuracy: About 15%, 73%, and 95% of simulated ΔG/G values were smaller than 0.1%, 1%, and 5%, respectively.About 13%, 54%, and 92% of Δa1/a1 values were smaller than the same thresholds.
- Influencing factors: Higher repeat, heterozygosity, and sequencing-error rates lowered genome-size accuracy, with sequencing error the most difficult factor, especially in combination with the others.Smaller k-mer sizes were suggested for data with high heterozygosity and sequencing-error rates.
- Repeat-structure estimation: Estimated and theoretical a1 and b1 values were consistent across k-mer sizes, supporting accurate repeat-structure estimation and helping select k-mer sizes for de novo assembly.The a1 and b1 values relate to unique k-mer species and the genome fraction covered by unique k-mers.
- Heterozygosity estimation: Estimated heterozygous-rate magnitudes matched simulations except at extremely low heterozygosity, while repeat content and sequencing errors reduced accuracy.The effects were especially pronounced when the theoretical heterozygous rate was very low.
- Real-data evaluation: Over 80% of ΔG/G or Δa1/a1 values from real sequencing data were smaller than 5%, and the estimates retained the correct trend for guiding de novo assembly.Estimated genome sizes from error-corrected data tended to be smaller than those from raw and low-quality-filtered data.
Conclusion
The framework uses k-mer frequency analysis of raw sequencing data to estimate genome size, repeat structure, and heterozygous rate, while also informing sequencing characteristics. It provides practical programs and supports sequencing-strategy and assembly-algorithm decisions, although accurate estimation for complex real genomes remains challenging.
- Contributions: The framework estimates genome size, repeat structure, and heterozygous rate from raw sequencing data, with potential extensions to polyploidy and DNA contamination.K-mer frequency curves also provide information about sequencing error rate and coverage bias.
- Applications and limitations: The programs are freely accessible, and the method can help determine sequencing strategies and guide assembly-algorithm development.Accurate estimation for complex genomes from real sequencing data remains a challenge requiring better understanding and processing of sequencing characteristics.
- Accuracy: K-mer frequency analysis is more accurate than traditional approaches when sequencing coverage exceeds 30X rather than remaining below 10X.The conclusion contrasts higher coverage commonly produced by next-generation sequencing with the lower coverage required by traditional Sanger projects.
- Accuracy: New k-mer individuals and float precision methods could increase estimation accuracy by one or two magnitudes over rough integer-precision estimation.The model achieved very high accuracy on simulated sequencing data, while coverage, repeats, heterozygosity, errors, and coverage bias affected accuracy.
- Real sequencing data: Both standard and continuous models can estimate genome size and ai values under coverage bias, and real-data accuracy is sufficient for many applications.Error reduction lowers memory use and simplifies estimation, but preprocessing can introduce systematic bias that affects accuracy.
Authors' contributions
WF and BL designed the study and drafted the manuscript, while BL, YS, and JY performed the statistical analysis and wrote the programs. Other authors contributed to discussion, result confirmation, manuscript revision, and final approval.
- WF and BL designed the study and drafted the manuscript.
- BL, YS, and JY performed the statistical analysis and wrote the programs.
- XH, YT, HZ, NL, ZL, YC, JL, DM, and SL participated in discussion, confirmed the results, and revised the manuscript.
- All authors read and approved the manuscript.
Figures
The figures illustrate k-mer frequency models for ideal, repetitive, heterozygous, and error-affected genomes, and examine how sequencing conditions and model precision affect genomic estimation. Simulated and real-data figures compare estimated curves and accuracy across genomic structures, k-mer sizes, and sequencing biases.
- Ideal genome: Figure 1 shows base, k-mer species, and k-mer individuals distributions from 23.8X simulated sequencing of a 10-Mb ideal genome.The read length was 100 bp, k-mer size was 17, and expected k-mer coverage depth was 20.
- Ideal genome: The ideal-genome curves matched theoretical distributions, while k-mer individuals equaled k-mer species multiplied by corresponding depth.A Poisson distribution with expected coverage depth 12.6 produced an observed peak at 12, approximately 5% lower.
- Heterozygous model: Figure 3 varies heterozygous rates from zero to 5% for ideal and human repetitive genomes using k-mer size 17 and coverage depth 20.Both k-mer species and individuals curves were plotted from simulated reads under these heterozygous conditions.
Tables
The tables summarize reference and estimated genome sizes for four species and document how reference construction and reporting reliability affect interpretation. They also compare GCE and GSP estimation accuracy using simulated data.
- Reference genomes: Table 1 lists genome sizes and reference a_i values for four reference species.The accompanying note specifies that gap sequences were removed and reports the first five 17-mer a_i values without heterozygosity.
- Real sequencing data: Table 2 reports reference and estimated genome sizes for real sequencing data across species, datasets, and models.The a_1 values were calculated with k-mer size 17 from incomplete assembled genomes that usually lack repeat sequences.
- Interpretation: Previously reported Ecoli and Yanhuang genome sizes are considered most reliable, but all reported values should still be treated cautiously.The reliability assessment reflects the near-completeness of their reference genomes and the stated inaccuracy of reported values generally.
- Accuracy comparison: Table 3 compares estimation accuracy between GCE and GSP using simulated 25X reads from four reference genomes without heterozygosity.The simulations used k-mer size 17bp; one condition had no sequencing errors, while another used a 1% error rate.
Supplementary figures
The supplementary figures evaluate k-mer curve agreement and show how estimation errors vary across simulated genomic and sequencing conditions. They also compare estimation curves for real sequencing data.
- K-mer curve validation: For 48 data sets, estimated k-mer species and individuals curves are well consistent with real curves in most subfigures.Figure S1 compares 17-mer and 25-mer species and individuals curves using real and estimated data.
- Simulation-based error distributions: Figures S2 and S3 show ΔG/G and Δa1/a1 distributions for k-mer sizes 17 and 25 across heterozygous-rate and error-rate combinations.The distributions cover four reference species using simulated data.
- Factors affecting accuracy: Figure S4 separates ΔG/G distributions across k-mer size, heterozygous rate, error rate, and repeat content for 96 analysis sets.The X-axis uses a logarithmic scale, while the Y-axis shows the accumulated number of data points.
- Real-data model comparison: Figure S6 compares real, standard-model, and continuous-model k-mer species and individuals curves for five species with real sequencing data.The curves use K=17.
Supplementary tables
The supplementary tables provide reference a_i values, simulated heterozygous-rate estimates, and detailed real-data results for genome-size and k-mer parameter estimation. They also document data-processing categories and compare GSP with GCE.
- Reference values: Table S1 reports reference a_i values calculated from genomic sequences for genomes with various heterozygous rates.Only a_i values for i<=5 are shown; a[i] is equivalent to a_i and is used throughout the other tables.
- Simulated-data estimation: Table S2 presents heterozygous-rate estimates from simulated sequencing data using a[1/2] values under the heterozygous model.The estimated heterozygous rate is calculated from a[1/2] using formula (11) in the main text.
- Real-data estimation: Table S3 gives detailed real-data estimates of genome size G and a1 obtained with different methods.It includes three sequencing-data types: all, corrected, and filtered.
- Real-data estimation: The real-data categories distinguish raw k-mer counting, error-corrected counting, and low-quality k-mer filtering.“All” uses all sequencing data, “corrected” uses error-corrected data, and “filtered” removes low-quality k-mers from raw data.
- Method comparison: Table S5 compares estimation accuracy between GSP and GCE for genome size, a1 values, and intermediate results.The table contains detailed results for both methods and their middle results.
Supplementary methods … 3 Exploring the repetitive genomes
The supplementary methods describe k-mer counting, estimation of expected k-mer coverage depth, and repeat-genome analysis. They recommend k-mer-space sizing and use k-mer frequency curves and a compound Poisson model to characterize repetitive genomes and estimate genome size.
- 1 k-mer frequency counting: K-mer frequency can be counted from either a reference genome sequence or sequencing reads as a fundamental bioinformatics analysis.
- 1 k-mer frequency counting: The k-mer size K should provide a k-mer space 4^K several times larger than genome size G so most genomic k-mers remain unique.In practice, the method requires 4^K > 5*G, with larger spaces preferred.
- 1 k-mer frequency counting: The jellyfish algorithm is reported to offer faster running speed, smaller memory requirements, and k-mer sizes up to 31 bp.
- 1 k-mer frequency counting: K-mer counting can require very large memory, particularly for large genomes and data containing many sequencing errors.The authors note that memory-saving algorithms such as bloom filters may be needed when available memory is insufficient.
- 1 k-mer frequency counting: The authors developed kmerfreq with array-based and hash-based counting methods, matching array counting to smaller K and hash counting to larger K.“kmer_freq_array” uses 4^K computer memory, whereas “kmer_freq_hash” scales with the number of k-mer species; the text recommends array counting for K<=17 and hash counting for K<=27.
- 2.1 Estimating c based on definition: Expected coverage depth c is defined separately for genomes without repeats and genomes with repeats using k-mer species and individual counts.
- 2.2 Estimating c from Poisson distribution: For repeat-free genomes, k-mer species frequencies follow a Poisson distribution, providing a formula-based route to estimate c.
- 3 Exploring the repetitive genomes: Repeats create multiple peaks in k-mer species and individual coverage curves, while the ai and bi values jointly provide a more quantitative repeat characterization.Genome size for repetitive genomes is calculated using a compound Poisson model.
4 The standard model · 5 Exploring the heterozygous genomes
The standard model represents k-mer depth as a compound distribution over repeat copy numbers and sequencing coverage, then iteratively estimates repeat composition and unique-region depth. For heterozygous genomes, hybrid peaks at fractional coverage reveal heterozygosity, but distinguishing peaks and estimating heterozygous rates becomes harder in repetitive genomes.
- 4 The standard model: The standard model partitions the genome into repeat families by copy number and models their k-mer depth distributions under random sequencing.For k-mers with the same genomic frequency, sequencing depth is assumed to follow a Poisson distribution.
- 4 The standard model: It defines a_i as the fraction of k-mer species at genomic frequency i and b_i as the fraction of k-mer individuals at that frequency.Coverage-depth distributions are summarized separately for k-mer species and k-mer individuals.
- 4 The standard model: Uncovered k-mers are estimated with the Lander-Waterman model, using the expected depth of each repeat region to calculate its gap probability.The zero-coverage component cannot be counted directly from reads because some genomic regions remain uncovered.
- 4 The standard model: Because repeat copy number i scales with unique-region depth c, the model uses this relationship to estimate genomic characteristics from k-mer depth.The unique-region depth is denoted c, and repeat-associated depth is proportional to copy number.
- 4 The standard model: The estimation starts from prior and posterior probabilities for repeat frequencies and updates a_i iteratively while removing repeat components from the raw depth distribution.After repeat removal, the remaining unique curve is used to estimate the unique depth c.
- 5 Exploring the heterozygous genomes: Heterozygous sites generate hybrid peaks at i*c/2, and in non-repetitive genomes the c/2 peak rises with heterozygous rate and can replace the c peak as major.These peaks arise in the k-mer coverage-depth distribution.
- 5 Exploring the heterozygous genomes: The k-mer individual curve better exposes the homozygous c peak at relatively high heterozygosity, whereas the k-mer species curve better exposes the heterozygous c/2 peak at relatively low heterozygosity.As heterozygosity increases, hybrid peaks may become difficult to distinguish from repeat peaks, requiring additional biological background.
- 5 Exploring the heterozygous genomes: Formula (11) roughly estimates heterozygous rate in diploid genomes, but it is less accurate for repetitive genomes because a_i lacks a clear relationship with heterozygosity.The method remains usable with repeats in practice, though its accuracy is lower than for non-repetitive genomes.
6 Sequencing error model · 7 Continuous compound model
The sequencing-error model distinguishes erroneous k-mers that are absent from the genome from those that merge with correct k-mers, but is not used because key parameters are difficult to estimate. The continuous compound model represents coverage-biased k-mer depths with Poisson-based distributions and estimates genomic coverage around the unique peak.
- 6 Sequencing error model: Sequencing errors produce low-depth k-mers absent from the genome and erroneous k-mers that exist genomically and increase correct k-mer depth.The model separates these two error classes by their effects on expected depth.
- 6 Sequencing error model: The error model’s parameters depend on sequencing error rate, k-mer size, and genome size, making them difficult to estimate in practice.Consequently, the sequencing-error model was not used later.
- 6 Sequencing error model: The method reduces erroneous k-mers by filtering counted k-mers according to their sequencing-quality-based probability of being correct.Sequencing quality is represented on the Phred scale.
- 6 Sequencing error model: K-mers are retained only when P(Kmer, correct) exceeds a threshold such as 0.95.This threshold is intended to filter most erroneous k-mers.
- 7 Continuous compound model: Coverage depth depends on both genomic k-mer frequency and sequencing bias, including uneven GC content and specific sequence structure.Therefore, observed genomic frequency is treated as continuous in real sequencing data.
- 7 Continuous compound model: Assuming independent sequencing within groups sharing an observed genomic frequency, the model assigns each group a Poisson depth distribution.Integrating these distributions yields continuous compound models for total k-mer and k-mer-individual coverage depth.
- 7 Continuous compound model: In practice, a dense discrete model approximates the continuous compound model, using peak-order terms and iteration formulas.The symbol k denotes peak order here and is unrelated to kc.
- 7 Continuous compound model: Genomic coverage depth c can be estimated from kc at the unique peak or from the average kc when most unique k-mers are bias-free.The relationship between ka and ia remains difficult to estimate because sequencing coverage bias cannot currently be modeled.
8 Dealing with sequencing error and coverage bias … 2 Analyzing real sequencing data from de novo genome projects
The framework addresses sequencing errors and coverage bias with an experience-based filtering method, while simulations and real-data analyses characterize estimation accuracy and practical limitations. Results support assembly-independent genomic characterization but show reduced accuracy under challenging sequencing conditions.
- 8 Dealing with sequencing error and coverage bias: A united model represents erroneous and correct k-mers together, but genomic and sequencing characteristics make the model difficult to resolve.The method therefore uses an experience-based procedure to exclude erroneous k-mers caused by sequencing errors.
- 8 Dealing with sequencing error and coverage bias: Coverage bias may alter peak position and distribution shape, yet the major peak still roughly reflects unique k-mer coverage depth in most real data.The analysis recalculates low-depth k-mer missing ratios using a compound Poisson model during Bayes iteration.
- 1 Estimating with various types of simulated data: 8.11% is the largest ΔG/G among 96 simulated analyses, observed for maize with 1% heterozygosity, 1% error, and 25 k-mer size.Sequencing error was the most difficult factor because higher error rates hindered exclusion of erroneous k-mers.
- 1 Estimating with various types of simulated data: Larger k-mers increased ΔG/G with sequencing error or heterozygosity but improved estimation when neither problem was present by better handling repeats.Larger k-mers also lower ck-mer, which tends to reduce estimation accuracy under error or heterozygosity.
- 1 Estimating with various types of simulated data: 12.6% is the largest Δa1/a1 among 96 simulated analyses, observed for maize with 0.1% heterozygosity, 1% error, and 17 k-mer size.Estimating ai is more complex because repeat and heterozygosity effects interact with the chosen standard or heterozygous model.
- 2 Analyzing real sequencing data from de novo genome projects: Filtering low-quality k-mers reduces erroneous k-mers, facilitates estimation of correct k-mer individuals, and significantly decreases computer-memory consumption.For real datasets, approximately 22% and 80% of ΔG/G values were below 1% and 5%, respectively; the largest was 10% for ant using raw k-mers.
- 2 Analyzing real sequencing data from de novo genome projects: Approximately 17% and 80% of real-data Δa1/a1 values were below 1% and 5%, respectively, with a largest value of 11% for ant low-quality filtered data using the standard model.Differences between estimated and reference a1 (b1) values were attributed mainly to coverage bias and uncertainty relationships between observed and genomic ai values.