GlycoForge is a simulation tool for generating glycomic relative-abundance datasets with customizable biological group differences and controllable batch-effect injection.

Key Features

• Two simulation modes: Fully synthetic or templated (extract factor from input reference data + simulate batch effect)
• Paired multi-glycome simulation: simulate_paired() generates two glycomic datasets (e.g., N- and O-glycomics) from the same biological samples, with shared batch labels and optional controllable cross-class coupling
• Controllable effects injection: Systematic grid search over biological effect or batch effect strength parameters
• Motif-level effects: For both bio and batch effects, desired motif differences (e.g., Neu5Ac: down) can be introduced. These are propagated in a dynamically constructed biosynthetic network to ensure physiological glycomics data (e.g., corresponding increase in desialylated glycans in the example of Neu5Ac: down)
• Missing data simulation: MNAR (left-censored, biased toward low-abundance glycans), MCAR, MAR, and a realistic mixed composite, all calibrated on the glycomics datasets in glycowork
• Known motif-level truth: simulate_known_truth() injects effects on the carriers of chosen motifs of a real glycome, or along its real group difference, and states, for every motif, whether it truly changed

Quick Start

Installation

• Python 3.11–3.15 required (>=3.11,<3.16). We recommend creating a dedicated virtual environment:

python3.11 -m venv .venv
source .venv/bin/activate   # Windows: .venv\Scripts\activate

• Core dependency: glycowork>=1.12.0

pip install glycoforge

OR

git clone https://github.com/BojarLab/GlycoForge.git
cd GlycoForge
python3.11 -m venv .venv
source .venv/bin/activate
pip install -e .

Usage

See run_simulation.ipynb for interactive simulation examples, or use_cases/batch_correction/ for batch correction workflows, and benchmarking_batch_effect_removal.ipynb for the full six-method benchmark.

How the simulator works

All simulation operates in CLR (centered log-ratio) space using a Gaussian copula sampler:

