注意：我对此进行了重大编辑，从而将模拟时间从14小时缩短到了14分钟。
我刚开始学习编程，但我已经做出了一个模拟程序，试图模拟生物中的无性复制并量化母细胞和子细胞之间染色体数目的差异。但模拟运行极慢，需要约6个小时才能完成。因此，我想知道如何最好地加速模拟的运行速度。
这些数字生物有x条染色体。与大多数生物不同的是，这些染色体都是相互独立的，因此它们被传递到子代的机会是相等的。
在这种情况下，染色体进入子细胞的分布遵循概率为0.5的二项式分布。
函数sim_repo接受一个数字生物矩阵，其中包含已知染色体数量的数字生物，并让它们经过12代的复制。该函数复制这些染色体，然后使用rbinom函数随机生成一个数字。这个数字然后被分配给一个子细胞。由于在无性生殖过程中没有染色体丢失，因此另一个子细胞获得剩余的染色体。这个过程再经过G代重复，然后从矩阵中的每一行中随机抽取一个值。
sim_repo = function( x1, G=12, k=1, t=25, h=1000 ) {
# x1 is the list of copy numbers for a somatic chromosome
# G is the number of generations, default is 12
# k is the transfer size, default is 1
# t is the number of transfers, default is 25
# h is the number of times to replicate, default is 1000
dup <- x1 * 2 # duplicate the initial somatic chromosome copy number for replication
pop <- 1 # set generation time
set.seed(11)
z <- matrix(rbinom(n=rep(1,length(dup)),size = as.vector(dup),prob = 0.5),nrow = nrow(dup)) # amount of somatic chromosome is distributed to one of the daughter cells
z1 <- dup - z # as no somatic chromosomes are lost, the other daughter cells receives the remainder somatic chromosomes
x1 <- cbind(z, z1) # put both in a matrix
for ( pop in 1:G ) { # this loop does the replication for each cell in each generation
pop <- 1 + pop # number of generations. This is a count for the for loop
dup <- x1 * 2 # double the somatic chromosomes for replication
set.seed(11)
z <- matrix(rbinom(n=rep(1,length(dup)),size = as.vector(dup),prob = 0.5),nrow = nrow(dup)) # amount of somatic c hromosomes distributed to one of the daughter cells
z1 <- dup - z # as no somatic chromosomes are lost, the other daughter cells receives the remainder somatic chromosomes
x1 <- cbind(z, z1) # put both in a matrix
}
# the following for loop randomly selects one cell in the population that was created
# the output is a matrix of 1 column
x1 <- matrix(apply(x1, 1, sample, size=k), ncol=1)
x1
}
在我的研究中，我对最初祖先生物染色体方差和模拟结束时点的变化感兴趣。以下函数表示将细胞转移到新的生存环境中。它使用从函数sim_rep的输出来生成更多的世代。然后找到第一列和最后一列矩阵行之间的方差，并找到它们之间的差异。
# The following function is mostly the same as I talked about in the description.
# The only difference is I changed some aspects to take into account I am using
# matrices and not lists.
# The function outputs the difference between the intial variance component between
# 'cell lines' with the final variance after t number of transfers
sim_exp = function( x1, G=12, k=1, t=25, h=1000 ) {
xn <- matrix(NA, nrow(x1), t)
x <- x1
xn[,1] <- x1
for ( l in 2:t ) {
x <- sim_repo( x, G, k, t, h )
xn[, l] <- x
}
colvar <- matrix(apply(xn,2,var),ncol=ncol(xn))
ivar <- colvar[,1]
fvar <- colvar[,ncol(xn)]
deltavar <- fvar - ivar
deltavar
}
我需要复制这个模拟h次。因此，我编写了以下函数，该函数将调用函数sim_exph次。
sim_1000 = function( x1, G=12, k=1, t=25, h=1000 ) {
xn <- vector(length=h)
for ( l in 2:h ) {
x <- sim_exp( x1, G, k, t, h )
xn[l] <- x
}
xn
}
当我调用带有6个值的sim_exp函数时，它需要大约52秒才能完成。
x1 <- matrix(data=c(100,100,100,100,100,100),ncol=1)
system.time(sim_1000(x1,h=1))
user system elapsed
1.280 0.105 1.369
如果我能更快地获得它，那么我可以完成更多的模拟，并在模拟中应用选择模型。
我的输入将如下所示：
x1，一个矩阵，每个祖先生物都在其自己的行中。x1 <- matrix(data=c(100,100,100,100,100,100),ncol=1) # a matrix of 6 organisms
当我运行以下命令时：
a <- sim_repo(x1, G=12, k=1)
我的期望输出将是：
a
[,1]
[1,] 137
[2,] 82
[3,] 89
[4,] 135
[5,] 89
[6,] 109
system.time(sim_repo(x1))
user system elapsed
1.969 0.059 2.010
当我调用sim_exp函数时，
b <- sim_exp(x1, G=12, k=1, t=25)
它会调用sim_repo函数G次并输出：
b
[1] 18805.47
当我调用
sim_1000函数时，通常会将h设置为1000，但在这里我将其设置为2。因此，在这里，sim_1000将调用sim_exp并复制它两次。c <- sim_1000(x1, G=12, k=1, t=25, h=2)
c
[1] 18805.47 18805.47
sim_exp()内部的cbind()和sim_1000()内部的c()可能非常耗费资源。 - flodelsim_exp()中，我会制作一个与最终输出相同列数和行数的矩阵，但将值填充为NULL吗？ - Kevinsim_exp效果很好。速度快了很多。你会如何在sim_repo函数中预分配矩阵？由于我的循环是按照世代数进行的，所以我无法像你在sim_exp中建议的那样将下一列设置为输入。 - Kevin