Quick start (end to end)

This walkthrough goes from raw files to mapped single-cell eQTL effects (main, interaction, or total), using only the command line. It runs on a small simulated demo dataset (download below) – 300 donors, 4,050 cells, and 3 genes (CTSS, ACTB, and TSPAN32) with 500 cis-variants each.

For your own study, the same four inputs are needed – just swap the file names:

Genotypes: a coordinate-sorted VCF with imputed dosages (DS) and/or genotype probabilities (GP) – here genotypes.vcf.gz.
Expression: a Market Exchange Format (MEX) count-matrix (as exported by Seurat/scanpy/10x pipelines) – here expr.mtx.gz + expr.features.gz + expr.barcodes.gz.
Metadata: one row per cell with the cell id, donor id, contexts (e.g. cell states), and covariates – here meta.tsv, with donor column donor_id, cell-id column cell_id, contexts cytotoxicity,treg_activation,central_memory, and normalized/scaled covariates scaled_age,sex,scaled_log_nUMI,percent_mito,gPC1-5,ePC1-5.
Gene coordinates: chromosome, start, end, and strand of the variants to test – here genes.bed (normally derived from GTF annotation).

Step 0: Get Julia, Dynema and demo data

Install Julia >= 1.9. Super easy!

Note

To use Dynema, no Julia coding is needed! Everything is command-line and simple!

Download Dynema:

curl -L https://github.com/immunogenomics/Dynema.jl/archive/refs/heads/main.tar.gz | tar -xz
mv Dynema.jl-main Dynema
cd Dynema

Download the demo dataset (~1 MB) from Dynema_datasets:

input=quick_start_data
mkdir -p "$input"

curl -L https://github.com/immunogenomics/Dynema_datasets/archive/refs/heads/main.tar.gz | \
  tar -xz -C "$input" --strip-components=3 "Dynema_datasets-main/data/quick_start"

Note

No need to worry about installing Dynema either! Once you use it on the command line, Dynema and all dependencies are automatically installed. This step happens only once and takes just a few minutes.

Step 1: Index the VCF (once per VCF)

dynema-map extracts each gene's variants to test (e.g. cis-window) from the VCF from a VCF file (BGZF-compressed and tabix-indexed). We provide the tool dynema-prepare-vcf in case your input file is not already BGZF-compressed and/or tabix-indexed.

dynema-prepare-vcf automatically BGZF-compresses (if required) and indexes the VCF file:

./bin/dynema-prepare-vcf --vcf "$input/genotypes.vcf.gz"

Your input file is never modified: any needed compression writes a new file alongside it, and the command reports the ready file to use. The demo VCF is plain gzip (not BGZF), so here that ready file is $input/genotypes_dynema.vcf.bgz – the file all the commands below pass to --vcf.

Note

This uses Dynema's bundled htslib (no system bgzip/tabix needed!) and handles plain .vcf, plain-gzip .vcf.gz, and already-BGZF files alike – see the prepare-vcf extra. If your VCF already ships with a .tbi index, skip this step and pass it to --vcf directly.

Step 2: Index the expression matrix (once per matrix)

Extracting the expression of a single gene normally requires scanning the whole Matrix Market file. Building a gene index .dgx once turns that into a milliseconds-scale read for every later run:

./bin/dynema-prepare-expr --expr-prefix "$input/expr"

Note

On a real study-scale matrix (thousands of genes, hundreds of thousands of cells, a multi-GB .mtx.gz) the benefit of indexing becomes obvious (each gene loading in ~0.5 s instead of minutes).

Step 3: Specify variants to test

A small bed-like file tells dynema-map what to map and where: one row per gene with chr, start, end, gene, strand. The TSS is derived FastQTL-style: the start position for + strand, and the end position for the- strand. The demo ships three genes as examples:

cat "$input/genes.bed"

chr1	150730706	150765778	CTSS	-
chr7	5527151	5530601	ACTB	-
chr11	2323216	2339430	TSPAN32	+

For your own study, write the file for a batch of genes you want to analyze. Preparing bed files as batches facilitates parallelization in HPC infrastructures.

Note

We also provide a tool to build bed files straight from your GTF annotation with dynema-prepare-bed – see the prepare-bed extra, facilitating parallelization.

Step 4: Map single-cell eQTL effects

Main effect: A convential eQTL test (context-independent). We are trying to answer the question: is there a constant/average effect of a variant on regulating gene expression?

The cell-state contexts enter as ordinary covariates here:

./bin/dynema-map \
  --expr-prefix "$input/expr" \
  --meta "$input/meta.tsv" \
  --vcf "$input/genotypes_dynema.vcf.bgz" \
  --bed "$input/genes.bed" \
  --covariates scaled_age,sex,scaled_log_nUMI,percent_mito,gPC1,gPC2,gPC3,gPC4,gPC5,ePC1,ePC2,ePC3,ePC4,ePC5,cytotoxicity,treg_activation,central_memory \
  --interaction-with cytotoxicity,treg_activation,central_memory \
  --donor-col donor_id --cell-id-col cell_id \
  --effect main \
  --out "$input/main"

