生物信息维护的种群动态模拟

Question

背景

通过这个模拟，我研究了一个酶在细胞中增殖的系统。在酶的复制过程中，寄生虫可能是由突变引起的。它们会使系统灭绝。我感兴趣的是参数空间共存在哪里是可能的。

在这个程序中，系统是一个列表，细胞是带有两个键的字典:酶是"e"，寄生虫是"p"。键的值是两个变体的数目。

我们的参数是：

• pop_size：细胞的数量
• cell_size：细胞分裂的最大分子数(enzymes+parasites)
• a_p：寄生虫的适应度相对于酶的适应度(例如，如果是a_p = 2，寄生虫的适应度是酶的两倍)。
• mutation_rate：复制事件中发生突变的概率
• gen_max：最大的世代数(一个世代对应于一个世代)
• while循环；如果系统灭绝，程序直到gen_max才能运行)

我们从pop_size细胞开始，用cell_size // 2酶和0寄生虫。在每个细胞中，分子不断增殖，直到其数量达到cell_size。每个细胞分裂，分子的分类根据二项分布(p=0.5)进行。带有"e" < 2的单元格作为死单元被丢弃。在此之后，如果存活细胞的数量大于pop_size，我们根据细胞适应度("e"/("e"+"p"))选择细胞中的pop_size，并将其转移到下一代。另一方面，如果存活细胞的数量是pop_size或更少，它们都会转移到下一代。

我的请求

我从来没有在学校学过编程。这个程序是大量搜索的结果。现在我已经到了需要有经验的人给我建议的地步。在某些参数值下，程序变得相当慢。

1. 有什么更好的解决方案存在性能比我的解决方案，在整个程序的操作列表的项目和写入数据文件？以及算法设计？
2. 为了有效地实现这类模型，我应该在哪些方向上改进Python的编程技能？或者，我在这方面是否接近Python能力的极限？
3. 我是否应该改为一种更合适的编程语言，以便在这类任务中获得更好的性能？如果是，我应该考虑哪种语言？(我猜是C.)

该程序由两个功能组成。simulation()进行模拟，writeoutfile()将数据写入文件。

# -*- coding: utf-8 -*-
from random import choices, random
import csv
import time
import numpy as np


