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)

threshold <- 2.22

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

#-----------------------------------------
# Stop and look at results: how many QTLs?
#-----------------------------------------


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

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

#--------------------
# STOP HERE AND WAIT FOR EXPLANATION
#

# fit two markers separately
anova(lm(Phenotype ~ m32, data=ped.f12))
anova(lm(Phenotype ~ m37, data=ped.f12))

# fit m32 conditional on m37
# (ie, fit m37 first, then see if m32 adds anything more)
fit0 <- lm(Phenotype ~ m37, data=ped.f12)
fit1 <- lm(Phenotype ~ m37 + m32, data=ped.f12)
anova(fit0, fit1)
# or more simply
anova(fit1)


# 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)

#---------------------------------------------------------------
# Stop and look at results: why has the second peak disappeared?
#---------------------------------------------------------------

# 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)

#-------------------------
# Stop and look at results
#-------------------------

# 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?