如何实现依赖前次结果的N次循环？以HWE杂合子模拟为例

我需要重复以下流程100次以获取不同的heterozygousHWE值，但每次循环必须基于前一次的结果执行——即第二次循环要使用第一次的结果，确保heterozygousHWE在每次循环中递减。目前我能实现100次循环，但每次都是独立执行，未关联前次结果，不知如何修改代码。

原始代码 #

n11 <- 10
n12 <- 12
n22 <- 3
#total de individuos
N <- n11 + n12 + n22

#calculos de frecuecnias alelicas
n1 <- (2*n11) + n12
n2 <- (2*n22) + n12

p1 = n1 /(2*N)
p1

p2 = 1- p1
p2
#TODO: Calcular p-value prueba chi cuadrada

#ESPERADOS 
n11e= N * (p1^2)
n12e = 2*N *p1*p2
n22e= N*((p2)^2)

#observados - esperados 
restn11 = n11- n11e
restn11

restn12 = n12- n12e
restn12

restn22 = n22- n22e
restn22

# (O-E)^2/2
x11= (restn11^2)/n11e
x11
x12= (restn12^2)/n12e
x12
x22= (restn22^2)/n22e
x22

#suma total
TOTAL = x11 + x12 + x22

#x^2
TOTAL

#valor P (p-value)
vp1gl= 0.4549
vp2gl= 1.3863

pchisq(TOTAL ,df=1)

#Simular poblacion HWE
#numbers GENERA 200 n aleatorios
#A1 cuantas personas son mayores a p1
numbers <- runif(N)
numbers

a1 <- numbers>p1
a1

numbers <- runif(N)
a2 <- numbers>p1
a2

genotypes <- a1 + a2
heterozygousHWE <- sum (genotypes==1)
heterozygousHWE
    

get.bias <- function(i) {
  numbers <- runif(N)
  numbers
  
  a1 <- numbers>p1
  a1
  
  numbers <- runif(N)
  a2 <- numbers>p1
  a2
  
  genotypes <- a1 + a2
  heterozygousHWE <- sum (genotypes==1)
  heterozygousHWE
}
set.seed(1)
result <- t(sapply(1:1000,get.bias))
m <- head(result)
m

修改方案 #

原来的sapply调用是每次独立执行模拟，没有传递前一次的群体状态，要改成迭代式循环，每次基于上一轮的基因型结果更新等位基因频率，再进行下一轮模拟。具体步骤：

1. 初始化存储结果的向量，记录每一轮的heterozygousHWE值
2. 第一轮用初始参数完成模拟并记录结果
3. 后续99轮循环中，每一轮都基于上一轮的基因型重新计算基因型计数、等位基因频率，再执行模拟并记录结果

修改后的完整代码：

set.seed(1) # 固定随机种子，保证结果可复现

# 初始参数
n11 <- 10
n12 <- 12
n22 <- 3
N <- n11 + n12 + n22

# 存储100轮的heterozygousHWE结果
hwe_results <- numeric(100)

# 第一轮模拟：基于初始参数
n1 <- (2*n11) + n12
p1 <- n1 / (2*N)
a1 <- runif(N) > p1
a2 <- runif(N) > p1
genotypes <- a1 + a2
hwe_results[1] <- sum(genotypes == 1)

# 后续99轮迭代：完全依赖前一轮结果
for (i in 2:100) {
  # 从当前基因型计算新的纯合/杂合计数
  current_n11 <- sum(genotypes == 0)
  current_n12 <- sum(genotypes == 1)
  current_n22 <- sum(genotypes == 2)
  
  # 计算新的等位基因频率p1
  new_n1 <- (2*current_n11) + current_n12
  new_p1 <- new_n1 / (2*N)
  
  # 基于新p1模拟下一轮基因型
  a1 <- runif(N) > new_p1
  a2 <- runif(N) > new_p1
  genotypes <- a1 + a2
  
  # 记录当前轮的杂合子数量
  hwe_results[i] <- sum(genotypes == 1)
}

# 查看前10轮结果验证递减趋势
head(hwe_results, 10)

关键说明 #

• 每次循环通过genotypes传递前一轮的群体状态，确保模拟的连续性
• 若需要严格保证heterozygousHWE递减，可在循环中加入判断逻辑（比如当前值不小于前一个则重新模拟），但基于群体遗传逻辑，只要每一轮p1向纯合方向偏移，就能实现递减

内容的提问来源于stack exchange，提问作者Viviana Robayo