Несколько раз запустите симуляцию pop-gen и сохраните результаты каждого в новом столбце во фрейме данных. - PullRequest
Новая
       17

Несколько раз запустите симуляцию pop-gen и сохраните результаты каждого в новом столбце во фрейме данных.

1 голос
/ 21 марта 2019

У меня есть базовая симуляция Райта-Фишера для двух аллелей, которая прекрасно работает и дает хорошо выглядящий график, на котором аллель фиксируется или вымирает случайно, как и ожидалось. Он экспортирует каждое поколение в расчете на фрейм данных d, поэтому у меня есть значения для каждого поколения. То, что я хочу сделать, это запустить все это, скажем, 20 раз и автоматически сохранить каждую полную симуляцию в новом столбце, чтобы я мог построить их все на графике ggplot с цветами и всем этим хорошим материалом. В основном я заинтересован в том, чтобы получить аккуратную рамку для создания привлекательных сюжетов для проекта, а не для головокружительной эффективности. 

#Wright Fisher model Mk1

#Simulation Parameters
# n = pop.size
# f = frequency of focal allele
# x = number of focal allele, do not set by hand
# y = number of the other allele, do not set by hand
# g = generations desired
n = 200
f = 0.6
x = (n*f)
y = (n-x)
g = 200

#This creates a data frame of the correct size to store each generation

d = data.frame(f = rep(0,g))

#Creates the graph.
plot(1,0, type = "n", xlim = c(1,200), ylim = c(0,n),
     xlab = "Generation", ylab = "Frequency A")

#Creates the population, this model is limited to only two alleles, and can only plot one
alleles<- c(rep("A",x), rep("a",y))

#this is the loop that actually simulates the population
#It has code for plotting each generation on the graph as a point 
#Exports the number of focal allele A to the data frame
for (i in 1:g){ 
  alleles <- sample(alleles, n, replace = TRUE)
points(i, length(alleles[alleles=="A"]), pch = 19, col= "red")
F = sum(alleles == "A")
d[i, ] = c(F)
}

Итак, я хочу запустить этот последний бит несколько раз и как-то сохранить каждую полную итерацию. Я знаю, что мог бы зациклить функцию, вложив ее, даже если это быстро и грязно, но при этом сохраняются только значения последней итерации внешнего цикла.

1 Ответ

0 голосов
/ 21 марта 2019

Здесь много возможностей для улучшения, но это должно помочь вам.Я показываю только пять симуляций, но вы должны быть в состоянии расширить.По сути, поместите большую часть вашего кода в функцию, и тогда вы сможете использовать либо map функции из пакета purrr, либо вы также можете сделать что-то с replicate: 

library(tidyverse)

n = 200
f = 0.6
x = (n*f)
y = (n-x)
g = 200

d = data.frame(f = rep(0,g))

run_sim <- function() {
  alleles <- c(rep("A", x), rep("a", y))

  for (i in 1:g) { 
    alleles <- sample(alleles, n, replace = TRUE)
    cnt_A = sum(alleles == "A")
    d[i, ] = c(cnt_A)
  }

  return(d)
}

sims <- paste0("sim_", 1:5)

set.seed(4) # for reproducibility

sims %>%
  map_dfc(~ run_sim()) %>%
  set_names(sims) %>%
  gather(simulation, results) %>%
  group_by(simulation) %>%
  mutate(period = row_number()) %>%
  ggplot(., aes(x = period, y = results, group = simulation, color = simulation)) +
  geom_line()

Создано в 2019-03-21 с помощью пакета prex (v0.2.1)

ПРИМЕЧАНИЕ. Вы также можете добавить аргументы в run_sim функция, скажем, x и y (т.е. run_sim <- function(x, y) { ... }), которая позволит вам изучить другие возможности.

...