API Reference

map_locus(f::FormulaTerm; pheno::AbstractVector, geno::Union{AbstractDataFrame, AbstractVector},
          meta::AbstractDataFrame, groups::Union{AbstractDataFrame, AbstractVector},
          termtest::Union{String, Vector{String}}) -> Dynema_struct

Performs a single-cell eQTL mappping for a gene. All data arguments are defined at the single-cell level observation.

Mandatory arguments

• f: Modelling formula focused on predictors. Should include all the terms to be tested and covariates. As genotype information

can change when iterating across variants, G is reserved to represent genotype information for each genetic variant. For example, a minimal formula can be defined as 0 ~ 1 + G, where 1 is the intercept term and G represents the genetic variant being tested. Gene expression can be represented with no name simply as 0.

• pheno: Gene expression counts for a particular gene. A vector of length n, where n is the number of cells
• geno: Genotype information in allele dosage format or genotype probabilities (e.g. [0, 2, 1, 0, 1])

expanded at the single-cell level (length n). If a DataFrame is provided, each column should correspond to an individual genetic variant and all variants will be tested. The dimensions of geno should be (n x gv), where gv is the total number of variants

• meta: Single-cell level metadata including single-cell cell states of interest. It also should include

single-cell and donor-level covariates expanded at the single-cell level. Dimensions n x p, where p is the number of variables

• groups: A Vector or DataFrame indicating how cells should be grouped. At least one grouping variable should be

specified (e.g. donor structure). For example ["donor1", "donor2", "donor1", "donor3"].

• termtest: Term included in the formula that is tested. In the case of dynamic effects, this should be the interaction

G and a single-cell covariate (e.g. "G & CV1")

Optional arguments

• parallel: Distributes contiguous blocks of variants across worker processes

(Distributed.jl; start workers first, e.g. julia -p 4). Each worker builds its own locus workspace, so the fast path applies in parallel too. Warm-start chains restart at block boundaries, so interaction-test results are reproducible for a fixed worker count but can differ at IRLS convergence tolerance (~1e-7) across worker counts or vs. a serial run. Without parallel, if the process has more than one Julia thread (start with julia -t N), variants are threaded across contiguous blocks instead – same block semantics and reproducibility caveat, lower overhead than workers (no serialization, one process). Threading is not used when boot = true (the bootstrap library manages its own threading)

• H0: Null hypothesis value. By default 0

Testing always uses a cluster-robust (CRVE) score/Lagrange-multiplier test at the restricted (null) fit. If the null model contains no genotype-dependent terms (as in main/total effect tests, where every tested term involves G), it is fit only once and reused across all variants.

• betas: Which variants get unrestricted-model coefficient estimates attached as extra summary columns (they are not

needed for the score test itself): :lead (default) fits only the lead variant(s) – smallest analytical p-value, including exact ties – leaving the estimate columns missing for all other variants; :all fits every variant (slowest); :none skips estimates entirely (fastest). Each genotype estimate column (G and every G & context) is accompanied by an se_wald_<term> column: the one-way cluster-robust (CRVE) standard error of that coefficient at the unrestricted fit (see crve_se); covariate and context main-effect coefficients are written without SEs. Estimates and SEs are on the natural-log scale and give a Wald z = estimate / SE (t on G - 1 df with small = true). They are distinct from the score-test z/t and p, which are computed at the null fit, so |estimate| / SE generally differs from |z|; the score-test p (or p_boot) remains the primary inference. SEs carry a small-sample correction only with small = true, and are NaN under multiway clustering. Covariate/context names in the formula must not coincide with the summary columns Dynema writes (variant, chr, pos, z, t, χ², F, p, p_boot, p_boot_approx, p_<context>, or the se_wald_<term> of another term): map_locus stops with an error naming the offending column(s) – rename them in the metadata

• boot: Apply score bootstrapping? If false, analytical p-values using a CRVE are returned. Bootstrapping

requires the optional WildBootTests package (Pkg.add("WildBootTests"); using WildBootTests alongside using Dynema)

