How to calculate the memory usage of storing kmers in RAM

Question

I want to write a program in C++ that stores kmers in a hash or in a trie. How can I calculate how much RAM I would need for each type of data structure? For this application my kmers are strand-specific, so I cannot reduce the size complexity by considering canonical kmers.

For a hash, for example, I know that there are $4^k$ possible kmers of length $k$. I also am assuming that DNA character with an alphabet $\in \{G,A,T,C\}$ can be stored in two bits ($2^2 = 4$). Using this logic the kmer hash should occupy $(4^k) \cdot 2 \cdot k$ bits. Plugging this into Wolfram Alpha ((4^13) * 2 * 13 bits to gigabytes), I calculated the amount of RAM that some small kmer sizes would require to store the entire hash in memory in GB.

k = 13 requires    0.2181 GB
k = 15 requires    4.027  GB
k = 17 requires   73.01  GB
k = 19 requires 1306   GB

This escalated rather quickly, and I think that my calculations must be wrong. Otherwise, how could programs count kmers in RAM?

The kmers that I want to store are sometimes sparse compared to the entire space (for k=13, only 29 million or so out of 67 million possibilities).

Constraints

• I do not need to store the number of occurences of a kmer, just its existence
• I cannot risk false positives. Bloom filters won't work in this case. :^/

Answer 1

To follow up on Devon Ryan's answer, I thought it would be a little fun to write a Python script that demonstrates using a bit array to maintain a presence/absence table.

Note: I wrote a C++ port that includes a custom bitset implementation that can be sized at runtime. This and the Python script are available on Github: https://github.com/alexpreynolds/kmer-boolean

#!/usr/bin/env python

import sys
import os
import bitarray

# read FASTA
def read_sequences():
    global seqs
    seqs = []
    seq = ""
    for line in sys.stdin:
        if line.startswith('>'):
            if len(seq) > 0:
                seqs.append(seq)
                seq = ""
        else:
            seq += line.strip()
    seqs.append(seq)

# build and initialize bit array
def initialize_bitarray():
    global ba 
    ba = bitarray.bitarray(4**k)
    ba.setall(False)
    sys.stderr.write("Memory usage of bitarray.bitarray instance is [%ld] bytes\n" % (ba.buffer_info()[1]))

# process sequences
def process_sequences():
    global observed_kmers
    observed_kmers = {}
    for seq in seqs:
        for i in range(0, len(seq)):
            kmer = seq[i:i+k]
            if len(kmer) == k:
                observed_kmers[kmer] = None
                idx = 0
                for j in range(k-1, -1, -1):
                    idx += 4**(k-j-1) * bm[kmer[j]]
                ba[idx] = True

def test_bitarray():
    test_idx = 0
    for j in range(k-1, -1, -1):
        test_idx += 4**(k-j-1) * bm[test_kmer[j]]
    test_result = ba[test_idx]
    if test_result:
        sys.stdout.write("%s found\n" % (test_kmer))
        sys.exit(os.EX_OK)
    else: 
        sys.stdout.write("%s not found\n" % (test_kmer))
        sys.exit(os.EX_DATAERR)

def main():
    global k
    k = int(sys.argv[1])

    global bm 
    bm = { 'A' : 0, 'C' : 1, 'T' : 2, 'G' : 3 }

    read_sequences()
    initialize_bitarray()
    process_sequences()

    try:
        global test_kmer
        test_kmer = sys.argv[2]
        if len(test_kmer) == k:
            test_bitarray()
        else:
            raise ValueError("test kmer (%s) should be of length k (%d)" % (test_kmer, k))
    except IndexError as err:
        keys = list(observed_kmers.keys())
        for i in range(0, len(keys)):
            sys.stdout.write("%s found\n" % (keys[i]))

    sys.exit(os.EX_OK)

if __name__== "__main__":
    main()

Note that this doesn't look at canonical kmers, e.g., AG is considered a distinct 2mer from its reverse complement CT.

To use this script, you pipe in your FASTA, specify the k, and an optional kmer that you want to test for presence/absence, e.g.:

$ echo -e ">foo\nCATTCTC\nGGGAC\n>bar\nTTATAT\n>baz\nTTTATTAG\nACCTCT" | ./kmer-bool.py 2 CG
Memory usage of bitarray.bitarray instance is [2] bytes
CG found

Or:

$ echo -e ">foo\nCATTCTC\nGGGAC\n>bar\nTTATAT\n>baz\nTTTATTAG\nACCTCT" | ./kmer-bool.py 3 AAA
Memory usage of bitarray.bitarray instance is [8] bytes
AAA not found

Or if the optional test kmer is left out:

$ echo -e ">foo\nCATTCTC\nGGGAC\n>bar\nTTATAT\n>baz\nTTTATTAG\nACCTCT" | ./kmer-bool.py 5
Memory usage of bitarray.bitarray instance is [128] bytes
CATTC found
ATTCT found
TTCTC found
TCTCG found
CTCGG found
TCGGG found
CGGGA found
GGGAC found
TTATA found
TATAT found
TTTAT found
TTATT found
TATTA found
ATTAG found
TTAGA found
TAGAC found
AGACC found
GACCT found
ACCTC found
CCTCT found

Or for the ~67M kmers in a 13-mer set, for which a roughly 8.4MB bit array is reserved:

$ echo -e ">foo\nCATTCTC\nGGGAC\n>bar\nTTATAT\n>baz\nTTTATTAG\nACCTCT" | ./kmer-bool.py 13
Memory usage of bitarray.bitarray instance is [8388608] bytes
TTTATTAGACCTC found
TTATTAGACCTCT found

Answer 2

I think that my calculations must be wrong. Otherwise, how could programs count kmers in RAM?

Hash table based k-mer counters only keep k-mers seen in data. For $16<k\le32$, you need 64-bit integers to store a k-mer. Given $n$ distinct k-mers, a naive implementation with open addressing hash table would roughly take $2n\cdot 64/8=16n$ bytes. We assume the table is half full – this is where factor 2 comes from. There are various tricks to reduce the memory footprint. You'd better read papers. On a side note, don't use std::unordered_map. It wastes extra memory.

A bitmap described by other answers takes $4^k/8=2^{2k-3}$ bytes, which is a better choice sometimes, but rarely used by k-mer counters as it does not work with long k, either. Trie or full-text indices are overkilling for fixed k-mers and are usually slower.

By the way, this question reminds me of your question on efficient reverse complement. Note that for k-mer counting, you traverse both strands at the same time. You should not reverse complement the whole k-mer.

Answer 3

If you simply need a set of present/absent k-mers and your kmer length is short, then convert the sequence length to an integer (a 13-mer would require 26 bits or 4 bytes) and use that as an index into a bit array. Initialize the array to 0 and use 1 for "k-mer present". If you have ~67 million k-mers then your array is 67 million bits, which is only about 8.5 megs.

This will scale "OK", but eventually balloon to a huge size if you use larger k-mers. At that point you simply have to use other data structures. The reality is that most applications are fairly tolerant to small amounts of error.