Multi-context Interaction effect: This test leverages the heterogeneity in single-cell data to infer dynamic effect of variants on gene expression. We are trying to answer the question: does the variant's effect on gene expression change across the contexts (e.g. cytotoxicity,treg_activation,central_memory)?

The context's main effects are added to the model automatically via --interaction-with:

./bin/dynema-map \
  --expr-prefix "$input/expr" \
  --meta "$input/meta.tsv" \
  --vcf "$input/genotypes_dynema.vcf.bgz" \
  --bed "$input/genes.bed" \
  --covariates scaled_age,sex,scaled_log_nUMI,percent_mito,gPC1,gPC2,gPC3,gPC4,gPC5,ePC1,ePC2,ePC3,ePC4,ePC5 \
  --interaction-with cytotoxicity,treg_activation,central_memory \
  --donor-col donor_id --cell-id-col cell_id \
  --effect interaction \
  --out interaction

Note

Alongside the joint (multi-context) p-value p, Dynema automatically returns the single-context interaction p-values as one extra column per context (here p_cytotoxicity, p_treg_activation, p_central_memory) – computed from the same model fit, at no extra cost. These show which context(s) drive a joint signal; the joint p remains the primary test (pass --per-context false to omit them).

Total effect: this test captures both main effect and its interactions jointly. We are trying to answer the question: is there any genetic effect, constant or context-dependent regulating gene expression?:

./bin/dynema-map \
  --expr-prefix "$input/expr" \
  --meta "$input/meta.tsv" \
  --vcf "$input/genotypes_dynema.vcf.bgz" \
  --bed "$input/genes.bed" \
  --covariates scaled_age,sex,scaled_log_nUMI,percent_mito,gPC1,gPC2,gPC3,gPC4,gPC5,ePC1,ePC2,ePC3,ePC4,ePC5 \
  --interaction-with cytotoxicity,treg_activation,central_memory \
  --donor-col donor_id --cell-id-col cell_id \
  --effect total \
  --out "$input/total"

Note

You can also run multiple tests with a single command (e.g. --effect main,interaction). Each gene's data is then loaded and extracted only once, shared across the tests. The effect name is appended to each output file – <out>_<effect>_<gene>.tsv, e.g. results_main_CTSS.tsv and results_interaction_CTSS.tsv with --out results – and the summary (one row per gene × effect) stays in a single <out>_summary.tsv. Note that the main test in this mode omits the untested interaction terms from its model, exactly matching the standalone command above.

Step 6: Read the results

Each command-line run writes one file per gene plus a summary:

<out>_<gene>.tsv (e.g. main_CTSS.tsv; <out>_<effect>_<gene>.tsv in a multi-effect run) – the per-variant results, one row per cis variant.
<out>_summary.tsv – the batch summary, one row per gene × effect with its lead variant.
<out>.log – the full log messages.

The per-gene table (`<out>_<gene>.tsv`)

One row per cis variant, with these columns:

variant – variant identifier, as given in the VCF.
chr, pos – the variant's chromosome and position, copied from the VCF.
z, t, χ² or F – the cluster-robust score statistic. Tests of a single term (--effect main, or interaction with one context) report a signed z. Joint tests (interaction with 2+ contexts, total) report a χ² with as many degrees of freedom as tested terms. With --small (see below) these become t and F respectively.
p – the analytical p-value of the score test with cluster-robust (donor-level) variance: two-sided normal for z, upper-tail χ² for joint tests; with --small, two-sided Student's t with G − 1 degrees of freedom for t and an upper-tail F with (q, G − 1) degrees of freedom for F (G = number of donors, q = number of tested terms).
p_<context> (e.g. p_CV1, p_CV2, ...) – only for multi-context interaction tests (on by default; --per-context false to omit): per-context 1-df p-values decomposing the joint statistic, computed from the same null fit and cluster-robust covariance. They show which context(s) drive a joint signal.
p_boot, p_boot_approx – only with --boot: the wild score bootstrap p-value and its beta-distribution-smoothed version (useful when the true p is below the resolution of the bootstrap replicates).
Effect-size columns – one column per model term, named as in the fitted formula: (Intercept), G (the genotype main effect), one column per covariate and context, and G & <context> for each interaction term. These come from refitting the full unrestricted GLM at the variant, so they are exact maximum-likelihood estimates, not score-test approximations. The model is Poisson with a log link, so estimates are on the natural-log scale: G is the log fold change in expression per copy of the counted allele, and G & <context> is the change in that per-allele log fold change per unit of the context variable. With --betas lead (the default) only the lead variant's row – smallest analytical p, including exact ties – is filled and all other rows are missing; --betas all fits every variant (slower); --betas none omits the columns entirely.
Standard-error columns – one se_wald_<term> column per genotype effect-size column (se_wald_G, se_wald_G & <context>; covariate and context main-effect coefficients carry no SE), filled for exactly the same rows. Each is the one-way cluster-robust (CRVE, donor-clustered) standard error of the corresponding coefficient, computed as a sandwich estimator at the unrestricted fit, on the same natural-log scale as the estimate. So estimate / se is a Wald z-statistic (a t on G − 1 df with --small), and estimate ± 1.96 × se an approximate 95% confidence interval (use the t quantile with --small). Two things to keep in mind:
- These are not the quantities behind z/t and p. Those come from a score test at the null fit, so |estimate| / se will generally differ from |z| (closely for small effects; more for large effects, rare variants, or few donors). The score-test p (or p_boot) remains the primary inference.
- Without --small no small-sample correction is applied (consistent with the score test), so SEs are optimistic with few donors. With --small the variance carries the same G/(G−1) factor as the score test. Under multiway clustering (Julia API only) the SE columns are NaN.

