Re: [Denovoassembler-users] Decreasing assembly quality with increasing kmer values

[Rick W. <wes...@pu...>]

   A peak in effective kmer size for a quality assembly is common in almost all assembly programs.  I see it in ABySS, in Ray, in SOAPdenovo.  This is not to say that you are doing everything correctly -- you may be making an error since you are suspiciously peaking out at your memory constraints -- but rather it is to say that longer kmers simply may not be useful.  In my current project Ray's best assembly is at kmer41.  It is worse at lower kmers and higher kmers.  There is no hard-and-fast rule.
    
   The above is the important point. Below is a long winded explanation which I am sure that other people could do much better than I.

   I recall reading in the Ray manual that longer kmers are simply too divergent to be useful but can not find that passage at the moment.  Let's see if I can explain -- if nothing else the practice is good for me.  Basically for any given kmer the assembly programs are looking for a peak number of reads versus occurrence.  This is roughly the depth of coverage and will be used for further refinement.  As an example, for kmer41 I might get:

kmers occurring 1 time:      20   [basically noise; in the ideal world]
kmers occurring 2 times:     30
kmers occurring 3 times:     20
...
kmers occurring 20 times:   300
kmers occurring 21 times: 2,500
kmers occurring 22 times: 9,000
kmers occurring 23 times: 7,500
kmers occurring 24 times: 1,500
...

Thus my peak is about 22 times.   For a genome assembly program it is assumed that the entire genome will be at this coverage.    The assembly program can then throw away -- or treat differently -- the parts of reads that have kmers in the low occurrence category (these are due to sequencing error or random noise) as well as those with high occurrence (these are repeat areas).   In the above example the assembly program might discards reads with kmers outside the '10 to 40' range.

The above concept works well if there are no sequencing errors.  But of course there are thus kmers with errors become unique and we see the initial part of the occurrence list look more like :

kmers occurring 1 time:  10,000  
kmers occurring 2 times:  3,000
kmers occurring 3 times:    500
...

As we increase the kmer length then the number of kmers with an error in them, and thus unique, increases. More and more kmers get put into the 'occurring 1 time' (or 2 times) category and are thus discarded or at least not used in the initial stages of assembly.  Additionally the peak gets flattened and harder to determine.   Thus we want want to balance kmer size with the need not to have too short of a kmer.  

As an example, let's say that sequencing error is 1 in 100 bases (with some distribution function).  If read sizes are 100 bases and the kmer size is 31 then There will be 69 kmers from this reads. What is the chance that any given kmer will have an error in it and thus be unique?  It does depend on the distribution but we can say (very!) roughly about 25% (if I have that right).   How about if the kmer size is 90?  Then roughly 90% of the kmers will have an error and thus be unique.  At this point we might as well not try to assemble anything. 

Bottom line: longer kmers are not always better.  


----- Original Message -----
> Just wondering if you might be able to offer some insight as to what
> might be going with our runs of ray. We are assembly a large
> eukaryotic genome and have been using both Ray and Celera to do so.
> Initially, Ray gave us some great results and increasing kmer values
> up to 31-kmer readily improved our assemblies. We then ran into a
> memory wall and got some folks to help us out running things off site.
> We were able to run a 33-kmer run on the same set, but it looked worse
> than the 31-kmer run. Additionally, longer runs of say 61-kmers have
> looked even worse. Genome should be around 30-40x coverage, I think,
> with Hi-seq and Mi-seq data (2x150 and 2x250 paired-end libraries,
> each heavily trimmed to get rid of anything suspicious) and we thought
> we'd still have some room to increase kmer sizes >31 to improve the
> assembly. Any thoughts on what might be going on here? Did we possibly
> compile the program incorrectly for larger kmer sizes? Ever see
> something similar?
> 
> Thanks in advance for any thoughts on the matter.
> 
> Cheers,
> Nate

-- 
Rick Westerman 
wes...@pu...

Bioinformatics specialist at the Genomics Facility.
Phone: (765) 494-0505           FAX: (765) 496-7255
Department of Horticulture and Landscape Architecture
625 Agriculture Mall Drive
West Lafayette, IN 47907-2010
Physically located in room S049, WSLR building