<p>I wrote a command-line k-mer counter called <code>kmer-counter</code> that will output results in a form that your Python script can consume: <a href="https://github.com/alexpreynolds/kmer-counter" rel="noreferrer">https://github.com/alexpreynolds/kmer-counter</a></p>

<p>You can grab, build and install it like so:</p>

<pre><code>$ git clone https://github.com/alexpreynolds/kmer-counter.git
$ cd kmer-counter
$ make
$ cp kmer-counter /usr/local/bin
</code></pre>

<p>Once the binary is in your path, you might use it in Python like so:</p>

<pre><code>k = 6
fastaFile = '/path/to/some/seqs.fa'
kmerCmd = 'kmer-counter --fasta --k=%d %s' % (k, fastaFile)
try:
    output = subprocess.check_output(kmerCmd, shell=True)
    result = {}
    for line in output.splitlines():
        (header, counts) = line.strip().split('\t')
        header = header[1:]
        kmers = dict((k,int(v)) for (k,v) in [d.split(':') for d in counts.split(' ')])
        result[header] = kmers
    sys.stdout.write("%s" % (str(result)))
except subprocess.CalledProcessError as error:
    sys.stderr.write("%s" % (str(error)))
</code></pre>

<p>Given example FASTA like this:</p>

<pre><code>&gt;foo
TTAACG
&gt;bar
GTGGAAGTTCTTAGGGCATGGCAAAGAGTCAGAATTTGAC
</code></pre>

<p>For k=6, you would get an iterable Python dictionary like this:</p>

<pre><code>{'foo': {'TTAACG': 1, 'CGTTAA': 1}, 'bar': {'GTTCTT': 1, 'AGAACT': 1, 'GAGTCA': 1, 'ATGGCA': 1, 'GAACTT': 1, 'ATTCTG': 1, 'CTAAGA': 1, 'CTTCCA': 1, 'ATTTGA': 1, 'GGAAGT': 1, 'AGGGCA': 1, 'CCTAAG': 1, 'CTCTTT': 1, 'AATTTG': 1, 'TCTGAC': 1, 'TTTGCC': 1, 'CTTAGG': 1, 'TTTGAC': 1, 'GAAGTT': 1, 'CCCTAA': 1, 'AGAATT': 1, 'AGTCAG': 1, 'CTGACT': 1, 'TCTTAG': 1, 'CGTTAA': 1, 'GTGGAA': 1, 'TGCCAT': 1, 'ACTCTT': 1, 'GGGCAT': 1, 'TTAGGG': 1, 'CTTTGC': 1, 'TGGAAG': 1, 'GACTCT': 1, 'CATGCC': 1, 'GCAAAG': 1, 'AAATTC': 1, 'GTCAAA': 1, 'TGACTC': 1, 'TAGGGC': 1, 'AAGTTC': 1, 'ATGCCC': 1, 'TCAAAT': 1, 'CAAAGA': 1, 'AACTTC': 1, 'GTCAGA': 1, 'CAAATT': 1, 'TAAGAA': 1, 'CATGGC': 1, 'AAGAAC': 1, 'AAGAGT': 1, 'TCTTTG': 1, 'TTCCAC': 1, 'TGGCAA': 1, 'GGCAAA': 1, 'AGTTCT': 1, 'AGAGTC': 1, 'TCAGAA': 1, 'GAATTT': 1, 'AAAGAG': 1, 'TGCCCT': 1, 'CCATGC': 1, 'GGCATG': 1, 'TTGCCA': 1, 'CAGAAT': 1, 'AATTCT': 1, 'GCATGG': 1, 'ACTTCC': 1, 'TTCTTA': 1, 'GCCATG': 1, 'GCCCTA': 1, 'TTCTGA': 1}}
</code></pre>

<p>You can use standard Python calls to manipulate this dictionary object and get sums of counts per record, for sequence, etc. which seems to answer your question. Please feel free to clarify what you're looking for if this object representation is not clear.</p>

