Are there any rolling hash functions that can hash a DNA sequence and its reverse complement to the same value?

Question

A common bioinformatics task is to decompose a DNA sequence into its constituent k-mers and compute a hash value for each k-mer. Rolling hash functions are an appealing solution for this task, since they can be computed very quickly. A rolling hash does not compute each the hash value from scratch with each k-mer: rather it updates a running hash value using an update strategy and a sliding window over the data.

It's also very useful for many applications to have a k-mer hash to the same value as its reverse complement. Unless the data were generated using a strand-specific sample prep, it's impossible to distinguish a k-mer from its reverse complement, and they should be treated as the same sequence.

Are there any rolling hashes that will map reverse complements to the same value? If not, how would we develop such an algorithm?

UPDATE: Ideally the hash function would be able to support k > 32, which would be lossy unless using something larger than a 64-bit integer.

ANOTHER UPDATE: I don't think it's necessary to store both the running k-mer and its reverse complement in a single value. If storing two k-mer strings and/or two hash values makes this easier, I'm totally cool with that.

Answer 1

Let f and r be two integers. They always keep the k-mer on the forward and reverse strand, respectively. At a new base c in a proper 2-bit encoding, we update the two integers as follows:

f = (f<<2|c) & ((1ULL<<2*k) - 1)
r = r>>2 | (3ULL-c)<<2*(k-1)

With this updating rule, f keeps the forward strand k-mer ending at c and r is f's reverse complement. The hash of the k-mer can be min(f,r). If integer operations take constant time, computing the hashes of all k-mers of a sequence of length L takes O(L) time.

The above only works if k-mer fits a word, which means k<=32 on x86. For long k-mers, you can overload bit operators in C++. In case of 32<k<=64, a simpler solution is to split into lower and higher bits:

c1 = c&1, c2 = c>>1
f1 = (f1<<1|c1) & ((1ULL<<k) - 1)
f2 = (f2<<1|c2) & ((1ULL<<k) - 1)
r1 = r1>>1 | (1ULL-c1)<<(k-1)
r2 = r2>>1 | (1ULL-c2)<<(k-1)

This is twice as slow as the k<=32 version. It is possible to implement these lines with SSE2 intrinsics, but that is overkilling.

When you want to hash a long k-mer into fewer bits, there are a few options:

xor:      f^r
min:      min(f,r)
min+hash: hash(min(f,r))

where hash() is a generic randomization hash function, which could be Murmur3, Wang's integer hash function, etc.

The following is a microbenchmark. Here, we squeeze all k-mers into a 32-bit hash table. collision_rate = 1 - #distinct_hash/#distinct_kmer. This is not the best metric, but it somehow measures the randomness, which is better than nothing.

====================================================
 Algorithm    data  max-k  k   %collision  CPU time
----------------------------------------------------
 xor          chr11  32    31     8.5%       0.9s
 xor+Wang     chr11  32    31     9.6%       1.3s
 min+Wang     chr11  32    31     1.4%       1.3s
 min+Wang     chr11  64    31     1.4%       1.9s
 min+murmur3  chr11  inf   31     1.4%       5.7s
 min+Wang     chr11  64    51     1.5%       2.1s
 min+murmur3  chr11  inf   51     1.5%       6.8s
 min+Wang     chr11  64    63     1.5%       2.1s
 min+murmur3  chr11  inf   63     1.5%       7.5s
 radix sort   chr11  -     -      -        5.6-7.3s
----------------------------------------------------
 min+Wang     hg38   64    63    26.2%        37s
 min+murmur3  hg38   inf   63    26.2%       192s
 radix sort   hg38   -     -      -         ~138s
====================================================

Observations:

• Don't use XOR. In fact, XOR hashes all palindromes to 0xffff.... This is already worrying enough.

• Unless you want to work with very long k-mers, don't use generic string hash functions like FNV and Murmur. Murmur is even slower than radix sorting all hashes.

• For k<=64, min+Wang is the best here. It is fast and simple to compute and has randomness nearly as good as murmur.

Answer 2

