Robust Estimator of Gene-ENvironment Architectures
Whether a gene-by-environment (GxE) interaction is present can depend on the scale the phenotype is measured on. REGENA fits a variance-component model on a grid of monotone transformed phenotypes (Box-Cox as an example) in a single pass over the genotypes. At each scale it reports:
- the environment-specific genetic variances σ²_g0 and σ²_g1 (amplification when they differ)
- the cross-environment genetic correlation ρ_g (effect heterogeneity when ρ_g < 1)
- the environment-specific residual variances (NxE)
A companion script uses the scan to separate scale-dependent GxE, which some Box-Cox scale removes, from scale-independent GxE, which no scale removes. It also picks the scale to report.
REGENA is a standalone C++ program built from source on Linux or macOS.
| Dependency | Version | Needed for |
|---|---|---|
C++ compiler with C++11 support (g++ or clang++) |
tested: g++ 12.2.0, clang++ 14, Apple clang | building REGENA |
| CMake | ≥ 3.10 (tested: 4.3.2) | building REGENA |
| make | tested: GNU Make 4.3 | building REGENA |
| POSIX threads | provided by the system | multithreading |
| Eigen | 3.3.3, bundled in include/Eigen |
nothing to install |
Python 3 with numpy and pandas |
tested: Python 3.10 to 3.12, numpy 1.26.4 and 2.0.1, pandas 2.2.2 and 3.0.2 | scripts/regena_select_scale.py and the simulation scripts only |
scripts/extract_wide_table.py needs only the Python standard library.
| Operating system | Architecture | Compiler |
|---|---|---|
| Debian GNU/Linux 12 (bookworm) | x86-64 | g++ 12.2.0 |
| Ubuntu 22.04 (GitHub Actions) | x86-64 | g++ 12 and clang++ 14 |
| macOS 14 and macOS 15 (GitHub Actions) | arm64 | Apple clang (Xcode default) |
The Ubuntu and macOS rows are built and run on the demo at every push by the workflows in
.github/workflows.
No non-standard hardware is required. The demo runs on a laptop with a single core and
less than 100 MB of memory. Biobank-scale analyses (hundreds of thousands of individuals
and SNPs) are meant for a multi-core server or cluster node with tens of GB of memory;
-t sets the number of threads.
git clone https://github.com/sriramlab/REGENA.git
cd REGENA
mkdir build
cd build/
cmake ..
make
This produces the executable build/REGENA. There is nothing else to install: the Python
scripts in scripts/ are run in place.
Typical install time: about 1 minute on a normal desktop computer. The build itself
(cmake .. and make) took 42 seconds on a single 2.3 GHz core; make -j4 is faster.
The example directory holds a small, fully synthetic dataset: 5,000 individuals and
10,000 SNPs, with covariates, a binary environment, annotation files and two phenotypes.
No real genotype or phenotype data is involved. Every file is written from fixed seeds by
example/simulate_inputs.sh (independent SNPs in Hardy-Weinberg equilibrium with minor
allele frequency between 0.05 and 0.5), so the dataset can be regenerated byte for byte.
scan_demo.pheno is a trait with h² = 0.5 that is additive on the log scale and that the
environment multiplies (scripts/simulate_scale_gxe.py). On the observed scale this
becomes amplification GxE, which the log scale (λ = 0) removes.
After building REGENA:
cd example
bash test.sh
test.sh runs four REGENA jobs (a 13-scale Box-Cox scan, a single scale, a model
partitioned over two annotations, and the scan again from a config file), then writes the
per-λ table scan_demo.tsv and the scale selection scan_demo.scale.tsv. The selection
step is skipped with a message if numpy and pandas are not installed. Every job fixes
the seed of REGENA's random vectors (-s 1, or seed=1 in the config file), so the
output is the same on every run; without -s the seed is random and the estimates vary
slightly between runs.
Each REGENA job ends with REGENA ran successfully. The scan (scan_demo.tsv; tests from
scan_demo.trajectory.tsv) should read as follows, SEs in parentheses:
| λ | σ²_g, env0 (E = 1) | σ²_g, env1 (E = 0) | ρ_g | p, amplification | p, ρ_g = 1 |
|---|---|---|---|---|---|
| −0.5 | 0.437 (0.065) | 0.585 (0.052) | 0.73 (0.12) | 0.076 | 0.023 |
| −0.25 | 0.515 (0.078) | 0.552 (0.049) | 0.77 (0.13) | 0.69 | 0.081 |
| 0 | 0.586 (0.091) | 0.506 (0.046) | 0.81 (0.14) | 0.43 | 0.17 |
| 0.25 | 0.646 (0.105) | 0.451 (0.043) | 0.83 (0.15) | 0.086 | 0.26 |
| 0.5 | 0.691 (0.118) | 0.391 (0.040) | 0.84 (0.15) | 0.016 | 0.31 |
| 1 | 0.722 (0.145) | 0.268 (0.032) | 0.83 (0.16) | 0.0022 | 0.28 |
On the observed scale (λ = 1) the two genetic variances differ 2.7-fold (p = 0.002), the
factor exp(2 × 0.5) that was simulated. They agree near λ = 0 (p = 0.43), and ρ_g is not
distinguishable from 1 at either scale. scan_demo.scale.tsv therefore reports
classification scale-dependent
lambda_selected 0.25
p_amp_at_observed 0.00221276
p_rho_at_observed 0.282836
among its fields: three scales (λ = −0.25, 0 and 0.25) remove the GxE, and the rule in
Scale selection picks λ = 0.25 among them. The full expected files are
in example/expected/. With g++ on Linux the numbers match to the printed precision for
any number of threads. Another compiler or C++ standard library can draw different random
vectors from the same seed, which moves the estimates by a small fraction of their SE.
About 10 seconds for all of test.sh on a normal desktop computer (2 to 3 seconds per
REGENA job; measured on a single 2.3 GHz core).
To run REGENA on your own data, prepare the input files described in
File formats, run a Box-Cox scan as in Quick start with
one phenotype and one binary environment per job, and pass the output to
scripts/regena_select_scale.py (Scale selection). The phenotype must
be on its measured scale and strictly positive; add a constant shift first if it is not.
./REGENA -g <genotype prefix> -p <phenotype> -c <covariates> -e <binary env> -a <annotation> \
-m G+GxE+NxE -k 10 -jn 100 -t 8 \
-bx -bx-min -3 -bx-max 3 -bx-n 100 \
-o trait.env.out
python3 scripts/regena_select_scale.py --out trait.env.out --pheno <phenotype> --env <binary env>
Alternatively, run REGENA with a newline-separated config file. The config file uses full key names
(e.g. genotype=, use-boxcox=1) instead of flags; see example/config.txt.
./REGENA --config <config file>
-h, --help Show this message and exit
-g, --genotype The path of PLINK BED genotype file (prefix)
-p, --phenotype The path of phenotype file
-c, --covariate The path of covariate file
-a, --annot The path of input annotation for partitioned heritability
-o, --output Output file path
-e, --environment The path of the binary environment file
-m, --model Specification of the model. Currently there are 3 options:
1. additive genetic (G) effects only (arg: G)
2. additive genetic (G) and gene-environment (GxE) effects (arg: G+GxE)
3. additive genetic (G), gene-environment (GxE) and heterogeneous
noise (NxE) effects (arg: G+GxE+NxE)
-k, --num-vec The number of random vectors (10 is recommended)
-jn, --num-jack The number of jackknife blocks (100 is recommended)
-t, --nthreads The number of threads for multithreading
-s, --seed Seed for the random vectors (default: random)
-v, --verbose Verbose mode; output extra information (normal equations, sample sizes, ...)
-eXa, --eXannot Partition the GxE component with respect to the annotation file
(by default a single GxE component is fitted)
-np, --norm-proj-pheno Standardize the phenotype after regressing out covariates (default: 1)
-i, --cov-add-intercept Append an intercept to the covariates (default: 1)
Box-Cox scan
-bx, --use-boxcox Fit the model on a grid of Box-Cox transformed phenotypes
-bx-min, --boxcox-min Smallest lambda in the grid (default: -2)
-bx-max, --boxcox-max Largest lambda in the grid (default: 2)
-bx-n, --boxcox-n Number of evenly spaced lambda values (default: 100)
At each λ the phenotype is transformed with (y^λ − 1)/λ (log y at λ = 0), then
standardised to mean 0 and variance 1 across all non-missing individuals. This puts every
scale on one variance scale, so λ = 1 is the observed scale. -bx-n 1 with
-bx-min λ -bx-max λ fits a single scale.
Genotype: PLINK BED format (.bed/.bim/.fam). SNPs with MAF = 0 must be excluded.
Phenotype: header "FID IID name". With -bx the values must be strictly positive.
Covariate: header "FID IID name_of_cov_1 ... name_of_cov_n".
Environment: header "FID IID name", one column coded 0/1.
Annotation: M rows (M = number of SNPs) and K columns (K = number of annotations), no header.
If SNP i belongs to annotation j there is "1" in row i and column j, otherwise "0".
(delimiter is " ")
1) The number and order of individuals must be the same in the phenotype, genotype,
environment and covariate files.
2) The number and order of SNPs must be the same in the bim file and the annotation file.
3) Individuals with NA or -9 in the phenotype or environment are excluded.
4) Run one phenotype per job. With -bx the phenotype is expanded into one column per λ,
and every column shares its missingness.
REGENA writes the variance components, heritabilities and jackknife SEs for every fitted
scale to <output>. The same results appear in a tab-separated block between a WIDE_TABLE line
and a ***** line, one row per λ:
| column | meaning |
|---|---|
phenotype, lambda |
row index and its Box-Cox λ |
gamma_bin<b> |
genetic covariance between environments (annotation b) |
sigma2_g_env<e>_bin<b> |
genetic variance in environment e |
sigma2_e_env<e> |
residual variance in environment e |
rho_bin<b>_env0-1 |
genetic correlation between environments, γ / √(σ²_g0 σ²_g1) |
h2_g_env<e>_bin<b>, h2_e_env<e> |
heritability and residual fraction within environment e |
pval.rho, pval.sigma.diff, pval.h2.diff |
Wald tests of ρ_g = 1, σ²_g0 = σ²_g1 and h²_g0 = h²_g1 (annotation 0) |
Every estimate has a matching _se column. Environment labels: a 0/1 environment E is
expanded to the indicator columns [E, 1 − E]. So env0 is the stratum with E = 1 and env1 is
the stratum with E = 0.
scripts/extract_wide_table.py <output> prints this block as a TSV.
scripts/regena_select_scale.py reads the scan and the phenotype REGENA was run on:
- Valid region. Assume z_i(λ) = f_λ(y_i) ~ N(μ, σ²). The profile log-likelihood is ℓ_p(λ) = −(n/2) log σ̂²(λ) + (λ − 1) Σ log y_i. A λ is kept if (ℓ_p(λ) − max ℓ_p)/n ≥ −0.5 and it lies in the contiguous block of the grid around the maximum. This excludes transformations with implausible tail behaviour. Scales where either σ²_g is not positive are also dropped.
- GxE tests at every kept λ: amplification (H0: σ²_g0 = σ²_g1) and effect
heterogeneity (H0: ρ_g = 1). S is the set of λ where neither test rejects at
--alpha. - Scale-dependent (S non-empty): a Box-Cox scale removes the GxE. The reporting scale is the member of S that maximises (1 − w) σ²_g0 + w σ²_g1, where w is the mean of the environment column.
- Scale-independent (S empty): no scale removes the GxE. The optimal scale λ* is the
member of {λ : p_amp > 0.05 or p_rho > 0.05, min z(σ²_g) ≥ 2} closest to λ = 1, and
λ* = 1 if that set is empty.
remaining_gxenames the component still present at λ*. If σ²_g in one environment is not distinguishable from zero at λ = 1, S can be empty for reasons unrelated to scale;weak_sigma2_g_at_observedflags such pairs.
--alpha should be corrected for the number of trait × environment pairs tested.
--trajectory writes every per-λ statistic and test to a TSV.
To simulate phenotypes with GxE effects on real genotypes, see https://github.com/sriramlab/Simulator.
REGENA is released under the MIT License, an OSI-approved open-source license.
The bundled Eigen library (include/Eigen) keeps its own license, the Mozilla Public
License 2.0, with the exceptions noted in its file headers.
REGENA builds on GENIE and RHE-mc:
1. Ali Pazokitoroudi, Zhengtong Liu, Andrew Dahl, Noah Zaitlen, Saharon Rosset, Sriram Sankararaman.
AJHG (2024); doi: 10.1016/j.ajhg.2024.05.015
2. Ali Pazokitoroudi, Yue Wu, Kathryn S. Burch, Kangcheng Hou, Aaron Zhou, Bogdan Pasaniuc, Sriram Sankararaman.
Nature Communications (2020); doi: 10.1038/s41467-020-17576-9
v1.0.0