Gene-set and pathway analysis

gsea() connects prepared gene-level evidence to the native gsea library. The same named sets and explicit gene background support several questions:

Method	Question	Main result
ora	Are selected genes over-represented in a set?	Hypergeometric upper-tail p-value and fold enrichment
preranked	Does a set concentrate towards either end of a gene ranking?	ES, NES, leading edge and gene-label permutation p-value
competitive	Is a gene score associated with a pathway after covariate adjustment?	GLS coefficient, standard error and upper-tail enrichment p-value
magma	What does the bounded MAGMA summary workflow estimate from prepared gene p-values?	Probit-score preprocessing and competitive set regression
bayesc, bayesr	Which pathway effects are supported under a joint shrinkage model?	Posterior means, variances, PIPs and chain summaries

The last two choices also accept several traits fitted jointly. Bayesian PIPs are probabilities under the specified model and priors; they are not p-values. Attach gsea and use its re-exported vegas() for the separate marker-set association test below. The implementation and native integration remain in glma; users do not need to attach both packages. Use gbase::checkStat() when summary statistics need mapping or harmonisation first. Gene-ID conversion and allele harmonization remain explicit external preparation. The capability overview describes the separate native gene-evidence operations.

From GWAS summary statistics to pathways

Download the complete GWAS-to-pathway example. Run it with installed gsuite, gsea and gsim packages from the repository root:

source("tools/examples/gwas_pathways.R")

The script simulates two continuous traits in 10,000 individuals, with 240 markers assigned to 80 artificial genes across four artificial chromosomes. It prepares LD, computes ordinary GWAS Z statistics, derives gene evidence and runs pathway analyses. Genotypes are only needed to construct this demonstration; vegas() itself consumes ordinary checked summary-statistics data frames and the sparse-LD resource attached to a Glist. Outputs stay together in build/examples/gwas-pathways/: a saved workflow, one pathway table, one figure and the reference files.

The analysis boundary is explicit:

library(gsea)

evidence <- vegas(
  statistics, Glist, sets = genes,
  metadata = metadata,
  control = list(tail_method = "controlled_series"),
  threads = 4, verbose = TRUE
)

input <- gsea_input(evidence)

fit <- gsea(input$stat, sets, method = "preranked",
            control = list(seed = 31, replicates = 999, adjustment = "BH"))
fit$pathways
plot(fit, trait = "trait_A")

Each statistics entry is an ordinary data frame in complete sparse-LD marker and allele order, using qgg-style names such as rsids, ea, nea, stat, b, and seb. gbase::checkStat() can establish that order and orientation from a Glist. Z-score, effect/standard-error and score/information forms are accepted. A named list supplies traits; native set preparation is reused across them. When the same reference and gene definitions serve multiple analyses, retain that work explicitly:

prepared_genes <- prepare_vegas(Glist, genes, threads = 4)
evidence_A <- vegas(statistics_A, prepared = prepared_genes)
evidence_B <- vegas(statistics_B, prepared = prepared_genes)

The prepared object is an ordinary R object and can be saved with saveRDS(). sets = genes supplies the marker membership for each gene, for example list(gene1=c("m1","m2"), gene2=c("m3","m4")). Nonnegative quadratic weights are supplied separately through weights. The same vegas() interface accepts regions or other marker sets. Its set column is mapped explicitly to the gene-labelled schema expected by gsea() by gsea_input() as shown above. The adapter rejects unavailable p-values and inconsistent traits. The helper verifies identities and allele coding but does not harmonize them.

sets alone defines the regions. Native code opens one validated sparse-LD reader and prepares eigenvalue spectra in bounded parallel batches across genes. Trait evaluation reuses those spectra and does not reread LD. The default maximum is 2,048 markers per gene.

The returned evidence$stat includes the standardized quadratic gene score, p_value, tail diagnostics, marker count and SNP-wise mean metadata. VEGAS does not construct cross-gene covariance. ORA and preranked enrichment can consume the gene statistics directly. Competitive, MAGMA and Bayesian pathway models require a separately prepared, labelled gene covariance; marker-set membership is not treated as evidence of independence.

