R/build_scan_regions.R
build_scan_regions.RdRegion dispatch for the unified GLOW scan. Returns the
per-chromosome regions table plus the variant filter(s) for a single
scan unit (a chromosome for gene/coding, a genome chunk for window). This is
the one structural fork of the scan: gene/window run a single
filter over many regions, while coding runs many filters (one per
category) over the same regions, so the return carries either one filter
object or a named category -> filter map.
build_scan_regions(
region_type,
chr,
config,
chunk = NULL,
variant_index = NULL,
gds_path = NULL
)Character: one of "gene", "window",
"coding", "noncoding".
Integer or character chromosome.
An environment or named list of resolved config values. Fields
used per type: filter_spec (all); resolved_gene_list
(optional, gene); window_size, step_size (window);
categories (coding).
A single-row data.frame (one row of the genome-chunk table) for
region_type = "window"; ignored otherwise. Must carry
chunk_start, chunk_end, span_end. Default NULL.
A glow_variant_index of the unit's GDS from
GLOWr::build_variant_index(), or NULL (default). Used by
region_type = "gene", for the large-gene segmentation below, and by
"noncoding", whose candidate lists are row numbers in it (built
from gds_path when NULL).
Character path of the unit's GDS/aGDS, or NULL (default).
Required for region_type = "noncoding"; ignored by the others.
A named list:
A data.frame with at least chr, start,
end, label (the define_regions_* shape); for
noncoding, a NAMED list of such tables, one per mask, each with
the n_records, category, gene_rank and
variant_rows columns of GLOWr::define_regions_noncoding.
A single glow_variant_filter (gene/window) OR a
named list mapping each category label to a glow_variant_filter
(coding, noncoding).
region_type = "gene". regions is
GLOWr::define_regions_gene(chr = chr) — the bundled STAARpipeline gene
table for the chromosome — unless config$resolved_gene_list is present
(a user-supplied custom gene list), in which case that list is subset to
chr. filters is the single filter
do.call(GLOWr::variant_filter, config$filter_spec).
region_type = "window". Requires a chunk row (one row
of the genome-chunk table from define_genome_chunks, with
chunk_start, chunk_end, span_end). regions is
GLOWr::define_regions_window(chr, chunk$chunk_start, chunk$span_end,
config$window_size, config$step_size) trimmed by the sliding boundary rule
(keep start <= chunk$chunk_end): windows are generated to the
chromosome span_end so boundary-spanning windows are kept, then only
those owned by this chunk (start within [chunk_start, chunk_end]) are
retained, so the union of chunks equals one full-span slide. filters is
the single variant_filter from config$filter_spec.
region_type = "coding". regions is the chromosome gene
table (as for gene). filters is a NAMED list mapping each
category label to its resolved glow_variant_filter:
config$categories is either a character vector of builtin coding
category names (each resolved via GLOWr::coding_filter) or a named list
mapping labels to a builtin name OR a list of custom DNF annotation clauses
(resolved via GLOWr::variant_filter with the run's
filter_spec). The MAF/MAC/min-variant axes come from
config$filter_spec so a category never hardcodes them.
region_type = "noncoding". The gene-centric non-coding masks
(see GLOWr::define_regions_noncoding). Needs gds_path, because
the regions come from the chromosome's annotations, not from gene
coordinates: the chromosome's QC, type, category and presence vectors are
read once, and for every mask in config$noncoding_categories (default
GLOWr::noncoding_masks(), the eight masks) one regions table is built
whose rows are the genes with linked records and whose variant_rows
column lists each gene's candidate records. regions is then a NAMED
list of those tables, one per mask, and filters the matching map of
GLOWr::noncoding_filter(mask, run filter) (the run's filter with no
annotation clause). The GDS handle is closed on return.
At whole-genome-sequencing density a gene body can hold tens of thousands
of records, and the per-gene cost and memory grow with the square of the
variant count. When variant_index is supplied and
config$gene_max_records is a finite number, every gene's GDS records
inside its span are counted from the index (before any filter, so the count
depends only on the file), and a gene with more records than that is
replaced by non-overlapping segments of config$gene_segment_bp base
pairs (default 10,000) that start at the gene start; the last segment ends
at the gene end. Each segment is its own region row with the gene's
region_id, its own start and end, and
label = <gene label>_seg<k>, and it is tested and reported like a
gene; there is no gene-level combination of the segments. A gene whose span
fits one segment is left as it is. The returned table then also carries
n_records (the index count of each row) and the attributes
n_genes_segmented and n_segments. With variant_index
NULL, or gene_max_records NULL or infinite, nothing is segmented.
if (FALSE) { # \dontrun{
env <- new.env(); source("runs/myrun/config.R", local = env)
# Gene: one filter, the chr gene table.
g <- build_scan_regions("gene", chr = 22, config = env)
# Window: trimmed to the chunk boundary.
chunk <- define_genome_chunks(22, window_size = env$window_size,
step_size = env$step_size)[1, ]
w <- build_scan_regions("window", chr = 22, config = env, chunk = chunk)
# Coding: a category -> filter map.
co <- build_scan_regions("coding", chr = 22, config = env)
names(co$filters)
# Non-coding: one regions table per mask, from the chromosome's annotations.
nc <- build_scan_regions("noncoding", chr = 22, config = env,
gds_path = "data/chr22.gds")
sapply(nc$regions, nrow)
} # }