Skip to contents

Fine-Mapping Pipeline

Source: R/fineMappingPipeline.R
fineMappingPipeline.Rd

S4-dispatched per-region fine-mapping entry point that replaces the deprecated univariateAnalysisPipeline, multivariateAnalysisPipeline, rssAnalysisPipeline, and susieRssPipeline pipelines. Accepts:

  • a QtlDataset for individual-level cohort fits (per-context / per-trait univariate SuSiE; joint multi-trait or multi-context mvSuSiE; joint multi-trait fSuSiE per context);

  • a MultiStudyQtlDataset which recurses through each embedded QtlDataset per study and processes the optional embedded QtlSumStats via the sumstat method;

  • a QtlSumStats for per-trait SuSiE-RSS fits and per-(study, trait) multi-context mvSuSiE-RSS fits;

  • a GwasSumStats for per-(study, LD-block) SuSiE-RSS fine-mapping (used by qtlEnrichmentPipeline downstream).

Method tokens are unified across input classes; auto-dispatch picks the individual-level vs RSS implementation based on the input class. The supported tokens are:

susie

susieR::susie with unmappable_effects = "none" on individual-level input; susieR::susie_rss (same) on RSS.

susieInf

unmappable_effects = "inf" variant of the same.

susieAsh

unmappable_effects = "ash" variant of the same.

ser

Single-effect regression via susieR::susie_ser on summary statistics (z, n); LD-free (no R, no L), so distinct from susie with L = 1. Sumstat input only (QtlSumStats / GwasSumStats, or the sumstat side of a MultiStudyQtlDataset); rejected on individual-level QtlDataset.

mvsusie

mvsusieR::mvsusie on individual- level input (requires multi-trait OR multi-context Y), mvsusieR::mvsusie_rss on sumstat input (requires multi-context within a single (study, trait) group). Errors on GwasSumStats input.

fsusie

fsusieR::susiF joint multi-trait fit per context. Individual-level only; errors on any SumStats input.

mrmash

Always rejected here. mr.mash is a TWAS weight-oriented method and lives in twasWeightsPipeline.

Usage

fineMappingPipeline(data, ...)

# S4 method for class 'QtlDataset'
fineMappingPipeline(
  data,
  methods,
  contexts = NULL,
  traitId = NULL,
  region = NULL,
  cisWindow = NULL,
  genotypeFilterParam = GenotypeFilterParam(),
  jointRegions = FALSE,
  jointSpecification = NULL,
  initializeWithSusieInf = TRUE,
  credibleSetParam = CredibleSetParam(),
  fineMappingResult = NULL,
  crossValidationParam = CrossValidationParam(),
  residualizationParam = ResidualizationParam(),
  signalScreenParam = SignalScreenParam(),
  usePCA = FALSE,
  nPCs = 10L,
  seed = NULL,
  mrmashPrior = NULL,
  dataDrivenPriorWeightsCutoff = 1e-10,
  naAction = c("drop", "impute"),
  verbose = 1,
  fitRetention = c("slim", "full"),
  dataType = NULL
)

# S4 method for class 'MultiStudyQtlDataset'
fineMappingPipeline(
  data,
  methods,
  contexts = NULL,
  traitId = NULL,
  region = NULL,
  cisWindow = NULL,
  genotypeFilterParam = GenotypeFilterParam(),
  panelFilterParam = PanelFilterParam(),
  jointRegions = FALSE,
  jointSpecification = NULL,
  initializeWithSusieInf = TRUE,
  credibleSetParam = CredibleSetParam(),
  fineMappingResult = NULL,
  mrmashPrior = NULL,
  dataDrivenPriorWeightsCutoff = 1e-10,
  crossValidationParam = CrossValidationParam(),
  residualizationParam = ResidualizationParam(),
  signalScreenParam = SignalScreenParam(),
  susieRssParam = SusieRssParam(),
  usePCA = FALSE,
  nPCs = 10L,
  seed = NULL,
  naAction = c("drop", "impute"),
  verbose = 1,
  fitRetention = c("slim", "full"),
  dataType = NULL
)

# S4 method for class 'QtlSumStats'
fineMappingPipeline(
  data,
  methods,
  contexts = NULL,
  traitId = NULL,
  jointSpecification = NULL,
  initializeWithSusieInf = TRUE,
  credibleSetParam = CredibleSetParam(),
  fineMappingResult = NULL,
  mrmashPrior = NULL,
  dataDrivenPriorWeightsCutoff = 1e-10,
  verbose = 1,
  fitRetention = c("slim", "full"),
  susieRssParam = SusieRssParam(),
  panelFilterParam = PanelFilterParam(),
  crossValidationParam = CrossValidationParam(),
  residualizationParam = ResidualizationParam(),
  dataType = NULL
)

