###################################################
### chunk number 1: setup
###################################################
#line 60 "vignettes/vsn/inst/doc/vsn.Rnw"
options(width=41, signif=3, digits=3)
set.seed(0xdada)

## To create bitmap versions of plots with many dots, circumventing
##   Sweave's fig=TRUE mechanism...
##   (pdfs are too large)
openBitmap = function(nm, rows=1, cols=1) {
  png(paste("vsn-", nm, ".png", sep=""), 
       width=600*cols, height=700*rows, pointsize=14)
  par(mfrow=c(rows, cols), cex=2)
}


###################################################
### chunk number 2: overv1
###################################################
#line 80 "vignettes/vsn/inst/doc/vsn.Rnw"
require("vsn")
data("kidney")
xnorm = justvsn(kidney)


###################################################
### chunk number 3: overv2
###################################################
#line 91 "vignettes/vsn/inst/doc/vsn.Rnw"
M = exprs(xnorm)[,1] - exprs(xnorm)[,2]


###################################################
### chunk number 4: overv3
###################################################
#line 101 "vignettes/vsn/inst/doc/vsn.Rnw"
fit = vsn2(kidney)
ynorm = predict(fit, kidney)


###################################################
### chunk number 5: overv4
###################################################
#line 105 "vignettes/vsn/inst/doc/vsn.Rnw"
stopifnot(
  identical(exprs(xnorm), exprs(ynorm)), 
  identical(exprs(xnorm), exprs(fit)))


###################################################
### chunk number 6: nkid-scp1
###################################################
#line 207 "vignettes/vsn/inst/doc/vsn.Rnw"
openBitmap("nkid-scp", cols=2)


###################################################
### chunk number 7: nkid-scp2
###################################################
#line 210 "vignettes/vsn/inst/doc/vsn.Rnw"
select = (0==rowSums(exprs(kidney)<=0))
plot(log2(exprs(kidney)[select, ]), 
  main = "a) raw", pch = ".", asp=1)
plot(exprs(xnorm), main = "b) vsn", 
  pch = ".", asp=1,
  col=ifelse(select, "black", "orange"))


###################################################
### chunk number 8: nkid-scp3
###################################################
#line 218 "vignettes/vsn/inst/doc/vsn.Rnw"
dev.off()


###################################################
### chunk number 9: nkid-sdmean1
###################################################
#line 233 "vignettes/vsn/inst/doc/vsn.Rnw"
openBitmap("nkid-sdmean", cols=2)


###################################################
### chunk number 10: nkid-sdmean2
###################################################
#line 236 "vignettes/vsn/inst/doc/vsn.Rnw"
meanSdPlot(xnorm, ranks=TRUE)
meanSdPlot(xnorm, ranks=FALSE)


###################################################
### chunk number 11: nkid-sdmean3
###################################################
#line 240 "vignettes/vsn/inst/doc/vsn.Rnw"
dev.off()


###################################################
### chunk number 12: nkid-histM
###################################################
#line 270 "vignettes/vsn/inst/doc/vsn.Rnw"
hist(M, breaks=33, col="#d95f0e")


###################################################
### chunk number 13: lymphoma
###################################################
#line 292 "vignettes/vsn/inst/doc/vsn.Rnw"
data("lymphoma")
dim(lymphoma)


###################################################
### chunk number 14: pDatalym
###################################################
#line 299 "vignettes/vsn/inst/doc/vsn.Rnw"
pData(lymphoma)


###################################################
### chunk number 15: lymjustvsn
###################################################
#line 314 "vignettes/vsn/inst/doc/vsn.Rnw"
lym = justvsn(lymphoma)


###################################################
### chunk number 16: lym-sdmean1
###################################################
#line 317 "vignettes/vsn/inst/doc/vsn.Rnw"
openBitmap("lym-sdmean")


###################################################
### chunk number 17: lym-sdmean2
###################################################
#line 320 "vignettes/vsn/inst/doc/vsn.Rnw"
meanSdPlot(lym)


###################################################
### chunk number 18: lym-sdmean3
###################################################
#line 323 "vignettes/vsn/inst/doc/vsn.Rnw"
dev.off()


