Predicting genomic breeding values

Practical 5: Genomic relationships, GBLUP and held-out prediction

Learning objectives

Align genotype and phenotype identities; define a genomic relationship scale; compare pedigree BLUP and GBLUP; and distinguish model agreement, in-sample fit and held-out phenotype prediction.

Continue from Practical 4. Use base R, qgg and the shared helper script. The analyses use the existing small mouse teaching dataset; they are not a benchmark or breeding-strategy qualification.

1. Inspect and align imputed genotypes

source("../scripts/mouse_practical_helpers.R")
data <- mouse_pedigree_data()
mouse <- data$mouse
ids <- data$ids
genotypes <- mouse_file("genotypes_imputed",
  "cd21078d4def9df35e97789396c730e1e73c3958",
  "c3b04ffca04a05f9b8da248873eca689")
stopifnot(is.matrix(genotypes), !is.null(rownames(genotypes)),
          !is.null(colnames(genotypes)), !anyDuplicated(rownames(genotypes)),
          !anyDuplicated(colnames(genotypes)), all(ids %in% rownames(genotypes)),
          all(is.finite(genotypes)), all(genotypes >= 0 & genotypes <= 2))
dim(genotypes)

[1] 1267 1813

range(genotypes)

[1] 0 2

M <- genotypes[ids,,drop=FALSE]
stopifnot(identical(rownames(M), ids))

The source contains 1,267 individuals and 1,813 marker columns. Values are supplied imputed dosages on a 0–2 scale; missing calls have already been handled upstream. This file does not establish the imputation procedure, physical marker metadata or allele orientation. Column labels are identifiers, not verified genome positions. We use this supplied representation for a teaching calculation without inferring those missing details.

Question 1. Why can matching dimensions still produce an incorrect analysis?

Suggested answer

Row orders may differ. A phenotype can be assigned another individual’s relationships even when matrix dimensions agree. We select rows by checked character IDs and verify the resulting order. Marker column ordering must also be retained when applying a training transformation to validation genotypes.

2. Define the genomic relationship matrix

For this tutorial, standardize each retained marker using the phenotype sample’s mean and sample standard deviation, and define

\[W_{ij}=\frac{M_{ij}-\bar M_j}{s_j},\qquad G=WW^T/m.\]

Markers with zero variance are excluded. Every retained standardized marker receives equal weight. This is an explicit standardized-marker kernel, not the same normalization as a pooled \(2p_j(1-p_j)\) denominator. Its base is the specified reference sample, rather than the pedigree founder base.

marker_sd <- apply(M,2,sd)
retained <- is.finite(marker_sd) & marker_sd > 0
W <- scale(M[,retained,drop=FALSE])
G <- tcrossprod(W)/ncol(W)
dimnames(G) <- list(ids,ids)
A <- mouse_relationship(data$pedigree)
K <- A[ids,ids,drop=FALSE]
c(retained_markers=ncol(W), excluded_markers=sum(!retained))

retained_markers excluded_markers 
            1813                0 

summary(diag(G))

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
 0.2624  0.9351  0.9967  0.9992  1.0615  1.4202 

stopifnot(isTRUE(all.equal(G,t(G),tolerance=1e-10)))
# A transparent tiny kernel check: two individuals, one standardized marker.
stopifnot(max(abs(tcrossprod(matrix(c(-1,1),2,1))-
                   matrix(c(1,-1,-1,1),2,2))) < 1e-12)

This dense calculation is appropriate for the teaching file. Larger workflows should use supported genomic resources and bounded computation. A centered \(G\) can be singular; positive residual variance makes \(V=G\sigma_G^2+I\sigma_E^2\) invertible without adding an arbitrary ridge to \(G\).

local({
  old <- par(mfrow=c(1,2)); on.exit(par(old))
  image(1:40,1:40,K[1:40,1:40],main="Pedigree A",
        xlab="Record index",ylab="Record index",col=hcl.colors(20,"Blues"))
  image(1:40,1:40,G[1:40,1:40],main="Genomic G",
        xlab="Record index",ylab="Record index",col=hcl.colors(20,"Blue-Red 3"))
})

Figure 1: Pedigree and genomic relationships among the same first 40 records. Each panel has its own color scale.

Question 2. Why can \(G\) contain negative off-diagonal entries?

Suggested answer

Genomic similarities are centered relative to the chosen sample. A negative value indicates less similarity than that reference, not negative biological kinship. Pedigree \(A\) measures expected relationships relative to assumed founders; its entries and \(G\) entries are not automatically on interchangeable bases.

3. Fit and compare evaluations

GBLUP uses \(u\sim N(0,G\sigma_G^2)\) with the same sex/replicate fixed effects and independent residual model as before. Pedigree BLUP uses \(A\) instead. Fit each separately, on identical records.

X <- model.matrix(~sex+reps,data=mouse)
y <- mouse$BW
names(y) <- ids
fitA <- mouse_fit(y,X,K,"A")
fitG <- mouse_fit(y,X,G,"G")
knitr::kable(rbind(pedigree=fitA$theta, genomic=fitG$theta),digits=3,
             col.names=c("Relationship scale","Residual variance"))

	Relationship scale	Residual variance
pedigree	6.570	13.384
genomic	8.255	11.106

c(pedigree_iterations=fitA$niter, genomic_iterations=fitG$niter)

pedigree_iterations  genomic_iterations 
                  9                  15 

uA <- drop(fitA$g[,"A"])
uG <- drop(fitG$g[,"G"])
cor(uA,uG)

[1] 0.7519603

