Fine-mapping with pecotmr
pecotmr authors
2026-11-02
Source:vignettes/fine-mapping.Rmd
fine-mapping.RmdOverview
fineMappingPipeline() is a single pipeline that can
accept any of the following inputs:
| Input class | What it represents | Methods supported |
|---|---|---|
QtlDataset |
Single-study individual-level data | susie, susieInf, susieAsh, mvsusie, fsusie |
MultiStudyQtlDataset |
Multiple QtlDataset objects (optionally with an
embedded QtlSumStats) |
susie, susieInf, susieAsh, mvsusie, fsusie (individual-level only),
ser (via an embedded QtlSumStats) |
QtlSumStats |
QTL summary statistics annotated by
(study, context, trait)
|
susie, susieInf, susieAsh, mvsusie, ser |
GwasSumStats |
GWAS summary statistics annotated by study
|
susie, susieInf, susieAsh, ser |
Method arguments are the same across input classes;
fineMappingPipeline() routes methods = "susie"
to susieR::susie for individual-level data and to
susieR::susie_rss for summary-statistics inputs
automatically.
The output is a QtlFineMappingResult (a
GwasFineMappingResult for GwasSumStats input):
a GRangesList-subclass collection, one element per
fine-mapped tuple.
Bundled inputs
data(
qtlDatasetExample,
qtlSumStatsExample,
gwasSumStatsS4Example,
multiStudyQtlDatasetExample
)
qtlDatasetExample## QtlDataset for study 'study1'
## 1 context(s): brain
## 1 unique traits across contexts
## Genotypes: plink1 @ pecotmr://extdata/toy_canonical
## Genotype covariates: 0 cols
## Samples: 165The *_sumstats_example objects ship with non-empty
qcOptions() so they bypass the
summaryStatsQc() gate. In a real analysis you’d construct
the raw QtlSumStats / GwasSumStats from your
data, then run summaryStatsQc() to populate that slot (see
the summary-statistics QC vignette).
Individual-level fine-mapping
fineMappingPipeline() always needs to know which
genotype window to fit on. With an individual-level
QtlDataset, pass cisWindow (in basepairs
around each trait’s TSS), supply an explicit region, or set
both per-call.
fmr <- fineMappingPipeline(
qtlDatasetExample,
methods = "susie",
cisWindow = 1e6
)
fmr## QtlFineMappingResult: 1 entries
## 1 studies, 1 contexts, 1 traits, 1 methods
## LD reference: NULL (individual-level fit)Per-tuple PIPs:
pipAll <- pip(
fmr,
studyName = "study1",
context = "brain",
trait = "ENSG_example",
method = "susie"
)
head(sort(pipAll, decreasing = TRUE), 5)## chr22:15611133:C:G chr22:15570730:C:T chr22:15569194:A:G chr22:15567276:C:T
## 0.2111065 0.1427403 0.1361272 0.1270685
## chr22:15594669:T:C
## 0.1169922Credible sets:
credibleSets(
fmr,
studyName = "study1",
context = "brain",
trait = "ENSG_example",
method = "susie",
coverage = 0.95
)## GRanges object with 0 ranges and 12 metadata columns:
## seqnames ranges strand | variant_id chrom pos A1
## <Rle> <IRanges> <Rle> | <character> <character> <integer> <character>
## A2 N af posterior_effect posterior_se theta
## <character> <numeric> <numeric> <numeric> <numeric> <numeric>
## pip logBF
## <numeric> <numeric>
## -------
## seqinfo: no sequencesThe top_loci table consolidates per-method PIPs, CS
membership, and sumstats into one long-format data.frame:
head(
topLoci(
fmr,
studyName = "study1",
context = "brain",
trait = "ENSG_example",
method = "susie"
),
5
)## GRanges object with 5 ranges and 24 metadata columns:
## seqnames ranges strand | variant_id chrom pos
## <Rle> <IRanges> <Rle> | <character> <character> <integer>
## [1] chr22 15567276 * | chr22:15567276:C:T chr22 15567276
## [2] chr22 15569194 * | chr22:15569194:A:G chr22 15569194
## [3] chr22 15570730 * | chr22:15570730:C:T chr22 15570730
## [4] chr22 15583103 * | chr22:15583103:G:A chr22 15583103
## [5] chr22 15594252 * | chr22:15594252:T:C chr22 15594252
## A1 A2 N af posterior_effect posterior_se
## <character> <character> <numeric> <numeric> <numeric> <numeric>
## [1] T C 165 0.199387 0.03363552 0.0918419
## [2] G A 165 0.200617 0.03621764 0.0950491
## [3] T C 165 0.200000 0.03810962 0.0972983
## [4] A G 165 0.263975 0.00630325 0.0395730
## [5] C T 165 0.173780 0.02151834 0.0741082
## theta pip logBF cs_95 cs_70 cs_50
## <numeric> <numeric> <numeric> <character> <character> <character>
## [1] 0 0.1270685 5.41846 susie_0 susie_0 susie_0
## [2] 0 0.1361272 5.48732 susie_0 susie_0 susie_0
## [3] 0 0.1427403 5.53476 susie_0 susie_0 susie_0
## [4] 0 0.0271347 3.87454 susie_0 susie_0 susie_0
## [5] 0 0.0839251 5.00365 susie_0 susie_0 susie_0
## cs_95_purity cs_70_purity cs_50_purity within_cs_pip method
## <numeric> <numeric> <numeric> <numeric> <character>
## [1] 0 0 0 NA susie
## [2] 0 0 0 NA susie
## [3] 0 0 0 NA susie
## [4] 0 0 0 NA susie
## [5] 0 0 0 NA susie
## gene event grange_start grange_end
## <character> <character> <integer> <integer>
## [1] <NA> <NA> 14680768 16681267
## [2] <NA> <NA> 14680768 16681267
## [3] <NA> <NA> 14680768 16681267
## [4] <NA> <NA> 14680768 16681267
## [5] <NA> <NA> 14680768 16681267
## -------
## seqinfo: 1 sequence from an unspecified genome; no seqlengthsMulti-study fine-mapping
A MultiStudyQtlDataset runs the pipeline per-study and
concatenates the results into a single
QtlFineMappingResult:
msFmr <- fineMappingPipeline(
multiStudyQtlDatasetExample,
methods = "susie",
cisWindow = 1e6
)
table(msFmr$study)##
## study1 study2
## 1 1QTL summary-statistics fine-mapping
fmrSs <- fineMappingPipeline(qtlSumStatsExample, methods = "susie")
fmrSs## QtlFineMappingResult: 1 entries
## 1 studies, 1 contexts, 1 traits, 1 methods
## LD reference: plink1 @ pecotmr://extdata/toy_canonicalEach row of the QtlSumStats collection names its LD
reference (here a genotype panel); the pipeline pulls each region’s LD
from its row’s reference automatically.
GWAS summary-statistics fine-mapping
fmrGwas <- fineMappingPipeline(gwasSumStatsS4Example, methods = "susie")
fmrGwas## GwasFineMappingResult: 1 entries
## 1 studies, 1 methods
## LD reference: plink1 @ pecotmr://extdata/toy_canonicalThe GWAS path is similar to the QTL-sumstats path but only annotates by study.
Common parameters
fineMappingPipeline(
qtlDatasetExample,
methods = c("susie", "susieInf"),
contexts = "brain",
traitId = "ENSG_example",
cisWindow = 5e5,
credibleSetParam = CredibleSetParam(
coverage = 0.95,
secondaryCoverage = c(0.7, 0.5),
signalCutoff = 0.025,
minAbsCorr = 0.8,
maxNumSingleEffects = 10L
),
initializeWithSusieInf = TRUE,
phenotypeFilterParam = PhenotypeFilterParam(naAction = "drop"),
verbose = 1
)-
phenotypeFilterParamis the run’s phenotype QC. ItsnaActionsays how phenotype NAs are handled:"drop"(the default) removes samples with any NA across the requested traits,"impute"mean-imputes each trait independently.outlierAction = "drop"also removes samples flagged as multivariate outliers. The QtlDataset accessors apply none of this; it belongs to the run. -
initializeWithSusieInf = TRUEruns SuSiE-inf first and uses its result to initialise the susie fit, which improves credible-set purity in regions with background polygenic signal. -
CredibleSetParam(maxNumSingleEffects = ...)bounds how many single effects the fit may carry, and so how many credible sets can exist. -
CredibleSetParam(secondaryCoverage = ...)adds secondary CS columns at the requested coverages (alongside the primarycoverage); useful when you want both narrow and broad CS reports. - Each of these settings is a constructor rather than a loose argument, so a misspelled option is an error instead of being silently ignored.
Joint multi-context / multi-trait fits
The joint methods (mvsusie, fsusie) fit
across traits when the input has two or more, else across contexts. To
choose the axes yourself, pass a jointSpecification: a
character vector of axes, or list(axes = , scope = ) to
limit the fit to listed studies, contexts or traits:
# `qd`: a QtlDataset whose traits are measured in brain, liver and more.
fineMappingPipeline(
qd,
methods = "mvsusie",
jointSpecification = list(list(
axes = "context",
scope = list(context = c("brain", "liver"))
))
)The jointSpecification argument in
?fineMappingPipeline gives the full grammar, including the
study axis for summary statistics. On summary statistics a
joint fit uses one LD matrix, so the rows it combines must share an LD
reference; a group whose rows do not is an error.
Reading the result
QtlFineMappingResult is a GRangesList
subclass, one element per (study, context, trait, method)
tuple. The element is the fitted variant set as a GRanges,
so the result is addressable by genomic range as well as by tuple; the
tuple itself and the per-tuple payloads live in mcols(),
reachable with $:
fmr## QtlFineMappingResult: 1 entries
## 1 studies, 1 contexts, 1 traits, 1 methods
## LD reference: NULL (individual-level fit)
colnames(fmr)## [1] "study" "context" "trait" "method"
## [5] "susieFit" "cvResult" "derivedColumns" "traitPos"
## [9] "region" "ldReference" "cvPartition" "dataSource"
## [13] "standardized"
lengths(fmr) # variants per tuple## study1|brain|ENSG_example|susie
## 198Rather than reaching into the payload, use the view accessors.
topLoci() returns a GRanges over the whole
collection with the tuple columns carried in mcols(); the
others return tidy tables:
## GRanges object with 3 ranges and 31 metadata columns:
## seqnames ranges strand | study context trait
## <Rle> <IRanges> <Rle> | <character> <character> <character>
## [1] chr22 15567276 * | study1 brain ENSG_example
## [2] chr22 15569194 * | study1 brain ENSG_example
## [3] chr22 15570730 * | study1 brain ENSG_example
## blockId method dataType dataSource standardized
## <character> <character> <character> <character> <character>
## [1] <NA> susie <NA> individual TRUE
## [2] <NA> susie <NA> individual TRUE
## [3] <NA> susie <NA> individual TRUE
## variant_id chrom pos A1 A2
## <character> <character> <integer> <character> <character>
## [1] chr22:15567276:C:T chr22 15567276 T C
## [2] chr22:15569194:A:G chr22 15569194 G A
## [3] chr22:15570730:C:T chr22 15570730 T C
## N af posterior_effect posterior_se theta pip
## <numeric> <numeric> <numeric> <numeric> <numeric> <numeric>
## [1] 165 0.199387 0.0336355 0.0918419 0 0.127069
## [2] 165 0.200617 0.0362176 0.0950491 0 0.136127
## [3] 165 0.200000 0.0381096 0.0972983 0 0.142740
## logBF cs_95 cs_70 cs_50 cs_95_purity cs_70_purity
## <numeric> <character> <character> <character> <numeric> <numeric>
## [1] 5.41846 susie_0 susie_0 susie_0 0 0
## [2] 5.48732 susie_0 susie_0 susie_0 0 0
## [3] 5.53476 susie_0 susie_0 susie_0 0 0
## cs_50_purity within_cs_pip gene event grange_start
## <numeric> <numeric> <character> <character> <integer>
## [1] 0 NA <NA> <NA> 14680768
## [2] 0 NA <NA> <NA> 14680768
## [3] 0 NA <NA> <NA> 14680768
## grange_end
## <integer>
## [1] 16681267
## [2] 16681267
## [3] 16681267
## -------
## seqinfo: 1 sequence from an unspecified genome; no seqlengths
head(credibleSets(fmr), 3) # credible-set membership## GRanges object with 0 ranges and 8 metadata columns:
## seqnames ranges strand | study context trait blockId
## <Rle> <IRanges> <Rle> | <character> <character> <character> <character>
## method dataType dataSource standardized
## <character> <character> <character> <character>
## -------
## seqinfo: no sequences
head(credibleSetSummary(fmr), 3) # one row per credible set## DataFrame with 0 rows and 8 columnsPassing a tuple selects a single fit and returns that tuple’s values on their own — a bare vector rather than a table:
key <- list(
studyName = fmr$study[[1L]],
context = fmr$context[[1L]],
trait = fmr$trait[[1L]],
method = fmr$method[[1L]]
)
pipOne <- do.call(pip, c(list(fmr), key))
length(pipOne)## [1] 198The underlying SuSiE fit is available with susieFit()
and the cross-validation record with cvResult(), both keyed
the same way. Range-based subsetting uses the standard verbs: subset the
collection by its rows’ analysis windows, or narrow a ranged view
directly:
window <- GenomicRanges::GRanges(
"chr22",
IRanges::IRanges(20000000, 21000000)
)
# the fits whose analysis window overlaps it
nrow(fmr[IRanges::overlapsAny(genomicRegion(fmr), window)])## [1] 0
# the top loci inside it
length(IRanges::subsetByOverlaps(topLoci(fmr), window))## [1] 0Calling a fitter directly
Each methods = value dispatches to an exported fitting
function, usable on its own when you have matrices rather than a
QtlDataset:
methods = |
Individual-level | Summary-statistics |
|---|---|---|
susie |
susieWeights() |
susieRssWeights() |
susieInf |
susieInfWeights() |
susieInfRssWeights() |
susieAsh |
susieAshWeights() |
susieAshRssWeights() |
mvsusie |
mvsusieWeights() / fitMvsusie()
|
mvsusieRssWeights() / fitMvsusieRss()
|
fsusie |
fsusieWeights() / fitFsusie()
|
— |
The *Weights() functions return a weight vector, which
is what makes a fine-mapper usable as a TWAS
weight method; the fit*() functions return the fit
object itself when you want the posterior rather than the
coefficients.
Two are worth calling out:
-
fsusieWrapper()wrapsfitFsusie()with the credible-band and affected-region post-processing, andfsusieGetCs(),fsusieCredibleBand()andfsusieAffectedRegions()read those out. fSuSiE needs a functional phenotype — a curve per sample, not a scalar — so it does not accept the sameYthe other methods do. -
fitSusieInfThenSusieRss()chains the two:susie_inffirst to estimate the infinitesimal component, then SuSiE-RSS on the residual signal. That is whatmethods = "susieInf"runs on summary statistics.
Next steps
-
QtlFineMappingResultis used bycolocPipelinefor QTL-GWAS colocalization and bytwasWeightsPipelineso that fine-mapping weights are contributed to the TWAS ensemble weights model alongside regularized regression methods. - Run the full
causalInferencePipeline(TWAS + MR) with aQtlFineMappingResultandGwasSumStats.
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.4
##
## 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 dbplyr_2.6.0
## [19] backports_1.5.1 XVector_0.50.0
## [21] Rsamtools_2.26.0 rmarkdown_2.32
## [23] tzdb_0.5.0 UCSC.utils_1.6.1
## [25] ragg_1.5.2 purrr_1.2.2
## [27] bit_4.6.0 xfun_0.61
## [29] Rfast_2.1.5.2 MultiAssayExperiment_1.36.2
## [31] cachem_1.1.0 GenomeInfoDb_1.46.2
## [33] jsonlite_2.0.0 blob_1.3.0
## [35] DelayedArray_0.36.1 reshape_0.8.10
## [37] tictoc_1.2.1 BiocParallel_1.44.0
## [39] irlba_2.4.1 parallel_4.5.3
## [41] R6_2.6.1 bslib_0.12.0
## [43] stringi_1.8.9 RColorBrewer_1.1-3
## [45] GenomicRanges_1.62.1 jquerylib_0.1.4
## [47] numDeriv_2016.8-1.1 Rcpp_1.1.2
## [49] Seqinfo_1.0.0 SummarizedExperiment_1.40.0
## [51] knitr_1.52 readr_2.2.0
## [53] IRanges_2.44.0 BiocBaseUtils_1.12.0
## [55] Matrix_1.7-6 tidyselect_1.2.1
## [57] abind_1.4-8 yaml_2.3.12
## [59] codetools_0.2-20 curl_8.0.0
## [61] metafor_5.0-1 lattice_0.23-1
## [63] tibble_3.3.1 plyr_1.8.9
## [65] withr_3.0.3 Biobase_2.70.0
## [67] S7_0.2.2 evaluate_1.0.5
## [69] archive_1.1.13 desc_1.4.3
## [71] RcppParallel_5.1.11-2 BiocFileCache_3.0.0
## [73] Biostrings_2.78.0 filelock_1.0.3
## [75] pillar_1.11.1 MatrixGenerics_1.22.0
## [77] metadat_1.6-0 checkmate_2.3.3
## [79] stats4_4.5.3 generics_0.1.4
## [81] vroom_1.7.1 mathjaxr_2.0-0
## [83] S4Vectors_0.48.1 hms_1.1.4
## [85] ggplot2_4.0.3 scales_1.4.0
## [87] glue_1.8.1 tools_4.5.3
## [89] robustbase_0.99-7 fs_2.1.0
## [91] grid_4.5.3 PlinkMatrix_1.0.0
## [93] tidyr_1.3.2 nlme_3.1-170
## [95] cli_3.6.6 zigg_0.0.2
## [97] textshaping_1.0.5 mixsqp_0.3-54
## [99] S4Arrays_1.10.1 dplyr_1.2.1
## [101] gtable_0.3.6 DEoptimR_1.2-1
## [103] sass_0.4.10 digest_0.6.39
## [105] BiocGenerics_0.56.0 SparseArray_1.10.10
## [107] htmlwidgets_1.6.4 farver_2.1.2
## [109] memoise_2.0.1 htmltools_0.5.9
## [111] pkgdown_2.2.1 lifecycle_1.0.5
## [113] httr_1.4.9 bit64_4.8.6
## [115] MASS_7.3-66