Pairwise colocalization with colocPipeline
pecotmr authors
2026-10-22
Source:vignettes/coloc-pipeline.Rmd
coloc-pipeline.RmdOverview
colocPipeline() runs pairwise QTL ↔︎ GWAS colocalization
using coloc::coloc.bf_bf on the SuSiE log Bayes-factor
(LBF) matrices carried by each fine-mapping entry. For multi-trait
colocalization across many QTL outcomes in one model, use colocboostPipeline
instead.
Inputs:
-
QTL input — a
QtlFineMappingResultproduced byfineMappingPipeline. -
GWAS input — either a
GwasSumStats(fine-mapping run internally) or aGwasFineMappingResult(used directly).
The output is a ColocResult: one element per tested
(QTL credible set, GWAS credible set, block) pair, holding
that pair’s aligned variants. Pair-level posteriors
(PP.H0.abf … PP.H4.abf) and the per-variant
SNP.PP.H4 are both kept, and four accessors project the
object to whichever granularity an analysis needs — see Reading the result.
Bundled inputs
data(qtlDatasetExample, gwasSumStatsS4Example)
# Build a QtlFineMappingResult from the bundled QtlDataset.
qtlFmr <- fineMappingPipeline(
qtlDatasetExample,
methods = "susie",
cisWindow = 1e6
)QTL FMR + GWAS sumstats (inline GWAS fine-mapping)
The recommend approach is to provide colocPipeline() a
QtlFineMappingResult and a GwasSumStats. The
pipeline runs fineMappingPipeline() internally on the GWAS
sumstats (configured by GwasFineMappingParam()) over the
exact window used in the QtlFineMappingResult, then runs
coloc.bf_bf per (QTL, GWAS) pair.
res <- colocPipeline(
qtlFineMappingResult = qtlFmr,
gwasInput = gwasSumStatsS4Example,
gwasFineMappingParam = GwasFineMappingParam(methods = "susie")
)
res## ColocResult with 10 colocalized pair(s)
## studies : study1
## paired with : trait1
## variants : 1980 across all pairs
## max PP.H4 : 0.9166Use returnGwasFineMapping = TRUE to get the inline GWAS
FMR back as an attribute on the result:
res2 <- colocPipeline(
qtlFineMappingResult = qtlFmr,
gwasInput = gwasSumStatsS4Example,
returnGwasFineMapping = TRUE
)
attr(res2, "gwasFineMapping")If the QtlFineMappingResult has an LD reference because
it was run on summary statistics, it must be the same LD reference as
the GwasSumStats carries.
QTL FMR + pre-computed GWAS FMR
If you already have a GwasFineMappingResult (typically
computed over LD blocks), you can use that fine-mapping result
directly:
gwasFmr <- fineMappingPipeline(gwasSumStatsS4Example, methods = "susie")
colocPipeline(qtlFineMappingResult = qtlFmr, gwasInput = gwasFmr)The disadvantage of this approach is that LD blocks will not align with cis-QTL windows. Ultimately this means that the variable selection in the QTL and GWAS fine-mapping will have been done over different variant sets, which may reduce the sensitivity for colocalizations, particuarly in areas with complex LD.
Filtering effects fed into coloc
coloc::coloc.bf_bf accepts one LBF matrix for the QTL
and one for hte GWAS. By default, every effect with non-trivial prior
variance (V > priorTol) gets included. Two opt-in
filters can adjust this:
| Filter | Description | Default |
|---|---|---|
ColocLbfFilterParam(filterLbfCs=) |
TRUE keeps only effects that produced a credible set
(trimmedFit$sets$cs_index). |
FALSE |
ColocLbfFilterParam(secondary=) |
Numeric coverage in (0,1). Runs a CS-concentration filter at that coverage. | NULL |
ColocLbfFilterParam(concentration=) |
Numeric in (0,1). Concentration factor: a CS at coverage X is kept
only if it spans fewer than nVariants * X * concentration
variants. |
0.5 |
ColocLbfFilterParam(priorTol=) |
Drop effects whose estimated prior variance is at or below this threshold. | 1e-9 |
colocPipeline(
qtlFineMappingResult = qtlFmr,
gwasInput = gwasSumStatsS4Example,
colocLbfFilterParam = ColocLbfFilterParam(
secondary = 0.5, # 50% CS concentration filter
concentration = 0.5
)
) # span < 25% of locus variantsEnrichment-aware colocalisation (enloc mode)
Passing a per-(gwasStudy, qtlContext) enrichment table —
the output of qtlEnrichmentPipeline() — switches the pair
from coloc to enloc. This is a port of fastenloc, not a
reweighting of coloc: the priors are derived from the
enrichment estimate rather than the baseline p12 being
scaled by it.
The estimator returns the enrichment parameters (a_0, a_1) on the log-odds scale, and the
priors follow as fastenloc’s set_enrich_params
writes them, with P_{\text{eQTL}} the
per-variant QTL prior:
p_1 = (1 - P_{\text{eQTL}})\,\sigma(a_0), \quad p_2 = P_{\text{eQTL}}\,\bigl(1 - \sigma(a_0 + a_1)\bigr), \quad p_{12} = P_{\text{eQTL}}\,\sigma(a_0 + a_1)
qtlEnrichmentPipeline() carries all three out as the
colocP1, colocP2 and colocP12
columns, so the table you pass already holds them. Because the
enrichment determines the priors, supplying priors as well
is an error rather than a silent precedence rule.
What makes the scoring enloc rather than coloc is that the prior depends on whether the GWAS hit is the causal QTL — \sigma(a_0 + a_1) when it is, \sigma(a_0) when it is not. coloc applies one p_{12} to every shared configuration and cannot express that.
Pairs without a matching enrichment row fall back to the baseline
priors with a warning. Four columns are added: enrichment
and p12Used (the per-pair factor and the derived
shared-signal prior), plus RCP and LCP,
enloc’s regional and locus colocalisation probabilities. Those last two
are the same quantities as PP.H4.abf and
PP.H3.abf + PP.H4.abf, reported under enloc’s names as
well.
enrichmentTable <- qtlEnrichmentPipeline(
gwasFineMappingResult = gwasFineMappingExample,
qtlFineMappingResult = qtlFineMappingExample,
verbose = FALSE
)
colocPipeline(
qtlFineMappingResult = qtlFmr,
gwasInput = gwasSumStatsS4Example,
enrichment = enrichmentTable
)PIP renormalization on variant intersection
When the QTL and GWAS fine-mapping results were fit on different
variant sets (common when the user declines RAISS imputation, or when
the GWAS FMR carries variants the QTL FMR doesn’t),
adjustPips = TRUE (default) renormalizes the PIPs for each
objectl to the cross-FMR variant intersection before any per-pair
inference. Pass FALSE to use the FMRs as supplied.
colocPipeline(
qtlFineMappingResult = qtlFmr,
gwasInput = gwasSumStatsS4Example,
adjustPips = TRUE
) # defaultCommon parameters
colocPipeline(
qtlFineMappingResult = qtlFmr,
gwasInput = gwasSumStatsS4Example,
colocLbfFilterParam = ColocLbfFilterParam(
filterLbfCs = FALSE,
secondary = NULL, # e.g. 0.5
concentration = 0.5,
priorTol = 1e-9
),
priors = ColocPriorParam(p1 = 1e-4, p2 = 1e-4, p12 = 5e-6),
gwasFineMappingParam = GwasFineMappingParam( # inline GWAS FM when needed
methods = "susie",
credibleSetParam = CredibleSetParam(coverage = 0.95)
),
returnGwasFineMapping = FALSE,
enrichment = NULL, # data.frame for enloc mode
adjustPips = TRUE
)Reading the result
Nothing is filtered when the result is written, so every narrowing below is a view decision the caller makes — and can revisit — rather than something baked in at run time.
Pairs
The native coloc.bf_bf granularity: one row per tested
pair. Every testable pair is reported verbatim, including both halves of
a credible set that straddles a block boundary.
pairs <- colocPairs(res)
head(pairs[order(pairs$PP.H4.abf, decreasing = TRUE), ])## # A tibble: 6 × 23
## study context trait method gwasStudy gwasContext gwasTrait gwasMethod blockId
## <chr> <chr> <chr> <chr> <chr> <chr> <chr> <chr> <chr>
## 1 study1 brain ENSG… susie trait1 NA NA susie chr22_…
## 2 study1 brain ENSG… susie trait1 NA NA susie chr22_…
## 3 study1 brain ENSG… susie trait1 NA NA susie chr22_…
## 4 study1 brain ENSG… susie trait1 NA NA susie chr22_…
## 5 study1 brain ENSG… susie trait1 NA NA susie chr22_…
## 6 study1 brain ENSG… susie trait1 NA NA susie chr22_…
## # ℹ 14 more variables: qtlCs <int>, gwasCs <int>, nSnps <int>, hit1 <chr>,
## # hit2 <chr>, PP.H0.abf <dbl>, PP.H1.abf <dbl>, PP.H2.abf <dbl>,
## # PP.H3.abf <dbl>, PP.H4.abf <dbl>, idx1 <int>, idx2 <int>,
## # qtlRetainedMass <dbl>, gwasRetainedMass <dbl>qtlCs / gwasCs are the fine-mapping effect
indices, which are stable across blocks; coloc’s own idx1 /
idx2 number the rows of a single call and are not
comparable between them. qtlRetainedMass /
gwasRetainedMass report how much of each effect’s posterior
survived reconciliation to the shared variant set — a low value means
the effect rested largely on variants the other side does not carry.
Variants
The per-variant layer: colocPp is
PP.H4.abf * SNP.PP.H4, the posterior that this particular
variant is the shared causal one.
variants <- colocVariants(res)
head(variants[order(variants$colocPp, decreasing = TRUE), ])## # A tibble: 6 × 26
## study context trait method gwasStudy gwasContext gwasTrait gwasMethod blockId
## <chr> <chr> <chr> <chr> <chr> <chr> <chr> <chr> <chr>
## 1 study1 brain ENSG… susie trait1 NA NA susie chr22_…
## 2 study1 brain ENSG… susie trait1 NA NA susie chr22_…
## 3 study1 brain ENSG… susie trait1 NA NA susie chr22_…
## 4 study1 brain ENSG… susie trait1 NA NA susie chr22_…
## 5 study1 brain ENSG… susie trait1 NA NA susie chr22_…
## 6 study1 brain ENSG… susie trait1 NA NA susie chr22_…
## # ℹ 17 more variables: qtlCs <int>, gwasCs <int>, nSnps <int>, hit1 <chr>,
## # hit2 <chr>, PP.H0.abf <dbl>, PP.H1.abf <dbl>, PP.H2.abf <dbl>,
## # PP.H3.abf <dbl>, PP.H4.abf <dbl>, idx1 <int>, idx2 <int>,
## # qtlRetainedMass <dbl>, gwasRetainedMass <dbl>, variant_id <chr>,
## # SNP.PP.H4 <dbl>, colocPp <dbl>Credible sets
The coloc credible set is the smallest set of variants whose
cumulative SNP.PP.H4 reaches coverage. It is a
third variant set — not guaranteed to be a subset of either
input credible set — so its purity is recomputed from LD rather than
inherited. Filtering happens before any LD extraction, so a stricter
threshold costs strictly less.
head(colocCredibleSets(res, coverage = 0.95, minPp4 = 0.1))## # A tibble: 6 × 29
## study context trait method gwasStudy gwasContext gwasTrait gwasMethod blockId
## <chr> <chr> <chr> <chr> <chr> <chr> <chr> <chr> <chr>
## 1 study1 brain ENSG… susie trait1 NA NA susie chr22_…
## 2 study1 brain ENSG… susie trait1 NA NA susie chr22_…
## 3 study1 brain ENSG… susie trait1 NA NA susie chr22_…
## 4 study1 brain ENSG… susie trait1 NA NA susie chr22_…
## 5 study1 brain ENSG… susie trait1 NA NA susie chr22_…
## 6 study1 brain ENSG… susie trait1 NA NA susie chr22_…
## # ℹ 20 more variables: qtlCs <int>, gwasCs <int>, nSnps <int>, hit1 <chr>,
## # hit2 <chr>, PP.H0.abf <dbl>, PP.H1.abf <dbl>, PP.H2.abf <dbl>,
## # PP.H3.abf <dbl>, PP.H4.abf <dbl>, idx1 <int>, idx2 <int>,
## # qtlRetainedMass <dbl>, gwasRetainedMass <dbl>, csSize <int>,
## # csCoverage <dbl>, leadVariant <chr>, leadPp <dbl>, csVariants <list>,
## # purity <dbl>Genes
Pooled to one row per gene. Within a QTL credible set the per-block
and per-GWAS-credible-set results are mutually exclusive — the shared
causal variant is in one block or the other — so they are summed;
distinct QTL credible sets are independent signals and combine as
1 - prod(1 - p).
colocGenes(res)## # A tibble: 1 × 11
## study context trait method gwasStudy gwasContext gwasTrait gwasMethod PP.H4
## <chr> <chr> <chr> <chr> <chr> <chr> <chr> <chr> <dbl>
## 1 study1 brain ENSG_e… susie trait1 NA NA susie 1
## # ℹ 2 more variables: nQtlCs <int>, nPairs <int>The warning above is the diagnostic that rule produces: a within-credible-set sum above 1 means several GWAS credible sets are competing for one QTL signal. The value is clipped to 1, and the pair view shows which pairs were involved.
If you need the flat table the pipeline used to return,
as.data.frame(res) gives the pair view.
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] coloc_5.2.3 DBI_1.3.0
## [3] bitops_1.1-0 gridExtra_2.3.1
## [5] httr2_1.3.0 rlang_1.3.0
## [7] magrittr_2.0.5 otel_0.2.0
## [9] matrixStats_1.5.0 susieR_0.16.6
## [11] compiler_4.5.3 RSQLite_3.53.3
## [13] systemfonts_1.3.2 vctrs_0.7.3
## [15] quadprog_1.5-8 stringr_1.6.0
## [17] pkgconfig_2.0.3 crayon_1.5.3
## [19] fastmap_1.2.0 dbplyr_2.6.0
## [21] backports_1.5.1 XVector_0.50.0
## [23] utf8_1.2.6 Rsamtools_2.26.0
## [25] rmarkdown_2.32 tzdb_0.5.0
## [27] ragg_1.5.2 purrr_1.2.2
## [29] bit_4.6.0 xfun_0.61
## [31] Rfast_2.1.5.2 MultiAssayExperiment_1.36.2
## [33] cachem_1.1.0 jsonlite_2.0.0
## [35] blob_1.3.0 DelayedArray_0.36.1
## [37] reshape_0.8.10 tictoc_1.2.1
## [39] BiocParallel_1.44.0 irlba_2.4.1
## [41] parallel_4.5.3 R6_2.6.1
## [43] bslib_0.12.0 stringi_1.8.9
## [45] RColorBrewer_1.1-3 GenomicRanges_1.62.1
## [47] jquerylib_0.1.4 numDeriv_2016.8-1.1
## [49] Rcpp_1.1.2 Seqinfo_1.0.0
## [51] SummarizedExperiment_1.40.0 knitr_1.52
## [53] readr_2.2.0 IRanges_2.44.0
## [55] BiocBaseUtils_1.12.0 Matrix_1.7-6
## [57] tidyselect_1.2.1 viridis_0.6.5
## [59] abind_1.4-8 yaml_2.3.12
## [61] codetools_0.2-20 curl_8.0.0
## [63] metafor_5.0-1 lattice_0.23-1
## [65] tibble_3.3.1 plyr_1.8.9
## [67] withr_3.0.3 Biobase_2.70.0
## [69] S7_0.2.2 evaluate_1.0.5
## [71] archive_1.1.13 desc_1.4.3
## [73] RcppParallel_5.1.11-2 BiocFileCache_3.0.0
## [75] Biostrings_2.78.0 filelock_1.0.3
## [77] pillar_1.11.1 MatrixGenerics_1.22.0
## [79] metadat_1.6-0 checkmate_2.3.3
## [81] stats4_4.5.3 generics_0.1.4
## [83] vroom_1.7.1 mathjaxr_2.0-0
## [85] S4Vectors_0.48.1 hms_1.1.4
## [87] ggplot2_4.0.3 scales_1.4.0
## [89] glue_1.8.1 tools_4.5.3
## [91] data.table_1.18.6.1 robustbase_0.99-7
## [93] fs_2.1.0 grid_4.5.3
## [95] PlinkMatrix_1.0.0 tidyr_1.3.2
## [97] nlme_3.1-170 cli_3.6.6
## [99] zigg_0.0.2 textshaping_1.0.5
## [101] viridisLite_0.4.3 mixsqp_0.3-54
## [103] S4Arrays_1.10.1 dplyr_1.2.1
## [105] gtable_0.3.6 DEoptimR_1.2-1
## [107] sass_0.4.10 digest_0.6.39
## [109] BiocGenerics_0.56.0 SparseArray_1.10.10
## [111] htmlwidgets_1.6.4 farver_2.1.2
## [113] memoise_2.0.1 htmltools_0.5.9
## [115] pkgdown_2.2.1 lifecycle_1.0.5
## [117] bit64_4.8.6 MASS_7.3-66