### R code from vignette source 'Hiiragi2013.Rnw'

###################################################
### code chunk number 1: style
###################################################
BiocStyle::latex(short.fignames=TRUE, use.unsrturl=FALSE)


###################################################
### code chunk number 2: packages
###################################################
library("Hiiragi2013")
set.seed(2013)


###################################################
### code chunk number 3: loadx
###################################################
data("x")
x


###################################################
### code chunk number 4: showWT
###################################################
with(subset(pData(x), genotype=="WT"),
     addmargins(table(Embryonic.day, Total.number.of.cells), 2))


###################################################
### code chunk number 5: showFGF4KO
###################################################
with(subset(pData(x), genotype=="FGF4-KO"), table(Embryonic.day))


###################################################
### code chunk number 6: sampleColours
###################################################
groups = split(seq_len(ncol(x)), pData(x)$sampleGroup)
sapply(groups, length)


###################################################
### code chunk number 7: sampleColours
###################################################
sampleColourMap = setNames(unique(pData(x)$sampleColour), unique(pData(x)$sampleGroup))
sampleColourMap


###################################################
### code chunk number 8: FGF4probes
###################################################
FGF4probes = (fData(x)$symbol == "Fgf4")
stopifnot(sum(FGF4probes)==4)


###################################################
### code chunk number 9: expressed1
###################################################
selectedSamples = with(pData(x), genotype=="WT") 
xe = x[, selectedSamples]
stopifnot(ncol(xe)==66)


###################################################
### code chunk number 10: figHistSds
###################################################
sdxe = rowSds(exprs(xe))
thresh = 0.5
hist(sdxe, 100, col = "skyblue")
abline(v = thresh)


###################################################
### code chunk number 11: anaHistSds
###################################################
table(sdxe>=thresh)
length(unique(fData(xe)$ensembl[ sdxe>=thresh ]))


###################################################
### code chunk number 12: mas5calls
###################################################
data("a")
stopifnot(nrow(pData(a))==ncol(x))
mas5c = mas5calls(a[, selectedSamples])


###################################################
### code chunk number 13: numCalls1
###################################################
myUnique = function(x) setdiff(unique(x), "")
allEnsemblIDs = myUnique(fData(xe)$ensembl)

callsPerGenePerArray = matrix(0, nrow = length(allEnsemblIDs), ncol = ncol(mas5c)+1,
                              dimnames = list(allEnsemblIDs, NULL))

for(j in seq_len(ncol(mas5c))) {
  for(k in 1:2) {
    ids = myUnique(fData(xe)$ensembl[ exprs(mas5c)[, j]==c("M","P")[k] ])
    callsPerGenePerArray[ids, j] = k
  }
}
 
fractionOfArrays = 0.1
for(k in 1:2) {
  ids = myUnique(fData(xe)$ensembl[ apply(exprs(mas5c)==c("M","P")[k], 1, 
                                          function(v) (mean(v)>fractionOfArrays)) ])
  callsPerGenePerArray[ids, ncol(mas5c)+1] = k
}

numCalls = apply(callsPerGenePerArray, 2, table)
numCalls = numCalls[rev(seq_len(nrow(numCalls))), ]


###################################################
### code chunk number 14: numCalls2
###################################################
numCalls[, 67]


###################################################
### code chunk number 15: figmas5calls
###################################################
barplot2(numCalls, names.arg = paste(seq_len(ncol(callsPerGenePerArray))),
         col = c(brewer.pal(8, "Paired")[2:1], "#e8e8e8"), ylab = "number of genes")


###################################################
### code chunk number 16: clusterResampling
###################################################
clusterResampling = function(x, ngenes, k = 2, B = 250, prob = 0.67) {
  mat = exprs(x)
  ce = cl_ensemble(list = lapply(seq_len(B), function(b) {
    selSamps = sample(ncol(mat), size = round(prob*ncol(mat)), replace = FALSE)
    submat = mat[, selSamps, drop = FALSE]
    selFeats = order(rowVars(submat), decreasing = TRUE)[seq_len(ngenes)]
    submat = submat[selFeats,, drop = FALSE]
    pamres = pam(t(submat), k = k, metric = "euclidean")
    pred = cl_predict(pamres, t(mat[selFeats, ]), "memberships")
    as.cl_partition(pred)
  }))
  cons = cl_consensus(ce)
  ag = sapply(ce, cl_agreement, y = cons)
  return(list(agreements = ag, consensus = cons))
}


###################################################
### code chunk number 17: ce1
###################################################
ce = list(
  "E3.25" = clusterResampling(x[, unlist(groups[c("E3.25")])], ngenes = 20),
  "E3.5"  = clusterResampling(x[, unlist(groups[c("E3.5 (EPI)", "E3.5 (PE)")])],
                              ngenes = 20))


###################################################
### code chunk number 18: figClue1
###################################################
par(mfrow = c(1,2))
colours  = c(sampleColourMap["E3.25"], brewer.pal(9,"Set1")[9])

boxplot(lapply(ce, `[[`, "agreements"), ylab = "agreement probabilities", col = colours)

mems = lapply(ce, function(x) sort(cl_membership(x$consensus)[, 1]))
mgrp = lapply(seq(along = mems), function(i) rep(i, times = length(mems[[i]])))
myjitter = function(x) x+seq(-.4, +.4, length.out = length(x))

plot(unlist(lapply(mgrp, myjitter)), unlist(mems),
     col = colours[unlist(mgrp)], ylab = "membership probabilities",
     xlab = "consensus clustering", xaxt = "n", pch = 16)
text(x = 1:2, y = par("usr")[3], labels = c("E3.25","E3.5"), adj = c(0.5, 1.4), xpd = NA)


###################################################
### code chunk number 19: test
###################################################
wilcox.test(ce$E3.25$agreements, ce$E3.5$agreements)


###################################################
### code chunk number 20: ce2
###################################################
sampleSets = list(
  `E3.5` = unlist(groups[c("E3.5 (EPI)", "E3.5 (PE)", "E4.5 (FGF4-KO)")]),
  `E4.5` = unlist(groups[c("E4.5 (EPI)", "E4.5 (PE)", "E4.5 (FGF4-KO)")]))
k = 3
csa = lapply(sampleSets, function(samps) list(
  colours = x$sampleColour[samps], 
  r = clusterResampling(x[!FGF4probes, samps], ngenes = 20, k = k)))


###################################################
### code chunk number 21: figClue2
###################################################
par(mfrow = c(2,3))
for(i in seq(along = csa))
  for(j in seq_len(k))
    plot(cl_membership(csa[[i]]$r$consensus)[, j], 
         ylab = paste0("P(cluster ", j, ")"), xlab = "samples",
         main = names(sampleSets)[i], col = csa[[i]]$colours, pch = 16, cex = 1.5)