def simulation(pop_size, cell_size, a_p, mutation_rate, gen_max):
    def fitness(pop):
        return [i["e"] / (i["e"] + i["p"]) for i in pop]

    def output(pop, gen, pop_size, cell_size, mutation_rate, a_p, boa_split):
        if pop:
            gyaklist_e = [i["e"] for i in pop]
            gyaklist_p = [i["p"] for i in pop]
            fitnesslist = fitness(pop)
            return (
                gen,
                sum(gyaklist_e), sum(gyaklist_p),
                sum([1 for i in pop if i["e"] > 1]),
                np.mean(gyaklist_e), np.var(gyaklist_e),
                np.percentile(gyaklist_e, 25),
                np.percentile(gyaklist_e, 50),
                np.percentile(gyaklist_e, 75),
                np.mean(gyaklist_p), np.var(gyaklist_p),
                np.percentile(gyaklist_p, 25),
                np.percentile(gyaklist_p, 50),
                np.percentile(gyaklist_p, 75),
                np.mean(fitnesslist), np.var(fitnesslist),
                np.percentile(fitnesslist, 25),
                np.percentile(fitnesslist, 50),
                np.percentile(fitnesslist, 75),
                pop_size, cell_size, mutation_rate, a_p, boa_split
                )
        return (
            gen,
            0, 0,
            0,
            0, 0,
            0, 0, 0,
            0, 0,
            0, 0, 0,
            0, 0,
            0, 0, 0,
            pop_size, cell_size, mutation_rate, a_p, boa_split
            )

    pop = [{"e": cell_size // 2, "p": 0} for _ in range(pop_size)]
    gen = 0
    yield output(
        pop,
        gen, pop_size, cell_size, mutation_rate, a_p, boa_split="aft"
        )
    print(
        "N = {}, rMax = {}, aP = {}, U = {}".format(
            pop_size, cell_size, a_p, mutation_rate
            )
        )

    while pop and gen < gen_max:
        gen += 1

        for i in pop:
            while not i["e"] + i["p"] == cell_size:
                luckyreplicator = choices(
                    ["e", "p"], [i["e"], a_p*i["p"]]
                    )
                if luckyreplicator[0] == "e" and random() < mutation_rate:
                    luckyreplicator[0] = "p"
                i[luckyreplicator[0]] += 1

        if gen % 100 == 0:
            yield output(
                pop,
                gen, pop_size, cell_size, mutation_rate, a_p, boa_split="bef"
                )

        newpop = [
            {"e": np.random.binomial(i["e"], 0.5),
             "p": np.random.binomial(i["p"], 0.5)}
            for i in pop
            ]
        for i in zip(pop, newpop):
            i[0]["e"] -= i[1]["e"]
            i[0]["p"] -= i[1]["p"]

        pop += newpop
        newpop = [i for i in pop if i["e"] > 1]

        if newpop:
            fitnesslist = fitness(newpop)
            fitness_sum = np.sum(fitnesslist)
            fitnesslist = fitnesslist / fitness_sum
            pop = np.random.choice(
                newpop, min(pop_size, len(newpop)),
                replace=False, p=fitnesslist
                ).tolist()
        else:
            pop = newpop
            for i in range(2):
                yield output(
                    pop,
                    gen+i, pop_size, cell_size, mutation_rate, a_p, boa_split="aft"
                    )
            print("{} generations are done. Cells are extinct.".format(gen))

        if gen % 100 == 0 and pop:
            yield output(
                pop,
                gen, pop_size, cell_size, mutation_rate, a_p, boa_split="aft"
                )

        if gen % 1000 == 0 and pop:
            print("{} generations are done.".format(gen))


def writeoutfile(simulationresult, runnumber):
    localtime = time.strftime(
        "%m_%d_%H_%M_%S_%Y", time.localtime(time.time())
        )
    with open("output_data_" + localtime + ".csv", "w", newline="") as outfile:
        outfile.write(
            "gen"+";" +
            "eSzamSum"+";"+"pSzamSum"+";" +
            "alive"+";" +
            "eSzamAtl"+";"+"eSzamVar"+";" +
            "eSzamAKv"+";" +
            "eSzamMed"+";" +
            "eSzamFKv"+";" +
            "pSzamAtl"+";" + "pSzamVar" + ";" +
            "pSzamAKv"+";" +
            "pSzamMed"+";" +
            "pSzamFKv"+";" +
            "fitAtl"+";"+"fitVar"+";" +
            "fitAKv"+";" +
            "fitMed"+";" +
            "fitFKv"+";" +
            "N"+";"+"rMax"+";"+"U"+";"+"aP"+";"+"boaSplit"+"\n"
            )
        outfile = csv.writer(outfile, delimiter=";")
        counter = 0
        print(counter, "/", runnumber)
        for i in simulationresult:
            outfile.writerows(i)
            counter += 1
            print(counter, "/", runnumber)


RESULT = [simulation(100, 20, 1, 0, 10000)]
RESULT.append(simulation(100, 20, 1, 1, 10000))
N_RUN = 2
writeoutfile(RESULT, N_RUN)
# Normally I call the functions from another script,
# these last 4 lines are meant to be an example.

关于参数值

到目前为止，对这些数值的组合进行了审查：

• pop_size：100；200；500；1000
• cell_size：20；50；100；200；500；1000
• a_p：0.75；1；1.25；1.5；1.75；2；3
• mutation_rate：0-1
• gen_max：10000

主要是我想增加pop_size和1000以上的细胞，程序比我希望的慢。当然，这有点主观，但举例来说，一百万个细胞将是一个完全合理的假设，在这个数量级上，我认为它在客观上是不可能缓慢的。

随着cell_size的增加，程序也会变慢，a_p也会稍微慢一些，但目前我对前者的值感到满意，而后者的效果是可以容忍的。

突变率对速度的影响也是可以容忍的。

除了pop_size之外，gen_max还应该增加，并且对运行时有显著的影响。我知道我并不能捕捉到每一次10000代的灭绝事件。20000更好，50000就足够了，100000就像用大锤敲螺母一样。

Answer 1

Numpy可以非常快，接近C或其他低级语言的速度(因为它使用C!)。但条件是慢的东西实际上是在Numpy做的。我的意思是，你不能一直在列表和字典中循环，然后在Numpy中选择操作，你必须坚持Numpy数组和元素级操作。

我会给出一些关于风格的评论，然后再回到这一点。

• 首先，在整个代码中没有任何注释。我建议在函数开始时使用"""docstrings"""，在代码有点混乱的行之间使用简短的# Comments。
• F-字符串是python 3.6+特性，它极大地提高了可读性。它们用于代替.format()和字符串连接。例如：

print(f'{gen} generations are done. Cells are extinct.')

• 当更长的代码行变得更干净时，可以将大量代码分散在几行代码上。您没有非常高度嵌套的代码，所以行甚至不会那么长。
• 很好地利用了yield。这是一些新程序员经常跳过的东西，很高兴看到它被用于这里的效果。
• 您的导入是干净的，最小的，并且与代码的其余部分有很好的分离。
• 一些命名可能需要一些工作来帮助清晰。只需将您的键命名为enzyme和parasite，而不是e和p。a_p是什么？尽量不要使用内置函数名作为参数名称(pop)，因为它可能会导致问题和混淆。在这里，它显然是人口的缩写，但要小心。使用snake_case命名小写对象( ratherthanthis ).
• 您经常返回大量的值。如果您总是将0打印到文件中，则不需要将其返回，每次只需将其写入文件，然后编写其余的返回值。有些事情，如gen，应该保持外部跟踪，而不是每次都被返回。如果某物是静态的，您可能不需要将它输入函数，然后将其吐回未咀嚼的状态。
• 多行字符串可以通过三重引号来实现：

example = """
          Like
          This
          """

回到Numpy

• 正如我说的，要快速，您需要使用Numpy开始-在您的慢部分完成。如果使用纯python生成一个列表，然后将其转换为一个数组，然后将其返回到纯python，则通常不会节省时间。它甚至可以慢于单纯的蟒蛇。
• 例如，您的健身函数应该使用元素级操作。
• 如果将纯python中最慢的部分替换为纯Numpy，则应该会看到一些很好的改进。您可以尝试一个代码剖析器，以找到确切的挂起位置。