######################################################################
# Code for "A guide to QTL mapping with R/qtl"
# Karl W Broman and Saunak Sen
#
# Chapter 7: Working with covariates
#######################################################################

######################################################################
# 7.1 Additive covariates
######################################################################

# Access the gutlength data
library(qtl)
library(qtlbook)
data(gutlength)

# Omit chromosome 15
gutlength <- subset(gutlength, chr = -15)

# Box plots of phenotype by sex and cross
boxplot(gutlength ~ sex*cross, data=gutlength$pheno, 
        horizontal=TRUE, xlab="Gut length (cm)",
        col=c("red","blue"))

# ANOVA relating gut length to sex and cross
anova(aov(gutlength ~ sex*cross, data=gutlength$pheno))

# Create indicators frommom and direction
cross <- as.numeric(pull.pheno(gutlength, "cross"))
frommom <- as.numeric(cross < 3)
forw <- as.numeric(cross == 1 | cross == 3)
gutlength$pheno$frommom <- frommom
gutlength$pheno$forw <- forw

# ANOVA relating gut length to sex and frommom, and forw
anova(aov(gutlength ~ sex*frommom*forw, data=gutlength$pheno))

# Make sex numeric 
sex <- as.numeric(pull.pheno(gutlength, "sex") == "M")

# Create a matrix for the cross covariate
crossX <- cbind(frommom, forw, frommom*forw)

# Paste sex and cross together
x <- cbind(sex, crossX)

# Genome scan with and without the covariates
gutlength <- calc.genoprob(gutlength, step=1, 
                           error.prob=0.001)
out.0 <- scanone(gutlength)
out.a <- scanone(gutlength, addcovar=x)

# Plot of the lod scores with and without the covariates
plot(out.0, out.a, col=c("blue", "red"), lty=1:2, 
     ylab="LOD score", alternate.chrid=TRUE)

# Plot of the lod score differences, with and without the covariates
plot(out.a - out.0, ylab="LOD w/ covar - LOD w/o covar", 
     ylim=c(-1, 1), alternate.chrid=TRUE)
abline(h=0, lty=2)

# Create strata regarding amount of genotyping
strat <- (nmissing(gutlength) < 50)

# Permutations with and without covariates, separately for autosomes and X chr
operm.0 <- scanone(gutlength, n.perm=1000, perm.Xsp=TRUE,
                   perm.strata=strat)
operm.a <- scanone(gutlength, addcovar=x, n.perm=1000, 
                   perm.Xsp=TRUE, perm.strata=strat)

# Combine results
out.both <- c(out.0, out.a, labels=c("nocovar", "covar"))
operm.both <- cbind(operm.0, operm.a, 
                    labels=c("nocovar", "covar"))

# LOD thresholds
summary(operm.both, 0.05)

# Summary of the main results
summary(out.both, perms=operm.both, format="allpeaks",
        alpha=0.2, pvalues=TRUE)

######################################################################
# 7.2 QTL $\times$ covariate interactions
######################################################################

# Genome scan with sex as an interactive covariate.
out.i <- scanone(gutlength, addcovar=x, intcovar=sex)

# plot of LODf and LODi
plot(out.i, out.i - out.a, ylab="LOD score", 
     col=c("blue", "red"), alternate.chrid=TRUE)

# Permutations with and without sex as interactive covariate
set.seed(54955149)
operm.a <- scanone(gutlength, addcovar=x, n.perm=1000, 
                   perm.Xsp=TRUE, perm.strata=strat)
set.seed(54955149)
operm.i <- scanone(gutlength, addcovar=x, intcovar=sex,
                   n.perm=1000, perm.Xsp=TRUE, 
                   perm.strata=strat)

# Combine the results
out.ia <- c(out.i, out.i - out.a, labels=c("f","i"))
operm.ia <- cbind(operm.i, operm.i - operm.a, 
                  labels=c("f","i"))

# Loci with large LOD_f 
summary(out.ia, perms=operm.ia, alpha=0.2, pvalues=TRUE)

