Model Fit vs k-mer size about genomescope HOT 4 OPEN

Hi Ben, Thanks for your interest. I would say this is great consistency between the two runs: a 1% variance in genome size and a 5% variance in the expected heterozygosity rate is within the limits of the data. I noticed your kmer peaks are truncated at 100,000 which probably contributes to some of this variance, especially for the genome size. If you increase this to 1,000,000 I would be the values to be even closer. For most genomes, we recommend a kmer length of 21 to balance out repetitive elements vs the impact of sequencing errors and heterozygosity. This is more-or-less in agreement with the value Merqury suggests over typical genome sizes (a few megabases to a few gigabases). If you use a shorter kmer, then more of the genome becomes repetitive until the peaks are poorly resolved. If you use a longer kmer, more of the kmers will be unique, although this will increasingly obscure how the heterozygosity is distributed within the genome. For a given amount of heterozygosity observed in the kmers, GenomeScope infers the level of heterozygosity at the nucleotide level by assuming the heterozygosity is distributed at a uniform rate. This was necessary for modeling purposes, but in practice, heterozygosity is not uniformly distributed due to variable selective pressures (coding/non-coding etc) so there can be some variability depending on the kmer length used. As an extreme example, in a typical human genome there is a heterozygous variant every 1000bp. So, if you used 1,000 mers or longer mers almost every kmer in a human genome would look heterozygous but it is ambiguous as to if every nucleotide is heterozygous or only 1 in 1000. As you shrink the kmer length, fewer and fewer kmers will overlap a heterozyous nucleotide until the true value can be found of about 0.1%. Hope this helps! Mike

Since jellyfish imposes a maximum count on the histogram, there might be kmers that occur more than 100,000 times, just they are reported as 100,000 times in the file. The number (and sequence) of kmers that exceed this limit will depend in part on the kmer length -- short kmers are less impacted by the ends of the reads and sequencing errors to tend to have higher coverage. So you might partially stabilize the results by increasing the maximum kmer count to around 500,000 times so all kmers are considered in both runs. And for your second question, the relative heights of the peaks directly depends on the choice of k since a single heterozygous base leads to 2*k heterozygous kmers (k from the maternal haplotype and k from the paternal). Consequently, the heterozygous peak for 31mers will be about 50% taller than it will be for 21mers from the same data (it is actually a little bit more subtle than this since changing k will also change the average kmer coverage). The GenomeScope model knows how to correct for this, which is why the estimates different by only a few percent and not 50%. I would encourage you to read the GenomeScope 1.0 supplement where these factors are discussed - you may even want to implement Supplemental equations 2 and 3 in R/python so you can play with how different parameters impact the shapes of the curves. The nls_4peak function here will get you started: https://github.com/schatzlab/genomescope/blob/master/genomescope.R Hope this helps Mike