A rolling hash function for DNA sequences called ntHash has recently been published in Bioinformatics and the authors dealt with reverse complements:

Using this table, we can easily compute the hash value for the reverse-complement (as well as the canonical form) of a sequence efficiently, without actually reverse- complementing the input sequence, as follows: ...

EDIT (by @user172818): I will add more details about how ntHash works. The notations used in its paper are somewhat uncommon. The source code is more informative.

Let's first define rotation functions for 64-bit integers:

rol(x,k) := x << k | x >> (64-k)
ror(x,k) := x >> k | x << (64-k)

We then define a hash function h() for each base. In the implementation, the authors are using:

h(A) = 0x3c8bfbb395c60474
h(C) = 0x3193c18562a02b4c
h(G) = 0x20323ed082572324
h(T) = 0x295549f54be24456
h(N) = 0

The rolling hash function of a forward k-mer s[i,i+k-1] is:

f(s[i,i+k-1]) := rol(h(s[i]),k-1) ^ rol(h(s[i+1]),k-2) ^ ... ^ h(s[i+k-1])

where ^ is the XOR operator. The hash function of its reverse complement is:

r(s[i,i+k-1]) := f(~s[i,i+k-1])
               = rol(h(~s[i+k-1]),k-1) ^ rol(h(~s[i+k-2]),k-2) ^ ... ^ h(~s[i])

where ~ gives the reverse complement of a DNA sequence. Knowing f(s[i,i+k-1]) and r(s[i,i+k-1]), we can compute their values for the next k-mer:

f(s[i+1,i+k]) = rol(f(s[i,i+k-1]),1) ^ rol(h(s[i]),k)  ^ h(s[i+k])
r(s[i+1,i+k]) = ror(r(s[i,i+k-1]),1) ^ ror(h(~s[i]),1) ^ rol(h(~s[i+k]),k-1)

This works because rol, ror and ^ can all be switched in order. Finally, for a k-mer s, the hash function considering both strands is the smaller between f(s) and r(s):

h(s) = min(f(s),r(s))

This is a linear algorithm regardless of the k-mer length. It only uses simple arithmetic operations, so should be fairly fast. I have briefly tested its randomness. It seems comparable to murmur. ntHash is probably the best algorithm so far if you want to hash an arbitrarily long k-mer into 64 bits.

Answer 3

If your goal is to minimise storage by just having one hash per kmer and its reverse complement, there is a simple solution for non-rolling hashes. For any sequence S, you compute and store the hash of the smaller of S and its reverse complement. A simple lexicographical comparison for "smaller" is enough. In programming terms:

hash=computeHash(min(S,rev(S));

If your goal is to minimise computing time via rolling hashes AND store only one kmer value, this will be hard. Hashes are normally designed to minimise collisions and you are asking it not only to collide, but also to collide in very, very specific circumstances.

Question is: is the computation of a hash a bottleneck for an application? Maybe, maybe not.

This post on StackExchange seems to imply that, in 2012, MurmurHash was able to compute de-novo a hash for a UUID (similar to a standard kmer in size and complexity) in ~250ns, i.e., ~4 million hashes per second. That's already quite fast.

My guess is that most applications will use a hash to look up things in memory or disk. And here, even with look-ups in memory, real life data will let you run into CPU cache miss speed penalties very quickly and I expect the impact of this to be higher than the speed of the hash computation itself.

Answer 4

At its easiest, you just store the forward ($F$) and reverse ($R$) hash value.

You update the forward hash value by conventional means, e.g. bit-shifting the base value into its lower bits:

$$ F_{n + 1} = ((F_n \ll B \mid x) \mathop\& M, $$

$B$ is the bit size of the base encoding, $M$ is the word mask for a word of length $W$ bits, and can be omitted if the width of the hash is the width of the machine word; and the reverse hash value by doing the mathematical inverse, e.g. bit-shifting the value of the reverse complement of the base into the upper bits:

$$ R_{n + 1} = ((R_n \gg B) \mid (\operatorname{compl}x \ll (W - B))) \mathop\& M. $$

And then you compute a combined hash value by xoring the two hashes, $H = F \oplus R$.