######################################################################
#
# hort.r
#
# copyright (c) 2004, Brian S. Yandell, University of Wisconsin-Madison
# Licensed under the GNU General Public License version 2 (June, 1991)
# 
# Relies on R/qtl and R/bim, both available from cran.r-project.org
#
######################################################################


library(qtlbim)
################################################################################
## data setup
data(Bnapus)
## 1 = Major has longer flowering time (biennial)
## 2 = Stellar does not require OW (annual)

## two selected markers (known to be near QTL)
mpair <- find.marker(Bnapus,chr=rep("N2",2),pos=c(54,110))

markers <- as.data.frame(argmax.geno(Bnapus,1)$geno[[2]]$argmax[,mpair])
for(i in 1:2)
  markers[[i]] <- factor(markers[[i]],labels=c("M","S"))

fl <- data.frame(f4 = Bnapus$pheno$flower4, f8 = Bnapus$pheno$flower8, f0 = Bnapus$pheno$flower0)
fl1 <- fl[markers[[1]] == "M",]
fl2 <- fl[markers[[1]] == "S",]

################################################################################

plot.map(Bnapus)
points(rep(2,2),pull.map(Bnapus)[[2]][mpair],col="red",lwd=2,cex=2)

## counts: recombination (filling in most probably genotype for missing values
tbl <- table(markers,exclude=NULL)
print(tbl)
four <- c(tbl)
names(four) <- c(outer(dimnames(tbl)[[1]],dimnames(tbl)[[2]],paste,sep=""))
print(four)

################################################################################
## Figure 1: Pie Charts and Bar Plots
jpeg("fig1.jpg",width=960,height=960)
par(mar=c(2.5,2.5,0.1,0.1),mfrow=c(2,2), lwd = 2)
pie(four[c(2,1,3,4)],names(four)[c(2,1,3,4)],col="white",cex=1.5)
text(-.4,.9,"M1 = Major",cex=1.5)
text(.4,.9,"M1 = Stellar",cex=1.5)
text(.9,.4,"M2 = Major",srt=270,cex=1.5)
text(.9,-.4,"M2 = Stellar",srt=270,cex=1.5)
text(-1,1,"(a)",cex=1.5)
text(-1,-1,"(c)",cex=1.5)
pie(four[c(1,4,2,3)],names(four)[c(1,4,2,3)],col="black",
  density=c(10,20,30,40),angle=c(45,90,135,180),cex=1.5)
text(0,.9,"non-recombinants",cex=1.5)
text(0,-.9,"recombinants",cex=1.5)
text(-1,1,"(b)",cex=1.5)
text(-1,-1,"(d)",cex=1.5)

barplot(four,angle=c(45,135,45,135),density=c(20,40,40,20),col="black")
barplot(tbl,angle=c(135,45,135,45),density=c(20,0),col="black")
#  legend=TRUE)
mtext(mpair[1],1,1,cex=1.25)
#text(1.4,55,mpair[2],cex=1.5, adj = 0)
dev.off()

################################################################################
## Figure 2: Recombination frequency and LOD score
Bnapus <- est.rf(Bnapus)
print(Bnapus$rf[mpair,mpair])

## NB: Figure uses as customized version of plot.rf().

jpeg("fig2.jpg",width=960,height=480)
par(mar=c(2.5,2.5,2.5,2.5),mfrow=c(1,2), lwd = 2)
plot.rf(Bnapus,chr="N2",col="gray", main = "")
abline(0,1,col="white")
mtext("(a)",1,1,at=0,cex=1.5)
plot.rf(Bnapus,col="gray",cex=0.75, chr = c(2,7))
abline(0,1,col="white")
mtext("(b)",1,1,at=0,cex=1.5)
dev.off()

################################################################################
## Figure 3. Histograms
## sumaries
stem(fl2$f4,2)
stem(fl2$f8)