###################################################
### code chunk number 22: ce3
###################################################
sampleSets = unlist(groups[c("E3.25", "E3.5 (FGF4-KO)")])
selSamples = x[!FGF4probes, sampleSets]
resampledSampleSet = clusterResampling(selSamples, ngenes = 20)


###################################################
### code chunk number 23: figClue3
###################################################
par(mfrow = c(1,2))
  for(j in seq_len(2))
    plot(cl_membership(resampledSampleSet$consensus)[, j], 
         ylab = paste0("P(cluster ", j, ")"), xlab = "Samples",
         col = x$sampleColour[sampleSets], pch = 16, cex = 1.5)


###################################################
### code chunk number 24: figHeatmapClue
###################################################
ngenes = 100
selFeats = order(rowVars(exprs(selSamples)), decreasing = TRUE)[seq_len(ngenes)]
myHeatmap(selSamples[selFeats, ], collapseDuplicateFeatures = TRUE, haveColDend = TRUE)


###################################################
### code chunk number 25: doCluster
###################################################
ngenes = c(10, 25, 50, 100, 250, 500, 1000, 2500, 5000, 10000)


###################################################
### code chunk number 26: cluster
###################################################
xForClustering = x[, x$Embryonic.day=="E3.5" & x$genotype=="WT"]
clusters = sapply(ngenes, pamCluster, x = xForClustering)


###################################################
### code chunk number 27: figClusters
###################################################
image(x = seq_len(nrow(clusters)), y = seq_len(ncol(clusters)), z = clusters,
      col = c("#f0f0f0", "#000000"), ylab = "ngenes", xlab = "samples", yaxt = "n")
text(x = 0, y = seq_len(ncol(clusters)), paste(ngenes), xpd = NA, adj = c(1, 0.5))


###################################################
### code chunk number 28: choosengenes
###################################################
i = which(ngenes==1000); stopifnot(length(i)==1)
ngenes = ngenes[i]
clusters = factor(clusters[, i])


###################################################
### code chunk number 29: compareClusteringAnnotation
###################################################
table(clusters, pData(x)[names(clusters), "lineage"])


###################################################
### code chunk number 30: compareClusteringAnnotationCheck
###################################################
stopifnot(all(
 table(clusters, pData(x)[names(clusters), "lineage"])==cbind(c(0,11),c(11,0))
))


###################################################
### code chunk number 31: print
###################################################
cbind(unlist(groups[c("E3.5 (EPI)", "E3.5 (PE)")]), ce[[2]]$consensus$.Data[, 1])


###################################################
### code chunk number 32: deBetweenClusters
###################################################
deCluster = rowttests(xForClustering, fac = clusters)


###################################################
### code chunk number 33: figIndepFiltSetup
###################################################
varianceRank = rank(-rowVars(exprs(xForClustering)))
plot(varianceRank, deCluster$p.value, pch = ".", log = "y",
     main = "Parameters for the independent filtering",
     xlab = "variance rank", ylab = "p-value")
nFilt = 20000
smallpValue = 1e-4
abline(v = nFilt, col = "blue")
abline(h = smallpValue, col = "orange")


###################################################
### code chunk number 34: checkFilterClaim
###################################################
stopifnot(!any((deCluster$p.value<smallpValue)&(varianceRank>nFilt)))


###################################################
### code chunk number 35: BH
###################################################
passfilter = which(varianceRank<=nFilt)
adjp = rep(NA_real_, nrow(x))
adjp[passfilter] = p.adjust(deCluster$p.value[passfilter], method = "BH")
ord = order(adjp)
numFeaturesReport = 200
differentially = ord[seq_len(numFeaturesReport)]
length(unique(fData(x)$symbol[differentially]))


###################################################
### code chunk number 36: write.csv
###################################################
deFeat35 = cbind(deCluster[differentially,], `FDR-adjusted p-value` = adjp[differentially], 
                 fData(x)[differentially,])
write.csv(deFeat35, file = "differentially-expressed-features-3.5.csv")


###################################################
### code chunk number 37: figHeatmapAll
###################################################
myHeatmap(x[differentially, x$genotype=="WT"])


###################################################
### code chunk number 38: figHeatmapAllColl
###################################################
myHeatmap(x[differentially, x$genotype=="WT"], collapseDuplicateFeatures = TRUE)


###################################################
### code chunk number 39: figHeatmap35
###################################################
myHeatmap(xForClustering[differentially, ])


###################################################
### code chunk number 40: figHeatmap35Coll
###################################################
myHeatmap(xForClustering[differentially, ], collapseDuplicateFeatures = TRUE)


###################################################
### code chunk number 41: de45
###################################################
x45 = x[, x$Embryonic.day=="E4.5" & x$genotype=="WT"]
de45 = rowttests(x45, fac = "lineage")


###################################################
### code chunk number 42: figIF45
###################################################
varianceRank = rank(-rowVars(exprs(x45)))
plot(varianceRank, de45$p.value, pch = ".", log = "y",
     main = "Parameters for the independent filtering",
     xlab = "variance rank", ylab = "p-value")
abline(v = nFilt, col = "blue")
abline(h = smallpValue, col = "orange")


###################################################
### code chunk number 43: BH45
###################################################
passfilter = which(varianceRank<=nFilt)
adjp = rep(NA_real_, nrow(x))
adjp[passfilter] = p.adjust(de45$p.value[passfilter], method = "BH")
ord = order(adjp)
differentially = ord[seq_len(numFeaturesReport)]
length(unique(fData(x)$symbol[differentially]))


###################################################
### code chunk number 44: write.csv
###################################################
deFeat45 = cbind(de45[differentially, ], `FDR-adjusted p-value` = adjp[differentially], 
                 fData(x)[differentially, ])
write.csv(deFeat45, file = "differentially-expressed-features-4.5.csv")


###################################################
### code chunk number 45: differentiallyE325toE35
###################################################
samples = unlist(groups[c("E3.25", "E3.5 (EPI)", "E3.5 (PE)")])
deE325toE35 = union(
  getDifferentialExpressedGenes(x, groups, "E3.25", "E3.5 (EPI)"),
  getDifferentialExpressedGenes(x, groups, "E3.25", "E3.5 (PE)"))


###################################################
### code chunk number 46: figdifferentiallyE325toE35
###################################################
myHeatmap(x[deE325toE35, samples], collapseDuplicateFeatures = TRUE)


###################################################
### code chunk number 47: E325samples
###################################################
projection = deCluster$dm[differentially] %*% exprs(x)[differentially, ]


