## Adapted from Karl Broman's R/qtl Tour
## See www.rqtl.org for most recent version.

## Load package as R library.
library(qtl)

## Get hyper data.
data(hyper)

## Summary of hyper data.
summary(hyper)

## Plot hyper data.
plot(hyper)

## Enlarged plot of missing data only, reordered.
plot.missing(hyper, reorder = TRUE)

## Calculate LOD for genotyping errors. This takes time...
hyper <- calc.errorlod(hyper, error.prob=0.01)

## Show top genotyping errors.
top.errorlod(hyper)

## Plot genotypes on chr 1 for selected individuals.
plot.geno(hyper, chr=1, ind=c(117:119,137:139,157:184))

## Calculate genotype probabilities; larger steps are faster.
hyper <- calc.genoprob(hyper, step=2, error.prob=0.01)

## Scan to get interval mapping (em) and Haley-Knott (hk) LOD curves.
out.em <- scanone(hyper)
out.hk <- scanone(hyper, method="hk")

## Summaries: max per chromosome above threshold.
summary(out.em, threshold=3)
summary(out.hk, threshold=3)

## Plots
plot(out.em, chr = c(1,4,6,15))
plot(out.hk, chr = c(1,4,6,15), add = TRUE, lty = 2)

## Permutation thresholds (much faster with hk).
operm.hk <- scanone(hyper, method="hk", n.perm=1000)
summary(operm.hk, alpha=c(0.01,0.05))
abline(h = summary(operm.hk, alpha=c(0.01,0.05)), lty = 3)

## Summary with p-value column.
summary(out.hk, perms=operm.hk, alpha=0.05, pvalues=TRUE)

## Calculate genoprob on coarser grid (faster) for 2-D scan.
hyper <- calc.genoprob(hyper, step=5, error.prob=0.01)
## 2-D scan with Haley-Knott (faster than em).
out2.hk <- scantwo(hyper, method="hk")

## 2-D Summaries.
summary(out2.hk, thresholds=c(6.0, 4.7, 4.4, 4.7, 2.6))
summary(out2.hk, thresholds=c(6.0, 4.7, 4.4, 4.7, 2.6), what = "int")

## 2-D Plot.
plot(out2.hk, chr=c(1,4,6,15))

## ANOVA for multiple QTL
hyper <- sim.geno(hyper, step=2, n.draws=16, error.prob=0.01)
qtl <- makeqtl(hyper, chr = c(1, 1, 4, 6, 15), pos = c(50, 76, 30, 70, 20))

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

## Effect plots and interaction plot.
hyper <- sim.geno(hyper, step=2, n.draws=16, error.prob=0.01)
effectplot(hyper, pheno.col = 1, mname1 = "D1Mit334")
effectplot(hyper, pheno.col = 1, mname1 = "D4Mit164")
markers <- find.marker(hyper, chr = c(6,15), pos = c(70,20))
markers
effectplot(hyper, pheno.col = 1, mname1 = markers[1], mname2 = markers[2])
effectplot(hyper, pheno.col = 1, mname1 = markers[2], mname2 = markers[1])

## Strip plot of data
plot.pxg(hyper, "D1Mit334")
plot.pxg(hyper, "D4Mit164")
plot.pxg(hyper, markers)