汉明德斯坦逆

Question

*这是一个简短的介绍，具体问题在最后一段用粗体字表示。在

我正在尝试生成所有具有给定Hamming距离的字符串，以有效地解决生物信息分配问题。在

我们的想法是，给定一个字符串（即'acgttgcatcgcatgagct'）、要搜索的单词的长度（即4）以及在字符串中搜索该单词时可接受的不匹配（即1），返回最频繁的单词或“变种”单词。在

明确地说，给定字符串中长度为4的单词可以是（在“[]”之间）：

[ACGT]TGCATGTCGCATGATGCATGAGAGCT #ACGT

这个

^{pr2}$

或者这个

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)

我想做的是，不是生成所有可能的单词，而是只生成出现在给定字符串中的单词的突变，具有给定数量的不匹配，换句话说，给定一个Hamming距离和一个单词，返回具有该距离（或更少）的所有可能的变异单词，然后使用它们进行搜索给定的字符串。在

我希望我是清楚的。谢谢您。在

Answer 1

设给定的汉明距离为D，设w为“word”子串。从w，您可以通过深度限制depth-first search生成距离≤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))

（这绝不会很快，但可能会给人以启发。）

Answer 2

如果我正确地理解了你的问题，你想确定基因组中得分最高的k-mers。k-mer的得分是它在基因组中出现的次数加上汉明距离小于m的任何k-mer在基因组中出现的次数。请注意，这假设您只对出现在您的基因组中的k-mer感兴趣（正如@j_random_hacker所指出的）。在

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

1. 确定基因组中的所有k-mer。在
2. 计算每个k-mer在G中出现的次数。在
3. 对于每对k-mers（K1，K2），如果K1和{}的汉明距离小于m，则增加它们的计数。在
4. 找到maxk-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)

Answer 3

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)

此函数生成汉明距离小于或等于给定数字的单词的所有突变。它是相对有效的，并且不检查无效的单词。然而，有效突变可以出现多次。如果希望每个元素都是唯一的，请使用集合。在