A PolyGenius Study

Where samples come from

The sample list is the genotype files' sample list. sample.names is built from them, in file order, and nothing else defines it.

study$samples$<name> <- value places each item of value on a sample. Items that carry sample IDs (a vector's names, a table's row names, or an ID column named through align()) are matched to sample.names. A sample the value doesn't mention gets NA, and an item for an unknown sample is dropped. Items without IDs are taken in order, and there must be exactly one per sample.

pheno <- read.csv("cohort-phenotypes.csv")          # own row order, an IID column, some samples missing
study$samples$phenotypes <- align(pheno, by = IID)  # matched by IID; the missing samples get NA

The case to watch is a table with no IDs at all. A read.csv() result with the right number of rows but a different order than the genotype files is accepted as it is, and every sample silently gets someone else's values. Give a table IDs, through align() or row names, whenever its order is not the genotype files' order.

In a study built from several genotypes, the same ID can appear in more than one of them. A value matched by ID aborts on such an ID, because it could mean either sample. Give one value per genotype instead, as a list named by genotype: list(gsa = pheno.gsa, omni = pheno.omni). The Assigning sample-level data section of the PolyGeniusStudy() reference lists every accepted form.

To add several values at once, assign to study$samples itself, or pass the same value to PolyGeniusStudy(samples = ). A table's columns become columns of the study's sample table. A named list assigns each element under its own name. A list named by genotype takes a table, or a named list of values, for each genotype:

study$samples <- align(pheno, by = IID)                   # every column of pheno joins the table
study$samples <- list(gsa = pheno.gsa, omni = pheno.omni)  # a column one table lacks is NA for that genotype

A name that already exists aborts, and nothing is written. study$samples$<name> <- NULL removes a value.

One study, one genotype set, or many

GenotypeSource names one study's files. That can be many files — a genome split by chromosome is still one study.

GenotypeSourceSet is different: it concatenates different sample sets along the sample axis. Registering chromosomes as a set rather than as one genotype multiplies your sample count instead of combining variants.

Sample IDs are unique within one genotype but may collide across a set, which is why character-based subsetting can select more rows than you intended, and why relatedness is only ever evaluated within one fileset, never across two.

Constructing the object

PolyGeniusStudy(name, library, genotypes, samples = NULL, sample.pairs = NULL, uns = NULL)

library is a PGS or PGSLibrary; genotypes is a GenotypeSource or GenotypeSourceSet. A bare PGS/GenotypeSource is coerced automatically. If you don't have either yet: a GenotypeSource points at PLINK files you already have on disk (see Your First Analysis for the minimal call), and a PGS/PGSLibrary comes from importing a published score or building one — see Model I/O and the PGS Catalog and Construction Algorithms.

It validates one thing you cannot recover from later: the models' build must equal the genotypes' build. That check happens here and nowhere else.

Some things are filled in for you. Each sample's study$samples$genotype column records which genotype (fileset) it came from. Model variant tables and their weights are not copied onto the study one model at a time — the study narrows the library you passed in down to what its own models reference and keeps that as study$library, so a model's variants and weight column are always read through it, never duplicated per model. GWAS and algorithm provenance work the same way, one level further down: they read through to study$library$backbone$models$annotations.

That narrowing is a copy, and a deliberate one: the study owns its library from then on. Absorbing more models into the PGSLibrary you still hold a binding to, renaming something in it, or subsetting it does not reach into a study already built from it — and nothing the study does reaches back out.

What it does not check: whether your phenotype table overlaps your genotypes. That happens when you assign it, as above. Nor does it check sample identity or order against the model library — there is no such axis to check, since a model library has no samples of its own, but it is worth naming as one more thing construction leaves entirely to you, in the same spirit as the gotcha above.

Built for scale

