Dfyjkt

I have a large array with millions of DNA sequences which are all 24 characters long. The DNA sequences should be random and can only contain A,T,G,C,N. I am trying to find strings that are within a certain hamming distance of each other.

My first approach was calculating the hamming distance between every string but this would take way to long.

My second approach used a masking method to create all possible variations of the strings and store them in a dictionary and then check if this variation was found more then 1 time. This worked pretty fast(20 min) for a hamming distance of 1 but is very memory intensive and would not be viable to use for a hamming distance of 2 or 3.

Python 2.7 implementation of my second approach.

sequences = 
masks = 
for sequence in sequences:
 for i in range(len(sequence)):
 try:
 masks[sequence[:i] + '?' + sequence[i + 1:]].append(sequence[i])
 except KeyError:
 masks[sequence[:i] + '?' + sequence[i + 1:]] = [sequence[i], ]

matches = 
for mask in masks:
 if len(masks[mask]) > 1:
 matches[mask] = masks[mask]

I am looking for a more efficient method. I came across Trie-trees, KD-trees, n-grams and indexing but I am lost as to what will be the best approach to this problem.

One approach is Locality Sensitive Hashing

First, you should note that this method does not necessarily return all the pairs, it returns all the pairs with a high probability (or most pairs).

Locality Sensitive Hashing can be summarised as: data points that are located close to each other are mapped to similar hashes (in the same bucket with a high probability). Check this link for more details.

Your problem can be recast mathematically as:

Given N vectors v ∈ R^24, N<<5^24 and a maximum hamming distance d, return pairs which have a hamming distance atmost d.

The way you'll solve this is to randomly generates K planes P_1,P_2,...,P_K in R^24; Where K is a parameter you'll have to experiment with. For every data point v, you'll define a hash of v as the tuple Hash(v)=(a_1,a_2,...,a_K) where a_i∈0,1 denotes if v is above this plane or below it. You can prove (I'll omit the proof) that if the hamming distance between two vectors is small then the probability that their hash is close is high.

So, for any given data point, rather than checking all the datapoints in the sequences, you only check data points in the bin of "close" hashes.

Note that these are very heuristic based and will need you to experiment with K and how "close" you want to search from each hash. As K increases, your number of bins increase exponentially with it, but the likelihood of similarity increases.

Judging by what you said, it looks like you have a gigantic dataset so I thought I would throw this for you to consider.

Found my solution here: http://www.cs.princeton.edu/~rs/strings/

This uses ternary search trees and took only a couple of minutes and ~1GB of ram. I modified the demo.c file to work for my use case.