Covariance: The Heart of Genomic Prediction

Every breeder who uses genomic selection has seen it happen. A line that was never planted receives a breeding value, and real decisions are made on it: which lines to advance, which to cross, which to drop.

Most explanations start with the software, the markers or ridge regression. We want to start one step earlier, with the question underneath: what makes it possible to predict anything at all?

The answer is covariance, and covariance has a history. It begins as alleles copied from parents to offspring, becomes a probability, then a matrix, and finally the engine of a mixed model. Here is that story in three acts.

1. Act I · Figure 1Resemblance is shared historyShared alleles, coancestry, Wright's paths, A
2. Act II · Figure 2Markers read what really happenedThe same pedigree, now G from DNA
3. Act III · Figure 3One letter changes everythingData, BLUEs and BLUPs, then GEBVs

Resemblance is shared history

Start with one locus and a small pedigree (Fig. 1a). Grandparents A and B have two children, C and D, who are full sibs. Each marries into an unrelated family, E and F, and each has one child: i and j, first cousins.

Every allele in this picture is a physical copy of an ancestral allele. Follow A1. Grandparent A passed it to both C and D, and each of them passed it on. The cousins now carry copies of the same piece of ancestral DNA. This is identity by descent (IBD). Two plants can also carry the same allele by coincidence, from unrelated ancestors. That is identity by state, and it creates no lasting resemblance.

Now zoom out from one locus to whole chromosomes (Fig. 1b). Each bar is one chromosome copy, labelled by the grandparental chromosome it came from. Where the copy i received from C and the copy j received from D descend from the same grandparental chromosome, the cousins carry identical DNA (outlined in pink). On chromosome 1 they share a long stretch of B1; on chromosome 3, a piece of A1; on chromosome 4, nothing. Across the genome they share 17.8%, where the pedigree expected 25%.

Any gene for yield inside a pink segment is the same copy in both plants. Whatever it does to i, it does to j. Their genetic values rise and fall together: that is covariance.

Every pair has a relationship, not only the cousins. The coefficient of coancestry, f, is the probability that one allele drawn at random from each of two plants is identical by descent (Malécot, 1948). Twice that is the additive relationship, a = 2f: the share of alleles the two hold by descent. Sewall Wright (1922) showed how to read it straight from a pedigree: find every path that joins the two plants through a common ancestor, and add (1/2)ⁿ for each path of n steps (Fig. 1c). Parent A and child C are one step apart: 1/2. Full sibs C and D are joined through A and through B, two steps each: 1/4 + 1/4 = 1/2. The cousins are joined by two paths of four steps: 1/16 + 1/16 = 1/8. E and F have no common ancestor: 0.

A is a covariance matrix. Write a for every pair, with 1 on the diagonal for each plant with itself, and you have the numerator relationship matrix A (Fig. 1d). Each cell is a covariance in waiting. Multiply by the additive genetic variance and it becomes the covariance of breeding values:

Covariance is not a property of a trait, or of a plant. It is a property of shared history. The pedigree supplies the history; the trait supplies the scale.

Markers read the history that actually happened

A is an expectation. It gives every pair of first cousins 1/8 and every pair of full sibs 1/2, because that is what happens on average. Real pairs drift away from those values, because each parent passes a random half of its genome to each child. This is Mendelian sampling, and the pedigree cannot see it.

So go back to the same pedigree and genotype every plant (Fig. 2a). Each plant now carries a barcode of 2,000 SNPs, scored as 0, 1 or 2 copies of an allele. From these barcodes we build a second relationship matrix, G (Fig. 2b). Its recipe, from VanRaden (2008), compares the barcodes of every pair of plants and adds up the alleles they share over all markers, then scales the result so that it can be read on the same scale as A.

The important point is not the arithmetic but the identity of the object. G is a relationship matrix, and therefore a covariance matrix, exactly like A. Only the source has changed: pedigree records for A, DNA for G.

Compare the two matrices pair by pair (Fig. 2c). Parent and offspring stay at about one half (0.51). The full sibs C and D share 0.46 rather than 0.50. The cousins share 0.09 rather than 0.125, which is exactly the shortfall we saw on their chromosomes in Figure 1b. Unrelated E and F stay at zero. G keeps the family structure of A and replaces each expectation with a measurement.

For a breeder, this is where markers earn their cost. Full sibs from one cross all share the same A, so the pedigree cannot tell them apart. G can, before any of them has been planted.

A says who should be alike. G says who is alike.

One letter changes everything: from I to G

Now put the covariance to work in a familiar setting (Fig. 3a). Two inbred parents are crossed, and six lines are advanced to F5 by single-seed descent. L1 and L5 come from the same F3 plant, as do L3 and L6; L2 and L4 come from other families. All six are genotyped, giving the 6 × 6 G matrix in panel a. The trial is a randomized complete block design with three blocks (replications).

Route 1, the classical trial (Fig. 3b). All six lines are planted, 18 plots in all. Before any model runs, the data are arranged as three objects: the vector y of plot yields, the design matrix X that records which block each plot sits in, and the design matrix Z that records which line each plot belongs to. Mixed-model software such as ASReml-R or sommer then fits y = Xβ + Zu + e. The block effects β are fixed, and their estimates are BLUEs: 4.78, 5.18 and 5.35 t/ha. The line effects u are random with u ~ N(0, Iσ²_u), and their predictions are BLUPs.

The identity matrix I makes a strong statement: every line is unrelated to every other line. Each BLUP is built only from that line's own plots, shrunk toward the mean by how much of its advantage could be plot noise.

