genomescope | Fast genome analysis from unassembled short reads | Genomics library

A: The most common problem is you have too low of coverage for the model to confidently identify the peaks corresponding to the homozygous kmers. This can be recognized by a lack of any peaks in the kmer plots or that the model fit doesnt match the observed kmer profile very well. To correct this problem, first make sure that you have used the cannonically kmer counting mode (-C if you are using jellyfish). If this still fails, you can try slightly decreasing the kmer size used to 17 or 19. If all of these attempts still fail you will unfortunately need to generate additional sequencing data. A: Thank you for writing. I think there is a bit of confusion over the variable names and how they relate to each other. The first thing to note is λ and kcov refer to the same value, just that we use λ in the written document and kcov in the code. The modeling tries to identify 4 peaks centered at λ, 2λ, 3λ, and 4λ. These 4 peaks correspond to the mean coverage levels of the unique heterozygous, unique homozygous, repetitive heterozygous and repetitive homozygous sequences, respectively. So when it estimates the haploid genome size, it divides by the 2λ, which is the average homozygous coverage, not "2 times the estimated coverage for homozygous k-mers" as you write. The other potentially confusing aspect is what is meant by haploid genome size versus diploid genome size. We consider the haploid genome size to mean the span of one complete set of haploid chromosomes and the diploid genome size to be the span of both haploid copies (total DNA content in one diploid cell). In particular, in a human cell, the haploid genome size is about 3Gbp and the diploid genome size is about 6Gbp. If you sequence a total of 300Gbp for a human genome, that would be about 150Gbp (50x coverage) of the maternal haplotype and about 150Gbp (50x coverage) of the paternal haplotype. But since the heterozygosity rate in humans is so low, the main peak in the distribution would be centered around 100x. However, GenomeScope will still try to fit the 4 peaks, so should set the heterozygous kmer coverage λ equal to 50x, and thus the homozygous coverage to 2*λ = 100x. From this GenomeScope will compute the haploid genome size as the total amount of sequence data (300GB) divided by the homozygous coverage (100x) to report 3Gbp as expected. Kmers with higher coverage are naturally scaled as well: kmers that occur 200 or 300 times in the kmer profile (and thus are 2 or 3 copy repeats in the diploid genome, 4 or 6 times in the haploid sequences) are still scaled by 100x to contribute 2 or 3 copies to the estimate. Finally, note that if the two haplotypes have significantly different lengths, then the reported haploid genome size will be the average of the two. A: No, GenomeScope is only appropriate for diploid genomes. In principle the model could be extend to higher levels of polyploidy by extending the model to consider additional peaks in the k-mer profile, but this is not currently supported. GenomeScope also does not support genomes that have uneven copy number of their chromosomes, such as aneuploid cancer genomes or even unequal numbers of sex chromosomes. In these scenarios the reported heterozygosity rate will represent the fraction of bases that are haploid (copy number 1) versus diploid (copy number 2) as well as any heterozygous positions in the other chromosomes. While GenomeScope can automate most of the analysis, it does have one critical parameter on which high frequency kmers should be filtered out. This parameter is needed because we often see in real samples that kmers that occur 10,000 times or 100,000 times or greater are artifacts such as phiX sequencing, or organelle sequencing (or other contamination). To account for these, GenomeScope will by default exclude any kmers that occur more than 1000 times from the analysis, which lead to a genome size estimate of 649Mbp. If you increase this cutoff to 10,000, then the new size estimate is 697Mbp. You could perhaps raise this limit even higher, but it looks like your histogram is truncated at 10,000 (which is the default for jellyfish). If you want to include these ultrahigh frequency kmers, then you will have to regenerate the histogram from jellyfish and then set the max coverage threshold to 100,000 or perhaps even 1,000,000. This will likely increase the estimated genome size, but also make the analysis even more sensitive to any artifacts in the data. Unfortunately every project is a little bit different on how to best remove those artifacts. Please see the paper for details (especially the section in the supplement on characterizing the high frequency kmers in Arabidopsis). The short answer is by mixing the samples together, this becomes a tetraploid sample, but GenomeScope currently doesnt work with ploidy > 2. Ive been doing some simulations to further understand this: I randomly generate a 1Mbp maternal sequence m1, and then introduce heterozygous variants to create the paternal sequence p1 at a given rate of heterozygosity r1. In parallel, I create a second maternal genome m2 that differs from m1 by a separate rate of heterozygosity rM. From m2, I derive the paternal genome p2 using the rate of heterozygosity r2. Note that there are 3 rates of heterozygosity that can be adjusted (r1, r2, and rM), and from your data we know that r1 and r2 are around 0.1% but rM is unknown, although is assumed to be very low because the mitochondrial genomes are identical. From this framework, I tried several values of rM, and see that when rM is set to 0.02% then GenomeScope infers the overall rate of heterozygosity at about 0.05%, similar to your data (see attached results). In the joint GenomeScope plot, there is a new peak of heterozygous kmers centered at around 10x coverage (below the expected coverage for heterozygous variants), but GenomeScope incorrectly assumes these are errors just because it doesnt understand how to deal with higher levels of ploidy. We are working on extended the model to consider situations like this, but for now if you want to demonstrate this further, you'd have to align the reads to your assembly, and then you should detect variants that occur at about 50% allele frequency (where m1 and m2 are different), plus variants that occur at about 25% allele frequency (where there is further heterozygosity in p1 or p2). If you would like to run additional simulations, the code is available in the repo at https://github.com/schatzlab/genomescope/tree/master/analysis/genomesim/polyploid.