Compute mvSuSiE-RSS TWAS weights from summary statistics

Multi-context summary-statistics analog of mvsusieWeights: extracts coefficients from an existing mvsusieR::mvsusie_rss fit. It never runs fine-mapping itself — mvsusieRssFit is required.

Usage

mvsusieRssWeights(
  stat,
  ldReference = NULL,
  mvsusieRssFit = NULL,
  fitRetention = c("none", "slim", "full")
)

Arguments

stat

A list with z (matrix variants x conditions), n (numeric vector or scalar) and, to check against an ldReference, variantNames, describing the block the fit came from. The fit must cover the same number of variants; a mismatch is an error, since weights from one block are not comparable against another's.

ldReference

Optional LD reference for that block: a genotype panel, a GenotypeMixture or an LdMatrix. With one, the statistics' variants must also be in it. No LD is read.

mvsusieRssFit

A fitted mvsusieRss object. Required: this extracts weights and never runs fine-mapping itself.

fitRetention

How much of the fit is kept on the returned weights: "none" attaches nothing, "slim" and "full" attach it as the "fit" attribute. These extractors hold no trimmable intermediates, so the two retaining levels behave alike here.

Value

A numeric matrix of per-variant per-context weights (variants x conditions).

Details

Follows the RSS weights contract: stat (and ldReference, when given) describe the block the fit is being extracted against, and the fit is checked against them rather than refitted.

Examples

data(eqtlRegionExample)
X <- eqtlRegionExample$X[, 1:30]
y <- eqtlRegionExample$yRes
ss <- lapply(seq_len(ncol(X)), function(j) {
  coef(summary(lm(y ~ X[, j])))[2, 1:2]
})
stat <- list(
  bhat = vapply(ss, `[`, numeric(1), 1L),
  shat = vapply(ss, `[`, numeric(1), 2L),
  z = vapply(ss, function(s) s[1] / s[2], numeric(1)),
  n = rep(nrow(X), ncol(X)))
LD <- cor(X)
fit <- fitMvsusieRss(Z = stat$z, R = LD, N = nrow(X),
  prior_variance = 1)
#> Eigendecomposition cache: K=1, common_cov=TRUE [mem: 1.01 GB]
#> Model initialized: J=30, R=1, L=10, K=1 [mem: 1.01 GB]
#> iter          ELBO       delta   sigma2      mem      V
#>    1     -588.3595           -   diag[1,1]   1.01 GB  [0 x 10]
#>    2     -588.3595    0.00e+00   diag[1,1]   1.01 GB  [0 x 10]  converged
mvsusieRssWeights(stat, mvsusieRssFit = fit)
#>       [,1]
#>  [1,]    0
#>  [2,]    0
#>  [3,]    0
#>  [4,]    0
#>  [5,]    0
#>  [6,]    0
#>  [7,]    0
#>  [8,]    0
#>  [9,]    0
#> [10,]    0
#> [11,]    0
#> [12,]    0
#> [13,]    0
#> [14,]    0
#> [15,]    0
#> [16,]    0
#> [17,]    0
#> [18,]    0
#> [19,]    0
#> [20,]    0
#> [21,]    0
#> [22,]    0
#> [23,]    0
#> [24,]    0
#> [25,]    0
#> [26,]    0
#> [27,]    0
#> [28,]    0
#> [29,]    0
#> [30,]    0