Inference vs. effect size: which columns to use

Use z (or χ²) and p for inference – deciding whether an association exists. They come from the score test evaluated at the null fit with donor-clustered variance, which is the statistic whose calibration Dynema guarantees, and it is computed for every variant.

Use G, G & <context> and their se_wald_* columns for describing the effect – its size, direction and uncertainty in tables, confidence intervals and forest plots. These come from refitting the full model at the lead variant(s), so the estimate is the exact maximum-likelihood value and the SE is the cluster-robust sandwich evaluated at β̂.

The two are asymptotically equivalent under the null but will not agree exactly: G / se_wald_G is a Wald statistic and differs from the score z, increasingly so for larger effects. Do not derive a second p-value from the Wald ratio – report the score-test p and the Wald-type G ± se_wald_G. Dynema deliberately writes no first-order (one-step, null-fit) effect-size approximation: every estimate and SE in the output comes from the exact unrestricted refit.

Small-sample correction (`--small`)

By default Dynema uses the asymptotic reference distributions (normal / χ²) with the plain cluster-robust variance. With few donors this is anti-conservative, so dynema-map --small (default: off) switches all analytical and bootstrap results to the small-sample versions, mirroring WildBootTests.jl's small = true:

	single-term test (1 df)	joint test (q df)
default	`z`, normal	`χ²`, χ²(q)
`--small`	`t`, Student's t (G − 1 df)	`F`, F(q, G − 1)

With --small, the cluster-robust variance is inflated by G/(G − 1) (G = number of donors), the statistic is rescaled by a factor n/(n − 1) (n = number of cells; essentially 1 for single-cell data), the score bootstrap uses the same corrected statistic, and every standard error (se_wald_<term>) carries the same factor. Because bootstrap p-values compare the observed statistic with a bootstrap distribution rescaled identically, the empirical p_boot is unchanged by the rescaling itself, while p_boot_approx maps statistics to nominal p-values with the t/F reference.

The batch summary (`<out>_summary.tsv`)

One row per gene × effect, concatenated across chunks, this is the table to scan for a study's top associations:

gene, effect – the gene label and which test the row summarizes (main, interaction, or total).
status – ok, skipped (output already existed under --skip-existing; lead statistics are re-read from the existing file when possible), or failed (one failed gene × effect does not stop the rest of the batch; check the log for the error).
chr, tss – the gene's chromosome and TSS from the bed file.
n_variants – number of cis variants tested.
lead_variant, lead_pos – the variant with the smallest analytical p for this gene × effect.
stat_type, lead_stat – which statistic the effect uses (z, t, χ² or F; the latter two joint) and its value at the lead variant.
lead_p – the lead variant's analytical p-value.
lead_p_boot, lead_p_boot_approx – the lead variant's bootstrap p-values (empty unless the run used --boot).
lead_G, lead_se_wald_G – the lead variant's genotype main-effect estimate and its Wald CRVE standard error (from the per-variant table's G / se_wald_G). For interaction tests, also lead_G & <context> and lead_se_wald_G & <context> for every tested context. Columns that do not apply to an effect (e.g. interaction terms on a main row) are empty, as are all of them with --betas none.
out_file – path to the gene's per-variant table, for drill-down.

Because the data are simulated with known effects, the results should tell a clean story: TSPAN32 has a hit near its TSS in all three runs (a main effect plus a G × cytotoxicity interaction); CTSS has no main effect but has considerable interaction and total effects (its simulated effect exists only through treg_activation); ACTB stays null everywhere.

The dataset included in this demo is simulated to demonstrate Dynema's application. It does not represent a real biological scenario.

Where to next

Map thousands of genes by chunking the bed file and submitting one job per chunk – see batching, including --check and --skip-existing for safe HPC runs.
Add bootstrap p-values for small cohorts with --boot.
Every option: ./bin/dynema-map --help and the CLI overview.