随机DNA突变生成器

Question

我想为一系列变异的DNA链创建一个字典的字典，每个字典都展示原始碱基和它变异后的碱基。

具体来说，我想做的是创建一个生成器，允许用户输入一个特定的DNA链，然后生成100条随机生成的链，这些链的变异频率是0.66%（这个频率适用于每个碱基，并且每个碱基可以变异成任何其他碱基）。接下来，我想创建一系列字典，每个字典详细说明在特定随机生成的链中发生的变异。我希望字典的键是原始碱基，值是新的变异碱基。有没有简单的方法可以做到这一点？到目前为止，我一直在尝试一个循环，代码大概是这样的：

#yields a strand with an A-T mutation frequency of 0.066%
def mutate(string, mutation, threshold):
    dna = list(string)
    for index, char in enumerate(dna):
        if char in mutation:
            if random.random() < threshold:
                dna[index] = mutation[char]

    return ''.join(dna)

dna = "ATGTCGTACGTTTGACGTAGAG"
print("DNA first:", dna)
newDNA = mutate(dna, {"A": "T"}, 0.0066)
print("DNA now:", newDNA)

但是我用这段代码只能生成一条链，而且只关注T变成A的变异。我也不太确定怎么把字典和这个结合起来。有没有人能给我一个更好的方法？谢谢。

Answer 1

听起来你的问题可以分成两个部分。第一部分是你想多次改变你的DNA序列，第二部分是你想收集一些关于这些变化的额外信息，放在某种数据结构里。我会分别来讲这两部分。

从同一个源字符串生成100个随机结果其实很简单。你可以用一个明确的循环（比如在一个生成器函数里），但你也可以用列表推导式来反复运行一个单一的变异函数：

results = [mutate(original_string) for _ in range(100)]

当然，如果你把mutate函数做得更复杂，这段简单的代码可能就不太合适了。如果它返回的是某种更复杂的数据结构，而不仅仅是一个字符串，你可能需要做一些额外的处理，把数据组合成你想要的格式。

至于如何构建这些数据结构，我觉得你现在的代码已经是个不错的开始。你需要决定你将如何访问这些数据，然后根据这个决定选择合适的容器。

例如，如果你只是想简单记录所有对字符串的变异，我建议用一个基本的列表，里面包含变异前后的元组。另一方面，如果你想高效地查找某个基因变异成什么，使用一个字典，值是列表，可能更合适。如果你想的话，还可以包含变异基因的索引。

这里有一个简单的函数尝试，它返回变异后的字符串以及记录所有变异的元组列表：

bases = "ACGT"

def mutate(orig_string, mutation_rate=0.0066):
    result = []
    mutations = []
    for base in orig_string:
        if random.random() < mutation_rate:
            new_base = bases[bases.index(base) - random.randint(1, 3)] # negatives are OK
            result.append(new_base)
            mutations.append((base, new_base))
        else:
            result.append(base)
    return "".join(result), mutations

这段代码中最棘手的部分是我如何选择当前基因的替代品。表达式bases[bases.index(base) - random.randint(1, 3)]一次性搞定了所有。我们来分解一下。bases.index(base)给出了在代码顶部的全局bases字符串中，前一个基因的索引。然后我从这个索引中减去一个随机偏移量（random.randint(1, 3)）。新的索引可能是负数，但没关系，因为当我们用它去索引bases字符串时（bases[...]），负索引是从右边开始计算的，而不是从左边。

这是你可以如何使用它：

string = "ATGT"
results = [mutate(string) for _ in range(100)]
for result_string, mutations in results:
    if mutations: # skip writing out unmutated strings
        print(result_string, mutations)

对于短字符串，比如"ATGT"，你很可能只会得到一次变异，甚至一次都不太常见。上面的循环每次运行通常会打印出2到4个结果（也就是说，超过95%的四个字符长度的字符串根本没有变异）。较长的字符串会更频繁地发生变异，而且在一个字符串中看到多个变异的可能性更大。