Pairwise colocalization with colocPipeline
pecotmr authors
2026-09-04
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 (with finemappingMethods) over the exact window
used in the QtlFineMappingResult, then runs
coloc.bf_bf per (QTL, GWAS) pair.
res <- colocPipeline(
qtlFineMappingResult = qtlFmr,
gwasInput = gwasSumStatsS4Example,
finemappingMethods = "susie"
)
res## ColocResult with 10 colocalized pair(s)
## QTL studies : study1
## GWAS studies: 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 QTL FineMappingResult has an LD sketch because it
was run on summary statistics it must match the LD sketch for the
GWASSumStats.
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 |
|---|---|---|
filterLbfCs |
TRUE keeps only effects that produced a credible set
(trimmedFit$sets$cs_index). |
FALSE |
filterLbfCsSecondary |
Numeric coverage in (0,1). Runs a CS-concentration filter at that coverage. | NULL |
filterLbfCsConcentration |
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 |
priorTol |
Drop effects whose estimated prior variance is at or below this threshold. | 1e-9 |
colocPipeline(
qtlFineMappingResult = qtlFmr,
gwasInput = gwasSumStatsS4Example,
filterLbfCsSecondary = 0.5, # 50% CS concentration filter
filterLbfCsConcentration = 0.5
) # span < 25% of locus variantsEnrichment-aware coloc (enloc mode)
Pass a per-(gwasStudy, qtlContext) enrichment table (the
output of qtlEnrichmentPipeline()) to scale the
shared-signal prior:
p_{12}^{\text{used}} = \min\bigl(p_{12}\,(1 + \text{enrichment}),\, p_{12,\max}\bigr)
Pairs without a matching enrichment row fall back to the baseline
p12 with a warning. Two extra columns are added to the
output: enrichment (the lookup) and p12Used
(what was actually passed to coloc.bf_bf).
enrichmentTable <- data.frame(
gwasStudy = "trait1",
qtlContext = "brain",
enrichment = 2.5,
stringsAsFactors = FALSE
)
colocPipeline(
qtlFineMappingResult = qtlFmr,
gwasInput = gwasSumStatsS4Example,
enrichment = enrichmentTable,
p12Max = 1e-3
)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,
filterLbfCs = FALSE,
filterLbfCsSecondary = NULL, # e.g. 0.5
filterLbfCsConcentration = 0.5,
priorTol = 1e-9,
p1 = 1e-4,
p2 = 1e-4,
p12 = 5e-6,
finemappingMethods = "susie", # inline GWAS FM when needed
returnGwasFineMapping = FALSE,
enrichment = NULL, # data.frame for enloc mode
p12Max = 1e-3,
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 <- getColocPairs(res)
head(pairs[order(pairs$PP.H4.abf, decreasing = TRUE), ])## # A tibble: 6 × 21
## study context trait method gwasStudy gwasMethod blockId qtlCs gwasCs nSnps
## <chr> <chr> <chr> <chr> <chr> <chr> <chr> <int> <int> <int>
## 1 study1 brain ENSG_ex… susie trait1 susie chr22_… 1 2 198
## 2 study1 brain ENSG_ex… susie trait1 susie chr22_… 1 4 198
## 3 study1 brain ENSG_ex… susie trait1 susie chr22_… 1 6 198
## 4 study1 brain ENSG_ex… susie trait1 susie chr22_… 1 8 198
## 5 study1 brain ENSG_ex… susie trait1 susie chr22_… 1 10 198
## 6 study1 brain ENSG_ex… susie trait1 susie chr22_… 1 1 198
## # ℹ 11 more variables: 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 <- getColocVariants(res)
head(variants[order(variants$colocPp, decreasing = TRUE), ])## # A tibble: 6 × 24
## study context trait method gwasStudy gwasMethod blockId qtlCs gwasCs nSnps
## <chr> <chr> <chr> <chr> <chr> <chr> <chr> <int> <int> <int>
## 1 study1 brain ENSG_ex… susie trait1 susie chr22_… 1 2 198
## 2 study1 brain ENSG_ex… susie trait1 susie chr22_… 1 4 198
## 3 study1 brain ENSG_ex… susie trait1 susie chr22_… 1 6 198
## 4 study1 brain ENSG_ex… susie trait1 susie chr22_… 1 8 198
## 5 study1 brain ENSG_ex… susie trait1 susie chr22_… 1 10 198
## 6 study1 brain ENSG_ex… susie trait1 susie chr22_… 1 1 198
## # ℹ 14 more variables: 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(getColocCredibleSets(res, coverage = 0.95, minPp4 = 0.1))## # A tibble: 6 × 27
## study context trait method gwasStudy gwasMethod blockId qtlCs gwasCs nSnps
## <chr> <chr> <chr> <chr> <chr> <chr> <chr> <int> <int> <int>
## 1 study1 brain ENSG_ex… susie trait1 susie chr22_… 1 1 198
## 2 study1 brain ENSG_ex… susie trait1 susie chr22_… 1 2 198
## 3 study1 brain ENSG_ex… susie trait1 susie chr22_… 1 3 198
## 4 study1 brain ENSG_ex… susie trait1 susie chr22_… 1 4 198
## 5 study1 brain ENSG_ex… susie trait1 susie chr22_… 1 5 198
## 6 study1 brain ENSG_ex… susie trait1 susie chr22_… 1 6 198
## # ℹ 17 more variables: 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).
getColocGenes(res)## # A tibble: 1 × 9
## study context trait method gwasStudy gwasMethod PP.H4 nQtlCs nPairs
## <chr> <chr> <chr> <chr> <chr> <chr> <dbl> <int> <int>
## 1 study1 brain ENSG_example susie trait1 susie 1 1 10The 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.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.8
##
## loaded via a namespace (and not attached):
## [1] coloc_5.2.3 bitops_1.0-9
## [3] gridExtra_2.3.1 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 systemfonts_1.3.2
## [11] vctrs_0.7.3 quadprog_1.5-8
## [13] stringr_1.6.0 pkgconfig_2.0.3
## [15] crayon_1.5.3 fastmap_1.2.0
## [17] XVector_0.50.0 utf8_1.2.6
## [19] Rsamtools_2.26.0 rmarkdown_2.32
## [21] tzdb_0.5.0 ragg_1.5.2
## [23] purrr_1.2.2 bit_4.6.0
## [25] xfun_0.60 Rfast_2.1.5.2
## [27] MultiAssayExperiment_1.36.1 cachem_1.1.0
## [29] jsonlite_2.0.0 DelayedArray_0.36.0
## [31] reshape_0.8.10 tictoc_1.2.1
## [33] BiocParallel_1.44.0 irlba_2.3.7
## [35] parallel_4.5.3 R6_2.6.1
## [37] bslib_0.12.0 stringi_1.8.9
## [39] RColorBrewer_1.1-3 GenomicRanges_1.62.1
## [41] jquerylib_0.1.4 numDeriv_2016.8-1.1
## [43] Rcpp_1.1.2 Seqinfo_1.0.0
## [45] SummarizedExperiment_1.40.0 knitr_1.51
## [47] readr_2.2.0 IRanges_2.44.0
## [49] BiocBaseUtils_1.12.0 Matrix_1.7-6
## [51] splines_4.5.3 tidyselect_1.2.1
## [53] abind_1.4-8 yaml_2.3.12
## [55] viridis_0.6.5 codetools_0.2-20
## [57] metafor_5.0-1 lattice_0.23-1
## [59] tibble_3.3.1 plyr_1.8.9
## [61] Biobase_2.70.0 withr_3.0.3
## [63] S7_0.2.2 evaluate_1.0.5
## [65] archive_1.1.13 desc_1.4.3
## [67] survival_3.8-11 RcppParallel_5.1.11-2
## [69] snpStats_1.60.0 Biostrings_2.78.0
## [71] pillar_1.11.1 MatrixGenerics_1.22.0
## [73] metadat_1.6-0 stats4_4.5.3
## [75] generics_0.1.4 vroom_1.7.1
## [77] mathjaxr_2.0-0 S4Vectors_0.48.0
## [79] hms_1.1.4 ggplot2_4.0.3
## [81] scales_1.4.0 glue_1.8.1
## [83] tools_4.5.3 data.table_1.18.6.1
## [85] fs_2.1.0 grid_4.5.3
## [87] tidyr_1.3.2 nlme_3.1-170
## [89] cli_3.6.6 zigg_0.0.2
## [91] textshaping_1.0.5 mixsqp_0.3-54
## [93] S4Arrays_1.10.1 viridisLite_0.4.3
## [95] dplyr_1.2.1 gtable_0.3.6
## [97] sass_0.4.10 digest_0.6.39
## [99] BiocGenerics_0.56.0 SparseArray_1.10.8
## [101] htmlwidgets_1.6.4 farver_2.1.2
## [103] htmltools_0.5.9 pkgdown_2.2.1
## [105] lifecycle_1.0.5 bit64_4.8.6
## [107] MASS_7.3-66