<p>I have recently come across this submission from <a href="https://arxiv.org/abs/2003.11817" rel="nofollow noreferrer">Wang et al. (2020)</a>, which is a method using the tool <a href="https://github.com/fanagislab/kmerfreq" rel="nofollow noreferrer">kmerfreq</a> on error-corrected long reads (e.g. generated using Canu, FLAS, or LoRMA). It looks to be a promising approach.</p>

CC BY-SA 4.0
2020-04-19T18:34:18.857Z
2020-04-19T18:34:18.857Z
Tane
<p>I don't think there is a method that would estimate a genome size using raw long reads.</p>

<p>The genome size estimates based on raw reads are done by fitting a model to kmer spectra (for instance <a href="https://github.com/schatzlab/genomescope" rel="nofollow noreferrer">Genomescope</a>). 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.</p>

<p>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.</p>

<p>Another option would be to use a trick that they used to generate solid kmers for <a href="https://www.ncbi.nlm.nih.gov/pmc/articles/PMC5206522/" rel="nofollow noreferrer">genome assembly using generalized de brujin graph using noisy reads</a>. Maybe it would be better to take a look on the newer version of the assembler, that is called <a href="https://github.com/fenderglass/Flye" rel="nofollow noreferrer">Flye</a>.</p>

<p>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.</p>

<p>-- edit --</p>

<p>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. <em>Ammopiptanthus nanus</em> was assembled to 823.74 Mb while kmer estimate was 889 Mb (<a href="https://academic.oup.com/gigascience/article/7/7/giy074/5039704" rel="nofollow noreferrer">ref</a>). <em>Plasmodium falciparum</em> was assembled to <a href="https://wellcomeopenresearch.org/articles/3-52/v1" rel="nofollow noreferrer">23.5 Mb</a> 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.</p>

<p>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 <a href="https://github.com/blaxterlab/blobology" rel="nofollow noreferrer">Blobology</a>).</p>

CC BY-SA 4.0
2018-07-24T11:45:29.373Z
ok, I can get genome size from assembly, as discussed. thanks for the input
aechchiki
CC BY-SA 4.0
2018-08-25T06:32:44.430Z
Using a genome assembly as a way to estimate the genome size is not reliable due to potential contamination causing assembly errors, possible failure to collapse some haplotypes, or assembly errors from inconsistent coverage.
conchoecia
CC BY-SA 4.0
2018-08-25T10:11:11.377Z
@conchoecia I do agree, but the contamination problem is in kmer spectra analysis as well. But in my experience, long read assemblies I have seen so far were quite close to flow cytometry estimates and consistent with kmer spectra analysis. I will do some more research. Do you have any example where the assembly size was very biased?
Kamil S Jaron
CC BY-SA 4.0
2018-08-26T20:17:52.023Z
@conchoecia better?
Kamil S Jaron
CC BY-SA 4.0
2018-08-27T20:14:13.660Z
Hi @KamilSJaron, the only caveat that I had in mind is that both canu and FALCON do not handle haploid assemblies well if the animal is diploid and highly heteryzygous, like over 5% in the case of FALCON. For example, I am assembling ctenophore genomes which regularly have uncollapsed haplotypes due to their high heterozygosity. As a result, the genome assembly size is approximately twice what is expected based on a kmer spectrum histogram from Illumina reads. Most people probably won't be working with high-heterozygosity animals, but it is important to point out. That's all!
conchoecia
CC BY-SA 4.0
2018-08-27T20:58:12.670Z
Ah, got you, I ve been there. I was assembling amphioxus genome that IS extremely heterozygous. I ended up setting the parameters so both haplotypes got assembled separately and it worked very well (assembled precisely the diploid genome size)
Kamil S Jaron
CC BY-SA 4.0
2018-07-23T11:12:42.180Z
2018-08-26T20:17:35.090Z
2018-08-26T20:17:35.090Z
Kamil S Jaron
<p>Can anyone suggest a software/method for kmer analysis using PacBio reads (RSII)? </p>

