Shared genomic simulation design

This design fixes the inputs for ten focused software validation runs across glma, gcorr, gbayes, gscore, gmap and gsea. The retained compact results were produced from this fixed design. The current LOCO scan and report driver use the current package interfaces. Website rendering never runs these analyses.

Population and truth

Each seed generates 10,000 unrelated individuals, 50,000 biallelic markers in 500 independent artificial LD regions, and three continuous traits. Disjoint groups of 3,000 supply one GWAS per trait; 1,000 other individuals are reserved for scoring. The pooled 9,000 GWAS individuals supply the matched LD reference. The workflow writes the genotypes blockwise to a packed BED file, then uses gsim to generate phenotypes from declared effects and the prepared panel. It records causal effects, genetic values, residual covariance and sample assignments. Population LD and the generating effects determine the heritability, covariance, correlation and marginal-effect truth; fitted LD does not define truth.

One data-generating process contains three named evaluation patterns:

  1. Distributed effects: effects across regions permit genome-wide GWAS, heritability, genetic correlation, Bayesian regression and prediction checks.
  2. Shared and distinct regional signals: the first two regions have known shared and distinct causal variants for fine-mapping and colocalisation.
  3. Annotation enrichment: causal effects are more frequent in the first 80 regions, with predefined artificial genes and pathways.

These are patterns within each simulated dataset, not three independent phenotype-generating scenarios. Ten fixed seeds are prescribed; failed seeds are recorded and never replaced. The run reuses the established connected workflow and saves compact results by seed, with hashes, package versions, warnings, status and stage timings. Large genotypes, fits and LD resources remain in the ignored build cache. Two R workers process odd and even seeds concurrently; elapsed times describe this machine and workload, not an isolated benchmark.

Generate a replicate in R

The code below shows the data-generation part of the connected workflow. Run it from the gsuite repository root with development installations of gsuite and gsim. Change replicate from 1 through 10 to use the prescribed seeds. Genotypes are written one marker at a time to PLINK BED files; gsim then generates the phenotypes from known effects. The recorded design also includes population truth, prepared LD and all six analysis outputs.

.libPaths(c("build/r-library", .libPaths()))
library(gsuite)
library(gsim)

replicate <- 1L
seed <- 20262000L + 100L * replicate
set.seed(seed)
n <- 10000L; m <- 50000L; size <- 100L; nr <- m %/% size
traits <- c("A", "B", "C")
ids <- sprintf("m%05d", seq_len(m))
individuals <- sprintf("id%05d", seq_len(n))
region <- rep(sprintf("region%03d", seq_len(nr)), each = size)
rho <- seq(.2, .8, length.out = nr)
rows <- setNames(lapply(0:2, function(k) k * 3000L + seq_len(3000L)), traits)
holdout <- 9001:10000

out <- file.path("build/examples/simulated-genomics",
                 sprintf("replicate-%02d", replicate))
dir.create(out, recursive = TRUE, showWarnings = FALSE)
prefix <- file.path(out, "genotypes")
con <- file(paste0(prefix, ".bed"), "wb")
writeBin(as.raw(c(0x6c, 0x1b, 0x01)), con)
for (j in seq_len(m)) {
  w <- rbinom(n, 2, .3)
  r <- (j - 1L) %/% size + 1L
  if ((j - 1L) %% size) {
    keep <- runif(n) < rho[r]
    w[keep] <- previous[keep]
  }
  previous <- w
  codes <- c(3L, 2L, 0L)[w + 1L]
  writeBin(as.raw(colSums(matrix(codes, 4L) * c(1, 4, 16, 64))), con)
}
close(con)
write.table(data.frame(rep(seq_len(nr), each = size), ids, 0,
                       as.character(rep(seq_len(size) * 1000L, nr)), "A", "G"),
            paste0(prefix, ".bim"), quote = FALSE, row.names = FALSE,
            col.names = FALSE)
write.table(data.frame(individuals, individuals, 0, 0, 0, -9),
            paste0(prefix, ".fam"), quote = FALSE, row.names = FALSE,
            col.names = FALSE)
Glist <- gbase::gprep(bedfiles = paste0(prefix, ".bed"))

active <- runif(m) < ifelse(rep(seq_len(nr), each = size) <= 80, .08, .008)
B <- matrix(0, m, 3, dimnames = list(ids, traits))
correlation <- matrix(.6, 3, 3); diag(correlation) <- 1
B[active, ] <- matrix(rnorm(sum(active) * 3, sd = .035), sum(active), 3) %*%
  chol(correlation)
B[1:200, ] <- 0
B[50, ] <- .5
B[125, 1] <- .5; B[175, 2:3] <- .5

simulation <- gsim::gsim(
  Glist = Glist, nt = 3, architecture = "fixed", beta = B,
  h2 = c(.3, .4, .35), vg = c(.3, .4, .35), re = 0,
  standardize_W = FALSE, scale_effects = TRUE, seed = seed + 1L,
  compute_sumstats = FALSE, chunk_size = 64
)
colnames(simulation$Y) <- colnames(simulation$G) <-
  colnames(simulation$B) <- traits
saveRDS(list(seed = seed, Glist = Glist, simulation = simulation,
             gwas_rows = rows, holdout_rows = holdout, region = region, rho = rho),
        file.path(out, "simulation.rds"))

The first 80 regions have a higher causal-marker probability; the first two regions also contain the declared shared and distinct signals. The GWAS groups, holdout and pooled LD reference are fixed by the row indices above. The saved object includes generating effects and individual genetic values. The connected workflow additionally computes population heritability and genetic-correlation truth from those effects and the known LD construction.

Checks and interpretation

Exact linear GWAS results are checked against selected R ordinary least-squares calculations within \(10^{-9}\). Applying the true marker effects as a score must reproduce centered simulated genetic values within \(10^{-9}\). All requested results, unavailable estimates and numerical statuses are retained. Point estimates, intervals, prediction, posterior diagnostics, causal-marker coverage and pathway results are summarized across the ten attempts without replacing unfavorable outcomes.

The run status and numerical checks list each prescribed seed. The library pages link the compact result tables used in their summaries.

The glma Manhattan and effect plots use per-marker results retained from the prescribed tenth seed. They illustrate one dataset; the ten-run summaries remain the repeated evidence.

The connected 50K-marker BayesC and BayesR fits hold residual variance, effect-variance scale and mixture weights at declared values. They check posterior effects, heritability summaries, chain behavior and held-out scores under that contract. The separate reference-LD parameter-learning checks is the evidence for learning these parameters with a matched independent reference and its prescribed sparse-LD construction.

One additional glma LOCO timing run uses the genotypes and trait A phenotype from seed 2, grouping the 500 artificial regions into 20 artificial chromosomes. It tests all 50,000 markers against a 500-marker background with two operation-local threads. This one run checks execution and timing. The unrelated-individual design and subset background do not qualify mixed-model null-tail calibration or heritability estimation. LOCO is not part of the ten repeated library comparisons.

Ten replicates support descriptive truth comparisons, not precise coverage, power, false-positive-rate or method-ranking claims. The LD reference overlaps the GWAS individuals. Binary traits, population structure, ancestry mismatch, independent external LD panels and complete package/platform validation are outside this design.