Route 2, the genomic trial (Fig. 3c). To save field space, only L1 to L4 are planted. The vector y now has 12 plots and X still has three blocks. Z still has six columns, but the columns for L5 and L6 are empty: not a single plot carries their data. The model is the same, with one change: g ~ N(0, Gσ²_g). Henderson's mixed-model equations show exactly where G enters:

The solution contains a value for all six lines. These are genomic estimated breeding values, GEBVs. L5 receives +0.54 t/ha and L6 −0.48 t/ha, although neither was planted. In the simulation their true genetic values were +0.54 and −0.52.

Figure 3d opens the box. The GEBV of an untested line is a weighted sum of the GEBVs of the tested lines, and the weights come from G. L5 takes 0.74 of the value of L1, its sib line, and almost nothing from the others: 0.74 × 0.73 + 0.08 × 0.10 + 0.05 × (−0.64) − 0.10 × (−0.31) = +0.54. L6 leans on its own sib line, L3, with a weight of 0.78. Under I every one of these weights is zero, and both lines are predicted at the population mean.

The whole story in one animated figure

Press play and follow a single allele from the grandparents to the cousins, watch Wright's paths add up to a relationship, see A fill in and turn into G, and then see the mixed model use G to predict a seedling that has never been planted. Use the I / G switch in the last step to see what happens when the relatives are taken away.

Open the animation full screen · Five steps, about 50 seconds. Relationships are from our simulation; the yields in step 5 are illustrative.

Route 1 and Route 2 by hand, in base R

No packages. We build X and Z, write Henderson's equations and solve them, first with I and then with G, using the data in Figure 3.

# Six inbred lines (L1-L6) in an RCBD with 3 blocks. Yield in t/ha.
y     <- c(5.7, 6.4, 6.3,   5.1, 4.5, 6.3,   3.5, 5.3, 4.0,    # L1, L2, L3
           4.8, 4.4, 4.9,   5.9, 5.9, 5.7,   3.7, 4.6, 4.9)    # L4, L5, L6
line  <- factor(rep(paste0("L", 1:6), each = 3))   # line of each plot
block <- factor(rep(1:3, times = 6))               # block of each plot
X <- model.matrix(~ block - 1)     # plots x blocks: the fixed effects
Z <- model.matrix(~ line - 1)      # plots x lines:  the random effects
lambda <- 1                        # sigma2_e / sigma2_g  (plot-level h2 = 0.5)

# ROUTE 1 - lines treated as unrelated: var(u) = I
I6  <- diag(6)
LHS <- rbind(cbind(t(X) %*% X, t(X) %*% Z),
             cbind(t(Z) %*% X, t(Z) %*% Z + solve(I6) * lambda))
RHS <- c(t(X) %*% y, t(Z) %*% y)
sol <- solve(LHS, RHS)             # Henderson's mixed-model equations
round(sol[1:3], 2)                 # BLUEs: the three block means
round(sol[4:9], 2)                 # BLUPs: line effects (deviations from the mean)

# ROUTE 2 - only L1-L4 grown; lines related through markers: var(g) = G
G <- matrix(c( 0.97,-0.12, 0.01, 0.17, 0.70, 0.02,
              -0.12, 0.95, 0.10,-0.18, 0.01, 0.00,
               0.01, 0.10, 0.95, 0.05, 0.06, 0.74,
               0.17,-0.18, 0.05, 1.00, 0.01,-0.01,
               0.70, 0.01, 0.06, 0.01, 0.96, 0.03,
               0.02, 0.00, 0.74,-0.01, 0.03, 0.90), nrow = 6)
G  <- G + diag(0.01, 6)            # tiny boost on the diagonal so G can be inverted
tr <- line %in% c("L1", "L2", "L3", "L4")   # plots of the four tested lines
Xt <- X[tr, ];  Zt <- Z[tr, ];  yt <- y[tr]
LHS2 <- rbind(cbind(t(Xt) %*% Xt, t(Xt) %*% Zt),
              cbind(t(Zt) %*% Xt, t(Zt) %*% Zt + solve(G) * lambda))
RHS2 <- c(t(Xt) %*% yt, t(Zt) %*% yt)
sol2 <- solve(LHS2, RHS2)
round(sol2[1:3], 2)                # BLUEs: block means from the tested plots
round(sol2[4:9], 2)                # GEBVs for all six lines, L5 and L6 included
# Now replace solve(G) with solve(I6) and run Route 2 again: L5 and L6 get exactly 0.

Three papers behind this edition

1. Wright S (1922). Coefficients of inbreeding and relationship. The American Naturalist 56: 330–338.The path rule in Figure 1c, still the cleanest way to read relationship from a pedigree.
2. Henderson CR (1975). Best linear unbiased estimation and prediction under a selection model. Biometrics 31: 423–447.BLUE and BLUP from one set of mixed-model equations: the engine of Figure 3.
3. VanRaden PM (2008). Efficient methods to compute genomic predictions. Journal of Dairy Science 91: 4414–4423.The G matrix in Figure 2b, still the default in most software.

Next time you see a genomic prediction, ask two questions: what is in its covariance, and how related is the training set to the lines being predicted? Together they explain most of what the model can and cannot do.

Your turn: in your programme, which lines are the least related to your training set? Tell us in the comments.

If this was useful, share it with one colleague and subscribe to The Crop Science Brief.

#CropScience #PlantBreeding #GenomicSelection #QuantitativeGenetics #PlantNura