Molecular QTL Enrichment Analysis for Complex Traits
Ru Feng
2026-09-04
Source:vignettes/xqtl-enrichment.Rmd
xqtl-enrichment.RmdThis vignette demonstrates the enrichment analysis part of the
pecotmr package, which largely follows from
fastENLOC (https://github.com/xqwen/fastenloc) but uses
susieR fitted objects as input to estimate priors for use
with the coloc package (coloc v5, aka SuSiE-coloc). The
main differences are:
- enrichment is based on all QTL variants whether or not they are inside signal clusters;
- causal QTL are sampled from SuSiE single effects, not signal clusters;
- a variant may be a QTL for not only multiple conditions (e.g. cell types) but also multiple regions (e.g. genes).
The example uses de-identified data shipped with the package. All sample names, variant positions, and identifiers are synthetic.
Set up environment and load data
The enrichment pipeline consumes fine-mapping results as
pecotmr S4 objects: a GwasFineMappingResult
for the GWAS and a QtlFineMappingResult for the molecular
QTLs. Both are shipped as example datasets.
## GwasFineMappingResult: 1 entries
## 1 studies, 1 methods
## LD sketch: plink1 @ pecotmr://extdata/toy_canonical
qtlFineMappingExample## QtlFineMappingResult: 1 entries
## 1 studies, 1 contexts, 1 traits, 1 methods
## LD sketch: NULL (individual-level fit)Input data
Each object is a GRangesList-backed collection, one
element per tuple; use the accessors rather than reaching into slots.
getTopLoci() returns the per-variant posterior table
(identity columns, pip, and credible-set membership):
head(getTopLoci(gwasFineMappingExample))## # A tibble: 6 × 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 study_1 NA NA region… susie chr22:323… chr22 3.23e7 T C NA
## 2 study_1 NA NA region… susie chr22:323… chr22 3.24e7 G A NA
## 3 study_1 NA NA region… susie chr22:323… chr22 3.24e7 G A NA
## 4 study_1 NA NA region… susie chr22:323… chr22 3.24e7 C T NA
## 5 study_1 NA NA region… susie chr22:323… chr22 3.24e7 G A NA
## 6 study_1 NA NA region… susie chr22:323… chr22 3.24e7 A C NA
## # ℹ 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(getTopLoci(qtlFineMappingExample))## # A tibble: 6 × 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 study_1 contex… gene… NA susie chr22:325… chr22 3.26e7 C T 415
## 2 study_1 contex… gene… NA susie chr22:327… chr22 3.27e7 C T 415
## 3 study_1 contex… gene… NA susie chr22:327… chr22 3.27e7 A C 415
## 4 study_1 contex… gene… NA susie chr22:327… chr22 3.27e7 A G 415
## 5 study_1 contex… gene… NA susie chr22:327… chr22 3.27e7 C T 415
## 6 study_1 contex… gene… NA susie chr22:327… chr22 3.27e7 A G 415
## # ℹ 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>Posterior inclusion probabilities across all variants in each region, taken from the underlying SuSiE fit:
gwas_pip <- getSusieFit(gwasFineMappingExample)$pip
qtl_pip <- getSusieFit(qtlFineMappingExample)$pip
op <- par(mfrow = c(1, 2))
plot(
gwas_pip,
ylab = "PIP",
xlab = "Variant index",
main = "GWAS PIPs",
pch = 16,
cex = 0.6
)
plot(
qtl_pip,
ylab = "PIP",
xlab = "Variant index",
main = "QTL PIPs",
pch = 16,
cex = 0.6
)
par(op)Perform enrichment analysis
qtlEnrichmentPipeline() takes the two fine-mapping
collections directly and estimates, for every
(gwasStudy, qtlStudy, qtlContext) combination, how enriched
the molecular QTLs are among the GWAS causal variants. The GWAS
collection must be RSS-derived (it carries an ldSketch);
gwasFineMappingExample already does.
enrich_res <- qtlEnrichmentPipeline(
gwasFineMappingResult = gwasFineMappingExample,
qtlFineMappingResult = qtlFineMappingExample,
verbose = FALSE
)## Warning: numGwas is not provided. Estimating piGwas from the data. Note that
## this estimate may be biased if the input gwasPip does not contain genome-wide
## variants.## Warning: piQtl is not provided. Estimating piQtl from the data. Note that this
## estimate may be biased if either 1) the input susieQtlRegions does not have
## enough data, or 2) the single effects only include variables inside of credible
## sets or signal clusters.## Fine-mapped GWAS and QTL data loaded successfully for enrichment analysis!
## Proportion of xQTL missing from GWAS variants: 0 in MI round 0
## Proportion of xQTL missing from GWAS variants: 0 in MI round 1
## Proportion of xQTL missing from GWAS variants: 0 in MI round 2
## Proportion of xQTL missing from GWAS variants: 0 in MI round 3
## Proportion of xQTL missing from GWAS variants: 0 in MI round 4
## Proportion of xQTL missing from GWAS variants: 0 in MI round 5
## Proportion of xQTL missing from GWAS variants: 0 in MI round 6
## Proportion of xQTL missing from GWAS variants: 0 in MI round 7
## Proportion of xQTL missing from GWAS variants: 0 in MI round 8
## Proportion of xQTL missing from GWAS variants: 0 in MI round 9
## Proportion of xQTL missing from GWAS variants: 0 in MI round 10
## Proportion of xQTL missing from GWAS variants: 0 in MI round 11
## Proportion of xQTL missing from GWAS variants: 0 in MI round 12
## Proportion of xQTL missing from GWAS variants: 0 in MI round 13
## Proportion of xQTL missing from GWAS variants: 0 in MI round 14
## Proportion of xQTL missing from GWAS variants: 0 in MI round 15
## Proportion of xQTL missing from GWAS variants: 0 in MI round 16
## Proportion of xQTL missing from GWAS variants: 0 in MI round 17
## Proportion of xQTL missing from GWAS variants: 0 in MI round 18
## Proportion of xQTL missing from GWAS variants: 0 in MI round 19
## Proportion of xQTL missing from GWAS variants: 0 in MI round 20
## Proportion of xQTL missing from GWAS variants: 0 in MI round 21
## Proportion of xQTL missing from GWAS variants: 0 in MI round 22
## Proportion of xQTL missing from GWAS variants: 0 in MI round 23
## Proportion of xQTL missing from GWAS variants: 0 in MI round 24
## EM updates completed!
## Outlier filtering removed 1 MI round(s)
enrich_res## # A tibble: 1 × 6
## gwasStudy qtlStudy qtlContext enrichment enrichmentSe enrichmentLogOdds
## <chr> <chr> <chr> <dbl> <dbl> <dbl>
## 1 study_1 study_1 context_1 NA NA NAThe result is a data frame with one row per
(gwasStudy, qtlStudy, qtlContext) combination, giving the
enrichment estimate (enrichment), its standard error
(enrichmentSe), and the log-odds
(enrichmentLogOdds). These enrichment priors can then be
supplied to a colocalization analysis. This synthetic single-region
example is for illustration only; stable enrichment estimates require
many fine-mapped regions.
Session information
Here are some details about the computing environment, including the versions of R and the R packages used to generate these results.
## 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] 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 magrittr_2.0.5
## [13] compiler_4.5.3 rlang_1.3.0
## [15] sass_0.4.10 tools_4.5.3
## [17] utf8_1.2.6 yaml_2.3.12
## [19] knitr_1.51 S4Arrays_1.10.1
## [21] htmlwidgets_1.6.4 bit_4.6.0
## [23] DelayedArray_0.36.0 plyr_1.8.9
## [25] RColorBrewer_1.1-3 abind_1.4-8
## [27] BiocParallel_1.44.0 withr_3.0.3
## [29] purrr_1.2.2 numDeriv_2016.8-1.1
## [31] BiocGenerics_0.56.0 desc_1.4.3
## [33] grid_4.5.3 stats4_4.5.3
## [35] susieR_0.16.6 ggplot2_4.0.3
## [37] scales_1.4.0 MASS_7.3-66
## [39] MultiAssayExperiment_1.36.1 SummarizedExperiment_1.40.0
## [41] cli_3.6.6 rmarkdown_2.32
## [43] metafor_5.0-1 crayon_1.5.3
## [45] ragg_1.5.2 generics_0.1.4
## [47] otel_0.2.0 RcppParallel_5.1.11-2
## [49] tzdb_0.5.0 BiocBaseUtils_1.12.0
## [51] cachem_1.1.0 stringr_1.6.0
## [53] metadat_1.6-0 parallel_4.5.3
## [55] XVector_0.50.0 matrixStats_1.5.0
## [57] vctrs_0.7.3 Matrix_1.7-6
## [59] jsonlite_2.0.0 IRanges_2.44.0
## [61] hms_1.1.4 S4Vectors_0.48.0
## [63] bit64_4.8.6 mixsqp_0.3-54
## [65] irlba_2.3.7 archive_1.1.13
## [67] systemfonts_1.3.2 tidyr_1.3.2
## [69] jquerylib_0.1.4 glue_1.8.1
## [71] pkgdown_2.2.1 codetools_0.2-20
## [73] stringi_1.8.9 gtable_0.3.6
## [75] GenomicRanges_1.62.1 quadprog_1.5-8
## [77] tibble_3.3.1 pillar_1.11.1
## [79] htmltools_0.5.9 Seqinfo_1.0.0
## [81] R6_2.6.1 zigg_0.0.2
## [83] textshaping_1.0.5 vroom_1.7.1
## [85] evaluate_1.0.5 lattice_0.23-1
## [87] Biobase_2.70.0 readr_2.2.0
## [89] tictoc_1.2.1 Rsamtools_2.26.0
## [91] Rfast_2.1.5.2 bslib_0.12.0
## [93] Rcpp_1.1.2 SparseArray_1.10.8
## [95] nlme_3.1-170 xfun_0.60
## [97] fs_2.1.0 MatrixGenerics_1.22.0
## [99] pkgconfig_2.0.3