plot(uA,uG,pch=16,cex=0.5,col="#247BA0",
     xlab="Pedigree breeding-value prediction",
     ylab="Genomic breeding-value prediction")
abline(0,1,lty=2)

Figure 2: Agreement between pedigree and genomic evaluations fitted using the same phenotypes; this is not independent validation.

Question 3. Does higher agreement with pedigree predictions prove better genomic prediction?

Suggested answer

No. Both evaluations use the same records, and neither is observed truth. Differences reflect relationship information, reference scale and fitted model components. Compare prediction performance on records excluded from fitting, with a split that reflects the intended use.

4. Predict a held-out replicate

Use replicates 1–2 for training and replicate 3 for validation. This illustrates prediction of phenotypes in a held-out replicate, with relatives potentially present in training. It is not validation in unrelated families or a future-generation experiment.

Replicate 3 has no training observations from which to estimate its fixed effect. For a declared baseline, fit sex only in this validation exercise. An unmodelled replicate shift may affect calibration. Recompute marker means and scales using training genotypes only, and fit variance components using training phenotypes only.

training <- as.character(mouse$reps) != "3"
validation <- !training
Xv <- model.matrix(~sex,data=mouse)
Mt <- M[training,,drop=FALSE]
center_t <- colMeans(Mt)
scale_t <- apply(Mt,2,sd)
keep_t <- is.finite(scale_t) & scale_t > 0
Wt <- sweep(sweep(Mt[,keep_t,drop=FALSE],2,center_t[keep_t],"-"),
            2,scale_t[keep_t],"/")
Wv <- sweep(sweep(M[validation,keep_t,drop=FALSE],2,center_t[keep_t],"-"),
            2,scale_t[keep_t],"/")
Gtt <- tcrossprod(Wt)/ncol(Wt)
Gvt <- Wv %*% t(Wt)/ncol(Wt)
dimnames(Gtt) <- list(ids[training],ids[training])
fit_train <- mouse_fit(y[training],Xv[training,,drop=FALSE],Gtt,"G")
VG_t <- unname(fit_train$theta["G"])
VE_t <- unname(fit_train$theta["E"])
Vt <- VG_t*Gtt + diag(VE_t,sum(training))
residual_t <- y[training]-drop(Xv[training,,drop=FALSE] %*% fit_train$b)
uhat_validation <- drop(VG_t*Gvt %*% solve(Vt,residual_t))
predicted_validation <- drop(Xv[validation,,drop=FALSE] %*% fit_train$b) +
                        uhat_validation
baseline_validation <- predict(lm(BW~sex,data=mouse[training,]),
                               newdata=mouse[validation,])
validation_results <- data.frame(
  model=c("Training sex means","Training GBLUP"),
  RMSE=c(sqrt(mean((y[validation]-baseline_validation)^2)),
         sqrt(mean((y[validation]-predicted_validation)^2))),
  correlation=c(cor(y[validation],baseline_validation),
                cor(y[validation],predicted_validation)))
knitr::kable(validation_results,digits=3)

model	RMSE	correlation
Training sex means	5.421	0.676
Training GBLUP	4.762	0.788

c(training_records=sum(training), validation_records=sum(validation),
  training_fit_iterations=fit_train$niter)

       training_records      validation_records training_fit_iterations 
                    842                     335                      15 

plot(predicted_validation,y[validation],pch=16,cex=0.6,col="#247BA0",
     xlab="Held-out predicted body weight",ylab="Observed body weight")
abline(0,1,lty=2)

Figure 3: Body-weight predictions for replicate 3 using only replicate 1–2 phenotypes for training.

These metrics describe held-out phenotype prediction, including sex differences. They are not correlations with true breeding values. A single split does not establish general superiority or calibrated uncertainty. Upstream imputation was supplied; we cannot verify whether that preprocessing was independent of the split.

Question 4. Why report RMSE as well as correlation?

Suggested answer

Correlation describes association and can remain high despite an offset or incorrect scale. RMSE measures errors in trait units. Family overlap, replicate shifts and the target population determine what this exercise can support.

Exercises

1. Compare the genomic and pedigree prediction rankings. Are the highest-ranked mice identical?
2. Count how many validation sires also occur in training. What deployment question does this split address?
3. Design a validation split for predicting new families. State what preprocessing and model fitting must be repeated, without running a broad validation campaign.

Suggested solutions

evaluations <- data.frame(id=ids,pedigree=uA,genomic=uG)
head(evaluations[order(-evaluations$genomic),])

       id pedigree  genomic
212   212 5.514750 8.527394
447   447 2.558604 7.452242
254   254 4.671404 7.131224
204   204 2.253301 6.983472
1144 1144 2.395455 6.490302
420   420 3.739838 6.368773

validation_sires <- unique(as.character(mouse$sire[validation]))
training_sires <- unique(as.character(mouse$sire[training]))
c(validation_sires=length(validation_sires),
  shared_sires=sum(validation_sires %in% training_sires))

validation_sires     shared_sires 
              11               11 

Rank differences follow differences in information and model bases. With shared sires, the split examines prediction across replicates in related material. For new-family prediction, assign whole families to validation, define treatment of connected relatives, recompute preprocessing on training information and re-estimate variance components. Compare appropriate baselines and retain the intended deployment assumptions.

From evaluation to breeding decisions

The five practicals progress from phenotype exploration to marker effects, pedigree variance estimation, BLUP and genomic evaluation. Selection objectives, economic weights and mating decisions are separate tasks for gsel; biological offspring simulation belongs in gsim. Real-data predicted breeding values remain distinct from simulated truth.