fl.hist <- list()
jpeg("fig3.jpg",width=960,height=960)
par(mfrow=c(2,2),mar=c(2.1,2.1,2.1,0.1), lwd = 2)
tmp <- hist(fl1$f4,breaks=seq(0,100,by=5),xlab="",main="(a) Major 4-week")
fl.hist$days <- tmp$breaks[-1] -2.5
fl.hist$m4 <- tmp$counts
tmp <- hist(fl1$f8,breaks=seq(0,100,by=5),xlab="",main="(b) Major 8-week")
fl.hist$s4 <- tmp$counts
tmp <- hist(fl2$f4,breaks=seq(0,100,by=5),xlab="",main="(c) Stellar 4-week")
fl.hist$m8 <- tmp$counts
tmp <- hist(fl2$f8,breaks=seq(0,100,by=5),xlab="",main="(d) Stellar 8-week")
fl.hist$s8 <- tmp$counts
fl.hist <- as.data.frame(fl.hist)
dev.off()

################################################################################
## Figure 4. Q-Q plots

myfun <- function(data, ylim = range(data, na.rm = TRUE), main = "")
{
  tmp <- qqnorm(subset(data, !is.na(data)), cex = 1.5, lwd = 2, type = "n",
    main = main, ylim = ylim)
  lines(sort(tmp$x), sort(tmp$y), type = "b", cex = 1.5, lwd = 2)
  tmpx <- quantile(tmp$x, c(.25, .75))
  tmpy <- quantile(tmp$y, c(.25, .75))
  abline(mean(tmpy),diff(tmpy)/diff(tmpx), col = "darkgray", lwd = 2)
}

jpeg("fig4.jpg",width=960,height=960)
par(mfrow = c(2,2), lwd = 2, mar = c(2.5,2.5,3.5,0.1))
myfun(fl1$f4, range(fl1$f4, fl2$f4, na.rm = TRUE), "(a) Major 4-week")
myfun(fl2$f4, range(fl1$f4, fl2$f4, na.rm = TRUE), "(b) Stellar 4-week")
myfun(fl1$f8, range(fl1$f8, fl2$f8, na.rm = TRUE), "(c) Major 8-week")
myfun(fl2$f8, range(fl1$f8, fl2$f8, na.rm = TRUE), "(d) Stellar 8-week")
dev.off()

################################################################################
## Figure 5. pheno by geno plot
jpeg("fig5.jpg",width=960,height=480)
par(yaxt="n",mar=c(4.5,4.5,0.1,0.1),mfrow=c(1,2), lwd = 2, cex = 1.5)
tmp2 <- plot.pxg(Bnapus,mpair,11,ylab="8-week",xlab="", col = "black")
par(yaxt="s")
axis(2,log10(c(10,20,30,40,50)),c(10,20,30,40,50))
mtext("(a)",1,1,at=0.25,cex=1.5)
#tmp2$p
par(yaxt="n")
tmp2 <- plot.pxg(Bnapus,mpair,10,ylab="4-week",xlab="", col = "black")
par(yaxt="s")
axis(2,log10(c(10,20,30,40,50,75,100)),c(10,20,30,40,50,75,100))
mtext("(b)",1,1,at=0.25,cex=1.5)
#text(3.5,log10(100),"no flowering")
#tmp2$p
dev.off()

################################################################################
## Figure 6. scatter plots
jpeg("fig6.jpg",width=960,height=480)
par(mar=c(3.1,3.1,0.1,0.1),mfrow=c(1,2),pty="s")
plot(fl$f8,fl$f4,xlab="",ylab="",cex=1.75)
mtext("8-week vern flower time",1,2)
mtext("4-week vern flower time",2,2)
text(45,15,"no labels\nskew to right",adj = 1,cex=1.5)
mtext("(a)",1,1,at=5,cex=1.5)
plot(fl$f8,fl$f4,xlab="",ylab="",log="xy",type="n")
mtext("8-week vern flower time",1,2)
mtext("4-week vern flower time",2,2)
text(fl$f8,fl$f4,as.character(markers[[1]]),cex=1,
  )
#  col = c("red","blue")[unclass(markers[[1]])])
text(45,10,"group labels\nmore symmetric",adj = 1,cex=1.5)
mtext("(b)",1,1,at=7.5,cex=1.5)
dev.off()