MAGMA’s sample-size and MAC adjustments additionally require explicit metadata columns gene, sample_size, and mean_minor_allele_count in gene order. A named trait list can provide different GWAS metadata per trait. The example calculates these from its GWAS individuals. The LD-reference sample count is provenance and is never substituted for GWAS sample size. Weighted genes do not acquire unit-weight SNP-wise mean metadata automatically.

Unavailable native tails remain visible through p_available, unavailable_reason, tail_error_bound and tail_terms. They are not silently replaced or filtered. The MAGMA route rejects nonfinite or endpoint gene p-values. A score-based model can be selected explicitly, but answers a different question. Singular gene dependence also remains an error in downstream GLS.

The script plots gene evidence together with ORA and preranked pathway results. Pathway A contains effects in both traits; pathway B contains effects only in trait B; pathway C contains none. This is one reproducible functional example, not a calibration study. Covariance-based pathway methods remain available through their explicit sampling contract, but are not included in this direct VEGAS handoff.

A small gene-level example

Download the runnable R example. It contains all six methods and uses 120 artificial genes and three overlapping pathways. This is an interface demonstration, not a comparative calibration study.

library(gsuite)
set.seed(2401)
genes <- paste0("gene", seq_len(120))
sets <- list(pathway_A=genes[1:25], pathway_B=genes[21:45],
             pathway_C=genes[71:95])
X <- vapply(sets, function(s) as.numeric(genes %in% s), numeric(120))
noise <- matrix(rnorm(240),120,2)
stat <- list(
  trait_A=data.frame(gene=genes, score=2*X[,1]+noise[,1]),
  trait_B=data.frame(gene=genes, score=1.5*X[,1]+X[,2]+
                      .3*noise[,1]+sqrt(.91)*noise[,2])
)

stat is a data frame for one trait or a named list for several traits. Every frame must contain the same unique character gene IDs in the same order. All genes in the declared background remain in the analysis, including genes outside every set. Unknown or duplicate set members are errors. Set overlap is allowed. For competitive and Bayesian regression, a labelled numeric feature matrix can replace binary set membership.

Selected genes or a ranked list

ora <- gsea(transform(stat$trait_A, selected=score>2), sets, "ora",
            control=list(adjustment="BH"))
ranked <- gsea(stat$trait_A, sets, "preranked",
               control=list(replicates=999, seed=2401, adjustment="BH"))
ranked$pathways
ranked$details$trait$leading_edge
plot(ranked)

ORA requires a logical selected column over the complete background; the selection rule is the caller’s decision. Preranked GSEA uses a finite score, a stable decreasing ranking and whole-background gene-label permutations. It does not correct for LD or perform phenotype permutations. Excluded sets and unavailable directional inference remain in the result with diagnostics. adjustment="BH" or "BY" adjusts the tested p-values, not pooled NES values as in Broad GSEA’s FDR calculation. Defaults are no adjustment, exponent 1, set sizes 1–500, seed 1, 999 permutations and one thread.

Competitive regression

fit <- gsea(stat, sets, "competitive",
             sampling=list(gene_covariance="independent"),
             control=list(adjustment="BH"))
fit$pathways
summary(fit)
plot(fit, trait="trait_A")

Classical methods analyse traits separately, with adjustment within each trait. Competitive regression fits pathways jointly by default; joint=FALSE fits each pathway separately with the same covariates and background. Covariates are a numeric matrix with gene row names or a data frame with gene and numeric columns; an intercept is added internally. Encode factors explicitly.

Regression always requires sampling$gene_covariance. Supply "independent" only when that working assumption is intended, or provide a labelled covariance matrix / named list of labelled blocks that partition the complete background. Block labels are used explicitly to map to gene order. Blocks must be symmetric positive definite; singular designs fail rather than silently drop columns. Student-t inference is a working approximation for non-Gaussian gene scores.

The bounded MAGMA route

gene_p <- transform(stat$trait_A,
  p_value=pnorm(score,lower.tail=FALSE), marker_count=10,
  effective_parameters=8, sample_size=10000, mean_minor_allele_count=2000)
magma <- gsea(gene_p, sets, "magma",
               sampling=list(gene_covariance="independent"),
               control=list(adjustment="BH"))
magma$pathways
magma$preparation$trait

