GLOWpipeline: Orchestrating a Genome-Wide GLOW Scan

Source: vignettes/glowpipeline-overview.Rmd

Introduction

GLOWpipeline is the orchestration-and-summary layer of the GLOW software family. The methods layer, GLOWr, provides the single-region variant-set test (GLOWr::glow_region_test()); GLOWpipeline scales that single-region test up to a genome-wide scan and turns the per-region output into genome-wide tables and diagnostics. It is the STAARpipeline analog in the GLOWpipeline -> GLOWr -> GFisher stack.

A genome-wide GLOW scan has four conceptual phases:

  1. Set up the shared inputs once (prepare_scan_run()).
  2. Chunk the genome into array tasks (define_genome_chunks()).
  3. Scan each chunk’s windows with GLOW (glow_scan_chunk() / run_scan_unit()), typically one compute-cluster array task per chunk.
  4. Aggregate and summarize the per-chunk outputs into a genome-wide table (aggregate_scan_results()) with diagnostics (summarize_qq_and_lambdas()).

GLOWpipeline is deliberately cohort- and cluster-agnostic: resolved objects go in, tables come out. A workflow script binds it to a dataset and a scheduler. This vignette shows the call shapes; the data-dependent steps are not evaluated here (they need a GDS file and a fitted null model), while the two self-contained steps run for real.

1. Set up the run

prepare_scan_run() resolves the parts of a run that every chunk task shares – the fitted null model, the optimal-weight inputs, and the output directories – so the work is done once rather than per task. It takes a resolved config list and an output_dir.

spine <- prepare_scan_run(
  config     = config,        # resolved analysis config (cohort + method knobs)
  output_dir = "scan_out",
  chroms     = 1:22
)

The returned setup carries the null model and weight inputs that the scan consumes in phase 3.

2. Chunk the genome

define_genome_chunks() splits a set of chromosome spans into contiguous sliding-window chunks. Each row is one chunk – one compute-cluster array task – mapping an array index to a chromosome and a contiguous run of window start positions. The boundary rule guarantees the union of all chunks’ windows equals one full-span slide with no duplicates and no gaps.

This step is self-contained, so we run it for real on a small two-chromosome example:

chrom_spans <- data.frame(
  chr   = c("21", "22"),
  start = c(5010000L, 10510000L),
  end   = c(48120000L, 50810000L)
)

chunks <- define_genome_chunks(
  chrom_spans,
  window_size    = 2000,
  step_size      = 1000,
  chunk_strategy = "fixed_mb",
  chunk_size_mb  = 5
)

nrow(chunks)
#> [1] 18
head(chunks)
#>   chunk_id chr chunk_start chunk_end span_end est_n_windows
#> 1        1  21     5010000  10009999 48120000          5000
#> 2        2  21    10010000  15009999 48120000          5000
#> 3        3  21    15010000  20009999 48120000          5000
#> 4        4  21    20010000  25009999 48120000          5000
#> 5        5  21    25010000  30009999 48120000          5000
#> 6        6  21    30010000  35009999 48120000          5000

Each chunk’s chunk_id is the value an array task reads (e.g. from SLURM_ARRAY_TASK_ID) to learn which chromosome and window range it owns.

3. Scan one chunk

An output policy declares which artifacts the scan writes. The flat per-window results table is always on; the heavier per-window evidence, STAAR detail, and result objects default off. This constructor is self-contained:

policy <- glow_output_policy()
policy
#> <glow_output_policy>
#>   flat_table     : on (always; the load-bearing output)
#>   evidence       : off 
#>   staar_detail   : off (needs a staar_context at the engine)
#>   result_objects : off

With the spine, a chunk’s window range, and the policy in hand, run_scan_unit() builds the chunk’s regions and filters and runs glow_scan_chunk() over them, persisting the results. In a cluster job this is the per-array-task call:

# `gds_path` and the null-model / weight inputs come from the shared setup;
# one of spine$b_func and spine$b_model is NULL, the other the B source.
k <- 1L                                   # this task's chunk (SLURM_ARRAY_TASK_ID)
run_scan_unit(
  gds_path          = gds_path,
  null_model        = spine$null_model,
  pi_models         = spine$pi_models,
  pi_features       = spine$pi_features,
  reference_medians = spine$reference_medians,
  b_func            = spine$b_func,
  b_model           = spine$b_model,
  region_type       = "window",
  config            = config,
  chunk             = chunks[k, ],        # one row of define_genome_chunks()
  results_path      = file.path("scan_out", "results", sprintf("scan_chunk_%d.csv", k)),
  output_dir        = "scan_out",
  output_policy     = policy
)

