Skip to content

Latest commit

 

History

History
194 lines (141 loc) · 9.36 KB

File metadata and controls

194 lines (141 loc) · 9.36 KB

Simulation design — Methods reference

Reference for writing the "Simulation design" Methods subsection. Every formula and parameter below is taken directly from the generator scripts (scripts/simulation/simulate_stage1_map.py, simulate_stage1_map_v2.py, and simulate_stage3_map_cider_<form>.R / simulate_stage4_map_cider_<form>.R), verified 2026-07-17. Notation: X = a TF's expression, T = the quantitative trait, β = effect_size, σ = noise_std, ε ~ N(0, σ²).


1. Overview: a 2×2 stage design

The benchmark crosses two factors, giving four stages:

Stage TF expression TF–TG relationship CIDER model
1 Simulated (independent Gaussian) Linear interaction CIDER-linear (OLS)
2 Real (GTEx whole blood) Linear interaction CIDER-linear (OLS)
3 Simulated (independent Gaussian) Nonlinear interaction (5 forms) CIDER-GAM
4 Real (GTEx whole blood) Nonlinear interaction (5 forms) CIDER-GAM

Stages 1/2 test detection of a linear TF×trait interaction; stages 3/4 test nonlinear interaction surfaces invisible to correlation-based methods. Simulated-TF stages (1,3) isolate the statistical model from real-covariate structure; real-TF stages (2,4) confirm the results under empirical whole-blood expression.

Each stage is additionally run at two trait-separability settings, v1 and v2 (§3), which control how well the dichotomized trait (used by the benchmark methods) recovers the latent structure.

Fixed dimensions in every stage: 20 TFs, 80 target genes (TGs), 10 true-positive TF–TG pairs, giving 1600 candidate TF–TG pairs of which 10 are positive.


2. TF expression

2.1 Simulated TFs (stages 1, 3)

The 20 TFs are drawn independently from a standard multivariate normal:

$$\mathbf{X}_i \sim \mathcal{N}(\mathbf{0}, \mathbf{I}_{20}), \qquad i = 1,\dots,n,$$