A PGS library is rarely one model. It is often a grid — the same algorithm swept across many parameter settings, or hundreds of published scores pulled from the PGS Catalog for the same trait — and those models overlap heavily in which variants they use. PolyGenius is built so that a library like this costs roughly what its distinct variants cost, not what models-times-variants would cost, because every model in the library shares one variant identity list instead of each carrying its own copy. Ten thousand C+T models swept across a parameter grid mostly reuse the same handful of underlying GWAS variants at different thresholds; PolyGenius stores that variant identity once and lets every model point at the rows it needs.

This is a deliberate design choice, not an incidental optimization, and it is worth knowing about even if you never look at the internals: growing your library doesn't quietly turn a fast job into one that runs out of memory or takes all night. A ten-model pilot and a ten-thousand-model production sweep behave the same way, qualitatively, because the cost tracks the shared variant set, not the model count.

For the actual representation behind this — the shared dictionary, how a model's weights are stored, and how it all persists to disk — see Data Architecture and Persistence and the Backbone.

The backbone and the variant dictionary, in plain terms

Two ideas explain most of what follows on this page.

The variant dictionary is the master list of every distinct variant that appears anywhere in the model library — one shared list, not one per model. A model does not carry its own private copy of "which variants, at which positions, with which alleles"; it points into this one shared list and adds its own weight for each variant it uses. Think of the dictionary as the filing cabinet the models' variant sheets are pulled from: the cabinet holds one master set of folders, and each model just says which folders are its.

The backbone is the object holding that dictionary, plus the roster of which models are registered against it. The word is not a metaphor here — "backbone" is literally what the object is called elsewhere in PolyGenius, including later in this page. Each study owns one. The sharing is within a study's library, across the models in it — not across studies: build two cohorts' studies from the same PGSLibrary and each ends up with its own backbone, carrying the same models over the same variant identities. That is what makes their results comparable, and it is a claim about content rather than about two studies holding the same object in memory.

Again, see Data Architecture and Persistence and the Backbone for what the dictionary and backbone actually look like underneath.

The slots

There is a third axis — variants — but it is not a slot on the study. It lives one level down, on the model library: study$library$backbone owns the variant dictionary every model's rows resolve against, and study$library$backbone$variants$annotations is where attached per-variant metadata ends up. Nothing in the table above is aligned to variants, so there is no study$variants to assign into.

The samples plane: one noun, three kinds of content

study$samples is not a plain data frame. study$samples$<name> <- looks at what you assign and routes it:

• a plain vector becomes a column of the primary table
• a table or matrix with one row per sample becomes a table member (a data frame is never silently squeezed into a single list-column — a rectangular value always becomes a member, even a one-column one)
• a matrix wrapped in reduction() becomes an embedding

Reading is uniform regardless of which of the three it is: study$samples$PCA just works. This is why compute$populationStructure() in Your First Analysis can write straight into data$samples$PCA — the function already returns a reduction()-tagged matrix, and the plane recognises it as an embedding rather than a column. reduction() itself is exported mainly so a producer can tag its own output; you rarely call it directly.

Models and the study's own library

study$library is the read-only handle to the library the study owns. study$model.names and study$n.models read straight through it, so they always report what that library currently holds rather than a count captured at construction. What they do not track is the PGSLibrary you passed in: that object and the study's library parted ways the moment the study was built, so absorbing more models into the one you kept a binding to adds nothing to the study. study$models is a third, separate thing again: a per-study overlay you populate yourself, not a copy of the registry.

This matters when you subset. A model-side subset returns a new study whose backbone is narrowed to exactly what the retained models reference, so the dictionary shrinks with the model list instead of dragging the parent's full variant set along. The parent keeps its own backbone, whole and unchanged.

Reading things out

$fetch() resolves unquoted names and $-paths across samples, sample.pairs, models, model.pairs and scores, so you rarely need to know which slot something lives in. It does not search uns or the three result slots (associations/evaluations/signals) — fetch those by name directly.

It reports rather than raises: an unresolved or ambiguous name comes back in an errors attribute, and the return type changes shape when that happens. Print a fetch result before building a model on it.

