library(IRkernel) library(IRdisplay) library(repr) library(base64enc) library(ggplot2) # would shadow Biobase::exprs suppressPackageStartupMessages({ library(readxl) library(destiny) library(Biobase) }) 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') # doesn’t work with spheres 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(readxl) raw_ct <- read_xls('mmc4.xls', 'Sheet1') raw_ct[1:9, 1:9] #preview of a few rows and columns library(destiny) library(Biobase) # as.data.frame because ExpressionSets # do not support readxl’s “tibbles” ct <- as.ExpressionSet(as.data.frame(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(smp) all(smp <= 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) dm <- DiffusionMap(guo_norm) plot(dm) palette(cube_helix(6)) # configure color palette plot(dm, pch = 20, # pch for prettier points col_by = 'num_cells', # or “col” with a vector or one color legend_main = 'Cell stage') plot(dm, 1:2, col_by = 'num_cells', legend_main = 'Cell stage') library(rgl) plot3d( eigenvectors(dm)[, 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 = dm, colour = factor(num_cells)) + scale_color_cube_helix() # or alternatively: #ggplot(dif, aes(DC1, DC2, colour = factor(num.cells))) + ... qplot(y = eigenvalues(dm)) + theme_minimal() + labs(x = 'Diffusion component (DC)', y = 'Eigenvalue') oh <- options('repr.plot.height') options(repr.plot.height = 3) plot(dm, 3:4, col_by = 'num_cells', draw_legend = FALSE) plot(dm, 19:20, col_by = 'num_cells', draw_legend = FALSE) options(oh) hist(exprs(cleaned_ct)['Aqp3', ], breaks = 20, xlab = 'Ct of Aqp3', main = 'Histogram of Aqp3 Ct', col = palette()[[4]], border = 'white') dilutions <- read_xls('mmc6.xls', 1L) dilutions$Cell <- NULL # remove annotation column get_lod <- function(gene) gene[[max(which(gene != 28))]] lods <- apply(dilutions, 2, get_lod) 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_dm <- DiffusionMap(guo, censor_val = lod_norm, censor_range = c(lod_norm, max_cycles_norm), verbose = FALSE) plot(thresh_dm, 1:2, col_by = 'num_cells', 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', legend_main = 'Cell stage') ct64 <- guo[, guo$num_cells == 64] dm64 <- DiffusionMap(ct64) ct32 <- guo[, guo$num_cells == 32] pred32 <- dm_predict(dm64, ct32) plot(dm64, 1:2, col = palette()[[6]], new_dcs = pred32, col_new = palette()[[4]])