###################################################
### code chunk number 48: figProjection
###################################################
plotProjection(projection, label = sampleNames(x),
               col = x$sampleColour, colourMap = sampleColourMap)


###################################################
### code chunk number 49: safeSelect-1
###################################################
safeSelect = function(grpnames){
  stopifnot(all(grpnames %in% names(groups)))
  unlist(groups[grpnames])
}
g = safeSelect(c("E3.25",
                 "E3.5 (EPI)", "E3.5 (PE)",
                 "E4.5 (EPI)", "E4.5 (PE)"))


###################################################
### code chunk number 50: varianceOrder
###################################################
nfeatures = 100
varianceOrder = order(rowVars(exprs(x[, g])), decreasing = TRUE)
varianceOrder = setdiff(varianceOrder, which(FGF4probes))
selectedFeatures = varianceOrder[seq_len(nfeatures)]
xwt = x[selectedFeatures, g]


###################################################
### code chunk number 51: doPCA1
###################################################
tab = table(xwt$sampleGroup)
sp = split(seq_len(ncol(xwt)), xwt$sampleGroup)
siz = max(listLen(sp))
sp = lapply(sp, sample, size = siz, replace = (siz>length(x)))
xwte = xwt[, unlist(sp)]


###################################################
### code chunk number 52: doPCA2
###################################################
thepca = prcomp(t(exprs(xwte)), center = TRUE)
pcatrsf = function(x) scale(t(exprs(x)), center = TRUE, scale = FALSE) %*% thepca$rotation
stopifnot(all( abs(pcatrsf(xwte) - thepca$x) < 1e-6 ))


###################################################
### code chunk number 53: myPCA
###################################################
myPCAplot = function(x, labels, ...) {
  xy = pcatrsf(x)[, 1:2]
  plot(xy, pch = 16, col = x$sampleColour, cex = 1, xlab = "PC1", ylab = "PC2", ...)
  if(!missing(labels))
    text(xy, labels, cex = 0.35, adj = c(0.5, 0.5))
}


###################################################
### code chunk number 54: figmds-1
###################################################
myPCAplot(xwt)


###################################################
### code chunk number 55: figpcaloadings
###################################################
par(mfrow = c(1,2))
for(v in c("PC1", "PC2")) {
  loading = thepca$rotation[, v]
  plot(sort(loading), main = v, ylab = "")
  sel = order(loading)[c(1:10, (-9:0)+length(loading))]
  print(data.frame(
   symbol = fData(x)$symbol[selectedFeatures[sel]],
   probe = names(loading)[sel],
   loading = loading[sel], stringsAsFactors = FALSE
  ))
}


###################################################
### code chunk number 56: figmds-2
###################################################
myPCAplot(x[selectedFeatures, ])


###################################################
### code chunk number 57: figmds-3
###################################################
myPCAplot(x[selectedFeatures, ], labels = paste(seq_len(ncol(x))))


###################################################
### code chunk number 58: figmds-4
###################################################
myPCAplot(x[selectedFeatures, ], labels = paste(x$Total.number.of.cells))


###################################################
### code chunk number 59: hmprep
###################################################
mat = exprs(x[selectedFeatures, ])
rownames(mat) = fData(x)[selectedFeatures, "symbol"]


###################################################
### code chunk number 60: fighmko
###################################################
heatmap.2(mat, trace = "none", dendrogram = "none", scale = "row",
          col = bluered(100), keysize = 0.9,
          ColSideColors = x$sampleColour, margins = c(7,5))


###################################################
### code chunk number 61: figFGF4expression
###################################################
x325  = x[, with(pData(x), Embryonic.day=="E3.25")]
rv325 = rowVars(exprs(x325))
featureColours = brewer.pal(sum(FGF4probes), "Dark2")
py = t(exprs(x325)[FGF4probes, ])
matplot(py, type = "p", pch = 15, col = featureColours,
        xlab = "arrays", ylab = expression(log[2] ~ signal),
        ylim = range(py) + c(-0.7, 0))
legend("topright", legend = rownames(fData(x325))[FGF4probes], fill = featureColours)
points(seq_len(nrow(py)), rep(par("usr")[3]+0.2, nrow(py)),
       pch = 16, col = x325$sampleColour)


###################################################
### code chunk number 62: checkClaim1
###################################################
stopifnot(all(c("1420085_at", "1420086_x_at") %in% rownames(fData(x325))[FGF4probes]),
          all(x325$genotype[1:36]=="WT"),
          all(x325$genotype[37:44]=="FGF4-KO"))


###################################################
### code chunk number 63: FGF4expression2
###################################################
fgf4Expression = colMeans(exprs(x325)[c("1420085_at", "1420086_x_at"), ])
fgf4Genotype = factor(x325$genotype, 
                      levels = sort(unique(x325$genotype),decreasing = TRUE))


###################################################
### code chunk number 64: figFGF4jitter
###################################################
plot(x = jitter(as.integer(fgf4Genotype)),
     y = fgf4Expression,
     col = x325$sampleColour, xlim = c(0.5, 2.5), pch = 16,
     ylab = expression(FGF4~expression~(log[2]~units)),
     xlab = "genotype", xaxt = "n")
cm = sampleColourMap[sampleColourMap %in% x325$sampleColour]
legend("topright", legend = names(cm), fill = cm)


###################################################
### code chunk number 65: Fgf4MDS
###################################################
zero2one = function(x) (x-min(x))/diff(range(x))
rgb2col = function(x) {x = x/255; rgb(x[, 1], x[, 2], x[, 3])}

colours = x325$sampleColour
wt325 = x325$genotype=="WT"
colourBar = function(x) rgb2col(colorRamp(c("yellow", "blue"))(zero2one(x)))
colours[wt325] = colourBar(fgf4Expression)[wt325]

selMDS = order(rv325, decreasing = TRUE)[seq_len(100)]


###################################################
### code chunk number 66: figFgf4MDS
###################################################
MDSplot(x325[selMDS, ], col = colours)


###################################################
### code chunk number 67: figFgf4MDSColourBar
###################################################
atColour = seq(min(fgf4Expression), max(fgf4Expression), length = 20)
image(z = rbind(seq(along = atColour)), col = colourBar(atColour),
      xaxt = "n", y = atColour, ylab = "")


###################################################
### code chunk number 68: corFgf4Dist
###################################################
KOmean = rowMeans(exprs(x325)[selMDS, x325$genotype=="FGF4-KO"])
dists = colSums((exprs(x325)[selMDS, wt325] - KOmean)^2)^0.5
ct = cor.test(fgf4Expression[wt325], dists, method = "spearman")
ct


