Multiple paths per variant

[standage]

Currently kevlar find reports:

• "interesting" k-mers
• reads containing "interesting" k-mers
• linear paths from the assembly graph baited by the "interesting" k-mers

In an ideal world, each variant would be associated with a single contig/path that is long enough to determine its position in the genome uniquely. Practically, that's not happening.

Consider the following example, newly introduced in #16. The case has a 5bp deletion relative to the controls. The kevlar find script found 11 "interesting" k-mers (k=13) in 16 reads, and these 11 k-mers pulled out 4 unique linear paths from the assembly graph.

        GGGTCCCTATTGACCTCTTTACCA
        GGGTCCCTATTGACC
              CTATTGACCTCTTTACCA
TGGGCTCGGGGTCCCTATTGACC

A couple of things to note.

• The 2nd and 3rd contigs are subsequences of the 1st. These are easy to collapse, which is also introduced in Collapsing linear paths, new test #16 (so that the final output of kevlar find shows only two linear paths).
• The 1st and 4th contigs overlap, but for some reason they do not exist as a single linear path in the assembly graph.

This second bullet point is what perplexes me now. All along when I've observed this in the human data, I assumed it was the result of sequencing errors introducing structure in the assembly graph and truncating linear paths. But for this example, I simulated reads with no errors or mutations. Why do either of these cases occur?

[standage]

Increasing k ≥ 15 solves the second case, but not the first.

[ctb]

I wonder if this is caused by false positives? You might try calculating cg.kmer_degree and then taking a look at @camillescott's JunctionCountAssembler code in dib-lab/khmer#1502

[standage]

If it's a CountMinSketch FP, this should be handled by loading reads with interesting k-mers into a separate node graph.

[standage]

Reviving this dormant thread.

As far as I've been able to tell this issue has been caused primarily by sequencing errors. khmer has plenty of scripts and functions to deal with this, but we for a while we couldn't come up with the secret QC sauce. Trimming or error correction was never feasible for the raw data. Variable-coverage trimming had weird and unintended results on the boring-reads-dumped data. Digital normalization would provide at best a marginal improvement with a risk of loss of sensitivity (lower risk but also lower benefit by raising normalization threshold). So alas, the last few months I've devoted time to, among other things, implementing a greedy assembly algorithm from scratch.

In implementing the assembler, we came up with an approach for partitioning the reads based on shared interesting k-mers. This turned out to be much more useful than khmer's partitioning code, and led to a great discovery: running abundance trimming (not variable coverage trimming!) on each individual partition, followed by greedy assembly, produced reliable contigs!

This morning, I was curious and revisited kevlar collect which I had all but deprecated (this uses khmer's linear path assembly). Turns out, running this on individual partitions that have been abundance trimmed also seems to produce reliable contigs. Sigh.

So it looks like the breakthrough, ridiculously simple in retrospect, is to partition interesting reads by shared interesting k-mers, and then do abundance trimming and assembly on individual partitions. Assembly details may or may not make much of a difference.