A bare PRS model name resolves to a column of the current scores layer (scores.layer, defaulting to X); a bare sample name resolves to a row of it. Model names may repeat, in which case fetching by bare name is refused as ambiguous.

Slot accessors

Most results carry more than their numbers. provenance() gives the record of what produced them — the call, the arguments you supplied, and what the operation resolved. artifacts() gives the derived tables that go with a result, and diagnostics() the counts and warnings. The same three verbs work whatever you hand them: a score matrix, a kinship matrix, an embedding, a model library, a study, an association or an evaluation.

print() on the study object marks each slot with what it carries — ¹ a provenance record, ² diagnostics, ³ artifacts — and a legend under it names the verb for each, so you can see what is inspectable.

Subsetting

study$subsetSamples(expr), study$subsetModels(expr), study[i, j] and subset(study, samples = ..., models = ...) all return a new object with a narrower working set. None of them mutate in place, and none of them rewrite your genotype files.

The index vocabulary is shared everywhere: names, integer positions, a logical vector, tidyselect helpers, or a predicate resolved via $fetch().

Two costs worth knowing. Subsetting rebuilds each GenotypeSource, which means one PLINK call per fileset — cheap once, not something to put in a loop. And associations, evaluations and signals are dropped whenever the subset actually removes a sample or a model — a fit of models on samples is invalidated by touching either axis, and neither can be repaired by filtering the result, so the child gets NULL instead of a silently stale number. A subset that removes nothing (subset(study), or a model subset naming every model) is a genuine no-op and keeps everything. uns keys are kept except reserved, dot-prefixed ones, which are dropped whenever the subset actually changes something. When something is dropped, the names are logged at info level.

Several cohorts

A PolyGeniusStudy holds one cohort. When you have several — a discovery and a replication set, or every cohort you plan to meta-analyse — you keep them as separate studies and call each step on each one. There is no container that groups them, and nothing in the package needs one: every verb takes a study as its first argument and writes back onto that study.

Build each one from the same model library, so every cohort is scored on the same models over the same variants:

models      <- generate$from.pgs.catalog("PGS002280")
 
discovery   <- PolyGeniusStudy(name = "discovery",   library = models, genotypes = gs.discovery)
replication <- PolyGeniusStudy(name = "replication", library = models, genotypes = gs.replication)

Passing the same library to both is fine — each study narrows its own library out of it, so neither can disturb the other's afterwards. What you should avoid is calling generate$from.pgs.catalog() once per cohort: two calls give you two libraries that agree today, with nothing checking they still agree after you filter or extend one of them.

From there, each step is one call per study:

discovery$scores$X   <- compute$scores(discovery,   maf.thr = 0.01)
replication$scores$X <- compute$scores(replication, maf.thr = 0.01)

The sample axes stay independent — nothing requires two cohorts to share a sample count, an order, or even a phenotype table with the same columns.

Pooling itself happens on results, not on studies. Fit each cohort separately, then hand the association objects to associate$meta(), which combines them cell by cell — see Meta-Analysis.

Saving and loading

savePolyGenius(study, path) writes a study to a single .pgd file; loadPolyGenius(path) reads it back. The file is self-contained and portable: it carries the model library in full — dictionary, weights, registry, annotations — and the genotypes as pointers to the PLINK files they already live in, never as copies of the genotypes themselves. Move the file without moving those files and everything that does not need genotypes still loads.

If the genotype files moved as well, name each fileset's new directory: loadPolyGenius(path, genotypes = c(cohort = "/new/dir")). The load checks that the files there hold the same samples in the same order, and aborts otherwise.

A load reads and validates the whole payload before it returns, so a truncated or inconsistent file fails at the load call rather than on whatever line first touches the damage. See Persistence and the Backbone for what a round trip keeps and what it deliberately does not.

Before you analyse

Check that the axes are non-empty, that the score layer you are about to name exists and has the right shape, that every covariate resolves through $fetch(), and that the genotype paths are still reachable.