These artificial p-values illustrate the API. Real analyses require calibrated gene-level evidence. Native preprocessing applies the documented probit transformation, clipping and gene-property corrections; the default adjusts the outcome for internal covariates before fitting each pathway separately. Constant internal covariates are reported. correct_size, correct_density, correct_sample_size and correct_mac default to TRUE; their positive input columns are required whenever enabled. truncate and residualize_internal also default to TRUE; bounds are lower_z_bound=-3, upper_z_sd=6 and covariate_sd_bound=5. External custom covariates use competitive instead.

For this method, gene dependence means quadratic-score correlation, with nonnegative entries and unit diagonal. The native workflow applies the recorded working-correlation power (default 1.1). This is not exact covariance of the transformed probit scores. It is a defined subset of MAGMA’s workflow, not full MAGMA equivalence. SNP-to-gene tests, allele alignment and reference preparation are separate native operations and are not hidden inside this call.

PoPS-style gene ranking

pops() is separate from pathway testing. It learns a joint feature model of gene evidence and predicts every gene from a fit that excluded the gene’s chromosome. Feature selection, nuisance adjustment and ridge-penalty selection are therefore repeated for every chromosome.

genes <- paste0("gene", seq_len(60))
chromosome <- rep(as.character(1:3), each = 20)
features <- cbind(
  expression = seq(-2, 2, length.out = 60),
  constraint = rep(c(-1, 1), 30),
  distance = cos(seq_len(60) / 6)
)
rownames(features) <- genes

gene_stat <- data.frame(
  gene = genes,
  chromosome = chromosome,
  score = 1.5 * features[, "expression"] + sin(seq_len(60) / 4)
)

ranking <- pops(
  gene_stat,
  features,
  sampling = list(gene_covariance = "independent"),
  control_features = "expression"
)

ranking$predictions[order(ranking$predictions$genome_rank), ]
ranking$folds
ranking$features
plot(ranking)

Use gene_covariance="independent" only when that assumption is intended. For correlated gene evidence, pass named covariance blocks using the same sampling contract as gsea(); every block must stay within one chromosome. The input score may come from a checked MAGMA-style gene analysis, gsea::vegas(), or another checked gene-level method, but pops() does not perform SNP-to-gene mapping or infer chromosome labels.

Bayesian PoPS uses the same chromosome-held-out prediction rule and the existing native gsea/gbayes sampler:

ranking <- pops(gene_stat, features, method="bayesC",
  pi=.001, nit=5000, nburn=1000, nrun=2)
ranking$predictions
ranking$folds
ranking$fits[[1]]$diagnostics

ranking <- pops(gene_stat, features, method="bayesR",
  pi=.001, nit=5000, nburn=1000, nrun=2)

sets=sets can replace the feature matrix. The familiar magma(stat, sets, method="bayesC", pi=..., nit=..., nburn=...) adapter also uses the existing pathway engine; it fits pathways rather than held-out gene predictions. These familiar calls default to an independent-gene working covariance. Supply sampling when gene dependence must be modelled.

Bayesian PoPS removes the target chromosome before optional marginal selection, nuisance projection, response/feature scaling and fitting. Without threshold, all features that vary in that training fold enter the Bayesian model. Supplied covariance blocks must stay within a chromosome. updateB, updateE and updatePi default to true and map to native effect-variance, residual-variance and mixture-weight learning. Defaults start with residual variance .5 and effect variance .5 divided by feature count and weighted active prior mass; variance learning uses four-degree scaled inverse-chi-square priors and weight learning uses unit Dirichlet concentrations. Advanced prior and control remain available. Priors use prepared training coordinates; predictions use training-centered features and original-unit posterior mean effects, excluding nuisance effects. This is a single-response workflow; it differs from qgg’s ordinary full-data pops() prediction rule. Inspect the compact native posterior and chain diagnostics in each fold’s fits before interpreting rankings.

The result is a feature-based priority score with genome-wide and within-chromosome ranks. It is not a p-value, causal probability or locus definition, and the function does not automatically choose a top gene per locus. Ridge fits retain the penalty, effective degrees of freedom, selected features and coefficients for every held-out chromosome. Bayesian fits retain posterior means, PIPs and sampler summaries for each training fold.

Joint Bayesian pathways

