0
$\begingroup$

I implemented the following find_neighbors_with_expected_hamming_distance function that generates all k-mers of Hamming distance at most d from the given Pattern.

from .Hamming_Distance import find_hamming_distance

dp_map_that_stores_all_neighbors_with_expected_hamming_distance = {}

unique_nucleotides = {'A', 'C', 'G', 'T'}


def suffix(pattern):
    if len(pattern) <= 1:
        return None
    else:
        return pattern[1: (len(pattern))]


def first_symbol(pattern):
    if len(pattern) == 0:
        return None
    else:
        return pattern[0]


def find_neighbors_with_expected_hamming_distance(pattern, expected_hamming_distance):
    if expected_hamming_distance == 0:
        return {pattern}

    if len(pattern) == 1:
        return {'A', 'C', 'G', 'T'}

    if (pattern, expected_hamming_distance) in dp_map_that_stores_all_neighbors_with_expected_hamming_distance.keys():
        return dp_map_that_stores_all_neighbors_with_expected_hamming_distance[(pattern, expected_hamming_distance)]

    neighborhood_set = set()

    suffix_neighbors_set = find_neighbors_with_expected_hamming_distance(suffix(pattern), expected_hamming_distance)

    for suffix_neighbor in suffix_neighbors_set:

        if find_hamming_distance(suffix(pattern), suffix_neighbor) < expected_hamming_distance:

            for nucleotide in unique_nucleotides:
                neighborhood_set.add(nucleotide + suffix_neighbor)

        else:
            neighborhood_set.add(first_symbol(pattern) + suffix_neighbor)

    if (pattern, expected_hamming_distance) not in dp_map_that_stores_all_neighbors_with_expected_hamming_distance.keys():
        dp_map_that_stores_all_neighbors_with_expected_hamming_distance[
            (pattern, expected_hamming_distance)] = neighborhood_set

    return neighborhood_set
    

I want your suggestion on my time complexity calculation:

For a given pattern of length N and D as maximum expected hamming distance, I calculated time complexity as -

$\sum_{n=1}^N$ $\sum_{d=1}^D$ $n\choose d$ * $3^i$

Reasoning: For any length n, we will choose any d letters (nucleotide in this case) to be replaced with the other 3 nucleotides that is different from the nucleotide present at that particular position.

For example: Pattern = ATGCAT and D = 2, then suppose for one of the case, I select 2nd and 4th positions to be replaced.

Then, at the second position (current nucleotide = 'T'), I can substitute any nucleotide from the set {'A', 'C', 'G'} and

for fourth position (current nucleotide = 'C'), I can substitute any nucleotide from the set {'A', 'T', 'G'}.

Kindly tell me if this approach is correct?

EDIT:

I have provided the pseudocode for the Neighbors function as follows:

Neighbors(Pattern, d)
    if d = 0
        return {Pattern}
    if |Pattern| = 1 
        return {A, C, G, T}
    Neighborhood ← an empty set
    SuffixNeighbors ← Neighbors(Suffix(Pattern), d)
    for each string Text from SuffixNeighbors
        if HammingDistance(Suffix(Pattern), Text) < d
            for each nucleotide x
                add x • Text to Neighborhood
        else
            add FirstSymbol(Pattern) • Text to Neighborhood
    return Neighborhood

The symbol • means concatenation of 2 strings. For example, 'S' • 'V' = 'SV'.

If we remove the first symbol of Pattern (denoted FirstSymbol(Pattern)), then we will obtain a

(k − 1)-mer that we denote by Suffix(Pattern).

Illustration:

Given pattern = 'ACG'

Expected Hamming distance (D) = 1

Answer = {'CCG', 'TCG', 'GCG', 'AAG', 'ATG', 'AGG', 'ACA', 'ACC', 'ACT', 'ACG'}

Note that the Answer is a set of all strings who have hamming distance of 1 with respect to the given pattern 'ACG'.

$\endgroup$
3
  • $\begingroup$ I'm confused by the "expectation" part. Isn't this deterministic? $\endgroup$ Jul 1 at 4:02
  • $\begingroup$ Also, I'd recommend adding a pseudocode explanation of your approach in addition to the python code. $\endgroup$ Jul 1 at 4:10
  • $\begingroup$ @Throckmorton, thank you very much for having a look. I have added the pseudocode with an illustration. $\endgroup$ Jul 1 at 12:54
2
$\begingroup$

Let's first establish a lower bound on any solution:

Given a hamming distance $D$ and an input pattern $P$ of length $N$, there are $\sum_{k=0}^D{N \choose k} 3^k$ other strings within hamming distance $D$ of $P$. Of course in the worst case, when $D=N$, there are $4^n$ possible strings. We can also get a closed bound on the number of solutions. We can show that the sum of the first $D$ binomial coefficients for a fixed $N$ is bounded by ${N \choose D} {N-(D-1) \over N-(2D-1)}$ as shown here. Therefore we have the following closed upper bound on the number of strings within hamming distance $D$ of some pattern of length $N$

$$ \sum_{k=0}^D{N \choose k} 3^k \leq {N \choose D} {3^D (N-(D-1)) \over N-(2D-1)} \leq \left(\frac{3N}{D}\right)^D {N-(D-1) \over N-(2D-1)} $$

However, for your approach, you are going through recursively, so when the subproblem is on strings within $D$ of the suffix of length $n$, you have $\sum_{k=0}^D{n \choose k} 3^k$ of these. Given that you have $N$ subproblems (suffixes), at first glance it does indeed look like you have roughly $\sum_{n=1}^N\sum_{k=0}^D{n \choose k} 3^k$ iterations of the innermost for loop. However, in that for loop, you are also computing the Hamming Distance function, which is $O(n)$. This means that you have roughly

$$ \sum_{n=1}^Nn\sum_{k=0}^D{n \choose k} 3^k $$

computations.

I would suggest finding a way to get rid of having to compute the hamming function. Since you are constructing the sequences yourself, you can keep track of each sequences hamming distance from the pattern during construction.

$\endgroup$
2
  • $\begingroup$ Excellent! I missed the complexity for computing Hamming distance. Also, you can correct the first sentence - "Given a hamming distance D" $\endgroup$ Jul 1 at 17:25
  • $\begingroup$ Thank you very much @Throckmorton. $\endgroup$ Jul 1 at 17:27

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.