Post-processing fine-mapping results
pecotmr authors
2026-10-22
Source:vignettes/finemapping-postprocessing.Rmd
finemapping-postprocessing.RmdOverview
fineMappingPipeline()
returns a QtlFineMappingResult or
GwasFineMappingResult. This vignette covers what you do
with one afterwards: summarising credible sets, comparing results across
studies, intersecting QTL with GWAS, and reaching the underlying SuSiE
quantities when the summaries are not enough.
Every function here reads an existing result — none of them refit anything.
library(pecotmr)
data(qtlFineMappingExample, gwasFineMappingExample)
data(qtlFineMappingLbfExample, gwasFineMappingLbfExample)
data(qtlFineMappingPairedExample, qtlDatasetExample)
qtlFineMappingExample## QtlFineMappingResult: 1 entries
## 1 studies, 1 contexts, 1 traits, 1 methods
## LD reference: NULL (individual-level fit)The two tables
A fine-mapping result carries two flat summaries, and which one you want depends on whether your unit of interest is a variant or a credible set.
topLoci() is one range per retained
variant, with its PIP, effect size, and credible-set membership at each
coverage level carried in mcols(). Take a table with
as.data.frame():
## [1] 15
head(as.data.frame(loci)[, c(
"variant_id",
"pip",
"cs_95",
"cs_95_purity",
"posterior_effect",
"posterior_se"
)], 4)## variant_id pip cs_95 cs_95_purity posterior_effect
## 1 chr22:32596246:T:C 0.14449697 susie_0 0.0000000 -0.012846434
## 2 chr22:32724254:T:C 0.04426635 susie_1 0.6255032 0.006013377
## 3 chr22:32724508:C:A 0.04426635 susie_1 0.6255032 0.006013377
## 4 chr22:32728910:G:A 0.03949625 susie_1 0.6255032 0.005349899
## posterior_se
## 1 0.03244071
## 2 0.02838201
## 3 0.02838201
## 4 0.02680179posterior_effect is the posterior mean effect — the same
quantity a TWAS weight is, and literally what
susieWeights() returns — and posterior_se its
standard error. They are deliberately not called
beta/se: marginalEffects() publishes the
marginal single-variant estimate under those names, so
sharing them would let the two views be joined into a silent mix of
posterior and marginal effects.
Both are on the standardized genotype scale, because the
QtlDataset accessors centre and scale every genotype column
to unit variance before fitting — residualizedGenotypes()
does it, governed by
ResidualizationParam(scaleResiduals = FALSE) per call or
the dataset’s own scaleResiduals slot. The published
weight_standardized records whichever way it resolved. That
is what makes a posterior effect comparable with a TWAS weight from any
other method — they are all on the scale of the X the
fitter was given, and on the same scale as the LD correlations and
z-scores a TWAS Z is computed from. The published
weight_standardized column records it.
cs_95 names the credible set the variant belongs to, or
is NA when it belongs to none. The cs_70 /
cs_50 columns give the same membership at looser coverage,
so a variant can be in a 95% set but not a 50% one.
credibleSetSummary() is one row per
credible set instead:
as.data.frame(credibleSetSummary(qtlFineMappingExample))[, c(
"cs",
"coverage",
"n_variants",
"purity_min",
"cs_log10bf",
"lead_variant",
"lead_pip"
)]## cs coverage n_variants purity_min cs_log10bf lead_variant
## 1 susie_1 0.95 42 0.6255032 16.22962 chr22:32745710:T:C
## lead_pip
## 1 0.0769631purity_min is the smallest absolute correlation between
any pair of variants in the set — the standard measure of whether a
credible set is a single signal or an artefact of a few uncorrelated
variants being lumped together. cs_log10bf says how much
the data prefer that effect being present.
A credible set with low purity_min deserves suspicion
before interpretation, which is why the pipeline exposes
minPurity as a filter rather than dropping such sets
silently.
Merging overlapping credible sets
A multi-effect SuSiE fit can return credible sets that overlap,
especially when two effects are in strong LD.
mergeSusieCs() collapses those into one row per variant:
credibleSetNames labels every set the variant belongs to
(comma-joined when it belongs to more than one), and maxPip
/ medianPip summarise its PIP across them. The example fit
has a single set, so each label names just that one:
merged <- mergeSusieCs(qtlFineMappingExample, coverage = 0.95)
dim(merged)## [1] 14 4
head(as.data.frame(merged), 4)## variant_id credibleSetNames maxPip medianPip
## 1 chr22:32724254:T:C cs_1_1 0.04426635 0.04426635
## 2 chr22:32724508:C:A cs_1_1 0.04426635 0.04426635
## 3 chr22:32728910:G:A cs_1_1 0.03949625 0.03949625
## 4 chr22:32729588:T:C cs_1_1 0.03451015 0.03451015Use it when downstream code needs a variant to belong to at most one set — overlap-counting, for instance, would otherwise double-count.
Intersecting QTL and GWAS
overlapTopLoci() joins two fine-mapping results on
shared variants, prefixing each side’s columns. It is the cheap first
pass before colocalization: if the two
credible sets share no variants at all, there is nothing for a coloc
model to work with.
overlap <- overlapTopLoci(qtlFineMappingLbfExample, gwasFineMappingLbfExample)
nrow(overlap)## NULL
head(as.data.frame(overlap)[, c(
"variant_id",
"qtl_pip",
"gwas_pip",
"qtl_cs_95",
"gwas_cs_95"
)], 4)## variant_id qtl_pip gwas_pip qtl_cs_95 gwas_cs_95
## 1 chr22:10416669:A:T 0.03101916 0.03101916 susie_0 susie_0
## 2 chr22:10424460:A:T 0.02609678 0.02609678 susie_0 susie_0
## 3 chr22:10443846:A:T 0.02806086 0.02806086 susie_0 susie_0
## 4 chr22:10456378:A:T 0.02518179 0.02518179 susie_0 susie_0The variants worth attention are those in a credible set on both sides:
## NULL
table(both$qtl_cs_95, both$gwas_cs_95)##
## susie_0
## susie_0 100Here that filter removes nothing, because topLoci in
this result retains only credible-set variants — every shared variant is
in one set on each side. The filter earns its place on results where
topLoci also carries high-PIP variants outside any set,
which is the usual case once pipCutoff is in play.
A non-empty overlap is not evidence of a shared causal variant,
however large: two independent signals in the same LD block overlap too.
That question is what colocPipeline() answers.
Combining results
combineFineMappingResults() stacks several results into
one collection, which is how per-region or per-chromosome runs are
assembled after parallel execution:
combined <- combineFineMappingResults(
qtlFineMappingExample,
qtlFineMappingLbfExample
)
nrow(combined)## [1] 3
data.frame(
study = combined$study,
context = combined$context,
trait = combined$trait,
method = combined$method
)## study context trait method
## 1 study_1 context_1 gene_1 susie
## 2 test_study context1 ENSG00000283047 susie
## 3 test_study context2 ENSG00000283047 susieThe inputs keep their own identity columns, so a combined object
still distinguishes which study and method each row came from. Note that
the identity columns are read with $:
as.data.frame() on a collection flattens its
variants, not its per-row metadata.
Between-credible-set correlation
computeCsCorrelation() reports how correlated a fit’s
credible sets are with each other. It is computed on demand rather than
stored, because it needs LD — and the LD lives with the data, not with
the fit:
computeCsCorrelation(qtlFineMappingPairedExample, qtlDatasetExample)## L2 L3
## L2 1.00000000 -0.08404977
## L3 -0.08404977 1.00000000The ldSource argument takes a QtlDataset,
QtlSumStats or GwasSumStats — whichever object
the fit was derived from. It must cover the fit’s region: a source whose
variants do not reach the fit’s window cannot supply the correlations,
and the function says so rather than returning something empty.
Near-zero off-diagonal entries, as here, mean the credible sets are capturing distinct signals. Strongly correlated sets suggest the fit split one signal in two.
Reaching the SuSiE quantities
When the flat tables are not enough — you want the per-effect posterior, or to re-derive PIPs under a different assumption — the underlying fit is available.
lbf() returns the log Bayes factors, one row per variant
and one column per single effect:
lbf <- lbf(qtlFineMappingLbfExample)
unit <- lbf[lbf$trait == lbf$trait[1] & lbf$context == lbf$context[1], ]
lbfMatrix <- as.matrix(unit[, grep("^lbf_L", names(unit))])
dim(lbfMatrix)## [1] 200 5lbfToAlpha() converts those to single-effect inclusion
probabilities. It normalises across each row, so the
matrix has to be oriented with effects as rows — transpose the
variants-by-effects layout above:
alpha <- lbfToAlpha(t(lbfMatrix))
dim(alpha)## [1] 5 200
rowSums(alpha) # each single effect sums to 1 across variants## lbf_L1 lbf_L2 lbf_L3 lbf_L4 lbf_L5
## 1 1 1 1 1Getting that orientation wrong is quiet rather than loud: normalising across effects for each variant also produces a matrix whose rows sum to 1, just not one that means anything.
The LBF matrix is also what colocPipeline() consumes,
which is why fineMappingPipeline(..., keepLbf = TRUE) is
required before colocalization.
susieFit() returns the raw per-unit SuSiE object when
you need something the accessors do not expose:
## [1] "susie"
names(fit)## [1] "pip" "sets" "sets_secondary"
## [4] "alpha" "lbf_variable" "V"
## [7] "niter" "max_L" "n_effects"
## [10] "X_column_scale_factors" "mu" "mu2"
## [13] "sigma2" "null_index" "converged"
## [16] "pi"alpha and lbf_variable are the per-effect
matrices, sets the credible sets before any purity
filtering, and V the estimated prior variances — everything
the flat tables above are derived from.
How the tables are built
The summaries above are not stored by the fitter — they are derived
from the raw SuSiE fit by a short chain that
fineMappingPipeline() runs for you:
postprocessFinemappingFits() # per-method fits -> a post object
-> computeCsTables() # credible sets from alpha, with purity
-> buildTopLoci() # the flat one-row-per-variant table
-> formatFinemappingOutput() # assemble the result payloadCall them directly when you have SuSiE fits from elsewhere and want
pecotmr’s tables over them, or when you need to change how credible sets
are formed — computeCsTables() is where the coverage levels
and purity threshold are applied, so re-running it with a different
coverage gives a different top_loci without
refitting anything.
Two smaller extractors sit at the row level and are what
credibleSetSummary() is assembled from:
-
extractCsInfo(fmRow, csNames, topLociTable, ldSource)— per-credible-set detail for one result row, including the between-variant correlations that givepurity_min, which is why it needs anldSource. -
extractTopPipInfo(fmRow, sumstats)— the lead variant and its statistics, joined against summary statistics.
susieResult() unwraps a per-condition fit list into the
single SuSiE object the other helpers expect.
Writing results out
Two generics cover every published format. exportToVcf()
serialises a fine-mapping result — or summary statistics, or TWAS
weights — to GWAS-VCF, the interchange format most external tools read:
one record per variant, one sample per study or fitted model, and the
statistics as FORMAT fields.
exportToVcf(
qtlFineMappingExample,
outputPath = "finemapping.vcf.gz",
context = "context_1"
)It needs the VariantAnnotation package, which is a
suggested dependency rather than a hard one.
exportToBed() writes the published long-format tables
instead, one row per published observation, bgzipped and tabix-indexed
when the path ends in .gz:
# the top-loci table: the posterior view joined to the marginal one
exportToBed(qtlFineMappingExample, "top_loci.bed.gz")
# every remaining column the view carries, after the published schema
exportToBed(qtlFineMappingExample, "top_loci_full.bed.gz", columns = "all")
# the fSuSiE table: one row per (context, credible set), carrying the fitted
# curve at both resolutions -- the wavelet grid and its reinterpolation onto
# the trait positions. `levels` keeps only one of the two.
exportToBed(fsusieFineMappingExample, "fsusie.bed.gz", schema = "fsusie")
exportToBed(fsusieFineMappingExample, "fsusie_grid.bed.gz",
schema = "fsusie", levels = "wavelet")Both accept a list of results and fold it into one
sorted, indexed output, so a genome’s worth of per-region RDSes needs no
separate concatenate step.
Reading results back
readFineMappingResult() is the inverse, and it is also
the way into pecotmr for a result produced somewhere else: it reads
either format and rebuilds the collection, keyed by whatever
(study, context, trait, method) tuples the file holds.
res <- readFineMappingResult("top_loci.bed.gz")
# only the variants of one region, from a tabix-indexed file
res <- readFineMappingResult("top_loci.bed.gz", region = "chr22:32e6-33e6")It requires almost nothing: coordinate-encoding variant ids and a
PIP. Every other column is optional, because each consumer
needs a different few — Mendelian randomization the marginal effect and
its standard error, colocalization the log Bayes factors — and a table
converted from another tool will carry some and not others. A consumer
reached without what it needs names the column that is missing rather
than quietly producing nothing.
Two arguments exist for files that were assembled by hand.
studyName, context, trait and
method fill in identity columns the file does not record (a
value that disagrees with the file is an error, not an
override), and weightStandardized says which genotype scale
a published WGT is on — required when the file does not
say, because the two scales differ by a per-variant,
allele-frequency-dependent factor and neither is recognisable on its
own.
readTwasWeights() does the same for a weights file. Note
that a fine-mapping table reads as weights too: the posterior effect
is the prediction weight, which is why both schemas publish it
as WGT.
Next steps
- Fine-mapping — producing the results this vignette consumes.
- Colocalization — the model that answers whether an overlap reflects a shared causal variant.
- TWAS weights — reusing a fine-mapping fit as a weight method.
Session info
## R version 4.5.3 (2026-03-11)
## Platform: x86_64-conda-linux-gnu
## Running under: Ubuntu 24.04.5 LTS
##
## Matrix products: default
## BLAS/LAPACK: /home/runner/work/pecotmr/pecotmr/.pixi/envs/default/lib/libopenblasp-r0.3.34.so; LAPACK version 3.12.0
##
## locale:
## [1] LC_CTYPE=C.UTF-8 LC_NUMERIC=C LC_TIME=C.UTF-8
## [4] LC_COLLATE=C.UTF-8 LC_MONETARY=C.UTF-8 LC_MESSAGES=C.UTF-8
## [7] LC_PAPER=C.UTF-8 LC_NAME=C LC_ADDRESS=C
## [10] LC_TELEPHONE=C LC_MEASUREMENT=C.UTF-8 LC_IDENTIFICATION=C
##
## time zone: Etc/UTC
## tzcode source: system (glibc)
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] pecotmr_0.9.3
##
## loaded via a namespace (and not attached):
## [1] DBI_1.3.0 bitops_1.1-0
## [3] httr2_1.3.0 rlang_1.3.0
## [5] magrittr_2.0.5 otel_0.2.0
## [7] matrixStats_1.5.0 susieR_0.16.6
## [9] compiler_4.5.3 RSQLite_3.53.3
## [11] systemfonts_1.3.2 vctrs_0.7.3
## [13] quadprog_1.5-8 stringr_1.6.0
## [15] pkgconfig_2.0.3 crayon_1.5.3
## [17] fastmap_1.2.0 backports_1.5.1
## [19] dbplyr_2.6.0 XVector_0.50.0
## [21] Rsamtools_2.26.0 rmarkdown_2.32
## [23] tzdb_0.5.0 ragg_1.5.2
## [25] purrr_1.2.2 bit_4.6.0
## [27] xfun_0.61 Rfast_2.1.5.2
## [29] MultiAssayExperiment_1.36.2 cachem_1.1.0
## [31] jsonlite_2.0.0 blob_1.3.0
## [33] DelayedArray_0.36.1 reshape_0.8.10
## [35] tictoc_1.2.1 BiocParallel_1.44.0
## [37] irlba_2.4.1 parallel_4.5.3
## [39] R6_2.6.1 bslib_0.12.0
## [41] stringi_1.8.9 RColorBrewer_1.1-3
## [43] GenomicRanges_1.62.1 jquerylib_0.1.4
## [45] numDeriv_2016.8-1.1 Rcpp_1.1.2
## [47] Seqinfo_1.0.0 SummarizedExperiment_1.40.0
## [49] knitr_1.52 readr_2.2.0
## [51] IRanges_2.44.0 BiocBaseUtils_1.12.0
## [53] Matrix_1.7-6 tidyselect_1.2.1
## [55] abind_1.4-8 yaml_2.3.12
## [57] codetools_0.2-20 curl_8.0.0
## [59] metafor_5.0-1 lattice_0.23-1
## [61] tibble_3.3.1 plyr_1.8.9
## [63] Biobase_2.70.0 withr_3.0.3
## [65] S7_0.2.2 evaluate_1.0.5
## [67] archive_1.1.13 desc_1.4.3
## [69] RcppParallel_5.1.11-2 BiocFileCache_3.0.0
## [71] Biostrings_2.78.0 filelock_1.0.3
## [73] pillar_1.11.1 MatrixGenerics_1.22.0
## [75] metadat_1.6-0 checkmate_2.3.3
## [77] stats4_4.5.3 generics_0.1.4
## [79] vroom_1.7.1 mathjaxr_2.0-0
## [81] S4Vectors_0.48.1 hms_1.1.4
## [83] ggplot2_4.0.3 scales_1.4.0
## [85] glue_1.8.1 tools_4.5.3
## [87] fs_2.1.0 grid_4.5.3
## [89] PlinkMatrix_1.0.0 tidyr_1.3.2
## [91] nlme_3.1-170 cli_3.6.6
## [93] zigg_0.0.2 textshaping_1.0.5
## [95] mixsqp_0.3-54 S4Arrays_1.10.1
## [97] dplyr_1.2.1 gtable_0.3.6
## [99] sass_0.4.10 digest_0.6.39
## [101] BiocGenerics_0.56.0 SparseArray_1.10.10
## [103] htmlwidgets_1.6.4 farver_2.1.2
## [105] memoise_2.0.1 htmltools_0.5.9
## [107] pkgdown_2.2.1 lifecycle_1.0.5
## [109] bit64_4.8.6 MASS_7.3-66