##############################################################
# R code for "A brief tour of R/qtl"
#
# Karl W Broman
# Department of Biostatistics, Johns Hopkins University
#
# http://www.rqtl.org
#
# 24 August 2007
##############################################################

save.image()

library(qtl)

ls()

help(read.cross)
?read.cross

############################################################
# Example 1: Hypertension
############################################################
data(hyper)
ls()
?hyper

summary(hyper)

nind(hyper)
nphe(hyper)
nchr(hyper)
totmar(hyper)
nmar(hyper)

plot(hyper)

plot.missing(hyper)
plot.map(hyper)
plot.pheno(hyper, pheno.col=1)

plot.map(hyper, chr=c(1, 4, 6, 7, 15), show.marker.names=TRUE)

plot.missing(hyper, reorder=TRUE)

hyper <- drop.nullmarkers(hyper)
totmar(hyper)

hyper <- est.rf(hyper)
plot.rf(hyper)
plot.rf(hyper, chr=c(1,4))

plot.rf(hyper, chr=6)
plot.missing(hyper, chr=6)

newmap <- est.map(hyper, error.prob=0.01)
plot.map(hyper, newmap)

hyper <- replace.map(hyper, newmap)

hyper <- calc.errorlod(hyper, error.prob=0.01)

top.errorlod(hyper)

plot.geno(hyper, chr=16, ind=c(24:34, 71:81))

plot.info(hyper)
plot.info(hyper, chr=c(1,4,15))
plot.info(hyper, chr=c(1,4,15), method="entropy")
plot.info(hyper, chr=c(1,4,15), method="variance")

hyper <- calc.genoprob(hyper, step=1, error.prob=0.01)

out.em <- scanone(hyper)
out.hk <- scanone(hyper, method="hk")

hyper <- sim.geno(hyper, step=2, n.draws=16, error.prob=0.01)
out.imp <- scanone(hyper, method="imp")

summary(out.em)
summary(out.em, threshold=3)
summary(out.hk, threshold=3)
summary(out.imp, threshold=3)

max(out.em)
max(out.hk)
max(out.imp)

plot(out.em, chr=c(1,4,15))
plot(out.em, out.hk, out.imp, chr=c(1,4,15))
plot(out.em, chr=c(1,4,15))
plot(out.hk, chr=c(1,4,15), col="blue", add=TRUE)
plot(out.imp, chr=c(1,4,15), col="red", add=TRUE)

operm.hk <- scanone(hyper, method="hk", n.perm=1000)

summary(operm.hk, alpha=0.05)

summary(out.hk, perms=operm.hk, alpha=0.05, pvalues=TRUE)

save.image()

hyper <- calc.genoprob(hyper, step=5, error.prob=0.01)

out2.hk <- scantwo(hyper, method="hk")

summary(out2.hk, thresholds=c(6.0, 4.7, 4.4, 4.7, 2.6))

summary(out2.hk, thresholds=c(6.0, 4.7, Inf, 4.7, 2.6))

plot(out2.hk)
plot(out2.hk, chr=c(1,4,6,15))

max(out2.hk)

operm2.hk <- scantwo(hyper, method="hk", n.perm=100)

summary(operm2.hk)

summary(out2.hk, perms=operm2.hk, pvalues=TRUE,
        alphas=c(0.05, 0.05, 0, 0.05, 0.05))

chr <- c(1, 1, 4, 6, 15)
pos <- c(50, 76, 30, 70, 20)
qtl <- makeqtl(hyper, chr, pos)

my.formula <- y ~ Q1 + Q2 + Q3 + Q4 + Q5 + Q4:Q5
out.fitqtl <- fitqtl(hyper$pheno[,1], qtl, formula=my.formula)
summary(out.fitqtl)

ls()
rm(list=ls())

############################################################
# Example 2: Genetic mapping
############################################################
data(badorder)
summary(badorder)
plot(badorder)

badorder <- est.rf(badorder)
plot.rf(badorder)

plot.rf(badorder, chr=1)

newmap <- est.map(badorder, verbose=TRUE)
plot.map(badorder, newmap)

plot.rf(badorder, chr=2:3)

pull.map(badorder, chr=2)
pull.map(badorder, chr=3)

badorder <- movemarker(badorder, "D2M937", 3, 48)
badorder <- movemarker(badorder, "D3M160", 2, 28.8)

plot.rf(badorder, chr=2:3)

rip1 <- ripple(badorder, chr=1, window=6)
summary(rip1)

rip2 <- ripple(badorder, chr=1, window=3, err=0.01, method="likelihood")
summary(rip2)

badorder.rev <- switch.order(badorder, 1, rip1[2,])
rip1r <- ripple(badorder.rev, chr=1, window=6)
summary(rip1r)

badorder.rev <- switch.order(badorder.rev, 1, rip1r[2,])
rip2r <- ripple(badorder.rev, chr=1, window=3, err=0.01)
summary(rip2r)

badorder.rev <- est.rf(badorder.rev)
plot.rf(badorder.rev, 1)

