Molecular QTL Enrichment Analysis for Complex Traits

This 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: 1) now enrichment is based on all QTL variants whether or not they are inside signal clusters; 2) Causal QTL are sampled from SuSiE single effects, not signal clusters; 3) Allow a variant to be QTL for not only multiple conditions (eg cell types) but also multiple regions (eg 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

Load the pecotmr package for this enrichment analysis.

library(pecotmr)

Load fine-mapping results from the package example data and write them to temporary RDS files, since xqtlEnrichmentWrapper expects file paths.

data(gwas_finemapping_example)
data(qtl_finemapping_example)
gwas_path <- tempfile(fileext = ".rds")
qtl_path <- tempfile(fileext = ".rds")
saveRDS(gwas_finemapping_example, gwas_path)
saveRDS(qtl_finemapping_example, qtl_path)
gwas <- gwas_finemapping_example
qtl <- qtl_finemapping_example

Input data requirements

GWAS data requirements

In the GWAS data, the file is stored as a list, with the fine mapping results already stored under the first layer.

names(gwas[[1]])

## [1] "alpha" "pip"   "V"     "sets"

Posterior inclusion probabilities of GWAS data:

plot(gwas[[1]]$pip, ylab = "PIP", xlab = "Variant index", main = "GWAS PIPs", pch = 16, cex = 0.6)

QTL data requirements

In the QTL data, the file is stored as a list and is more complicated. The first layer is named as the molecular trait region ID, followed by different contexts as the second layer. The variant names and fine mapping results can be found under these contexts (variant_names and susie_result_trimmed).

names(qtl)

## [1] "region_1"

names(qtl[["region_1"]])

## [1] "context_1"

names(qtl[["region_1"]][["context_1"]])

## [1] "susie_result_trimmed" "variant_names"

names(qtl[["region_1"]][["context_1"]][["susie_result_trimmed"]])

## [1] "alpha" "pip"   "V"     "sets"

Posterior inclusion probabilities of QTL data:

plot(qtl[["region_1"]][["context_1"]][["susie_result_trimmed"]]$pip,
     ylab = "PIP", xlab = "Variant index", main = "QTL PIPs", pch = 16, cex = 0.6)

Perform enrichment analysis

The input for xqtlEnrichmentWrapper requires the paths for GWAS and QTL fine mapping data. You need to specify the list name of SuSiE results in each dataset as xqtl_finemapping_obj and the list name of variant names as xqtl_varname_obj. For example, if SuSiE results for QTL data are under qtl[["region_1"]][["context_1"]][["susie_result_trimmed"]], then xqtl_finemapping_obj would be c("context_1", "susie_result_trimmed"), and xqtl_varname_obj would be c("context_1", "variant_names") since the variant names are stored directly under context_1. For GWAS data, the fine mapping results and variant names are stored under the first layer directly, so you do not need to specify gwas_finemapping_obj and gwas_varname_obj. However, if your data does not follow this structure, you may need to specify these parameters accordingly.

enrich_res <- xqtlEnrichmentWrapper(
  gwas_files = gwas_path, xqtl_files = qtl_path,
  xqtl_finemapping_obj = c("context_1", "susie_result_trimmed"),
  xqtl_varname_obj = c("context_1", "variant_names")
)

## Warning in xqtlEnrichmentWrapper(gwas_files = gwas_path, xqtl_files = qtl_path,
## : xqtlEnrichmentWrapper() has been removed. Use qtlEnrichmentPipeline() with
## FineMappingResult collections for the GWAS and the QTLs.

print(enrich_res[[1]])

## NULL

print(head(enrich_res[["unused_xqtl_variants"]][[1]]))

## NULL

The enrichment analysis output is stored as a 2-layer list. The first layer contains the enrichment results, which will be used as priors in coloc analysis. The second layer stores the QTL variants that were not detected in the GWAS dataset.

Session information

Here are some details about the computing environment, including the versions of R, and the R packages, used to generate these results.

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.33.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.6.7
## 
## 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] yaml_2.3.12                 knitr_1.51                 
## [19] S4Arrays_1.10.1             htmlwidgets_1.6.4          
## [21] bit_4.6.0                   DelayedArray_0.36.0        
## [23] plyr_1.8.9                  RColorBrewer_1.1-3         
## [25] abind_1.4-8                 BiocParallel_1.44.0        
## [27] purrr_1.2.2                 numDeriv_2016.8-1.1        
## [29] BiocGenerics_0.56.0         desc_1.4.3                 
## [31] grid_4.5.3                  stats4_4.5.3               
## [33] susieR_0.16.4               ggplot2_4.0.3              
## [35] scales_1.4.0                MASS_7.3-66                
## [37] SummarizedExperiment_1.40.0 cli_3.6.6                  
## [39] rmarkdown_2.31              metafor_5.0-1              
## [41] crayon_1.5.3                ragg_1.5.2                 
## [43] generics_0.1.4              otel_0.2.0                 
## [45] RcppParallel_5.1.11-2       tzdb_0.5.0                 
## [47] cachem_1.1.0                stringr_1.6.0              
## [49] metadat_1.6-0               parallel_4.5.3             
## [51] XVector_0.50.0              matrixStats_1.5.0          
## [53] vctrs_0.7.3                 Matrix_1.7-5               
## [55] jsonlite_2.0.0              IRanges_2.44.0             
## [57] hms_1.1.4                   S4Vectors_0.48.0           
## [59] bit64_4.8.2                 mixsqp_0.3-54              
## [61] irlba_2.3.7                 systemfonts_1.3.2          
## [63] tidyr_1.3.2                 jquerylib_0.1.4            
## [65] glue_1.8.1                  pkgdown_2.2.1              
## [67] codetools_0.2-20            stringi_1.8.7              
## [69] gtable_0.3.6                GenomicRanges_1.62.1       
## [71] quadprog_1.5-8              tibble_3.3.1               
## [73] pillar_1.11.1               htmltools_0.5.9            
## [75] Seqinfo_1.0.0               R6_2.6.1                   
## [77] zigg_0.0.2                  textshaping_1.0.5          
## [79] vroom_1.7.1                 evaluate_1.0.5             
## [81] lattice_0.22-9              Biobase_2.70.0             
## [83] readr_2.2.0                 Rsamtools_2.26.0           
## [85] tictoc_1.2.1                Rfast_2.1.5.2              
## [87] bslib_0.11.0                Rcpp_1.1.2                 
## [89] SparseArray_1.10.8          nlme_3.1-170               
## [91] xfun_0.60                   fs_2.1.0                   
## [93] MatrixGenerics_1.22.0       pkgconfig_2.0.3