###################################################
### chunk number 19: lym-M
###################################################
#line 332 "vignettes/vsn/inst/doc/vsn.Rnw"
iref = seq(1, 15, by=2)
ismp = seq(2, 16, by=2)
M= exprs(lym)[,ismp]-exprs(lym)[,iref] 
A=(exprs(lym)[,ismp]+exprs(lym)[,iref])/2
colnames(M) = lymphoma$sample[ismp]
colnames(A) = colnames(M)

j = "DLCL-0032"
smoothScatter(A[,j], M[,j], main=j, 
        xlab="A", ylab="M", pch=".")
abline(h=0, col="red")


###################################################
### chunk number 20: affy1
###################################################
#line 367 "vignettes/vsn/inst/doc/vsn.Rnw"
require("affydata")
data("Dilution")
d_vsn = vsnrma(Dilution)


###################################################
### chunk number 21: affy2
###################################################
#line 373 "vignettes/vsn/inst/doc/vsn.Rnw"
d_rma = rma(Dilution)


###################################################
### chunk number 22: figaffy1
###################################################
#line 379 "vignettes/vsn/inst/doc/vsn.Rnw"
openBitmap("affy", cols=3, rows=1)


###################################################
### chunk number 23: figaffy2
###################################################
#line 382 "vignettes/vsn/inst/doc/vsn.Rnw"
par(pch=".")
ax = c(2, 16)

plot(exprs(d_vsn)[,c(1,3)], 
  main = "vsn: array 1 vs 3",
  asp=1, xlim=ax, ylim=ax)

plot(exprs(d_rma)[,c(1,3)], 
  main = "rma: array 1 vs 3",
  asp=1, xlim=ax, ylim=ax)

plot(exprs(d_rma)[,1], 
     exprs(d_vsn)[,1], 
     xlab="rma", ylab="vsn", 
     asp=1, xlim=ax, ylim=ax,
     main = "array 1")
abline(a=0, b=1, col="#ff0000d0")


###################################################
### chunk number 24: figaffy3
###################################################
#line 401 "vignettes/vsn/inst/doc/vsn.Rnw"
dev.off()


###################################################
### chunk number 25: limma
###################################################
#line 428 "vignettes/vsn/inst/doc/vsn.Rnw"
require("limma")
wg = which(lymphoma$dye=="Cy3")
wr = which(lymphoma$dye=="Cy5")

lymRG = new("RGList", list(
  R=exprs(lymphoma)[, wr],
  G=exprs(lymphoma)[, wg]))

lymNCS = justvsn(lymRG)


###################################################
### chunk number 26: addmeta
###################################################
#line 468 "vignettes/vsn/inst/doc/vsn.Rnw"
vmd = data.frame(
  labelDescription=I(c("array ID", 
      "sample in G", "sample in R")),
  channel=c("_ALL", "G", "R"),
  row.names=c("arrayID", "sampG", "sampR"))

arrayID = lymphoma$name[wr]
stopifnot(identical(arrayID, 
          lymphoma$name[wg])) 

## remove sample number suffix 
sampleType = factor(sub("-.*", "", 
                   lymphoma$sample)) 

v = data.frame(
  arrayID = arrayID,
  sampG   = sampleType[wg],
  sampR   = sampleType[wr])
v
              
adf = new("AnnotatedDataFrame", 
  data=v,
  varMetadata=vmd)

phenoData(lymNCS) = adf 


###################################################
### chunk number 27: width1
###################################################
#line 498 "vignettes/vsn/inst/doc/vsn.Rnw"
oo=options(width=75L)


###################################################
### chunk number 28: lymNCS
###################################################
#line 501 "vignettes/vsn/inst/doc/vsn.Rnw"
lymNCS


###################################################
### chunk number 29: width2
###################################################
#line 504 "vignettes/vsn/inst/doc/vsn.Rnw"
options(oo)


###################################################
### chunk number 30: lymM
###################################################
#line 519 "vignettes/vsn/inst/doc/vsn.Rnw"
lymM = (assayData(lymNCS)$R - 
        assayData(lymNCS)$G)


###################################################
### chunk number 31: design
###################################################
#line 526 "vignettes/vsn/inst/doc/vsn.Rnw"
design = model.matrix( ~ lymNCS$sampR)
lf = lmFit(lymM, design[, 2, drop=FALSE])
lf = eBayes(lf)