############################################################
# Example 3: Listeria susceptibility
############################################################
data(listeria)
summary(listeria)
plot(listeria)
plot.missing(listeria)

listeria$pheno$logSurv <- log(listeria$pheno[,1])
plot(listeria)

listeria <- est.rf(listeria)
plot.rf(listeria)
plot.rf(listeria, chr=c(5,13))

newmap <- est.map(listeria, error.prob=0.01)
plot.map(listeria, newmap)
listeria <- replace.map(listeria, newmap)

listeria <- calc.errorlod(listeria, error.prob=0.01)
top.errorlod(listeria)
top.errorlod(listeria, cutoff=3.5)
plot.geno(listeria, chr=13, ind=61:70, cutoff=3.5)

listeria <- calc.genoprob(listeria, step=2)
out.2p <- scanone(listeria, pheno.col=3, model="2part", upper=TRUE)

summary(out.2p)
summary(out.2p, threshold=4.5)

summary(out.2p, format="allpeaks", threshold=3)
summary(out.2p, format="allpeaks", threshold=c(4.5,3,3))

plot(out.2p)
plot(out.2p, lodcolumn=2)
plot(out.2p, lodcolumn=1:3, chr=c(1,5,13,15))

operm.2p <- scanone(listeria, model="2part", pheno.col=3,
                    upper=TRUE, n.perm=25)
summary(operm.2p, alpha=0.05)

summary(out.2p, format="allpeaks", perms=operm.2p,
        alpha=0.05, pvalues=TRUE)

y <- listeria$pheno$logSurv
my <- max(y, na.rm=TRUE)
z <- as.numeric(y==my)
y[y==my] <- NA
listeria$pheno$logSurv2 <- y
listeria$pheno$binary <- z
plot(listeria)

out.mu <- scanone(listeria, pheno.col=4)
plot(out.mu, out.2p, lodcolumn=c(1,3), chr=c(1,5,13,15), col=c("blue","red"))

out.p <- scanone(listeria, pheno.col=5, model="binary")
plot(out.p, out.2p, lodcolumn=c(1,2), chr=c(1,5,13,15), col=c("blue","red"))

out.np1 <- scanone(listeria, model="np", ties.random=TRUE)
out.np2 <- scanone(listeria, model="np", ties.random=FALSE)

plot(out.np1, out.np2, col=c("blue","red"))
plot(out.2p, out.np1, out.np2, chr=c(1,5,13,15))

############################################################
# Example 4: Covariates in QTL mapping
############################################################
data(fake.bc)
summary(fake.bc)
plot(fake.bc)

fake.bc <- calc.genoprob(fake.bc, step=2.5)
out.nocovar <- scanone(fake.bc, pheno.col=1:2)

sex <- fake.bc$pheno$sex
out.acovar <- scanone(fake.bc, pheno.col=1:2, addcovar=sex)

summary(out.nocovar, threshold=3, format="allpeaks")
summary(out.acovar, threshold=3, format="allpeaks")

plot(out.nocovar, out.acovar, chr=c(2, 5))
plot(out.nocovar, out.acovar, chr=c(2, 5), lodcolumn=2)

out.icovar <- scanone(fake.bc, pheno.col=1:2, addcovar=sex, intcovar=sex)

summary(out.icovar, threshold=3, format="allpeaks")

plot(out.acovar, out.icovar, chr=c(2,5), col=c("blue", "red"))
plot(out.acovar, out.icovar, chr=c(2,5), lodcolumn=2,
     col=c("blue", "red"))

out.sexint <- out.icovar - out.acovar
plot(out.sexint, lodcolumn=1:2, chr=c(2,5), col=c("green", "purple"))

seed <- ceiling(runif(1, 0, 10^8))
set.seed(seed)
operm.acovar <- scanone(fake.bc, pheno.col=1:2, addcovar=sex,
                        method="hk", n.perm=100)
set.seed(seed)
operm.icovar <- scanone(fake.bc, pheno.col=1:2, addcovar=sex,
                        intcovar=sex, method="hk", n.perm=100)

operm.sexint <- operm.icovar - operm.acovar

summary(operm.sexint, alpha=c(0.05, 0.20))

summary(out.sexint, perms=operm.sexint, alpha=0.1,
        format="allpeaks", pvalues=TRUE)

############################################################
# Example 5: Multiple QTL mapping
############################################################
rm(list=ls())
data(hyper)

source("http://www.rqtl.org/multqtlfunc.R")

hyper <- sim.geno(hyper, step=2.5, n.draws=16, err=0.01)

out1 <- scanone(hyper, method="imp")
plot(out1)

max(out1)

find.marker(hyper, 4, 29.5)

g <- pull.geno(hyper)[,"D4Mit164"]
mean(is.na(g))

g <- pull.geno(fill.geno(hyper))[,"D4Mit164"]

out1.c4 <- scanone(hyper, method="imp", addcovar=g)

plot(out1, out1.c4, col=c("blue", "red"))

plot(out1.c4 - out1, ylim=c(-3,3))
abline(h=0, lty=2, col="gray")

out1.c4i <- scanone(hyper, method="imp", addcovar=g, intcovar=g)

