查找 DNA 串中最常见的错配 k-mer

Question

我正在尝试解决在字符串中查找最常见的不匹配 k-mers 的问题。要求如下：

Frequent Words with Mismatches 问题：找出字符串中出现不匹配的最频繁的 k-mers。

输入：字符串 Text 以及整数 k 和 d。 （您可以假设 k ≤ 12 和 d ≤ 3。）

输出：所有最常见的 k-mer，文本中最多有 d 个不匹配。

这是一个例子：

示例输入：

ACGTTGCATGTCGCATGATGCATGAGAGCT

4 1

样本输出：

GATG ATGC ATGT

最简单但效率最低的方法是列出文本中的所有k-mers并计算它们之间的hamming_difference，并挑选出hamming_difference小于或等于d的模式，下面是我的代码：

kmer = 4
in_genome = "ACGTTGCATGTCGCATGATGCATGAGAGCT";
in_mistake = 1;
out_result = [];
mismatch_list = []

def hamming_distance(s1, s2):
    # Return the Hamming distance between equal-length sequences
    if len(s1) != len(s2):
        raise ValueError("Undefined for sequences of unequal length")
    else:
        return sum(ch1 != ch2 for ch1, ch2 in zip(s1, s2))

for i in xrange(len(in_genome)-kmer + 1):
    v = in_genome[i:i + kmer]
    out_result.append(v)

for i in xrange(len(out_result)):
    for j in xrange(i + 1, len(out_result)):
        if hamming_distance((out_result[i], out_result[j]) <= in_mistake:
            mismatch_list.append(out_result[i], out_result[j])
    print mismatch_list

当我运行代码时，它总是显示“无效语法”并指出错误来自“if hamming_distance((out_result[i], out_result[j])

另一件事是我知道肯定存在更好的算法。有什么线索吗？

Answer 1

我不确定你想用 mismatch_list 做什么：你不能同时附加两个东西。您是否需要extend 来添加它们中的每一个？引用的行中多余的左括号很容易修复；这需要更多的语义工作。用

替换这些行

if hamming_distance(out_result[i], out_result[j]) <= in_mistake:
    mismatch_list.extend([out_result[i], out_result[j]])

给出这个输出（不匹配列表）：

['ACGT', 'ATGT', 'GTTG', 'GTCG', 'GTTG', 'GATG', 'TTGC', 'TCGC', 'TTGC', 'ATGC', 'TGCA', 'CGCA', 'TGCA', 'TGCA', 'GCAT', 'GCAT', 'GCAT', 'GCAT', 'CATG', 'CATG', 'CATG', 'GATG', 'CATG', 'CATG', 'ATGT', 'ATGA', 'ATGT', 'ATGC', 'ATGT', 'ATGA', 'CGCA', 'TGCA', 'GCAT', 'GCAT', 'CATG', 'GATG', 'CATG', 'CATG', 'ATGA', 'ATGC', 'ATGA', 'ATGA', 'TGAT', 'TGAG', 'GATG', 'CATG', 'ATGC', 'ATGA', 'TGAG', 'AGAG', 'GAGA', 'GAGC']

---

我不确定你想从这里得到什么；您指定了 d=1，但您打印的项目各有 4 个错误。下面是一些代码，它计算 mismatch_list 中的出现次数，按频率对结果进行排序（未使用，但您以后可能需要它），然后打印所有与假设频率 4 匹配的结果。

mismatch_count = collections.Counter(mismatch_list)
sorted_match = sorted(list(mismatch_count), key=lambda pair: pair[1], reverse=True)
for item, value in mismatch_count.iteritems():
    if value == 4:
        print item