###################################################
### code chunk number 69: figCorFgf4Dist
###################################################
plot(fgf4Expression[wt325], dists, pch = 16, main = "E3.25 WT samples",
     xlab = "FGF4 expression", ylab = "Distance to FGF4-KO", col = colours)


###################################################
### code chunk number 70: Hiiragi2013.Rnw:894-895
###################################################
stopifnot(ct$p.value<0.01)


###################################################
### code chunk number 71: FGF4variab1
###################################################
varGroups = split(seq_len(ncol(x)), f = list(x$sampleGroup,x$Total.number.of.cells),
                  sep = ":", drop = TRUE)


###################################################
### code chunk number 72: ScanDate
###################################################
data.frame(
  `number arrays` = listLen(varGroups),
  `ScanDates` = sapply(varGroups, function(v) 
    paste(as.character(unique(x$ScanDate[v])), collapse = ", ")),
  stringsAsFactors = FALSE)


###################################################
### code chunk number 73: FGF4variab1
###################################################
sel   = varianceOrder[seq_len(nfeatures)]
myfun = function(x) median(apply(exprs(x), 1, mad))

sds = lapply(varGroups, function(j) myfun(x[sel, j]))
names(sds) = sprintf("%s (n=%d)", names(sds), listLen(varGroups))

varGroupX = factor(sapply(strsplit(names(varGroups), split = ":"),`[`, 1))


###################################################
### code chunk number 74: figFGF4variab
###################################################
op = par(mai = c(2,0.7,0.1,0.1))
plot(jitter(as.integer(varGroupX)), sds, xaxt = "n", xlab = "", ylab = "")
axis(1, las = 2, tick = FALSE, at = unique(varGroupX), labels = unique(varGroupX))
par(op)


###################################################
### code chunk number 75: sds
###################################################
gps = split(seq_len(ncol(x)), f = x$sampleGroup)[c("E3.25", "E3.25 (FGF4-KO)")]
sds = sapply(gps, function(j) apply(exprs(x)[sel, j], 1, mad))
summary(sds)
apply(sds, 2, function(x) c(`mean` = mean(x), `sd` = sd(x)))


###################################################
### code chunk number 76: numberOfCells
###################################################
for(n in c(100, 1000)) {
  sel = order(rv325, decreasing = TRUE)[seq_len(n)]
  KOmean = rowMeans(exprs(x325)[sel, x325$genotype=="FGF4-KO"])
  dists = colSums((exprs(x325)[sel, wt325] - KOmean)^2)^0.5

  pdf(file = sprintf("Hiiragi2013-figNumberOfCells-%d.pdf", n), width = 5, height = 10)
  par(mfrow = c(2,1))
  plot(x325$Total.number.of.cells[wt325], dists, pch = 16, main = "",
       xlab = "Total number of cells", ylab = "Distance to FGF4-KO")
  MDSplot(x325[sel, ], pointlabel = ifelse(x325$genotype=="WT",
                                           paste(x325$Total.number.of.cells), "KO"), cex = 1)
  dev.off()
}


###################################################
### code chunk number 77: selectedGroups
###################################################
selectedGroups = c("E3.25", "E3.25 (FGF4-KO)")
xKO = x[, safeSelect(selectedGroups)]
selectedFeatures = order(rowVars(exprs(xKO)), decreasing = TRUE)[seq_len(100)]


###################################################
### code chunk number 78: figHeatmapKO
###################################################
myHeatmap(xKO[selectedFeatures, ], collapseDuplicateFeatures = TRUE, haveColDend = TRUE)


###################################################
### code chunk number 79: checkClaim2
###################################################
stopifnot("Fgf4" %in% fData(xKO)$symbol[selectedFeatures])


###################################################
### code chunk number 80: dewtko
###################################################
x35 = x[, safeSelect(c("E3.5 (FGF4-KO)", "E3.5 (EPI)", "E3.5 (PE)"))]
f1 = f2 = x35$sampleGroup
f1[f1=="E3.5 (PE)" ] = NA
f2[f2=="E3.5 (EPI)"] = NA
x35$EPI = factor(f1, levels = c("E3.5 (FGF4-KO)", "E3.5 (EPI)"))
x35$PE  = factor(f2, levels = c("E3.5 (FGF4-KO)", "E3.5 (PE)"))
de = list(`EPI` = rowttests(x35, "EPI"),
          `PE`  = rowttests(x35, "PE"))
for(i in seq(along = de))
  de[[i]]$p.adj = p.adjust(de[[i]]$p.value, method = "BH")


###################################################
### code chunk number 81: figIndepFiltDeWtKo
###################################################
par(mfcol = c(3,2))
rkv = rank(-rowVars(exprs(x35)))
fdrthresh = 0.05
fcthresh =  1
for(i in seq(along = de)) {
  hist(de[[i]]$p.value, breaks = 100, col = "lightblue", main = names(de)[i], xlab = "p")
  plot(rkv, -log10(de[[i]]$p.value), pch = 16, cex = .25, main = "",
       xlab = "rank of overall variance", ylab = expression(-log[10]~p))
  plot(de[[i]]$dm, -log10(de[[i]]$p.value), pch = 16, cex = .25, main = "",
      xlab = "average log fold change", ylab = expression(-log[10]~p))
  abline(v = c(-1,1)*fcthresh[i], col = "red")
}


###################################################
### code chunk number 82: fig2DWTKO
###################################################
isSig = ((pmin(de$PE$p.adj, de$EPI$p.adj) < fdrthresh) &
         (pmax(abs(de$PE$dm), abs(de$EPI$dm)) > fcthresh))
table(isSig)
plot(de$PE$dm, de$EPI$dm, pch = 16, asp = 1,
     xlab = expression(log[2]~FGF4-KO / PE),
     ylab = expression(log[2]~FGF4-KO / EPI),
     cex  = ifelse(isSig, 0.6, 0.1), col = ifelse(isSig, "red", "grey"))


###################################################
### code chunk number 83: writetab35
###################################################
df = cbind(fData(x35)[isSig, ], `log2FC KO/PE` = de$PE$dm[isSig],
           `log2FC KO/EPI` = de$EPI$dm[isSig])
write.csv(df, file = "differentially_expressed_E3.5_vs_FGF4-KO.csv")
ctrls = grep("^AFFX", rownames(df), value = TRUE)


###################################################
### code chunk number 84: figFGF4ProbesInKO
###################################################
par(mfrow = c(2,2))
for(p in which(FGF4probes))
  plot(exprs(x35)[p, ], type = "p", pch = 16, col = x35$sampleColour, 
       main = rownames(fData(x325))[p], ylab = expression(log[2]~expression)) 


