Run ancestry-stratified SMMAT + Cauchy combination (Step 3)

Source: R/ancestry_smmat.R

Wraps GMMAT::SMMAT() across all ancestry-split GDS files (as produced by write_ancestry_gds), computes per-gene ancestry weights from Step 1's ancestry_counts (this replaces the retired src/step3_weight_finding.R script), and combines the per-gene p-values across populations with cauchy_combine. This is the package's Step 3 entry point; it does not run Step 1/2 for you (SMMAT runs can take hours, so this is always an explicit, separate call).

ancestry_smmat(
  gds_paths,
  pheno,
  formula,
  kinship,
  gene_group_file,
  ancestry_counts = NULL,
  variant_info = NULL,
  id_col = "id",
  family = stats::gaussian(),
  gene_col = "group",
  p_col = "E.pval",
  MAF.range = c(0, 0.5),
  miss.cutoff = 1,
  method = "davies",
  ncores = 1,
  verbose = TRUE,
  ...
)

Arguments

gds_paths

Named list/vector of GDS file paths, one per ancestry population, e.g. c(AFR = "afr.gds", EUR = "eur.gds") – as returned by write_ancestry_gds. Any K >= 1 populations are supported; add an extra named entry (e.g. OBSERVED = ...) to include an additional (e.g. non-ancestry-split) GDS in the same SMMAT + Cauchy combination run.

pheno

data.frame with an id_col column plus the response and covariates referenced in formula.

formula

Model formula for the null model, e.g. y ~ age + sex.

kinship

Numeric matrix (or object coercible via as.matrix) with row/column names equal to sample IDs.

gene_group_file

Either a path to the gene_group file passed to GMMAT::SMMAT() (no header; columns gene, chr, pos, ref, alt, weight), or a data.frame/matrix with those same six columns. Either way, it is read, has any leading "chr" stripped from the chr column (matching SeqArray::seqVCF2GDS()'s convention – see the "chr column" note below), and written to a fresh temporary file that is deleted with unlink() when ancestry_smmat() returns; the caller's own path/data.frame is never read-modified or deleted, and can safely be reused across repeated calls (e.g. one call per phenotype). chr column: pass either "19" or "chr19" – normalized automatically. Without this, a mismatch between the group file's chr convention and the GDS's (always bare, no "chr" prefix) would silently drop every affected variant, and the whole gene if all its variants are affected, since GMMAT::SMMAT() matches variants by exact chr:pos:ref:alt string equality with no fuzzy/positional alignment.

ancestry_counts, variant_info

Optional: pass ancestry_split()$ancestry_counts / $variant_info straight through to compute per-gene ancestry weights for the Cauchy combination. If omitted, equal weights are used.

id_col

Name of the sample ID column in pheno and formula's id argument to glmmkin (default "id").

family

A family object for GMMAT::glmmkin() (default gaussian()).

gene_col

Name of the gene/group ID column in SMMAT()'s output (default "group", matching GMMAT::SMMAT()).

p_col

Name of the p-value column in SMMAT()'s output to combine (default "E.pval", matching GMMAT::SMMAT()'s default tests = "E"; use "O.pval"/"S.pval"/ "B.pval" if you pass a different tests via ...).

MAF.range, miss.cutoff, method, ncores

Passed straight through to GMMAT::SMMAT() – same names, same meaning as in a direct SMMAT() call.

verbose

Print this function's own step-by-step progress messages (Steps 1-6 below). Distinct from GMMAT::SMMAT()'s own internal verbose (its per-call progress bar), which is not exposed separately and stays at SMMAT's default (FALSE); every other SMMAT() argument is available via ....

...

Any other GMMAT::SMMAT() argument, forwarded as-is – for example tests ("B", "S", "O", or the default "E"), rho, MAF.weights.beta, missing.method, use.minor.allele, auto.flip, Garbage.Collection, group.file.sep, or meta.file.prefix. null.obj, geno.file, group.file, and is.dosage are always supplied internally by ancestry_smmat() and cannot be overridden this way. See ?GMMAT::SMMAT for the full list and defaults.

Value

A list with two elements:

results

data.frame with one row per gene: gene, one p_<POP> column per population, one w_<POP> weight column per population (only if ancestry_counts/ variant_info were supplied), and p_cauchy.

smmat_results

Named list (one entry per population, same names as gds_paths) of the raw GMMAT::SMMAT() result data.frames.

Examples

if (requireNamespace("GMMAT", quietly = TRUE) &&
    requireNamespace("SeqVarTools", quietly = TRUE)) {
  ext <- system.file("extdata", package = "lantern")
  split <- ancestry_split(file.path(ext, "toy.vcf.gz"),
                          file.path(ext, "toy.msp.tsv"),
                          mode = "dosage", chrom = "chr19", verbose = FALSE)
  out_dir <- tempfile()
  gds <- write_ancestry_gds(split, out_dir, verbose = FALSE)

  pheno <- read.delim(file.path(ext, "toy_pheno.tsv"))
  kinship <- diag(nrow(pheno))   # unrelated samples
  dimnames(kinship) <- list(pheno$id, pheno$id)

  # gene_group_file: a path, or a data.frame with the same six columns
  out <- ancestry_smmat(gds, pheno, y ~ age + sex, kinship,
                        file.path(ext, "toy_genes.tsv"),
                        ancestry_counts = split$ancestry_counts,
                        variant_info    = split$variant_info,
                        verbose = FALSE)
  out$results              # per-gene p-values + Cauchy combination
  head(out$smmat_results$AFR)   # raw GMMAT::SMMAT() output for AFR
  unlink(out_dir, recursive = TRUE)
}