The gene-centric non-coding masks are the same call with region_type = "noncoding" and a chromosome in place of a chunk; config$noncoding_categories selects the masks (default all eight) and n_cores the forked workers:

run_scan_unit(
  gds_path          = gds_path,
  null_model        = spine$null_model,
  pi_models         = spine$pi_models,
  pi_features       = spine$pi_features,
  reference_medians = spine$reference_medians,
  b_func            = spine$b_func,
  b_model           = spine$b_model,
  region_type       = "noncoding",
  config            = config,
  chr               = 22,
  results_path      = file.path("scan_out", "results-glow", "glow_noncoding_chr22.csv"),
  output_dir        = "scan_out",
  output_policy     = policy,
  n_cores           = 4
)

glow_scan_chunk() returns a list whose results element is the flat per-window table (one row per non-skipped window, with the window coordinates, variant counts, and the GLOW test p-values).

region_type selects how the unit’s SNV sets are formed, and the same call serves all four types. A unit holds many regions, and GLOW tests one set per region: the region says which variants are candidates, and the filter (config$filter_spec: QC label, variant type, MAF and MAC cutoffs, minimum set size) says which candidates stay.

  • "gene": one set per gene of the bundled gene table (or of config$resolved_gene_list), every record in the gene’s span. A gene with more records than config$gene_max_records is tiled into segments.
  • "window": one set per sliding window of a genome chunk (chunk), every record in the window.
  • "coding": one set per gene and coding mask (config$categories), the records in the gene’s span whose coding annotation matches the mask.
  • "noncoding": one set per gene and gene-centric non-coding mask (config$noncoding_categories, the eight masks of GLOWr::noncoding_masks()), the records anywhere on the chromosome whose annotation links them to the gene. build_scan_regions() groups the chromosome’s records once per mask with GLOWr::define_regions_noncoding(), which needs the aGDS path. The rows of the per-unit table carry the mask as category.

With n_cores above 1 the unit’s regions run on forked workers, scheduled by their record counts. A region’s cost is its squared record count. A region whose cost reaches the unit’s total cost divided by 8 times n_cores runs alone, the rest are packed into blocks of at most 25 regions whose summed cost stays within that threshold, and the blocks run heaviest first. The written table is the single-core one.

4. Aggregate and summarize

After all chunk tasks finish, aggregate_scan_results() concatenates the per-chunk result files into a single genome-wide table (ordered by chromosome then position), flags significant windows against a Bonferroni-style threshold, and merges overlapping significant windows into loci.

agg <- aggregate_scan_results(
  chunk_files  = chunk_files,    # the per-chunk result files written in phase 3
  primary_test = "GLOW_Omni",
  alpha        = 0.05
)

agg$scan_results_all   # genome-wide flat table (all tested windows)
agg$scan_sig_windows   # the significant windows

For a gene-centric coding scan, aggregate_coding_results() is the analog. It reads results-glow/glow_coding_chr<N>.csv under output_dir. The gene-centric non-coding masks use it too: result_prefix = "glow_noncoding" reads their per-chromosome tables, type_label = "noncoding" names the returned tables, gene_anchor (from noncoding_gene_anchor()) orders a gene’s masks together although their spans differ, and family_alpha sets the across-mask cutoff, which is family_alpha divided by the number of tested (gene, mask) cells. Its thresholds table lists every cutoff with its denominator, and its completeness_table one row per requested chromosome. staar_modes = character(0) reads no native STAARpipeline tables; the comparator columns inside the GLOW tables are kept either way.

tbl <- do.call(rbind, lapply(1:22, function(k) read.csv(
  sprintf("scan_out/results-glow/glow_noncoding_chr%d.csv", k),
  stringsAsFactors = FALSE, colClasses = c(chr = "character"))))
nc <- aggregate_coding_results(
  output_dir        = "scan_out",
  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$noncoding_merged
nc$thresholds
nc$completeness_table

Genomic-control diagnostics come from summarize_qq_and_lambdas(), which computes per-test lambda values and writes QQ plots:

summarize_qq_and_lambdas(
  df           = agg$scan_results_all,
  label        = "genome",
  output_dir   = "scan_out",
  test_columns = c("GLOW_Burden", "GLOW_SKAT", "GLOW_Omni"),
  omni_test    = "GLOW_Omni"
)

For the gene and window routes, run completeness, run-log parsing and per-error variant drill-downs are available via check_scan_completeness() (chunks), parse_chr_log() / parse_chunk_log() and write_error_drilldowns(). The non-coding route writes its errors and its per-mask tallies beside each per-chromosome table instead (<stem>_errors.csv, <stem>_summary.csv), which aggregate_coding_results() reads into its completeness_table; a candidate list has no interval for a drill-down.

Where to go next