Skip to contents

S4-dispatched per-region pipeline for learning TWAS prediction weights. Accepts:

  • a QtlDataset for individual-level cohort fits;

  • a QtlSumStats for per-trait RSS fits;

  • a MultiStudyQtlDataset, whose individual-level studies and optional summary-statistic studies are each fitted, then combined.

Usage

twasWeightsPipeline(data, ...)

# S4 method for class 'QtlDataset'
twasWeightsPipeline(
  data,
  methods = "default",
  contexts = NULL,
  traitId = NULL,
  region = NULL,
  cisWindow = NULL,
  genotypeFilterParam = GenotypeFilterParam(),
  jointRegions = FALSE,
  jointSpecification = NULL,
  fineMappingResult = NULL,
  twasWeights = NULL,
  mashPrior = NULL,
  crossValidationParam = CrossValidationParam(),
  ensembleParam = EnsembleParam(),
  estimatePi = TRUE,
  usePCA = FALSE,
  nPCs = 10L,
  fitRetention = c("slim", "none", "full"),
  residualizationParam = ResidualizationParam(),
  weightsDataType = NULL,
  naAction = c("drop", "impute"),
  verbose = 1,
  seed = NULL
)

# S4 method for class 'QtlSumStats'
twasWeightsPipeline(
  data,
  methods = NULL,
  contexts = NULL,
  traitId = NULL,
  jointSpecification = NULL,
  fineMappingResult = NULL,
  twasWeights = NULL,
  fitRetention = c("slim", "none", "full"),
  weightsDataType = NULL,
  verbose = 1L,
  panelFilterParam = PanelFilterParam(),
  crossValidationParam = CrossValidationParam()
)

# S4 method for class 'MultiStudyQtlDataset'
twasWeightsPipeline(
  data,
  methods = "default",
  contexts = NULL,
  traitId = NULL,
  region = NULL,
  cisWindow = NULL,
  genotypeFilterParam = GenotypeFilterParam(),
  panelFilterParam = PanelFilterParam(),
  jointRegions = FALSE,
  jointSpecification = NULL,
  fineMappingResult = NULL,
  twasWeights = NULL,
  mashPrior = NULL,
  estimatePi = TRUE,
  usePCA = FALSE,
  nPCs = 10L,
  fitRetention = c("slim", "none", "full"),
  weightsDataType = NULL,
  naAction = c("drop", "impute"),
  verbose = 1,
  residualizationParam = ResidualizationParam(),
  crossValidationParam = CrossValidationParam(),
  ensembleParam = EnsembleParam(),
  seed = NULL
)

Arguments

data

A QtlDataset, MultiStudyQtlDataset, or QtlSumStats. The MultiStudyQtlDataset method iterates the embedded individual-level QtlDataset entries and the optional embedded QtlSumStats, then rbinds the results.

...

Reserved for method-specific arguments.

methods

A character vector of short method names, a preset string ("default" or "fastDefault"), or a named list of <method>_weights = args entries. For QtlSumStats input the default is c("lasso", "prsCs", "dprGibbs"), which there run lassosum, PRS-CS and SDPR; the fine-mapping methods are not in it, and are requested by name together with a fineMappingResult.

contexts

Optional context selection; NULL (the default) uses every context. A character vector applies to every study: a context no study has is an error, and one only some studies lack is skipped for those, with a warning. A list keyed by study name picks each study's contexts separately (list(A = "liver")); studies it does not name keep all of theirs.

traitId

Optional trait selection; NULL (the default) uses every trait, or those overlapping region when one is given. A character vector applies to every study. A list keyed by study name picks each study's traits separately (list(A = "ENSG1")), and a study's entry may itself be a list keyed by context (list(A = list(brain = "ENSG1"))). Studies and contexts the list does not name keep all of their traits. A joint fit spanning a study's contexts uses every trait selected in any of them, minus the (context, trait) pairs a per-context entry leaves out.

region

Optional variant window for QtlDataset trait selection: a GRanges, a "chr:start-end" string, or a one-row data.frame with chrom/start/end. Mutually exclusive with traitId.

cisWindow

For QtlDataset: cis-window (bp) around each trait's genomic position when extracting variants. Required when traitId is supplied. Mutually exclusive with region.

genotypeFilterParam

For QtlDataset: per-call genotype-filter overrides, built with GenotypeFilterParam. Each field that is set replaces the corresponding construct-time QtlDataset slot for this call only (applied to a validated copy); a field left unset leaves the stored value in place. Variant QC is a property of the data, so these are applied identically here and in fineMappingPipeline — there is deliberately no TWAS-specific variant filter.

