汉明距离的逆运算

Question

我正在尝试生成所有与给定的汉明距离相同的字符串，以高效地解决生物信息学作业。

具体思路是，给定一个字符串（例如“ACGTTGCATGTCGCATGATGCATGAGAGCT”），要搜索的单词长度（例如4）以及在搜索该单词时允许的不匹配数（例如1），返回最常见的单词或“突变”单词。

明确一下，在给定字符串中，长度为4的单词可以是以下格式之一（用方括号'[]'表示）：

请翻译以上内容

[ACGT]TGCATGTCGCATGATGCATGAGAGCT #ACGT

这个

A[CGTT]GCATGTCGCATGATGCATGAGAGCT #CGTT

或者这个

ACGTTGCATGTCGCATGATGCATGAG[AGCT] #AGCT

我所做的事情是（这个过程非常低效，当需要生成10个字符的单词时速度会变得很慢），生成所有与给定距离相同的可能单词：

itertools.imap(''.join, itertools.product('ATCG', repeat=wordSize))

然后搜索和比较给定字符串中的每个单词，如果生成的单词（或其变异体）出现在循环中：

wordFromString = givenString[i:i+wordSize]
mismatches = sum(ch1 != ch2 for ch1, ch2 in zip(wordFromString, generatedWord))
if mismatches <= d:
    #count that generated word in a list for future use
    #(only need the most repeated)

我想做的是，不是生成所有可能的单词，而是仅生成给定字符串中出现的单词的变异体，其具有给定数量的不匹配项，换句话说，给定一个海明距离和一个单词，返回所有可能的变异单词以及这个距离（或更少），然后将它们用于在给定的字符串中搜索。

希望我表述清楚了。谢谢。

- JackS

---

你会计算重叠吗？例如，对于字符串“CCACCAC”和单词“CCAC”，你会计算两个出现还是只有一个？因为“C”重叠了。 - גלעד ברקן

不可知的答案及其解释 https://dev59.com/f2Ij5IYBdhLWcg3wuXeN#34523345 - Salvador Dali

Answer 1

def mutations(word, hamming_distance, charset='ATCG'):
    for indices in itertools.combinations(range(len(word)), hamming_distance):
        for replacements in itertools.product(charset, repeat=hamming_distance):
            mutation = list(word)
            for index, replacement in zip(indices, replacements):
                mutation[index] = replacement
            yield "".join(mutation)

这个函数生成一个单词的所有汉明距离小于或等于给定数字的变异。它相对高效，但不检查无效的单词。然而，有效的变异可能会出现多次。如果您想要每个元素都是唯一的，请使用集合(set)。

- Joshua

Answer 2

假设给定的汉明距离为D，并且让w成为“单词”子字符串。从w开始，您可以通过深度限制的深度优先搜索生成所有距离≤D的单词：

def dfs(w, D):
    seen = set([w])      # keep track of what we already generated
    stack = [(w, 0)]

    while stack:
        x, d = stack.pop()
        yield x

        if d == D:
            continue

        # generate all mutations at Hamming distance 1
        for i in xrange(len(x)):
            for c in set("ACGT") - set(x[i])
                 y = x[:i] + c + x[i+1:]
                 if y not in seen:
                     seen.add(y)
                     stack.append((y, d + 1))

（这绝不会很快，但它可能作为灵感。）

- Fred Foo

Answer 3

如果我正确理解您的问题，您想要在基因组G中识别得分最高的k-mer。一个k-mer的得分是它在基因组中出现的次数加上任何Hamming距离小于m的k-mer也在基因组中出现的次数。请注意，这假定您只对出现在您的基因组中的k-mer感兴趣（如@j_random_hacker所指出的）。

您可以通过以下四个基本步骤来解决此问题：

1.识别基因组G中的所有k-mer。 2.计算每个k-mer在G中出现的次数。 3.对于每一对（K1，K2）k-mer，如果它们的Hamming距离小于m，则同时增加K1和K2的计数。 4.找到max k-mer及其计数。

以下是Python示例代码：

from itertools import combinations
from collections import Counter

# Hamming distance function
def hamming_distance(s, t):
    if len(s) != len(t):
        raise ValueError("Hamming distance is undefined for sequences of different lengths.")
    return sum( s[i] != t[i] for i in range(len(s)) )

# Main function
# - k is for k-mer
# - m is max hamming distance
def hamming_kmer(genome, k, m):
    # Enumerate all k-mers
    kmers = [ genome[i:i+k] for i in range(len(genome)-k + 1) ]

    # Compute initial counts
    initcount  = Counter(kmers)
    kmer2count = dict(initcount)

    # Compare all pairs of k-mers
    for kmer1, kmer2 in combinations(set(kmers), 2):
        # Check if the hamming distance is small enough
        if hamming_distance(kmer1, kmer2) <= m:
            # Increase the count by the number of times the other
            # k-mer appeared originally
            kmer2count[kmer1] += initcount[kmer2]
            kmer2count[kmer2] += initcount[kmer1]

    return kmer2count