###################################################
### chunk number 32: figlimma
###################################################
#line 536 "vignettes/vsn/inst/doc/vsn.Rnw"
par(mfrow=c(1,2))
hist(lf$p.value, 100, col="orange") 
pdat=t(lymM[order(lf$p.value)[1:5],])
matplot(pdat, 
  lty=1, type="b", lwd=2, 
  col=hsv(seq(0,1,length=5), 0.7, 0.8), 
  ylab="M", xlab="arrays")


###################################################
### chunk number 33: makebg
###################################################
#line 574 "vignettes/vsn/inst/doc/vsn.Rnw"
rndbg=function(x, off, fac)
 array(off+fac*runif(prod(dim(x))),
       dim=dim(x))

lymRGwbg = lymRG
lymRGwbg$Rb = rndbg(lymRG, 100, 30)
lymRGwbg$Gb = rndbg(lymRG,  50, 20)


###################################################
### chunk number 34: justvsnwbg
###################################################
#line 589 "vignettes/vsn/inst/doc/vsn.Rnw"
lymESwbg = justvsn(lymRGwbg[, 1:3], 
   backgroundsubtract=TRUE)


###################################################
### chunk number 35: pinId
###################################################
#line 622 "vignettes/vsn/inst/doc/vsn.Rnw"
ngr = ngc = 4L
nsr = nsc = 24L
arrayGeometry = data.frame(
  spotcol = rep(1:nsc, 
    times = nsr*ngr*ngc),
  spotrow = rep(1:nsr, 
    each = nsc, times=ngr*ngc),
  pin = rep(1:(ngr*ngc), 
    each = nsr*nsc))


###################################################
### chunk number 36: strata
###################################################
#line 636 "vignettes/vsn/inst/doc/vsn.Rnw"
EconStr = justvsn(lymRG[,1], 
     strata=arrayGeometry$pin)


###################################################
### chunk number 37: nostrata
###################################################
#line 644 "vignettes/vsn/inst/doc/vsn.Rnw"
EsenzaStr = justvsn(lymRG[,1])


###################################################
### chunk number 38: figstrata1
###################################################
#line 651 "vignettes/vsn/inst/doc/vsn.Rnw"
openBitmap("figstrata")


###################################################
### chunk number 39: figstrata2
###################################################
#line 654 "vignettes/vsn/inst/doc/vsn.Rnw"
j = 1L
plot(assayData(EsenzaStr)$R[,j],
     assayData(EconStr)$R[,j],
     pch = ".", asp = 1, 
     col = hsv(seq(0, 1, length=ngr*ngc), 
       0.8, 0.6)[arrayGeometry$pin], 
     xlab = "without strata", 
     ylab = "print-tip strata",
     main = sampleNames(lymNCS)$R[j])


###################################################
### chunk number 40: figstrata3
###################################################
#line 665 "vignettes/vsn/inst/doc/vsn.Rnw"
dev.off()


###################################################
### chunk number 41: miss1
###################################################
#line 687 "vignettes/vsn/inst/doc/vsn.Rnw"
lym2 = lymphoma
nfeat = prod(dim(lym2))
wh = sample(nfeat, nfeat/10)
exprs(lym2)[wh] = NA
table(is.na(exprs(lym2)))


###################################################
### chunk number 42: miss2
###################################################
#line 695 "vignettes/vsn/inst/doc/vsn.Rnw"
fit1 = vsn2(lymphoma, lts.quantile=1)
fit2 = vsn2(lym2, lts.quantile=1)


###################################################
### chunk number 43: figmiss
###################################################
#line 703 "vignettes/vsn/inst/doc/vsn.Rnw"
par(mfrow=c(1,2))
for(j in 1:2){
  p1 = coef(fit1)[,,j]
  p2 = coef(fit2)[,,j]
  d  = max(abs(p1-p2))
  stopifnot(d < c(0.05, 0.03)[j])
  plot(p1, p2, pch = 16, asp = 1,
    main = paste(letters[j], 
      ": max diff=", signif(d,2), sep = ""),
    xlab = "no missing data",
    ylab = "10% of data missing")
  abline(a = 0, b = 1, col = "blue")
}


###################################################
### chunk number 44: spikein
###################################################
#line 750 "vignettes/vsn/inst/doc/vsn.Rnw"
spikeins = 100:200
spfit = vsn2(kidney[spikeins,], 
             lts.quantile=1)
