##############################
# Preliminaries
##############################
library(qtl)
ls()
help(read.cross)
?read.cross
options(htmlhelp=TRUE)

##############################
# Example 1: Hypertension
###############################

# (1) Get access to the data
data(hyper)
ls()
?hyper

# (2) Summary information on the data
summary(hyper)

nind(hyper)
nphe(hyper)
nchr(hyper)
totmar(hyper)
nmar(hyper)

# (3) Plot summary information on the data
plot(hyper)

plot.missing(hyper)
plot.map(hyper)
hist(hyper$pheno[,1], breaks=30)

# (4) Plot of missing data, individuals in phenotype order
plot.missing(hyper,reorder=TRUE)

# (5) Drop markers with no data
hyper <- drop.nullmarkers(hyper)
totmar(hyper)

# (6) Estimate and plot pairwise recombination fractions
hyper <- est.rf(hyper)
plot.rf(hyper)
plot.rf(hyper,c(1,4))

plot.rf(hyper,6)
plot.missing(hyper,6)

# (7) Re-estimate genetic map
newmap <- est.map(hyper, error.prob=0.01, trace=TRUE)
plot.map(hyper, newmap)

hyper <- replace.map(hyper, newmap)

# (8) Identify potential genotyping errors
hyper <- calc.genoprob(hyper, error.prob=0.01)
hyper <- calc.errorlod(hyper, error.prob=0.01)

plot.errorlod(hyper)
top.errorlod(hyper)
plot.errorlod(hyper, chr=c(4,11,16))

# (9) Plot the questionable genotypes
plot.geno(hyper, chr=16, ind=71:90, min.sep=4)

# (10) Plot proporition of available marker information
hyper <- calc.genoprob(hyper, step=2, err=0.01)
plot.info(hyper)
plot.info(hyper, chr=c(1,4,15))
plot.info(hyper, chr=c(1,4,15), method="entropy")
plot.info(hyper, chr=c(1,4,15), method="variance")

# (11) Single-QTL genome scan by EM and Haley-Knott
out.em <- scanone(hyper)
out.hk <- scanone(hyper, method="hk")

# (12) Summary of the results
summary(out.em)
summary(out.em, 3)
summary(out.hk, 3)

# (13) Just the highest peaks
max(out.em)
max(out.hk)

# (14) Plot the results
plot(out.em, chr=c(1,4,15))
plot(out.hk, out.em, chr=c(1,4,15), col=c("red","blue"), lty=1)

plot(out.em, chr=c(1,4,15), col="blue")
plot(out.hk, chr=c(1,4,15), col="red", add=TRUE)

# (15) Multiple imputation
hyper <- sim.geno(hyper, step=2, n.draws=32, err=0.01)
out.imp <- scanone(hyper, method="imp")
plot(out.hk, out.em, out.imp, chr=c(1,4,15), col=c("red","blue","orange"), lty=1)

# (16) Permutation test
operm.hk <- scanone(hyper, method="hk", n.perm=25)
quantile(operm.hk, 0.95)

# (17) Save what you have done
save.image()

# (18)  Two-dimensional, 2-QTL genome scan
hyper.coarse <- calc.genoprob(hyper, step=10, error.prob=0.01)
out2.hk <- scantwo(hyper.coarse, method="hk")

# (19) Summary and plot of the results
summary(out2.hk, c(8,3,3))
summary(out2.hk, c(0,4,Inf))
summary(out2.hk, c(0,Inf,4))

plot(out2.hk)
plot(out2.hk,chr=c(1,4))

# (20) Peaks with maximum joint LOD and interactive LOD
max(out2.hk)

# (21) Permutation test for 2-dim scan
operm2 <- scantwo(hyper.coarse, method="hk", n.perm=3)
apply(operm2, 2, quantile, 0.95)

# (22) Fit of a multiple-QTL model by imputation
chr <- c(1, 1, 4, 6, 15)
pos <- c(50, 76, 30, 70, 20)
qtl <- makeqtl(hyper, chr, pos)

my.formula <- y ~ Q1 + Q2 + Q3 + Q4 + Q5 + Q4:Q5
out.fitqtl <- fitqtl(hyper$pheno[,1], qtl, formula=my.formula)
summary(out.fitqtl)

##############################
# Example 2: Genetic mapping
##############################

# (1) Get access to the data; get summaries
data(badorder)

summary(badorder)
plot(badorder)

# (2) Estamite pairwise recombination fractions
badorder <- est.rf(badorder)
plot.rf(badorder)

plot.rf(badorder, chr=1)

# (3) Re-estimate genetic map
newmap <- est.map(badorder, trace=TRUE)
plot.map(badorder, newmap)

# (4) Assess marker order on chromosome 1
rip1 <- ripple(badorder, chr=1, window=6)
summary(rip1)

# (5) Assess marker order with likelihood method
rip2 <- ripple(badorder, chr=1, window=3, err=0.01, method="likelihood")
summary(rip2)

# (6) Switch markers 9-11 and re-assess marker order
badorder.rev <- switch.order(badorder, 1, rip1[2,])
rip1r <- ripple(badorder.rev, chr=1, window=6)
summary(rip1r)

# (7) Switch markers and re-assess marker order
badorder.rev <- switch.order(badorder.rev, 1, rip1r[2,])
rip2r <- ripple(badorder.rev, chr=1, window=3, err=0.01)
summary(rip2r)

# (8) Re-estimate pairwise recombination fractions
badorder.rev <- est.rf(badorder.rev)
plot.rf(badorder.rev, 1)


##############################
# Internal data structure
##############################

# (1) Get access to data
data(fake.bc)

# (2) class
class(fake.bc)

# (3) geno and pheno compontents
names(fake.bc)

# (4) the phenotype data
fake.bc$pheno

# (5) The geno component
names(fake.bc$geno)
sapply(fake.bc$geno, class)

# (6) The geno component: data and map
names(fake.bc$geno[[3]])
fake.bc$geno[[3]]$data[1:5,]
fake.bc$geno[[3]]$map

# (7) The geno component after calculations
names(fake.bc$geno[[3]])

fake.bc <- calc.genoprob(fake.bc, step=10, err=0.01)
names(fake.bc$geno[[3]])

fake.bc <- sim.geno(fake.bc, step=10, n.draws=10, err=0.01)
names(fake.bc$geno[[3]])

fake.bc <- argmax.geno(fake.bc, step=10, err=0.01)
names(fake.bc$geno[[3]])

fake.bc <- calc.errorlod(fake.bc, err=0.01)
names(fake.bc$geno[[3]])

# (8) The rf component
names(fake.bc)
fake.bc <- est.rf(fake.bc)
names(fake.bc)