# S4 method for class 'GwasSumStats'
fineMappingPipeline(
  data,
  methods,
  initializeWithSusieInf = TRUE,
  credibleSetParam = CredibleSetParam(),
  fineMappingResult = NULL,
  verbose = 1,
  fitRetention = c("slim", "full"),
  susieRssParam = SusieRssParam(),
  panelFilterParam = PanelFilterParam(),
  crossValidationParam = CrossValidationParam(),
  residualizationParam = ResidualizationParam(),
  dataType = NULL
)

Arguments

data

A QtlDataset, MultiStudyQtlDataset, QtlSumStats, or GwasSumStats.

...

Reserved for future per-method arguments.

methods

Method specification. Accepts either:

  • A character vector of method tokens, e.g. c("susie", "susieInf", "mvsusie") (any subset of c("susie", "susieInf", "susieAsh", "mvsusie", "fsusie"), subject to per-class compatibility).

  • A FineMappingMethodsParam record keyed by method token, where each value is a list of per-method kwargs to splice into the underlying fitter, e.g. FineMappingMethodsParam(susie = list(L = 1, refine = FALSE), mvsusie = list(max_iter = 500)). Each entry is checked against the engine that receives it. Mirrors the convention of twasWeightsPipeline's methods argument. User-supplied kwargs override the capability-table defaults and any base / chained args set by the pipeline (e.g. you can override model_init even when fitting from a susieInf chain).

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

Per-call overrides of the QtlDataset's own filters, built with GenotypeFilterParam. A field left unset keeps the dataset's construct-time value, so GenotypeFilterParam(mafCutoff = 0) pins the cutoff at zero while GenotypeFilterParam() changes nothing. Applies to the QtlDataset method only.

jointRegions

For QtlDataset with a multi-range region: FALSE (default) fits each range independently and merges the per-range results into one entry per (study, context, trait, method) – the merged susieFit is a named list of per-region fits and credible-set labels are renumbered to stay unique. TRUE concatenates the ranges' genotypes into one joint fit. Ignored for a single-range / cis (traitId + cisWindow) request.

jointSpecification

Optional joint-fit specification saying which axes mvsusie / fsusie fit across. When NULL (the default) they fit across traits if the input has two or more, else across contexts. 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.

initializeWithSusieInf

Logical. When susieInf is in methods alongside susie and/or susieAsh, whether the SuSiE-inf fit initialises the chained downstream method(s). Default TRUE. A choice between methods rather than a property of the credible sets, which is why it stands alone.

credibleSetParam

How credible sets are built and reported, built with CredibleSetParam: coverage (default 0.95), secondaryCoverage (c(0.7, 0.5)), signalCutoff (the PIP cutoff for top-loci selection, 0.025), minAbsCorr (0.8) and medianAbsCorr (NULL) for purity — a set is kept if it passes either, OR-logic — and includeAllCs. perCsColumns adds the per-set variant-level columns to topLoci. maxNumSingleEffects bounds how many credible sets can exist (default 10) and maxNumSingleEffectsGreedy the greedy loop of the SuSiE-inf refinement.

Most fields are read by postprocessFinemappingFits, so a stored fit can be re-summarized under different settings without refitting. maxNumSingleEffects / maxNumSingleEffectsGreedy are the exception: they are seeded onto each SuSiE-family token's own arguments, so changing them needs a refit.

fineMappingResult

Optional existing QtlFineMappingResult or GwasFineMappingResult to use as a resume cache; tuples already present are not refit.

crossValidationParam

Cross-validation settings, built with CrossValidationParam: folds (default 0, no CV), numThreads for the per-fold refits, and samplePartition to reuse a fixed partition across methods. When folds > 1 each method is refit on every fold's training samples and used to predict the held-out ones; the partition, the out-of-fold predictions and the metrics are stored on each FineMappingRow's cvResult slot. fitFullData = FALSE keeps only that: each row carries its cvResult and no full-data fit, variants or top loci. A fineMappingResult resume cache is only reused where it carries what the run asks for – a full-data fit unless the run is cross-validation-only, and a cvResult when the run cross-validates – and the other tuples are refit.

