Marker effects recovered from a fitted relationship-based term

Back-transforms the BLUPs of a random term fitted with a genomic relationship matrix (GBLUP) into BLUPs of marker effects, optionally with standard errors and back-solved GWAS-type tests. This is the post-fit equivalent of ASReml's mef argument and does not require refitting the model.

Usage

meffects_mmes(object, term, M, method = c("A.mat", "custom"),
              min.MAF = 0, Z = NULL, scale = NULL, blend = 0, Gu = NULL,
              se = FALSE, chunk = 1000)

Arguments

object: A fitted mmes model.
term: Name (as in names(object$uList)) or position of the random term fitted with Gu.
M: Marker matrix (individuals by markers, coded -1/0/1 as for A.mat) with rownames matching the levels of the term. Use the same matrix (same rows) that was used to build the relationship matrix, because allele frequencies are computed from it. Rows that are not levels of the term are ignored.
method: "A.mat" reproduces exactly how A.mat builds G from M (imputation, min.MAF and monomorphic filters, centring by $2p$ and scaling by $2\sum p_jq_j$). "custom" uses Z and scale.
min.MAF: The min.MAF used in A.mat.
Z: For method="custom": transformed marker matrix with rownames such that $G = scale \cdot ZZ^\top$.
scale: For method="custom": the scalar $s$ in $G = sZZ^\top$.
blend: $\lambda$ if the fitted matrix was $(1-\lambda)G+\lambda I$. Small bending constants added only for invertibility can be left at 0.
Gu: Optional relationship matrix (or inverse, with attr(Gu,"inverse")=TRUE). By default the inverse used in the fit is rebuilt from the model inputs.
se: If TRUE, add standard errors, z statistics and p-values.
chunk: Number of markers processed per block when se=TRUE.

Details

If the fitted relationship matrix is $G=(1-\lambda)sZZ^\top+\lambda I$, the BLUP of the marker effects is $$\hat a=(1-\lambda)s\,Z^\top G^{-1}\hat u,$$ using the same $G^{-1}$ as the fit, so $Z\hat a=\hat u-\lambda G^{-1}\hat u$ (exactly $\hat u$ when $\lambda=0$). For structured terms such as vsm(dsm(env), ism(id), Gu=Gi) effects are returned per coordinate (environment, trait, ...).

With se=TRUE, following Gualdron-Duarte et al. (2014), $$\mathrm{Var}(\hat a)=H^\top(\sigma^2_c G - C^{uu})H,\qquad H=(1-\lambda)sG^{-1}Z,$$ where $\sigma^2_c$ is the variance of the coordinate and $C^{uu}$ the prediction error covariance of the term; $z=\hat a/\mathrm{SE}$ gives back-solved GWAS tests. Only diagonal elements are formed: the $C^{uu}$ part uses solves with the stored coefficient matrix in blocks of chunk markers. The cost grows as $O(mn^2)$ for m markers and n individuals.

The rebuild of $G^{-1}$ requires the data used in the original call; supply Gu otherwise. Markers removed by the filters are returned as NA.

Value

A data frame with columns marker, coordinate (structured terms only), effect and, if se=TRUE, se, z and p.value. Attributes scale and markersKept.

References

Gualdron-Duarte JL, Cantet RJC, Bates RO, Ernst CW, Raney NE, Steibel JP (2014). Rapid screening for phenotype-genotype associations by linear transformations of genomic evaluations. BMC Bioinformatics 15: 246.

Examples

set.seed(7)
M <- matrix(sample(c(-1, 0, 1), 40 * 80, TRUE), 40)
rownames(M) <- paste0("g", 1:40); colnames(M) <- paste0("snp", 1:80)
G <- 0.95 * A.mat(M) + 0.05 * diag(40)
Gi <- solve(G); attr(Gi, "inverse") <- TRUE
d <- data.frame(id = factor(rep(rownames(M), each = 2), levels = rownames(M)))
d$y <- rnorm(80)
fit <- mmes(y ~ 1, random = ~ vsm(ism(id), Gu = Gi), data = d, verbose = FALSE)
head(meffects_mmes(fit, 1, M, blend = 0.05, se = TRUE))