22

*这是一个简单的介绍，具体问题在最后一段以粗体显示。

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

这个想法是，给定一个字符串（即'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)

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

我希望我很清楚。谢谢你。

4

6 回答 6

18 
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)

This function generates all mutations of a word with a hamming distance less than or equal to a given number. It is relatively efficient, and does not check invalid words. However, valid mutations can appear more than once. Use a set if you want every element to be unique.

于 2013-11-12T22:54:31.137 回答
5

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

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

  1. 识别基因组中的所有 k-mer G。
  2. 计算每个 k-mer 在 中出现的次数G。
  3. 对于每对 ( K1, ) k-mers ，如果它们的汉明距离小于，则K2增加两者的计数。K1K2m
  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)
于 2013-11-12T22:38:24.607 回答
5

让给定的汉明距离为D，让w为“单词”子串。从w ，您可以通过深度限制的深度优先搜索生成距离 ≤<em>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))

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

于 2013-11-12T22:32:42.947 回答
2

这就是我认为您要解决的问题是：您有一个长度为 n 的“基因组”，并且您想找到在该基因组中近似出现最频繁的 k-mer，其中“近似出现”意味着出现汉明距离 <= d。这个k-mer实际上不需要精确地出现在基因组的任何地方（例如对于基因组ACCA，k=3和d=1，最好的k-mer是CCC，出现两次）。

如果您从字符串中的某个 k-mer 生成汉明距离 <= d 的所有 k-mers，然后搜索字符串中的每个 k-mers，就像您目前正在做的那样，那么您正在添加一个不必要的 O(n)搜索时间的因素（除非您使用Aho-Corasick 算法同时搜索它们，但这在这里有点过分了）。

您可以通过遍历基因组做得更好，在每个位置 i，从基因组中的位置 i 开始，生成与 k-mer 距离 <= d 的所有 k-mer 的集合，并为每个位置增加一个计数器一。

于 2013-11-12T23:19:25.563 回答
1 
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:]

将生成 dist 为 1 的排列（它也会重复包含重复）

得到一组在 2 内的所有...然后用第一个解决方案中的每一个作为初始猜测调用它

于 2013-11-12T22:31:51.987 回答
1

这里有正确的答案，它大量利用了 python 的神奇功能，几乎可以为你做所有事情。我将尝试用数学和算法来解释事物，以便您可以将其应用于您想要的任何语言。

因此，在您的情况下，您有一个字母表{a1, a2, ... a_a}（的基数a），{'A', 'C', 'G', 'T'}基数是 4。您有一个长度字符串，k并且您想要生成汉明距离小于或等于 的所有字符串d。

首先你有几个？答案不取决于您选择的字符串。如果您选择了一个弦，您将拥有与您的弦C(d, k)(a-1)^d有汉明距离d的弦。所以字符串的总数是：

在此处输入图像描述

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

那么你将如何推导出一个生成所有字符串的算法呢？请注意，很容易生成一个距您的 wold 最多一个汉明距离的字符串。您只需要遍历字符串中的所有字符，并为每个字符尝试字母表中的每个字母。如您所见，有些词是相同的。

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

所以这是一个伪代码：

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
}
于 2015-12-30T05:18:13.007 回答