## Data file: ftp://statgen.ncsu.edu/pub/qtlcart/data/zengetal99/
## Need to convert files to CSV.
## Change 9, 999.99900 to missing (empty value).

bm4zb <- read.csv("bm4zb.csv", header = TRUE)
bm6zb <- read.csv("bm6zb.csv", header = TRUE)
bs4zb <- read.csv("bs4zb.csv", header = TRUE)
bs6zb <- read.csv("bs6zb.csv", header = TRUE)

for(i in c("bm4zb","bm6zb","bs4zb","bs6zb")) {
  cat(i, "\n")
  tmp <- read.csv(paste(i, "csv", sep = "."), header = TRUE)

  ## Change X chr values.
  tmp[, 2:7] <- (tmp[, 2:7] / 2)

  ## Change chr 2,3 values to 1,2.
  if(substring(i, 1, 2) == "bm")
    tmp[, 8:47] <- 2 - tmp[, 8:47]

  ## Reorder for read.cross
  tmp <- tmp[, c(1, 47:51, 2:46)]
  cat(names(tmp), "\n")

  ## Add chromosome.
  tmp <- rbind(c(rep(NA, 6), rep(1, 6), rep(2, 16), rep(3, 23)), tmp)

  ## Write cross file.
  write.csv(tmp, file = paste("n", i, ".csv", sep = ""), row.names = FALSE, na = "")
}

bm4zb <- cbind(run = 4, read.csv("nbm4zb.csv", header = TRUE))
bm6zb <- cbind(run = 6, read.csv("nbm6zb.csv", header = TRUE))
bm4zb$label <- as.character(bm4zb$label)
bm6zb$label <- as.character(bm6zb$label)
bmzb <- rbind(bm4zb, bm6zb)
bmzb <- bmzb[-(nrow(bm4zb)+1),]
bmzb[1,1] <- NA
bmzb$label <- factor(bmzb$label)
write.csv(bmzb, file = "bmzb.csv", row.names = FALSE, na = "", quote = FALSE)

bs4zb <- cbind(run = 4, read.csv("nbs4zb.csv", header = TRUE))
bs6zb <- cbind(run = 6, read.csv("nbs6zb.csv", header = TRUE))
bs4zb$label <- as.character(bs4zb$label)
bs6zb$label <- as.character(bs6zb$label)
bszb <- rbind(bs4zb, bs6zb)
bszb <- bszb[-(nrow(bs4zb)+1),]
bszb[1,1] <- NA
bszb$label <- factor(bszb$label)
write.csv(bszb, file = "bszb.csv", row.names = FALSE, na = "", quote = FALSE)

######################################################################################3
library(qtl)
bmzb <- read.cross("csv",, "bmzb.csv", genotype = c(0,1))
######################################################################################3

bmzb <- calc.genoprob(bmzb, step = 2)
addcov <- as.matrix(bmzb$pheno$run)
plot(scanone(bmzb,, 3, addcov = addcov))

two <- scantwo(bmzb,, 3, addcov = addcov, type = "hk")
png("bmzbtwo.png", width = 960, height = 960)
plot(two, col = "gray")
dev.off()
######################################################################################3
library(qtlbim)
bmzb <- read.cross("csv",, "bmzb.csv", genotype = c(0,1))

bmzb <- qb.genoprob(bmzb, step = 2)
qb <- qb.mcmc(bmzb, pheno.col = 3, fixcov = 1)
save(qb, file = "qb.bmzb.RData", compress = TRUE)
qb <- qb.split.chr(qb)
one <- qb.scanone(qb, type = "2logBF")
lpd <- qb.scanone(qb, type = "LPD")

bszb <- read.cross("csv",, "bszb.csv", genotype = c(0,1))
bszb <- qb.genoprob(bszb, step = 2)
qb.bs <- qb.mcmc(bszb, pheno.col = 3, fixcov = 1)
save(qb.bs, file = "qb.bszb.RData", compress = TRUE)
one.bs <- qb.scanone(qb.bs, type = "2logBF")

png("bmzbqb.png", width = 960, height = 960)
par(mfrow = c(2,1), mar = c(3.1,4.1,3.1,0.1))
plot(one, col = c(sum="darkgray", epistasis = "black", main = "black"), lty = c(3,2,1),
     main = "mauritiana bc")
plot(one.bs, col = c(sum="darkgray", epistasis = "black", main = "black"), lty = c(3,2,1),
     main = "simulans bc")
dev.off()