diff --git a/Functions.R b/Functions.R new file mode 100644 index 0000000..6e35a1b --- /dev/null +++ b/Functions.R @@ -0,0 +1,1359 @@ + +##### Modified Seurat functions +##### SCeNGEA APP April 2020 +#options(repos = BiocManager::repositories()) +#options("repos") + +load("Dataset_20October_2020.rda") + +utr <- c("WBGene00023498","WBGene00023497","WBGene00004397","WBGene00006843", + "WBGene00004010","WBGene00006789","WBGene00001135","WBGene00001079", + "WBGene00001135","WBGene00006783","WBGene00000501","WBGene00006788", + "WBGene00001555") + +### GetXYAesthetics ---- +GetXYAesthetics <- function(plot, geom = 'GeomPoint', plot.first = TRUE) { + geoms <- sapply( + X = plot$layers, + FUN = function(layer) { + return(class(x = layer$geom)[1]) + } + ) + geoms <- which(x = geoms == geom) + if (length(x = geoms) == 0) { + stop("Cannot find a geom of class ", geom) + } + geoms <- min(geoms) + if (plot.first) { + x <- as.character(x = plot$mapping$x %||% plot$layers[[geoms]]$mapping$x)[2] + y <- as.character(x = plot$mapping$y %||% plot$layers[[geoms]]$mapping$y)[2] + } else { + x <- as.character(x = plot$layers[[geoms]]$mapping$x %||% plot$mapping$x)[2] + y <- as.character(x = plot$layers[[geoms]]$mapping$y %||% plot$mapping$y)[2] + } + return(list('x' = x, 'y' = y)) +} + +### LabelClusters2 ---- +LabelClusters2 <- function( + plot, + id, + clusters = NULL, + labels = NULL, + split.by = NULL, + repel = TRUE, + dims = dims, + l1, l2, + ... +) { + xynames <- unlist(x = GetXYAesthetics(plot = plot), use.names = TRUE) + if (!id %in% colnames(x = plot$data)) { + stop("Cannot find variable ", id, " in plotting data") + } + if (!is.null(x = split.by) && !split.by %in% colnames(x = plot$data)) { + warning("Cannot find splitting variable ", id, " in plotting data") + split.by <- NULL + } + data <- plot$data[, c(xynames, id, split.by)] + possible.clusters <- as.character(x = na.omit(object = unique(x = data[, id]))) + groups <- clusters %||% as.character(x = na.omit(object = unique(x = data[, id]))) + if (any(!groups %in% possible.clusters)) { + stop("The following clusters were not found: ", paste(groups[!groups %in% possible.clusters], collapse = ",")) + } + labels.loc <- lapply( + X = groups, + FUN = function(group) { + data.use <- data[data[, id] == group, , drop = FALSE] + data.medians <- if (!is.null(x = split.by)) { + do.call( + what = 'rbind', + args = lapply( + X = unique(x = data.use[, split.by]), + FUN = function(split) { + medians <- apply( + X = data.use[data.use[, split.by] == split, xynames, drop = FALSE], + MARGIN = 2, + FUN = median, + na.rm = TRUE + ) + medians <- as.data.frame(x = t(x = medians)) + medians[, split.by] <- split + return(medians) + } + ) + ) + } else { + as.data.frame(x = t(x = apply( + X = data.use[, xynames, drop = FALSE], + MARGIN = 2, + FUN = median, + na.rm = TRUE + ))) + } + data.medians[, id] <- group + return(data.medians) + } + ) + labels.loc <- do.call(what = 'rbind', args = labels.loc) + labels <- labels %||% groups + if (length(x = unique(x = labels.loc[, id])) != length(x = labels)) { + stop("Length of labels (", length(x = labels), ") must be equal to the number of clusters being labeled (", length(x = labels.loc), ").") + } + names(x = labels) <- groups + for (group in groups) { + labels.loc[labels.loc[, id] == group, id] <- labels[group] + } + #print(head(labels.loc)) + if (length(l1) > 1){ + labels.loc <- dplyr::filter(labels.loc, UMAP_1 > l1[1], UMAP_1l2[1], UMAP_2 0] + cutoff <- quantile(x = data, probs = this.quantile) + } + return(as.numeric(x = cutoff)) +} + + +### BlendMatrix ---- +BlendMatrix <- function(n = 10, col.threshold = 0.5) { + if (0 > col.threshold || col.threshold > 1) { + stop("col.threshold must be between 0 and 1") + } + return(outer( + X = 1:n, + Y = 1:n, + FUN = function(i, j) { + red <- 1 / (1 + exp(x = -(i - col.threshold * n) / 0.9)) + green <- 1 / (1 + exp(x = -(j - col.threshold * n) / 0.9)) + blue <- 0.2 + alpha <- lapply(X = list(red, green, blue), FUN = '^', 40) + alpha <- Reduce(f = '+', x = alpha) + alpha <- 0.99 - 0.1 * exp(x = -alpha / 1) + return(rgb(red = red, green = green, blue = blue, alpha = alpha)) + } + )) +} + +BlendExpression <- function(data) { + if (ncol(x = data) != 2) { + stop("'BlendExpression' only blends two features") + } + features <- colnames(x = data) + data <- as.data.frame(x = apply( + X = data, + MARGIN = 2, + FUN = function(x) { + return(round(x = 9 * (x - min(x)) / (max(x) - min(x)))) + } + )) + data[, 3] <- data[, 1] + data[, 2] * 10 + ## CHECK!! + aaa <- features[1] + bbb <- features[2] + features <- c(as.character(dplyr::filter(gene_list, gene_id == features[1])$gene_name), as.character(dplyr::filter(gene_list, gene_id == features[2])$gene_name) ) + if (aaa %in% utr) { features[1] <- paste0(features[1],' (Transgene!!)')} + if (bbb %in% utr) { features[2] <- paste0(features[2],' (Transgene!!)')} + colnames(x = data) <- c(features, paste(features, collapse = '_')) + for (i in 1:ncol(x = data)) { + data[, i] <- factor(x = data[, i]) + } + return(data) +} + +### SingleDimPlot2 ---- +SingleDimPlot2 <- function( + data, + dims, + col.by = NULL, + cols = NULL, + pt.size = NULL, + shape.by = NULL, + order = NULL, + label = FALSE, + repel = FALSE, + label.size = 4, + cells.highlight = NULL, + cols.highlight = 'red', + sizes.highlight = 1, + na.value = 'grey50' +) { + pt.size <- pt.size %||% AutoPointSize(data = data) + if (length(x = dims) != 2) { + stop("'dims' must be a two-length vector") + } + if (!is.data.frame(x = data)) { + data <- as.data.frame(x = data) + } + if (is.character(x = dims) && !all(dims %in% colnames(x = data))) { + stop("Cannot find dimensions to plot in data") + } else if (is.numeric(x = dims)) { + dims <- colnames(x = data)[dims] + } + if (!is.null(x = cells.highlight)) { + highlight.info <- SetHighlight( + cells.highlight = cells.highlight, + cells.all = rownames(x = data), + sizes.highlight = sizes.highlight %||% pt.size, + cols.highlight = cols.highlight, + col.base = cols[1] %||% 'black', + pt.size = pt.size + ) + order <- highlight.info$plot.order + data$highlight <- highlight.info$highlight + col.by <- 'highlight' + pt.size <- highlight.info$size + cols <- highlight.info$color + } + if (!is.null(x = order) && !is.null(x = col.by)) { + if (typeof(x = order) == "logical"){ + if (order) { + data <- data[order(data[, col.by]), ] + } + } else { + order <- rev(x = c( + order, + setdiff(x = unique(x = data[, col.by]), y = order) + )) + data[, col.by] <- factor(x = data[, col.by], levels = order) + new.order <- order(x = data[, col.by]) + data <- data[new.order, ] + if (length(x = pt.size) == length(x = new.order)) { + pt.size <- pt.size[new.order] + } + } + } + if (!is.null(x = col.by) && !col.by %in% colnames(x = data)) { + warning("Cannot find ", col.by, " in plotting data, not coloring plot") + col.by <- NULL + } else { + # col.index <- grep(pattern = col.by, x = colnames(x = data), fixed = TRUE) + col.index <- match(x = col.by, table = colnames(x = data)) + if (grepl(pattern = '^\\d', x = col.by)) { + # Do something for numbers + col.by <- paste0('x', col.by) + } else if (grepl(pattern = '-', x = col.by)) { + # Do something for dashes + col.by <- gsub(pattern = '-', replacement = '.', x = col.by) + } + colnames(x = data)[col.index] <- col.by + } + if (!is.null(x = shape.by) && !shape.by %in% colnames(x = data)) { + warning("Cannot find ", shape.by, " in plotting data, not shaping plot") + } + plot <- ggplot(data = data) + + geom_point( + mapping = aes_string( + x = dims[1], + y = dims[2], + color = paste0("`", col.by, "`"), + shape = shape.by + ), + size = pt.size + ) + + guides(color = guide_legend(override.aes = list(size = 3))) + + labs(color = NULL) + if (label && !is.null(x = col.by)) { + plot <- LabelClusters( + plot = plot, + id = col.by, + repel = repel, + size = label.size + ) + } + if (!is.null(x = cols)) { + plot <- plot + if (length(x = cols) == 1) { + scale_color_brewer(palette = cols, na.value = na.value) + } else { + scale_color_manual(values = cols, na.value = na.value) + } + } + plot <- plot + theme_minimal() + return(plot) +} + +Melt <- function(x) { + if (!is.data.frame(x = x)) { + x <- as.data.frame(x = x) + } + return(data.frame( + rows = rep.int(x = rownames(x = x), times = ncol(x = x)), + cols = unlist(x = lapply(X = colnames(x = x), FUN = rep.int, times = nrow(x = x))), + vals = unlist(x = x, use.names = FALSE) + )) +} + +`%iff%` <- function(lhs, rhs) { + if (!is.null(x = lhs)) { + return(rhs) + } else { + return(lhs) + } +} + +### FeaturePlot2 ---- +FeaturePlot2 <- function( + object, + features, + dims = c(1, 2), + cells = NULL, + cols = c("lightgrey", "blue"), + pt.size = NULL, + order = FALSE, + min.cutoff = NA, + max.cutoff = NA, + reduction = NULL, + split.by = NULL, + shape.by = NULL, + blend = FALSE, + blend.threshold = 0.5, + label = FALSE, + label.size = 4, + ncol = NULL, + combine = TRUE, + coord.fixed = FALSE, + by.col = TRUE, + ranges=NULL +) { + no.right <- theme( + axis.line.y.right = element_blank(), + axis.ticks.y.right = element_blank(), + axis.text.y.right = element_blank(), + axis.title.y.right = element_text( + face = "bold", + size = 14, + margin = margin(r = 7) + ) + ) + if (is.null(reduction)) { + default_order <- c('umap', 'tsne', 'pca') + reducs <- which(default_order %in% names(object@reductions)) + reduction <- default_order[reducs[1]] + } + if (length(x = dims) != 2 || !is.numeric(x = dims)) { + stop("'dims' must be a two-length integer vector") + } + if (blend && length(x = features) != 2) { + stop("Blending feature plots only works with two features") + } + dims <- paste0(Key(object = object[[reduction]]), dims) + cells <- cells %||% colnames(x = object) + data <- FetchData(object = object, vars = c(dims, features), cells = cells) + features <- colnames(x = data)[3:ncol(x = data)] + min.cutoff <- mapply( + FUN = function(cutoff, feature) { + return(ifelse( + test = is.na(x = cutoff), + yes = min(data[, feature]), + no = cutoff + )) + }, + cutoff = min.cutoff, + feature = features + ) + max.cutoff <- mapply( + FUN = function(cutoff, feature) { + return(ifelse( + test = is.na(x = cutoff), + yes = max(data[, feature]), + no = cutoff + )) + }, + cutoff = max.cutoff, + feature = features + ) + check.lengths <- unique(x = vapply( + X = list(features, min.cutoff, max.cutoff), + FUN = length, + FUN.VALUE = numeric(length = 1) + )) + if (length(x = check.lengths) != 1) { + stop("There must be the same number of minimum and maximum cuttoffs as there are features") + } + brewer.gran <- ifelse( + test = length(x = cols) == 1, + yes = brewer.pal.info[cols, ]$maxcolors, + no = length(x = cols) + ) + data[, 3:ncol(x = data)] <- sapply( + X = 3:ncol(x = data), + FUN = function(index) { + data.feature <- as.vector(x = data[, index]) + min.use <- SetQuantile(cutoff = min.cutoff[index - 2], data.feature) + max.use <- SetQuantile(cutoff = max.cutoff[index - 2], data.feature) + data.feature[data.feature < min.use] <- min.use + data.feature[data.feature > max.use] <- max.use + if (brewer.gran == 2) { + return(data.feature) + } + data.cut <- if (all(data.feature == 0)) { + 0 + } + else { + as.numeric(x = as.factor(x = cut( + x = as.numeric(x = data.feature), + breaks = brewer.gran + ))) + } + return(data.cut) + } + ) + colnames(x = data)[3:ncol(x = data)] <- features + rownames(x = data) <- cells + data$split <- if (is.null(x = split.by)) { + RandomName() + } else { + switch( + EXPR = split.by, + ident = Idents(object = object)[cells], + object[[split.by, drop = TRUE]][cells] + ) + } + if (!is.factor(x = data$split)) { + data$split <- factor(x = data$split) + } + if (!is.null(x = shape.by)) { + data[, shape.by] <- object[[shape.by, drop = TRUE]] + } + plots <- vector( + mode = "list", + length = ifelse( + test = blend, + yes = 4, + no = length(x = features) * length(x = levels(x = data$split)) + ) + ) + xlims <- ranges$x + ylims <- ranges$y + if (blend) { + ncol <- 4 + color.matrix <- BlendMatrix(col.threshold = blend.threshold) + colors <- list( + color.matrix[, 1], + color.matrix[1, ], + as.vector(x = color.matrix) + ) + } + for (i in 1:length(x = levels(x = data$split))) { + ident <- levels(x = data$split)[i] + data.plot <- data[as.character(x = data$split) == ident, , drop = FALSE] + if (blend) { + data.plot <- cbind(data.plot[, dims], BlendExpression(data = data.plot[, features[1:2]])) + features <- colnames(x = data.plot)[3:ncol(x = data.plot)] + } + for (j in 1:length(x = features)) { + feature <- features[j] + if (blend) { + cols.use <- as.numeric(x = as.character(x = data.plot[, feature])) + 1 + cols.use <- colors[[j]][sort(x = unique(x = cols.use))] + } else { + cols.use <- NULL + } + plot <- SingleDimPlot2( + data = data.plot[, c(dims, feature, shape.by)], + dims = dims, + col.by = feature, + order = order, + pt.size = pt.size, + cols = cols.use, + shape.by = shape.by, + label = label, + label.size = label.size + ) + + scale_x_continuous(limits = xlims) + + scale_y_continuous(limits = ylims) + + theme_minimal() + if (length(x = levels(x = data$split)) > 1) { + plot <- plot + theme(panel.border = element_rect(fill = NA, colour = 'black')) + plot <- plot + if (i == 1) { + labs(title = feature) + } else { + labs(title = NULL) + } + if (j == length(x = features) && !blend) { + suppressMessages( + expr = plot <- plot + + scale_y_continuous(sec.axis = dup_axis(name = ident)) + + no.right + ) + } + if (j != 1) { + plot <- plot + theme( + axis.line.y = element_blank(), + axis.ticks.y = element_blank(), + axis.text.y = element_blank(), + axis.title.y.left = element_blank() + ) + } + if (i != length(x = levels(x = data$split))) { + plot <- plot + theme( + axis.line.x = element_blank(), + axis.ticks.x = element_blank(), + axis.text.x = element_blank(), + axis.title.x = element_blank() + ) + } + } else { + plot <- plot + labs(title = feature) + } + if (!blend) { + plot <- plot + guides(color = NULL) + if (length(x = cols) == 1) { + plot <- plot + scale_color_brewer(palette = cols) + } else if (length(x = cols) > 1) { + plot <- suppressMessages( + expr = plot + scale_color_gradientn( + colors = cols, + guide = "colorbar" + ) + ) + } + } + if (coord.fixed) { + plot <- plot + coord_fixed() + } + plot <- plot + plots[[(length(x = features) * (i - 1)) + j]] <- plot + } + } + if (blend) { + blend.legend <- BlendMap(color.matrix = color.matrix) + for (i in 1:length(x = levels(x = data$split))) { + suppressMessages(expr = plots <- append( + x = plots, + values = list( + blend.legend + + scale_y_continuous( + sec.axis = dup_axis(name = ifelse( + test = length(x = levels(x = data$split)) > 1, + yes = levels(x = data$split)[i], + no = '' + )), + expand = c(0, 0) + ) + + labs( + x = features[1], + y = features[2], + title = if (i == 1) { + paste('Color threshold:', blend.threshold) + } else { + NULL + } + ) + + no.right + ), + after = 4 * i - 1 + )) + } + } + plots <- Filter(f = Negate(f = is.null), x = plots) + if (combine) { + if (is.null(x = ncol)) { + ncol <- 2 + if (length(x = features) == 1) { + ncol <- 1 + } + if (length(x = features) > 6) { + ncol <- 3 + } + if (length(x = features) > 9) { + ncol <- 4 + } + } + ncol <- ifelse( + test = is.null(x = split.by) || blend, + yes = ncol, + no = length(x = features) + ) + legend <- if (blend) { + 'none' + } else { + split.by %iff% 'none' + } + if (by.col & !is.null(x = split.by)) { + plots <- lapply(X = plots, FUN = function(x) { + suppressMessages(x + + theme_minimal() + ggtitle("") + + scale_y_continuous(sec.axis = dup_axis(name = "")) + no.right + ) + }) + nsplits <- length(x = levels(x = data$split)) + idx <- 1 + for (i in (length(x = features) * (nsplits - 1) + 1):(length(x = features) * nsplits)) { + plots[[i]] <- suppressMessages(plots[[i]] + scale_y_continuous(sec.axis = dup_axis(name = features[[idx]])) + no.right) + idx <- idx + 1 + } + idx <- 1 + for (i in which(x = 1:length(x = plots) %% length(x = features) == 1)) { + plots[[i]] <- plots[[i]] + ggtitle(levels(x = data$split)[[idx]]) + idx <- idx + 1 + } + idx <- 1 + if (length(x = features) == 1) { + for (i in 1:length(x = plots)) { + plots[[i]] <- plots[[i]] + ggtitle(levels(x = data$split)[[idx]]) + idx <- idx + 1 + } + } + plots <- plots[c(do.call( + what = rbind, + args = split(x = 1:length(x = plots), f = ceiling(x = seq_along(along.with = 1:length(x = plots))/length(x = features))) + ))] + plots <- CombinePlots( + plots = plots, + ncol = nsplits, + legend = legend + ) + } + else { + plots <- CombinePlots( + plots = plots[c(1,2,3)], + ncol = 3, + legend = legend, + nrow = split.by %iff% length(x = levels(x = data$split)) + ) + } + } + return(plots) +} + + +SingleExIPlot <- function( + data, + idents, + split = NULL, + type = 'violin', + sort = FALSE, + y.max = NULL, + adjust = 1, + pt.size = 0, + cols = NULL, + log = FALSE +) { + set.seed(seed = 42) + if (!is.data.frame(x = data) || ncol(x = data) != 1) { + stop("'SingleExIPlot requires a data frame with 1 column") + } + feature <- colnames(x = data) + data$ident <- idents + if ((is.character(x = sort) && nchar(x = sort) > 0) || sort) { + data$ident <- factor( + x = data$ident, + levels = names(x = rev(x = sort( + x = tapply( + X = data[, feature], + INDEX = data$ident, + FUN = mean + ), + decreasing = grepl(pattern = paste0('^', tolower(x = sort)), x = 'decreasing') + ))) + ) + } + if (log) { + noise <- rnorm(n = length(x = data[, feature])) / 200 + data[, feature] <- data[, feature] + 1 + } else { + noise <- rnorm(n = length(x = data[, feature])) / 100000 + } + if (all(data[, feature] == data[, feature][1])) { + warning(paste0("All cells have the same value of ", feature, ".")) + } else{ + data[, feature] <- data[, feature] + noise + } + axis.label <- ifelse(test = log, yes = 'Log Expression Level', no = 'Expression Level') + y.max <- y.max %||% max(data[, feature]) + if (is.null(x = split) || type != 'violin') { + vln.geom <- geom_violin + fill <- 'ident' + } else { + data$split <- split + vln.geom <- geom_split_violin + fill <- 'split' + } + switch( + EXPR = type, + 'violin' = { + x <- 'ident' + y <- paste0("`", feature, "`") + xlab <- 'Identity' + ylab <- axis.label + geom <- list( + vln.geom(scale = 'width', adjust = adjust, trim = TRUE), + theme(axis.text.x = element_text(angle = 45, hjust = 1)) + ) + jitter <- geom_jitter(height = 0, size = pt.size) + log.scale <- scale_y_log10() + axis.scale <- ylim + }, + 'ridge' = { + x <- paste0("`", feature, "`") + y <- 'ident' + xlab <- axis.label + ylab <- 'Identity' + geom <- list( + geom_density_ridges(scale = 4), + #theme_ridges(line_size = 0.001), + scale_y_discrete(expand = c(0.01, 0)), + scale_x_continuous(expand = c(0, 0)) + ) + jitter <- geom_jitter(width = 0, size = pt.size) + log.scale <- scale_x_log10() + axis.scale <- function(...) { + invisible(x = NULL) + } + }, + stop("Unknown plot type: ", type) + ) + plot <- ggplot( + data = data, + mapping = aes_string(x = x, y = y, fill = fill, color = fill)[c(3,4,2,1)] + ) + + labs(x = xlab, y = ylab, title = feature, fill = NULL) + + theme_minimal() + plot <- do.call(what = '+', args = list(plot, geom)) + plot <- plot + if (log) { + log.scale + } else { + axis.scale(min(data[, feature]), y.max) + } + if (pt.size > 0) { + plot <- plot + jitter + } + if (!is.null(x = cols)) { + if (!is.null(x = split)) { + idents <- unique(x = as.vector(x = data$ident)) + splits <- unique(x = as.vector(x = data$split)) + labels <- if (length(x = splits) == 2) { + splits + } else { + unlist(x = lapply( + X = idents, + FUN = function(pattern, x) { + x.mod <- gsub( + pattern = paste0(pattern, '.'), + replacement = paste0(pattern, ': '), + x = x, + fixed = TRUE + ) + x.keep <- grep(pattern = ': ', x = x.mod, fixed = TRUE) + x.return <- x.mod[x.keep] + names(x = x.return) <- x[x.keep] + return(x.return) + }, + x = unique(x = as.vector(x = data$split)) + )) + } + if (is.null(x = names(x = labels))) { + names(x = labels) <- labels + } + } else { + labels <- unique(x = as.vector(x = data$ident)) + } + plot <- plot + scale_fill_manual(values = cols, labels = labels) + scale_color_manual(values=rep("white",length(labels))) + } + return(plot) +} + +### DimPlot2 ---- +DimPlot2 <- function( + object, + dims = c(1, 2), + cells = NULL, + cols = NULL, + pt.size = NULL, + reduction = NULL, + group.by = NULL, + split.by = NULL, + shape.by = NULL, + order = NULL, + label = FALSE, + label.size = 4, + repel = FALSE, + cells.highlight = NULL, + cols.highlight = 'red', + sizes.highlight = 1, + na.value = 'grey50', + combine = TRUE, + ncol = NULL, + l1, + l2, + ... +) { + if (length(x = dims) != 2) { + stop("'dims' must be a two-length vector") + } + reduction <- reduction %||% { + default.reductions <- c('umap', 'tsne', 'pca') + object.reductions <- FilterObjects(object = object, classes.keep = 'DimReduc') + reduc.use <- min(which(x = default.reductions %in% object.reductions)) + default.reductions[reduc.use] + } + cells <- cells %||% colnames(x = object) + data <- Embeddings(object = object[[reduction]])[cells, dims] + data <- as.data.frame(x = data) + dims <- paste0(Key(object = object[[reduction]]), dims) + object[['ident']] <- Idents(object = object) + group.by <- group.by %||% 'ident' + data[, group.by] <- object[[group.by]][cells, , drop = FALSE] + for (group in group.by) { + if (!is.factor(x = data[, group])) { + data[, group] <- factor(x = data[, group]) + } + } + if (!is.null(x = shape.by)) { + data[, shape.by] <- object[[shape.by, drop = TRUE]] + } + if (!is.null(x = split.by)) { + data[, split.by] <- object[[split.by, drop = TRUE]] + } + plots <- lapply( + X = group.by, + FUN = function(x) { + plot <- SingleDimPlot2( + data = data[, c(dims, x, split.by, shape.by)], + dims = dims, + col.by = x, + cols = cols, + pt.size = pt.size, + shape.by = shape.by, + order = order, + label = FALSE, + cells.highlight = cells.highlight, + cols.highlight = cols.highlight, + sizes.highlight = sizes.highlight, + na.value = na.value + ) + if (label) { + plot <- LabelClusters2( + plot = plot, + id = x, + repel = repel, + size = label.size, + split.by = split.by, + l1 = l1, l2 = l2 + ) + } + if (!is.null(x = split.by)) { + #plot <- plot + FacetTheme() + + plot <- plot + theme_minimal() + + facet_wrap( + facets = vars(!!sym(x = split.by)), + ncol = if (length(x = group.by) > 1) { + length(x = unique(x = data[, split.by])) + } else { + NULL + } + ) + } + return(plot) + } + ) + if (combine) { + plots <- CombinePlots( + plots = plots, + ncol = if (!is.null(x = split.by) && length(x = group.by) > 1) { + 1 + } else { + ncol + }, + ... + ) + } + return(plots) +} + +FacetTheme <- function(...) { + return(theme( + strip.background = element_blank(), + strip.text = element_text(face = 'bold'), + # Validate the theme + validate = TRUE, + ... + )) +} + + +SetHighlight <- function( + cells.highlight, + cells.all, + sizes.highlight, + cols.highlight, + col.base = 'black', + pt.size = 1 +) { + if (is.character(x = cells.highlight)) { + cells.highlight <- list(cells.highlight) + } else if (is.data.frame(x = cells.highlight) || !is.list(x = cells.highlight)) { + cells.highlight <- as.list(x = cells.highlight) + } + cells.highlight <- lapply( + X = cells.highlight, + FUN = function(cells) { + cells.return <- if (is.character(x = cells)) { + cells[cells %in% cells.all] + } else { + cells <- as.numeric(x = cells) + cells <- cells[cells <= length(x = cells.all)] + cells.all[cells] + } + return(cells.return) + } + ) + cells.highlight <- Filter(f = length, x = cells.highlight) + names.highlight <- if (is.null(x = names(x = cells.highlight))) { + paste0('Group_', 1L:length(x = cells.highlight)) + } else { + names(x = cells.highlight) + } + sizes.highlight <- rep_len( + x = sizes.highlight, + length.out = length(x = cells.highlight) + ) + cols.highlight <- c( + col.base, + rep_len(x = cols.highlight, length.out = length(x = cells.highlight)) + ) + size <- rep_len(x = pt.size, length.out = length(x = cells.all)) + highlight <- rep_len(x = NA_character_, length.out = length(x = cells.all)) + if (length(x = cells.highlight) > 0) { + for (i in 1:length(x = cells.highlight)) { + cells.check <- cells.highlight[[i]] + index.check <- match(x = cells.check, cells.all) + highlight[index.check] <- names.highlight[i] + size[index.check] <- sizes.highlight[i] + } + } + plot.order <- sort(x = unique(x = highlight), na.last = TRUE) + plot.order[is.na(x = plot.order)] <- 'Unselected' + highlight[is.na(x = highlight)] <- 'Unselected' + highlight <- as.factor(x = highlight) + return(list( + plot.order = plot.order, + highlight = highlight, + size = size, + color = cols.highlight + )) +} + +AutoPointSize <- function(data) { + return(min(1583 / nrow(x = data), 1)) +} + +BlendMap <- function(color.matrix) { + color.heat <- matrix( + data = 1:prod(dim(x = color.matrix)) - 1, + nrow = nrow(x = color.matrix), + ncol = ncol(x = color.matrix), + dimnames = list( + 1:nrow(x = color.matrix), + 1:ncol(x = color.matrix) + ) + ) + xbreaks <- seq.int(from = 0, to = nrow(x = color.matrix), by = 2) + ybreaks <- seq.int(from = 0, to = ncol(x = color.matrix), by = 2) + color.heat <- Melt(x = color.heat) + color.heat$rows <- as.numeric(x = as.character(x = color.heat$rows)) + color.heat$cols <- as.numeric(x = as.character(x = color.heat$cols)) + color.heat$vals <- factor(x = color.heat$vals) + plot <- ggplot( + data = color.heat, + mapping = aes_string(x = 'rows', y = 'cols', fill = 'vals') + ) + + geom_raster(show.legend = FALSE) + + theme(plot.margin = unit(x = rep.int(x = 0, times = 4), units = 'cm')) + + scale_x_continuous(breaks = xbreaks, expand = c(0, 0), labels = xbreaks) + + scale_y_continuous(breaks = ybreaks, expand = c(0, 0), labels = ybreaks) + + scale_fill_manual(values = as.vector(x = color.matrix)) + + theme_cowplot() + return(plot) +} + + +FeaturePlot3 <- function( + object, + features, + dims = c(1, 2), + cells = NULL, + cols = c("lightgrey", "blue"), + pt.size = NULL, + order = FALSE, + min.cutoff = NA, + max.cutoff = NA, + reduction = NULL, + split.by = NULL, + shape.by = NULL, + blend = FALSE, + blend.threshold = 0.5, + label = FALSE, + label.size = 4, + ncol = NULL, + combine = TRUE, + coord.fixed = FALSE, + by.col = TRUE, + ranges=NULL +) { + no.right <- theme( + axis.line.y.right = element_blank(), + axis.ticks.y.right = element_blank(), + axis.text.y.right = element_blank(), + axis.title.y.right = element_text( + face = "bold", + size = 14, + margin = margin(r = 7) + ) + ) + if (is.null(reduction)) { + default_order <- c('umap', 'tsne', 'pca') + reducs <- which(default_order %in% names(object@reductions)) + reduction <- default_order[reducs[1]] + } + if (length(x = dims) != 2 || !is.numeric(x = dims)) { + stop("'dims' must be a two-length integer vector") + } + if (blend && length(x = features) != 2) { + stop("Blending feature plots only works with two features") + } + dims <- paste0(Key(object = object[[reduction]]), dims) + cells <- cells %||% colnames(x = object) + data <- FetchData(object = object, vars = c(dims, features), cells = cells) + features <- colnames(x = data)[3:ncol(x = data)] + min.cutoff <- mapply( + FUN = function(cutoff, feature) { + return(ifelse( + test = is.na(x = cutoff), + yes = min(data[, feature]), + no = cutoff + )) + }, + cutoff = min.cutoff, + feature = features + ) + max.cutoff <- mapply( + FUN = function(cutoff, feature) { + return(ifelse( + test = is.na(x = cutoff), + yes = max(data[, feature]), + no = cutoff + )) + }, + cutoff = max.cutoff, + feature = features + ) + check.lengths <- unique(x = vapply( + X = list(features, min.cutoff, max.cutoff), + FUN = length, + FUN.VALUE = numeric(length = 1) + )) + if (length(x = check.lengths) != 1) { + stop("There must be the same number of minimum and maximum cuttoffs as there are features") + } + brewer.gran <- ifelse( + test = length(x = cols) == 1, + yes = brewer.pal.info[cols, ]$maxcolors, + no = length(x = cols) + ) + data[, 3:ncol(x = data)] <- sapply( + X = 3:ncol(x = data), + FUN = function(index) { + data.feature <- as.vector(x = data[, index]) + min.use <- SetQuantile(cutoff = min.cutoff[index - 2], data.feature) + max.use <- SetQuantile(cutoff = max.cutoff[index - 2], data.feature) + data.feature[data.feature < min.use] <- min.use + data.feature[data.feature > max.use] <- max.use + if (brewer.gran == 2) { + return(data.feature) + } + data.cut <- if (all(data.feature == 0)) { + 0 + } + else { + as.numeric(x = as.factor(x = cut( + x = as.numeric(x = data.feature), + breaks = brewer.gran + ))) + } + return(data.cut) + } + ) + colnames(x = data)[3:ncol(x = data)] <- features + rownames(x = data) <- cells + data$split <- if (is.null(x = split.by)) { + RandomName() + } else { + switch( + EXPR = split.by, + ident = Idents(object = object)[cells], + object[[split.by, drop = TRUE]][cells] + ) + } + if (!is.factor(x = data$split)) { + data$split <- factor(x = data$split) + } + if (!is.null(x = shape.by)) { + data[, shape.by] <- object[[shape.by, drop = TRUE]] + } + plots <- vector( + mode = "list", + length = ifelse( + test = blend, + yes = 4, + no = length(x = features) * length(x = levels(x = data$split)) + ) + ) + xlims <- ranges$x + ylims <- ranges$y + if (blend) { + ncol <- 4 + color.matrix <- BlendMatrix(col.threshold = blend.threshold) + colors <- list( + color.matrix[, 1], + color.matrix[1, ], + as.vector(x = color.matrix) + ) + } + for (i in 1:length(x = levels(x = data$split))) { + ident <- levels(x = data$split)[i] + data.plot <- data[as.character(x = data$split) == ident, , drop = FALSE] + if (blend) { + data.plot <- cbind(data.plot[, dims], BlendExpression(data = data.plot[, features[1:2]])) + features <- colnames(x = data.plot)[3:ncol(x = data.plot)] + } + for (j in 1:length(x = features)) { + feature <- features[j] + if (blend) { + cols.use <- as.numeric(x = as.character(x = data.plot[, feature])) + 1 + cols.use <- colors[[j]][sort(x = unique(x = cols.use))] + } else { + cols.use <- NULL + } + plot <- SingleDimPlot2( + data = data.plot[, c(dims, feature, shape.by)], + dims = dims, + col.by = feature, + order = order, + pt.size = pt.size, + cols = cols.use, + shape.by = shape.by, + label = label, + label.size = label.size + ) + + scale_x_continuous(limits = xlims) + + scale_y_continuous(limits = ylims) + + theme_minimal() + theme(legend.position = c(0.95, 0.9)) ## Replace legend to have same scale + if (length(x = levels(x = data$split)) > 1) { + plot <- plot + theme(panel.border = element_rect(fill = NA, colour = 'black')) + plot <- plot + if (i == 1) { + labs(title = feature) + } else { + labs(title = NULL) + } + if (j == length(x = features) && !blend) { + suppressMessages( + expr = plot <- plot + + scale_y_continuous(sec.axis = dup_axis(name = ident)) + + no.right + ) + } + if (j != 1) { + plot <- plot + theme( + axis.line.y = element_blank(), + axis.ticks.y = element_blank(), + axis.text.y = element_blank(), + axis.title.y.left = element_blank() + ) + } + if (i != length(x = levels(x = data$split))) { + plot <- plot + theme( + axis.line.x = element_blank(), + axis.ticks.x = element_blank(), + axis.text.x = element_blank(), + axis.title.x = element_blank() + ) + } + } else { + plot <- plot + labs(title = feature) + } + if (!blend) { + plot <- plot + guides(color = NULL) + if (length(x = cols) == 1) { + plot <- plot + scale_color_brewer(palette = cols) + } else if (length(x = cols) > 1) { + plot <- suppressMessages( + expr = plot + scale_color_gradientn( + colors = cols, + guide = "colorbar" + ) + ) + } + } + if (coord.fixed) { + plot <- plot + coord_fixed() + } + plot <- plot + plots[[(length(x = features) * (i - 1)) + j]] <- plot + } + } + if (blend) { + blend.legend <- BlendMap(color.matrix = color.matrix) + for (i in 1:length(x = levels(x = data$split))) { + suppressMessages(expr = plots <- append( + x = plots, + values = list( + blend.legend + + scale_y_continuous( + sec.axis = dup_axis(name = ifelse( + test = length(x = levels(x = data$split)) > 1, + yes = levels(x = data$split)[i], + no = '' + )), + expand = c(0, 0) + ) + + labs( + x = features[1], + y = features[2], + title = if (i == 1) { + paste('Color threshold:', blend.threshold) + } else { + NULL + } + ) + + no.right + ), + after = 4 * i - 1 + )) + } + } + plots <- Filter(f = Negate(f = is.null), x = plots) + if (combine) { + if (is.null(x = ncol)) { + ncol <- 2 + if (length(x = features) == 1) { + ncol <- 1 + } + if (length(x = features) > 6) { + ncol <- 3 + } + if (length(x = features) > 9) { + ncol <- 4 + } + } + ncol <- ifelse( + test = is.null(x = split.by) || blend, + yes = ncol, + no = length(x = features) + ) + legend <- if (blend) { + 'none' + } else { + split.by %iff% 'none' + } + if (by.col & !is.null(x = split.by)) { + plots <- lapply(X = plots, FUN = function(x) { + suppressMessages(x + + theme_minimal() + ggtitle("") + + scale_y_continuous(sec.axis = dup_axis(name = "")) + no.right + ) + }) + nsplits <- length(x = levels(x = data$split)) + idx <- 1 + for (i in (length(x = features) * (nsplits - 1) + 1):(length(x = features) * nsplits)) { + plots[[i]] <- suppressMessages(plots[[i]] + scale_y_continuous(sec.axis = dup_axis(name = features[[idx]])) + no.right) + idx <- idx + 1 + } + idx <- 1 + for (i in which(x = 1:length(x = plots) %% length(x = features) == 1)) { + plots[[i]] <- plots[[i]] + ggtitle(levels(x = data$split)[[idx]]) + idx <- idx + 1 + } + idx <- 1 + if (length(x = features) == 1) { + for (i in 1:length(x = plots)) { + plots[[i]] <- plots[[i]] + ggtitle(levels(x = data$split)[[idx]]) + idx <- idx + 1 + } + } + plots <- plots[c(do.call( + what = rbind, + args = split(x = 1:length(x = plots), f = ceiling(x = seq_along(along.with = 1:length(x = plots))/length(x = features))) + ))] + plots <- CombinePlots( + plots = plots, + ncol = nsplits, + legend = legend + ) + } + else { + plots <- CombinePlots( + plots = plots, + ncol = ncol, + legend = legend, + nrow = split.by %iff% length(x = levels(x = data$split)) + ) + } + } + return(plots) +} diff --git a/server.R b/server.R new file mode 100644 index 0000000..3e89f05 --- /dev/null +++ b/server.R @@ -0,0 +1,1393 @@ +#### CenGenAPP 2018-2020 +#### Report bugs to: gabrielsantperebaro@gmail.com + + +library(BiocManager) +require(shiny) +library(shinyjs) +library(shinythemes) +library(ggridges) +library(DT) +library(cowplot) +require(dplyr) +require(Seurat) +library(ggplot2) +library(tidyverse) +library(pheatmap) +library(expss) +library(plotly) +#library(cairo) +#library(MAST) +options(repos = BiocManager::repositories()) +#setwd("/Users/gabrielsantpere/Documents/Scripts/RshinyDeveloping_April2020") +#source("Functions.R") + + +utr <- + c( + "WBGene00023498", + "WBGene00023497", + "WBGene00004397", + "WBGene00006843", + "WBGene00004010", + "WBGene00006789", + "WBGene00001135", + "WBGene00001079", + "WBGene00001135", + "WBGene00006783", + "WBGene00000501", + "WBGene00006788", + "WBGene00001555" + ) + +server <- function(input, output) { + ### Plot features ---- + observeEvent(input$PlotButton, { + withProgress(message = "Generating Plots...", + value = 0, + expr = { + if (exists("id1")) { + removeNotification(id1) + } + + if (input$dataset == "Neurons only") { + all <- allNeurons + } + if (input$dataset == "All cell types") { + all <- allCells + } + Idents(object = all) <- "Neuron.type" + ranges <- reactiveValues(x = NULL, y = NULL) + observeEvent(input$plotPermanent_dblclick, { + brush <- input$plotPermanent_brush + if (!is.null(brush)) { + ranges$x <- c(brush$xmin, brush$xmax) + ranges$y <- c(brush$ymin, brush$ymax) + } else { + ranges$x <- NULL + ranges$y <- NULL + } + }) + + + plots <- input$Plots + if (plots == 'Metadata') { + varToPlot <- input$featuresMetadata + featuresInput <- function() { + if (is.null(all)) { + return(NULL) + } else{ + withProgress(message = "Generating Features Plot...", value = 0, { + DimPlot2( + object = all, + reduction = "umap", + group.by = 'Neuron.type', + label = TRUE, + pt.size = 0.7, + repel = TRUE, + l1 = ranges$x, + l2 = ranges$y + ) + theme(legend.position = "none") + coord_cartesian(xlim = ranges$x, + ylim = ranges$y, + expand = FALSE) + }) + } + } + + if (input$featuresMetadata == 'Neuron.type' && + input$dataset == "Neurons only") { + Type <- input$Cluster + output$SinglePlot <- + renderPlot({ + featuresInput() + }, height = 400, width = 600) + fn1 <- + function() { + withProgress(message = "Generating Features Plot...", value = 0, { + DimPlot2( + object = all, + reduction = "umap", + group.by = "Neuron.type", + cells.highlight = names(all$Neuron.type[all$Neuron.type == Type]), + cols.highlight = "red", + pt.size = 0.7, + l1 = ranges$x, + l2 = ranges$y + ) + theme(legend.position = "none") + coord_cartesian(xlim = ranges$x, + ylim = ranges$y, + expand = FALSE) + }) + } + output$High <- + renderPlot({ + fn1() + }, height = 400, width = 600, execOnResize = FALSE) + output$downloadHigh <- + downloadHandler( + filename = function() { + filename = paste("NeuronType-", Type, ".png", sep = "") + }, + content = function(file) { + ggsave( + fn1(), + file = file, + height = 100, + width = 150, + units = "mm" , + limitsize = FALSE, + device = "png" + ) + } + ) + } + if (input$featuresMetadata == 'Detection' && + input$dataset == "Neurons only") { + Det <- input$Cluster3 + output$SinglePlot <- + renderPlot({ + featuresInput() + }, height = 400, width = 600) + fn3 <- + function() { + withProgress(message = "Generating Features Plot...", value = 0, { + DimPlot2( + object = all, + reduction = "umap", + group.by = "Detection" , + cells.highlight = names(all$Detection[all$Detection == Det]), + cols.highlight = "red", + pt.size = 0.7 + ) + theme(legend.position = "none") + coord_cartesian(xlim = ranges$x, + ylim = ranges$y, + expand = FALSE) + }) + } + output$High <- + renderPlot({ + fn3() + }, height = 400, width = 600, execOnResize = FALSE) + output$downloadHigh <- + downloadHandler( + filename = function() { + filename = paste("Detection-", Det, ".png", sep = "") + }, + content = function(file) { + ggsave( + fn3(), + file = file, + height = 100, + width = 150, + units = "mm" , + limitsize = FALSE, + device = "png" + ) + } + ) + } + if (input$featuresMetadata == 'Experiment' && + input$dataset == "Neurons only") { + Exp <- input$Cluster2 + output$SinglePlot <- + renderPlot({ + featuresInput() + }, height = 400, width = 600) + fn2 <- + function() { + withProgress(message = "Generating Features Plot...", value = 0, { + DimPlot2( + object = all, + reduction = "umap", + group.by = "Experiment", + cells.highlight = names(all$Experiment[all$Experiment == Exp]), + cols.highlight = "red", + pt.size = 0.7 + ) + theme(legend.position = "none") + coord_cartesian(xlim = ranges$x, + ylim = ranges$y, + expand = FALSE) + }) + } + output$High <- + renderPlot({ + fn2() + }, height = 400, width = 600, execOnResize = FALSE) + output$downloadHigh <- + downloadHandler( + filename = function() { + filename = paste("Experiment-", Exp, ".png", sep = "") + }, + content = function(file) { + ggsave( + fn2(), + file = file, + height = 100, + width = 150, + units = "mm" , + limitsize = FALSE, + device = "png" + ) + } + ) + } + if (input$featuresMetadata2 == 'Neuron.type' && + input$dataset == "All cell types") { + Type <- input$Cluster.2 + output$SinglePlot <- + renderPlot({ + featuresInput() + }, height = 400, width = 600) + fn1 <- + function() { + withProgress(message = "Generating Features Plot...", value = 0, { + DimPlot2( + object = all, + reduction = "umap", + group.by = "Neuron.type", + cells.highlight = names(all$Neuron.type[all$Neuron.type == Type]), + cols.highlight = "red", + pt.size = 0.7, + l1 = ranges$x, + l2 = ranges$y + ) + theme(legend.position = "none") + coord_cartesian(xlim = ranges$x, + ylim = ranges$y, + expand = FALSE) + }) + } + output$High <- + renderPlot({ + fn1() + }, height = 400, width = 600, execOnResize = FALSE) + output$downloadHigh <- + downloadHandler( + filename = function() { + filename = paste("NeuronType-", Type, ".png", sep = "") + }, + content = function(file) { + ggsave( + fn1(), + file = file, + height = 100, + width = 150, + units = "mm" , + limitsize = FALSE, + device = "png" + ) + } + ) + } + if (input$featuresMetadata2 == 'Detection' && + input$dataset == "All cell types") { + Det <- input$Cluster3.2 + output$SinglePlot <- + renderPlot({ + featuresInput() + }, height = 400, width = 600) + fn3 <- + function() { + withProgress(message = "Generating Features Plot...", value = 0, { + DimPlot2( + object = all, + reduction = "umap", + group.by = "Detection" , + cells.highlight = names(all$Detection[all$Detection == Det]), + cols.highlight = "red", + pt.size = 0.7 + ) + theme(legend.position = "none") + coord_cartesian(xlim = ranges$x, + ylim = ranges$y, + expand = FALSE) + }) + } + output$High <- + renderPlot({ + fn3() + }, height = 400, width = 600, execOnResize = FALSE) + output$downloadHigh <- + downloadHandler( + filename = function() { + filename = paste("Detection-", Det, ".png", sep = "") + }, + content = function(file) { + ggsave( + fn3(), + file = file, + height = 100, + width = 150, + units = "mm" , + limitsize = FALSE, + device = "png" + ) + } + ) + } + if (input$featuresMetadata2 == 'Experiment' && + input$dataset == "All cell types") { + Exp <- input$Cluster2.2 + output$SinglePlot <- + renderPlot({ + featuresInput() + }, height = 400, width = 600) + fn2 <- + function() { + withProgress(message = "Generating Features Plot...", value = 0, { + DimPlot2( + object = all, + reduction = "umap", + group.by = "Experiment", + cells.highlight = names(all$Experiment[all$Experiment == Exp]), + cols.highlight = "red", + pt.size = 0.7 + ) + theme(legend.position = "none") + coord_cartesian(xlim = ranges$x, + ylim = ranges$y, + expand = FALSE) + }) + } + output$High <- + renderPlot({ + fn2() + }, height = 400, width = 600, execOnResize = FALSE) + output$downloadHigh <- + downloadHandler( + filename = function() { + filename = paste("Experiment-", Exp, ".png", sep = "") + }, + content = function(file) { + ggsave( + fn2(), + file = file, + height = 100, + width = 150, + units = "mm" , + limitsize = FALSE, + device = "png" + ) + } + ) + } + if (input$featuresMetadata2 == 'Tissue.type' && + input$dataset == "All cell types") { + output$SinglePlot <- + renderPlot({ + featuresInput() + }, height = 400, width = 600) + fn4 <- + function() { + withProgress(message = "Generating Features Plot...", value = 0, { + DimPlot2( + object = all, + reduction = "umap", + group.by = "Tissue.type", + cells.highlight = names(all$Tissue.type[all$Tissue.type == input$Cluster4]), + cols.highlight = "red", + pt.size = 0.7 + ) + theme(legend.position = "none") + coord_cartesian(xlim = ranges$x, + ylim = ranges$y, + expand = FALSE) + }) + } + output$High <- + renderPlot({ + fn4() + }, height = 400, width = 600, execOnResize = FALSE) + output$downloadHigh <- + downloadHandler( + filename = function() { + filename = paste("Tissue-", input$Cluster4, ".png", sep = "") + }, + content = function(file) { + ggsave( + fn4(), + file = file, + height = 100, + width = 150, + units = "mm" , + limitsize = FALSE, + device = "png" + ) + } + ) + } + + output$downloadClust <- + downloadHandler( + filename = function() { + filename = paste("Clusters.png", sep = "") + }, + content = function(file) { + ggsave( + featuresInput(), + file = file, + height = 200, + width = 300, + units = "mm" , + limitsize = FALSE, + device = "png" + ) + } + ) + } + + ### Plot genes individual ---- + if (plots == 'Individual genes') { + g <- input$GeneName + print(g) + + if (g %in% gene_list$gene_name) { + varToPlot <- filter(gene_list, gene_name == g)$gene_id + } + + if (g %in% gene_list$gene_id) { + varToPlot <- g + } + + if (g %in% gene_list$seqnames) { + varToPlot <- filter(gene_list, seqnames == g)$gene_id + } + + output$geneUTR <- + isolate(renderText({ + shiny::validate(need( + !varToPlot %in% utr, + message = paste0( + "WARNING: ", + g, + " expression is unreliable as it has been overexpressed to generate transgenic strains." + ) + )) + })) + + + f <- function() { + DimPlot2( + object = all, + reduction = "umap", + group.by = "Neuron.type", + label = TRUE, + repel = TRUE, + l1 = ranges$x, + l2 = ranges$y + ) + theme(legend.position = "none") + coord_cartesian(xlim = ranges$x, + ylim = ranges$y, + expand = TRUE) + } + f2 <- function() { + RidgePlot( + object = all, + features = as.character(varToPlot), + ncol = 1, + group.by = "Neuron.type", + sort = "decreasing" + ) + theme( + legend.position = "none", + panel.grid = element_blank(), + panel.background = element_blank() + ) + } + f3 <- + function() { + FeaturePlot3( + object = all, + features = as.character(varToPlot), + reduction = "umap", + pt.size = 0.7, + combine = TRUE, + cols = c("beige", "darkred"), + coord.fixed = FALSE, + ranges = ranges + ) + } + + output$SinglePlot2 <- renderPlot({ + shiny::validate( + need(as.character(g) %in% gene_list$gene_name | as.character(g) %in% gene_list$gene_id | as.character(g) %in% gene_list$seqnames, message = "Gene name is not present in this dataset or wrong name") + ) + f() + }, height = 400, width = 600, execOnResize = FALSE) + output$FeaturePlot <- renderPlot({ + shiny::validate( + need(as.character(g) %in% gene_list$gene_name | as.character(g) %in% gene_list$gene_id | as.character(g) %in% gene_list$seqnames, message = "Gene name is not present in this dataset or wrong name") + ) + f3() + }, height = 400, width = 600, execOnResize = FALSE) + output$Violin <- renderPlot({ + shiny::validate( + need(as.character(g) %in% gene_list$gene_name | as.character(g) %in% gene_list$gene_id | as.character(g) %in% gene_list$seqnames, message = "Gene name is not present in this dataset or wrong name") + ) + f2() + }, height = 1300, width = 600) + + + output$downloadExp <- + downloadHandler( + filename = function() { + filename = paste("Dim-", varToPlot, ".png", sep = "") + }, + content = function(file) { + ggsave( + f3(), + file = file, + height = 100, + width = 150, + units = "mm" , + limitsize = FALSE, + device = "png" + ) + } + ) + output$downloadViolin <- + downloadHandler( + filename = function() { + filename = paste("Violin-", varToPlot, ".png", sep = "") + }, + content = function(file) { + ggsave( + f2(), + file = file, + height = 1300, + width = 600, + units = "mm" , + limitsize = FALSE, + device = "png" + ) + } + ) + } + + ### Plot genes colocoalization ---- + + if (plots == 'Colocalization') { + g1 <- input$ColocalizationGenes1 + g2 <- input$ColocalizationGenes2 + #gene1 <- filter(gene_list, gene_name == g1)$gene_id + #gene2 <- filter(gene_list, gene_name == g2)$gene_id + + if (g1 %in% gene_list$gene_name) { + varToPlot1 <- filter(gene_list, gene_name == g1)$gene_id + } + + if (g1 %in% gene_list$gene_id) { + varToPlot1 <- g1 + } + + if (g1 %in% gene_list$seqnames) { + varToPlot1 <- filter(gene_list, seqnames == g1)$gene_id + } + + if (g2 %in% gene_list$gene_name) { + varToPlot2 <- filter(gene_list, gene_name == g2)$gene_id + } + + if (g2 %in% gene_list$gene_id) { + varToPlot2 <- g2 + } + + if (g2 %in% gene_list$seqnames) { + varToPlot2 <- filter(gene_list, seqnames == g2)$gene_id + } + + + featuresInput <- function() { + if (is.null(all)) { + return(NULL) + } else{ + varBlend <- input$blend + withProgress(message = "Generating Features Plot...", value = + 0, { + FeaturePlot2( + object = all, + features = c( + #as.character(filter(gene_list, gene_name == g1)$gene_id), + #as.character(filter(gene_list, gene_name == g2)$gene_id) + varToPlot1, varToPlot2 + ), + reduction = "umap", + pt.size = 0.7, + combine = TRUE, + cols = c("beige", "darkred"), + blend = TRUE, + blend.threshold = varBlend, + coord.fixed = TRUE, + ranges = ranges + ) + }) + } + } + + output$SinglePlotDouble <- isolate(renderPlot({ + shiny::validate( + need(as.character(g1) %in% gene_list$gene_name | as.character(g1) %in% gene_list$gene_id | as.character(g1) %in% gene_list$seqnames, message = "Gene name 1 is not present in this dataset or wrong name") + ) + shiny::validate( + need(as.character(g2) %in% gene_list$gene_name | as.character(g2) %in% gene_list$gene_id | as.character(g2) %in% gene_list$seqnames, message = "Gene name 2 is not present in this dataset or wrong name") + ) + featuresInput() + }, height = 400, width = 1000)) + output$gene1utr <- + renderText({ + shiny::validate(need( + !varToPlot1 %in% utr , + message = paste0( + "WARNING: ", + g1, + " expression is unreliable as it has been overexpressed to generate transgenic strains." + ) + )) + }) + output$gene2utr <- + renderText({ + shiny::validate(need( + !varToPlot2 %in% utr , + message = paste0( + "WARNING: ", + g2, + " expression is unreliable as it has been overexpressed to generate transgenic strains." + ) + )) + }) + output$downloadCol <- + downloadHandler( + filename = function() { + filename = paste("Colocalization-", g1, "_", g2, ".png", sep = "") + }, + content = function(file) { + ggsave( + featuresInput(), + file = file, + height = 200, + width = 500, + units = "mm" , + limitsize = FALSE, + device = "png" + ) + } + ) + output$SinglePlotPermanent <- + renderPlot( + DimPlot2( + object = all, + reduction = "umap", + group.by = "Neuron.type", + label = TRUE, + repel = TRUE, + l1 = ranges$x, + l2 = ranges$y + ) + theme(legend.position = "none") + coord_cartesian( + xlim = ranges$x, + ylim = ranges$y, + expand = FALSE + ), + height = 400, + width = 600, + execOnResize = FALSE + ) + + } + }) + + }) + + ### Tables of markers ---- + markers$p_val <- + formatC(markers$p_val, format = "e", digits = 3) %>% gsub(" ", "", .) + markers$p_val_adj <- + formatC(markers$p_val_adj, format = "e", digits = 3) %>% gsub(" ", "", .) + markers$avg_logFC <- + formatC(markers$avg_logFC, digits = 3) %>% gsub(" ", "", .) + markersAllcells$p_val <- + formatC(markersAllcells$p_val, + format = "e", + digits = 3) %>% gsub(" ", "", .) + markersAllcells$p_val_adj <- + formatC(markersAllcells$p_val_adj, + format = "e", + digits = 3) %>% gsub(" ", "", .) + markersAllcells$avg_logFC <- + formatC(markersAllcells$avg_logFC, digits = 3) %>% gsub(" ", "", .) + + observeEvent(input$Markers, { + print(input$top2) + output$MarkTable <- DT::renderDataTable({ + DT::datatable( + filter(markers, cluster == input$Markers) %>% arrange(p_val_adj, desc(avg_logFC)) %>% head(as.numeric(input$top)) , + style = 'jQueryUI', + class = 'cell-border stripe', + rownames = FALSE + ) %>% formatStyle(c(1:8), color = "black") + }) + t1 <- + filter(markers, cluster == input$Markers) %>% arrange(p_val_adj, desc(avg_logFC)) + output$downloadMarkers <- + downloadHandler( + filename = function() { + paste("CellMarkers_NeuronsOnly-", input$Markers, ".csv", sep = "") + }, + content = function(file) { + write.csv(t1, file, sep = "\t") + } + ) + }) + observeEvent(input$Markers2, { + output$MarkTable2 <- DT::renderDataTable({ + DT::datatable( + filter(markersAllcells, cluster == input$Markers2) %>% arrange(p_val_adj, desc(avg_logFC)) %>% head(as.numeric(input$top2)), + style = 'jQueryUI', + class = 'cell-border stripe', + rownames = FALSE + ) %>% formatStyle(c(1:8), color = "black") + }) + t2 <- + filter(markersAllcells, cluster == input$Markers2) %>% arrange(p_val_adj, desc(avg_logFC)) + output$downloadMarkers2 <- + downloadHandler( + filename = function() { + paste("CellMarkers_AllCells-", input$Markers2, ".csv", sep = "") + }, + content = function(file) { + write.csv(t2, file, sep = "\t") + } + ) + }) + + + ### Tables of DEX ---- + observeEvent(input$DEXButton, { + print(input$ClusterCells1) + withProgress( + message = paste0( + 'Calculating DEX between ', + input$ClusterCells1, + ' and ', + input$ClusterCells2, + " using ", + input$test + ), + value = 0, + { + Idents(object = allCells) <- "Neuron.type" + tableDEX <- + FindMarkers( + allCells, + ident.1 = input$ClusterCells1, + ident.2 = input$ClusterCells2, + test.use = input$test + ) + tableDEX$gene <- rownames(tableDEX) + tableDEX <- + merge( + tableDEX, + gene_list, + by.x = "gene", + by.y = "gene_id", + all.x = TRUE + ) + tableDEX <- tableDEX %>% arrange(p_val_adj) + if (input$test != "roc") { + tableDEX$p_val <- + as.numeric(formatC(tableDEX$p_val, format = "e", digits = 3) %>% gsub(" ", "", .)) + tableDEX$p_val_adj <- + as.numeric(formatC( + tableDEX$p_val_adj, + format = "e", + digits = 3 + ) %>% gsub(" ", "", .)) + tableDEX$avg_logFC <- + as.numeric(formatC(tableDEX$avg_logFC, digits = 3) %>% gsub(" ", "", .)) + } + } + ) + if (nrow(tableDEX) > 0) { + output$MarkTable_ClusterCells <- DT::renderDataTable({ + DT::datatable( + tableDEX %>% head(input$topM2), + options = list(pageLength = input$topM2), + style = 'jQueryUI', + class = 'cell-border stripe', + rownames = FALSE + ) %>% formatStyle(c(1:8), color = "black") + }) + output$downloadDEX <- + downloadHandler( + filename = function() { + paste("DEXGens-", + input$ClusterCells1, + "", + input$ClusterCells2, + ".csv", + sep = "") + }, + content = function(file) { + write.csv(tableDEX, file, sep = "\t") + } + ) + } else { + output$MarkTable_ClusterCells <- + "No features pass logfc.threshold threshold" + } + }, ignoreNULL = TRUE) + + observeEvent(input$DEXButtonBatch, { + print(input$batch1) + print(input$batch2) + Idents(object = allCells) <- "Neuron.type" + b1 <- unlist(strsplit(input$batch1, split = ",")) + b1 <- gsub(" ", "", as.character(b1)) + b2 <- unlist(strsplit(input$batch2, split = ",")) + b2 <- gsub(" ", "", as.character(b2)) + + if (mean(c(b1, b2) %in% allCells$Neuron.type) == 1 | b2 == "ALL" | b2 == "NEURONS") { + output$text2 <- renderText({ + "" + }) + + if (!b2 %in% c("ALL", "NEURONS")){ + + tableDEX <- + FindMarkers( + allCells, + #ident.1 = unlist(strsplit(input$batch1, split = ",| ")), + #ident.2 = unlist(strsplit(input$batch2, split = ",| ")), + ident.1 = b1, + ident.2 = b2, + test.use = input$test + ) + tableDEX$gene <- rownames(tableDEX) + + } + + if (b2 == "ALL"){ + print(b2) + tableDEX <- + FindMarkers( + allCells, + #ident.1 = unlist(strsplit(input$batch1, split = ",| ")), + ident.1 = b1, + test.use = input$test + ) + + tableDEX$gene <- rownames(tableDEX) + } + + if (b2 == "NEURONS"){ + print(b2) + tableDEX <- + FindMarkers( + allCells, + #ident.1 = unlist(strsplit(input$batch1, split = ",| ")), + ident.1 = b1, + ident.2 = filter(allCells@meta.data, !Neuron.type %in% b1, Tissue.type == "Neuron")$Neuron.type %>% unique, + test.use = input$test + ) + + tableDEX$gene <- rownames(tableDEX) + } + + + tableDEX <- + merge( + tableDEX, + gene_list, + by.x = "gene", + by.y = "gene_id", + all.x = TRUE + ) + tableDEX <- tableDEX %>% arrange(p_val_adj) + if (input$test != "roc") { + tableDEX$p_val <- + as.numeric(formatC(tableDEX$p_val, format = "e", digits = 3) %>% gsub(" ", "", .)) + tableDEX$p_val_adj <- + as.numeric(formatC( + tableDEX$p_val_adj, + format = "e", + digits = 3 + ) %>% gsub(" ", "", .)) + tableDEX$avg_logFC <- + as.numeric(formatC(tableDEX$avg_logFC, digits = 3) %>% gsub(" ", "", .)) + } + if (nrow(tableDEX) > 0) { + output$MarkTable_Batch <- DT::renderDataTable({ + DT::datatable( + tableDEX %>% head(input$topM2), + options = list(pageLength = input$topM2), + style = 'jQueryUI', + class = 'cell-border stripe', + rownames = FALSE + ) %>% formatStyle(c(1:8), color = "black") + }) + output$downloadDEX <- + downloadHandler( + filename = function() { + paste( + "DEXGens-", + paste(b1, collapse = ","), + "-", + paste(b2, collapse = ","), + ".csv", + sep = "" + ) + }, + content = function(file) { + write.csv(tableDEX, file, sep = "\t") + } + ) + } else { + output$MarkTable_Batch <- + "No features pass logfc.threshold threshold" + } + } else { + output$MarkTable_Batch <- DT::renderDataTable({ + "NULL" + }) + output$text2 <- + renderText({ + "One or more cell types introduced are not correct" + }) + } + }, ignoreNULL = TRUE) + + ### React to thresholds ---- + + observeEvent(input$TCell, { + withProgress(message = "Obtaining information...", value = 0, { + print(input$Tcell_cut) + t4 <- + filter(ths, threshold == input$Tcell_cut) %>% dplyr::select(gene_name, input$Tcell_name) + t4d <- + filter(ths, threshold == input$Tcell_cut) %>% dplyr::select(gene_name, X, input$Tcell_name) + t4 <- t4[rev(order(t4[, 2])), ] + t4d <- t4d[rev(order(t4d[, 3])), ] + t4[, 2] <- + as.numeric(formatC(t4[, 2], digits = 3, format = "f") %>% gsub(" ", "", .)) + colnames(t4) <- c("Gene name", "Expression level") + t4d[, 3] <- + as.numeric(formatC(t4d[, 3], digits = 3, format = "f") %>% gsub(" ", "", .)) + colnames(t4d) <- c("Gene name", "Gene ID", "Expression level") + output$Tcell_name_table <- + DT::renderDataTable({ + DT::datatable( + t4, + options = list(pageLength = 10, autoWidth = TRUE), + rownames = FALSE, + style = 'jQueryUI', + class = 'cell-border stripe' + ) %>% formatStyle(c(1:2), color = "black") + }) + output$get_download_gene <- renderUI({ + req(input$TCell) + downloadButton('downloadGene', "Download table") + }) + + output$downloadGene <- + downloadHandler( + filename = function() { + paste( + "GenesExpressed_in_", + input$Tcell_name, + "-thrs", + input$Tcell_cut, + ".csv", + sep = "" + ) + }, + content = function(file) { + write.csv(t4d , file, dec = ".", sep = "\t") + } + ) + }) + }) + + observeEvent(input$TGene, { + + g <- input$Tgene_name + print(g) + + var <- "" + + if (g %in% gene_list$gene_name) { + var <- g + } + + if (g %in% gene_list$gene_id) { + var <- filter(gene_list, gene_id == g)$gene_name + } + + if (g %in% gene_list$seqnames) { + var <- filter(gene_list, seqnames == g)$gene_name + } + + withProgress(message = "Obtaining information...", value = 0, { + if (var %in% filter(ths, threshold == input$Tgene_cut)$gene_name) { + print("YES") + t3 <- + sort( + filter( + ths, + gene_name == var, + threshold == input$Tgene_cut + )[, 2:129], + decreasing = TRUE + ) + t3 <- data.frame(CellType = names(t3), expression = as.numeric(t3)) + t3 <- dplyr::filter(t3, expression > 0) + if (nrow(t3) > 0) { + t3[, 2] <- + as.numeric(formatC(t3[, 2], digits = 3, format = "f") %>% gsub(" ", "", .)) + } + colnames(t3) <- c("Cell type", "Expression level") + output$text1 <- + isolate(renderText({ + shiny::validate(need( + !filter(gene_list, gene_name == var)$gene_id %in% utr, + message = paste0( + "WARNING: ", + input$Tgene_name, + " expression is unreliable as it has been overexpressed to generate transgenic strains." + ) + )) + })) + #output$text1 <- renderText({""}) + } + else { + print("NO") + t3 <- NULL + output$text1 <- + renderText({ + "Gene is not expressed at any threshold or does not exist" + }) + } + + output$Tgene_name_table <- + DT::renderDataTable({ + DT::datatable( + t3, + options = list(pageLength = 10, autoWidth = TRUE), + rownames = FALSE, + style = 'jQueryUI', + class = 'cell-border stripe' + ) %>% formatStyle(c(1:2), color = "black") + }) + + output$get_download_cell <- renderUI({ + req(input$Tgene_name) + downloadButton('downloadCell', "Download table") + }) + + output$downloadCell <- + downloadHandler( + filename = function() { + paste( + "GenesExpressing-", + input$Tgene_name, + "-thrs", + input$Tgene_cut, + ".csv", + sep = "" + ) + }, + content = function(file) { + write.csv(t3, file, dec = ".", sep = "\t") + } + ) + }) + }) + + observeEvent(input$Tgene_name_batch, { + + + gns1 <- unlist(strsplit(input$Tgene_name_batch, split = ",")) + gns1 <- gsub(" ", "", as.character(gns1)) + gns <- unique(c(gns1, filter(gene_list, gene_id %in% gns1 | seqnames %in% gns1)$gene_name)) + print(gns) + req(input$Tgene_name_batch) + + if (length(which(gns %in% unique(ths$gene_name))) == length(gns1)) { + print("YES") + tb <- + filter(ths, + gene_name %in% gns, + threshold == input$Tgene_cut_batch)[, c(1, 131, 130, c(2:129))] + output$textb <- renderText({ + "" + }) + output$TGeneBatch <- + downloadHandler( + filename = function() { + paste( + "GenesExpressing-", + paste(gns1, collapse = ","), + "-thrs", + input$Tgene_cut_batch, + ".csv", + sep = "" + ) + }, + content = function(file) { + write.csv(tb, file, dec = ".", sep = "\t") + } + ) + + } else { + print("NO") + tb <- NULL + output$textb <- + renderText({ + "One or more genes is not expressed or does not exist" + }) + } + }) + + ### Percentages ---- + + observeEvent(input$Filter, { + s1 <- unlist(strsplit(as.character(input$String1), split = ",")) + s1 <- gsub(" ", "", as.character(s1)) + s2 <- unlist(strsplit(as.character(input$String2), split = ",")) + s2 <- gsub(" ", "", as.character(s2)) + expressed <- input$Expressed + not_expressed <- input$NonExpressed + + print(expressed) + yes <- + filter(pcttable, id %in% s1, pct.exp >= as.numeric(expressed)) + yes.genes <- names(which(table(yes$gene_name) == length(s1))) + no <- + filter(pcttable, id %in% s2, pct.exp < as.numeric(not_expressed)) + no.genes <- names(which(table(no$gene_name) == length(s2))) + tt3 <- intersect(yes.genes, no.genes) + tt3t <- unique(filter(pcttable, gene_name %in% tt3)[, c(1, 5)]) + tt1 <- filter(yes, gene_name %in% tt3) + tt2 <- filter(no, gene_name %in% tt3) + + + output$YesExpressed <- + DT::renderDataTable({ + DT::datatable( + yes, + options = list(pageLength = 10, autoWidth = TRUE), + rownames = FALSE, + style = 'jQueryUI', + class = 'cell-border stripe' + ) %>% formatStyle(c(1:6), color = "black") + }) + output$NoExpressed <- + DT::renderDataTable({ + DT::datatable( + no, + options = list(pageLength = 10, autoWidth = TRUE), + rownames = FALSE, + style = 'jQueryUI', + class = 'cell-border stripe' + ) %>% formatStyle(c(1:6), color = "black") + }) + output$Result <- + DT::renderDataTable({ + DT::datatable( + as.data.frame(tt3t), + options = list(pageLength = 10, autoWidth = TRUE), + rownames = FALSE, + style = 'jQueryUI', + class = 'cell-border stripe' + ) %>% formatStyle(c(1:2), color = "black") + }) + + output$downloadQuery <- + downloadHandler( + filename = function() { + paste( + "GenesProportion-", + paste(s1, collapse = ","), + "-thrs", + paste(s2, collapse = ","), + ".csv", + sep = "" + ) + }, + content = function(file) { + write.csv(rbind(tt1, tt2), file, sep = "\t") + } + ) + + }) + + +############################################################################# +########### Heatmaps ######################################################## +############################################################################# + + +observeEvent(input$PlotHeatmap, { + + #ss <- unlist(strsplit(as.character(input$genelist), split = ",")) + #ss <- gsub(" ", "", as.character(ss)) + ss <- strsplit(as.character(input$genelist), "\n| |\\,") + ss <- as.data.frame(ss)[,1] + + ss <- unique(c(ss, filter(gene_list, gene_id %in% ss | seqnames %in% ss)$gene_name)) + + missing <- ss[ss %in% c(gene_list$gene_id, gene_list$gene_name, gene_list$seqnames) & !ss %in% med.scaled.long$gene_name & ss %in% gene_list$gene_name] + + flp.neuron.scaled <- med.scaled.long[which(med.scaled.long$gene_name %in% ss),] + flp.ids <- as.character(vlookup(unique(flp.neuron.scaled$gene_name), gene_list, result_column = 1, lookup_column = 2)) + flp.expr <- L4.TPM.medium[flp.ids,, drop=FALSE] + if ( nrow(flp.expr) >1 ) { + flp.neuron.order <- pheatmap(flp.expr, scale = "row") + flp.neuron.order <- flp.neuron.order[["tree_row"]]$order + flp.neuron.order <- flp.ids[flp.neuron.order] + } else { + flp.neuron.order <- flp.ids[1] + } + flp.neuron.order <- as.character(vlookup(flp.neuron.order, gene_list)) + flp.neuron.scaled$gene_name <- factor(flp.neuron.scaled$gene_name, levels = c(rev(flp.neuron.order), missing)) + #flp.neuron.scaled$gene_name <- fct_rev(flp.neuron.scaled$gene_name) + + for( i in missing ){ + dff <- data.frame(gene_name=i, cell.type= colnames(ths)[-c(1,130,131)], scaled.expr=0, prop=0, Modality="NA") + flp.neuron.scaled <- rbind(flp.neuron.scaled, dff) + } + + fnh <- + function() { + withProgress(message = "Generating heatmap Plot...", value = 0, { + ggplot(flp.neuron.scaled, aes(y =gene_name, x = cell.type)) + + geom_point(aes(color = scaled.expr, size = prop)) + + #geom_point(data = NULL, aes(y = vec), pch = NA) + + theme(axis.text.x.bottom = element_text(angle = 90, vjust= 0.5, size = 7, color = "black", hjust = 1), + panel.background = element_blank(), axis.text.y.left = element_text(size = 7, color = "black")) + + scale_color_gradientn("Scaled TPM", colors = c("orange", "maroon", "navy")) + + scale_size_continuous(name = "Proportion", limit = c(0.5, 100), range = c(0,5)) + + labs(y = "Gene", x= "Neuron") + theme(panel.grid = element_line(size = 0.5, color = "grey85")) + }) + } + + + output$downloadheatmap <- + downloadHandler( + filename = function() { + filename = paste("Heatmap.png", sep = "") + }, + content = function(file) { + ggsave( + fnh(), + file = file, + height = 200, + width = 500, + units = "mm" , + limitsize = FALSE, + device = "png" + ) + } + ) + + output$heatmap <- + renderPlot( + ggplot(flp.neuron.scaled, aes(y =gene_name, x = cell.type)) + + geom_point(aes(color = scaled.expr, size = prop)) + + #geom_point(data = NULL, aes(y = vec), pch = NA) + + theme(axis.text.x.bottom = element_text(angle = 90, vjust= 0.5, size = 7, color = "black", hjust = 1), + panel.background = element_blank(), axis.text.y.left = element_text(size = 7, color = "black")) + + scale_color_gradientn("Scaled TPM", colors = c("orange", "maroon", "navy")) + + scale_size_continuous(name = "Proportion", limit = c(0.5, 100), range = c(0,5)) + + labs(y = "Gene", x= "Neuron") + theme(panel.grid = element_line(size = 0.5, color = "grey85")) + ) + + output$dynamic <- renderUI({ + #req(input$plot_hover) + verbatimTextOutput("vals", placeholder = TRUE) + }) + + output$vals <- renderPrint({ + hover <- input$plot_hover + #print(input$plot_hover) # list + y <- nearPoints(flp.neuron.scaled, input$plot_hover) + req(nrow(y) != 0) + y + }) +}) + +### From file + + observeEvent(input$PlotHeatmap2, { + inFile <- input$file1 + ss<-read.table(inFile$datapath, header=FALSE)$V1 + + ss <- strsplit(as.character(ss), "\n| |\\,") + + ss <- unique(c(ss, filter(gene_list, gene_id %in% ss | seqnames %in% ss)$gene_name)) + + missing <- ss[ss %in% c(gene_list$gene_id, gene_list$gene_name, gene_list$seqnames) & !ss %in% med.scaled.long$gene_name & ss %in% gene_list$gene_name] + + flp.neuron.scaled <- med.scaled.long[which(med.scaled.long$gene_name %in% ss),] + flp.ids <- as.character(vlookup(unique(flp.neuron.scaled$gene_name), gene_list, result_column = 1, lookup_column = 2)) + flp.expr <- L4.TPM.medium[flp.ids,, drop=FALSE] + if ( nrow(flp.expr) >1 ) { + flp.neuron.order <- pheatmap(flp.expr, scale = "row") + flp.neuron.order <- flp.neuron.order[["tree_row"]]$order + flp.neuron.order <- flp.ids[flp.neuron.order] + } else { + flp.neuron.order <- flp.ids[1] + } + flp.neuron.order <- as.character(vlookup(flp.neuron.order, gene_list)) + flp.neuron.scaled$gene_name <- factor(flp.neuron.scaled$gene_name, levels = c(rev(flp.neuron.order), missing)) + + for( i in missing ){ + dff <- data.frame(gene_name=i, cell.type= colnames(ths)[-c(1,130,131)], scaled.expr=0, prop=0, Modality="NA") + flp.neuron.scaled <- rbind(flp.neuron.scaled, dff) + } + + fnh <- + function() { + withProgress(message = "Generating heatmap Plot...", value = 0, { + ggplot(flp.neuron.scaled, aes(y =gene_name, x = cell.type)) + + geom_point(aes(color = scaled.expr, size = prop)) + + #geom_point(data = NULL, aes(y = vec), pch = NA) + + theme(axis.text.x.bottom = element_text(angle = 90, vjust= 0.5, size = 7, color = "black", hjust = 1), + panel.background = element_blank(), axis.text.y.left = element_text(size = 7, color = "black")) + + scale_color_gradientn("Scaled TPM", colors = c("orange", "maroon", "navy")) + + scale_size_continuous(name = "Proportion", limit = c(0.5, 100), range = c(0,5)) + + labs(y = "Gene", x= "Neuron") + theme(panel.grid = element_line(size = 0.5, color = "grey85")) + }) + } + + + output$downloadheatmap <- + downloadHandler( + filename = function() { + filename = paste("Heatmap.png", sep = "") + }, + content = function(file) { + ggsave( + fnh(), + file = file, + height = 200, + width = 500, + units = "mm" , + limitsize = FALSE, + device = "png" + ) + } + ) + + output$heatmap <- + renderPlot( + ggplot(flp.neuron.scaled, aes(y =gene_name, x = cell.type)) + + geom_point(aes(color = scaled.expr, size = prop)) + + #geom_point(data = NULL, aes(y = vec), pch = NA) + + theme(axis.text.x.bottom = element_text(angle = 90, vjust= 0.5, size = 7, color = "black", hjust = 1), + panel.background = element_blank(), axis.text.y.left = element_text(size = 7, color = "black")) + + scale_color_gradientn("Scaled TPM", colors = c("orange", "maroon", "navy")) + + scale_size_continuous(name = "Proportion", limit = c(0.5, 100), range = c(0,5)) + + labs(y = "Gene", x= "Neuron") + theme(panel.grid = element_line(size = 0.5, color = "grey85")) + ) + + output$dynamic <- renderUI({ + #req(input$plot_hover) + verbatimTextOutput("vals", placeholder = TRUE) + }) + + output$vals <- renderPrint({ + hover <- input$plot_hover + #print(input$plot_hover) # list + y <- nearPoints(flp.neuron.scaled, input$plot_hover) + req(nrow(y) != 0) + y + }) + }) + +} \ No newline at end of file diff --git a/spinnerManage.js b/spinnerManage.js new file mode 100644 index 0000000..3c7df07 --- /dev/null +++ b/spinnerManage.js @@ -0,0 +1,8 @@ +$(document).ready(function() { + $('#klickButton').click(function() { + $(".loading-spinner").show(); + }); + }); + $(document).on("shiny:connected", function(e) { + $(".loading-spinner").hide(); + }); diff --git a/ui.R b/ui.R new file mode 100644 index 0000000..140d1a8 --- /dev/null +++ b/ui.R @@ -0,0 +1,615 @@ +#### CenGenAPP 2018-2020 +#### Report bugs to: gabrielsantperebaro@gmail.com + +library(BiocManager) +require(shiny) +library(shinyjs) +library(shinythemes) +library(ggridges) +library(DT) +library(cowplot) +require(dplyr) +require(Seurat) +library(ggplot2) +library(shinybusy) +library(tidyverse) +library(pheatmap) +library(expss) +library(plotly) +#library(cairo) +#library(MAST) +options(repos = BiocManager::repositories()) +options("repos") +options(shiny.maxRequestSize = 400 * 1024 ^ 2) +source("Functions.R") + +options(scipen = 0) +options(digits = 2) + +utr <- + c( + "WBGene00023498", + "WBGene00023497", + "WBGene00004397", + "WBGene00006843", + "WBGene00004010", + "WBGene00006789", + "WBGene00001135", + "WBGene00001079", + "WBGene00001135", + "WBGene00006783", + "WBGene00000501", + "WBGene00006788", + "WBGene00001555" + ) + + + +## UI ---- +ui <- fluidPage( + theme = "Theme.min.css", + #shinytheme("flatly"), + #background-color: #2f2d2d; + #shinythemes::themeSelector(), + tags$head(tags$style( + HTML(".shiny-output-error-validation {color: red;}") + )), + + # App title ---- + titlePanel( + "Discovery and analysis of the C. elegans Neuronal Gene Expression Network -- CeNGEN" + ), + hr(), + add_busy_spinner( + spin = "double-bounce", + color = "orange", + timeout = 100, + position = c("top-right"), + onstart = TRUE, + margins = c(10, 10), + height = "50px", + width = "50px" + ), + #add_busy_bar(color = "darkorange"), + # Side panel entering what to plot ---- + + + # Main panel showing plots ---- + tabsetPanel( + type = "pills", + id = "tabs", + ### Cell types Panel ---- + tabPanel( + "Gene expression by cell type", + fluidPage( + hr(), + h6( + "Find which genes are expressed in each cell type ordered by expression level or which cell type express a particular gene." + ), + h6("Select the threshold of expression or unfiltered data"), + h6("4 More stringent, 1 less stringent"), + hr(), + fluidRow( + column(1), + column( + 4, + selectInput( + inputId = "Tcell_name", + label = "Select cell type", + choices = colnames(ths)[2:129], + selected = "ADA" + ), + selectInput( + inputId = "Tcell_cut", + label = "Select threshold", + choices = c(1:4, "Unfiltered"), + selected = 2 + ), + actionButton("TCell", "Expressed genes", icon = icon("hand-o-right")) + + ), + #column(width = 1, offset = 0, style='padding:5px;'), + column( + 4, + textInput( + inputId = "Tgene_name", + label = "Type gene name", + value = "zig-4" + ), + selectInput( + inputId = "Tgene_cut", + label = "Select threshold", + choices = c(1:4, "Unfiltered"), + selected = 2 + ), + actionButton("TGene", "Which cell types", icon = icon("hand-o-right")) + ), + column( + 2, + offset = 0, + style = 'padding:5px;', + textInput( + inputId = "Tgene_name_batch", + label = "Query multiple genes for download", + value = "zig-4,ins-26" + ), + selectInput( + inputId = "Tgene_cut_batch", + label = "Select threshold", + choices = c(1:4, "Unfiltered"), + selected = 2 + ), + downloadButton("TGeneBatch", "Download batch"), + span(textOutput("textb"), style = + "color:red") + + ) + ), + br(), + fluidRow( + column(1), + column( + 4, + br(), + DT::dataTableOutput("Tcell_name_table"), + br(), + uiOutput("get_download_gene") + #downloadButton('downloadGene', "Download table") + ), + #column(width = 1, offset = 0, style='padding:5px;'), + column( + 4, + br(), + DT::dataTableOutput("Tgene_name_table"), + br(), + #downloadButton('downloadCell', "Download table"), + uiOutput("get_download_cell"), + span(textOutput("text1"), style = + "color:red") + ) + ) + ) + ), + + tabPanel( + "Find markers based on percentage of expression", + fluidPage( + hr(), + h6( + "Find which genes are expressed in a group of cell types and not in other cell types." + ), + fluidRow( + #column(1), + column( + 4, + textInput( + inputId = "String1", + label = "Group 1", + value = "AWC_ON,AWC_OFF" + ), + textInput( + inputId = "Expressed", + label = "Minimum percentage of cells expressing the gene", + value = 65 + ), + actionButton("Filter", "Run query", icon = icon("hand-o-right")) + + ), + #column(width = 1, offset = 0, style='padding:5px;'), + column( + 4, + textInput( + inputId = "String2", + label = "Group 2", + value = "SMD,SIB" + ), + textInput( + inputId = "NonExpressed", + label = "Maximum percentage of cells expressing the gene", + value = 2 + ), + downloadButton('downloadQuery', "Download table") + ) + ), + fluidRow( + #column(1), + column(4, + br(), DT::dataTableOutput("YesExpressed")), + column(4, + br(), + DT::dataTableOutput("NoExpressed")), + + column(4, + br(), + DT::dataTableOutput("Result") + #span(textOutput("text3"), style="color:red")) + ) + + ) + ) + ), + ### Enriched types Panel ---- + tabPanel( + "Enriched Genes by Cell Type", + fluidPage( + hr(), + h6( + "Find genes overexpressed in one cell type compared to the rest of cells in the dataset (neurons only or all cell types)." + ), + textOutput("TopMarkers cell plot"), + selectInput( + inputId = "dataset2", + label = "Choose dataset", + choices = c("All cell types", "Neurons only") + ), + conditionalPanel( + "input.dataset2 == 'Neurons only'", + selectInput( + inputId = "Markers", + label = "Select cluster", + choices = sort(unique(allNeurons$`Neuron.type`)), + selected = "ADA" + ), + textInput( + inputId = "top", + label = "Show top X genes", + value = "30" + ), + DT::dataTableOutput("MarkTable"), + downloadButton('downloadMarkers', "Download table"), + h6("HEADER LEGEND:"), + h6( + "p-val and p_val_adj: nominal and adjusted P-values of the test, respectively." + ), + h6( + "pct.1 and pct.2: The percentage of cells where the gene is detected in the first or second group" + ), + h6( + "avg_logFC: Log of the expression fold change between group 1 and group 2." + ) + ), + conditionalPanel( + "input.dataset2 == 'All cell types'", + selectInput( + inputId = "Markers2", + label = "Select cluster", + choices = sort(unique(allCells$`Neuron.type`)), + selected = "Sperm" + ), + textInput( + inputId = "top2", + label = "Show top X genes", + value = "30" + ), + DT::dataTableOutput("MarkTable2"), + downloadButton('downloadMarkers2', "Download table"), + h6("HEADER LEGEND:"), + h6( + "p-val and p_val_adj: nominal and adjusted P-values of the test, respectively." + ), + h6( + "pct.1 and pct.2: The percentage of cells where the gene is detected in the first or second group" + ), + h6( + "avg_logFC: Log of the expression fold change between group 1 and group 2." + ) + ) + ) + ), + ### DEX panel ---- + tabPanel( + "Find Differential Expression between Cell Types", + fluidPage( + hr(), + h6("The calculation can take a few minutes."), + hr(), + fluidRow( + column( + 3, + selectInput( + inputId = "ClusterCells1", + label = "Select cluster 1", + choices = sort(unique(allCells$`Neuron.type`)), + selected = "Sperm" + ), + selectInput( + inputId = "ClusterCells2", + label = "Select cluster 2", + choices = sort(unique(allCells$`Neuron.type`)), + selected = "Intestine" + ), + actionButton("DEXButton", "Calculate DEX", icon = icon("hand-o-right")) + + ), + column( + 3, + selectInput( + inputId = "test", + label = "Select statistical test", + choices = c("wilcox", "bimod", "roc", "t", "LR") + ), + textInput( + inputId = "topM2", + label = "Show top X genes", + value = "30" + ) + + ), + column( + 3, + textInput( + inputId = "batch1", + label = "Group 1: Introduce cell types separated with ,", + value = "AVM,AWA,AWB" + ), + textInput( + inputId = "batch2", + label = "Group 2: Introduce cell types, NEURONS or ALL", + value = "AWC_ON,AWC_OFF" + ), + actionButton( + "DEXButtonBatch", + "Calculate DEX in groups", + icon = icon("hand-o-right") + ) + ) + ), + br(), + h6("HEADER LEGEND:"), + h6( + "p-val and p_val_adj: nominal and adjusted P-values of the test, respectively." + ), + h6( + "pct.1 and pct.2: The percentage of cells where the gene is detected in the first or second group" + ), + h6( + "avg_logFC: Log of the expression fold change between group 1 and group 2." + ), + br(), + + DT::dataTableOutput("MarkTable_ClusterCells"), + span(textOutput("text2"), style = + "color:red"), + br(), + DT::dataTableOutput("MarkTable_Batch"), + downloadButton('downloadDEX', "Download table") + ) + ), + ### Single cell panel ---- + tabPanel("Single cell plot", fluidPage( + hr(), + column( + 3, + h4('Integrated single-cell (10X) analysis'), + wellPanel( + selectInput( + inputId = "dataset", + label = "Choose dataset", + choices = c("All cell types", "Neurons only") + ), + selectInput( + inputId = "Plots", + label = "Features to plot", + choices = c("Metadata", "Individual genes", "Colocalization"), + selected = "Individual genes" + ), + # One feature ---- + conditionalPanel( + ### METADATA PANELS neurons ---- + condition = "input.Plots == 'Metadata' && input.dataset == 'Neurons only'", + selectInput( + "featuresMetadata", + "Metadata", + colnames(allNeurons@meta.data) + ), + + ### CLUSTERS PANELS + conditionalPanel( + condition = "input.featuresMetadata == 'Neuron.type' && input.dataset == 'Neurons only'", + selectInput( + "Cluster", + "Highlight cells", + choices = unique(allNeurons$Neuron.type) + ) + ), + conditionalPanel( + condition = "input.featuresMetadata == 'Experiment' && input.dataset == 'Neurons only'", + selectInput( + "Cluster2", + "Highlight cells", + choices = unique(allNeurons$Experiment) + ) + ), + conditionalPanel( + condition = "input.featuresMetadata == 'Detection' && input.dataset == 'Neurons only'", + selectInput("Cluster3", "Highlight cells", choices = unique(allNeurons$Detection)) + ) + ), + + conditionalPanel( + ### METADATA PANELS all cells ---- + condition = "input.Plots == 'Metadata' && input.dataset == 'All cell types'", + selectInput( + "featuresMetadata2", + "Metadata", + colnames(allCells@meta.data) + ), + ### CLUSTERS PANELS + conditionalPanel( + condition = "input.featuresMetadata2 == 'Neuron.type' && input.dataset == 'All cell types'", + selectInput( + "Cluster.2", + "Highlight cells", + choices = unique(allCells$Neuron.type) + ) + ), + conditionalPanel( + condition = "input.featuresMetadata2 == 'Experiment' && input.dataset == 'All cell types'", + selectInput( + "Cluster2.2", + "Highlight cells", + choices = unique(allCells$Experiment) + ) + ), + conditionalPanel( + condition = "input.featuresMetadata2 == 'Detection' && input.dataset == 'All cell types'", + selectInput("Cluster3.2", "Highlight cells", choices = unique(allCells$Detection)) + ), + conditionalPanel( + condition = "input.featuresMetadata2 == 'Tissue.type' && input.dataset == 'All cell types'", + selectInput("Cluster4", "Highlight cells", choices = unique(allCells$Tissue.type)) + ) + ), + # One gene ---- + conditionalPanel( + condition = "input.Plots == 'Individual genes'", + textInput( + inputId = "GeneName", + label = "Gene name (Symbol)", + value = "ric-4" + ) + ), + # Two gene colocalization ---- + conditionalPanel( + condition = "input.Plots == 'Colocalization'", + textInput( + inputId = "ColocalizationGenes1", + label = "Gene 1", + value = "zig-4" + ), + hr(), + textInput( + inputId = "ColocalizationGenes2", + label = "Gene 2", + value = "F23D12.3" + ) + ), + fluidRow(column( + 6, actionButton("PlotButton", "Plot data", icon = icon("hand-o-right")) + )) + ) + ), + + column( + 9, + textOutput("Single cell plot"), + conditionalPanel( + condition = "input.Plots == 'Metadata'", + downloadLink("downloadClust", "Download Plot"), + plotOutput( + "SinglePlot", + width = "100%", + dblclick = "plotPermanent_dblclick", + brush = brushOpts(id = "plotPermanent_brush", resetOnNew = TRUE) + ), + h5( + "Select region and do double-click to zoom in and double-click to zoom out." + ), + hr(), + downloadLink("downloadHigh", "Download Plot"), + plotOutput("High", width = "100%") + ), + conditionalPanel( + condition = "input.Plots == 'Colocalization'", + sliderInput( + "blend", + "Choose blend", + min = 0, + max = 1, + value = 0.25 + ), + plotOutput("SinglePlotDouble", width = + "50%"), + downloadLink("downloadCol", "Download plot"), + h4(span(textOutput("gene1utr"), style = + "color:red")), + h4(span(textOutput("gene2utr"), style = + "color:red")), + hr(), + h5( + "Select region in the panel below and do double-click to zoom in and double-click to zoom out." + ), + plotOutput( + "SinglePlotPermanent", + width = "50%", + dblclick = "plotPermanent_dblclick", + brush = brushOpts(id = "plotPermanent_brush", resetOnNew = TRUE) + ), + hr() + ), + conditionalPanel( + condition = "input.Plots == 'Individual genes'", + h5( + "Select region and do double-click to zoom in and double-click to zoom out." + ), + plotOutput( + "SinglePlot2", + width = "100%", + dblclick = "plotPermanent_dblclick", + brush = brushOpts(id = "plotPermanent_brush", resetOnNew = TRUE) + ), + h4(span(textOutput("geneUTR"), style = + "color:red")), + hr(), + downloadLink("downloadExp", "Download plot"), + plotOutput("FeaturePlot", width = "100%"), + hr(), + downloadLink("downloadViolin", "Download plot"), + plotOutput("Violin", width = "100%") + ) + ) + )), + tabPanel( + "Heatmaps of gene expression", + fluidPage( + hr(), + h6( + "Introduce a list of genes to display a heatmap showing expression levels and proportion of cells expressing the genes." + ), + textAreaInput( + inputId = "genelist", + label = "Introduce a list of genes", + value = "flp-1\nflp-2,flp-3,WBGene00001447\nWBGene00001448\nflp-6\nflp-7\nflp-8\nflp-9\nflp-10\nflp-11\nflp-12\nflp-13\nflp-14\nflp-15\nflp-16\nflp-17\nflp-18\nflp-19\nflp-20\nflp-21\nflp-22\nflp-23\nflp-24\nflp-25\nflp-26\nflp-27\nflp-28\nflp-32\nflp-33\nflp-34", + width = "500px", + height = "100px" + ), + actionButton( + "PlotHeatmap", + "Plot heatmap from list", + icon = icon("hand-o-right") + ), + hr(), + fluidRow( + column(3,fileInput("file1", NULL, + accept = c( + "text/csv", + "text/comma-separated-values,text/plain", + "txt") + )), + #column(3,actionButton("resetFile", "Clear uploaded file")), + column(3, actionButton( + "PlotHeatmap2", + "Plot heatmap from file", + icon = icon("hand-o-right") + )) + ), + + hr(), + h6( + "You can identify circles by clicking on them." + ), + + + uiOutput("dynamic"), + hr(), + plotOutput("heatmap", width = "100%", hover = "plot_hover"), + + hr(), + downloadLink("downloadheatmap", "Download plot") + + + ) + ) + ) + ) +