# Loci with large LOD_i 
summary(out.ia, perms=operm.ia, alpha=0.2, pvalues=TRUE, 
        lodcolumn=2)

# Separate genome scan within each sex
out.m <- scanone(subset(gutlength, ind = sex==1), 
                 addcovar=crossX[sex==1,])
out.f <- scanone(subset(gutlength, ind = sex==0), 
                 addcovar=crossX[sex==0,])

# plot of LOD scores in males and females separately.
plot(out.m, out.f, col=c("blue", "red"), ylab="LOD score", 
     alternate.chrid=TRUE)

# Plot LOD(males) + LOD(females) - LOD_f
plot(out.m + out.f - out.i, ylim=c(-0.5,0.5), 
     ylab="LOD(males) + LOD(females) - LODf", 
     alternate.chrid=TRUE)
abline(h=0, lty=2)

# Permutations in each of males and females
operm.m <- scanone(subset(gutlength, ind = sex==1), 
                   addcovar=crossX[sex==1,], n.perm=1000,
                   perm.strata=strat[sex==1], perm.Xsp=TRUE)
operm.f <- scanone(subset(gutlength, ind = sex==0), 
                   addcovar=crossX[sex==0,], n.perm=1000,
                   perm.strata=strat[sex==0], perm.Xsp=TRUE)

# Combine the results
out.sexsp <- c(out.m, out.f, labels=c("male","female"))
operm.sexsp <- cbind(operm.m, operm.f, 
                     labels=c("male","female"))

# LOD thresholds
summary(operm.sexsp, 0.05)

# LOD thresholds
summary(out.sexsp, perms=operm.sexsp, alpha=0.2, 
        pvalues=TRUE, format="allpeaks")

# Effect plots for four selected chromosomes
gutlength <- sim.geno(gutlength, n.draws=128, step=1,
                      error.prob=0.001)
par(mfrow=c(2,2))
effectplot(gutlength, mname1="sex", ylim=c(15.1, 17.2),
           mname2="4@65", main="Chromosome 4", 
           add.legend=FALSE)
effectplot(gutlength, mname1="sex", ylim=c(15.1, 17.2),
           mname2="5@22", main="Chromosome 5", 
           add.legend=FALSE)
effectplot(gutlength, mname1="sex", ylim=c(15.1, 17.2),
           mname2="18@48.2", main="Chromosome 18", 
           add.legend=FALSE)
effectplot(gutlength, mname1="sex", ylim=c(15.1, 17.2),
           mname2="X@58", main="X chromosome", 
           add.legend=FALSE)

######################################################################
# 7.3 Covariates with non-normal phenotypes
######################################################################

# load the nf1 data set
data(nf1)
nf1 <- drop.nullmarkers(nf1)

# Proportion affected
mean(pull.pheno(nf1,"affected"))
tapply(pull.pheno(nf1,"affected"), 
       pull.pheno(nf1,"from.mom"), mean)

# chi-square test
chisq.test(pull.pheno(nf1,"affected"), 
           pull.pheno(nf1,"from.mom"))

# genome scans for nf1
nf1 <- calc.genoprob(nf1, step=1, error.prob=0.001)
from.mom <- pull.pheno(nf1,"from.mom")
out.a <- scanone(nf1, model="binary", addcovar=from.mom)
out.i <- scanone(nf1, model="binary", addcovar=from.mom,
                 intcovar=from.mom)

# Permutation tests
set.seed(1310709)
operm.a <- scanone(nf1, model="binary", addcovar=from.mom,
                   n.perm=1000)
set.seed(1310709)
operm.i <- scanone(nf1, model="binary", addcovar=from.mom,
                   intcovar=from.mom, n.perm=1000)

# Combine the results
out.all <- c(out.i, out.a, out.i-out.a, labels=c("f","a","i"))
operm.all <- cbind(operm.i, operm.a, operm.i - operm.a,
                   labels=c("f","a","i"))

# Plot nf1 results
plot(out.all, lod=1:3, ylab="LOD score")

summary(out.all, perms=operm.all, alpha=0.2, pvalues=TRUE)