Cross-validation holds out samples, so it is not possible on QtlSumStats / GwasSumStats input. Those methods accept the argument and error if it is set, rather than accepting it and quietly returning results that were never cross-validated.

maxVariants and weightMethods are honoured by twasWeightsPipeline only and are ignored here.

residualizationParam

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

Only meaningful for QtlDataset / MultiStudyQtlDataset input. A QtlSumStats / GwasSumStats 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.

The marginal univariate effects stored on each FineMappingRow obey the same choice as the SuSiE fit itself — they are computed against the same residualized X / Y.

signalScreenParam

Individual-level single-effect (SER) pre-screen applied to each residualized (X, y) block before a full fit, built with SignalScreenParam. A susie model with L = 1 is fit and the block is skipped unless the chosen metric exceeds its cutoff:

  • pip — the maximum PIP. A negative cutoff uses the adaptive 3 / nVariants threshold.

  • absZ — the maximum marginal |z|, which needs no fit.

  • bf / logBf — the maximum per-variant Bayes factor or log Bayes factor from the L = 1 fit.

Only one metric may be enabled, which SignalScreenParam enforces when you build it. Unset (the default) screens nothing. The summary-statistics analog lives in summaryStatsQc.

usePCA

Logical (length 1). Individual-level data only: a QtlDataset, or the QtlDataset components of a MultiStudyQtlDataset. When TRUE (default FALSE), each multi-trait context's PCA-reduced phenotype is fine-mapped with univariate SuSiE on its top principal components (ports the legacy fsusie.R susie_on_top_pc). Each PC becomes a pseudo-trait row keyed trait = "topPC\{i\}", method = "susie". Single-trait contexts have no PCA and are skipped.

nPCs

Integer (length 1). Individual-level data only. Caps the number of top principal components fine-mapped per context when usePCA = TRUE (default 10). The effective count is min(nPCs, usable traits).

seed

Optional integer. When non-NULL, withr::local_seed(seed) is called at the start of the call for reproducible fits; the global RNG state is restored on return. Default NULL (no seeding).

mrmashPrior

Optional TwasWeights from a previous mr.mash twasWeightsPipeline run, mined for the mvSuSiE data-driven prior: the retained mrmash fits supply the prior mixture, and their cross-validation payload supplies honest per-fold priors. NULL uses the canonical prior instead.

dataDrivenPriorWeightsCutoff

Numeric or NULL. Cutoff below which data-driven prior weights are zeroed; NULL disables the cutoff.

naAction

Character. How to handle missing values in the extracted phenotype/genotype data.

verbose

Verbosity (0 silent, 1 default). Default 1.

fitRetention

How much of each fit the susieFit slot keeps: "slim" (default) a trimmed view, or "full" the whole susie() return, so susieFit() and non-default-coverage credibleSets() queries can read the full posterior matrices. The per-variant topLoci table is fully populated either way — its per-credible-set columns are governed by credibleSetParam's perCsColumns, not by this.

dataType

Optional assay / modality label for the phenotype the fit was run on ("bulk_rnaseq", "methylation"), recorded as a per-row provenance column and published by exportToBed / exportToVcf. NULL (default) omits the column. It cannot be derived — the input carries the values, not what they measure.

panelFilterParam

LD-reference-panel filters for the summary-statistics and GWAS methods, built with PanelFilterParam. See Panel filters on the RSS path above.

susieRssParam

Summary-statistics solver settings for QtlSumStats / GwasSumStats input, built with SusieRssParam:

  • serFallback (default FALSE) — after each multi-effect SuSiE-RSS fit, read susieR's finite-sample R diagnostics and, when fit$R_finite_diagnostics$R_reliability_flag is TRUE, report the single-effect (ser_model) result for that region instead.

  • keepFullFit (default "fallback") — which pre-fallback multi-effect fits to retain: only the regions that fell back, "all", or "none". Meaningful only with serFallback. The retained fit and the decision are read via susieFit(res)$multiEffectFit, $R_reliability_flag and $serFallbackUsed.

  • rFinite — finite-sample size for susieR's R_finite correction. NULL (default) uses susieR's own, except when a finite/EB mode is active (serFallback or rMismatch != "none"), where it falls back to the LD panel's nSamples(ldReference).

  • rMismatch (default "none") — LD-mismatch correction forwarded as R_mismatch: also "eb" or "eb_mix".

  • control — options for susieR::susie_rss_control(), built with SusieRssControlOptions. This was a bare named list, which meant a misspelled option was dropped silently; the constructor checks the names against that function's live formals.