Omega <- matrix(c(1,.3,.3,1),2,dimnames=list(names(stat),names(stat)))
G <- diag(2);dimnames(G)<-dimnames(Omega)
fit <- gsea(stat, sets, "bayesr",
  sampling=list(gene_covariance="independent",trait_covariance=Omega),
  prior=list(effect_covariance=G),
  control=list(burnin=100,sampling_sweeps=1000,seeds=c(11,29)))
fit$pathways                 # trait-specific effects and PIPs
fit$joint                    # probability of activity in any trait
fit$parameters$joint         # prepared-scale parameter summaries
fit$diagnostics$joint        # conditioning, threads and chain PIP spread
plot(fit, trait="trait_A")
saveRDS(fit, "pathway-fit.rds")

Bayesian models fit all pathways jointly and adjust both response and features for the shared nuisance design. They support overlapping pathways and proper priors when the number of features exceeds the number of genes; highly dependent sets can still make individual PIPs difficult to interpret.

For one trait, supply prior=list(residual_variance=1,effect_variance=1) with values chosen for the prepared score scale. BayesC defaults to inclusion probability 0.5. Scalar BayesR defaults to weights (0.7,0.1,0.1,0.1) and variance multipliers (0,0.01,0.1,1). Optional scalar learning flags are estimate_effect_variance, estimate_residual_variance and estimate_weights. The corresponding priors are effect_variance_prior=list(df=...,scale=...), residual_variance_prior=list(df=...,scale=...) and positive weight_prior counts.

Joint fits require a labelled effect_covariance prior and a separate, fixed residual trait_covariance on original response units. The latter describes remaining gene-score dependence between traits, not pathway-effect covariance. All traits share genes, features, covariates and gene dependence. The native likelihood uses a separable gene-by-trait covariance model. It does not infer sample overlap or substitute a genetic-correlation matrix for residual covariance.

By default, joint models enumerate activity patterns across 2–12 traits. For larger trait sets, explicitly provide restricted binary patterns with labelled trait columns, distinct rows and the null row first; the native maximum is 64 traits. A restricted pattern set changes the model. pattern_weights default to equal probabilities. Joint BayesC has one positive component scale, 1; BayesR uses (0.01,0.1,1) with equal conditional component weights. The null is represented by the null pattern, not a zero component scale.

Joint learning is optional: estimate_covariance requires covariance_prior=list(df=...,scale=...) (an inverse-Wishart scale matrix); estimate_pattern_weights requires positive pattern_prior counts; and BayesR’s estimate_component_weights requires positive component_prior counts. The residual trait covariance stays fixed.

Results, units and bounds

• pathways retains method-specific statistics; summary() also returns parameters, diagnostics and provenance. plot() requires a trait choice when several are present. No automatic significance or PIP cutoff is applied.
• Bayesian pathway means and variances use original response/feature units. Priors, native chain traces and parameter summaries use prepared coordinates. preparation records response and feature scales and nuisance rank.
• scale_response and scale_features default to TRUE. Chain defaults are 100 burn-in sweeps, 1,000 sampling sweeps, thinning 1, seeds (1,2) and one thread. retain_diagnostics=TRUE retains parameter traces. Joint retain_pattern_probabilities=TRUE additionally stores pathway-by-pattern probabilities. Chain PIP spread is not an ESS or R-hat convergence check.
• Regression defaults bound features (255 classical / 4,096 Bayesian), user covariates (128), each covariance block (2,048 genes) and preparation payload (512 MiB). Bayesian preparation checks its native payload estimate; the R classical check bounds the feature matrix, not total process memory. Joint sampler storage has a separate 256 MiB memory_limit_bytes estimate.
• Competitive fits permit at most 256 coefficients by default and require reciprocal condition at least 1e-10; both controls are configurable. The bounded MAGMA route retains those native limits. No dense whole-genome covariance is formed when explicitly partitioned blocks are supplied.

The focused R checks compare ORA to the hypergeometric calculation, GLS to a dense R reference, and unadjusted MAGMA preprocessing to the corresponding competitive fit. Bayesian checks cover chains, units, labels, joint probabilities, learning and serialization. The GWAS adapter additionally checks quadratic moments and dependence against a direct R calculation, equal-eigenvalue chi-square tails, explicit block assumptions, unavailable tails and the connected pathway fits. These integration checks complement the small native oracles; they do not establish broad genomic calibration or full package qualification.