################################################################################
## Figure 7. separate plots by genotype
tmpfn <- function(fl1)
{
  fit <- summary(lm(log10(f8)~log10(f4),fl1))
  std.dev <- fit$sigma
  x0 <- log10(10)+c(-.01,0,-.01)
  y0 <- log10(20)+c(0,0,std.dev)
  x1 <- log10(10)+c(.01,0,.01)
  y1 <- log10(20)+c(0,rep(std.dev,2))
  segments(10^x0,10^y0,10^x1,10^y1, lwd = 2)
  text(10^x1[1],10^mean(y0[2:3]),"SD",adj=0)
  print(std.dev)
}

jpeg("fig7.jpg",width=960,height=480)
par(mar=c(3.1,3.1,0.1,0.1),mfrow=c(1,2),pty="s")
plot(fl1$f4,fl1$f8,xlab="",ylab="",log="xy", lwd = 2, cex = 1.5, col = "gray",
  xlim = range(fl$f4, na.rm = TRUE ), ylim = range( fl$f8, na.rm = TRUE ))
mtext("4-week vern flower time",1,2)
mtext("8-week vern flower time",2,2)
r <- range( fl2$f4, na.rm = TRUE )
tmp1 <- predict(lm(log10(f8)~log10(f4),fl2),newdata=data.frame(f4=r))
lines(r,10^tmp1,col="black",lwd = 2, lty = 2)
r <- range( fl1$f4, na.rm = TRUE )
tmp2 <- predict(lm(log10(f8)~log10(f4),fl1),newdata=data.frame(f4=r))
lines(r,10^tmp2, col="black",lwd = 2, lty = 1)
r <- seq(r[1], r[2], length = 100)
tmp1 <- predict(lm(log10(f8)~log10(f4),fl2), se.fit = TRUE,
  newdata=data.frame(f4 = r))
tmp2 <- predict(lm(log10(f8)~log10(f4),fl1), se.fit = TRUE,
  newdata=data.frame(f4 = r))
tmp <- sqrt(tmp1$se.fit^2 + tmp2$se.fit^2)
lines(r,10^(tmp2$fit - 2 * tmp), col="black",lwd = 2)
lines(r,10^(tmp2$fit + 2 * tmp), col="black",lwd = 2)
text(50,10,"Major genotype")
tmpfn(fl1)

plot(fl2$f4,fl2$f8,xlab="",ylab="",log="xy", lwd = 2, cex = 1.5, col = "gray",
  xlim = range(fl$f4, na.rm = TRUE ), ylim = range( fl$f8, na.rm = TRUE ))
mtext("4-week vern flower time",1,2)
mtext("8-week vern flower time",2,2)
r <- range( fl1$f4, na.rm = TRUE )
lines(r,10^predict(lm(log10(f8)~log10(f4),fl1),newdata=data.frame(f4=r)),
  col="black",lwd = 2, lty = 2)
r <- range( fl2$f4, na.rm = TRUE )
lines(r,10^predict(lm(log10(f8)~log10(f4),fl2),newdata=data.frame(f4=r)),
  col="black",lwd = 2)
r <- seq(r[1], r[2], length = 100)
tmp1 <- predict(lm(log10(f8)~log10(f4),fl1), se.fit = TRUE,
  newdata=data.frame(f4 = r))
tmp2 <- predict(lm(log10(f8)~log10(f4),fl2), se.fit = TRUE,
  newdata=data.frame(f4 = r))
