Contents

compute$scores

Compute polygenic scores

Computes PRS values for every model in a PolyGeniusStudy by driving PLINK2 --score over the study's genotypes. Returns the score matrix and writes nothing into study; the caller assigns the result into a data$scores layer.

Usage

compute$scores(
  study,
  maf.thr = 0,
  model.filter = function(variants) rep(TRUE, nrow(variants)),
  models = NULL,
  simplify = FALSE,
  memory = NULL,
  nthreads = NULL
)

Arguments

Argument	Description
`study`	A study whose own model library is empty is fine as long as models supplies them; the call aborts only when neither does.
`maf.thr`	Numeric scalar in [0, 0.5], default 0. Variants whose in-sample frequency falls outside [threshold, 1 - threshold] are excluded and reported as maf.filtered; any value above 0 forfeits the fused scoring path.
`model.filter`	Function of one argument, default function(variants) rep(TRUE, nrow(variants)). Called once per model with that model's variant table and must return a logical vector of length nrow(variants) with no NA; anything else aborts, naming the offending model. A supplied filter also disables the file-backed regime for the call, and a non-identity filter that removes nothing across the whole library logs a warning. Check the chromosome spelling before comparing: a model's chr column may be "19" or "chr19", so v$chr == 19 can silently match nothing and leave a negated filter retaining every row.
`models`	A PGS, a PGSLibrary, a variant data.frame, a named list of variant data.frames, a list of PGS objects, or omitted (the default). Explicit models to score instead of the ones attached to study; omitting it scores every attached model. Supplying it as NULL aborts, because that is an expression that resolved to nothing rather than a request to score everything. Each variant table needs chr, position, ea, nea and beta columns, and a raw table additionally requires study's genotype to expose a genome build. Any other shape aborts. Score columns are named after the models. For a list, the list's own names win over each PGS's $name, and an entry left unnamed keeps its model's; a list of variant data.frames has no other name source, so an unnamed entry there aborts. Two models may share a name, and the returned matrix then carries that label twice. A name labels a column; position identifies it. Subset such a layer by column number, not by name -- scores[, "shared"] returns the first match twice.
`simplify`	Logical scalar, default FALSE. When TRUE and exactly one model is scored, the result is dropped to a named numeric vector instead of a one-column matrix.
`memory`	Positive numeric scalar (GiB), a named list, or NULL (default). Memory budget; NULL resolves from workspace$config$max.memory, and a bare scalar is a total GiB budget of its own. The list form overrides any subset of total (GiB), plink (MB), cells (integer cell count) and file.backed (integer cell count); any other key aborts. See Memory and batching.
`nthreads`	Positive integer scalar, or NULL (default). Passed to PLINK --threads; NULL omits the flag and leaves PLINK's own default.

Value

A numeric matrix of PRS values (samples x models), columns in the requested model order, or — with simplify = TRUE and exactly one model — a named numeric vector. Either carries a provenance record (provenance()) whose misc holds strategy ("fused" or "via.extract"), n.batches, n.models, n.variants.total and maf.thr, plus two diagnostics entries, coverage (diagnostics(x, "coverage")) and the variant fate record (diagnostics(x, "fate")), each described in its section above.

Details

The variants of every scored model are collated into one union and matched to the genotype records (PLINK .pvar/.bim) by chromosome and position, accepting a record when the alleles agree in either orientation (nea == ref & ea == alt, or nea == alt & ea == ref). Each matched record takes the canonical id chr:pos:ref:alt, the same id PLINK assigns via --set-all-var-ids, so the score file needs no rename step. Orientation is resolved by PLINK and never in R: the effect allele goes into the score file's EA column and --score counts that allele's dosage directly, so no dosage is ever flipped.

Genotypes split across several files are matched and scored per file and the per-file PRS values summed, which is exact because a PRS is additive. Files are scored sequentially; nthreads parallelises only within a single PLINK call.

A model variant matching more than one genotype record is excluded and reported as collision. Multi-allelic records (--max-alleles 2) and duplicated chr:pos:ref:alt records (--rm-dup exclude-all) are dropped inside PLINK and reported as missing.

A model that lists a variant more than once -- one loaded from an older library or .pgd, or built by a route that bypasses PGS() -- is scored under the rule PGS() applies: exact copies count once, and a variant listed with conflicting betas is left out. Either case warns, naming the models.

Memory and batching

memory resolves into three knobs at three layers: cells, this call's per-batch cap in cells (V_batch x M_batch), not bytes; plink, PLINK's own --memory flag in MB, a genuine concurrent reservation; and file.backed, the dense backbone$n.variants() x n.models cell count above which the scored union is read out-of-core from the model-weight store (backbone-backed input only). memory = NULL falls back to workspace$config$max.memory, whose own default Inf means no budget at all: cells = 200e6, PLINK uncapped, file.backed = 2.5e8. Under a total budget, plink is subtracted from total before cells is derived from the remainder, file.backed is derived from the full un-subtracted total, and any knob supplied explicitly wins verbatim over its derivation. cells is the master knob for peak memory: all models score in one fused genotype pass while the scored union times the model count fits it, and switch to a batched via-extract path once it does not. maf.thr > 0 also forfeits the fused path but adds no extra genotype read, since the frequency computation is folded into the same canonical-extraction pass.

Variant fate record

diagnostics(x, "fate") is a PolyGeniusVariantFate, holding what became of every requested variant: counts (named integers over scored, missing, maf.filtered, no.freq, collision, counted over unique variants), n.requested, n.flipped, n.strand.ambiguous, and counts.by.genotype when more than one genotype contributed.

Printing it shows those counts and a short head; as.data.table() and as.data.frame() expand it into one row per requested variant per genotype — chr, position, nea, ea, status, matched, flipped, strand.ambiguous, freq, n.records. It stores a byte per variant per genotype rather than that table, so the detail is always kept and there is no retention mode to choose. It sits beside coverage in diagnostics() rather than in provenance, whose misc must stay bounded.

provenance(x)$params records every supplied argument as the expression the caller wrote, so a model library or a filter closure is named rather than retained inside the study's serialized form. What the request resolved to is in misc: n.models and n.variants.total, both counted after model.filter, and the applied maf.thr. n.variants.total is the sum of per-model variant counts, so it exceeds the number of distinct variants whenever two models share one.

Coverage diagnostics

diagnostics(x, "coverage") is a data.table, one row per model in column order, with columns model, n.weights, n.scored, weight.total, weight.scored, coverage.var, coverage.count and coverage.palindromic. coverage.var is sum(beta^2 * 2p(1-p)) over scored variants divided by the same sum over every requested variant, with p the model's own source-plane eaf; it is NA for a model whose variant table carries no usable eaf, and coverage.count (n.scored / n.weights) is the always-computable naive fallback. A variant a model lists more than once counts once in n.weights; one listed with conflicting betas is never scored, so it counts in n.weights but not in n.scored, and adds nothing to the weight sums. coverage.palindromic is the strand-ambiguous share of sum(beta^2) over every requested variant, which coverage.var cannot supply: a model scored perfectly backwards still reads coverage.var = 1.0. The entry stays on this layer only — assigning the layer into data$scores copies it nowhere else; see diagnostics().

Examples

data$scores$X <- compute$scores(data)

# Cap peak memory at 8 GiB.
data$scores$X <- compute$scores(data, memory = 8)

# Score an ad hoc model library, dropping near-zero weights.
scores <- compute$scores(data, models = my.models,
                         model.filter = function(v) abs(v$beta) > 1e-4)
coverage <- diagnostics(scores, "coverage")