###################################################
### code chunk number 85: checkClaim3
###################################################
stopifnot(sum(FGF4probes)==4)


###################################################
### code chunk number 86: figCtrls
###################################################
stopifnot(length(ctrls)>=8)


###################################################
### code chunk number 87: figCtrls
###################################################
par(mfcol = c(4, 2))
for(p in ctrls[1:8])
  plot(exprs(x35)[p, ], type = "p", pch = 16, col = x35$sampleColour, 
       main = p, ylab = expression(log[2]~expression)) 


###################################################
### code chunk number 88: kegg (eval = FALSE)
###################################################
## allGenes = unique(fData(x35)$symbol)
## sigGenes = unique(df$symbol)
## ntot = length(allGenes)
## n1   = length(sigGenes)
## keggpw = as.list(mouse4302PATH2PROBE)
## fts = lapply(keggpw, function(ps) {
##    pwGenes = unique(as.character(mouse4302SYMBOL[ps]))
##    ovlGenes = intersect(pwGenes, sigGenes)
##    n2  = length(pwGenes)
##    nov = length(ovlGenes)
##    mat = matrix(c(ntot-n1-n2+nov, n2-nov, n1-nov, nov), ncol = 2, nrow = 2,
##                 dimnames = list(`pathway` = c("no","yes"), `diff. expressed` = c("no","yes")))
##    ft = fisher.test(mat)
##    list(mat = mat, p.value = ft$p.value, estimate = ft$estimate, 
##         genes = paste(sort(ovlGenes), collapse = ", ")) 
##    })


###################################################
### code chunk number 89: kegg
###################################################
keggpw = as.list(mouse4302PATH2PROBE)
dat  = de$PE$dm+de$EPI$dm
fts = lapply(keggpw, function(ps) {
  inpw = rownames(fData(x35)) %in% ps
  ft   = ks.test( dat[inpw], dat[!inpw] )
  list(p.value = ft$p.value, statistic = ft$statistic, n = sum(inpw)) 
})


###################################################
### code chunk number 90: pws1
###################################################
pws = c("04010", "04070", "04150", "04310", "04330", "04340", "04350", "04370", "04630")


###################################################
### code chunk number 91: pws2 (eval = FALSE)
###################################################
## pws = c(pws, names(fts)[ sapply(fts, function(x) 
##     with(x, (sum(mat[2, ])<=40) && (p.value<0.01) && (estimate>=2))) ] )


###################################################
### code chunk number 92: pws2
###################################################
pws = unique( c(pws, names(fts)[ sapply(fts, function(x) 
    with(x, (n<=100) && (p.value<0.01))) ] ) )


###################################################
### code chunk number 93: pwInfo
###################################################
query  = paste0("mmu", pws)
query  = split(query, seq(along = query) %/% 10)
pwInfo = unlist(lapply(query, keggGet), recursive = FALSE)


###################################################
### code chunk number 94: assertion
###################################################
stopifnot(all(pws %in% names(fts)))
stopifnot(length(pwInfo)==length(pws))


###################################################
### code chunk number 95: pwImgs (eval = FALSE)
###################################################
## pwImgs = lapply(query, keggGet, option = "image")
## stopifnot(length(pwImgs)==length(pws))


###################################################
### code chunk number 96: xtable (eval = FALSE)
###################################################
## df = data.frame(
##   `id`   = pws,
##   `name` = sub(" - Mus musculus (mouse)", "", sapply(pwInfo, `[[`, "NAME"), fixed = TRUE),
##   `n`    = sapply(pws, function(x) sum(fts[[x]]$mat[2, ])),
##   `differentially expressed genes` = sapply(pws, function(x) fts[[x]]$genes),
##   `p`          = as.character(signif(sapply(pws, function(x) fts[[x]]$p.value),  2)),
##   `odds-ratio` = as.character(signif(sapply(pws, function(x) fts[[x]]$estimate), 2)),
##   stringsAsFactors = FALSE, check.names = FALSE)


###################################################
### code chunk number 97: xtable
###################################################
df = data.frame(
  `id`   = pws,
  `name` = sub(" - Mus musculus (mouse)", "", sapply(pwInfo, `[[`, "NAME"), fixed = TRUE),
  `n`         =                     sapply(pws, function(x) fts[[x]]$n),
  `p`         = as.character(signif(sapply(pws, function(x) fts[[x]]$p.value),  2)),
  `statistic` = as.character(signif(sapply(pws, function(x) fts[[x]]$statistic), 2)),
  stringsAsFactors = FALSE, check.names = FALSE)


###################################################
### code chunk number 98: printxtable
###################################################
print(xtable(df,
  caption = paste("Gene set enrichment analysis of selected KEGG pathways, for the",
                   "differentially expressed genes between E3.5 FGF-4 KO and WT samples.",
                   "n: number of microarray features annotated to genes in the pathway."),
  label = "tab_KEGG", align = c("l", "l","p{4cm}","r","r","r")),
  type = "latex", file = "tab_KEGG.tex")


###################################################
### code chunk number 99: figShowPws
###################################################
par(mfrow = c(7,4)); stopifnot(length(pws)<=42)
for(i in seq(along = pws)) {
  inpw = factor(ifelse(rownames(fData(x35)) %in% keggpw[[pws[i]]],"in pathway","outside"))
  ord  = order(inpw)
  enr  = paste0("p=", df$p[i], if(as.numeric(df$p[i])<0.05)
                                  paste(" D=", df$statistic[i]) else "")
  multiecdf( (de$PE$dm+de$EPI$dm) ~ inpw, xlim = c(-1,1)*3,
             main = paste(df$name[i], enr, sep = "\n"),
             xlab = "mean difference between FGF4-KO and wildtype samples" )
}


###################################################
### code chunk number 100: defJSD
###################################################
H = function(p) -sum(p*log2(p))
JSD = function(m, normalize = TRUE) {
  stopifnot(is.matrix(m), all(dim(m)>1))
  if(normalize)
    m = m/rowSums(m)
  H(colMeans(m)) - mean(apply(m, 2, H))
}


###################################################
### code chunk number 101: defSDV
###################################################
SDV = function(m) sqrt(mean(rowVars(m)))


###################################################
### code chunk number 102: numberFeatures
###################################################
numberFeatures = c(50, 200, 1000, 4000)


