CIDER — simulation benchmark

Code for the simulation study behind CIDER, a method that tests, for each transcription-factor (TF) → target-gene (TG) pair, whether the TF→TG relationship varies with a quantitative trait.

• Linear regime — OLS TG ~ TF + trait + TF:trait; test the TF:trait coefficient.
• Nonlinear regime — GAM (mgcv) TG ~ s(TF) + s(trait) + ti(TF, trait); test the ti() tensor-product interaction smooth. ti() excludes the marginals, so it isolates trait-dependent modulation.

Per-pair interaction p-values are BH-adjusted and pairs are ranked by interaction significance; performance is average precision against the known true positives.

This repository contains the simulation and benchmark code only. The application of CIDER to cohort data (LLFS, FHS, HRS, WHI) is not included: those data are not licensed for redistribution, and neither the data nor the cohort-specific analysis code is part of this deposit.

---

What is benchmarked

CIDER is compared against four differential-correlation methods, all of which dichotomise the trait at its median:

Method	Statistic
DGCA	differential correlation z-score / permutation p (Pearson)
DiffCorr	parametric Fisher-z differential correlation
discordant	mixture-model differential correlation posterior
DINGO	partial-correlation differential network (iDINGO)
CIDER_discrete	ablation — the CIDER linear model on the same binarised trait

The CIDER_discrete ablation separates what the model buys from what the continuous trait encoding buys.

Design

Factorial over TF-expression regime × interaction form × (sample size, effect size, noise) × 10 replicates. Internally the four combinations of regime and regime-type are called "stages":

	simulated TFs	real GTEx whole-blood TFs
linear interaction	stage 1	stage 2
nonlinear interactions	stage 3	stage 4

and two trait-separability versions: v1 = high separability (mixture means ±1.5), v2 = low (±0.5).

• 20 TFs, 80 TGs → 1,600 pairs per replicate, of which 10 are true positives (a 0.625% base rate; random-baseline AP ≈ 0.006).
• Grid: 10 log-spaced sample sizes {10 … 1000}, 10 effect sizes {0.1 … 10}, 10 noise levels {0.1 … 10} = 1,000 settings, 10 replicates each. Real-TF stages cap the top tier at 803, the number of GTEx whole-blood donors, which is the only difference between the two grid files.
• Core subset for presentation: 90 settings = 10 sample sizes × the central 3×3 of (β, σ). See sim/grids/simulation_grid_subset.txt.

Interaction forms

form	g(TF, trait)	shape
linear	TF·trait (+ main effects)	plane
tanh	tanh(3·TF)·trait	monotonic in TF; slope sign flips with trait
sinusoidal	sin(2·TF·trait)	oscillatory; the trait sets the frequency
stoichiometry	exp(−(TF−trait)²/2)	diagonal ridge, peaks where TF ≈ trait
synergy	σ(TF)·σ(trait)	AND-gate corner
quadratic	trait·(TF²−1)	symmetric parabola; concavity flips with trait sign

The quadratic form is the one that separates the methods: within each trait group the TF–TG relationship is a symmetric parabola, so the linear correlation is ≈ 0 in both groups and every correlation-based method is blind to it by construction. Definitions live in one place, sim/R/interaction_forms.R.

---

Layout

config.sh                  every path the pipeline uses; edit config.local.sh for your site
sim/
  R/
    simulate_cider_gam.R        nonlinear generator + evaluator (stages 3, 4), all options explicit
    interaction_forms.R         the six interaction forms + the background menu
  python/
    simulate_cider_linear.py    linear generator + evaluator (stages 1, 2)
    build_gtex_wholeblood_input.py   rebuilds the GTEx input from public sources
  benchmarks/                   benchmark_{DGCA,DiffCorr,discordant,DINGO,CIDER_discrete*}.R
  grids/                        parameter grids + cells.tsv (every cell and how to run it)
slurm/                     four array-job templates covering the whole pipeline
eval/                      aggregation, coverage audit, significance tests, diagnostics
figures/                   every published figure
tests/verification/        reads mgcv's basis and penalty objects out of fitted models
data/aggregated/           per-cell average-precision tables (gzipped) — regenerate figures without re-running the grid
docs/                      methods references, findings, and the consolidation notes

---

Running it

1. Configure

cp config.local.sh.example config.local.sh
$EDITOR config.local.sh          # CIDER_OUT, CIDER_R_LIBS, CIDER_PYTHON, CIDER_GTEX_DIR, ...
source config.sh
mkdir -p logs                    # the SLURM templates write here

config.sh sources config.local.sh first, so anything you set there wins over the defaults. Nothing outside config.local.sh contains an absolute path — the shipped example records the original cluster's settings for reference.

2. Dependencies