i.e. zero mean and identity covariance (mutually independent TFs, unit variance). This makes each TF marginally N(0,1), so E[X²]=1 (used by the quadratic form's centering, §6).

2.2 Real TFs (stages 2, 4)

TFs are sampled from GTEx v10 whole-blood expression (gtex_v10_wholeblood_3cpm_25perc_log1tmmcpm.tsv; ≥3 CPM in ≥25% of samples, log1p(TMM-CPM)). For each replicate, simulate_tf_expression:

  1. restricts the gene annotation to the 834 transcription factors (is TF == True, from the Lambert TF list) and draws 20 TFs at random without replacement;
  2. draws n donors at random (without replacement) from the 803 available;
  3. transposes to a samples × TF matrix and standardizes each TF (scale(), z-score).

Because the draw of TFs, donors, and (below) true-positive assignments is repeated inside the replicate loop with no fixed seed, every replicate uses a fresh random 20-TF panel; TF-panel variability is therefore already reflected in replicate-to-replicate variation. (Real GTEx TFs are non-Gaussian, so after standardization E[X²] ≈ 1 only approximately — see the quadratic-form caveat in §6.)


3. Quantitative trait

The continuous trait is drawn from a balanced two-component Gaussian mixture:

$$T \sim \tfrac{1}{2},\mathcal{N}(-\mu,,1) ;+; \tfrac{1}{2},\mathcal{N}(+\mu,,1),$$

with equal mixing weights (0.5/0.5), unit component variance, and symmetric means ±μ (gmm_mean). The two components represent latent subpopulations; μ sets their separability:

  • v1 (high separability): μ = 1.5
  • v2 (low separability): μ = 0.5

The trait is symmetric about 0 and bimodal for large μ. It is used untransformed by CIDER; the benchmark methods use the dichotomized version (§4).


4. Trait dichotomization for the benchmark methods

The four benchmark methods (DiffCorr, DGCA, Discordant, DINGO) require two discrete groups, so the trait is split at its sample median:

$$g_i = \begin{cases} 2 & T_i > \operatorname{median}(T)\ 1 & \text{otherwise,}\end{cases}$$

giving approximately balanced groups (n₁ ≈ n₂). Because the median split is not the same as the latent GMM component assignment — the two diverge most under low separability (v2) — it represents what an analyst without oracle knowledge of the components would do, and is the fair input for the comparators. CIDER uses the continuous trait throughout.


5. Target gene expression — linear stages (1, 2)

TGs are generated as true positives or background. Ten TGs are chosen as true positives and paired 1:1 with ten randomly chosen TFs; background TFs are drawn with replacement.

True-positive TG (a genuine linear TF×trait interaction, with matched main effects):

$$\mathrm{TG} = \beta X + \beta T + \beta,(X \cdot T) + \varepsilon,$$

where all three terms share the coefficient β. The defining signal is the multiplicative X·T interaction; CIDER-linear tests the OLS TF:Trait term.

Background TG (no interaction) — one structure chosen uniformly at random per gene:

$$\mathrm{TG} \in {, \varepsilon,;; \beta X + \varepsilon,;; \beta T + \varepsilon,;; \beta X + \beta T + \varepsilon ,}.$$

Background genes may thus carry TF and/or trait main effects but never an interaction, so methods must distinguish interaction from marginal association. (In the R nonlinear stages the background set omits the additive-both option; see §6.)


6. Target gene expression — nonlinear stages (3, 4)

True-positive TGs are generated from a pure trait-modulated nonlinear interaction (no additive main-effect terms): TG = β · f(X, T) + ε. Five interaction surfaces f(X,T) are used (one form per run); each is designed so the TF→TG relationship's shape depends on the trait.

Form f(X, T) Note
tanh tanh(3X) · T monotone-saturating response, trait-scaled
sinusoidal sin(2·X·T) oscillatory tensor-style interaction
stoichiometry exp(−(X − T)² / 2) Gaussian "matching" bump (peaks when X ≈ T)
synergy sigmoid(X) · sigmoid(T) AND-gate; sigmoid(u)=1/(1+e^{−u})
quadratic T · (X² − 1) trait flips concavity; zero within-group linear correlation

The quadratic form is the key discriminating case: with standardized TF (E[X²]=1), the −1 centering makes it a pure interaction with no marginal quadratic main effect, and the within-group Pearson correlation is ≈0 in both trait halves — so correlation-based benchmarks are blind to it by construction while a GAM detects the concave↔convex flip. (For real TFs in stage 4, non-Gaussianity means X²−1 is not perfectly mean-zero, so a small amount of linear signal leaks and correlation methods creep up slightly — report honestly.)

Background TG (no interaction) — one of three structures chosen uniformly at random:

$$\mathrm{TG} \in {, \varepsilon,;; \beta X + \varepsilon,;; \beta T + \varepsilon ,}.$$

CIDER-GAM fits TG ~ s(TF, k=5) + s(Trait, k=5) + ti(TF, Trait, k=5) (REML) and tests the ti() interaction smooth.

Code note for the writer: in simulate_stage3_map_cider_{sinusoidal,stoichiometry}.R the descriptive writeLines(...) strings mention additive b1*TF + b2*trait terms, but the executed y <- ... expression contains only the interaction term + noise (verified). Describe the code as run (pure interaction), not the stale comment.


7. Simulation grid

Each (stage, form, separability) is evaluated over a grid of sample size × effect size (β) × noise (σ), with replicate simulations per cell.

  • Sample sizes (10 values, log-spaced, round(10^linspace(1,3,10))): 10, 16, 27, 46, 77, 129, 215, 359, 599, 1000. Real-TF stages (2, 4) cap the top tier at 803 (the number of GTEx whole-blood donors).
  • Effect size β and noise σ: each varied over 10 log-spaced values in [0.1, 10] for the full grid (10×10×10 = 1000 configurations).
  • Core subset used for the reported benchmark (simulation_grid_subset.txt): the central 3 β × 3 σ values — β, σ ∈ {0.464, 0.774, 1.292} — crossed with all 10 sample sizes = 90 cells per (stage, form, separability). This spans the informative regime while keeping compute tractable.
  • Replicates: 10 per cell originally; 50 in the bootstrap-CI re-run (_r50), each a fresh independent simulation (new TFs/trait/assignments).

Performance metric: average precision (AP) of the 1600 pairs against the 10-positive ground truth (rank-based; see the benchmark-methods doc for per-method score columns); MAP = mean AP across replicates.


8. Parameter summary

Parameter Value Source
# TFs 20 n_tfs
# TGs 80 n_tgs
# true positives 10 n_true_pos
candidate pairs 1600 (10 positive) 20 × 80
simulated TF N(0, I₂₀) simulate_tf_expression (py)
real TF 20 of 834 Lambert TFs, n of 803 donors, z-scored simulate_tf_expression (R)
trait ½N(−μ,1) + ½N(+μ,1) simulate_trait
separability μ 1.5 (v1) / 0.5 (v2) gmm_mean
dichotomization median split → groups 1/2 benchmark scripts
linear TP βX + βT + βXT + ε stage 1/2
nonlinear TP β·f(X,T) + ε, 5 forms stage 3/4
noise ε ~ N(0, σ²) noise_std
sample sizes 10…1000 (803 cap for real TF) grid
β, σ grid {0.464, 0.774, 1.292} core; [0.1,10] full grid
replicates 10 (50 in _r50) driver