• Reference construction: Zeros (the reference datasets' own missing values) are set to half of each glycan's minimum observed value before CLR, so the copula does not reproduce them as quasi-zero "clean" abundances. In synthetic mode, pool all class-matched glycowork datasets, compute per-dataset Ledoit-Wolf covariance, extract and average correlation matrices into a consensus R, and collect empirical CLR marginals from the pooled samples. In templated mode, compute a single Ledoit-Wolf covariance and extract empirical CLR marginals from the input dataset's samples.
• Clean data generation (simulate_clean_data): Draw Z ~ N(0, R_reg), apply the probability integral transform to get uniform marginals, then map to real CLR empirical quantiles via linear interpolation. This preserves both inter-feature correlation (Mantel r) and marginal realism (KS).
• Biological injection: In templated mode, inject real effect sizes in CLR space: z_U = z_H + m * lambda * d_robust, where m is the differential mask, lambda is bio_strength, and d_robust is the effect vector after robust_effect_size_processing. In synthetic mode, inject along the top-K eigenvectors of the pooled covariance with random signs, scaled by bio_strength * std * sqrt(n_glycans) to remain PVCA-meaningful across feature counts. The alpha_H and alpha_U parameters are computed to define the CLR-space injection direction (via clr(p_U) - clr(p_H))
• Batch effects: Two modes controlled by batch_mode. In additive mode (default, ComBat-correctable): Y_batch = Y_clean + kappa_mu * sigma * u_b + epsilon. In multiplicative mode (non-linear): Y_batch = Y_clean + kappa_mu * u_b * Y_clean + epsilon. Variance inflation uses batch-specific scale factors spread evenly around 1.0, controlled by var_b. Compositional pairing ensures substrate-product glycans receive correlated noise.
• Missingness (apply_missingness, chosen via missing_mechanism):
  • MNAR: logistic model in log-abundance space with per-sample intercept solved via Brent's method for exact target missingness fraction.
  • MCAR: p_ij = expit(c + s_i), independent of all values; s_i ~ N(0, sample_sd) models run-to-run variation (set sample_sd = 0 for textbook MCAR).
  • MAR: as MCAR, but s_i also rises with the share of the top_k most abundant glycans (never missing, so the driver is always observed) and with group/batch shifts.
  • mixed: shared detection limit (MNAR) with the MAR per-sample shifts, plus abundance-independent dropouts making up mcar_share of the missing values.
  • Defaults (missing_params: mcar_share = 0.13, sample_sd = 0.8, dominance = 0.35, top_k = 3, group_sd = 0.2, batch_sd = 0) are calibrated on the 18 glycowork glycomics datasets with sufficiently sized replicate groups, so that the same estimators run on templated GlycoForge data reproduce the real per-sample variation, dominance dependence, and between-group variation.

Simulation Modes

glycoforge.simulate() generates a single glycomic dataset per entity in two modes controlled by data_source. For paired multi-glycome data from the same samples, use glycoforge.simulate_paired() instead (see below). Configuration files are in sample_config/.

More than two groups (e.g., for benchmarking ANOVA and post-hoc tests): pass n_U as a list of group sizes. n_H=10, n_U=[10, 10, 10] simulates healthy plus three case groups (columns unhealthy1_*, unhealthy2_*, unhealthy3_*; bio_labels 0–3), each shifted from healthy by its own effect, with batches stratified within every group. group_effect_scale multiplies bio_strength per case group, so [1, 2, 0] gives a doubled effect and a null group with the same mean as healthy. In synthetic mode the first case group is drawn exactly as in the two-group setup and later ones along random rotations of the same top-eigenvector subspace (or, with motif_rules, from separate motif draws); in templated mode with use_real_effect_sizes=True, each case group gets the template's own effect sizes against its control group, taken from the glycowork contrasts in order (e.g. human_prostate_N_PMC8010466: control vs. stage1, stage2, ...) or from a list of prefixes in column_prefix['unhealthy'].

Synthetic mode (data_source="simulated") – Fully synthetic simulation (click to show detail introduction)

No real data dependency. Ideal for controlled experiments with known ground truth.

Pipeline steps:

1. Pools all class-matched glycowork datasets (filtered by glycan_class) to build a consensus Ledoit-Wolf correlation matrix and empirical CLR marginals via _build_copula_ref
2. Constructs alpha_H from the pooled mean abundance profile across all candidate datasets, scaled to sum proportional to 10 * n_glycans
3. For each random seed, generates alpha_U by randomly scaling alpha_H (30% up, 35% down by default), or via motif-guided compositional pairing if motif_rules is provided
4. Samples clean cohorts via Gaussian copula (LW correlation + empirical marginals) with biological injection along top-K eigenvectors of the pooled correlation matrix, each glycan shifted by bio_strength times its own SD (capped at 3 SD per direction)
5. Defines batch effect direction vectors u_dict once per simulation run (fixed seed ensures reproducible batch geometry across parameter sweep)
6. Applies batch effects controlled by kappa_mu (shift strength) and var_b (variance scaling)
7. Optionally applies missingness:
   • missing_fraction: proportion of missing values (0.0–1.0)
   • missing_mechanism: MNAR (default), MCAR, MAR or mixed; a list runs a grid over mechanisms
   • mnar_bias: intensity-dependent bias for MNAR/mixed (default 1.0, range 0.5–5.0)
   • missing_params: overrides for the calibrated mechanism defaults
8. Grid search over kappa_mu and var_b produces multiple datasets under identical batch effect structure

Key parameters: n_glycans, n_H, n_U, kappa_mu, var_b, batch_mode, missing_fraction, missing_mechanism, mnar_bias

Templated mode (data_source="real") – Extract biological effect from input reference data + simulate batch effect (click to show detail introduction)

Starts from real glycomics data to preserve biological signal structure. Accepts CSV file or glycowork.glycan_data datasets.

Pipeline steps:

1. Loads CSV and extracts healthy/unhealthy sample columns by prefix (configurable via column_prefix)
2. Runs CLR-based differential expression via glycowork.get_differential_expression to compute Cohen's d effect sizes
3. Reindexes effect sizes to match input glycan order (fills missing glycans with 0.0)
4. Applies differential_mask to select which glycans receive biological signal injection:
   • "All": inject into all glycans
   • "significant": only glycans marked significant by glycowork
   • "Top-N": top N glycans by absolute effect size (e.g., "Top-10")
5. Processes effect sizes through robust_effect_size_processing:
   • Centers effect sizes to remove global shift
   • Applies Winsorization to clip extreme outliers (auto-selects percentile 85–99, or uses winsorize_percentile)
   • Normalizes by baseline (baseline_method: median, MAD, or p75)
   • Returns normalized d_robust scaled by bio_strength
6. Injects effects in CLR space: z_U = z_H + mask * bio_strength * d_robust
7. Converts back to proportions: p_U = invclr(z_U)
8. Computes Ledoit-Wolf covariance from pooled CLR data
9. Samples clean cohorts via Gaussian copula using LW correlation and empirical CLR marginals, with the CLR injection vector shifting unhealthy samples
10. Defines batch effect direction vectors u_dict once per run (fixed seed ensures fair comparison across parameter combinations)
11. Applies batch effects: in additive mode, y_batch = y_clean + kappa_mu * sigma * u_b + epsilon; in multiplicative mode, y_batch = y_clean + kappa_mu * u_b * y_clean + epsilon, where variance inflation uses batch-specific scale factors controlled by var_b
12. Optionally applies missingness (same as Simplified mode)
13. Grid search over bio_strength, k_dir, variance_ratio, kappa_mu, var_b to systematically test biological signal and batch effect interactions

Key parameters: data_file, column_prefix, bio_strength, k_dir, variance_ratio, differential_mask, winsorize_percentile, baseline_method, kappa_mu, var_b, missing_fraction, missing_mechanism, mnar_bias, batch_mode

Paired mode (simulate_paired()) – Two glycomic classes from the same biological samples (click to show detail introduction)

Generates two glycomic datasets (e.g., N- and O-glycomics) that share sample identity: the same n_H healthy and n_U unhealthy individuals appear in both, so bio_labels, batch_labels, and column names are identical across glycomes. Each glycome is otherwise independently parameterised (different glycan counts, Dirichlet parameters, biological effect structures, and batch direction vectors).

Pipeline steps:

1. Draws independent log-normal healthy baselines alpha_H_A and alpha_H_B (fixed seeds 42/43 for reproducibility)
2. Generates per-glycome alpha_U using independent seeds so biological effects are not correlated by seed reuse
3. Simulates clean compositional data for both glycomes via Gaussian copula, using independently constructed pooled references from _build_copula_ref (class-filtered by glycan_class_A/B)
4. Optionally injects cross-class coupling in CLR space via shared latent factors Z ~ N(0, I):
   • Y_A_clr += coupling_strength * (Z @ U_A.T) * sigma_A
   • Y_B_clr += coupling_strength * (Z @ U_B.T) * sigma_B
   • At coupling_strength=0 the two glycomes are conditionally independent given sample labels; induced HSIC scales as coupling_strength²
   • Direction matrices U_A, U_B can be biased toward motif-matching glycans via coupling_motif_A/B
5. Round-trips through invclr to restore simplex validity after coupling injection
6. Applies shared batch labels with independent per-glycome direction vectors (same samples in the same batches, but different glycans affected)
7. Applies missingness independently per glycome (independent seeds prevent artificially correlated missing-value patterns; under MAR/mixed the shared bio and batch labels drive both)

Key parameters: n_glycans_A/B, bio_strength_A/B, k_dir_A/B, variance_ratio_A/B, coupling_strength, n_coupling_components, coupling_motif_A/B, kappa_mu, var_b, missing_fraction, missing_mechanism, mnar_bias

Known motif-level truth

glycoforge.simulate_known_truth() benchmarks motif-level tests on cohorts whose differentially abundant motifs are known in advance. A real glycome with its group signal removed is the template (motif_template(), which also holds glycowork's glycan x motif counts and containment DAG), effects are injected at the glycan level, and every motif's truth follows from the structures. These two functions need glycowork>=1.12.0; everything else in GlycoForge runs on 1.11.

from glycoforge import motif_template, simulate_known_truth
template = motif_template("human_serum_parkinson_GSL_PMID40379659")
sim = simulate_known_truth(template, scenario = "redistribution", n_H = 6, n_U = 6, log2fc = 1, missing_fraction = 0.2, seed = 0)
sim["table"], sim["truth"], sim["focus"], sim["children"]

Scenarios are null, inherited (a child's carriers move, and the parent only through them), own (carriers of a parent outside all its children move), redistribution (mass moves between two children, with the parent's total held), multi (the carriers of up to three DAG roots move), and cohort (the real difference between the template's two largest groups moves the glycans holding 90% of its energy, scaled so that its root mean square over them is log2fc). The first four are set-shaped, so a method built on carrier sets matches them by construction; cohort puts its signal on the noise axes as real cohorts do, and carrier_mix moves between the two. In the truth table, a motif is true if an injected glycan carries it and it changes by at least delta (0.2 log2 units, estimated from a large draw), grey if it carries the injection but changes less, closure if it changes without injected carriers, as closure moves every part of a composition, and null otherwise. Lists for n_U and log2fc give several case groups.

Use Cases

The use_cases/batch_correction/ directory demonstrates:

• Call glycoforge simulation, and then apply correction workflow
• Six-method batch correction benchmark (ComBat, Percentile, Ratio-ComBat, Harmony, limma-style, Stratified ComBat) across a parameter grid of biological signal strengths and batch effect severities
• Batch correction effectiveness metrics visualization

Limitation

More than two groups in simulate() only: simulate_paired(), simulate_circadian(), glycoforge_power() and glycoforge_reliability() still target the healthy/unhealthy setup.

Citation

If you use GlycoForge in your research, please cite:

Hu, S. and Bojar, D. (2026). GlycoForge generates realistic glycomics data under known ground truth for rigorous method benchmarking. bioRxiv, doi:10.64898/2026.02.20.707134

BibTeX:

@article{hu2026glycoforge,
  title   = {GlycoForge generates realistic glycomics data under known ground truth for rigorous method benchmarking},
  author  = {Hu, Siyu and Bojar, Daniel},
  journal = {bioRxiv},
  year    = {2026},
  doi     = {10.64898/2026.02.20.707134},
  url     = {https://www.biorxiv.org/content/10.64898/2026.02.20.707134v1}
}