control_for_qtl.R

library(qtl)
library(parallel)

# fake intercross data
data(fake.f2)

# single-QTL genome scan
fake.f2 <- calc.genoprob(fake.f2, step=1)
out <- scanone(fake.f2, method="hk")

# permutation test in parallel
operm <- scanone(fake.f2, method="hk", n.perm=1000,
                 n.cluster=detectCores())

# single-QTL scan summary
summary(out, perms=operm, alpha=0.05)

# maximum peak
mx <- max(out) # chr 13, 28 cM

# make a QTL object
qtl <- makeqtl(fake.f2, chr=mx$chr, pos=mx$pos, what="prob")

# scan for an additional QTL
out.c13 <- addqtl(fake.f2, qtl=qtl, method="hk")
summary(out.c13, perms=operm, alpha=0.05)

# scan for an additional QTL, interacting with c13 locus
out.c13i <- addqtl(fake.f2, qtl=qtl, method="hk",
                   formula=y~Q1*Q2)

# plot LOD scores for interaction
plot(out.c13i - out.c13, ylab="interaction LOD")