library(MASS)
library(lattice)
source("boulder_library.R")

#--------------------------------------------------------------------------
#  F12 designs (AILs)
#--------------------------------------------------------------------------

# F12 with 1 qtl

markers.f12 <- read.delim("f12.markers")
ped.f12 <- read.ped("f12.ped", markers.f12$name)
chrom.scan <- scan.markers(Phenotype ~ MARKER, data=ped.f12, marker.data=markers.f12)
head(chrom.scan)

# do this if there is time
max.logPs <- permute.scan.markers(perms=100, Phenotype ~ MARKER, data=ped.f12, marker.data=markers.f12)
threshold <- quantile(max.logPs, 1 - 0.05)
# or if not...
# threshold <- 2.22

plot.chrom.scan(chrom.scan, threshold=threshold)

# find the maximum marker
chrom.scan[which.max(chrom.scan$logP),]

# draw it on the graph
abline(v=72.0)

# redo the scan conditioning on the top marker
chrom.scan <- scan.markers(Phenotype ~ m37 + MARKER, h0=Phenotype ~ m37, data=ped.f12, marker.data=markers.f12)
plot.chrom.scan(chrom.scan, col="red", add=TRUE)

# Q: why has the second peak disappered?

# F12 with two QTLs
markers.f12a <- read.delim("f12_twoqtls.markers")
ped.f12a <- read.ped("f12_twoqtls.ped", markers.f12a$name)
chrom.scan <- scan.markers(Phenotype ~ MARKER, data=ped.f12a, marker.data=markers.f12a)
plot.chrom.scan(chrom.scan, threshold=threshold)

# find the maximum peak and condition on it
chrom.scan[which.max(chrom.scan$logP),]
abline(v=72.0)

chrom.scan <- scan.markers(Phenotype ~ m37 + MARKER, h0=Phenotype ~ m37, data=ped.f12a, marker.data=markers.f12a)
plot.chrom.scan(chrom.scan, col="red", add=TRUE)

# Q: how can you test whether there are further QTLs?

chrom.scan[which.max(chrom.scan$logP),]
abline(v=56.0)

chrom.scan <- scan.markers(Phenotype ~ m37 + m29 + MARKER, h0=Phenotype ~ m37 + m29, data=ped.f12a, marker.data=markers.f12a)
plot.chrom.scan(chrom.scan, col="blue", add=TRUE)

# Q: how many QTLs are there?