jointRegions

For QtlDataset with a multi-range region: FALSE (default) learns weights for each range independently and concatenates them into one entry per (study, context, trait, method); the per-region fits are kept as a named list and per-region CV is recorded as a flat cvResult data frame (one row per region). TRUE concatenates the ranges' genotypes into one joint fit. Ignored for a single-range / cis request.

jointSpecification

Optional joint-fit specification saying which axes mrmash fits across. When NULL (the default) it runs its implicit multi-trait / multi-context fit. Otherwise give a character vector of axes (one spec) or a list of specs. A spec is a character vector drawn from "study", "context" and "trait", or list(axes = <axes>, scope = <named list>), where scope (keys study, context, trait) limits the fit to the listed values, each of which must exist in data. A list-form spec always goes inside the outer list, even on its own. Other methods run per (study, context, trait) as usual. The "study" axis needs summary statistics: it is rejected on a QtlDataset, and on a MultiStudyQtlDataset only the summary-statistics studies take part.

fineMappingResult

Optional QtlFineMappingResult. When supplied, its SuSiE / SuSiE-inf / SuSiE-ash trimmed fits for the matching (study, context, trait) tuples are injected into learnTwasWeights via fittedModels so SuSiE-family weight methods reuse the prior fit instead of refitting.

twasWeights

Optional TwasWeights resume cache. For each requested (study, context, trait, method) tuple already present in this collection, the cached TwasWeightsRow is copied through and the corresponding weight fit is skipped. Only the un-cached method subset is fit; the cached and fresh entries are concatenated in the returned collection. Per-tuple matching mirrors the fineMappingResult cache in fineMappingPipeline. Multivariate dispatch (mvsusie, mr.mash) is unaffected because those methods produce one fit jointly across multiple (context, trait) columns.

mashPrior

Optional MashPrior (or coercible) supplying data-driven prior matrices for mr.mash weight methods; NULL to skip.

crossValidationParam

Cross-validation settings, built with CrossValidationParam: folds (default 0, no CV), threads, samplePartition, maxVariants (cap on the CV design matrix), weightMethods (which methods to cross-validate; unset means every method with non-zero weights) and fitFullData (FALSE keeps only the cross-validation results; default TRUE).

folds defaulted to 5 before, while fineMappingPipeline's defaulted to 0. The two now agree on 0, so an unspecified setting means one thing.

Cross-validation holds out samples, so it is not possible on QtlSumStats input; that method errors if it is set rather than returning results that were never cross-validated.

ensembleParam

SR-TWAS ensemble settings, built with EnsembleParam: enabled (default FALSE), r2Threshold, solver and alpha. Stacking reads out-of-fold predictions, so enabled = TRUE requires crossValidationParam with folds >= 2 and is an error otherwise — it previously returned a result quietly missing its ensemble row.

estimatePi

If TRUE, estimate spike-and-slab sparsity from mr.ash before BGLR / qgg spike-and-slab methods that consume it.

usePCA

Logical (length 1). Individual-level data only: a QtlDataset, or the QtlDataset components of a MultiStudyQtlDataset. When TRUE, each multi-trait context is additionally PCA-reduced and weights are fitted for each top principal component as a trait (topPC1, topPC2, ...), alongside the per-trait rows. The same reduction fineMappingPipeline applies, sharing its scores helper and producing groups the ordinary CV / ensemble / retention machinery handles unchanged. A single-trait context, or one with no usable component, contributes nothing.

nPCs

Integer (length 1). Caps the number of top principal components fitted per context when usePCA = TRUE (default 10). The effective count is min(nPCs, usable traits).

fitRetention

How much of each fit is kept on the entry: "slim" (default), "none" to keep no fit at all, or "full". The same setting fineMappingPipeline takes, which has no "none" — its susieFit slot is part of the returned object's contract.

"full" is honoured by the mr.mash engines and, on summary statistics, by lasso, scad, mcp and l0learn, whose full fit also keeps every candidate's coefficients; every other method keeps its fit whole or not at all. Requesting "full" for a run where no method can honour it is an error, and a mixed run warns rather than dropping it in silence.

residualizationParam

Covariate residualization settings, built with ResidualizationParam and forwarded to residualizedPhenotypes / residualizedGenotypes: phenotypeCovariates and genotypeCovariates name which covariates to regress out (NULL, the default, uses every available one), and residualizePhenotype / residualizeGenotype turn each side off.

