# R/binQTLScan.R In venyao/binQTL: QTL mapping with binmap

```binQTLScan <- function(phenotype="", genotype="", population="RIL") {
if (population == "MAGIC") {
result <- magicbinQTL(phenotype=phenotype, genotype=genotype)

return(result)
} else {
Z <- as.matrix(genotype[, -c(1:4)])
Y <- as.matrix(phenotype[, 2])
genotype.n <- genoConv(genotype)
kk <- kinship(genotype.n\$genotype.a)
m <- nrow(Z)
n <- ncol(Z)
X <- as.matrix(rep(1, n))

qq <- eigen(kk, symmetric = T)
uu <- qq\$vectors
delta <- qq\$values

r <- ncol(X)
x <- t(uu) %*% X
y <- t(uu) %*% Y
if (population=="RIL") {
L <- matrix(c(-0.5, 0.5), 1, 2)
} else {
L <- rbind(c(-0.5, 0, 0.5), c(-0.5, 1, -0.5))
}
invLL <- solve(L %*% t(L))
rnk <- nrow(L)

triplet <- TRUE
parr <- numeric()
for (k in 1:m) {
if ((triplet == FALSE) | (k == 1) | (k == m)) {
U <- as.factor(Z[k, ])
Zk <- model.matrix( ~ U - 1)
zk <- t(uu) %*% as.matrix(Zk)
fit <- random( x = x, z = zk, y = y, m = m, kk = qq, cov = "qq", delta = delta, n = n)
} else{
U1 <- as.factor(Z[k - 1, ])
U2 <- as.factor(Z[k, ])
U3 <- as.factor(Z[k + 1, ])
Z.1 <- model.matrix( ~ U1 - 1)
Z.2 <- model.matrix( ~ U2 - 1)
Z.3 <- model.matrix( ~ U3 - 1)
Zk <- cbind(Z.1, Z.2, Z.3)
zk <- t(uu) %*% as.matrix(Zk)
fit <- random1( x = x, z = zk, y = y, m = m, kk = qq, cov = "qq", delta = delta, n = n )
}
covparm <- fit[[1]]
vg <- covparm[1]
vk <- covparm[2]
ve <- covparm[3]
lambda <- covparm[4]
lambda.k <- covparm[5]
beta <- fit[[2]]
gamma <- fit[[3]]
stderr <- fit[[4]]
cov <- fit[[5]]
lrt <- fit[[6]]
effect <- L %*% gamma
vv <- L %*% cov %*% t(L)
wald <- t(effect) %*% solve(vv) %*% effect
dk <- rnk - sum(diag(invLL %*% vv)) / vk
if (dk < 1e-5) { pwald <- 1 } else { pwald <- 1 - pchisq(wald, dk) }
if (lrt < 1e-8) { plrt <- 1 } else { plrt <- 0.5 * (1 - pchisq(lrt, 1)) }
par <- c( k, vg, vk, ve, lambda, lambda.k, drop(beta), drop(gamma), drop(stderr),
drop(effect), wald, dk, pwald, -log10(pwald), lrt, plrt, -log10(plrt)  )
parr <- rbind(parr, par)
}

if (population=="RIL") {
varnames <- c("k", "vg", "vk", "ve", "lambda", "lambdak", "beta", "g11",
"g22", "sd11", "sd22")
varnames <- c(varnames, "additive", "wald", "dk", "pwald", "logpwald",
"lrt", "plrt", "logplrt")
} else if (population=="F2") {
varnames <- c("k", "vg", "vk", "ve", "lambda", "lambdak", "beta", "g11",
"g12", "g22", "sd11", "sd12", "sd22")
varnames <- c(varnames, "additive", "dominance", "wald", "dk", "pwald",
"logpwald", "lrt", "plrt", "logplrt")
}

out <- as.matrix(parr)
colnames(out) <- varnames
rownames(out) <- 1:nrow(out)
result <- data.frame(genotype[, 1:4], out, stringsAsFactors=FALSE)
result <- result[, c("Bin", "Chr", "Start", "Stop", "plrt", "logplrt")]
names(result)[5:6] <- c("p", "-logp")

return(result)
}
}
```
venyao/binQTL documentation built on June 9, 2019, 4:52 a.m.