###################################################
### code chunk number 103: computeDivergences
###################################################
computeDivergences = function(y, indices, fun) {
  y = y[indices, ]
  exprVals = y[, -1]
  strata = y[, 1]
  numStrata = max(strata)
  stopifnot(setequal(strata, seq_len(numStrata)))

  orderedFeatures = order(rowVars(t(exprVals)), decreasing = TRUE)
  sapply(numberFeatures, function(n) {
    selFeat = orderedFeatures[seq_len(n)]
    sapply(seq_len(numStrata), function(s)
             fun(exprVals[strata==s, selFeat, drop = FALSE]))
  })
}


###################################################
### code chunk number 104: strata
###################################################
x$sampleGroup = factor(x$sampleGroup)
x$strata = as.numeric(x$sampleGroup)


###################################################
### code chunk number 105: check
###################################################
stopifnot(!any(is.na(x$strata)))


###################################################
### code chunk number 106: doJSD
###################################################
bootJSD = boot(data = cbind(x$strata, t(exprs(x))),
               statistic = function(...) computeDivergences(..., fun = JSD),
               R = 100, strata = x$strata)
dim(bootJSD$t) = c(dim(bootJSD$t)[1], dim(bootJSD$t0))


###################################################
### code chunk number 107: doSD
###################################################
bootSDV = boot(data = cbind(x$strata, t(exprs(x))),
               statistic = function(...) computeDivergences(..., fun = SDV),
               R = 100, strata = x$strata)
dim(bootSDV$t) = c(dim(bootSDV$t)[1], dim(bootSDV$t0))


###################################################
### code chunk number 108: fig-JSD
###################################################
par(mfrow = c(length(numberFeatures),2))
colores = sampleColourMap[levels(x$sampleGroup)]
for(i in seq(along = numberFeatures)) {
  for(what in c("JSD", "SDV")){
    obj = get(paste("boot", what, sep = ""))
    boxplot(obj$t[, , i], main = sprintf("numberFeatures=%d", numberFeatures[i]),
            col = colores, border = colores, ylab = what, xaxtb = "n")
    px = seq_len(ncol(obj$t))
    text(x = px, y = par("usr")[3], labels = levels(x$sampleGroup),
         xpd = NA, srt = 45, adj = c(1,0.5), col = colores)
    points(px, obj$t0[, i], pch = 16)
  }
} 


###################################################
### code chunk number 109: xs
###################################################
xs = x[, pData(x)$sampleGroup %in%
   c("E3.25", "E3.5 (PE)", "E4.5 (PE)", "E3.5 (EPI)", "E4.5 (EPI)")]
xs$sampleGroup = factor(xs$sampleGroup)


###################################################
### code chunk number 110: myBoxplot
###################################################
myBoxplot = function(ps) {
  fac = factor(pData(xs)$sampleGroup)
  boxplot(exprs(xs)[ps, ] ~ fac, col = sampleColourMap[levels(fac)], lim = c(2,14),
          main = sprintf("%s (%s)", fData(x)[ps, "symbol"], ps), show.names = FALSE)
  text(seq(along = levels(fac)), par("usr")[3] - diff(par("usr")[3:4])*0.02,
       levels(fac), xpd = NA, srt = 90, adj = c(1,0.5), cex = 0.8)
}