<p>Something similar to <a href="http://www.cbcb.umd.edu/software/jellyfish/" rel="noreferrer">Jellyfish</a>, that I saw in a nice <a href="https://bioinformatics.uconn.edu/genome-size-estimation-tutorial/#" rel="noreferrer">tutorial</a> - but must be suitable for long, noisy reads. <a href="https://github.com/BioInfoTools/BBMap/blob/master/sh/kmercountexact.sh" rel="noreferrer">kmercounterexact from BBMapTools</a> might also be an option, but again tricky with long, noisy reads.</p>

<p>The closest solution <a href="http://pb-falcon.readthedocs.io/en/latest/tutorial.html#assembly-graph-and-pread-overlaps" rel="noreferrer">I could find</a> 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: </p>

<ol>
<li>running an initial assembly using a small seed read length</li>
<li>then plot the pre-assembled read overlap histogram.</li>
</ol>

<p>What do you think of this approach?</p>

<p>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 <a href="https://genome.cshlp.org/content/26/12/1710.full.pdf+html" rel="noreferrer">Carvalho et al. 2016</a></p>

<p>I have reads &amp; in-house genome assemblies (Canu) from different isolates of the same species of fungus to compare. </p>

<genome><software-recommendation><assembly><k-mer><pacbio>
estimate genome size: kmer-based approach from PacBio reads
# estimate genome size: kmer-based approach from PacBio reads

(Asked by: username_0 on 2018-07-23)

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

Something similar to [Jellyfish](http://www.cbcb.umd.edu/software/jellyfish/), that I saw in a nice [tutorial](https://bioinformatics.uconn.edu/genome-size-estimation-tutorial/#) - but must be suitable for long, noisy reads. [kmercounterexact from BBMapTools](https://github.com/BioInfoTools/BBMap/blob/master/sh/kmercountexact.sh) might also be an option, but again tricky with long, noisy reads.

The closest solution [I could find](http://pb-falcon.readthedocs.io/en/latest/tutorial.html#assembly-graph-and-pread-overlaps) 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](https://genome.cshlp.org/content/26/12/1710.full.pdf+html)

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

**Tags:** genome, software-recommendation, assembly, k-mer, pacbio

---

## 2 answers

### Answer by username_2 on 2018-07-23. Score: 5 | Accepted answer

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](https://github.com/schatzlab/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](https://www.ncbi.nlm.nih.gov/pmc/articles/PMC5206522/). Maybe it would be better to take a look on the newer version of the assembler, that is called [Flye](https://github.com/fenderglass/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](https://academic.oup.com/gigascience/article/7/7/giy074/5039704)). *Plasmodium falciparum* was assembled to [23.5 Mb](https://wellcomeopenresearch.org/articles/3-52/v1) 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](https://github.com/blaxterlab/blobology)).

> **Comments:**
> **username_0**: ok, I can get genome size from assembly, as discussed. thanks for the input  
> **username_3**: Using a genome assembly as a way to estimate the genome size is not reliable due to potential contamination causing assembly errors, possible failure to collapse some haplotypes, or assembly errors from inconsistent coverage.  
> **username_2**: @conchoecia I do agree, but the contamination problem is in kmer spectra analysis as well. But in my experience, long read assemblies I have seen so far were quite close to flow cytometry estimates and consistent with kmer spectra analysis. I will do some more research. Do you have any example where the assembly size was very biased?  
> **username_2**: @conchoecia better?  
> **username_3**: Hi @KamilSJaron, the only caveat that I had in mind is that both canu and FALCON do not handle haploid assemblies well if the animal is diploid and highly heteryzygous, like over 5% in the case of FALCON. For example, I am assembling ctenophore genomes which regularly have uncollapsed haplotypes due to their high heterozygosity. As a result, the genome assembly size is approximately twice what is expected based on a kmer spectrum histogram from Illumina reads. Most people probably won't be working with high-heterozygosity animals, but it is important to point out. That's all!  
> **username_2**: Ah, got you, I ve been there. I was assembling amphioxus genome that IS extremely heterozygous. I ended up setting the parameters so both haplotypes got assembled separately and it worked very well (assembled precisely the diploid genome size)  

---

### Answer by username_1 on 2020-04-19. Score: 1

I have recently come across this submission from [Wang et al. (2020)](https://arxiv.org/abs/2003.11817), which is a method using the tool [kmerfreq](https://github.com/fanagislab/kmerfreq) on error-corrected long reads (e.g. generated using Canu, FLAS, or LoRMA). It looks to be a promising approach.


---

