library(qtl)
data(hyper)
summary(hyper)
png("hyperplot.png",width=960,height=720)
plot(hyper)
dev.off()
png("hypermiss.png",width=960,height=720)
plot.missing(hyper, chr = c(1,4,6,15), reorder=TRUE)
dev.off()

hyper <- calc.errorlod(hyper, error.prob=0.01)
top.errorlod(hyper)
png("hypergeno.png",width=960,height=720)
plot.geno(hyper, chr=1, ind=c(117:119,137:139,157:184))
dev.off()

hyper <- calc.genoprob(hyper, step=1, error.prob=0.01)
out.em <- scanone(hyper)
out.hk <- scanone(hyper, method="hk")
summary(out.em, threshold=3)
summary(out.hk, threshold=3)

operm.hk <- scanone(hyper, method="hk", n.perm=1000)
summary(operm.hk, alpha=c(0.01,0.05))
summary(out.hk, perms=operm.hk, alpha=0.05, pvalues=TRUE)


png("hyperone.png",width=960,height=720)
plot(out.em, chr=c(1,4,6,15))
plot(out.hk, chr=c(1,4,6,15), col="blue", add=TRUE)
abline(h = summary(operm.hk, alpha=c(0.01,0.05)), col = "red", lwd = 2, lty = 2)
dev.off()

hyper <- calc.genoprob(hyper, step=5, error.prob=0.01)
out2.hk <- scantwo(hyper, method="hk")
png("hypertwo.png",width=960,height=720)
plot(out2.hk, chr=c(1,4,6,15))
dev.off()

hyper <- sim.geno(hyper, step=2, n.draws=16, error.prob=0.01)
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)
out.fitqtl$result.drop <- out.fitqtl$result.drop[,-6]
summary(out.fitqtl)

###############################################################3
library(qtlbim)
qb.load(hyper)
summary(qbHyper)

one <- qb.scanone(qbHyper, chr = c(1,4,6,15), type = "LPD")
summary(one)
plot(one)
plot(out.em,chr=c(1,4,6,15),add=T,col="red",lty = 2)

hpd <- qb.hpdone(qbHyper, profile = "2logBF")
summary(hpd)
plot(hpd)

tmp <- qb.BayesFactor(qbHyper)
summary(tmp)
plot(tmp)

two <- qb.scantwo(qbHyper, chr = c(6,15), type = "2logBF")
plot(two)
plot(two, chr = 6, slice = 15, show.locus = FALSE)
plot(two, chr = 15, slice = 6, show.locus = FALSE)
two <- qb.scantwo(qbHyper, chr = c(6,15), type = "LPD")
plot(two, chr = 6, slice = 15, show.locus = FALSE)
plot(two, chr = 15, slice = 6, show.locus = FALSE)

slice <- qb.sliceone(qbHyper, type = "cellmean",
  chr = 6, slice = 15)
summary(slice)
plot(slice)