Molecular QTL Enrichment Analysis for Complex Traits
Ru Feng
2026-10-22
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 reference: plink1 @ pecotmr://extdata/toy_canonical
qtlFineMappingExample## QtlFineMappingResult: 1 entries
## 1 studies, 1 contexts, 1 traits, 1 methods
## LD reference: 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.
topLoci() returns the per-variant posterior summary as a
GRanges, with the identity columns, pip, and
credible-set membership in mcols():
## GRanges object with 6 ranges and 31 metadata columns:
## seqnames ranges strand | study context trait
## <Rle> <IRanges> <Rle> | <character> <character> <character>
## [1] chr22 32337604 * | study_1 <NA> <NA>
## [2] chr22 32350609 * | study_1 <NA> <NA>
## [3] chr22 32354393 * | study_1 <NA> <NA>
## [4] chr22 32362470 * | study_1 <NA> <NA>
## [5] chr22 32362793 * | study_1 <NA> <NA>
## [6] chr22 32369273 * | study_1 <NA> <NA>
## blockId method dataType dataSource standardized
## <character> <character> <character> <character> <character>
## [1] region_1 susie <NA> <NA> <NA>
## [2] region_1 susie <NA> <NA> <NA>
## [3] region_1 susie <NA> <NA> <NA>
## [4] region_1 susie <NA> <NA> <NA>
## [5] region_1 susie <NA> <NA> <NA>
## [6] region_1 susie <NA> <NA> <NA>
## variant_id chrom pos A1 A2
## <character> <character> <integer> <character> <character>
## [1] chr22:32337604:C:T chr22 32337604 T C
## [2] chr22:32350609:A:G chr22 32350609 G A
## [3] chr22:32354393:A:G chr22 32354393 G A
## [4] chr22:32362470:T:C chr22 32362470 C T
## [5] chr22:32362793:A:G chr22 32362793 G A
## [6] chr22:32369273:C:A chr22 32369273 A C
## N af posterior_effect posterior_se theta pip
## <numeric> <numeric> <numeric> <numeric> <numeric> <numeric>
## [1] NA NA 0.002119882 0.00433735 0 0.1968901
## [2] NA NA 0.000792956 0.00281382 0 0.0751927
## [3] NA NA 0.001375308 0.00361922 0 0.1288942
## [4] NA NA 0.004781037 0.00552692 0 0.4367631
## [5] NA NA 0.001452388 0.00370652 0 0.1359622
## [6] NA NA 0.000270653 0.00166687 0 0.0262804
## logBF cs_95 cs_70 cs_50 cs_95_purity cs_70_purity
## <numeric> <character> <character> <character> <numeric> <numeric>
## [1] 16782.2 susie_rss_3 susie_rss_3 susie_rss_3 0.984661 0.984661
## [2] 16755.2 susie_rss_3 susie_rss_0 susie_rss_0 0.984661 0.000000
## [3] 16919.3 susie_rss_3 susie_rss_0 susie_rss_0 0.984661 0.000000
## [4] 16950.6 susie_rss_3 susie_rss_3 susie_rss_3 0.984661 0.984661
## [5] 16766.3 susie_rss_3 susie_rss_3 susie_rss_0 0.984661 0.984661
## [6] 21268.5 susie_rss_0 susie_rss_0 susie_rss_0 0.000000 0.000000
## cs_50_purity within_cs_pip gene event grange_start
## <numeric> <numeric> <character> <character> <integer>
## [1] 0.987238 0.1968901 <NA> <NA> <NA>
## [2] 0.000000 0.0751927 <NA> <NA> <NA>
## [3] 0.000000 0.1288942 <NA> <NA> <NA>
## [4] 0.987238 0.4367631 <NA> <NA> <NA>
## [5] 0.000000 0.1359622 <NA> <NA> <NA>
## [6] 0.000000 NA <NA> <NA> <NA>
## grange_end
## <integer>
## [1] <NA>
## [2] <NA>
## [3] <NA>
## [4] <NA>
## [5] <NA>
## [6] <NA>
## -------
## seqinfo: 1 sequence from an unspecified genome; no seqlengths## GRanges object with 6 ranges and 31 metadata columns:
## seqnames ranges strand | study context trait
## <Rle> <IRanges> <Rle> | <character> <character> <character>
## [1] chr22 32596246 * | study_1 context_1 gene_1
## [2] chr22 32724254 * | study_1 context_1 gene_1
## [3] chr22 32724508 * | study_1 context_1 gene_1
## [4] chr22 32728910 * | study_1 context_1 gene_1
## [5] chr22 32729588 * | study_1 context_1 gene_1
## [6] chr22 32732510 * | study_1 context_1 gene_1
## blockId method dataType dataSource standardized
## <character> <character> <character> <character> <character>
## [1] <NA> susie <NA> <NA> <NA>
## [2] <NA> susie <NA> <NA> <NA>
## [3] <NA> susie <NA> <NA> <NA>
## [4] <NA> susie <NA> <NA> <NA>
## [5] <NA> susie <NA> <NA> <NA>
## [6] <NA> susie <NA> <NA> <NA>
## variant_id chrom pos A1 A2
## <character> <character> <integer> <character> <character>
## [1] chr22:32596246:T:C chr22 32596246 C T
## [2] chr22:32724254:T:C chr22 32724254 C T
## [3] chr22:32724508:C:A chr22 32724508 A C
## [4] chr22:32728910:G:A chr22 32728910 A G
## [5] chr22:32729588:T:C chr22 32729588 C T
## [6] chr22:32732510:G:A chr22 32732510 A G
## N af posterior_effect posterior_se theta pip
## <numeric> <numeric> <numeric> <numeric> <numeric> <numeric>
## [1] 415 NA -0.01284643 0.0324407 0 0.1444970
## [2] 415 NA 0.00601338 0.0283820 0 0.0442663
## [3] 415 NA 0.00601338 0.0283820 0 0.0442663
## [4] 415 NA 0.00534990 0.0268018 0 0.0394963
## [5] 415 NA 0.00464131 0.0249437 0 0.0345102
## [6] 415 NA 0.00937127 0.0351058 0 0.0684514
## logBF cs_95 cs_70 cs_50 cs_95_purity cs_70_purity
## <numeric> <character> <character> <character> <numeric> <numeric>
## [1] 6.53015 susie_0 susie_0 susie_0 0.000000 0.000000
## [2] 15.67566 susie_1 susie_1 susie_1 0.625503 0.938229
## [3] 15.67566 susie_1 susie_1 susie_1 0.625503 0.938229
## [4] 15.56140 susie_1 susie_1 susie_0 0.625503 0.938229
## [5] 15.42615 susie_1 susie_1 susie_0 0.625503 0.938229
## [6] 16.11229 susie_1 susie_1 susie_1 0.625503 0.938229
## cs_50_purity within_cs_pip gene event grange_start
## <numeric> <numeric> <character> <character> <integer>
## [1] 0.000000 NA <NA> <NA> <NA>
## [2] 0.940595 0.0441714 <NA> <NA> <NA>
## [3] 0.940595 0.0441714 <NA> <NA> <NA>
## [4] 0.000000 0.0394021 <NA> <NA> <NA>
## [5] 0.000000 0.0344176 <NA> <NA> <NA>
## [6] 0.940595 0.0683547 <NA> <NA> <NA>
## grange_end
## <integer>
## [1] <NA>
## [2] <NA>
## [3] <NA>
## [4] <NA>
## [5] <NA>
## [6] <NA>
## -------
## seqinfo: 1 sequence from an unspecified genome; no seqlengthsPosterior inclusion probabilities across all variants in each region, taken from the underlying SuSiE fit:
gwas_pip <- susieFit(gwasFineMappingExample)$pip
qtl_pip <- susieFit(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 ldReference);
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 × 16
## gwasStudy gwasContext gwasTrait qtlStudy qtlContext enrichment enrichmentSe
## <chr> <chr> <chr> <chr> <chr> <dbl> <dbl>
## 1 study_1 NA NA study_1 context_1 -0.00000795 1.000
## # ℹ 9 more variables: enrichmentLogOdds <dbl>,
## # enrichmentLogOddsNoShrinkage <dbl>, enrichmentSeNoShrinkage <dbl>,
## # intercept <dbl>, interceptSe <dbl>, colocP1 <dbl>, colocP2 <dbl>,
## # colocP12 <dbl>, effectiveMiRounds <dbl>The 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.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] tidyselect_1.2.1 dplyr_1.2.1
## [3] farver_2.1.2 Biostrings_2.78.0
## [5] S7_0.2.2 bitops_1.1-0
## [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.52 S4Arrays_1.10.1
## [21] htmlwidgets_1.6.4 bit_4.6.0
## [23] DelayedArray_0.36.1 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.2 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.1
## [63] bit64_4.8.6 mixsqp_0.3-54
## [65] irlba_2.4.1 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 backports_1.5.1
## [91] Rsamtools_2.26.0 Rfast_2.1.5.2
## [93] bslib_0.12.0 Rcpp_1.1.2
## [95] SparseArray_1.10.10 nlme_3.1-170
## [97] checkmate_2.3.3 xfun_0.61
## [99] fs_2.1.0 MatrixGenerics_1.22.0
## [101] pkgconfig_2.0.3