nkid = predict(spfit, newdata=kidney)


###################################################
### chunk number 45: ref1
###################################################
#line 808 "vignettes/vsn/inst/doc/vsn.Rnw"
ref = vsn2(lymphoma[, ismp[1:7]])


###################################################
### chunk number 46: ref2
###################################################
#line 814 "vignettes/vsn/inst/doc/vsn.Rnw"
f8 = vsn2(lymphoma[, ismp[8]], 
          reference = ref)


###################################################
### chunk number 47: ref3
###################################################
#line 822 "vignettes/vsn/inst/doc/vsn.Rnw"
fall = vsn2(lymphoma[, ismp])


###################################################
### chunk number 48: ref4
###################################################
#line 824 "vignettes/vsn/inst/doc/vsn.Rnw"
coefficients(f8)[,1,]
coefficients(fall)[,8,]


###################################################
### chunk number 49: figref1
###################################################
#line 832 "vignettes/vsn/inst/doc/vsn.Rnw"
openBitmap("figref")


###################################################
### chunk number 50: figref2
###################################################
#line 835 "vignettes/vsn/inst/doc/vsn.Rnw"
plot(exprs(f8), exprs(fall)[,8], 
     pch=".", asp=1)
abline(a=0, b=1, col="red")


###################################################
### chunk number 51: figref3
###################################################
#line 840 "vignettes/vsn/inst/doc/vsn.Rnw"
dev.off()


###################################################
### chunk number 52: hiddenchecks
###################################################
#line 844 "vignettes/vsn/inst/doc/vsn.Rnw"
stopifnot(length(ismp)==8L)
maxdiff = max(abs(exprs(f8) - exprs(fall)[,8])) 
if(maxdiff>0.3)
  stop(sprintf("maxdiff is %g", maxdiff))


###################################################
### chunk number 53: nkid-calib1
###################################################
#line 903 "vignettes/vsn/inst/doc/vsn.Rnw"
coef(fit)[1,,]


###################################################
### chunk number 54: nkid-calib2
###################################################
#line 997 "vignettes/vsn/inst/doc/vsn.Rnw"
bkid = kidney
exprs(bkid)[,1]=0.25*(500+exprs(bkid)[,1])


###################################################
### chunk number 55: nkid-calib3
###################################################
#line 1004 "vignettes/vsn/inst/doc/vsn.Rnw"
bfit = vsn2(bkid)


###################################################
### chunk number 56: nkid-calib4
###################################################
#line 1008 "vignettes/vsn/inst/doc/vsn.Rnw"
oldwarn = options(warn=-1)
openBitmap("nkid-calib", cols=2)


###################################################
### chunk number 57: nkid-calib5
###################################################
#line 1012 "vignettes/vsn/inst/doc/vsn.Rnw"
plot(exprs(bkid), main="raw", 
     pch=".", log="xy")
plot(exprs(bfit), main="vsn", 
     pch=".")
coef(bfit)[1,,]


###################################################
### chunk number 58: nkid-calib6
###################################################
#line 1019 "vignettes/vsn/inst/doc/vsn.Rnw"
dev.off()
options(oldwarn)
rm(list="oldwarn")


###################################################
### chunk number 59: vsnQ
###################################################
#line 1042 "vignettes/vsn/inst/doc/vsn.Rnw"
lym_q = normalizeQuantiles(exprs(lymphoma))
lym_qvsn = vsn2(lym_q, calib="none")


###################################################
### chunk number 60: figqvsn1
###################################################
#line 1046 "vignettes/vsn/inst/doc/vsn.Rnw"
openBitmap("qvsn", cols=2, rows=1)


###################################################
### chunk number 61: figqvsn
###################################################
#line 1049 "vignettes/vsn/inst/doc/vsn.Rnw"
plot(exprs(lym_qvsn)[, 1:2], pch=".", 
   main="lym_qvsn")
plot(exprs(lym)[,1], exprs(lym_qvsn)[, 1], 
  main="lym_qvsn vs lym", pch=".", 
  ylab="lym_qvsn[,1]", xlab="lym[,1]")


###################################################
### chunk number 62: figqvsn3
###################################################
#line 1056 "vignettes/vsn/inst/doc/vsn.Rnw"
dev.off()


