8
$\begingroup$

I have a Fastq file and I want to get the exact count of each possible kmer from this file.

On a previous post called How to use khmer to count k-mers? Daniel Standage proposed a custom script based on khmer methods that you can see down below :

>>> outfile = open('outfile.txt', 'w') 
>>> seenkmers = set()  # Consumes a lot of memory for large input!!!
>>> for read in khmer.ReadParser('reads.fq.gz'):
...     for kmer in counts.get_kmers(read.sequence):
...         if kmer not in seenkmers:
...             print(kmer, counttable.get(kmer), file=outfile)
...             seenkmers.add(kmer)

The code works if I switch counttable.get(kmer) with counts.get(kmer) but the main problem is that when I try on a customized fastq file with 300 occurrences of one kmer, the count will still be stuck at 255 and can't go further, and for some reason I can't pass through this limit.

Is this coming from the minimum table size ?

Does anyone know if there is a parameter that can work on a custom script to disable this limit ?

$\endgroup$
5
$\begingroup$

255 is the default maximum size of a Counttable in khmer. You want to do the following:

import khmer
counts = khmer.Counttable(31, 1e7, 1)  # Feel free to change the parameters
counts.set_use_bigcount(True)

The last line increases the maximum value from 255 to 65535. I don't think there's any way to currently go above that without changing the khmer code (at some point you have to increase the size of the number used in the hash table).

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.