library(IRkernel)
library(IRdisplay)
library(RColorBrewer)
library(repr)
library(base64enc)

suppressPackageStartupMessages({
    library(xlsx)
    library(destiny)
})

options(device = function(...) png('/dev/null', 7, 6, 'in', res = 120))
options(repr.plot.width = 7, repr.plot.height = 6)
options(jupyter.plot_mimetypes = c('application/pdf', 'image/png'))

setHook('on.rgl.close', function(...) {
    name <- tempfile()
    par3d(windowRect = c(0, 0, 1200,1200))
    Sys.sleep(1)
    
    rgl.snapshot(  filename = paste0(name, '.png'))
    #rgl.postscript(filename = paste0(name, '.pdf'), fmt='pdf')
    
    res <- getOption('repr.plot.res')
    
    publish_mimebundle(list(
        'image/png'       = base64encode(paste0(name, '.png'))
     #, 'application/pdf' = base64encode(paste0(name, '.pdf'))
    ), list(
        width  = res * getOption('repr.plot.width'),
        height = res * getOption('repr.plot.height')
    ))
}, 'replace')

library(xlsx)
raw.ct <- read.xlsx('mmc4.xls', sheetName = 'Sheet1')

raw.ct[1:9, 1:9]  #preview of a few rows and columns

library(destiny)

ct <- as.ExpressionSet(raw.ct)
ct

num.cells <- gsub('^(\\d+)C.*$', '\\1', ct$Cell)
ct$num.cells <- as.integer(num.cells)

# cells from 2+ cell embryos
have.duplications <- ct$num.cells > 1
# cells with values ≤ 28
normal.vals <- apply(exprs(ct), 2, function(sample) all(sample <= 28))

cleaned.ct <- ct[, have.duplications & normal.vals]

housekeepers <- c('Actb', 'Gapdh')  # houskeeper gene names

normalizations <- colMeans(exprs(cleaned.ct)[housekeepers, ])

guo.norm <- cleaned.ct
exprs(guo.norm) <- exprs(guo.norm) - normalizations

library(destiny)
#data(guo.norm)
dif <- DiffusionMap(guo.norm)

plot(dif)

library(RColorBrewer)
palette(brewer.pal(6, 'Spectral')) # configure color palette

plot(dif, pch = 20,        # pch for prettier points
     col.by = 'num.cells', # or “col” to use a vector or a single color
     legend.main = 'Cell stage')

plot(dif, 1:2, pch = 20, col.by = 'num.cells',
     legend.main = 'Cell stage')

library(rgl)
plot3d(eigenvectors(dif)[, 1:3],
       col = log2(guo.norm$num.cells),
       type = 's', radius = .01)
view3d(theta = 10, phi = 30, zoom = .8)
# now use your mouse to rotate the plot in the window
rgl.close()

library(ggplot2)
qplot(DC1, DC2, data = dif, colour = factor(num.cells)) +
    scale_color_brewer(palette = 'Spectral')
# or alternatively:
#ggplot(dif, aes(DC1, DC2, colour = factor(num.cells))) + ...

plot(eigenvalues(dif), ylim = 0:1, pch = 20,
     xlab = 'Diffusion component (DC)', ylab = 'Eigenvalue')

oh <- options('repr.plot.height')
options(repr.plot.height = 3)

par(mfrow = c(1, 2), mar = c(2,2,2,2))

plot(dif, 3:4,   pch = 20, col.by = 'num.cells', draw.legend = FALSE)
plot(dif, 19:20, pch = 20, col.by = 'num.cells', draw.legend = FALSE)

options(oh)

sigmas <- find.sigmas(guo.norm, verbose = FALSE)
optimal.sigma(sigmas)

par(pch = 20, mfrow = c(2, 2), mar = c(3,2,2,2))

for (sigma in c(2, 5, optimal.sigma(sigmas), 100))
    plot(DiffusionMap(guo.norm, sigma), 1:2,
         main = substitute(sigma == s, list(s = round(sigma,2))),
         col.by = 'num.cells', draw.legend = FALSE)

hist(exprs(cleaned.ct)['Aqp3', ], breaks = 20,
     xlab = 'Ct of Aqp3', main = 'Histogram of Aqp3 Ct',
     col = 'slategray3', border = 'white')

dilutions <- read.xlsx('mmc6.xls', 1L)
dilutions$Cell <- NULL #remove annotation column

lods <- apply(dilutions, 2, function(col) col[[max(which(col != 28))]])
lod <- ceiling(median(lods))
lod

lod.norm <- ceiling(median(lods) - mean(normalizations))
max.cycles.norm <- ceiling(40 - mean(normalizations))

list(lod.norm = lod.norm, max.cycles.norm = max.cycles.norm)

guo <- guo.norm
exprs(guo)[exprs(cleaned.ct) >= 28] <- lod.norm

thresh.dif <- DiffusionMap(guo,
                           censor.val = lod.norm,
                           censor.range = c(lod.norm, max.cycles.norm),
                           verbose = FALSE)

plot(thresh.dif, 1:2, col.by = 'num.cells', pch = 20,
     legend.main = 'Cell stage')

# remove rows with divisionless cells
ct.w.missing <- ct[, ct$num.cells > 1L]
# and replace values larger than the baseline
exprs(ct.w.missing)[exprs(ct.w.missing) > 28] <- NA

housekeep <- colMeans(exprs(ct.w.missing)[housekeepers, ],
                      na.rm = TRUE)

w.missing <- ct.w.missing
exprs(w.missing) <- exprs(w.missing) - housekeep

exprs(w.missing)[is.na(exprs(ct.w.missing))] <- lod.norm

dif.w.missing <- DiffusionMap(w.missing,
                              censor.val = lod.norm,
                              censor.range = c(lod.norm,
                                               max.cycles.norm),
                              missing.range = c(1, 40),
                              verbose = FALSE)

plot(dif.w.missing, 1:2, col.by = 'num.cells', pch = 20,
     legend.main = 'Cell stage')

ct64 <- guo[, guo$num.cells == 64]

dif64 <- DiffusionMap(ct64)

ct32 <- guo[, guo$num.cells == 32]
pred32 <- dm.predict(dif64, ct32)

par(mar = c(2,2,1,5), pch = 20)
plot(dif64,  1:2,      col     = palette()[[6]],
     new.dcs = pred32, col.new = palette()[[5]])
colorlegend(c(32L, 64L), palette()[5:6])