estimate genome size: kmer-based approach from PacBio reads

Question

Can anyone suggest a software/method for kmer analysis using PacBio reads (RSII)?

Something similar to Jellyfish, that I saw in a nice tutorial - but must be suitable for long, noisy reads. kmercounterexact from BBMapTools might also be an option, but again tricky with long, noisy reads.

The closest solution I could find was to estimate the genome size from overlap statistics (for the raw reads, the error corrected preassembled reads or by mapping the raw reads to the assembled contigs) by:

1. running an initial assembly using a small seed read length
2. then plot the pre-assembled read overlap histogram.

What do you think of this approach?

I would like to use this kmer analysis to estimate genome size, but for now I could not find much. I do not have companion Illumina reads - if so, I could try the approach proposed by Carvalho et al. 2016

I have reads & in-house genome assemblies (Canu) from different isolates of the same species of fungus to compare.

Answer 1

I don't think there is a method that would estimate a genome size using raw long reads.

The genome size estimates based on raw reads are done by fitting a model to kmer spectra (for instance Genomescope). The kmer spectra built from long reads are really messy due to the high error rate of long reads. That makes fitting of a model quite difficult. These methods assume small error rates. If the error rate is as high as 0.1 (which is normal for long reads), then the probability of having one kmer of 23 bases correct is ~8%. Meaning that if the coverage is 100x, just ~8x of homozygous kmers would be the real genomic kmers and it's very likely that kmer peaks of the sequencing errors, heterozygous kmers and homozygous kmers will just blend into one.

I can imagine, that if your coverage would crazy high (like 600x or more) the kmer spectra could separate errors and genomic kmers. You could also try to decrease kmer size to increase the proportion of correct kmers.

Another option would be to use a trick that they used to generate solid kmers for genome assembly using generalized de brujin graph using noisy reads. Maybe it would be better to take a look on the newer version of the assembler, that is called Flye.

But to be honest, I think for a long reads, a genome assembly is much more reliable way how to predict a genome size than any of the kmer tricks. I would simply assemble the genome, remap back reads, check the uniformity of coverage and conclude the genome size out it.

-- edit --

I tried to find some comprehensive evaluation of the assembly sizes compared to genome sizes, but I could not, so I just checked two examples, one big, one small and they both match quite well. Ammopiptanthus nanus was assembled to 823.74 Mb while kmer estimate was 889 Mb (ref). Plasmodium falciparum was assembled to 23.5 Mb while the original genome size estimate is 22.8 Mb. I know it's just two, but there is also a good reason to think that the genome size will be reasonable, feel free to add here more examples.

If you would be worried about contamination in the assembly, you can first run it though one of the methods to detect a contamination (like Blobology).

Answer 2

I have recently come across this submission from Wang et al. (2020), which is a method using the tool kmerfreq on error-corrected long reads (e.g. generated using Canu, FLAS, or LoRMA). It looks to be a promising approach.