R : Copyright 2001, The R Development Core Team
Version 1.2.2  (2001-02-26)

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type `license()' or `licence()' for distribution details.

R is a collaborative project with many contributors.
Type `contributors()' for more information.

Type `demo()' for some demos, `help()' for on-line help, or
`help.start()' for a HTML browser interface to help.
Type `q()' to quit R.

> postscript()
> 
> source("demo.R",echo=TRUE)

> library(bqtl)

> cvl.map <- make.map.frame(read.table("../plant/cvl-data/cvl.map", 
    header = TRUE))

> plot(cvl.map)

> cvl.markers <- read.csv("../plant/cvl-data/cvl.markers")

> cvl.codes <- apply(cvl.markers, 2, function(x) ifelse(x == 
    "AA", 1, ifelse(x == "aa", 2, 3)))

> image(1:161, 1:163, cvl.codes, xlab = "RIL", ylab = "marker")

> abline(v = seq(1.5, by = 1, length = 160))

> abline(h = 0.5 + which(cvl.map$pos.type == "right"))

> increasing.cM <- function(x, extra = 1) {
    add.to.x <- cumsum(extra + c(0, x$cM)[which(x$pos.type == 
        "left")])
    add.to.x[x$chr.num] + x .... [TRUNCATED] 

> x.cM <- increasing.cM(cvl.map)

> x.cM <- c(x.cM, x.cM[163] + 1)

> image(1:161, x.cM, cvl.codes, xlab = "RIL", ylab = "cM")

> abline(v = seq(1.5, by = 1, length = 160))

> abline(h = x.cM[which(cvl.map$pos.type == "right")[1:4] + 
    1])

> pheno.dat <- read.csv("../plant/cvl-data/pheno.dat")

> summary(pheno.dat)
      blue            brz              dark           farred     
 Min.   : 3.49   Min.   : 4.800   Min.   :10.38   Min.   :2.700  
 1st Qu.: 5.66   1st Qu.: 6.370   1st Qu.:14.52   1st Qu.:4.370  
 Median : 6.46   Median : 7.520   Median :15.73   Median :4.870  
 Mean   : 6.66   Mean   : 7.662   Mean   :15.68   Mean   :4.975  
 3rd Qu.: 7.58   3rd Qu.: 8.820   3rd Qu.:17.01   3rd Qu.:5.430  
 Max.   :10.27   Max.   :12.640   Max.   :20.29   Max.   :7.760  
       ga              red             white             germ    
 Min.   : 4.510   Min.   : 5.150   Min.   :3.370   good    :148  
 1st Qu.: 6.730   1st Qu.: 7.590   1st Qu.:4.870   not.good: 13  
 Median : 7.700   Median : 8.750   Median :5.650                 
 Mean   : 7.882   Mean   : 8.886   Mean   :5.783                 
 3rd Qu.: 8.860   3rd Qu.:10.100   3rd Qu.:6.650                 
 Max.   :12.200   Max.   :14.080   Max.   :8.510                 

> pairs(pheno.dat)

> round(cor(data.matrix(pheno.dat)), 2)
        blue   brz  dark farred    ga   red white  germ
blue    1.00  0.64  0.68   0.73  0.76  0.70  0.74 -0.28
brz     0.64  1.00  0.62   0.72  0.58  0.60  0.49 -0.17
dark    0.68  0.62  1.00   0.61  0.62  0.51  0.56 -0.17
farred  0.73  0.72  0.61   1.00  0.54  0.70  0.49 -0.25
ga      0.76  0.58  0.62   0.54  1.00  0.76  0.89 -0.17
red     0.70  0.60  0.51   0.70  0.76  1.00  0.73 -0.19
white   0.74  0.49  0.56   0.49  0.89  0.73  1.00 -0.18
germ   -0.28 -0.17 -0.17  -0.25 -0.17 -0.19 -0.18  1.00

> image(1:161, 1:163, cvl.codes[order(pheno.dat$white), 
    ], xlab = "RIL", ylab = "marker")

> abline(v = seq(1.5, by = 1, length = 160))

> abline(h = 0.5 + which(cvl.map$pos.type == "right"))

> title("short             in white light left             long ")

> cvl.ana <- make.analysis.obj(pheno.dat, make.map.frame(cvl.map, 
    reso = 2), cvl.markers, method = "RI.self")

> plot(cvl.ana)

> image(1:161, 1:163, cvl.ana$state.matrix[, cvl.ana$map.frame$is.marker, 
    2], xlab = "RIL", ylab = "locus", zlim = c(0, 2))

> abline(v = seq(1.5, by = 1, length = 160))

> abline(h = 0.5 + which(cvl.map$pos.type == "right"))

> image(1:163, 1:163, cor(cvl.ana$state.matrix[, cvl.ana$map.frame$is.marker, 
    2]))

> contour(1:163, 1:163, cor(cvl.ana$state.matrix[, cvl.ana$map.frame$is.marker, 
    2]), levels = c(0.2, 0.4, 0.6, 0.8), add = T)

> abline(v = 0.5 + which(cvl.map$pos.type == "right"))

> abline(h = 0.5 + which(cvl.map$pos.type == "right"))

> x.cM <- increasing.cM(cvl.map)

> div.pt <- c(x.cM, x.cM[163] + 1)

> mid.pt <- (x.cM + div.pt[-1])/2

> marker.cors <- cor(cvl.ana$state.matrix[, cvl.ana$map.frame$is.marker, 
    2])

> image(div.pt, div.pt, marker.cors)

> contour(mid.pt, mid.pt, marker.cors, levels = c(0.2, 
    0.4, 0.6, 0.8), add = T)

> abline(v = div.pt[1 + which(cvl.map$pos.type == "right")])

> abline(h = div.pt[1 + which(cvl.map$pos.type == "right")])

> fit.null <- bqtl(white ~ 1, cvl.ana)

> summary(fit.null)
$coefficients
Intercept 
 5.782733 

$loglik
[1] -256.0794

$std.res
[1] 1.187223

$N
NULL


> fit.PVV4 <- bqtl(white ~ PVV4, cvl.ana)

> summary(fit.PVV4)
$coefficients
Intercept      PVV4 
5.7503866 0.4044548 

$loglik
[1] -246.3247

$std.res
[1] 1.116674

$N
     N N.omit N.used 
   161      0    161 


> scan.1 <- bqtl(white ~ locus(all), cvl.ana)

> lod.ratio <- log10(exp(1))

> lod.0 <- loglik(fit.null) * lod.ratio

> lod.1 <- loglik(scan.1) * lod.ratio

> par(mfrow = c(3, 2))

> plot(cvl.ana, lod.1 - lod.0)

> white.perm <- sample(1:161)

> perm.1 <- bqtl(white[white.perm] ~ locus(all), cvl.ana)

> scan.2 <- bqtl(white ~ locus(chromo = 1, cM = c(0, 
    30)) + locus(chromo = 1, cM = c(100, 150)), cvl.ana)

> scan.2.lod <- loglik(scan.2) * lod.ratio

> index.2 <- map.loc(scan.2)

> vec.2 <- unlist(index.2)

> ind.1 <- as.numeric(vec.2[names(vec.2) == "cM1"])

> ind.2 <- as.numeric(vec.2[names(vec.2) == "cM2"])

> image(unique(ind.1), unique(ind.2), matrix(scan.2.lod, 
    nr = length(unique(ind.1))), xlab = "cM", ylab = "cM")

> text(ind.1, ind.2, round(scan.2.lod - min(scan.2.lod), 
    1), cex = 0.6)

> lb.spec <- list(gene.number = 1:10, n.cycles = c(0, 
    400, rep(100, 8)))

> white.lb <- linear.bayes(white ~ locus(all), cvl.ana, 
    rp = 1, spec = lb.spec)

> par(mfrow = c(1, 1))

> plot(log10(white.lb$odds), type = "h")

> par(mfrow = c(3, 2))

> plot(cvl.ana, white.lb$loc.posterior, type = "h")
> 
> 
> dev.off()
null device 
          1 
>