S4-dispatched per-region pipeline for learning TWAS prediction weights. Accepts:
a
QtlDatasetfor individual-level cohort fits;a
QtlSumStatsfor 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, orQtlSumStats. TheMultiStudyQtlDatasetmethod iterates the embedded individual-levelQtlDatasetentries and the optional embeddedQtlSumStats, 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 = argsentries. For QtlSumStats input the default isc("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 afineMappingResult.- 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 overlappingregionwhen 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 withchrom/start/end. Mutually exclusive withtraitId.- cisWindow
For QtlDataset: cis-window (bp) around each trait's genomic position when extracting variants. Required when
traitIdis supplied. Mutually exclusive withregion.- genotypeFilterParam
For QtlDataset: per-call genotype-filter overrides, built with
GenotypeFilterParam. Each field that is set replaces the corresponding construct-timeQtlDatasetslot 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 infineMappingPipeline— 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 flatcvResultdata frame (one row per region).TRUEconcatenates the ranges' genotypes into one joint fit. Ignored for a single-range / cis request.- jointSpecification
Optional joint-fit specification saying which axes
mrmashfits 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", orlist(axes = <axes>, scope = <named list>), wherescope(keysstudy,context,trait) limits the fit to the listed values, each of which must exist indata. 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 aQtlDataset, and on aMultiStudyQtlDatasetonly 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 intolearnTwasWeightsviafittedModelsso SuSiE-family weight methods reuse the prior fit instead of refitting.- twasWeights
Optional
TwasWeightsresume cache. For each requested(study, context, trait, method)tuple already present in this collection, the cachedTwasWeightsRowis 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 thefineMappingResultcache infineMappingPipeline. 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;NULLto skip.- crossValidationParam
Cross-validation settings, built with
CrossValidationParam:folds(default0, 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) andfitFullData(FALSEkeeps only the cross-validation results; defaultTRUE).foldsdefaulted to5before, whilefineMappingPipeline's defaulted to0. The two now agree on0, so an unspecified setting means one thing.Cross-validation holds out samples, so it is not possible on
QtlSumStatsinput; 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(defaultFALSE),r2Threshold,solverandalpha. Stacking reads out-of-fold predictions, soenabled = TRUErequirescrossValidationParamwithfolds >= 2and 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 theQtlDatasetcomponents of aMultiStudyQtlDataset. WhenTRUE, 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 reductionfineMappingPipelineapplies, 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(default10). The effective count ismin(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 settingfineMappingPipelinetakes, which has no"none"— itssusieFitslot is part of the returned object's contract."full"is honoured by the mr.mash engines and, on summary statistics, bylasso,scad,mcpandl0learn, 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
ResidualizationParamand forwarded toresidualizedPhenotypes/residualizedGenotypes:phenotypeCovariatesandgenotypeCovariatesname which covariates to regress out (NULL, the default, uses every available one), andresidualizePhenotype/residualizeGenotypeturn each side off.QtlDataset/MultiStudyQtlDatasetinput only. AQtlSumStatsrun has no covariates to regress out, so the bundle is ignored there rather than refused — unlikecrossValidationParam, 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 (QtlDatasetpath), seeds both the main-process RNG and the parallel method-fitting / cross-validation RNG (BiocParallelRNGseed) 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 outerset.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,mvsusiefollow the multi-trait / multi-context rules of the fine-mappingmvsusiefamily (require at least two traits OR at least two contexts).prsCsruns only on summary statistics, so it is rejected onQtlDatasetinput.enet,dprVb,dprAdaptiveGibbsand the Bayesian alphabet (bayesA,bayesB,bayesC,bayesL,bayesN,bayesR,bLasso) run only on individual-level data, so they are rejected onQtlSumStatsinput.
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)