How to use Python to count k-mers?

Question

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

I would like:

1. Build an index for the FASTA file
2. 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

Answer 1

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)

Answer 2

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