# NOT RUN {
####=========================================####
#### For CRAN time limitations most lines in the
#### examples are silenced with one '#' mark,
#### remove them and run the examples using
#### command + shift + C |OR| control + shift + C
####=========================================####
data(wheatLines)
X <- wheatLines$wheatGeno; X[1:5,1:5]; dim(X)
Y <- wheatLines$wheatPheno
rownames(X) <- rownames(Y)
####=========================================####
#### select environment 1
####=========================================####
#y <- Y[,1] # response grain yield
#Z1 <- diag(length(y)) # incidence matrix
#K <- A.mat(X) # additive relationship matrix
####=========================================####
#### GBLUP pedigree-based approach
####=========================================####
#ETA <- list(line=list(Z=Z1, K=K))
#ans <- mmer(Y=y, Z=ETA, method="EMMA") # kinship based
#summary(ans)
####=========================================####
#### GBLUP marker based approach
####=========================================####
#ETA2 <- list(mar=list(Z=X))
#ans2 <- mmer(Y=y, Z=ETA2, method="EMMA") # marker based
#summary(ans2)
####=========================================####
#### compare and check that is the same result
####=========================================####
#plot(ans$u.hat$line, (X%*%ans2$u.hat$mar), xaxt="n", yaxt="n",
# ylab="Marker-based GBLUP", xlab="Pedigree-based GBLUP")
####=========================================####
#### PREDICT PROGENY
####=========================================####
#GEBV.pb <- ans$u.hat$line # this are the BV
#rownames(GEBV.pb) <- rownames(Y)
####=========================================####
#### all possible crosses = 179,101
####=========================================####
#crosses <- do.call(expand.grid, list(rownames(Y),rownames(Y))); dim(crosses)
#cross2 <- duplicated(t(apply(crosses, 1, sort)))
#crosses2 <- crosses[cross2,]; head(crosses2); dim(crosses2)
####=========================================####
#### match the possible crosses with the parental BV
####=========================================####
#GCA1 = GEBV.pb[match(crosses2[,1], rownames(GEBV.pb))] # get GCA1 BLUP of each hybrid
#GCA2 = GEBV.pb[match(crosses2[,2], rownames(GEBV.pb))] # get GCA1 BLUP of each hybrid
####=========================================####
#### join everything
####=========================================####
#BV <- data.frame(crosses2,GCA1,GCA2); head(BV)
#BV$BVcross <- apply(BV[,c(3:4)],1,mean); head(BV)
#plot(BV$BVcross)
# }
Run the code above in your browser using DataLab