Causal TWAS with ctwasPipeline
pecotmr authors
2026-10-22
Source:vignettes/ctwas-pipeline.Rmd
ctwas-pipeline.RmdOverview
A TWAS Z-score (see the TWAS Z-score + MR vignette) asks whether a gene’s predicted expression associates with the trait. It asks that question one gene at a time, so a gene can score highly purely because it sits near a strong GWAS signal that is really driven by a nearby variant.
ctwasPipeline() runs cTWAS, which
fine-maps genes and SNPs jointly within each LD block.
A gene’s posterior inclusion probability is therefore conditioned on the
SNP signal beside it, which is what separates a gene that carries the
association from one that merely neighbours it.
Inputs:
| Argument | What it is |
|---|---|
gwasSumStats |
one GwasSumStats whose elements are LD blocks |
twasWeights |
per-gene TwasWeights /
QtlFineMappingResult
|
twasZ |
optional pre-computed per-gene Z from
causalInferencePipeline()
|
fineMappingResult |
optional; enables the credible-set and PIP variant rescues |
The LD block convention
gwasSumStats is one
GwasSumStats object whose elements are LD blocks, keyed by
its blockId column — not a list of per-block objects. Build
it with
loadGwasSumStatsFromManifest(..., ldBlocks = <blocks>)
and pass it through summaryStatsQc().
A single-block input is rejected: cTWAS estimates its group priors by EM across regions, and one region gives it nothing to average over.
ldStem <- file.path(
system.file("extdata", "ld_reference", "chr22", package = "pecotmr"),
"protocol_example.LD.chr22"
)
gwasTsv <- system.file(
"extdata",
"manifests",
"protocol_example.twas.gwas_sumstats.chr22.tsv.gz",
package = "pecotmr"
)
blocks <- GenomicRanges::GRanges(
"chr22",
IRanges::IRanges(c(10000000, 15000001), c(15000000, 19000000)),
blockId = c("chr22_1", "chr22_2")
)
gwasSumStats <- loadGwasSumStatsFromManifest(
manifest = data.frame(study = "gwas1", sumStatsPath = gwasTsv),
genome = "hg38",
ldReference = ldStem,
region = "chr22:10000000-19000000",
ldBlocks = blocks
)
gwasByRegion <- summaryStatsQc(
gwasSumStats,
panelFilterParam = PanelFilterParam(mafCutoff = 0.0025)
)
data(ctwasWeightsExample)Running the pipeline
res <- ctwasPipeline(
gwasSumStats = gwasByRegion,
twasWeights = list(ctwasWeightsExample),
ctwasPriorParam = CtwasPriorParam(
thin = 1,
niterPrefit = 3,
niter = 10,
fallbackToPrefit = TRUE
),
methodArgs = CtwasOptions(
regionData = CtwasRegionDataOptions(min_group_size = 1),
estParam = CtwasEstParamOptions(
min_group_size = 1,
min_p_single_effect = 0
),
screen = CtwasScreenOptions(min_nonSNP_PIP = 0)
)
)
res## CtwasResult: 1 run(s)
## 1 GWAS stud(y/ies), 1 method(s): susieThe settings above are sized for a toy dataset and are not defaults you should carry into a real analysis:
-
thin = 1keeps every SNP (production runs thin to0.1). -
niterPrefit/niterare cut to keep the vignette fast. -
min_nonSNP_PIP = 0disables region screening. The bundled GWAS has no genome-wide-significant signal, so the default (0.5) selects no region at all and the result comes back empty. On real data, leave it alone. -
fallbackToPrefit = TRUErecovers from the accurate-EM step returningNaNon underpowered data by keeping the prefit estimates.
Where genes come from
twasWeights accepts two shapes:
- a flat
TwasWeightsorQtlFineMappingResult(or a homogeneous list of them) carryingregionprovenance. Each gene is assigned to its home LD block bystart(region), matching cTWAS’s own placement rule. - a pre-bucketed named list keyed by
region_id, used as-is.
Blocks with no gene weights still contribute their SNP-level signal to the group prior — the whole-chromosome pattern where only a few blocks carry genes.
A gene whose cis window straddles a block boundary is assigned by its
single anchor. The signal it loses across the boundary is a separate
concern, handled by mergeCtwasBoundaryRegions() below.
Trimming each gene’s variants
Weight matrices from a permissive QTL model can be large, and cTWAS’s cost scales with them. Four arguments trim each gene before cTWAS sees it:
| Argument | Effect |
|---|---|
weightCutoff |
drop variants with abs(weight) below the cutoff |
maxVariants |
cap variants per gene |
csMinCor |
keep variants in a credible set with purity at least this |
pipRescue |
keep variants with PIP above this |
csMinCor and pipRescue mark variants as
must-keep, so they survive the maxVariants cap. Both need a
fineMappingResult; without one they are no-ops and the cap
falls back to ordering by abs(weight).
ctwasPipeline(
gwasSumStats = gwasByRegion,
twasWeights = list(ctwasWeightsExample),
fineMappingResult = fmr,
variantPruningParam = VariantPruningParam(
weightCutoff = 1e-4,
maxVariants = 1000,
csMinCor = 0.8,
pipRescue = 0.1
)
)Reading the result
ctwasPipeline() returns a CtwasResult: one
row per (gwasStudy, study, context, method) run.
methodNames(res)## [1] "susie"
studyName(res)## [1] "protocol_example"
contexts(res)## [1] "bulk_rnaseq"The per-gene posteriors live in ctwasPosteriors():
fm <- ctwasPosteriors(res)
as.data.frame(fm)[, c(
"molecular_id",
"context",
"region_id",
"z",
"susie_pip",
"cs"
)]## molecular_id context region_id z susie_pip cs
## 1 ENSG00000130538 bulk_rnaseq chr22_1 1.826089 0.06802641 NAsusie_pip is the quantity cTWAS exists to produce: the
probability the gene is causal given the SNPs it competes with. Here it
is low (the toy GWAS has no real signal) and cs is
NA because the gene falls in no credible set.
susieAlpha() returns the per-single-effect alpha matrix
behind those PIPs, and the EM group priors sit on each row’s entry:
param <- ctwasGroupPriors(res$entry[[1]])
param$group_prior## bulk_rnaseq|bulk_rnaseq SNP
## 0.01326949 0.02400648
param$group_prior_var## bulk_rnaseq|bulk_rnaseq SNP
## 2.368016 4.549256
param$converged## [1] TRUEBy default the context-agnostic SNP background is dropped from the
result, since it is the null the genes are being scored against rather
than a finding. Pass keepSnps = TRUE to retain it as an
extra row (study = context = "SNP"), which makes the full
cTWAS output reconstructable.
The four chainable steps
ctwasPipeline() is a wrapper over four steps. Call them
directly when you want to inspect an intermediate payload, override the
estimated priors, or re-enter the chain without recomputing what came
before.
inputs <- assembleCtwasInputs(
gwasSumStats = gwasByRegion,
twasWeights = list(ctwasWeightsExample)
)
est <- estCtwasGroupPriors(
inputs,
ctwasPriorParam = CtwasPriorParam(
thin = 1,
niterPrefit = 3,
niter = 10,
fallbackToPrefit = TRUE
),
methodArgs = CtwasOptions(
regionData = CtwasRegionDataOptions(min_group_size = 1),
estParam = CtwasEstParamOptions(
min_group_size = 1,
min_p_single_effect = 0
)
)
)
screened <- screenCtwasRegions(
est,
methodArgs = CtwasOptions(
screen = CtwasScreenOptions(min_nonSNP_PIP = 0)
)
)
finemapped <- finemapCtwasRegions(screened)
nrow(finemapped$finemap_res)## [1] 20finemap_res here has more rows than
ctwasPosteriors(res) above: it still carries the SNP rows.
asCtwasResult() is the step that structures the payload and
drops them.
nrow(ctwasPosteriors(asCtwasResult(finemapped)))## [1] 1
nrow(ctwasPosteriors(asCtwasResult(finemapped, keepSnps = TRUE)))## [1] 20Each step’s payload is also bundled, so you can enter anywhere:
data(ctwasInputsExample) # after assembleCtwasInputs
data(ctwasEstExample) # after estCtwasGroupPriors
data(ctwasFinemapExample) # after finemapCtwasRegionsThe most common reason to use the granular path is to override the
priors. estCtwasGroupPriors()’s accurate-EM step can
diverge on underpowered data; fallbackToPrefit = TRUE
handles that automatically, but taking the steps apart lets you
substitute priors of your own into
screenCtwasRegions().
Boundary genes
A gene whose cis window crosses an LD block boundary is fine-mapped
in only one of the two blocks. mergeCtwasBoundaryRegions()
re-runs such genes over the merged pair:
merged <- mergeCtwasBoundaryRegions(
finemapped,
pipThresh = 0.5,
filterCs = FALSE,
maxSNP = Inf
)ctwasPipeline(boundaryMergeParam = BoundaryMergeParam(enabled = TRUE))
does the same thing inline, with mergePipThresh /
mergeFilterCs / mergeMaxSNP forwarding to
those arguments. It is off by default because it only matters when a
high-PIP gene actually straddles a boundary.
Common parameters
ctwasPipeline(
gwasSumStats = gwasByRegion,
twasWeights = list(ctwasWeightsExample),
method = NULL, # NULL: "ensemble" if present, else every method
maxNumSingleEffects = 5,
ncore = 1,
keepSnps = FALSE,
ctwasPriorParam = CtwasPriorParam(
thin = 0.1,
niterPrefit = 3,
niter = 30,
varStructure = "shared_type",
fallbackToPrefit = FALSE
),
variantPruningParam = VariantPruningParam(),
boundaryMergeParam = BoundaryMergeParam(enabled = FALSE)
)When method is NULL and several weight
methods are present without an "ensemble", each method runs
as its own independent cTWAS run and contributes its own row-set.
Next steps
- TWAS weights — training the weights cTWAS consumes.
-
TWAS Z-scores and MR — the per-gene Z
that
twasZoptionally supplies here. -
QC for summary statistics —
summaryStatsQc(), whichgwasSumStatsmust have been through.
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.3
##
## loaded via a namespace (and not attached):
## [1] mathjaxr_2.0-0 RColorBrewer_1.1-3
## [3] jsonlite_2.0.0 MultiAssayExperiment_1.36.2
## [5] magrittr_2.0.5 GenomicFeatures_1.62.0
## [7] farver_2.1.2 rmarkdown_2.32
## [9] fs_2.1.0 BiocIO_1.20.0
## [11] ragg_1.5.2 vctrs_0.7.3
## [13] memoise_2.0.1 Rsamtools_2.26.0
## [15] RCurl_1.98-1.19 mixsqp_0.3-54
## [17] htmltools_0.5.9 S4Arrays_1.10.1
## [19] BiocBaseUtils_1.12.0 curl_8.0.0
## [21] tictoc_1.2.1 SparseArray_1.10.10
## [23] sass_0.4.10 bslib_0.12.0
## [25] htmlwidgets_1.6.4 desc_1.4.3
## [27] plyr_1.8.9 plotly_4.12.1
## [29] zoo_1.9-0 cachem_1.1.0
## [31] gggrid_0.2-0 GenomicAlignments_1.46.0
## [33] lifecycle_1.0.5 pkgconfig_2.0.3
## [35] Matrix_1.7-6 R6_2.6.1
## [37] fastmap_1.2.0 MatrixGenerics_1.22.0
## [39] digest_0.6.39 numDeriv_2016.8-1.1
## [41] reshape_0.8.10 AnnotationDbi_1.72.0
## [43] S4Vectors_0.48.1 irlba_2.4.1
## [45] textshaping_1.0.5 GenomicRanges_1.62.1
## [47] RSQLite_3.53.3 metadat_1.6-0
## [49] pgenlibr_0.7.0 httr_1.4.9
## [51] abind_1.4-8 compiler_4.5.3
## [53] withr_3.0.3 bit64_4.8.6
## [55] S7_0.2.2 backports_1.5.1
## [57] metafor_5.0-1 BiocParallel_1.44.0
## [59] DBI_1.3.0 MASS_7.3-66
## [61] DelayedArray_0.36.1 rjson_0.2.23
## [63] tools_4.5.3 otel_0.2.0
## [65] glue_1.8.1 quadprog_1.5-8
## [67] restfulr_0.0.17 nlme_3.1-170
## [69] grid_4.5.3 checkmate_2.3.3
## [71] generics_0.1.4 LDlinkR_1.4.0
## [73] gtable_0.3.6 tzdb_0.5.0
## [75] susieR_0.16.6 locuszoomr_0.3.8
## [77] tidyr_1.3.2 ensembldb_2.34.0
## [79] data.table_1.18.6.1 hms_1.1.4
## [81] XVector_0.50.0 BiocGenerics_0.56.0
## [83] ggrepel_0.9.8 pillar_1.11.1
## [85] stringr_1.6.0 vroom_1.7.1
## [87] logging_0.10-111 dplyr_1.2.1
## [89] lattice_0.23-1 rtracklayer_1.70.1
## [91] bit_4.6.0 tidyselect_1.2.1
## [93] Biostrings_2.78.0 knitr_1.52
## [95] IRanges_2.44.0 Seqinfo_1.0.0
## [97] ProtGenerics_1.42.0 SummarizedExperiment_1.40.0
## [99] stats4_4.5.3 xfun_0.61
## [101] Biobase_2.70.0 matrixStats_1.5.0
## [103] stringi_1.8.9 UCSC.utils_1.6.1
## [105] lazyeval_0.2.3 yaml_2.3.12
## [107] evaluate_1.0.5 codetools_0.2-20
## [109] cigarillo_1.0.0 archive_1.1.13
## [111] tibble_3.3.1 cli_3.6.6
## [113] RcppParallel_5.1.11-2 systemfonts_1.3.2
## [115] jquerylib_0.1.4 Rcpp_1.1.2
## [117] GenomeInfoDb_1.46.2 zigg_0.0.2
## [119] png_0.1-9 Rfast_2.1.5.2
## [121] XML_3.99-0.25 parallel_4.5.3
## [123] pkgdown_2.2.1 ggplot2_4.0.3
## [125] readr_2.2.0 blob_1.3.0
## [127] AnnotationFilter_1.34.0 bitops_1.1-0
## [129] viridisLite_0.4.3 scales_1.4.0
## [131] ctwas_0.6.0 purrr_1.2.2
## [133] crayon_1.5.3 rlang_1.3.0
## [135] AMR_3.0.1 cowplot_1.2.0
## [137] KEGGREST_1.50.0