R 4.4.0 with mgcv, optparse, MASS, tidyverse, progress, yardstick, DGCA, DiffCorr, discordant, iDINGO, GGMridge, corrr, Biobase; Python 3 with numpy, pandas, scikit-learn, statsmodels, matplotlib, tqdm. See env/ for an installer and the exact recorded versions.

3. One cell, locally

Every option is explicit. sim/grids/cells.tsv records the exact option string behind each published cell, so the usual way to run one is to read it from there:

CELL=stage3_map_cider_quadratic
ARGS=$(awk -F'\t' -v c=$CELL '$1==c{print $9}' sim/grids/cells.tsv)

# generator (nonlinear, quadratic form, high trait separability)
Rscript sim/R/simulate_cider_gam.R $ARGS \
    --runid 0033 --samples 359 --effect_size 0.774 --noise 0.774 \
    --num_replicates 10 --output_dir $CIDER_OUT/runs/$CELL

# a benchmark consuming that generator's output
Rscript sim/benchmarks/benchmark_DGCA.R \
    --input_dir  $CIDER_OUT/runs/$CELL \
    --output_dir $CIDER_OUT/runs/stage3_map_DGCA_quadratic \
    --runid 0033 --num_replicates 10 \
    --tf_prefix TF --adjust perm --score raw_p

Generators must finish before their consumers start — the benchmarks read the per-replicate *_data.txt the generator writes, so that every method sees exactly the same simulated data.

4. The full grid, on SLURM

sim/grids/cells.tsv is the registry: all 256 cells with the runner, input cell, method arguments, grid file and replicate count for each. Every input_cell resolves to another row, so the generator/consumer dependencies are explicit. Drive a cell straight from it:

CELL=stage3_map_cider_quadratic
read -r RUNNER ARGS GRID <<< "$(awk -F'\t' -v c=$CELL '$1==c{print $7"\t"$9"\t"$10}' sim/grids/cells.tsv)"

sbatch --array=3-1001%190 \
       --export=ALL,CIDER_REPO=$PWD,CELL=$CELL,GRID=$GRID,GEN_ARGS="$ARGS" \
       slurm/run_cider_gam.sbatch

Array indices are 1-based line numbers in the grid file, whose line 1 is a header — runid N is on line N+2.

5. Aggregate and plot

python eval/aggregate_cell.py stage3_map_cider_quadratic   # -> $CIDER_OUT/aggregated/<cell>.txt
python eval/audit_coverage.py                              # completeness across all cells
python figures/plot_nonlinear_figures.py v1

To regenerate figures without re-running the grid, decompress the shipped tables into place:

mkdir -p $CIDER_OUT/aggregated
for f in data/aggregated/*.gz; do gunzip -c "$f" > "$CIDER_OUT/aggregated/$(basename "${f%.gz}")"; done

---

Reproducibility notes

• Defaults reproduce the paper, not best practice. Several behaviours in the original code look like oversights — a background menu whose fourth option is unreachable, GAM convergence failures being dropped rather than ranked last, a failed discordant replicate aborting a whole cell. They are the defaults because changing them changes the published numbers. Each has a flag that turns on the corrected behaviour, and new work should use it. See docs/CONSOLIDATION.md.
• Seeds. The published runs were unseeded. Both consolidated generators now accept --seed; supplying it changes nothing else.
• PYTHONHASHSEED matters. The linear-stage generators derive the background gene list with set(), whose iteration order over strings depends on Python's per-process hash seed — so two runs at the same numpy seed produce different data unless it is pinned. The SLURM template pins it. --deterministic_background removes the dependence, but changes which genes land in the background and so does not reproduce the published runs.
• Known non-convergence. discordant fails at high sample size combined with low noise (a mixture component collapses); those runs score AP = 0, with no survivorship credit. This is a property of the method on low-noise data, not a pipeline failure — see docs/DISCORDANT_CEILING.md.
• Cost. DINGO dominates: one dingo() call per replicate on the 100-gene network, which is why it has its own per-(runid, repid) grid and template.

---

Data

The GTEx whole-blood expression matrix used by the real-TF stages is not redistributed here. sim/python/build_gtex_wholeblood_input.py rebuilds it from the public GTEx v10, GENCODE v39 and Lambert et al. TF-list releases; see docs/DATA.md for the download locations and the exact filtering and normalisation applied.

---

Documentation

docs/simulation_study_reference.md	full design, parameters and findings
docs/simulation_design_methods.md	Methods-section reference for the design
docs/CONSOLIDATION.md	what was consolidated, the option axes, and why several defaults look wrong but are correct for reproducing the paper
docs/gam_smooth_specification.md	what the mgcv smooths actually are, verified against fitted objects
docs/benchmark_methods_code_check.md	cross-check of each benchmark's code against its published description
docs/DISCORDANT_CEILING.md	the discordant non-convergence ceiling
docs/DATA.md	how to rebuild the GTEx input

License

GPL-3.0-or-later. See LICENSE.