Causal inference with pecotmr: TWAS Z-scores and MR

Overview

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 fineMappingResult is 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 each `FineMappingEntry`’s `topLoci`	yes
both	weights from `TwasWeights`	yes (MR uses the FMR)

library(pecotmr)

Bundled inputs

data(qtl_dataset_example, gwas_sumstats_s4_example)

# Pre-compute the QTL-side products from the bundled QtlDataset.
fmr <- fineMappingPipeline(qtl_dataset_example, methods = "susie",
                            cisWindow = 1e6)
tw  <- twasWeightsPipeline(qtl_dataset_example,
                            methods           = "susie",
                            cisWindow         = 1e6,
                            fineMappingResult = fmr)

TWAS Z only (weights alone)

res <- causalInferencePipeline(
  gwasSumStats = gwas_sumstats_s4_example,
  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        NA

The 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      = gwas_sumstats_s4_example,
  twasWeights       = tw,
  fineMappingResult = fmr,
  mrMethod          = "ivwPerVariant",
  mrPipCutoff       = 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         0

MR variants: per-variant vs CS-aware IVW

`mrMethod`	What it does	Extra filters
`"ivwPerVariant"` (default)	Treats every `topLoci` variant with `pip > mrPipCutoff` as an IV; pools their Wald ratios via inverse-variance weighting.	`mrPipCutoff`
`"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 `mrCpipCutoff`), then IVW-pools across CSs. Reports Cochran’s Q and I-squared.	`mrCpipCutoff`

causalInferencePipeline(
  gwasSumStats      = gwas_sumstats_s4_example,
  twasWeights       = tw,
  fineMappingResult = fmr,
  mrMethod          = "csAware",
  mrCpipCutoff      = 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(qtl_dataset_example,
                                 methods = c("susie", "lasso"))
causalInferencePipeline(
  gwasSumStats   = gwas_sumstats_s4_example,
  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).

LD-sketch consistency

Both twasWeights and fineMappingResult carry an ldSketch slot. When provided, they must reference the same LD panel as gwasSumStats. The pipeline errors on mismatch.

A QTL input with ldSketch = NULL (e.g. weights learned from individual-level data without any RSS step) skips the panel check for that input — useful for analyses that mix individual-level QTL fits with sumstats GWAS.

Common parameters

causalInferencePipeline(
  gwasSumStats      = gwas_sumstats_s4_example,
  twasWeights       = tw,
  fineMappingResult = fmr,
  mrPipCutoff       = 0.5,            # ivwPerVariant
  mrMethod          = "ivwPerVariant",
  mrCpipCutoff      = 0.5,            # 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), ]

pecotmr authors

2026-07-30