Reads genotype data from a GDS/aGDS file, applies QC, variant type, positional, and annotation-based filters, and returns a prepared variant set for downstream processing and testing.
extract_variant_set(
gds,
region,
filter_spec,
sample_id = NULL,
annotation_names = NULL,
Annotation_dir = "annotation/info/FunctionalAnnotation",
Annotation_name_catalog = NULL,
impute_method = "mean",
variant_index = NULL,
verbose = 1
)An open SeqVarGDSClass object or character file path.
If a path, the file is opened and closed internally.
A list or single-row data.frame with elements/columns:
chr, start, end. Optional: region_id, label, and variant_rows, an
integer vector of row numbers in the chromosome's variant index (the
candidate list of define_regions_noncoding): unique,
strictly increasing, whole numbers within the index, also accepted as
the one-row list column of a regions table. When present it replaces the
interval lookup by start and end; integer(0) is an
empty set and returns NULL, and an invalid list stops with the
region's label. Everything after the lookup (QC, variant type, MAF, MAC,
the annotations) applies as for an interval.
A glow_variant_filter from variant_filter()
or coding_filter().
Character vector of sample IDs for extraction and ordering. Must match a subset of GDS sample.id. If NULL, all samples are used in GDS order.
Character vector of annotation names to extract for PI computation (e.g., c("CADD", "LINSIGHT")). Uses the annotation name catalog to resolve GDS paths. NULL means no annotations extracted.
Character path prefix for annotations in GDS. Default "annotation/info/FunctionalAnnotation".
Data frame mapping annotation names to GDS sub-paths. NULL uses the bundled catalog.
Character: "mean" or "minor". Default "mean".
"mean" replaces NA with the column mean (2 * AF), the expected
dosage under Hardy-Weinberg. "minor" replaces NA with 0 when
AF <= 0.5, or 2 when AF > 0.5; in both cases the imputed value assumes
the individual carries zero copies of the minor allele. This matches
the legacy STAARpipeline matrix_flip_minor() behavior.
A glow_variant_index from
build_variant_index for this GDS, or NULL (default). With an
index the region's records are found by binary search on the index and
the GDS filter is set by file position, so a scan over many regions does
not re-read the chromosome's variant table for each one. Without it the
function builds an index for this call, which costs one read of the
chromosome, as before. The selection and the returned object are the
same either way.
Integer verbosity level.
A glow_variant_set S3 object, or NULL if the region
has fewer than filter_spec$min_variants qualifying variants.
Cost: with a prebuilt index, locating a region costs two binary searches on
the chromosome's positions, so the per-region cost is set by the filtered
records of the region alone. Without an index, each call reads the
chromosome's variant table once (linear in the chromosome's record count).
The orchestration layer (GLOWpipeline::glow_scan_chunk) builds the
index once per chunk and passes it to every region.
if (FALSE) { # \dontrun{
library(SeqArray)
gds <- seqOpen("chr22_essentialdb.gds")
# Extract plof variants for a gene
region <- list(chr = "22", start = 17565844L, end = 17596578L,
label = "CECR1")
spec <- coding_filter("plof")
vset <- extract_variant_set(gds, region, spec,
sample_id = sample_ids,
annotation_names = c("CADD", "LINSIGHT"))
dim(vset$G) # n_samples x n_variants
vset$variant_info # variant metadata
vset$annotations # annotation score matrix
seqClose(gds)
} # }