Single-step SNP-BLUP: markers as equations, no G, conjugate gradients

Fits the equivalent model of the single step with the marker effects as unknowns: u_g = Z g + p_g, with g ~ N(0, I (1-rpg) s2u / k) and the residual polygenic p_g ~ N(0, rpg s2u A22), k = 2 sum p(1-p). The implied genomic relationship is G* = (1-rpg) Z Z' / k + rpg A22, WITHOUT the affine adjustment the genotypes= path applies; in populations far from the base the two paths differ by construction, and that difference is declared, not hidden. Every application of A22^-1 uses the identity A22^-1 v = A^22 v - A^21 (A^11)^-1 A^12 v over the sparse blocks of A^-1, with one sparse factorization of the non-genotyped block.

Usage

snp_blup(
  formula,
  data,
  pedigree,
  genotypes,
  theta,
  rpg = 0.05,
  missing_code = NULL,
  tol = 1e-08,
  maxiter = 2000L,
  metafounders = NULL,
  gamma = NULL,
  verbose = interactive()
)

Arguments

formula

as in model(). The markers enter every relationship group and every component of it, as H^-1 does on the genotypes= path: direct-maternal in one group, a reaction norm, direct and maternal in separate groups, the indirect effect (with its dilution =, the mate weights (n_i - 1)^(-d) described in model()). A group of dimension q with covariance K0 gets q sets of marker effects with covariance K0 (1-rpg) / k between them, so the precision is K0^-1 times the scalar one

data

data.frame

pedigree

data.frame animal, sire, dam (required: the model is single step). The ids are matched to the data by the rule of model() (see pedigree there)

genotypes

list with ids and m (0/1/2 matrix); NA is imputed with the mean. Checked as in model(): a repeated id, rownames that contradict ids or an animal below a 0.90 call rate are refused, and a parent-offspring conflict with the pedigree is a warning

theta

variance components, in the order model() reports for this formula (each random group, then the residual). This is a solver at fixed components: estimate them with model() (pedigree or single step) and solve here at scale

rpg

residual polygenic proportion, in (0, 1); plays the role the blend plays in the genotypes= path

missing_code

missing-value code for the trait

tol

relative residual of the conjugate gradients at which to stop

maxiter

maximum conjugate-gradient iterations

metafounders

as in model()

gamma

as in model()

verbose

print the solve as it walks, the relative residual of every PCG iteration, so a long solve is a progress report instead of silence. Defaults to interactive(), live in a session and quiet in scripts and checks. Every fitter also honors Ctrl+C now

Value

list with b (fixed-effect solutions, named term=level; the parametrization note of model() applies), ebv (per covariance group, all animals; ebv() works on it), g (marker effects in trait units per allele dose, NA for monomorphic markers: a vector with one relationship component, a marker x component matrix with more, the columns named as the components in theta), converged, iters, resnorm, message

Details

The structure checks of model() (see its section "Structure checks") run here as warnings: the components are given, so the solve exists.

References

Liu, Z., Goddard, M.E., Reinhardt, F. & Reents, R. (2014). A single-step genomic model with direct estimation of marker effects. Journal of Dairy Science 97:5833-5850.

Masuda, Y. et al. (2017). Avoiding the direct inversion of the numerator relationship matrix... Journal of Animal Science 95:49-52.

Vandenplas, J., Eding, H., Calus, M.P.L. & Vuik, C. (2018). Genetics Selection Evolution 50:51; Vandenplas, J. et al. (2019) 51:30.