tmp <- sqrt(tmp1$se.fit^2 + tmp2$se.fit^2)
lines(r,10^(tmp2$fit - 2 * tmp), col="black",lwd = 2)
lines(r,10^(tmp2$fit + 2 * tmp), col="black",lwd = 2)
text(50,10,"Stellar genotype")
tmpfn(fl2)
dev.off()
################################################################################
## Figure 8: interaction plot
jpeg("fig8.jpg",width=960,height=480)
par(mfrow=c(2,1),mar=c(2.1,2.1,0.1,0.1))
tmp=Bnapus
names(tmp$pheno)=c(names(tmp$pheno)[-(10:11)],paste(c(4,8),"-week flowering time",sep=""))
par(yaxt="n",mar=c(4.5,4.5,0.5,0.5),mfrow=c(1,2))
effectplot(tmp,10,mpair[2],mname2=mpair[1],main="")
par(yaxt="s")
axis(2,log10(seq(20,60,10)),seq(20,60,10))
mtext("(a)",1,1,at=0.65,cex=1.5)
par(yaxt="n")
effectplot(tmp,11,mpair[2],mname2=mpair[1],main="")
par(yaxt="s")
axis(2,log10(seq(10,30,2)),seq(10,30,2))
mtext("(b)",1,1,at=0.65,cex=1.5)
dev.off()
################################################################################
## other miscellaneous code
## box plot
par(mfrow=c(2,1),mar=c(2.1,2.1,0.1,0.1))
hist(fl2$f8,nclass=20,xlab="",main="")
text(40,10,"histogram",cex=1.5)
boxplot(fl2$f8,horizontal=TRUE)
text(40,1.4,"boxplot",cex=1.5)

x11(width=10,height=5)
par(mfrow=c(1,1),mar=c(2.1,2.1,0.1,0.1))
plot(fl2$f8,jitter(rep(1,nrow(fl2)),,0.1),xlab="",ylab="",yaxt="n",ylim=c(0.5,2.5),cex=1.5,lwd=2)
boxplot(fl2$f8,horizontal=TRUE,add=TRUE,at=2,cex=1.5,lwd=2)

tmp <- sample(fl2$f8,5)
plot(tmp,jitter(rep(1,length(tmp)),,0.1),xlab="",ylab="",yaxt="n",ylim=c(0.5,2.5),cex=1.5,lwd=2)
boxplot(tmp,horizontal=TRUE,add=TRUE,at=2,cex=1.5,lwd=2)


## box plots
fl0.box <- data.frame(day=Bnapus$pheno$flower0,marker = marker)
fl0.box$vern <- rep(0,nrow(fl0.box))
fl0.box <- fl0.box[marker==1|marker==2,]
fl4.box <- data.frame(day=fl$f4,marker = marker)
fl4.box$vern <- rep(4,nrow(fl4.box))
fl4.box <- fl4.box[marker==1|marker==2,]
fl8.box <- data.frame(day=fl$f8,marker = marker)
fl8.box$vern <- rep(8,nrow(fl8.box))
fl8.box <- fl8.box[marker==1|marker==2,]
fl.box <- rbind(fl0.box,fl4.box,fl8.box)
fl.box$marker <- factor(fl.box$marker)
fl.box$vern <- ordered(fl.box$vern)
x11(width=10,height=5)
par(mar=c(3.1,3.1,0.1,0.1))
boxplot(day~marker*vern,fl.box,names=rep(c("Stellar","Major"),3))
mtext(c("no vern","4-week","8-week"),1,2,at=c(1.5,3.5,5.5))
mtext("days",2,2)
text(5,100,"no flowering",adj=0)

## adjust for skew
boxplot(log10(day)~marker*vern,fl.box,names=rep(c("Stellar","Major"),3),yaxt="n")
axis(2,log10(c(10,20,50,100)),c(10,20,50,100))
mtext(c("no vern","4-week","8-week"),1,2,at=c(1.5,3.5,5.5))
mtext("days",2,2)
text(5,2,"no flowering",adj=0)

x11(width=10,height=5)
jpeg("jitter.jpg",width=960,height=480)
par(mar=c(3.1,3.1,0.1,0.1))
tmp <- interaction(fl.box[,c("marker","vern")])
plot(c(0.5,6.5),range(fl.box$day,na.rm=T),type="n",xlab="",ylab="",xaxt="n")
axis(1,1:6,rep(c("Stellar","Major"),3))
tmp2 <- tapply(fl.box$day,tmp,median,na.rm=T)
segments((1:6)-.35,tmp2,(1:6)+.35,tmp2,col="blue",lwd=3)
tmp2 <- fl.box$day
tmp2[tmp2==100] <- jitter(tmp2[tmp2==100],,2)
points(jitter(unclass(tmp)),tmp2,lwd=2,cex=1.5)
text(5,100,"no flowering",adj=0)
mtext("days",2,2)
mtext(c("no vern","4-week","8-week"),1,2,at=c(1.5,3.5,5.5))
dev.off()

