# NOT RUN {
####=========================================####
#### For CRAN time limitations most lines in the
#### examples are silenced with one '#' mark,
#### remove them and run the examples
####=========================================####
####=========================================####
####=========================================####
#### EXAMPLE 1
#### breeding values with 1 variance component
####=========================================####
####=========================================####
####=========================================####
#### simulate genotypic data
#### random population of 200 lines with 1000 markers
####=========================================####
M <- matrix(rep(0,200*1000),1000,200)
for (i in 1:200) {
M[,i] <- ifelse(runif(1000)<0.5,-1,1)
}
colnames(M) <- paste("geno",1:dim(M)[2], sep="")
rownames(M) <- paste("mark",1:dim(M)[1], sep="")
####=========================================####
#### simulate phenotypes
####=========================================####
QTL <- 100*(1:5) #pick 5 QTL
u <- rep(0,1000) #marker effects
u[QTL] <- 1
g <- as.vector(crossprod(M,u))
h2 <- 0.5
y <- g + rnorm(200,mean=0,sd=sqrt((1-h2)/h2*var(g)))
M <- t(M); M[1:5,1:5]
####=========================================####
#### fit the model using mmer (matrix based)
####=========================================####
Z1 <- diag(length(y))
K1 <-A.mat(M)
ETA <- list( add=list(Z=Z1, K=K1) )
#ans <- mmer(Y=y, Z=ETA)
#summary(ans)
#### run the same but as GWAS
#### just add the marker matrix in the argument W
#### markers are fixed effects
#ans <- mmer(Y=y, Z=ETA, W=M)
#summary(ans)
####=========================================####
#### fit the model using mmer2 (formula-based)
####=========================================####
K<-A.mat(M) # additive relationship matrix
dat <- data.frame(y=y,id=rownames(M));head(dat)
#ans <- mmer2(y~1, random=~g(id), G=list(id=K), data=dat)
#summary(ans)
#### run the same but as GWAS
#### just add the marker matrix in the argument W
#### markers are fixed effects
#ans <- mmer2(y~1, random=~g(id), G=list(id=K), W=M, data=dat)
#summary(ans)
############################################
############################################
############################################
############################################
############################################
############################################
############################################
############################################
####=========================================####
####=========================================####
#### EXAMPLE 2
#### breeding values with 3 variance components
####=========================================####
####=========================================####
####=========================================####
#### Hybrid prediction using mmer (matrix-based)
####=========================================####
data(cornHybrid)
hybrid2 <- cornHybrid$hybrid # extract cross data
A <- cornHybrid$K
y <- hybrid2$Yield
X1 <- model.matrix(~ Location, data = hybrid2);dim(X1)
Z1 <- model.matrix(~ GCA1 -1, data = hybrid2);dim(Z1)
Z2 <- model.matrix(~ GCA2 -1, data = hybrid2);dim(Z2)
Z3 <- model.matrix(~ SCA -1, data = hybrid2);dim(Z3)
colnames(Z1) <- levels(hybrid2$GCA1)
colnames(Z2) <- levels(hybrid2$GCA2)
colnames(Z3) <- levels(hybrid2$SCA)
####=========================================####
#### Realized IBS relationships for set of parents 1
####=========================================####
#K1 <- A[levels(hybrid2$GCA1), levels(hybrid2$GCA1)]; dim(K1)
####=========================================####
#### Realized IBS relationships for set of parents 2
####=========================================####
#K2 <- A[levels(hybrid2$GCA2), levels(hybrid2$GCA2)]; dim(K2)
####=========================================####
#### Realized IBS relationships for cross
#### (as the Kronecker product of K1 and K2)
####=========================================####
#S <- kronecker(K1, K2) ; dim(S)
#rownames(S) <- colnames(S) <- levels(hybrid2$SCA)
#ETA <- list(GCA1=list(Z=Z1, K=K1),
# GCA2=list(Z=Z2, K=K2),
# SCA=list(Z=Z3, K=S)
# )
#ans <- mmer(Y=y, X=X1, Z=ETA)
#ans$var.comp
#summary(ans)
####=========================================####
#### Hybrid prediction using mmer2 (formula-based)
####=========================================####
# data(cornHybrid)
# hybrid2 <- cornHybrid$hybrid # extract cross data
# A <- cornHybrid$K
# ####=========================================####
# #### Realized IBS relationships for set of parents 1
# ####=========================================####
# K1 <- A[levels(hybrid2$GCA1), levels(hybrid2$GCA1)]; dim(K1)
# ####=========================================####
# #### Realized IBS relationships for set of parents 2
# ####=========================================####
# K2 <- A[levels(hybrid2$GCA2), levels(hybrid2$GCA2)]; dim(K2)
# ####=========================================####
# #### Realized IBS relationships for cross
# #### (as the Kronecker product of K1 and K2)
# ####=========================================####
# S <- kronecker(K1, K2) ; dim(S)
# rownames(S) <- colnames(S) <- levels(hybrid2$SCA)
#
# ans <- mmer2(Yield~Location, random=~g(GCA1)+ g(GCA2) + g(SCA),
# G=list(GCA1=K1, GCA2=K2, SCA=S), data=hybrid2)
# summary(ans)
############################################
############################################
############################################
############################################
############################################
############################################
############################################
############################################
####=========================================####
####=========================================####
#### EXAMPLE 3
#### PREDICTING SPECIFIC PERFORMANCE
#### within biparental population
####=========================================####
####=========================================####
####=========================================####
#### using mmer (matrix-based)
####=========================================####
data(CPdata)
CPpheno <- CPdata$pheno[,-c(1:4)]
CPgeno <- CPdata$geno
#### look at the data
head(CPpheno)
CPgeno[1:5,1:5]
#### fit a model including additive and dominance effects
y <- CPpheno$color
Za <- diag(length(y))
Zd <- diag(length(y))
Ze <- diag(length(y))
#A <- A.mat(CPgeno) # additive relationship matrix
#D <- D.mat(CPgeno) # dominant relationship matrix
####================####
#### ADDITIVE MODEL ####
####================####
#ETA.A <- list(add=list(Z=Za,K=A))
#ans.A <- mmer(Y=y, Z=ETA.A)
#summary(ans.A)
####=========================####
#### ADDITIVE-DOMINANCE MODEL ####
####=========================####
#ETA.AD <- list(add=list(Z=Za,K=A), dom=list(Z=Zd,K=D))
#ans.AD <- mmer(Y=y, Z=ETA.AD)
#summary(ans.AD)
### greater accuracy !!!! 4 percent increment!!
### we run 100 iterations, 4 percent increment in general
####=========================================####
#### same using mmer2 (formula-based)
####=========================================####
# data(CPdata)
# CPpheno <- CPdata$pheno
# CPgeno <- CPdata$geno
# ### look at the data
# head(CPpheno); CPgeno[1:5,1:5]
# ###=========================================####
# ### fit a model including additive and dominance effects
# ###=========================================####
# A <- A.mat(CPgeno)
# D <- D.mat(CPgeno)
# ####=========================================####
# #### ADDITIVE MODEL
# ####=========================================####
# ans.A <- mmer2(Yield~1,random=~ g(id), G=list(id=A),
# data=CPpheno)
# summary(ans.A)
# ####=========================================####
# #### ADDITIVE-DOMINANCE MODEL
# ####=========================================####
# CPpheno$idd <- CPpheno$id
# ans.AD <- mmer2(Yield~1,random=~ g(id) + g(idd),
# G=list(id=A, idd=D), data=CPpheno)
# summary(ans.AD)
############################################
############################################
############################################
############################################
############################################
############################################
############################################
############################################
####=========================================####
####=========================================####
#### EXAMPLE 4
#### PARALLEL RESPONSES
#### accelerating analysis
#### using 'n.cores' argument for parallelization
####=========================================####
####=========================================####
# data(CPdata)
# CPpheno <- CPdata$pheno; CPpheno$idd <- CPpheno$id
# CPgeno <- CPdata$geno
# ###=========================================####
# ### Do incidence and relationship matrices
# ###=========================================####
# A <- A.mat(CPgeno)
# D <- D.mat(CPgeno)
# ###=========================================####
# ### Only additive model using all traits
# ###=========================================####
# head(CPpheno)
# ans.A <- mmer2(cbind(color,Yield,FruitAver)~1, random = ~ g(id),
# G=list(id=A), n.cores = 3, data=CPpheno)
# summary(ans.A)
# ###=========================================####
# ### Additive + dominance model
# ###=========================================####
# head(CPpheno)
# ans.AD <- mmer2(cbind(color,Yield,FruitAver)~1, random=~g(id)+g(idd),
# G=list(id=A, idd=D), n.cores = 3, data=CPpheno)
# summary(ans.AD)
############################################
############################################
############################################
############################################
############################################
############################################
####=========================================####
####=========================================####
#### EXAMPLE 5
#### MULTIVARIATE MODEL
####=========================================####
####=========================================####
# data(CPdata)
# CPpheno <- CPdata$pheno
# CPgeno <- CPdata$geno
# ### look at the data
# head(CPpheno)
# CPgeno[1:5,1:5]
# ## fit a model including additive effects
# A <- A.mat(CPgeno)
# ### MAKE SURE YOU SET 'MVM=TRUE'
# ans.A <- mmer2(cbind(color,Yield,FruitAver)~1, random = ~ g(id),
# G=list(id=A),MVM=TRUE, data=CPpheno)
# summary(ans.A)
############################################
############################################
############################################
############################################
############################################
############################################
####=========================================####
####=========================================####
#### EXAMPLE 6
#### OVERLAY IN MODELS
####=========================================####
####=========================================####
####=========================================####
#### using mmer (matrix-based)
####=========================================####
data(HDdata)
head(HDdata)
#### GCA matrix for half diallel using male and female columns
Z1 <- overlay(HDdata[,c("male","female")])
#### SCA matrix
Z2 <- model.matrix(~as.factor(geno)-1, data=HDdata)
#### Define response variable and run
y <- HDdata$sugar
ETA <- list(GCA=list(Z=Z1), SCA=list(Z=Z2)) # Zu component
modHD <- mmer(Y=y, Z=ETA)
summary(modHD)
####=========================================####
#### using overlay with mmer2 function (formula-based)
####=========================================####
# data(HDdata)
# head(HDdata)
# HDdata$female <- as.factor(HDdata$female)
# HDdata$male <- as.factor(HDdata$male)
# HDdata$geno <- as.factor(HDdata$geno)
# #### model using overlay
# modh <- mmer2(sugar~1, random=~female + and(male) + geno,
# data=HDdata)
# summary(modh)
# #### model using overlay [and(.)] and covariance structures [g(.)]
# A <- diag(7); A[1,2] <- 0.5; A[2,1] <- 0.5 # fake covariance structure
# colnames(A) <- as.character(1:7); rownames(A) <- colnames(A);A
#
# modh2 <- mmer2(sugar~1, random=~g(female) + and(g(male)) + geno,
# G=list(female=A, male=A),data=HDdata)
# summary(modh2)
# data(HDdata)
# head(HDdata)
# HDdata$female <- as.factor(HDdata$female)
# HDdata$male <- as.factor(HDdata$male)
# HDdata$geno <- as.factor(HDdata$geno)
# #### model using overlay
# modh <- mmer2(sugar~1, random=~female + and(male) + geno,
# data=HDdata)
# summary(modh)
# }
Run the code above in your browser using DataLab