power_paper_sim_eqtl.R

library(dplyr)
library(tidyr)
library(purrr)

simulate_genotype <- function(maf, n) rbinom(n, 2, maf)

simulate_data <- function(beta, maf, n) {
  data.frame(
    genotype = simulate_genotype(maf, n)
  ) |>
    mutate(
      expression = genotype * beta + rnorm(n, mean = 0, sd = 1)
    )
}

n_per_combination <- 2000
beta_values <- c(0.05, 0.1, 0.2)
maf_values <- c(0.1, 0.2, 0.3)

combinations <- expand.grid(beta = beta_values, maf = maf_values)

simulated_data <- combinations %>%
  mutate(data = map2(beta, maf, ~simulate_data(.x, .y, n_per_combination))) %>%
  unnest(data)

means <- combinations %>%
  mutate(data = map2(
    beta, 
    maf, 
    ~ data.frame(
      genotype=factor(0:2), 
      expression=0:2 * .x
      )
    )) %>%
  unnest(data)

combo_with_r <- combinations |>
  mutate(r = round(beta * sqrt(2 * maf * (1 - maf)), 3),
         label = paste("r =", r))

library(ggplot2)
library(ggforce)
simulated_data |>
  mutate(genotype = factor(genotype)) |>
  ggplot(aes(genotype, expression)) +
  geom_violin() +
  geom_sina(alpha=.1) +
  geom_point(
    data = means,
    color = "blue",
    size = 4
  ) +
  geom_text(
    data = combo_with_r, 
    aes(label = label, x = Inf, y = Inf), 
    hjust = 1.1, vjust = 1.5, size = 8,
    inherit.aes = FALSE
  ) +
  coord_cartesian(ylim=c(-2,2)) +
  facet_grid(
    beta ~ maf,  
    labeller = labeller(beta=label_both, maf=label_both)) + 
  theme(
    axis.title = element_text(size = 20),
    strip.text = element_text(size = 20)
  )