Learning TWAS weights with pecotmr
pecotmr authors
2026-09-08
Source:vignettes/twas-weights.Rmd
twas-weights.RmdOverview
TWAS weights are per-variant coefficients that predict a molecular phenotype (gene expression, splicing, protein abundance) from genotype. Once learned, the weights are combined with GWAS summary statistics to compute per-gene TWAS Z-scores (covered in the TWAS Z-score + causal inference vignette).
twasWeightsPipeline() is a single that accepts multiple
input types. The output is always a TwasWeights collection
(a GRangesList subclass keyed by
(study, context, trait, method), one element per
tuple).
| Input class | Methods supported |
|---|---|
QtlDataset |
Univariate: LASSO, elastic net, mr.ash, MCP, SCAD, L0Learn, Bayesian alphabet, DPR. Multivariate: mr.mash, mvsusie |
MultiStudyQtlDataset |
Iterates over embedded studies; same methods as
QtlDataset
|
QtlSumStats |
RSS-only: LASSO, mr.ash, MCP, SCAD, L0Learn, DPR, PRS-CS, P+T. Multivariate: mr.mash, mvsusie |
Bundled inputs
data(qtlDatasetExample, qtlSumStatsExample, multiStudyQtlDatasetExample)Individual-level weights
The simplest individual-level call learns a single regression-based
weight set from a QtlDataset. Cross-validation is enabled
by default (cvFolds = 5). Like
fineMappingPipeline(), you need to tell it which genotype
window to extract — cisWindow (per trait) is the most
common choice.
Fine-mapping methods (susie, susieInf,
susieAsh, mvsusie, fsusie) are
not re-fit by twasWeightsPipeline() - instead run
fine-mapping models with fineMappingPipeline() first and
pass the result via fineMappingResult = … (see the next
section).
tw <- twasWeightsPipeline(qtlDatasetExample, methods = "lasso", cisWindow = 1e6)
tw## TwasWeights: 1 entries
## 1 studies, 1 contexts, 1 traits, 1 methods
## LD sketch: NULL (individual-level fit)Extract the per-variant weights for one tuple:
entry <- getTwasWeights(
tw,
study = "study1",
context = "brain",
trait = "ENSG_example",
method = "lasso"
)
head(getWeights(entry))## chr22:14850625:T:G chr22:14870204:T:C chr22:14878387:G:A chr22:14880040:A:G
## 0 0 0 0
## chr22:14884399:T:C chr22:15063831:A:G
## 0 0CV performance:
getCvResult(entry)## $samplePartition
## # A tibble: 165 × 2
## Sample Fold
## <chr> <int>
## 1 NA12005 1
## 2 NA11919 1
## 3 NA12842 1
## 4 NA12843 1
## 5 NA10836 1
## 6 NA10831 1
## 7 NA12145 1
## 8 NA12761 1
## 9 NA12286 1
## 10 NA07349 1
## # ℹ 155 more rows
##
## $predictions
## NA06989 NA11891 NA11843 NA12341 NA12739 NA10850
## -0.049846523 0.071101616 0.000000000 0.000000000 0.000000000 0.000000000
## NA06984 NA12877 NA12275 NA06986 NA12272 NA10845
## 0.000000000 -0.173111627 0.055076902 0.153867473 -0.011121075 0.000000000
## NA10852 NA07051 NA12400 NA12344 NA12777 NA12287
## -0.050745478 -0.067440733 0.142604947 0.369153869 0.000000000 -0.045874154
## NA10837 NA12383 NA12340 NA12708 NA12273 NA11892
## 0.025284031 0.000000000 0.000000000 0.000000000 0.000000000 -0.069997549
## NA12546 NA12843 NA12766 NA12348 NA12817 NA10840
## -0.146018540 -0.241263924 -0.049865421 0.000000000 -0.082630013 -0.035957498
## NA06997 NA11917 NA12718 NA12282 NA11920 NA12776
## 0.000000000 -0.117916165 -0.050317441 -0.253274971 0.056337923 -0.035570735
## NA12283 NA07435 NA12828 NA07045 NA07031 NA12336
## -0.084527847 -0.030963247 -0.083894490 -0.051725971 -0.048466511 -0.261175909
## NA07349 NA12827 NA12375 NA12343 NA12335 NA12778
## -0.284366433 0.000000000 -0.049903000 0.000000000 0.057216412 0.000000000
## NA12832 NA11930 NA12890 NA07037 NA07347 NA12829
## 0.056966310 0.034463620 0.133021411 0.000000000 -0.074719952 -0.203697418
## NA12749 NA11894 NA12286 NA10865 NA10864 NA10853
## -0.230264107 -0.034478174 -0.300843685 0.071803864 -0.091158897 0.376585622
## NA11918 NA12830 NA12818 NA10836 NA11893 NA11919
## 0.055161238 0.000000000 -0.242491863 0.118303283 0.572821006 0.148243677
## NA12489 NA12399 NA12413 NA10843 NA12842 NA12347
## -0.077314964 0.049445481 0.000000000 0.000000000 0.062382803 0.369890590
## NA07346 NA12775 NA07014 NA12767 NA12889 NA12386
## -0.048693439 -0.010530936 -0.250633007 -0.033740477 0.143913300 0.000000000
## NA06995 NA12376 NA12342 NA12748 NA11931 NA12045
## -0.081141571 0.050337080 0.161343431 0.072509753 -0.050462020 0.048452362
## NA12750 NA11831 NA12146 NA11882 NA07056 NA12707
## 0.000000000 -0.061883388 0.333170268 0.156807277 0.056531721 -0.067586148
## NA12154 NA12753 NA11839 NA10859 NA12875 NA07348
## 0.055500400 -0.055547196 -0.052454445 0.160684270 -0.019321064 -0.034752152
## NA12156 NA12044 NA11992 NA11829 NA12239 NA12762
## -0.255985044 -0.071319524 -0.050392820 -0.057446110 -0.064590106 0.145925671
## NA12716 NA12878 NA10856 NA12874 NA12760 NA06985
## 0.000000000 -0.229220988 0.000000000 -0.034969029 0.274506212 -0.157753584
## NA12003 NA10835 NA07022 NA12813 NA10839 NA07055
## -0.050512630 0.000000000 -0.104197463 0.132423240 0.139063582 0.070882514
## NA12056 NA10863 NA12145 NA12814 NA10847 NA12006
## 0.000000000 0.112059304 -0.041282450 0.323978451 0.116310179 0.298310330
## NA12763 NA07357 NA12144 NA10831 NA07000 NA11832
## 0.000000000 -0.172491997 0.000000000 -0.044011218 -0.035746359 -0.036139064
## NA06991 NA11840 NA12802 NA12761 NA10830 NA10855
## -0.228693356 -0.146148277 0.240825549 -0.245041556 -0.263423707 -0.035585131
## NA06994 NA11993 NA11995 NA12891 NA12864 NA12751
## 0.000000000 -0.034642781 -0.072197205 0.071112478 -0.079599677 -0.094899456
## NA10861 NA12005 NA12234 NA07345 NA07029 NA12892
## -0.296682755 -0.256612254 -0.054552829 -0.034484954 0.537202025 -0.036115202
## NA12248 NA10846 NA12801 NA12872 NA12155 NA06993
## 0.000000000 0.000000000 -0.073877474 -0.004869147 -0.009002766 -0.069321802
## NA11830 NA10838 NA12249 NA12057 NA12812 NA11881
## -0.053694877 0.000000000 0.104125173 -0.066953967 -0.084943961 -0.254625615
## NA11994 NA12873 NA12815 NA12740 NA12752 NA12043
## -0.124120295 -0.072020180 0.273167259 -0.069612209 -0.220653149 -0.057500318
## NA12264 NA10854 NA12865
## 0.000000000 0.020780484 -0.049812528
##
## $metrics
## corr rsq adj_rsq pval RMSE MAE
## 0.051390763 0.002641010 -0.003477756 0.512119380 0.999951442 0.790754739
##
## $foldFits
## NULLThe entry’s variantIds slot lists the LD-aligned variant
set the weights are defined on:
length(getVariantIds(entry))## [1] 198Including a fine-mapping susie fit in the TWAS ensemble
The posterior effect estimates of SuSiE fine-mapping models are
themselves a set of per-variant weight estimates. When you pass an
existing FineMappingResult to
twasWeightsPipeline(), the trimmed susie fit for each
matching (study, context, trait) tuple is added to the
ensemble as an additional weights row instead of being re-fit.
fmr <- fineMappingPipeline(
qtlDatasetExample,
methods = "susie",
cisWindow = 1e6
)
tw2 <- twasWeightsPipeline(
qtlDatasetExample,
methods = c("susie", "lasso"),
cisWindow = 1e6,
fineMappingResult = fmr
)Multi-study weights
msTw <- twasWeightsPipeline(
multiStudyQtlDatasetExample,
methods = "lasso",
cisWindow = 1e6
)
table(msTw$study)##
## study1 study2
## 1 1The call for iterating over the studies of a MultiStudyQtlDataset is identical to the call for a single QtlDataset.
Multiple weight methods at once
Pass a character vector to fit several methods on the same fold split
and return them as separate (method) rows in the
collection. Fine-mapping methods (susie,
susieInf, susieAsh, mvsusie) can
only be included if their results are present in the object passed to
the fineMappingResult argument:
# TWAS-regression only — no fineMappingResult needed.
twasWeightsPipeline(
qtlDatasetExample,
methods = c("lasso", "enet", "ridge", "mrash")
)
# Add a fine-mapping susie row by providing the FineMappingResult.
fmr <- fineMappingPipeline(
qtlDatasetExample,
methods = "susie",
cisWindow = 1e6
)
twasWeightsPipeline(
qtlDatasetExample,
methods = c("susie", "lasso", "enet", "ridge", "mrash"),
cisWindow = 1e6,
fineMappingResult = fmr
)RSS / sumstats weights
twSs <- twasWeightsPipeline(qtlSumStatsExample, methods = "lasso")
twSs## TwasWeights: 1 entries
## 1 studies, 1 contexts, 1 traits, 1 methods
## LD sketch: plink1 @ pecotmr://extdata/toy_canonicalThe RSS path uses the ldSketch genotype panel inside the
QtlSumStats collection to compute LD as needed. The
same
Multivariate (multi-context) weights
When you have multiple contexts for the same trait (e.g. multi-tissue
expression), pass methods = "mrmash" to do a joint
multivariate fit:
twasWeightsPipeline(
qtlDatasetExample,
methods = "mrmash",
contexts = c("brain", "liver")
)For an explicit cross-context joint fit driven by a
jointSpecification:
twasWeightsPipeline(
multiStudyQtlDatasetExample,
methods = "mrmash",
jointSpecification = list(axis = "context", contexts = c("brain", "liver"))
)Common parameters
twasWeightsPipeline(
qtlDatasetExample,
methods = "susie",
contexts = "brain",
traitId = "ENSG_example",
region = "chr22:25000000-26000000",
cisWindow = 5e5,
cvFolds = 5,
samplePartition = NULL,
maxCvVariants = -1,
cvThreads = 1,
ensemble = TRUE,
ensembleR2Threshold = 0.01,
ensembleSolver = "quadprog",
estimatePi = TRUE,
phenotypeCovariatesToResidualize = NULL,
genotypeCovariatesToResidualize = NULL,
dataType = NULL,
naAction = "drop"
)-
ensemble = TRUEadds an additionalensembleWeightsrow to the collection — a weighted combination of the individual methods’ CV predictions filterd byensembleR2Threshold. -
naActionis forwarded togetResidualizedPhenotypes()(see the fine-mapping vignette for semantics).
Reading the result
TwasWeights is a GRangesList subclass keyed
by the 4-tuple (study, context, trait, method) — one
element per tuple, holding the variants its weights are aligned to. The
tuple and the per-tuple payloads are mcols() columns,
reached with $:
tw## TwasWeights: 1 entries
## 1 studies, 1 contexts, 1 traits, 1 methods
## LD sketch: NULL (individual-level fit)
colnames(tw)## [1] "study" "context" "trait" "method" "fits"
## [6] "cvResult" "standardized" "dataType" "traitPos"
tw$method## [1] "lasso"as.data.frame() gives the long per-variant view across
every method at once, which is usually what you want for inspection:
head(as.data.frame(tw), 3)## group group_name seqnames start end width strand A1 A2 weight
## 1 1 <NA> chr22 14850625 14850625 1 * G T 0
## 2 1 <NA> chr22 14870204 14870204 1 * C T 0
## 3 1 <NA> chr22 14878387 14878387 1 * A G 0To pull one method’s weights, pass the tuple to the accessor:
key <- list(
study = tw$study[[1L]],
context = tw$context[[1L]],
trait = tw$trait[[1L]],
method = tw$method[[1L]]
)
w <- do.call(getWeights, c(list(tw), key))
length(w)## [1] 198## [1] FALSE## NULLgetWeights() returns a named numeric vector for
univariate methods and a matrix for multivariate ones. The remaining
payloads are keyed the same way: getVariantIds() for the
aligned variant set, getCvResult() for the cross-validation
record (rsq, r, predicted values), and
getFits() for the full fit object when the pipeline was run
with retainFit = TRUE. resolveWeights()
returns the weights and their variant IDs together, which is the form
the downstream TWAS and MR steps consume.
Next steps
- Combine the learned weights with a GWAS via
causalInferencePipelineto compute per-gene TWAS Z-scores.
Session info
## 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.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.7.10
##
## 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 survival_3.8-11
## [13] magrittr_2.0.5 compiler_4.5.3
## [15] rlang_1.3.0 sass_0.4.10
## [17] tools_4.5.3 utf8_1.2.6
## [19] yaml_2.3.12 knitr_1.51
## [21] S4Arrays_1.10.1 htmlwidgets_1.6.4
## [23] bit_4.6.0 DelayedArray_0.36.0
## [25] plyr_1.8.9 RColorBrewer_1.1-3
## [27] abind_1.4-8 BiocParallel_1.44.0
## [29] withr_3.0.3 purrr_1.2.2
## [31] numDeriv_2016.8-1.1 BiocGenerics_0.56.0
## [33] desc_1.4.3 grid_4.5.3
## [35] stats4_4.5.3 susieR_0.16.6
## [37] ggplot2_4.0.3 iterators_1.0.14
## [39] scales_1.4.0 MASS_7.3-66
## [41] MultiAssayExperiment_1.36.1 SummarizedExperiment_1.40.0
## [43] cli_3.6.6 rmarkdown_2.32
## [45] metafor_5.0-1 crayon_1.5.3
## [47] ragg_1.5.2 generics_0.1.4
## [49] otel_0.2.0 RcppParallel_5.1.11-2
## [51] tzdb_0.5.0 BiocBaseUtils_1.12.0
## [53] cachem_1.1.0 stringr_1.6.0
## [55] splines_4.5.3 metadat_1.6-0
## [57] parallel_4.5.3 XVector_0.50.0
## [59] matrixStats_1.5.0 vctrs_0.7.3
## [61] glmnet_5.0 Matrix_1.7-6
## [63] jsonlite_2.0.0 IRanges_2.44.0
## [65] hms_1.1.4 S4Vectors_0.48.0
## [67] bit64_4.8.6 mixsqp_0.3-54
## [69] irlba_2.3.7 archive_1.1.13
## [71] systemfonts_1.3.2 foreach_1.5.2
## [73] tidyr_1.3.2 jquerylib_0.1.4
## [75] snpStats_1.60.0 glue_1.8.1
## [77] pkgdown_2.2.1 codetools_0.2-20
## [79] shape_1.4.6.1 stringi_1.8.9
## [81] gtable_0.3.6 GenomicRanges_1.62.1
## [83] quadprog_1.5-8 tibble_3.3.1
## [85] pillar_1.11.1 htmltools_0.5.9
## [87] Seqinfo_1.0.0 R6_2.6.1
## [89] zigg_0.0.2 textshaping_1.0.5
## [91] vroom_1.7.1 evaluate_1.0.5
## [93] lattice_0.23-1 Biobase_2.70.0
## [95] readr_2.2.0 tictoc_1.2.1
## [97] Rsamtools_2.26.0 Rfast_2.1.5.2
## [99] bslib_0.12.0 Rcpp_1.1.2
## [101] SparseArray_1.10.8 nlme_3.1-170
## [103] xfun_0.60 fs_2.1.0
## [105] MatrixGenerics_1.22.0 pkgconfig_2.0.3