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'.