One genomic dataset, seven R packages

Data design, reproducibility and shared objects for the package tutorials.

These tutorials follow one connected analysis. gbase prepares a PLINK BED/BIM/FAM panel and a reusable LD reference. The other packages perform association testing, Bayesian regression, genetic correlation, fine-mapping, scoring and gene-level analysis. The pages show actual tables and figures from one seeded simulation, along with the R calls that produced them.

The simulated study

Component	Design
Genotypes	20,000 individuals and 50,000 biallelic markers across 500 artificial chromosome labels, one per 100-marker region
Traits	Three quantitative traits, A, B and C
GWAS	Three disjoint groups of 6,000 individuals, one group per trait
LD reference	The first 18,000 individuals; 500 independent 100-marker regions
Prediction	The remaining 2,000 individuals, excluded from GWAS and LD preparation
Gene sets	1,000 artificial genes, two nonoverlapping marker sets per region
Pathways	Three gene collections labelled enriched, background and mixed

The simulation uses deliberately strong signals and simple independent LD regions to make the outputs readable. Effect sizes, p-values and run times are properties of this one simulated dataset. They do not establish calibration or performance on real cohorts.

Reproduce the complete run

Install the seven source packages and gsim. Download the complete R driver, place it in a new working directory and run:

Rscript run.R --threads=4

The driver creates the PLINK files and phenotypes, computes LD once, and runs all seven package workflows. It saves stage objects under build/benchmarks/genomic-20k-50k-4cores/ and small result tables under benchmarks/genomic_20k_50k/results/. This is a substantial simulation; allow disk space for the packed genotypes and LD files. The tutorial plots and tables were generated from its saved four-thread stage objects. Package installation and compilation are outside the driver.

To inspect a saved stage without repeating the analysis:

cache <- file.path("build", "benchmarks", "genomic-20k-50k-4cores")
stage <- function(name) {
  readRDS(file.path(cache, paste0(name, ".rds")))$value
}

reference_Glist <- stage("sparse_ld_and_scores")
statistics <- stage("gwas_and_checkstat")$statistics
bayes_fit <- stage("bayesian_regression")

base_cache <- file.path("build", "benchmarks", "genomic-20k-50k")
simulation_cache <- if (file.exists(file.path(cache, "data_simulation.rds")))
  cache else base_cache
simulation <- readRDS(file.path(simulation_cache,
  "data_simulation.rds"))$value
prefix <- simulation$prefix
Y <- simulation$simulation$Y
rownames(Y) <- sprintf("id%05d", seq_len(nrow(Y)))

marker_ids <- sprintf("m%05d", 1:50000)
regions <- split(marker_ids,
  rep(sprintf("region%03d", 1:500), each = 100))
gene_sets <- list()
for (r in seq_along(regions)) {
  cut <- 20L + r %% 60L
  gene_sets[[sprintf("gene%04d", 2L * r - 1L)]] <- regions[[r]][1:cut]
  gene_sets[[sprintf("gene%04d", 2L * r)]] <- regions[[r]][(cut + 1L):100L]
}

regions and gene_sets are named marker vectors; pathway_sets in the gsea tutorial contains gene names. These are different units: a gene set identifies markers in the Glist, while a pathway identifies genes already summarized by VEGAS. In an ordinary study, replace prefix, Y, the sample IDs and set definitions with your own data.

Follow the analysis

Step	Tutorial	Read the result as
1	gbase: data and LD	A checked Glist, sparse LD, LD scores and aligned statistics
2	glma: association and VEGAS	Marker effects, p-values and LD-aware gene evidence
3	gbayes: Bayesian regression	Posterior marker effects, diagnostics and scoring weights
4	gcorr: heritability and correlation	Genome-wide and partitioned estimates with uncertainty
5	gmap: fine-mapping	Regional PIPs, posterior sets and fit diagnostics
6	gscore: prediction	Individual scores and marker coverage
7	gsea: genes and pathways	Pathway tests and feature-based gene ranks

Each page also points to other supported routes. Those calls illustrate the package interface; only results explicitly labelled as part of this run appear in the displayed tables and plots.