Causal inference with pecotmr: TWAS Z-scores and MR
pecotmr authors
2026-11-02
Source:vignettes/twas-zscore.Rmd
twas-zscore.RmdOverview
causalInferencePipeline() pairs QTL fine-mapping and/or
regularized regression weights one or more GWAS studies to produce, per
(qtlStudy, context, trait, method, gwasStudy) tuple:
- a TWAS Z-score (and p-value), and
- when
fineMappingResultis supplied, a Mendelian Randomization Wald-ratio estimate over the QTL credible sets.
Three valid input combinations:
| What you supply | TWAS-Z source | MR? |
|---|---|---|
twasWeights only |
per-tuple weights from TwasWeights
|
no |
fineMappingResult only |
SuSiE coefficients from topLoci()
|
yes |
| both | weights from TwasWeights
|
yes (MR uses the FMR) |
Bundled inputs
data(qtlDatasetExample, gwasSumStatsS4Example)
# Pre-compute the QTL-side products from the bundled QtlDataset.
# crossValidation carries the per-fold susie fits to the weights pipeline,
# which never fine-maps for itself.
fmr <- fineMappingPipeline(
qtlDatasetExample,
methods = "susie",
cisWindow = 1e6,
crossValidationParam = CrossValidationParam(folds = 3)
)
tw <- twasWeightsPipeline(
qtlDatasetExample,
methods = "susie",
cisWindow = 1e6,
fineMappingResult = fmr
)TWAS Z only (weights alone)
res <- causalInferencePipeline(
gwasSumStats = gwasSumStatsS4Example,
twasWeights = tw
)
head(S4Vectors::mcols(res))## DataFrame with 1 row and 14 columns
## qtlStudy context trait method gwasStudy twasZ
## <character> <character> <character> <character> <character> <numeric>
## 1 study1 brain ENSG_example susie trait1 13.8179
## twasPval waldRatio waldRatioSe mrPval nIV Q I2
## <numeric> <numeric> <numeric> <numeric> <integer> <numeric> <numeric>
## 1 1.98845e-43 NA NA NA NA NA NA
## nCs
## <integer>
## 1 NAThe result is a GRanges object positioned at the variant
span of each QTL weight set, with identity columns
(qtlStudy, context, trait,
method, gwasStudy) and the TWAS columns
(twasZ, twasPval). MR columns are
NA in this mode.
TWAS Z + MR (weights + fine-mapping)
When you also pass fineMappingResult, the pipeline adds
Wald-ratio MR on top of the TWAS computation. The MR uses each
fine-mapping entry’s topLoci rows as candidate instrumental
variables.
res2 <- causalInferencePipeline(
gwasSumStats = gwasSumStatsS4Example,
twasWeights = tw,
fineMappingResult = fmr,
mrParam = MrParam(method = "ivwPerVariant", pipCutoff = 0.5)
)
S4Vectors::mcols(res2)[, c(
"twasZ",
"twasPval",
"waldRatio",
"waldRatioSe",
"mrPval",
"nIV"
)]## DataFrame with 1 row and 6 columns
## twasZ twasPval waldRatio waldRatioSe mrPval nIV
## <numeric> <numeric> <numeric> <numeric> <numeric> <integer>
## 1 13.8179 1.98845e-43 NA NA NA 0MR variants: per-variant vs CS-aware IVW
MrParam(method=) |
What it does | Extra filters |
|---|---|---|
"ivwPerVariant" (default) |
Treats every topLoci variant with
pip > pipCutoff as an IV; pools their Wald ratios via
inverse-variance weighting. |
pipCutoff |
"csAware" |
Groups variants by credible set (column cs in
topLoci), forms a PIP-weighted composite Wald ratio per CS
(subject to per-CS cumulative-PIP cutoff cpipCutoff), then
IVW-pools across CSs. Reports Cochran’s Q and I-squared. |
cpipCutoff |
causalInferencePipeline(
gwasSumStats = gwasSumStatsS4Example,
twasWeights = tw,
fineMappingResult = fmr,
mrParam = MrParam(method = "csAware", cpipCutoff = 0.5)
)Cross-method p-value combination
When a (qtlStudy, context, trait, gwasStudy) group has
TWAS results from multiple weight methods, combineMethods
runs combinePValues() to produce a single combined p-value
per group.
# Multi-method weights so combination is meaningful
tw_multi <- twasWeightsPipeline(
qtlDatasetExample,
methods = c("susie", "lasso")
)
causalInferencePipeline(
gwasSumStats = gwasSumStatsS4Example,
twasWeights = tw_multi,
combineMethods = c("acat", "hmp")
)Accepted values for combineMethods: acat,
hmp, bonferroni, fisher,
stouffer, invchisq, gbj,
bj, hc, ghc, minp,
gbj_omni, aspu, gates (see
?combinePValues).
The per-tuple p-value itself comes from waldTestPval(),
which is exported so the same test can be applied to any effect estimate
and standard error:
waldTestPval(beta = 2.5, se = 0.8, n = 10000)## [1] 0.00178314n sets the degrees of freedom, so a small study gives a
more conservative p-value than the normal approximation would.
LD reference consistency
Both twasWeights and fineMappingResult
record an LD reference per row. Where a row has one, it must be the same
LD reference as the gwasSumStats rows it overlaps carry.
The pipeline errors on mismatch.
A QTL row with no LD reference (e.g. weights learned from individual-level data without any RSS step) skips the panel check — useful for analyses that mix individual-level QTL fits with sumstats GWAS.
Common parameters
causalInferencePipeline(
gwasSumStats = gwasSumStatsS4Example,
twasWeights = tw,
fineMappingResult = fmr,
mrParam = MrParam(
method = "ivwPerVariant",
pipCutoff = 0.5 # cpipCutoff instead, for method = "csAware"
),
combineMethods = NULL
) # e.g. c("acat", "hmp")Reading the result
The GRanges output indexes results by the QTL variant
span (the union of variants used in each tuple’s weights). Extract
per-tuple values with S4Vectors::mcols():
mc <- as.data.frame(S4Vectors::mcols(res2))
# Significant TWAS hits at Bonferroni alpha = 0.05 across rows
mc[mc$twasPval < 0.05 / nrow(mc), ]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.4
##
## loaded via a namespace (and not attached):
## [1] DBI_1.3.0 bitops_1.1-0
## [3] httr2_1.3.0 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 RSQLite_3.53.3
## [11] systemfonts_1.3.2 vctrs_0.7.3
## [13] quadprog_1.5-8 stringr_1.6.0
## [15] pkgconfig_2.0.3 crayon_1.5.3
## [17] fastmap_1.2.0 dbplyr_2.6.0
## [19] backports_1.5.1 XVector_0.50.0
## [21] Rsamtools_2.26.0 rmarkdown_2.32
## [23] tzdb_0.5.0 UCSC.utils_1.6.1
## [25] ragg_1.5.2 purrr_1.2.2
## [27] bit_4.6.0 xfun_0.61
## [29] Rfast_2.1.5.2 MultiAssayExperiment_1.36.2
## [31] cachem_1.1.0 GenomeInfoDb_1.46.2
## [33] jsonlite_2.0.0 blob_1.3.0
## [35] DelayedArray_0.36.1 reshape_0.8.10
## [37] tictoc_1.2.1 BiocParallel_1.44.0
## [39] irlba_2.4.1 parallel_4.5.3
## [41] R6_2.6.1 bslib_0.12.0
## [43] stringi_1.8.9 RColorBrewer_1.1-3
## [45] GenomicRanges_1.62.1 jquerylib_0.1.4
## [47] numDeriv_2016.8-1.1 Rcpp_1.1.2
## [49] Seqinfo_1.0.0 SummarizedExperiment_1.40.0
## [51] knitr_1.52 readr_2.2.0
## [53] IRanges_2.44.0 BiocBaseUtils_1.12.0
## [55] Matrix_1.7-6 tidyselect_1.2.1
## [57] abind_1.4-8 yaml_2.3.12
## [59] codetools_0.2-20 curl_8.0.0
## [61] metafor_5.0-1 lattice_0.23-1
## [63] tibble_3.3.1 plyr_1.8.9
## [65] withr_3.0.3 Biobase_2.70.0
## [67] S7_0.2.2 evaluate_1.0.5
## [69] archive_1.1.13 desc_1.4.3
## [71] RcppParallel_5.1.11-2 BiocFileCache_3.0.0
## [73] Biostrings_2.78.0 filelock_1.0.3
## [75] pillar_1.11.1 MatrixGenerics_1.22.0
## [77] metadat_1.6-0 checkmate_2.3.3
## [79] stats4_4.5.3 generics_0.1.4
## [81] vroom_1.7.1 mathjaxr_2.0-0
## [83] S4Vectors_0.48.1 hms_1.1.4
## [85] ggplot2_4.0.3 scales_1.4.0
## [87] glue_1.8.1 tools_4.5.3
## [89] robustbase_0.99-7 fs_2.1.0
## [91] grid_4.5.3 PlinkMatrix_1.0.0
## [93] tidyr_1.3.2 nlme_3.1-170
## [95] cli_3.6.6 zigg_0.0.2
## [97] textshaping_1.0.5 mixsqp_0.3-54
## [99] S4Arrays_1.10.1 dplyr_1.2.1
## [101] gtable_0.3.6 DEoptimR_1.2-1
## [103] sass_0.4.10 digest_0.6.39
## [105] BiocGenerics_0.56.0 SparseArray_1.10.10
## [107] htmlwidgets_1.6.4 farver_2.1.2
## [109] memoise_2.0.1 htmltools_0.5.9
## [111] pkgdown_2.2.1 lifecycle_1.0.5
## [113] httr_1.4.9 bit64_4.8.6
## [115] MASS_7.3-66