R/aggregate_coding_results.R
aggregate_coding_results.RdCombines the per-chromosome output of a gene-centric coding-mask scan into one
genome-wide table and computes its significance. It reads up to three result
sources a coding run produces – GLOW's own results, and (optionally)
STAARpipeline's STAAR-O results run without and with the saddlepoint
approximation (SPA) – joins them by gene and variant category into a
side-by-side method comparison, and derives per-category and genome-wide
significance from the p-values. This is the coding-mask counterpart of
aggregate_scan_results.
aggregate_coding_results(
output_dir,
chroms,
gene_num_in_array,
staar_modes = c("spa", "nospa"),
genes_info = NULL,
alpha = 0.05,
gw_alpha = 2.5e-06,
primary_test = "GLOW_Omni",
allow_partial = FALSE,
result_prefix = "glow_coding",
type_label = "coding",
family_alpha = NULL,
gene_anchor = NULL,
reference_cutoffs = NULL
)Character. The run's outputs/ directory; the function
reads results-glow/, results-staar-<mode>/ under it.
Integer vector of the chromosomes the run covered (the GLOW
completeness denominator, e.g. config$chroms).
Integer. Genes per native STAAR array task (the
native array-task mapping, e.g. config$gene_num_in_array); used only
to reconcile native completeness.
Character vector of native modes to read (default
c("spa", "nospa"); pass character(0) for a GLOW-only run).
Optional data frame whose second column gives each gene's
chromosome (the STAARpipeline::genes_info layout). Used ONLY to derive
the native array-task completeness denominator when staar_modes is
non-empty; pass STAARpipeline::genes_info (or an equivalent table) to
enable native reconciliation. Not needed for GLOW-only aggregation
(staar_modes = character(0)) – the package has no hard dependency on
STAARpipeline.
Family-wise error rate for the per-category Bonferroni cutoff
(default 0.05).
Genome-wide significance threshold (default 2.5e-6,
the STAARpipeline convention).
Column ranking the significant-cells table and driving the
within-category order (default "GLOW_Omni").
Logical (default FALSE). When FALSE, a
missing GLOW chromosome file or native array task is an error; when TRUE it is a
warning and the assembly proceeds on the present files.
Character. The stem of the per-chromosome GLOW files
under results-glow/: "glow_coding" (default) reads
glow_coding_chr<N>.csv; "glow_noncoding" reads the
gene-centric non-coding masks' tables.
Character. Names the returned tables: "coding"
(default) gives glow_coding_all, coding_merged,
coding_sig_by_category, coding_sig_by_category_cells;
"noncoding" the same with that word.
Numeric or NULL (default). When given, the across-family
cutoff is computed as family_alpha / n_all, with n_all the
tested cells of the method over all categories, in place of the fixed
gw_alpha; both denominators are reported in thresholds.
A data.frame with gene, chr and
anchor (an integer position per gene), or NULL (default). When
given, rows are ordered by (chr, anchor, gene, category) instead of
(chr, start, category), so that a gene's categories stay together
when their spans differ (the non-coding masks); a gene absent from the
table falls back to its row's start. See
noncoding_gene_anchor.
A data.frame with category, cutoff
and derivation, or NULL (default): published constants printed
beside the computed cutoffs in thresholds (the category
"all" for the across-family one). They are never applied.
A named list (the four table names below take type_label in
place of coding):
Genome-wide GLOW table (gene x category), with the
*_glowG rename and label, chr, category fronted.
Genome-wide native tables per mode (0-row schema-only data.frame if a mode found nothing or is not requested).
One row per (gene, category): the four method
columns plus each source's sub-tests and count columns, sorted gene-first.
Per (category, method): n_tested,
thresholds, and significant-cell counts.
One row per significant
(gene, category) cell, gene-first, with every method's p-value, the
count chain, a global GLOW rank, and the significance flags.
The method columns present in coding_merged.
Character vector reconciling expected vs present GLOW chromosomes and native array tasks (the completeness-report lines).
One row per (method, category) and per
(method, "all"): n_tested, the computed cutoff, how it
was computed (cutoff_rule), and the reference constant with its
derivation when reference_cutoffs gives one.
One row per requested chromosome: its
status ("present", "empty results file" when every
cell skipped or errored, "no results file" under
allow_partial), the cells (rows) and genes, the tested cells of the
primary test and the rows whose primary p-value is missing, the cells
attempted, ok, skipped (n_skip_empty, n_staar_skip_single)
and errored read from <prefix>_chr<N>_summary.csv when that
sidecar exists (the non-coding type; NA otherwise), and the error rows
read from <prefix>_chr<N>_errors.csv when that file exists (NA
otherwise).
Three sources.
GLOW – the per-chromosome glow_coding_chr*.csv
row-bound genome-wide. The two baked-in STAAR omnibus columns on GLOW's
exact post-filter genotype matrix (STAAR_O, ACAT_O) are
renamed STAAR_O_glowG / ACAT_O_glowG so the merge can tell
them from the native engine's STAAR_O_native.
native non-SPA – the results-staar-nospa/staar_coding_*.csv
per array-task tables row-bound (columns suffixed _native).
native SPA – the results-staar-spa/staar_coding_*.csv
per array-task tables row-bound (columns suffixed _spa).
Merge. GLOW is the base: it emits a cell for every \(\ge 1\)-variant gene; native only for \(\ge 2\)-variant cells, so native \(\subseteq\) GLOW. A left join from GLOW captures every real cell; native cells with no GLOW row are reported category-aware – a cell in a GLOW-tested category is WARNED (a genuine extraction mismatch), one in a standard category this run's GLOW set does not test is INFO (expected: native always emits all 7 standard masks).
Significance (computed at runtime, never hardcoded). For each
(category, method) with \(n_{tested} > 0\) finite p-values, the
per-category Bonferroni cutoff is \(\alpha / n_{tested}\) and the genome-wide
cutoff is gw_alpha; counts of cells clearing each are tallied. The
significant-cells table lists every (gene, category) cell that any
present method makes significant (its per-(category, method) Bonferroni
cutoff, looked up from the counts table, or gw_alpha).
Completeness. The number of expected vs present GLOW chromosomes and
native array tasks is reconciled; allow_partial = FALSE (default) stops on
any gap (a missing chr/task distorts the significance denominators),
TRUE downgrades it to a warning and proceeds on what is present.
Gene-first output. coding_merged and
coding_sig_by_category_cells are sorted by
(chr, gene_start, category) – one block per gene at its genomic
position with its categories underneath. glow_coding_all,
staar_coding_nospa_all and staar_coding_spa_all keep their
(chr, start[/gene], category) order.
Li, X., Li, Z., Zhou, H., et al. (2020). Dynamic incorporation of multiple in silico functional annotations empowers rare variant association analysis of large whole-genome sequencing studies at scale. Nature Genetics, 52, 969-983. doi:10.1038/s41588-020-0676-4
Li, Z., Li, X., Zhou, H., et al. (2022). A framework for detecting noncoding rare-variant associations of large-scale whole-genome sequencing studies. Nature Methods, 19, 1599-1611. doi:10.1038/s41592-022-01640-x
if (FALSE) { # \dontrun{
agg <- aggregate_coding_results(
output_dir = "runs/maf0.5_snv_7cat/outputs",
chroms = 1:22,
gene_num_in_array = 500,
staar_modes = c("spa", "nospa"))
head(agg$coding_merged)
agg$coding_sig_by_category
# The gene-centric non-coding masks: the same function pointed at the
# glow_noncoding_chr*.csv tables, the across-family cutoff computed from the
# tested cells (family_alpha / n_all) and the rows ordered by the gene anchor.
# staar_modes = character(0) reads no native STAARpipeline tables; the
# comparator columns inside the GLOW tables (STAAR_O_glowG, ACAT_O_glowG)
# are kept either way.
tbl <- do.call(rbind, lapply(1:22, function(k) utils::read.csv(
sprintf("runs/myrun/outputs/results-glow/glow_noncoding_chr%d.csv", k),
stringsAsFactors = FALSE, colClasses = c(chr = "character"))))
nc <- aggregate_coding_results(
output_dir = "runs/myrun/outputs",
chroms = 1:22,
gene_num_in_array = 500,
staar_modes = character(0),
result_prefix = "glow_noncoding",
type_label = "noncoding",
family_alpha = 0.05,
gene_anchor = noncoding_gene_anchor(tbl))
nc$thresholds
nc$completeness_table
} # }