<p>For a fully-fleshed out demonstration, see: <a href="https://github.com/alexpreynolds/kmer-counter/blob/master/test/kmer-test.py" rel="noreferrer">https://github.com/alexpreynolds/kmer-counter/blob/master/test/kmer-test.py</a></p>

CC BY-SA 3.0
2017-06-08T22:04:16.040Z
2017-07-27T20:08:44.917Z
2017-07-27T20:08:44.917Z
Alex Reynolds
<p>If you want to count the number of unique k-mers that occur in your data set, you should use the <a href="http://khmer.readthedocs.io/en/v2.1.1/user/scripts.html#unique-kmers-py" rel="noreferrer">unique-kmers.py</a> script, which implements a HyperLogLog-based cardinality estimator.</p>

<p>If the number of unique k-mers is not what you're after, please clarify.</p>

<p><strong>EDIT</strong>:</p>

<p>khmer uses probabilistic data structures internally, which store k-mers as hashed values that cannot be un-hashed back into k-mer sequences uniquely. Querying a k-mer's abundance requires you to know which k-mer(s) you're looking for, which in this case would require a second pass over the reads.</p>

<p>The closest thing to what you're asking provided in khmer's command-line scripts is <a href="http://khmer.readthedocs.io/en/v2.1.1/user/scripts.html#abundance-dist-py" rel="noreferrer">abundance-dist.py</a> (or the alternative abundance-dist-single.py), which will produce a k-mer abundance histogram but not per-kmer abundance.</p>

<p>The <a href="http://khmer.readthedocs.io/en/v2.1.1/user/scripts.html#load-into-counting-py" rel="noreferrer">load-into-counting.py</a> scripts computes the k-mer abundances and stores them in a probabilistic data structure, which is written to disk and can subsequently be re-loaded into memory quickly. Many of the scripts in khmer expect the reads to have been pre-processed by <code>load-into-counting.py</code>, while others will invoke the counting routines directly themselve. It's also fairly easy to do this with the Python API. In every case, the accuracy of the k-mer counts is dependant on <a href="http://khmer.readthedocs.io/en/v2.1.1/user/choosing-table-sizes.html" rel="noreferrer">memory considerations</a>.</p>

<pre><code>&gt;&gt;&gt; import khmer
&gt;&gt;&gt; counts = khmer.Counttable(31, 1e7, 4)
&gt;&gt;&gt; counts.consume_seqfile('reads.fq.gz')
(25000, 1150000)
&gt;&gt;&gt; counts.get('TGACTTTCTTCGCTTCCTGACGGCTTATGCC')
117
</code></pre>

<p>To check the false positive rate:</p>

<pre><code>&gt;&gt;&gt; fpr = khmer.calc_expected_collisions(counts)
&gt;&gt;&gt; fpr
</code></pre>

<p>If you <em>must</em> have the counts of each k-mer, it's not much more work...although reporting each k-mer's abundance only once would consume a lot of memory with a naive approach.</p>

<pre><code>&gt;&gt;&gt; outfile = open('outfile.txt', 'w')
&gt;&gt;&gt; seenkmers = set()  # Consumes a lot of memory for large input!!!
&gt;&gt;&gt; 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)
</code></pre>

