Bug in JellyFish and DSK k-mer counting tool?
0
0
Entering edit mode
5.5 years ago
SmallChess ▴ 580

I have some k-mers that neither Jellyfish nor DSK was able to count. Is that a bug in the programs?

This is my script for JellyFish

1. grep TTACATAACACCCATTGTGGCGGCTGCAAGT ABCD.fq | wc

41 41 5166

The 31 length k-mer is in the sample file and there are many of them.

2. jellyfish count -m 31 -s 100M -C ABCD.fq -o 1.jf

3. jellyfish dump 1.jf -c > 1.txt
4. grep TTACATAACACCCATTGTGGCGGCTGCAAGT 1.txt

returns nothing

Then I tried DSK,

1. dsk -kmer-size 31 -abundance-min 0 -file ABCD.fq -out-dir ABCD -out ABCD.h5
2. dsk2ascii -file ABCD.h5 -out ABCD.txt
3. grep TTACATAACACCCATTGTGGCGGCTGCAAGT *.txt

still returns nothing

As you can see, the 31 k-mer is in the input file (ABCD.fq) and there are many of them not being filtered away. I don't see a problem in the input file - it came out from a simulator. It is a valid FASTQ file. Why both programs weren't able to count my k-mer?

dsk jellyfish k-mer • 2.4k views
2
Entering edit mode

Have you tried grepping for the reverse-complement?

Also, maybe you should try the BBMap package's KmerCountExact tool. It will, coincidentally, report that specific kmer by default, because it stores the higher, rather than lower, of each kmer/reverse-kmer pair. The other tools store the lower one.

2
Entering edit mode

DSK indeed will use a canonical form of the kmer that corresponds to the minimum value of the kmer and its reverse-complement. So in the example one should grep ACTTGCAGCCGCCACAATGGGTGTTATGTAA.

0
Entering edit mode

Brian. Thanks for the tip. I will give a try. I will try the BBMap tool and report back.

0
Entering edit mode

Why would both JellyFish and DSK count on reverse-complement?

3
Entering edit mode

It's faster and uses less memory to store either a kmer OR its reverse-complement instead of both, and for most purposes, it makes no difference because if a kmer exists in a genome, its reverse-complement must also exist because DNA is double-stranded. So, counting only one or the other cuts memory use in half.

2
Entering edit mode

Counting 'canonical' k-mers makes sense if your input genome assembly is incomplete..

jellyfish count --help
..
-C, --canonical                          Count both strand, canonical representation (false)
..