# Count the occurrences of each mismatched k-mer
genome = 'ACGTTGCATGTCGCATGATGCATGAGAGCT'
kmer2count = hamming_kmer(genome, 4, 1)

# Print the max k-mer and its count
print max(kmer2count.items(), key=lambda (k,v ): v )
# Output => ('ATGC', 5)

- mdml

---

最佳的k-mer不一定实际出现在基因组中。例如，假设基因组为ACCA，k=3，距离阈值d=1：那么3-mer CCC 出现了两次（一次从第一个位置开始，一次从第二个位置开始，每个的汉明距离为1），而在基因组中出现的两个3-mer (ACC 和 CCA) 每个只出现一次。 - j_random_hacker

@j_random_hacker：好观点，我更新了我的答案以注意到那个警告。 - mdml

这恐怕还有其他严重的问题。如果某个k-mer出现m次，则在自我比较中它将被计算m *（m-1）次 - 因此例如，如果基因组为AAAAA而k = 2，则您的代码将显示该2-merAA出现了12次！ - j_random_hacker

实际上它将被计算2m(m-1)次，这仍然是错误的答案--在上面的例子中为24。 - j_random_hacker

@j_random_hacker：更多好点子，谢谢建议！我更新了我的回答，使用了k-mer的“set”并存储了初始计数。 - mdml

不客气 :) 对于我们只关心在基因组中恰好出现一次的k-mer版本的问题，这看起来现在是正确的。它还具有一个很好的特点，即对于长度为n且具有r个不同k-mer的基因组，它仅需要O(r^2)时间运行（一个朴素算法仍需要O(n^2)时间）。 - j_random_hacker

Answer 4

这是我认为你试图解决的问题：你有一个长度为n的“基因组”，想要找到在该基因组中近似出现最频繁的k-mer，其中“近似出现”意味着与Hamming距离<= d。这个k-mer不必实际上在基因组中出现（例如对于基因组ACCA，k=3和d=1，最好的k-mer是CCC，出现两次）。

如果你从字符串中的某个k-mer生成所有Hamming距离<= d的k-mers，然后搜索每一个，在搜索时间上就会添加一个不必要的O(n)因子（除非你使用 Aho-Corasick算法同时搜索所有k-mers，但这在这里是过度的）。

你可以通过遍历基因组，在每个位置i处生成所有与基因组中起始于位置i的k-mer距离<= d的k-mer集合，并对每个k-mer增加计数器来实现更好的效果。

- j_random_hacker

Answer 5

这里有一些正确的答案，它们大量使用了 Python 的神奇函数，几乎为您做了所有事情。我将尝试用数学和算法来解释问题，以便您可以将其应用于任何您想要的语言。

---

所以你有一个字母表{a1，a2，... a_a}（a的基数），在这种情况下为{'A'，'C'，'G'，'T'}，基数为4。您有一个长度为k的字符串，并希望生成所有汉明距离小于或等于d的字符串。

首先，您有多少个字符串？答案不取决于您选择的字符串。如果您选择了一个字符串，则会有C(d，k)（a-1）^ d个字符串与您的字符串的汉明距离d相同。因此，字符串的总数是：

它在几乎每个参数方面呈指数增长，因此您将无法找到所有单词的任何快速算法。

---

那么，你如何推导出一个可以生成所有字符串的算法呢？注意，生成一个距离您的单词最多只有一个汉明距离是很容易的。您只需要迭代字符串中的所有字符，并为每个字符尝试字母表中的每个字母。正如您将看到的那样，其中一些单词将是相同的。

现在要生成所有与您的字符串相隔两个汉明距离的字符串，您可以对上一次迭代中的每个单词应用生成一个汉明距离单词的相同函数。

因此，这里是伪代码：

function generateAllHamming(str string, distance int): 
    alphabet = ['A', ...]// array of letters in your alphabet
    results = {} // empty set that will store all the results
    prev_strings = {str} // set that have strings from the previous iterations
    // sets are used to get rid of duplicates

    if distance > len(str){ distance = len(str)} // you will not get any new strings if the distance is bigger than length of the string. It will be the same all possible strings.

    for d = 0; d < distance; d++ {
        for one_string in prev_strings {
            for pos = 0; pos < len(one_string); pos++ {
                for char in alphabet {
                    new_str = substitute_char_at_pos(one_string, char, pos)
                    add new_str to set results 
                }
            }
        }

        populate prev_strings with elements from results
    }

    return your results
}

- Salvador Dali

Answer 6

def generate_permutations_close_to(initial = "GACT",charset="GACT"):
    for i,c in enumerate(initial):
         for x in charset:
             yield initial[:i] + x + inital[i+1:]

这将生成距离为1的排列（它也可能包含重复）

要获得所有距离为2以内的集合...然后使用每个第一个解作为初始猜测来调用此函数。

- Joran Beasley