• B`: Number of bootstrap iterations to apply. By default 39999 iterations at apply to achieve a p-value of 5 x 10^-5 for a two-tail test.

For adaptive bootstrapping, the number of iterations might be specified as a vector indicating the number of iterations to perform in each step

• ptype: Type of bootstrap p-value to return (:equaltail, :symmetric). Only applicable when one term is being tested.
• rboot: Whether to return the bootstrap distributions for each variant. Useful for direct assessment of bootstrap statistics and

internal debugging

• pos: A numeric value specifying a genomic location for each genetic variant. Stored in final output for convenience
• gene: Name of the gene being tested. Stored in final output for convenience
• chr: Chromosome position of gene being tested. Stored in final output for convenience.
• percontext: For multi-context interaction tests (per-variant null with more than one tested term, one-way

clustering), additionally report one 1-df p-value column per tested context (p_<context>), computed from the same restricted fit and CRVE covariance as the joint test (see crve_percomponent) – no extra model fitting. They decompose the joint statistic; exact under the joint null, but with correlated contexts a true effect in one context can partially project onto another's test, so the joint p-value remains the primary inference. Ignored for main/total effect tests (shared null) and single-context interactions. Default true

• small: Apply the small-sample correction of WildBootTests.jl's small = true to all analytical and bootstrap

results. Default false (normal/χ² reference, no correction). With small = true: the cluster-robust denominator is multiplied by G/(G-1) (G = number of donors/clusters), the statistic is rescaled by the nobs/(nobs-1) factor (see smallsample_factors), and the reference distribution becomes Student's t with G-1 df for a single tested term (summary column t instead of z) or F with (q, G-1) df for a joint test of q terms (column F instead of χ²). The bootstrap uses the same corrected statistic (p_boot is unchanged by this rescaling, being rank-based; p_boot_approx uses the t/F reference), and all standard errors below are on the same footing. Multiway clustering is handled by the WildBootTests implementation (needs the optional package) and has no score-based or Wald SEs (NaN).

Standard errors

Inference comes from the score test (z/t/χ²/F and p, evaluated at the restricted fit). Effect sizes come only from the unrestricted refit selected with betas: the genotype coefficient columns (G, G & context) are accompanied by se_wald_<term> columns, the Wald-type CRVE standard error at the unrestricted fit (see crve_se); covariate coefficients are written without SEs. G / se_wald_G is a Wald statistic and will not reproduce the score-test z; report the score-test p for inference and G ± se_wald_G for the effect. No one-step (null-fit) estimate or SE is written. Joint (multi-term) tests have no single SE; their per-context p-values are p_<context>.

expand_geno(mat::AbstractMatrix, rows::AbstractVector;
             rownames=nothing, colnames=nothing)

Memory-efficient wrapper for genotype matrices that allows lazy row and column subsetting. Instead of building a full matrix at the single-cell level, it provides a lazy genotype matrix representation compatible with downstream operations. Genotype data at the single-cell level is only materialized when subsetting a single column or row.

Arguments

• mat::AbstractMatrix: donor × SNP genotype matrix (e.g., Matrix{Float64} or SparseMatrixCSC).
• rows::AbstractVector: mapping from single-cell to donor rows.
• rownames::Union{Nothing,String,AbstractVector}=nothing: optional cell ids
• colnames::Union{Nothing,String,AbstractVector}=nothing: optional SNP names.

Indexing

• 2D indexing only: E[i, j], E[i, :], E[:, j], E[I, J].
• Single-index access is forbidden: E[i] throws an error.
• Materializes single rows or single columns only.
• Submatrices return ExpandedGenoView (lazy, no copying).

Examples

# Construct an ExpandedGeno
E = expand_geno(geno, indices;
                 rownames = cell_id,
                 colnames = snp_names)

# Access a single element
val = E[10, 5]        # Float64

# Access a single row
row = E[10, :]         # Vector{Float64}

# Access a single column
col = E[:, 5]          # Vector{Float64}

# Lazy submatrix (returns ExpandedGenoView)
view = E[1:100, 5:10]

# Access row and column names
rnames = names(E, 1)
cnames = names(E, 2)

expand_geno(geno::AbstractMatrix, geno_donor_ids::AbstractVector;
             sc_donor_ids=nothing, variant_names=nothing,
             cell_names=nothing)

Memory-efficient wrapper for genotype matrices that allows lazy row and column subsetting. Instead of building a full matrix at the single-cell level, it provides a lazy genotype matrix representation compatible with downstream operations. Genotype data at the single-cell level is only materialized when subsetting a single column or row.

Arguments

• geno::AbstractMatrix: donor × variant matrix containing allele dosages (e.g., Matrix{Float64} or SparseMatrixCSC).
• geno_donor_ids::AbstractVector: A vector of donor ids, in the same order as the rows in the genotype matrix.
• sc_donor_ids::AbstractVector: A vector of donor ids for each cell/barcode. This vector should be in the same order as the metadata and gene expression at the single-cell level.
• variant_names::Union{Nothing,String,AbstractVector}=nothing: A vector of optional variant names. Strongly encouraged.
• cell_names::Union{Nothing,String,AbstractVector}=nothing: (Optional) A vector of barcode/cell names.

Indexing

• 2D indexing only: E[i, j], E[i, :], E[:, j], E[I, J].
• Single-index access is forbidden: E[i] throws an error.
• Materializes single rows or single columns only.
• Submatrices return ExpandedGenoView (lazy, no copying).

Examples

# Construct an ExpandedGeno
ex_geno = expand_genotypes(Matrix(geno), donor_ids, sc_meta.donor_id, variant_names)


# Access a single element
val = ex_geno[10, 5]        # Float64

# Access a single row
row = ex_geno[10, :]         # Vector{Float64}

# Access a single column
col = ex_geno[:, 5]          # Vector{Float64}

# Lazy submatrix (returns ExpandedGenoView)
view = ex_geno[1:100, 5:10]

# Access row and column names
rnames = names(ex_geno, 1)
cnames = names(ex_geno, 2)

get_f(::Dynema.DynemaModel)

Extract formula used for a DynemaModel

get_termtest(::Dynema.DynemaModel)

Extract formula used for a DynemaModel

get_ncell(::Dynema.DynemaModel)

Extract number of cells used for a DynemaModel

get_ndonor(::DynemaModel)

Extract number of donors/individuals for a DynemaModel

get_summary(::Dynema.DynemaModel)

Extract all summary statistics for a DynemaModel

get_stat(::Dynema.DynemaModel)

Extract the test statistic for every variant (the summary column named by get_stattype: z (t with small = true) for a single tested term, χ² (F with small = true) for a joint test)

get_p(::Dynema.DynemaModel)

Extract the analytical CRVE p-values for a DynemaModel (bootstrap p-values, when computed, are the p_boot/p_boot_approx columns of get_summary)

get_variant(::Dynema.DynemaModel)

Extract variant names provided as column names in genotypying data from a DynemaModel

get_B(::Dynema.DynemaModel)

Extract number of bootstrap iterations applied iteratively for a DynemaModel

get_bootdists(::Dynema.DynemaModel)

Extract bootstrap stat distributions for each variant

get_time(::Dynema.DynemaModel)

Extract total elapsed time in seconds

get_stattype(::Dynema.DynemaModel)

Extract statistic type: z or χ² (default), or t or F when the small-sample correction (small = true) is used

get_testtype(::Dynema.DynemaModel)

Whether a score or wald test were performed

get_boot(::Dynema.DynemaModel)

Determine whether bootstrapping was perfomed

get_pos(::Dynema.DynemaModel)

Extract genomic position for each variant

get_gene(::Dynema.DynemaModel)

Extract name for tested gene

get_chr(::Dynema.DynemaModel)

Extract name for tested gene

set_pos!(::Dynema.DynemaModel)

Sets positions for all variants for a DynemaModel

set_gene!(::Dynema.DynemaModel)

Sets gene name for a DynemaModel

set_chr!(::Dynema.DynemaModel)

Sets chromosome name for gene tested

extract_geno_dataframe(; vcf, chr=nothing, tss=nothing, bed=nothing,
                        window=500_000, field="auto", samples_file=nothing,
                        donor_col="donor_id", maf=0.0, max_missing=0.1, verbose=true)

Core VCF cis-window extraction: resolves the region around a TSS (from a single-gene bed-like bed file – columns chr, start, end, gene, strand – FastQTL-style: gene start on the + strand, gene end on the - strand; or given directly via chr/tss), queries a bgzipped, tabix-indexed VCF for that region, converts each variant's GP/DS FORMAT field to an expected allele dosage (preferring GP, per field), optionally remaps VCF sample names to donor ids via samples_file (columns vcf_id, donor_id), and drops variants that are multiallelic, lack a usable dosage field, or exceed max_missing, or fall below maf. Remaining missing dosages are mean-imputed.

Returns a NamedTuple with:

• geno::DataFrame: one donor_col column plus one column per retained variant, ready to hand to map_locus (via expand_genotypes) directly.
• positions::Vector{Int}: each retained variant's genomic position (the VCF POS field), aligned with geno's variant columns.
• chr, tss, start_pos, end_pos: the resolved cis-window.
• n_samples_vcf, n_samples_matched: sample counts before/after samples_file matching.
• n_variants_total, n_multiallelic, n_no_field, n_high_missing, n_retained: how many variants were found in the region and why any were skipped.

Set verbose=false to suppress the built-in progress printlns – e.g. for a caller (like the dynema-map CLI) that wants to render its own progress messages from the returned counts instead.

extract_gene_expression(; mtx, features, barcodes, gene, id_col="cell_id", verbose=true)

Core Matrix Market single-gene extraction: streams through mtx once, extracting the expression of gene (matched against features, one entry per matrix row) across every cell (one entry per matrix column, named from barcodes), without ever loading the whole sparse matrix into memory (see module-level comment above). Matrix entries not present in the sparse file are zero.

gene is a gene name/symbol or a gene id. The features file may be a single column of gene names (e.g. a Seurat export), in which case gene must match that column, or the two/three-column 10x layout (geneid, genename[, modality]), in which case gene is searched against genename (column 2) first and, failing that, against geneid (column 1). Ensembl id version suffixes are ignored; no match is an error.

Returns a NamedTuple with:

• expr::DataFrame: two columns (id_col, then gene as provided), shaped exactly like a plain dynema-map --expr TSV would be.
• target_row, n_genes: which row gene was found at, and how many genes/rows features has in total.
• n_cells: number of cells/columns (from barcodes).
• n_found: number of nonzero entries read for gene.

If a gene-major indexed sidecar (.dgx, built once with prepare_gene_expression or the dynema-prepare-expr CLI) exists next to mtx, the gene is loaded from it in milliseconds; otherwise the whole file is scanned once (streaming, roughly decompression-bound).

Set verbose=false to suppress the built-in progress printlns – e.g. for a caller (like the dynema-map CLI) that wants to render its own progress messages from the returned counts instead.

resolve_mtx_triplet(prefix) -> (mtx, features, barcodes)

Given a filename prefix (e.g. "expr_0.05"), resolves the three companion files a Matrix Market export usually comes as: <prefix>.mtx, <prefix>.features, <prefix>.barcodes – each tried gzipped (<prefix>.mtx.gz, ...) first, then plain. Errors up front, naming exactly which file(s) are missing and every path that was tried for each, rather than failing deep inside the (potentially minutes-long) matrix scan.

prepare_gene_expression(; mtx, out = nothing, verbose = true) -> String

One-time conversion of a (possibly gzipped) cell-major Matrix Market count matrix into a gene-major indexed binary sidecar (.dgx, written next to mtx unless out is given), after which extract_gene_expression – and therefore dynema-map --expr-prefix – loads any single gene in milliseconds instead of rescanning the whole file.

Runs two streaming passes over mtx (count entries per gene, then place them), so it needs roughly 8 bytes of RAM per nonzero entry (~8 GB per billion entries) plus the file's decompression, and takes about twice one full scan. Values are stored as Float32 (exact for counts below 2^24).

prepare_vcf(; vcf, force = false, verbose = true) -> String

One-time preparation of a VCF for Dynema's indexed region queries, using the htslib binaries bundled with the package (htslib_jll) – no system bgzip/tabix install needed. The input file is never modified: any needed compression writes a new, Dynema-namespaced file alongside it, <name>_dynema.vcf.bgz (.bgz is htslib's extension for BGZF files) – a name no user file plausibly occupies, so nothing of the user's can be overwritten. The three states a VCF arrives in:

• plain .vcf: BGZF-compressed to a new <name>_dynema.vcf.bgz;
• .vcf.gz that is plain gzip rather than BGZF (tabix rejects these): recompressed to a new <name>_dynema.vcf.bgz;
• already BGZF: used as is (no new file).

On a rerun, an existing <name>_dynema.vcf.bgz is reused when its decompressed contents still match the input, and regenerated when they don't (a stale copy of an input that has since changed). The resulting file is then tabix-indexed (skipped if an up-to-date .tbi/.csi already exists, unless force = true; an index older than the file it indexes is rebuilt). Returns the path of the ready, indexed file – which is what to pass to dynema-map --vcf (it differs from the input when a compression step ran).

The VCF must be coordinate-sorted: Dynema deliberately never modifies or reorders a VCF's records, so an unsorted file errors with instructions to sort it with external tooling (e.g. bcftools sort) first.