Fine-mapping with pecotmr

Overview

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, susieSer
`MultiStudyQtlDataset`	Multiple `QtlDataset` objects (optionally with an embedded `QtlSumStats`)	susie, susieInf, susieAsh, mvsusie, fsusie (individual-level only), susieSer (individual-level only)
`QtlSumStats`	QTL summary statistics annotated by `(study, context, trait)`	susie, susieInf, susieAsh, mvsusie
`GwasSumStats`	GWAS summary statistics annotated by `study`	susie, susieInf, susieAsh, mvsusie

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 always a FineMappingResult (a GRangesList-subclass collection, one element per fine-mapped tuple).

library(pecotmr)

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: 165
##   Scale residuals: TRUE

The *_sumstats_example objects ship with a non-empty qcInfo slot 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 sketch: NULL (individual-level fit)

Per-tuple PIPs:

pip <- getPip(
    fmr,
    study = "study1",
    context = "brain",
    trait = "ENSG_example",
    method = "susie"
)
head(sort(pip, 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.1169922

Credible sets:

getCs(
    fmr,
    study = "study1",
    context = "brain",
    trait = "ENSG_example",
    method = "susie",
    coverage = 0.95
)

## # A tibble: 0 × 11
## # ℹ 11 variables: variant_id <chr>, chrom <chr>, pos <int>, A1 <chr>, A2 <chr>,
## #   N <dbl>, af <dbl>, beta <dbl>, se <dbl>, pip <dbl>, logBF <dbl>

The top_loci table consolidates per-method PIPs, CS membership, and sumstats into one long-format data.frame:

head(
    getTopLoci(
        fmr,
        study = "study1",
        context = "brain",
        trait = "ENSG_example",
        method = "susie"
    ),
    5
)

## # A tibble: 5 × 23
##   variant_id    chrom    pos A1    A2        N    af    beta     se    pip logBF
##   <chr>         <chr>  <int> <chr> <chr> <dbl> <dbl>   <dbl>  <dbl>  <dbl> <dbl>
## 1 chr22:155672… chr22 1.56e7 T     C       165 0.199 0.0336  0.0918 0.127   5.42
## 2 chr22:155691… chr22 1.56e7 G     A       165 0.201 0.0362  0.0950 0.136   5.49
## 3 chr22:155707… chr22 1.56e7 T     C       165 0.2   0.0381  0.0973 0.143   5.53
## 4 chr22:155831… chr22 1.56e7 A     G       165 0.264 0.00630 0.0396 0.0271  3.87
## 5 chr22:155942… chr22 1.56e7 C     T       165 0.174 0.0215  0.0741 0.0839  5.00
## # ℹ 12 more variables: cs_95 <chr>, cs_70 <chr>, cs_50 <chr>,
## #   cs_95_purity <dbl>, cs_70_purity <dbl>, cs_50_purity <dbl>,
## #   within_cs_pip <dbl>, method <chr>, gene <chr>, event <chr>,
## #   grange_start <int>, grange_end <int>

Multi-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      1

QTL summary-statistics fine-mapping

fmrSs <- fineMappingPipeline(qtlSumStatsExample, methods = "susie")
fmrSs

## QtlFineMappingResult: 1 entries
##   1 studies, 1 contexts, 1 traits, 1 methods
##   LD sketch: plink1 @ pecotmr://extdata/toy_canonical

The QtlSumStats collection carries an ldSketch genotype panel; the pipeline pulls the per-region LD from there automatically.

GWAS summary-statistics fine-mapping

fmrGwas <- fineMappingPipeline(gwasSumStatsS4Example, methods = "susie")
fmrGwas

## GwasFineMappingResult: 1 entries
##   1 studies, 1 methods
##   LD sketch: plink1 @ pecotmr://extdata/toy_canonical

The 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",
    region = "chr22:25000000-26000000",
    cisWindow = 5e5,
    coverage = 0.95,
    secondaryCoverage = c(0.7, 0.5),
    signalCutoff = 0.025,
    minAbsCorr = 0.8,
    addSusieInf = TRUE,
    naAction = "drop", # "drop" (default) or "impute"
    verbose = 1
)

naAction controls how phenotype NAs are handled by getResidualizedPhenotypes(): "drop" removes samples with any NA across the requested traits, "impute" mean-imputes each trait independently.
addSusieInf = TRUE runs SuSiE-inf first and uses its result to initialise the susie fit, which improves credible-set purity in regions with background polygenic signal.
secondaryCoverage adds secondary CS columns at the requested coverages (alongside the primary coverage); useful when you want both narrow and broad CS reports.

Joint multi-context / multi-trait fits

jointSpecification triggers the dispatch for cross-axis joint methods (mvsusie for cross-context, mvsusie or fsusie for multi-trait):

fineMappingPipeline(
    multiStudyQtlDatasetExample,
    methods = "mvsusie",
    jointSpecification = list(axis = "context", contexts = c("brain", "liver"))
)

See ?jointSpecification for the specification grammar.

Reading the result

FineMappingResult 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 sketch: NULL (individual-level fit)

colnames(fmr)

## [1] "study"    "context"  "trait"    "method"   "susieFit" "cvResult" "traitPos"

lengths(fmr) # variants per tuple

## [1] 198

Rather than reaching into the payload, use the view accessors. Each returns a tidy table over the whole collection, with the tuple columns carried along:

head(getTopLoci(fmr), 3) # per-variant posterior summary

## # A tibble: 3 × 27
##   study  context trait  blockId method variant_id chrom    pos A1    A2        N
##   <chr>  <chr>   <chr>  <chr>   <chr>  <chr>      <chr>  <int> <chr> <chr> <dbl>
## 1 study1 brain   ENSG_… NA      susie  chr22:155… chr22 1.56e7 T     C       165
## 2 study1 brain   ENSG_… NA      susie  chr22:155… chr22 1.56e7 G     A       165
## 3 study1 brain   ENSG_… NA      susie  chr22:155… chr22 1.56e7 T     C       165
## # ℹ 16 more variables: af <dbl>, beta <dbl>, se <dbl>, pip <dbl>, logBF <dbl>,
## #   cs_95 <chr>, cs_70 <chr>, cs_50 <chr>, cs_95_purity <dbl>,
## #   cs_70_purity <dbl>, cs_50_purity <dbl>, within_cs_pip <dbl>, gene <chr>,
## #   event <chr>, grange_start <int>, grange_end <int>

head(getCs(fmr), 3) # credible-set membership

## # A tibble: 0 × 5
## # ℹ 5 variables: study <chr>, context <chr>, trait <chr>, blockId <chr>,
## #   method <chr>

head(getCredibleSetSummary(fmr), 3) # one row per credible set

## # A tibble: 0 × 5
## # ℹ 5 variables: study <chr>, context <chr>, trait <chr>, blockId <chr>,
## #   method <chr>

Passing a tuple selects a single fit and returns that tuple’s values on their own — a bare vector rather than a table:

key <- list(
    study = fmr$study[[1L]],
    context = fmr$context[[1L]],
    trait = fmr$trait[[1L]],
    method = fmr$method[[1L]]
)
pip <- do.call(getPip, c(list(fmr), key))
length(pip)

## [1] 198

The underlying SuSiE fit is available with getSusieFit() and the cross-validation record with getCvResult(), both keyed the same way. Range-based subsetting works directly:

window <- GenomicRanges::GRanges(
    "chr22",
    IRanges::IRanges(20000000, 21000000)
)
sum(lengths(subsetRegion(fmr, window)))

## [1] 0

Next steps

QtlFineMappingResult is used by colocPipeline for QTL-GWAS colocalization and by twasWeightsPipeline so that fine-mapping weights are contributed to the TWAS ensemble weights model alongside regularized regression methods.
Run the full causalInferencePipeline (TWAS + MR) with a FineMappingResult and GwasSumStats.

Session info

sessionInfo()

## R version 4.5.3 (2026-03-11)
## Platform: x86_64-conda-linux-gnu
## Running under: Ubuntu 24.04.4 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.7.10
## 
## loaded via a namespace (and not attached):
##   [1] tidyselect_1.2.1            dplyr_1.2.1                
##   [3] farver_2.1.2                Biostrings_2.78.0          
##   [5] S7_0.2.2                    bitops_1.0-9               
##   [7] fastmap_1.2.0               reshape_0.8.10             
##   [9] mathjaxr_2.0-0              digest_0.6.39              
##  [11] lifecycle_1.0.5             survival_3.8-11            
##  [13] magrittr_2.0.5              compiler_4.5.3             
##  [15] rlang_1.3.0                 sass_0.4.10                
##  [17] tools_4.5.3                 utf8_1.2.6                 
##  [19] yaml_2.3.12                 knitr_1.51                 
##  [21] S4Arrays_1.10.1             htmlwidgets_1.6.4          
##  [23] bit_4.6.0                   DelayedArray_0.36.0        
##  [25] plyr_1.8.9                  RColorBrewer_1.1-3         
##  [27] abind_1.4-8                 BiocParallel_1.44.0        
##  [29] withr_3.0.3                 purrr_1.2.2                
##  [31] numDeriv_2016.8-1.1         BiocGenerics_0.56.0        
##  [33] desc_1.4.3                  grid_4.5.3                 
##  [35] stats4_4.5.3                susieR_0.16.6              
##  [37] ggplot2_4.0.3               scales_1.4.0               
##  [39] MASS_7.3-66                 MultiAssayExperiment_1.36.1
##  [41] SummarizedExperiment_1.40.0 cli_3.6.6                  
##  [43] rmarkdown_2.32              metafor_5.0-1              
##  [45] crayon_1.5.3                ragg_1.5.2                 
##  [47] generics_0.1.4              otel_0.2.0                 
##  [49] RcppParallel_5.1.11-2       httr_1.4.9                 
##  [51] tzdb_0.5.0                  BiocBaseUtils_1.12.0       
##  [53] cachem_1.1.0                stringr_1.6.0              
##  [55] splines_4.5.3               metadat_1.6-0              
##  [57] parallel_4.5.3              XVector_0.50.0             
##  [59] matrixStats_1.5.0           vctrs_0.7.3                
##  [61] Matrix_1.7-6                jsonlite_2.0.0             
##  [63] IRanges_2.44.0              hms_1.1.4                  
##  [65] S4Vectors_0.48.0            bit64_4.8.6                
##  [67] mixsqp_0.3-54               irlba_2.3.7                
##  [69] archive_1.1.13              systemfonts_1.3.2          
##  [71] tidyr_1.3.2                 jquerylib_0.1.4            
##  [73] snpStats_1.60.0             glue_1.8.1                 
##  [75] pkgdown_2.2.1               codetools_0.2-20           
##  [77] stringi_1.8.9               gtable_0.3.6               
##  [79] GenomeInfoDb_1.46.2         UCSC.utils_1.6.1           
##  [81] GenomicRanges_1.62.1        quadprog_1.5-8             
##  [83] tibble_3.3.1                pillar_1.11.1              
##  [85] htmltools_0.5.9             Seqinfo_1.0.0              
##  [87] R6_2.6.1                    zigg_0.0.2                 
##  [89] textshaping_1.0.5           vroom_1.7.1                
##  [91] evaluate_1.0.5              lattice_0.23-1             
##  [93] Biobase_2.70.0              readr_2.2.0                
##  [95] tictoc_1.2.1                Rsamtools_2.26.0           
##  [97] Rfast_2.1.5.2               bslib_0.12.0               
##  [99] Rcpp_1.1.2                  SparseArray_1.10.8         
## [101] nlme_3.1-170                xfun_0.60                  
## [103] fs_2.1.0                    MatrixGenerics_1.22.0      
## [105] pkgconfig_2.0.3

pecotmr authors

2026-09-08