snp <- c("AA Ct GG Ag AA Ct Ga AA tt CC GG AA")
ss <- strsplit(snp,split = " ")
snpM <- matrix(ss[[1]],3,4,byrow = T)
row.names(snpM) <- c("ID1","ID2","ID3")
colnames(snpM) <- c("snp1","snp2","snp3","snp4")
# Genotype:AA Aa aa
snpM
library(synbreed)
gp <- create.gpData(geno = snpM)
gpb <- codeGeno(gp,label.heter = "alleleCoding", maf = 0.01, nmiss = 0.1,
impute = TRUE, impute.type = "random", verbose = TRUE)
agenotype <- gpb$geno
# Genotype 0 1 2
MAF <- agenotype
MAF
# Genotype -1 0 2
M <- MAF-1
M
#对角线为每个个体含有纯合的位点，比如ID1有2个纯合位点
M%*%t(M)
t(M)%*%M
p1 <- c(0.383,0.244,0.167,0.067)
names(p1) <- paste("snp",1:4,sep = "")
p1
p=2*(p1-0.5)
P=matrix(p,byrow=T,nrow=nrow(M),ncol=ncol(M))
P
Z=as.matrix(M-P)
Z
b=1-p1
c=p1*b
d=2*(sum(c))
ZZt=Z%*%t(Z)
G1=(ZZt/d)
G1 = round(G1,3)
G1
D=1/(ncol(M)*(2*p1*(1-p1)))
G2 = Z%*%(D*t(Z))
G2 =round(G2,3)
G2
p1=array(0.5,ncol(M))
p=p1
P=matrix(p,byrow=T,nrow=nrow(M),ncol=ncol(M))
Z=as.matrix(M-P)
b=1-p1
c=p1*b
d=2*(sum(c))
ZZt=Z%*%t(Z)
G3=(ZZt/d)
G3
library(GeneticsPed)
data(Mrode3.1)
x <- Pedigree(x=Mrode3.1, subject="calf", ascendant=c("sire", "dam"),ascendantSex=c("Male", "Female"), sex="sex")
x
## Observed/measured phenotype records:
y <- x$pwg
y
## Design matrix(X) for sex effect:
X <- model.matrix(~x$sex -1)
X
t(X)
## Design matrix (Z) for additive genetic effect. Note that first three columns do not have indicators since
## these columns are for individuals without phenotype records and apear in the model only through the
## pedigree
Z <- model.matrix(object = x, y=x$pwg, id=x$calf)
Z
## Left hand side (LHS) of MME without G −1 :
LHS <- rbind(cbind(t(X) %*% X, t(X) %*% Z),cbind(t(Z) %*% X, t(Z) %*% Z))
## or more efficiently
#LHS <- rbind(cbind(crossprod(X), crossprod(X, Z)),+ cbind(crossprod(Z, X), crossprod(Z)))
LHS
## We want Ainv for all individuals in the pedigree not only individuals with records
x <- extend(x)
x
Ainv <- inverseAdditive(x=x)
Ainv
sigma2a <- 20
sigma2e <- 40
alpha <- sigma2e / sigma2a
q <- nIndividual(x)
q
p <- nrow(LHS) - q
p
#LHS[(p + 1):(p + q), (p + 1):(p + q)] <- LHS[(p + 1):(p + q), (p + 1):(p + q)] + Ainv * alpha
LHS[(p + 1):(p + q), (p + 1):(p + q)]
LHS[(p + 1):(p + q), (p + 1):(p + q)] <- LHS[(p + 1):(p + q), (p + 1):(p + q)] + Ainv * alpha
LHS
## Right hand side (RHS) of MME:
RHS <- rbind(t(X) %*% y, t(Z) %*% y)
## or more efficiently
#RHS <- rbind(crossprod(X, y),crossprod(Z, y))
t(RHS)
## Solution
sol <- solve(LHS) %*% RHS
## or more efficiently
sol <- solve(LHS, RHS)
sol
###################################################################################
library(GeneticsPed)
data(Mrode3.1)
x <- Pedigree(x=Mrode3.1, subject="calf", ascendant=c("sire", "dam"),ascendantSex=c("Male", "Female"), sex="sex")
x
## Observed/measured phenotype records:
y <- x$pwg
y
## Design matrix(X) for sex effect:
X <- model.matrix(~x$sex -1)
X
t(X)
## Design matrix (Z) for additive genetic effect. Note that first three columns do not have indicators since
## these columns are for individuals without phenotype records and apear in the model only through the
## pedigree
Z <- model.matrix(object = x, y=x$pwg, id=x$calf)
Z
## Left hand side (LHS) of MME without G −1 :
tzz <- diag(0.001, nrow = nrow(t(Z) %*% Z), ncol = ncol(t(Z) %*% Z)) + t(Z) %*% Z
tzz
LHS <- rbind(cbind(t(X) %*% X, t(X) %*% Z),cbind(t(Z) %*% X, tzz))
## or more efficiently
#LHS <- rbind(cbind(crossprod(X), crossprod(X, Z)),+ cbind(crossprod(Z, X), crossprod(Z)))
LHS
## Right hand side (RHS) of MME:
RHS <- rbind(t(X) %*% y, t(Z) %*% y)
## or more efficiently
#RHS <- rbind(crossprod(X, y),crossprod(Z, y))
t(RHS)
## Solution
sol <- solve(LHS) %*% RHS
## or more efficiently
#sol <- solve(LHS, RHS)
sol
library(lme4)
library(dplyr)
library(tidyr)
x
## We want Ainv for all individuals in the pedigree not only individuals with records
x <- extend(x)
x
xx <- x %>% gather(3:4, key = parents, value = val) %>% arrange(calf, sex) %>% filter(!is.na(pwg))
xx
str(xx)
temp <- lmer(pwg ~ 0 + sex + (1 | calf), xx)
print(temp)
anova(temp) ####求方差
ranef(temp) ####求随机效应的BLUP值
fixef(temp) ####求固定效应的BLUE值
model.matrix(temp)
## X and Z are the fixed- and random-effects model matrices respectively
X1 = getME(temp, "X")
X1
Z1 = getME(temp, "Z")
Z1
LHS <- rbind(cbind(t(X1) %*% X1, t(X1) %*% Z1),cbind(t(Z1) %*% X1, t(Z1) %*% Z1))
LHS
## Left hand side (LHS) of MME without G −1 :
tzz <- diag(0.001, nrow = nrow(t(Z1) %*% Z1), ncol = ncol(t(Z1) %*% Z1)) + t(Z1) %*% Z1
tzz
LHS <- rbind(cbind(t(X1) %*% X1, t(X1) %*% Z1),cbind(t(Z1) %*% X1, tzz))
LHS
## Observed/measured phenotype records:
y <- xx$pwg
y
## Right hand side (RHS) of MME:
RHS <- rbind(t(X1) %*% y, t(Z1) %*% y)
## or more efficiently
#RHS <- rbind(crossprod(X, y),crossprod(Z, y))
t(RHS)
## Solution
sol <- solve(LHS) %*% RHS
## or more efficiently
#sol <- solve(LHS, RHS)
sol
VarCorr(temp)
##############################################################################
summary(temp)