###################################################
### chunk number 63: calcshrink
###################################################
#line 1130 "vignettes/vsn/inst/doc/vsn.Rnw"
log2.na = function(x){
  w = which(x>0)
  res = rep(as.numeric(NA), length(x))
  res[w] = log2(x[w])
  return(res)
}
  
fc = 0.5                    ## true fold change
x2 = seq(0.5, 15, by=0.5)   ## 'true values' in sample 1
x1 = x2/fc                  ## 'true values' in sample 2
m = s = numeric(length(x1))
n  = 10000
sa = 1
b  = 1
sb = 0.1
sdeta = 0.1
for(i in 1:length(x1)){
  z1 = sa*rnorm(n)+x1[i]*b*exp(sb*rnorm(n))
  z2 = sa*rnorm(n)+x2[i]*b*exp(sb*rnorm(n))
  if(i%%2==1) {
    q = log2.na(z1/z2)
    m[i] = mean(q, na.rm=TRUE)
    s[i] = sd(q, na.rm=TRUE)
  } else {
    h = (asinh(z1/(sa*b/sb))-asinh(z2/(sa*b/sb)))/log(2)
    m[i] = mean(h)
    s[i] = sd(h)
  }
}

colq = c("black", "blue")
lty  = 1
pch  = c(20,18)
cex  = 1.4
lwd  = 2


###################################################
### chunk number 64: figshrink
###################################################
#line 1168 "vignettes/vsn/inst/doc/vsn.Rnw"
mai=par("mai"); mai[3]=0; par(mai)
plot(x2, m, ylim=c(-2, 3.5), type="n", xlab=expression(x[2]), ylab="")
abline(h=-log2(fc), col="red", lwd=lwd, lty=1)
abline(h=0, col="black", lwd=1, lty=2)
points(x2, m, pch=pch, col=colq, cex=cex)
segments(x2, m-2*s, x2, m+2*s, lwd=lwd, col=colq, lty=lty)
legend(8.75, -0.1, c("q","h"), col=colq, pch=pch, lty=lty, lwd=lwd)


###################################################
### chunk number 65: figgraph
###################################################
#line 1179 "vignettes/vsn/inst/doc/vsn.Rnw"
par(mai = c(0.8, 0.6, 0.01, 0.01))
x = seq(-70.5, 170.5, by=1)
cols = c("black", "blue", "darkgrey")
xoff = cc = 50
ymat = cbind(log2.na(x), log2( (x+sqrt(x^2+cc^2))/2 ), log2.na(x+xoff))
ylim = range(ymat, na.rm=TRUE)
matplot(x, ymat, 
        lty=c(1,1,2), col=cols, lwd=3, type="l", ylim=ylim, ylab="")
abline(v = 0, col = "#80808080", lty = 2)
legend(x = par("usr")[2], y = par("usr")[3], 
       legend = c(expression(y = log[2](x)), expression(y = glog[2](x,c)), expression(y = log[2](x+x[off]))), 
       fill=cols, density = c(NA, NA, 50), xjust=1.1, yjust=-0.1)


###################################################
### chunk number 66: lymbox
###################################################
#line 1282 "vignettes/vsn/inst/doc/vsn.Rnw"
colours = hsv(seq(0,1,length=nsr),0.6,1)
j = "CLL-13"
boxplot(A[, j] ~ arrayGeometry$spotrow, 
   col=colours, main=j, 
   ylab="A", xlab="spotrow")


###################################################
### chunk number 67: lymquscp1
###################################################
#line 1311 "vignettes/vsn/inst/doc/vsn.Rnw"
openBitmap("lymquscp")


###################################################
### chunk number 68: lymquscp2
###################################################
#line 1314 "vignettes/vsn/inst/doc/vsn.Rnw"
plot(A[,j], M[,j], pch=16, cex=0.3,
col=ifelse(arrayGeometry$spotrow%in%(22:23),
           "orange", "black"))
abline(h=0, col="blue")


###################################################
### chunk number 69: lymquscp3
###################################################
#line 1320 "vignettes/vsn/inst/doc/vsn.Rnw"
dev.off()


###################################################
### chunk number 70: sessionInfo
###################################################
#line 1462 "vignettes/vsn/inst/doc/vsn.Rnw"
toLatex(sessionInfo())