Value

A QtlFineMappingResult keyed by (study, context, trait, method) – for GwasSumStats input, a GwasFineMappingResult keyed by (study, method) and range. The ldReference slot is set automatically: NULL for individual-level (QtlDataset / all-individual-level MultiStudyQtlDataset) fits, the input's ldReference for RSS-derived fits.

Chained initialisation

When susieInf is requested alongside susie and/or susieAsh and initializeWithSusieInf = TRUE (the default), the SuSiE-inf fit is computed first and used as initialisation for the SuSiE / SuSiE-ash fits, mirroring the legacy univariateAnalysisPipeline / susieRssPipeline chained init behaviour. SuSiE-inf is dropped from the final result when the caller did not explicitly request it (only used as init).

QC contract

All QtlSumStats and GwasSumStats inputs must have been QC'd via summaryStatsQc; the pipeline errors on inputs where length(qcInfo(x)) == 0L. summaryStatsQc also drops variants absent from the ldReference, so by the time per-entry processing runs every variant is guaranteed to be present in the LD panel and a local LD matrix can be built with extractBlockGenotypes + computeLd("sample").

Optional resume cache

Supplying fineMappingResult of an existing fine-mapping result skips re-fitting any (study, context, trait, method) tuple that already has a matching row; cached entries are merged with the newly-fit entries in the returned collection.

Intentional behaviours dropped from the pre-stub pipelines

The four pre-stub pipelines (univariateAnalysisPipeline / multivariateAnalysisPipeline / rssAnalysisPipeline / susieRssPipeline) carried several behaviours that are deliberately not ported here:

  • TWAS weights computation (twasWeights = TRUE path): lives in twasWeightsPipeline now.

  • Filtering knobs (now genotypeFilterParam / panelFilterParam, ldReferenceMetaFile): individual- level QC lives on the QtlDataset constructor; sumstat QC lives in summaryStatsQc(). No filtering happens inside this pipeline.

  • Diagnostic re-analysis paths (singleEffect / bayesianConditionalRegression reanalysis on the RSS path): these are not exposed as dedicated method tokens. Callers who want a single-effect fit can request it via per-method kwargs, e.g. methods = FineMappingMethodsParam(susie = list(L = 1)) methods parameter).

  • loadRssData and explicit ldReferenceMetaFile arguments: the new QtlSumStats / GwasSumStats carry the (already-QC'd) sumstats and ldReference directly.

  • Verbose methodName suffixing (e.g. "susie_rss_NO_QC", "susie_rss_SLALOM_RAISS_imputed"): the method column on the returned collection carries the bare token ("susie", "susieInf", "mvsusie", ...) only. QC provenance is recorded on the sumstats' qcInfo.

  • An entry that summaryStatsQc(pipCutoffToSkip = ...) screened out (recorded as qcInfo$entryAudit[[i]]$pipScreenSkipped, and emptied to 0 variants) is skipped, not fit: it produces no row and a message with the screen reason. An all-screened collection yields a valid empty result rather than an error.

Panel filters on the RSS path

On QtlSumStats / GwasSumStats 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. Use it when the panel cannot support the LD estimate a rare variant would need.

The genotype scale

Individual-level genotypes reach every fit centred and scaled to unit variance, so a posterior effect — which is the prediction weight — is on the standardized scale, and the published weight_standardized records it. That scaling belongs to the QtlDataset accessors, not to this pipeline: it is residualizedGenotypes that applies it, governed by ResidualizationParam(scaleResiduals = ) per call with the dataset's own scaleResiduals slot as the default. One mechanism, one place, shared with twasWeightsPipeline and colocboostPipeline. An RSS fit needs nothing done: its likelihood is written in standardized units already.

It is what makes weights comparable across methods. Measured through pecotmr's own entry points, every TWAS method returns coefficients on the scale of the X it was given, and the DPR samplers are not even scale-equivariant — a uniform rescale of X changes their fit rather than just its units. One input scale removes the question, and it is the scale the LD correlations and z-scores behind a TWAS Z are already on.

Examples

data(qtlDatasetExample)
fineMappingPipeline(qtlDatasetExample, methods = "susie", cisWindow = 1e6)
#> Fitting susie for (context='brain', trait='ENSG_example') ...
#> QtlFineMappingResult: 1 entries
#>   1 studies, 1 contexts, 1 traits, 1 methods
#>   LD reference: NULL (individual-level fit)