## log
tmp <- interaction(fl.box[,c("marker","vern")])
plot(c(0.5,6.5),range(fl.box$day,na.rm=T),type="n",xlab="",ylab="",xaxt="n",log="y")
axis(1,1:6,rep(c("Stellar","Major"),3))
tmp2 <- tapply(fl.box$day,tmp,median,na.rm=T)
segments((1:6)-.35,tmp2,(1:6)+.35,tmp2,col="blue",lwd=3)
tmp2 <- fl.box$day
tmp2[tmp2==100] <- jitter(tmp2[tmp2==100],,2)
points(jitter(unclass(tmp)),tmp2,lwd=2,cex=1.5)
text(5,100,"no flowering",adj=0)
mtext("days",2,2)
mtext(c("no vern","4-week","8-week"),1,2,at=c(1.5,3.5,5.5))

## pheno by geno plot
x11(width=10,height=6)
par(yaxt="n",xaxt="n")
tmp2 <- plot.pxg(Bnapus,mpair,11,ylab="",xlab="")
mtext("genotype",1,2,cex=1.5)
mtext("phenotype",2,2,cex=1.5)
par(yaxt="s")
axis(2,log10(c(10,20,30,40,50)),c(10,20,30,40,50))
tmp2$p


################################################################################
## pairs plot
lower <- function(x,y) {
text(x,y,c("SS","SM","MS","MM")[marker],cex=0.65,
  col = c("red","green","purple","blue")[marker])
}
x11(height=7,width=7)
par(pty="s")
tmp <- as.data.frame(lapply(Bnapus$pheno[,6:8],jitter))
pairs(tmp,lower.panel=lower,labels=c("no vern","4-week","8-week"))




## 3 comps
jpeg("three.jpg",width=960,height=960)
tmp2 <- fl$f0
tmp2[tmp2==100] <- jitter(tmp2[tmp2==100],,5)
par(mar=c(3.1,3.1,0.1,0.1),mfrow=c(2,2),pty="s")
plot(fl$f4,tmp2,xlab="",ylab="",log="xy",type="n")
mtext("4-week",1,2)
mtext("no vern",2,2)
text(fl$f4,tmp2,as.character(markers[[1]]),cex=1,
  col = c("red","blue")[unclass(markers[[1]])])
abline(h=90,v=90,lty=2,lwd=2)
plot(fl$f8,tmp2,xlab="",ylab="",log="xy",type="n")
mtext("no vern",2,2)
mtext("8-week",1,2)
text(fl$f8,tmp2,as.character(markers[[1]]),cex=1,
  col = c("red","blue")[unclass(markers[[1]])])
abline(h=90,lty=2,lwd=2)
plot(fl$f4,fl$f8,xlab="",ylab="",log="xy",type="n")
mtext("4-week",1,2)
mtext("8-week",2,2)
text(fl$f4,fl$f8,as.character(markers[[1]]),cex=1,
  col = c("red","blue")[unclass(markers[[1]])])
abline(v=90,lty=2,lwd=2)
dev.off()



## regression diagnostics
fl.fit <- lm(log10(f8)~log10(f4),fl,subset=marker==1)

## 2-D QTL scan
Bnapus=calc.genoprob(Bnapus,2,error=0.01)
fl4 <- scantwo(Bnapus,pheno=10,chr=c(2,3))
summary(fl4)
fl8 <- scantwo(Bnapus,pheno=11,chr=2)
summary(fl8)

###################################################################################
## lattice (pick out some examples!)
lset(theme = col.whitebg())
par(ask=T)
example(xyplot)
## need to figure out how to get reasonable headings on this!
xyplot(f4~f8|m1*m2,fl,cex=1.5)

###################################################################################