plot(out1.c4i - out1.c4)

out2 <- scantwo(hyper, method="imp")

summary(out2, thr=c(6.0, 4.7, Inf, 4.7, 2.6))

summary( subset(out2, chr=1) )

summary( subset(out2, chr=c(7,15)) )

plot(out2, chr=c(1,4,6,7,15))

plot(out2, chr=1, lower="cond-add")
plot(out2, chr=c(6,15), lower="cond-int")
plot(out2, chr=c(7,15), lower="cond-int")

out2.c4 <- scantwo(hyper, method="imp", addcovar=g, chr=c(1,6,7,15))

summary(out2.c4, thr=c(6.0, 4.7, Inf, 4.7, 2.6))
summary( subset(out2.c4, chr=1) )
summary( subset(out2.c4, chr=c(7,15)) )

plot(out2.c4, chr=c(1,4,6,7,15))
plot(out2.c4, chr=1, lower="cond-int")
plot(out2.c4, chr=c(6,15), lower="cond-int")
plot(out2.c4, chr=c(7,15), lower="cond-int")

out2sub <- subset(out2, chr=c(1,6,7,15))
plot(out2.c4 - out2sub, allow.neg=TRUE, lower="cond-int")

qc <- c(1, 1, 4, 6, 15)
qp <- c(43.3, 78.3, 30.0, 62.5, 18.0)
qtl <- makeqtl(hyper, chr=qc, pos=qp)
phe <- hyper$pheno[,1]

myformula <- y ~ Q1+Q2+Q3+Q4+Q5 + Q4:Q5

out.fq <- fitqtl(phe, qtl, formula = myformula)
summary(out.fq)

out.fq <- fitqtl(phe, qtl, formula = myformula, drop=FALSE, get.ests=TRUE)
summary(out.fq)

out.rq <- refineqtl(hyper, chr=qc, pos=qp, formula = myformula)

qp - out.rq[,2]

qp2 <- out.rq[,2]
qtl2 <- makeqtl(hyper, chr=qc, pos=qp2)
out.fq2 <- fitqtl(phe, qtl2, formula=myformula)
summary(out.fq2)

out1.sq <- scanqtl(hyper, chr=c(1,4), pos = list( c(-Inf,Inf), 29.5) )

null <- scanqtl(hyper, chr=4, pos=list(29.5))
out1.c4r <- convert.scanqtl(out1.sq, null)

plot(out1.c4, out1.c4r, col=c("blue", "red"), chr=1)

out2.sq.add <- scanqtl(hyper, chr=c(1,1,4),
                       pos=list(c(-Inf,Inf), c(-Inf,Inf), 29.5))
out2.sq.full <- scanqtl(hyper, chr=c(1,1,4),
                        pos=list(c(-Inf,Inf), c(-Inf,Inf), 29.5),
                        formula=y~Q1+Q2+Q3+Q1:Q2)

out2.c4r <- convert.scanqtl(out2.sq.full, null, out2.sq.add)

out2.c4sub <- subset(out2.c4, chr=1)
plot(out2.c4sub - out2.c4r, lower="cond-add", allow.neg=TRUE)

newpos <- c( as.list(qp2), list(c(-Inf, Inf)) )
out.sq <- NULL
for(i in 1:19) {
  cat("Chr ", i, "\n")
  temp <- scanqtl(hyper, chr=c(qc,i), pos=newpos,
                  formula = y ~ Q1+Q2+Q3+Q4+Q5 + Q4:Q5 + Q6)
  out.sq <- rbind(out.sq, convert.scanqtl(temp, out.fq2))
}

plot(out.sq)

out.sqi <- NULL
for(i in 1:19) {
  cat("Chr ", i, "\n")
  temp <- scanqtl(hyper, chr=c(qc,i), pos=newpos,
                  formula = y ~ Q1+Q2+Q3+Q4+Q5 + Q4:Q5 + Q6 + Q5:Q6)
  out.sqi <- rbind(out.sqi, convert.scanqtl(temp, out.fq2))
}

plot(out.sqi)
plot(out.sqi - out.sq)

############################################################
# Example 6: Internal data structure
############################################################
data(fake.bc)

class(fake.bc)

names(fake.bc)

fake.bc$pheno[1:10,]

names(fake.bc$geno)
sapply(fake.bc$geno, class)

names(fake.bc$geno[[3]])
fake.bc$geno[[3]]$data[1:5,]
fake.bc$geno[[3]]$map

names(fake.bc$geno[[3]])
fake.bc <- calc.genoprob(fake.bc, step=10, err=0.01)
names(fake.bc$geno[[3]])
fake.bc <- sim.geno(fake.bc, step=10, n.draws=8, err=0.01)
names(fake.bc$geno[[3]])
fake.bc <- argmax.geno(fake.bc, step=10, err=0.01)
names(fake.bc$geno[[3]])
fake.bc <- calc.errorlod(fake.bc, err=0.01)
names(fake.bc$geno[[3]])

names(fake.bc)
fake.bc <- est.rf(fake.bc)
names(fake.bc)

# end of rqtltour.R