使用 Python 在 DNA 链上生成突变频率

Question

我想输入一个 DNA 序列并制作某种生成器，该生成器会产生具有一定突变频率的序列。例如，假设我有 DNA 链“ATGTCGTCACACACCGCAGATCCGTGTTTGAC”，我想创建 T->A 频率为 5% 的突变。我将如何去创造这个？我知道可以使用这样的代码来创建随机突变：

import random
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)

但我真的不知道该怎么做是固定突变频率。有人知道该怎么做吗？谢谢。

编辑：

如果我使用字节数组，格式应该是这样的，因为我收到了一个错误：

import random

dna = "ATGTCGTACGTTTGACGTAGAG"

def mutate(dna, mutation, threshold):
    dna = bytearray(dna) #if you don't want to modify the original
    for index in range(len(dna)):
        if dna[index] in mutation and random.random() < threshold:
                dna[index] = mutation[char]

    return dna

变异（DNA，{“A”：“T”}，0.05） print("我现在的 dna:", dna)

错误：“TypeError：没有编码的字符串参数”

编辑 2：

import random

myDNA = bytearray("ATGTCGTCACACACCGCAGATCCGTGTTTGAC")

def mutate(dna, mutation, threshold):
    dna = myDNA # if you don't want to modify the original
    for index in range(len(dna)):
        if dna[index] in mutation and random.random() < threshold:
                dna[index] = mutation[char]

    return dna

mutate(dna, {"A": "T"}, 0.05)
print("my dna now:", dna)

产生错误

Answer 1

你问我一个打印所有可能突变的函数，就是这样。输出的数量随着输入数据的长度呈指数增长，因此该函数仅打印可能性并且不会以某种方式存储它们（这可能会消耗大量内存）。我创建了一个递归函数，这个函数不应该与非常大的输入一起使用，我还将添加一个非递归函数，应该没有问题或限制。

def print_all_possibilities(dna, mutations, index = 0, print = print):
    if index < 0: return #invalid value for index

    while index < len(dna):
        if chr(dna[index]) in mutations:
            print_all_possibilities(dna, mutations, index + 1)
            dnaCopy = bytearray(dna)
            dnaCopy[index] = ord(mutations[chr(dna[index])])
            print_all_possibilities(dnaCopy, mutations, index + 1)

            return
        index += 1
    print(dna.decode("ascii"))

# for testing
print_all_possibilities(bytearray(b"AAAATTTT"), {"A": "T"})

这适用于我在 python 3 上，如果你愿意，我也可以解释代码。
注意：此函数需要函数测试中给出的字节数组。

说明：
该函数在 dna 中搜索可能发生突变的位置，它从索引开始，因此它通常从 0 开始并到末尾。这就是为什么每次执行循环时都会增加索引的while循环是for（它基本上是像for循环一样的正常迭代）。如果函数找到一个可能发生突变的地方（if chr(dna[index]) in mutations:），那么它会复制 dna 并让第二个 dna 发生突变（dnaCopy[index] = ord(mutations[chr(dna[index])])，注意 bytearray 是一个数值数组，所以我使用 chr 和 ord一直在 string 和 int 之间切换）。之后再次调用该函数以寻找更多可能的突变，因此这些函数再次在两个可能的 dna 中寻找可能的突变，但它们跳过了它们已经扫描的点，因此它们从 index + 1 开始。之后打印的命令被传递给调用的函数 print_all_possibilities，所以我们不必再做任何事情并使用return 退出执行。如果我们不再发现任何突变，我们将打印我们可能的 dna，因为我们不再调用该函数，所以没有其他人会这样做。
这听起来可能很复杂，但它或多或少是一个优雅的解决方案。此外，要理解递归，您必须了解递归，所以如果您现在不理解它，请不要打扰自己。如果您在一张纸上尝试一下，可能会有所帮助：取一个简单的 dna 字符串“TTATTATTA”，其中可能发生突变“A”->“T”（所以我们有 8 个可能的突变），然后执行以下操作：遍历字符串从左到右，如果你找到一个序列可以变异的位置（这里只是“A”），再次写下这个字符串，这次让字符串在给定位置变异，这样你的第二个字符串与原始字符串略有不同。在原件和副本中，标记你走了多远（可能在你让变异的字母后面加一个“|”），然后将副本作为新的原件重复此过程。如果您没有发现任何可能的突变，则在字符串下划线（这相当于打印它）。最后，您应该有 8 个不同的字符串，所有字符串都带有下划线。我希望这可以帮助理解它。

编辑：这是非递归函数：