CC BY-SA 3.0
2017-06-08T05:06:44.933Z
Not exactly what I'm looking for. EDITED.
SmallChess
CC BY-SA 3.0
2018-03-08T09:03:08.703Z
See this [question](https://bioinformatics.stackexchange.com/q/3794/48) if you run into problems not counting kmers above 255
llrs
CC BY-SA 3.0
2017-06-08T05:03:42.780Z
2017-06-08T06:04:19.577Z
2017-06-08T06:04:19.577Z
Daniel Standage
<p>I have some FASTQ sequence files and a FASTA file for some regions I'm interested in.</p>

<p>I would like:</p>

<ol>
<li>Build an index for the FASTA file</li>
<li>Use the index to count number of k-mers occurred in my sequence files</li>
</ol>

<p>I know how to do this in many k-mer counting tools, such as Kallisto, Jellyfish, DSK, Salmon. However, I'm more interested in Python (khmer is a Python package, all the other k-mer counting tools are C++).</p>

<p>Unfortunately, I can't figure out how to do that in khmer from the documentation:</p>

<blockquote>
  <p><a href="http://khmer.readthedocs.io/en/v2.1.1/introduction.html" rel="noreferrer">http://khmer.readthedocs.io/en/v2.1.1/introduction.html</a></p>
</blockquote>

<p>It talks about RNA-Seq, metagenomics, de Bruijn graph etc... All I need are just counts. Khmer is a popular k-mer tool, so it should be able to do it, but how?</p>

<p><strong>Q:</strong> How to build an index and count k-mers using khmer?</p>

<p><strong>EDIT:</strong></p>

<p>I don't need number of unique k-mers, I'm looking for a solution like Jellyfish where we build an index and count k-mers to the index.</p>

<p>Something like:</p>

<pre><code>MY_SEQ_1  100 counts
MY_SEQ_2  134 counts
</code></pre>

CC BY-SA 3.0
2017-07-27T17:04:35.153Z
Does this accomplish what you're looking for: https://github.com/dib-lab/khmer/blob/master/examples/python-api/exact-counting.py
Jon Deaton
<genome><k-mer><metagenome><khmer>
How to use Python to count k-mers?
# How to use Python to count k-mers?

(Asked by: username_0 on 2017-06-08)

I have some FASTQ sequence files and a FASTA file for some regions I'm interested in.
I would like:

Build an index for the FASTA file
Use the index to count number of k-mers occurred in my sequence files

I know how to do this in many k-mer counting tools, such as Kallisto, Jellyfish, DSK, Salmon. However, I'm more interested in Python (khmer is a Python package, all the other k-mer counting tools are C++).
Unfortunately, I can't figure out how to do that in khmer from the documentation:

http://khmer.readthedocs.io/en/v2.1.1/introduction.html

It talks about RNA-Seq, metagenomics, de Bruijn graph etc... All I need are just counts. Khmer is a popular k-mer tool, so it should be able to do it, but how?
Q: How to build an index and count k-mers using khmer?
EDIT:
I don't need number of unique k-mers, I'm looking for a solution like Jellyfish where we build an index and count k-mers to the index.
Something like:
MY_SEQ_1  100 counts
MY_SEQ_2  134 counts

**Tags:** genome, k-mer, metagenome, khmer
> **Comments:**
> **username_1**: Does this accomplish what you're looking for: https://github.com/dib-lab/khmer/blob/master/examples/python-api/exact-counting.py  

---

## 2 answers

### Answer by username_3 on 2017-06-08. Score: 5 | Accepted answer

If you want to count the number of unique k-mers that occur in your data set, you should use the unique-kmers.py script, which implements a HyperLogLog-based cardinality estimator.
If the number of unique k-mers is not what you're after, please clarify.
EDIT:
khmer uses probabilistic data structures internally, which store k-mers as hashed values that cannot be un-hashed back into k-mer sequences uniquely. Querying a k-mer's abundance requires you to know which k-mer(s) you're looking for, which in this case would require a second pass over the reads.
The closest thing to what you're asking provided in khmer's command-line scripts is abundance-dist.py (or the alternative abundance-dist-single.py), which will produce a k-mer abundance histogram but not per-kmer abundance.
The load-into-counting.py scripts computes the k-mer abundances and stores them in a probabilistic data structure, which is written to disk and can subsequently be re-loaded into memory quickly. Many of the scripts in khmer expect the reads to have been pre-processed by load-into-counting.py, while others will invoke the counting routines directly themselve. It's also fairly easy to do this with the Python API. In every case, the accuracy of the k-mer counts is dependant on memory considerations.
>>> import khmer
>>> counts = khmer.Counttable(31, 1e7, 4)
>>> counts.consume_seqfile('reads.fq.gz')
(25000, 1150000)
>>> counts.get('TGACTTTCTTCGCTTCCTGACGGCTTATGCC')
117

To check the false positive rate:
>>> fpr = khmer.calc_expected_collisions(counts)
>>> fpr

If you must have the counts of each k-mer, it's not much more work...although reporting each k-mer's abundance only once would consume a lot of memory with a naive approach.
>>> 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)

