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

#-------------------------------------------------------
# Looking at GxE: effects of Sex x Genotype interactions
#-------------------------------------------------------

threshold <- 2.2 # from before

# load data from before
marker.data <- read.delim("f2.markers")
ped <- read.ped("f2.ped", marker.names = marker.data$name)

# look at effect of sex
plot(Phenotype ~ Sex, data=ped)
# is it significant?
anova(lm(Phenotype ~ Sex, data=ped))

chrom.scan <- scan.markers(Phenotype ~ MARKER, data=ped, marker.data=marker.data)
plot.chrom.scan(chrom.scan, threshold=threshold)

# load data of an F2 with a sex effect independent of the QTL effect
markers.gender0 <- read.delim("f2_gender0.markers")
ped.gender0  <- read.ped("f2_gender0.ped", marker.names = markers.gender0$name)

plot(Phenotype ~ Sex, data=ped.gender0)
anova(lm(Phenotype ~ Sex, data=ped.gender0))

chrom.scan <- scan.markers(Phenotype ~ MARKER, data=ped.gender0, marker.data=markers.gender0)
plot.chrom.scan(chrom.scan, threshold=threshold)

# regress out the effects of sex
chrom.scan <- scan.markers(Phenotype ~ Sex + MARKER, h0=Phenotype ~ Sex, data=ped.gender0, marker.data=markers.gender0)
plot.chrom.scan(chrom.scan, col="green", add=TRUE)

# load data of an F2 with an interaction of sex
markers.gender1 <- read.delim("f2_gender1.markers")
ped.gender1  <- read.ped("f2_gender1.ped", marker.names = markers.gender1$name)

plot(Phenotype ~ Sex, data=ped.gender1)
anova(lm(Phenotype ~ Sex, data=ped.gender1))

chrom.scan <- scan.markers(Phenotype ~ MARKER, data=ped.gender1, marker.data=markers.gender1)
plot.chrom.scan(chrom.scan, threshold=threshold)

# regress out the effects of sex
chrom.scan <- scan.markers(Phenotype ~ Sex + MARKER, h0=Phenotype ~ Sex, data=ped.gender1, marker.data=markers.gender1)
plot.chrom.scan(chrom.scan, col="green", add=TRUE)

# do separate genome scans for males and females

females <- "F"==ped.gender1$Sex
males   <- "M"==ped.gender1$Sex

female.scan <- scan.markers(Phenotype ~ MARKER, data=ped.gender1[females,], marker.data=markers.gender1)
plot.chrom.scan(female.scan, col="magenta", add=TRUE)

male.scan <- scan.markers(Phenotype ~ MARKER, data=ped.gender1[males,], marker.data=markers.gender1)
plot.chrom.scan(male.scan, col="blue", add=TRUE)

# fit interaction on full data set

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

# what's going on? Find a marker near to the QTL
# ie, the marker with the highest p-value
chrom.scan[order(chrom.scan$logP),]

bwplot(Phenotype ~ as.factor(m12) | Sex, data=ped.gender1)
anova(lm(Phenotype ~ Sex * m13, data=ped.gender1))


#-------------------------------------
# Sex effects 2 : (answers come later)

markers.gender2 <- read.delim("f2_gender2.markers")
ped.gender2 <- read.ped("f2_gender2.ped", marker.names = markers.gender2$name)

chrom.scan <- scan.markers(Phenotype ~ MARKER, data=ped.gender2, marker.data=markers.gender2)
plot.chrom.scan(chrom.scan, threshold=threshold, ylim=c(0, 15))

# try regressing out the effect of sex
chrom.scan <- scan.markers(Phenotype ~ Sex + MARKER, h0=Phenotype ~ Sex, data=ped.gender2, marker.data=markers.gender2)
plot.chrom.scan(chrom.scan, col="green", add=TRUE)

# scan males and females separately

females <- "F"==ped.gender2$Sex
males   <- "M"==ped.gender2$Sex

male.scan <- scan.markers(Phenotype ~ MARKER, data=ped.gender2[males,], marker.data=markers.gender2)
plot.chrom.scan(male.scan, col="blue", add=TRUE)

female.scan <- scan.markers(Phenotype ~ MARKER, data=ped.gender2[females,], marker.data=markers.gender2)
plot.chrom.scan(female.scan, col="magenta", add=TRUE)

# include interaction

chrom.scan <- scan.markers(Phenotype ~ Sex * MARKER, h0 = Phenotype ~ Sex, data=ped.gender2, marker.data=markers.gender2)
plot.chrom.scan(chrom.scan, col="red", add=TRUE)

chrom.scan[order(chrom.scan$logP),]

bwplot(Phenotype ~ as.factor(m8) | Sex, data=ped.gender2)
anova(lm(Phenotype ~ Sex * m8, data=ped.gender2))