def print_all_possibilities(dna, mutations, printings = -1, print = print):
    mut_possible = []
    for index in range(len(dna)):
        if chr(dna[index]) in mutations: mut_possible.append(index)

    if printings < 0: printings = 1 << len(mut_possible)
    for number in range(min(printings, 1 << len(mut_possible)):
        dnaCopy = bytearray(dna) # don't change the original
        counter = 0
        while number:
            if number & (1 << counter):
                index = mut_possible[counter]
                dnaCopy[index] = ord(mutations[chr(dna[index])])
                number &= ~(1 << counter)
            counter += 1

        print(dnaCopy.decode("ascii"))

# for testing
print_all_possibilities(bytearray(b"AAAATTTT"), {"A": "T"})

这个函数自带一个附加参数，可以控制最大输出的个数，例如

print_all_possibilities(bytearray(b"AAAATTTT"), {"A": "T"}, 5)

只会打印 5 个结果。

说明：
如果你的 dna 有 x 个可能发生突变的位置，那么你就有 2 ^ x 个可能的突变，因为在每个地方 dna 都可以发生突变。这个函数找到你的 dna 可以变异的所有位置，并将它们存储在mut_possible（这是 for 循环的代码）中。现在mut_possible 包含 dna 可以突变的所有位置，因此我们有 2 ^len(mut_possible)（len(mut_possible) 是 mut_possible 中的元素数量）可能的突变。我写了1 << len(mut_possible)，它是一样的，但更快。如果printings 是负数，该函数将决定打印所有可能性并将打印设置为可能性的数量。如果printings 是正数，但低于可能性的数量，那么该函数将只打印printings 突变，因为min(printings, 1 << len(mut_possible)) 将返回较小的数字，即printings。否则，该功能将打印出所有可能性。现在我们有number 可以通过range(...)，所以这个每次打印一个突变的循环将执行所需的次数。此外，数字每次都会增加一。 （例如，range(4) 与 [0, 1, 2, 3] 类似！）。接下来我们使用数字来创建一个突变。要了解此步骤，您必须了解二进制数。如果我们的数字是 10，它是二进制 1010。这些数字告诉我们在哪些地方我们必须修改 dna 的代码 (dnaCopy)。第一位是0，所以我们不修改第一个可能发生突变的位置，下一位是1，所以我们修改这个位置，之后是0等等...... " 我们使用变量counter 的位。如果设置了counterth 位，number & (1 << counter) 将返回一个非零值，因此如果设置了该位，我们将在可能发生突变的counterth 位置修改我们的 dna。这是用 mut_possible 写的，所以我们想要的位置是mut_possible[counter]。在我们对那个位置的 dna 进行变异后，我们将该位设置为 0，以表明我们已经修改了这个位置。这是通过number &= ~(1 << counter) 完成的。之后我们增加计数器来查看其他位。只有当 number 不为 0 时，while 循环才会继续执行，因此如果 number 至少设置了一位（如果我们必须修改至少一个 dna 的位置）。修改 dnaCopy 后，while 循环结束，我们打印结果。

我希望这些解释能有所帮助。我知道您是 python 新手，所以请花点时间了解一下，如果您有任何其他问题，请与我联系。

Answer 2

读完后，这个问题似乎很容易回答。我误解的可能性很高，所以如果我错了，请纠正我。

如果您希望有 5% 的机会将 T 更改为 A，那么您应该编写

mutate(yourString, {"A": "T"}, 0.05)

我还建议您使用字节数组而不是字符串。字节数组类似于字符串，它只能包含字节（值从 0 到 255），而字符串可以包含更多字符，但字节数组是可变的。通过使用字节数组，您无需创建临时列表或最终加入它。如果这样做，您的代码将如下所示：

import random
def mutate(dna, mutation, threshold):
    if isinstance(dna, str):
        dna = bytearray(dna, "utf-8")
    else:
        dna = bytearray(dna)

    for index in range(len(dna)):
        if chr(dna[index]) in mutation and random.random() < threshold:
                dna[index] = ord(mutation[chr(dna[index])])

    return dna

dna = "ATGTCGTACGTTTGACGTAGAG"
print("DNA first:", dna)
newDNA = mutate(dna, {"A": "T"}, 0.05)
print("DNA now:", newDNA.decode("ascii")) #use decode to make newDNA a string

在 bytearray 版本遇到了所有愚蠢的问题之后，这里是对字符串进行操作的版本：

import random
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.05)
print("DNA now:", newDNA)

如果您使用具有较大输入的字符串版本，则计算时间和使用的内存都会更大。当您想使用更大的输入执行此操作时，bytearray-version 将是最好的。