> **Comments:**
> **username_0**: Not exactly what I'm looking for. EDITED.  
> **username_4**: See this [question](https://bioinformatics.stackexchange.com/q/3794/48) if you run into problems not counting kmers above 255  

---

### Answer by username_2 on 2017-06-08. Score: 7

I wrote a command-line k-mer counter called kmer-counter that will output results in a form that your Python script can consume: https://github.com/alexpreynolds/kmer-counter
You can grab, build and install it like so:
$ git clone https://github.com/alexpreynolds/kmer-counter.git
$ cd kmer-counter
$ make
$ cp kmer-counter /usr/local/bin

Once the binary is in your path, you might use it in Python like so:
k = 6
fastaFile = '/path/to/some/seqs.fa'
kmerCmd = 'kmer-counter --fasta --k=%d %s' % (k, fastaFile)
try:
    output = subprocess.check_output(kmerCmd, shell=True)
    result = {}
    for line in output.splitlines():
        (header, counts) = line.strip().split('\t')
        header = header[1:]
        kmers = dict((k,int(v)) for (k,v) in [d.split(':') for d in counts.split(' ')])
        result[header] = kmers
    sys.stdout.write("%s" % (str(result)))
except subprocess.CalledProcessError as error:
    sys.stderr.write("%s" % (str(error)))

Given example FASTA like this:
>foo
TTAACG
>bar
GTGGAAGTTCTTAGGGCATGGCAAAGAGTCAGAATTTGAC

For k=6, you would get an iterable Python dictionary like this:
{'foo': {'TTAACG': 1, 'CGTTAA': 1}, 'bar': {'GTTCTT': 1, 'AGAACT': 1, 'GAGTCA': 1, 'ATGGCA': 1, 'GAACTT': 1, 'ATTCTG': 1, 'CTAAGA': 1, 'CTTCCA': 1, 'ATTTGA': 1, 'GGAAGT': 1, 'AGGGCA': 1, 'CCTAAG': 1, 'CTCTTT': 1, 'AATTTG': 1, 'TCTGAC': 1, 'TTTGCC': 1, 'CTTAGG': 1, 'TTTGAC': 1, 'GAAGTT': 1, 'CCCTAA': 1, 'AGAATT': 1, 'AGTCAG': 1, 'CTGACT': 1, 'TCTTAG': 1, 'CGTTAA': 1, 'GTGGAA': 1, 'TGCCAT': 1, 'ACTCTT': 1, 'GGGCAT': 1, 'TTAGGG': 1, 'CTTTGC': 1, 'TGGAAG': 1, 'GACTCT': 1, 'CATGCC': 1, 'GCAAAG': 1, 'AAATTC': 1, 'GTCAAA': 1, 'TGACTC': 1, 'TAGGGC': 1, 'AAGTTC': 1, 'ATGCCC': 1, 'TCAAAT': 1, 'CAAAGA': 1, 'AACTTC': 1, 'GTCAGA': 1, 'CAAATT': 1, 'TAAGAA': 1, 'CATGGC': 1, 'AAGAAC': 1, 'AAGAGT': 1, 'TCTTTG': 1, 'TTCCAC': 1, 'TGGCAA': 1, 'GGCAAA': 1, 'AGTTCT': 1, 'AGAGTC': 1, 'TCAGAA': 1, 'GAATTT': 1, 'AAAGAG': 1, 'TGCCCT': 1, 'CCATGC': 1, 'GGCATG': 1, 'TTGCCA': 1, 'CAGAAT': 1, 'AATTCT': 1, 'GCATGG': 1, 'ACTTCC': 1, 'TTCTTA': 1, 'GCCATG': 1, 'GCCCTA': 1, 'TTCTGA': 1}}

You can use standard Python calls to manipulate this dictionary object and get sums of counts per record, for sequence, etc. which seems to answer your question. Please feel free to clarify what you're looking for if this object representation is not clear.
For a fully-fleshed out demonstration, see: https://github.com/alexpreynolds/kmer-counter/blob/master/test/kmer-test.py


---