pchisq(0.703 * 2 * log(10), 1, lower=FALSE)

# nf1 scans split by parent-of-origin
out.frommom <- scanone(subset(nf1, ind=(from.mom==1)),
                       model="binary")
out.fromdad <- scanone(subset(nf1, ind=(from.mom==0)),
                       model="binary")

# Permutation tests
operm.frommom <- scanone(subset(nf1, ind=(from.mom==1)),
                       model="binary", n.perm=1000)
operm.fromdad <- scanone(subset(nf1, ind=(from.mom==0)),
                       model="binary", n.perm=1000)

# Combine the results
out.bypoo <- c(out.frommom, out.fromdad, 
               labels=c("mom","dad"))
operm.bypoo <- cbind(operm.frommom, operm.fromdad,
                     labels=c("mom","dad"))

# Plot nf1, split by parent-of-origin, results
plot(out.bypoo, lod=1:2, col=c("red", "blue"), 
     ylab="LOD score")

# Summary of the results
summary(out.bypoo, perms=operm.bypoo, alpha=0.2, pvalues=TRUE,
        format="allpeaks")

# Plot proportion affected as a function of two-locus genotypes
nf1 <- sim.geno(nf1, n.draws=128, step=1, error.prob=0.001)
par(mfrow=c(1,2))
effectplot(nf1, mname1="NPcis", mark1=1-from.mom, 
           geno1=c("Mom", "Dad"), mname2="15@13", 
           ylim=c(0,1))
effectplot(nf1, mname1="NPcis", mark1=1-from.mom, 
           geno1=c("Mom", "Dad"), mname2="19@0", ylim=c(0,1))

######################################################################
# 7.4 Composite interval mapping
######################################################################

# Reload hyper and rerun genome scan
data(hyper)
hyper <- calc.genoprob(hyper, step=1, error.prob=0.001)
out <- scanone(hyper)

# Genotype data for marker near peak LOD score
mar <- find.marker(hyper, 4, 29.5)
g <- pull.geno(hyper)[,mar]
sum(is.na(g))

# Imputate genotype data
g <- pull.geno(fill.geno(hyper))[,mar]
sum(is.na(g))

# How to recode the genotype data for an intercross [don't run here]
g <- cbind(as.numeric(g==1), as.numeric(g==2))

# genome scan with the marker as an additive covariate
out.ag <- scanone(hyper, addcovar=g)

# Plot the results
plot(out, out.ag, col=c("blue", "red"), ylab="LOD score")

# Permutation test
strat <- (ntyped(hyper) > 100)
operm.ag <- scanone(hyper, addcovar=g, chr=-4, 
                    perm.strata=strat, n.perm=1000)

# LOD thresholds
summary(operm.ag, alpha=c(0.2, 0.05))

# Summary of results
summary(out.ag, perms=operm.ag, alpha=0.2, pvalues=TRUE)

# genome scan with the marker as an interactive covariate
out.ig <- scanone(hyper, addcovar=g, intcovar=g)

# Plot the interaction LOD scores
plot(out.ig - out.ag, ylab="interaction LOD score")

# run CIM with different window sizes
out.cim.20 <- cim(hyper, n.marcovar=3, window=20)
out.cim.40 <- cim(hyper, n.marcovar=3, window=40)
out.cim.inf <- cim(hyper, n.marcovar=3, window=Inf)

# Plot the CIM results
chr <- c(1, 2, 4, 6, 15)
par(mfrow=c(3,1))
plot(out, out.cim.20, chr=chr, ylab="LOD score",
     col=c("blue", "red"), main="window = 20 cM")
add.cim.covar(out.cim.20, chr=chr, col="green")
plot(out, out.cim.40, chr=chr, ylab="LOD score",
     col=c("blue", "red"), main="window = 40 cM")
add.cim.covar(out.cim.40, chr=chr, col="green")
plot(out, out.cim.inf, chr=chr, ylab="LOD score",
     col=c("blue", "red"), main="window = Inf")
add.cim.covar(out.cim.inf, chr=chr, col="green")

# end of chap07.R