QtlDataset / MultiStudyQtlDataset input only. A QtlSumStats run has no covariates to regress out, so the bundle is ignored there rather than refused — unlike crossValidationParam, which errors. Residualization is on by default, so refusing a non-default value would reject the default bundle; CV is off by default, so a non-default value there is an explicit request for something the input cannot do.

weightsDataType

Optional data-type label recorded on every TwasWeightsRow$dataType (e.g. "expression").

naAction

Character. How to handle missing values in the extracted data.

verbose

Verbosity (0 silent, 1 default, 2 includes external package messages).

seed

Integer or NULL. When supplied (QtlDataset path), seeds both the main-process RNG and the parallel method-fitting / cross-validation RNG (BiocParallel RNGseed) for reproducibility under multi-threading. The main-process seed is scoped to the call, so the session RNG is left as it was found. NULL (default) does not seed at all, so an outer set.seed() still governs the main-process draws.

panelFilterParam

For QtlSumStats: LD-reference-panel filters, built with PanelFilterParam. See Panel filters on the RSS path above.

Value

A TwasWeights collection keyed by (study, context, trait, method). The ldReference slot is NULL for individual-level fits and equals the input's ldReference for RSS-derived fits.

Details

Method-restriction rules (enforced):

  • mr.mash, mvsusie follow the multi-trait / multi-context rules of the fine-mapping mvsusie family (require at least two traits OR at least two contexts).

  • prsCs runs only on summary statistics, so it is rejected on QtlDataset input.

  • enet, dprVb, dprAdaptiveGibbs and the Bayesian alphabet (bayesA, bayesB, bayesC, bayesL, bayesN, bayesR, bLasso) run only on individual-level data, so they are rejected on QtlSumStats input.

A QtlSumStats input must have been QC'd via summaryStatsQc first; otherwise an error is raised pointing at that function.

The returned TwasWeights collection's ldReference slot is set automatically: NULL for individual-level fits, the input's ldReference for RSS-derived fits.

Optionally a QtlFineMappingResult may be supplied as a source of pre-fit SuSiE / SuSiE-inf / SuSiE-ash objects; their stored fits are passed through to learnTwasWeights / the RSS sub-pipelines via the fittedModels slot, avoiding a re-fit.

When the supplied result was produced with cross-validation (a CrossValidationParam with folds >= 2), each matching (study, context, trait) entry's cvResult is reused: its fold partition becomes the CV partition (unless samplePartition is given explicitly) and its per-fold out-of-fold predictions/metrics are fed directly into the SR-TWAS ensemble in place of re-fitting those fine-mapping methods here. Non-fine-mapping methods (lasso, enet, ...) are still cross-validated on the same shared partition.

Panel filters on the RSS path

On QtlSumStats input there is no genotype matrix to filter, so panelFilterParam's mafCutoff / macCutoff / imissCutoff are measured against the LD reference panel instead: a variant whose panel genotypes fall below the cutoffs is dropped before the z-scores and LD matrix are built. The thresholds mean the same thing as genotypeFilterParam's on the QtlDataset path (MAC is converted to a MAF equivalent and the stricter of the two applies), so one number carries across input types. The defaults filter nothing.

This discards observed variants, unlike summaryStatsQc(raissParam = ...), which only bounds which variants RAISS will impute.

The genotype scale

Individual-level genotypes reach every method centred and scaled to unit variance, so all the weights are on one scale and the recorded weight_standardized says so. The scaling belongs to the QtlDataset accessors rather than to this pipeline: residualizedGenotypes applies it, governed by ResidualizationParam(scaleResiduals = ) per call with the dataset's own scaleResiduals slot as the default.

It is not cosmetic: measured through pecotmr's own entry points, every method returns coefficients on the scale of the X it was given — exactly so for glmnet / ncvreg / L0Learn / mr.ash / the Bayesian alphabet and for susieWeights() — while the DPR samplers are not scale-equivariant at all, so a uniform rescale of X changes their fit and not merely its units. The chosen scale is also the one the LD correlations and z-scores behind a TWAS Z already use. A summary-statistics fit needs nothing done.

learnTwasWeights still takes weightStandardized: there the caller supplies X themselves, so only they can say.

Examples

data(qtlDatasetExample)
twasWeightsPipeline(qtlDatasetExample, methods = "lasso", cisWindow = 1e6)
#>   Fitting lasso ...
#>   Fitting lasso done in 0.3s
#> TwasWeights: 1 entries
#>   1 studies, 1 contexts, 1 traits, 1 methods
#>   LD reference: NULL (individual-level fit)