The previous methods have focused on algorithms to generate the inverse of the numerator relationship matrix (A) in order to
generate breeding values in the mixed model equations. This section focuses on methods to generate A for a subset of individuals
without having to set up A across all of the individuals in a pedigree (i.e. outlined in the
Recursive Method to Create A section ). Generating subsets of A are useful when generating mating designs on
a given list of sire and dams or to generate the H matrix when generating breeding values with single-step genomic BLUP. The algorithm is
outlined in Colleau (2002), although I found it difficult to figure out how to write up the algorithm. I found manuscripts by Aguilar et al.
(2011) and Misztal et al. (2009) to be helpfull and are a great place to go after reading over the initial Colleau (2002) manuscript.
The method generates subsets of a relationship matrix across the desired individuals based on solving the following equation:
$$ w = Av $$
where w is a vector of the relationship between a desired individual and all of the animals in the pedigree, A is the full relationship
and v is a vector with a value of '1' for the desired individual and a '0' for all remaining animals in the pedigree. The first part of the algorithm
generates inbreeding values across all individuals prior to generating relationships for a subset of individuals. Once inbreeding estimates are generated,
a vector, q is initialized that loops through all individuals in the pedigree to calculate the fraction that derivies from its ancestors and is similar to
code from a previous section (i.e. Computation of Inbreeding Coeffecients ). Lastly,
the final steps involve calculation of D (i.e. square root of diagonals of L) and the numerator relationship values between an animals and the other
animals in the subset.
Similar to previous methods, the pedigree has to be sorted so that parents come before progeny. Lastly, if animals are numbered from 1 to the
very last animal then the sire and dam values can be used to index where the respective elements are located within the algorithm. Algorithms to sort
and renumber a pedigree is outlined in the following section . It is assumed that the inbreeding
coefficients are known across all animals in the pedigree prior to actually starting to compute relationships. The first part of the code below computes
the inbreeding coefficients using the method outlined in the Quaas (1976) section or the
Meuwissen & Luo (1992) section.
A_Colleau2002 <- function(animal,sire,dam,subsetids,inbreeding = "Meuwissen_Luo")
{
if(inbreeding != "Meuwissen_Luo" & inbreeding != "Quass")
{
stop("Inbreeding method is not an available option.")
}
if(is.vector(animal) == FALSE){stop("The animal input variable needs to be a vector!!")}
if(is.vector(sire) == FALSE){stop("The sire input variable needs to be a vector!!")}
if(is.vector(dam) == FALSE){stop("The dam input variable needs to be a vector!!")}
if(is.vector(subsetids) == FALSE){stop("The subsetids input variable needs to be a vector!!")}
##########################
## Calculate Inbreeding ##
##########################
if(inbreeding == "Quass")
{
## Look at Quass Method to Generate Ainv for Comments ##
m <- length(animal); u = numeric(m); v = numeric(m)
for(k in 1:m)
{
if(sire[k] == 0 & dam[k] = =0){v[k] <- 1;}
if(sire[k] == 0 & dam[k] != 0){v[k] = sqrt(1-0.25*u[dam[k]])}
if(sire[k] != 0 & dam[k] == 0){v[k] = sqrt(1-0.25*u[sire[k]])}
if(sire[k] != 0 & dam[k] != 0){v[k] = sqrt(1-0.25*u[sire[k]]-0.25*u[dam[k]])}
if(k < m)
{
for(i in (k+1):m)
{
if(sire[i] >= k & dam[i] >= k)
{
v[i] = 0.5*v[sire[i]] + 0.5*v[dam[i]]
}
if(sire[i] < k & dam[i] >= k) { v[i] = 0.5*v[dam[i]]}
if(dam[i] < k & sire[i] >= k) { v[i] = 0.5*v[sire[i]]}
if(sire[i] < k & dam[i] >= k) { v[i] = 0.5*v[dam[i]]}
if(sire[i] < k & dam[i] < k) { v[i] = 0}
}
}
u[k:m] = u[k:m] + v[k:m]^2
}
F <- u-1; rm(m,u,v)
}
if(inbreeding == "Meuwissen_Luo")
{
## Look at Meuwissen &Luo Method to Generate Ainv for Comments ##
m <- length(animal);
L_ml <- rep(0,(length(animal))); D_ml <- rep(0,length(animal));
F_ml <- rep(0,(length(animal)+1)); F_ml[1] = -1;
for(k in 1:m)
{
L_ml <- rep(0,length(animal));
AN <- c(k);
ai = 0.0;
L_ml[k] = 1.0;
D_ml[k] = 0.5 - (0.25* (F_ml[sire[k]+1] + F_ml[dam[k]+1]))
j = k;
while(length(AN) > 0)
{
if(sire[j] != 0)
{
AN <- c(AN,sire[j])
L_ml[sire[j]] <- L_ml[sire[j]] + (0.5 * L_ml[j]);
{
if(dam[j] != 0)
{
AN <- c(AN,dam[j])
L_ml[dam[j]] <- L_ml[dam[j]] + (0.5 * L_ml[j]);
}
ai = ai + (L_ml[j] * L_ml[j] * D_ml[j]);
found = 1;
repeat
{
if(AN[found] == j){AN <- AN[-c(found)]; break;}
if(AN[found] != j){found <- found + 1;}
}
if(length(AN) > 0){j = -1; j <- max(AN); AN <- unique(AN)}
}
F_ml[k+1] = ai - 1;
}
F <- F_ml[2:length(F_ml)]; rm(ai,AN,D_ml,F_ml,found,j,k,L_ml,m)
}
##########################################
## Generate Subset of Inbreeding Matrix ##
##########################################
nGen <- length(subsetids); # Number of animals for subset ##
n <- length(animal); # Size of full pedigree ##
A22 <- matrix(data=NA,nrow=nGen,ncol=nGen,byrow=TRUE)
w <- c(rep(0,n)) # Solution to Av ##
v <- c(rep(0,n)) # Vector multiplied by A ##
for(i in 1:nGen)
{
v <- c(rep(0,n))
v[subsetids[i]] <- 1.0
q <- rep(0,n) # row for animal in T #
for(j in n:1)
{
q[j] = q[j] + v[j]
s <- sire[j]; d <- dam[j]
if(s != 0){q[s] = q[s] + q[j] * 0.5}
if(d != 0){q[d] = q[d] + q[j] * 0.5}
}
for(j in 1:n)
{
s <- sire[j]; d <- dam[j]; di = 0;
if(s != 0 && d != 0){di <- (2 / 4) - (0.25 * (F[s] + F[d]))}
if(s != 0 && d == 0){di <- (3 / 4) - (0.25 * (F[s]))}
if(s == 0 && d != 0){di <- (3 / 4) - (0.25 * (F[d]))}
if(s == 0 && d == 0){di <- 1}
tmp = 0;
if(s != 0){tmp = tmp + w[s]}
if(d != 0){tmp = tmp + w[d]}
w[j] = 0.5 * tmp;
w[j] = w[j] + (di * q[j]);
}
A22[,i] = w[subsetids]
}
return(A22)
}
The full pedigree file from Henderson (1976) can be utilized to compute relationships between the last
three animals (5, 6 and 7) utilizing the R function above. The columns are animal, sire and dam and the pedigree is already ordered so that parents come
before progeny. Lastly the animals go from 1 to the total number of animals. Outlined below is what A for the three animals, v and w look like at
the end of each iteration of the for loop.
- Colleau, J.J. 2002. An indirect approach to the extensive calculation of relationship coefficients. Genet Sel Evol. 34(4):409-21.
- Misztal, I., A. Legarra & I. Aguilar. 2009. Computing procedures for genetic evaluation including phenotypic, full pedigree, and genomic information.
J Dairy Sci. 92(9):4648-55.
- Aguilar, I., I. Misztal, A. Legarra & S. 2011. Tsuruta. Efficient computation of the genomic relationship matrix and other matrices used in single-step
evaluation. J Anim Breed Genet. 128(6):422-8.