###################################################
### code chunk number 111: showGenes
###################################################
exemplaryGenes = read.table(header = TRUE, stringsAsFactors = FALSE, text = "
symbol   thclass  probeset
Gata6      C-1 1425464_at
Tom1l1     C-1 1439337_at
Sox2       C-1 1416967_at
Pdgfra    PE-2 1438946_at
Sox17     PE-2 1429177_x_at
P4ha2     PE-2 1417149_at
Gata4     PE-3 1418863_at
Aldh18a1  PE-3 1415836_at
Col4a1    PE-3 1452035_at
Col4a2    PE-3 1424051_at
Cubn      PE-3 1426990_at
Lamb1     PE-3 1424113_at
Dab2      PE-4 1423805_at
Lrp2      PE-4 1427133_s_at
Amn       PE-4 1417920_at
Fgf4     EPI-2 1420085_at
Nanog    EPI-2 1429388_at
Tdgf1    EPI-2 1450989_at
Cldn4    EPI-4 1418283_at
Enox1    EPI-4 1436799_at")

pdf(file = "exemplaryGenes.pdf", width = 8, height = 11)
par(mfrow = c(5,3))
for(i in seq_len(nrow(exemplaryGenes))) {
  wh = featureNames(x)[fData(x)[, "symbol"]==exemplaryGenes[i, "symbol"]]
  stopifnot(length(wh)>=1)
  sapply(wh, myBoxplot)
}
dev.off()


###################################################
### code chunk number 112: fig-ExemplaryGenes
###################################################
layout(rbind(c(1, 4,   7, 13),
             c(2, 5,   8, 14),
            c(21, 6,   9, 15),
             c(3, 16, 10, 19),
            c(22, 17, 11, 20),
            c(23, 18, 12, 24)))

xs = x
xs$sampleGroup = factor(xs$sampleGroup)
xs$sampleGroup = relevel(xs$sampleGroup, "E4.5 (PE)")
xs$sampleGroup = relevel(xs$sampleGroup, "E4.5 (EPI)")
xs$sampleGroup = relevel(xs$sampleGroup, "E3.5 (PE)")
xs$sampleGroup = relevel(xs$sampleGroup, "E3.5 (EPI)")
xs$sampleGroup = relevel(xs$sampleGroup, "E3.25")

for(i in seq_len(nrow(exemplaryGenes))) 
  myBoxplot(exemplaryGenes[i, "probeset"])


###################################################
### code chunk number 113: avglfc
###################################################
greater = function(cond1, cond2, thresh = 2) {
  stopifnot(all(cond1 %in% names(groups)), all(cond2 %in% names(groups)))
  rm1 = lapply(groups[cond1], function(j) rowMeans(exprs(x)[, j]))
  rm2 = lapply(groups[cond2], function(j) rowMeans(exprs(x)[, j]))
  res = rep(TRUE, nrow(x))
  for(v1 in rm1) for(v2 in rm2) 
      res = res & ((v1-v2) > thresh)
  return(res)
}


###################################################
### code chunk number 114: minLevel
###################################################
xquantiles = apply(exprs(x), 1, quantile, probs = c(0.025, 0.975))
minLevel = xquantiles[1, ]
maxLevel = xquantiles[2, ]
trsf2Zero2One = function(x) {
  x = (x-minLevel)/(maxLevel-minLevel)
  x[x<0] = 0
  x[x>1] = 1
  return(x)
}
xsc = x
exprs(xsc) = trsf2Zero2One(exprs(x))


###################################################
### code chunk number 115: isOff
###################################################
isOff = function(cond, thresh = 2) {
  stopifnot(all(cond %in% names(groups)))
  samps = unlist(groups[cond])
  rowSums(exprs(x)[, samps] > minLevel+thresh) <= ceiling(length(samps) * 0.1)
}


###################################################
### code chunk number 116: classification
###################################################
  `C-1` =  greater("E3.25", c("E3.5 (EPI)", "E4.5 (EPI)", "E3.5 (PE)", "E4.5 (PE)"))
`EPI-2` =  greater("E3.5 (EPI)", c("E3.25")) &
           greater("E3.5 (EPI)", "E4.5 (EPI)", thresh = 0.5) &
           isOff("E3.5 (PE)")  & isOff("E4.5 (PE)" )
 `PE-2` =  greater("E3.5 (PE)",  c("E3.25")) &
           greater("E3.5 (PE)",  "E4.5 (PE)", thresh = 0.5)  & 
           isOff("E3.5 (EPI)") & isOff("E4.5 (EPI)")
`EPI-4` = (greater("E4.5 (EPI)", c("E3.25", "E3.5 (EPI)")) &
           isOff(c("E3.25", "E3.5 (EPI)")) & 
           isOff("E4.5 (PE)" ))
 `PE-4` = (greater("E4.5 (PE)", c("E3.25", "E3.5 (PE)")) &
           isOff(c("E3.25", "E3.5 (PE)")) & 
           isOff("E4.5 (EPI)"))
`EPI-3` = (greater("E3.5 (EPI)", "E3.25", thresh = 0.5) &
           greater("E4.5 (EPI)", "E3.5 (EPI)", thresh = 0.5) &
           greater("E4.5 (EPI)", "E3.25", thresh = 3)   & 
           isOff("E4.5 (PE)") & !`EPI-4`)
 `PE-3` = (greater("E3.5 (PE)",  "E3.25", thresh = 0.5) &
           greater("E4.5 (PE)",  "E3.5 (PE)",  thresh = 0.5) &
           greater("E4.5 (PE)",  "E3.25", thresh = 3) & 
           isOff("E4.5 (EPI)") & !`PE-4`)
thclasses = cbind(`C-1`, `EPI-2`, `EPI-3`, `EPI-4`, `PE-2`,  `PE-3`,  `PE-4`)


###################################################
### code chunk number 117: onegroup
###################################################
multiclass = thclasses[rowSums(thclasses)>1,, drop = FALSE]
stopifnot(nrow(multiclass)==0)


###################################################
### code chunk number 118: cs
###################################################
colSums(thclasses)


###################################################
### code chunk number 119: figHeatmapthclasses
###################################################
agr = groups[c("E3.25", "E3.5 (EPI)", "E4.5 (EPI)", "E3.5 (PE)", "E4.5 (PE)")]
fgr = apply(thclasses, 2, which)
xsub = xsc[unlist(fgr), unlist(agr)]
mat = exprs(xsub)
rownames(mat) = fData(xsub)[, "symbol"]
myHeatmap2(mat, keeprownames = TRUE,
           rowGroups = factor(rep(seq(along = fgr), listLen(fgr))),
           colGroups = factor(rep(seq(along = agr), listLen(agr))))


###################################################
### code chunk number 120: writecsv
###################################################
out = do.call(rbind, lapply(seq(along = fgr), function(i) 
    cbind(`class` = names(fgr)[i], fData(xsc)[fgr[[i]], ])))
write.csv(out, file = "featureclassification.csv")


###################################################
### code chunk number 121: fig2bGenes
###################################################
fig2bGenes = c("Aldh18a1", "Col4a1", "Cubn", "Foxq1", "Gata4", "Lamb1", "Serpinh1")


###################################################
### code chunk number 122: rbGenes
###################################################
rbGenes = sort(unique(fData(x)$symbol[thclasses[, "PE-3"]]))
rbGenes


###################################################
### code chunk number 123: setdiff
###################################################
setdiff(fig2bGenes, rbGenes)


###################################################
### code chunk number 124: load
###################################################
data("xq")


###################################################
### code chunk number 125: hmprep
###################################################
xq$sampleColours = sampleColourMap[sub("[[:digit:]]+ ", "", sampleNames(xq))]
stopifnot(!any(is.na(xq$sampleColours)))


###################################################
### code chunk number 126: myHeatmap3
###################################################
myHeatmap3 = function(x, log = TRUE, 
  col = colorRampPalette(brewer.pal(9, "Blues"))(ncol), ncol = 100, ...) {
  mat = exprs(x)
  if(log){
    mat[mat<1] = 1
    mat = log10(mat)
  }
  rownames(mat) = fData(x)[, "symbol"]
  heatmap.2(mat, trace = "none", dendrogram = "none", scale = "none", col = col,
            keysize = 0.9, ColSideColors = x$sampleColour, margins = c(5,7), ...)
}


###################################################
### code chunk number 127: fighm1
###################################################
myHeatmap3(xq)


###################################################
### code chunk number 128: fighm2
###################################################
selectedGenes = c("Col4a1", "Lamab1", "Cubn", "Gata4", "Serpinh1", "Foxq1", "Aldh18a1")
myHeatmap3(xq[selectedGenes, ])


###################################################
### code chunk number 129: checkAssertion
###################################################
stopifnot(length(selectedGenes)==7)


###################################################
### code chunk number 130: selectedPE
###################################################
selectedSamples = (xq$Cell.type %in% c("ICM", "PE"))
table(xq$sampleGroup[selectedSamples])
sxq = xq[selectedGenes, selectedSamples]
groups = c("E3.25", "E3.5 (PE)", "E4.5 (PE)")
sxq$fsampleGroup = factor(sxq$sampleGroup, levels = groups)
stopifnot(!any(is.na(sxq$fsampleGroup)))


###################################################
### code chunk number 131: fighm3
###################################################
myHeatmap3(sxq)


###################################################
### code chunk number 132: discr
###################################################
groupmedians = t(apply(exprs(sxq), 1, function(v) tapply(v, sxq$sampleGroup, median)))
stopifnot(identical(colnames(groupmedians), groups), length(groups)==3)
discrthreshs = (groupmedians[, 2:3] + groupmedians[, 1:2]) / 2
stst = sapply(split(seq(along = sxq$sampleGroup), sxq$sampleGroup), 
    function(v) {i1 = head(v,1); i2 = tail(v,1); stopifnot(identical(v, i1:i2)); c(i1,i2)})


###################################################
### code chunk number 133: show
###################################################
groupmedians
discrthreshs


###################################################
### code chunk number 134: figthr
###################################################
op = par(mfrow = c(nrow(sxq), 1), mai = c(0.1, 0.7, 0.01, 0.01))
for(j in seq_len(nrow(sxq))) {
  plot(exprs(sxq)[j, ], type = "n", xaxt = "n", ylab = fData(sxq)$symbol[j])
  segments(x0 = stst[1, ], x1 = stst[2, ], y0 = groupmedians[j, ], col = "#808080", lty=3)
  segments(x0 = stst[1,1:2], x1 = stst[2,2:3], y0 = discrthreshs[j, ], col = "#404040")
  points(exprs(sxq)[j, ], pch = 16, col = sxq$sampleColour)
}
par(op)


###################################################
### code chunk number 135: discretisation1
###################################################
discretize = function(x) {
  exprs(x) = t(sapply(seq_len(nrow(exprs(x))), function(r) {
      as.integer(cut(exprs(x)[r, ], breaks = c(-Inf, discrthreshs[r, ],+Inf)))
  }))
  return(x)
}


###################################################
### code chunk number 136: discretisation2
###################################################
sxqd = discretize(sxq)


###################################################
### code chunk number 137: fighmd
###################################################
myHeatmap3(sxqd, log = FALSE, Colv = FALSE, Rowv = FALSE, key = FALSE, ncol = 3)


###################################################
### code chunk number 138: bruteForceOptimisation1
###################################################
stopifnot(all(exprs(sxqd)%in%(1:3)))

costfun1 = function(x, sg) {
  k = (sg=="E3.5 (PE)")
  mean( (x[-1,k]==1) & (x[-nrow(x),k]>1 ) )
}
costfun2 = function(x, sg) {
  k = (sg=="E4.5 (PE)")
  mean( (x[-1,k]<3 ) & (x[-nrow(x),k]==3) )
}

perms = permutations(nrow(sxq), nrow(sxq))
bruteForceOptimisation = function(fun, samps) {
  if (missing(samps)) samps = rep(TRUE, ncol(sxqd))
  apply(perms, 1, function(i) fun(exprs(sxqd)[i, samps], sxq$sampleGroup[samps]))
}

costs = list(
 `E3.25 -> E3.5` = bruteForceOptimisation(costfun1),
  `E3.5 -> E4.5` = bruteForceOptimisation(costfun2))
whopt  = lapply(costs, function(v) which(v==min(v)))


###################################################
### code chunk number 139: fighmd
###################################################
outf = file("fighmd.tex", open="w")
for(i in seq(along = costs)) {
  sxqdsub = switch(i,
     { 
       rv = sxqd[, sxqd$sampleGroup %in% groups[2] ]
       exprs(rv)[ exprs(rv)>2 ] = 2
       rv
     }, {
       rv = sxqd[, sxqd$sampleGroup %in% groups[3] ]
       exprs(rv)[ exprs(rv)<2 ] = 2
       exprs(rv) = exprs(rv)-1
       rv
     })
  columnOrder = order(sxqdsub$sampleGroup, colSums(exprs(sxqdsub)))
  for(w in whopt[[i]]) {
    fn = sprintf("fighmd-%d-%d.pdf", i, w)
    pdf(fn, width = 7, height = 3)
    myHeatmap3(sxqdsub[perms[w, ], columnOrder], log = FALSE, Colv = FALSE, Rowv = FALSE,
               key = FALSE, ncol = 2, breaks = seq(0.5, 2.5, by = 1))
    dev.off()
    cat("\\includegraphics[width=0.49\\textwidth]{", fn, "}\n", file = outf, sep = "")
  }
}
close(outf)


###################################################
### code chunk number 140: printcosts
###################################################
for(i in seq(along = costs)) {
  k = unique(costs[[i]][whopt[[i]]])
  stopifnot(length(k)==1)
  cat(sprintf("%14s: cost %g\n", names(costs)[[i]], k))
}


###################################################
### code chunk number 141: boot-takes-a-long-time
###################################################
bopt = lapply(list(costfun1, costfun2), function(fun) {
  boot(data = seq_len(ncol(sxqd)), 
       statistic = function(dummy, idx) min(bruteForceOptimisation(fun, idx)),
       R = 250) 
})


###################################################
### code chunk number 142: figboot
###################################################
names(bopt) = names(costs)
lapply(bopt, `[[`, "t0")
multiecdf(lapply(bopt, `[[`, "t"))
t.test(x = bopt[[1]]$t, y = bopt[[2]]$t)


###################################################
### code chunk number 143: mdsplot
###################################################
data("xql")

labeledSampleColourMap = c(brewer.pal(5, "RdGy")[c(1,2)], brewer.pal(5, "BrBG")[c(4,5)])
names(labeledSampleColourMap) = c("E4.5_high", "E3.5_high", "E3.5_low",  "E4.5_low")
labeledGroups = with(pData(xql), list(
  `E3.5_high` = which(Label=="High" & Embryonic.day=="E3.5"),
  `E3.5_low`  = which(Label=="Low"  & Embryonic.day=="E3.5"),
  `E4.5_high`	= which(Label=="High" & Embryonic.day=="E4.5"),
  `E4.5_low`  = which(Label=="Low"  & Embryonic.day=="E4.5")))

labeledSampleColours = rep(NA_character_, ncol(xql))
for(i in seq(along = labeledGroups))
 labeledSampleColours[labeledGroups[[i]]] = labeledSampleColourMap[names(labeledGroups)[i]]
xql$sampleColours = labeledSampleColours

md = isoMDS( dist(t(log(exprs(xql),2))))$points
plot(md, col = xql$sampleColours, pch = 16, asp = 1, xlab = "",ylab = "")


###################################################
### code chunk number 144: figure5b
###################################################
xs = x[, (x$genotype %in% "WT") & (x$Embryonic.day %in% c("E3.25","E3.5","E4.5"))]

Fgf3 = c("1441914_x_at")
Fgf4 = c("1420086_x_at")
Fgfr2 = c("1433489_s_at")
Fgfr3 = c("1421841_at")
Fgfr4 = c("1418596_at")

correlationPlot = function(ID){
  xsl = xs[ID, ]
  mat = exprs(xsl)
  rownames(mat) = fData(xsl)[, "symbol"]
  plot(t(mat), pch = 16, asp = 1, cex = 1.25, cex.lab = 1.2, col = xs$sampleColour,
       xlim = c(2,12), ylim = c(3,11))
}

par(mfrow = c(1,4))
correlationPlot(c(Fgf4,Fgfr2))
correlationPlot(c(Fgf4,Fgfr3))
correlationPlot(c(Fgf3,Fgfr3))
correlationPlot(c(Fgf3,Fgfr4))


###################################################
### code chunk number 145: sessionInfo
###################################################
toLatex(sessionInfo())