#code to replicate simulations based on hypertension data set in Sugiyama et al. (2001) 
#run with R version 2.4.1 and R/qtl version 1.07 
library(pLODsims)
data(map10)

#This is the first of sixteen simulations run in parallel
#For the full set of simulations, replace 01 with the values 02, ..., 16
set.seed(52967665+01)

n.sim <- 250

simres.tgb <-  list()

#specify model parameters
RMSE<-sqrt(0.9363305)
scale <- 0.6 / RMSE
aeffects <- c(0.08177, 0.41753, -0.56086, 0.24533, 0.17272, -0.51780, 0.55867) * scale # scaled to have sd=1
deffects <- c(-0.28760, -0.24667, -0.09485, -0.46404, -0.13206, -0.19185, -0.33115) * scale
i1effects <- c(0.11936, -1.23483, 0.06423, 1.04768) * scale #first interaction aa, ad, da, dd
i2effects <- c(-0.14228, 1.06596, 0.33489, -0.60775) * scale #first interaction aa, ad, da, dd
i3effects <- c(0.49512, -0.47932, 0.11983, -1.18010) * scale #first interaction aa, ad, da, dd

qtl.pos <- cbind(c(5,7,8,10,16,17,18), c(80.0, 63.5, 14.0, 6.0, 48.2, 11.75, 30.0))
main.effects <- cbind(1:7, aeffects, deffects)

intxns <- rbind(c(2,4, i1effects), c(3,7,i2effects), c(1,5,i3effects))

for(i in 1:n.sim) {
  tgbsim <- pLODsims.sim.cross(map10[1:19], qtl.pos=qtl.pos, qtl.main=main.effects, qtl.epis=intxns, n.ind=500, type="f2")
  tgbsim <- calc.genoprob(tgbsim, step=2.5)

  #pLOD_a with 5% thresholds
  simres.tgb[[i]] <- search.hk(tgbsim, T1=3.5808, markers.only=1, refinement=1, 
	max.steps=100, max.df=50, fs2=1, verbose=0, pheno.col=1, restrict.diag=2, incl.intxns=0)
                   
  if(i==floor(i/10)*10){
    cat(i,"\n")
    save(simres.tgb, file="tgb01.RData", compress=T)
 }
}

save(simres.tgb, file="tgb01.RData", compress=T)