R/graphics.R

In specmine: Metabolomics and Spectral Data Analysis and Mining

### BOXPLOT of each variable ####

"boxplot_variables" = function(dataset, variables = NULL, samples = NULL, horizontal = TRUE,
                               col = "lightblue", nchar.label = 10, cex.axis = 0.8, ...) {
  
  if (is.null(variables)) {
    variables = rownames(dataset$data)
  }
  
  if (is.null(samples)) {
    samples = colnames(dataset$data)
  }
  if (is.numeric(variables))
    names.short = substr(get_x_values_as_text(dataset)[variables], 1, nchar.label)
  else
    names.short = substr(variables, 1, nchar.label)
  
  if (length(variables) > 1) {
    boxplot(t(dataset$data[variables, samples]), names = names.short,
            horizontal = horizontal, las = 2,
            col = col, cex.axis = cex.axis, ...)
  }
  else {
    boxplot(dataset$data[variables, samples], xlab = names.short,
            horizontal = FALSE, las = 2,
            col = col, cex.axis = cex.axis, ...)
  }
}

boxplot_vars_factor = function(dataset, meta.var, variables = NULL, samples = NULL,
                               horizontal = FALSE, nchar.label = 10, col = NULL,
                               vec.par = NULL, cex.axis = 0.8, ylabs = NULL, ...)
{
  if (is.null(variables)) {
    variables = rownames(dataset$data)
  }
  
  if (is.null(samples)) {
    samples = colnames(dataset$data)
  }
  
  if (is.numeric(variables))
    names.short = substr(get_x_values_as_text(dataset)[variables], 1, nchar.label)
  else
    names.short = substr(variables, 1, nchar.label)
  
  if (is.null(vec.par))
    vec.par = c(length(variables), 1)
  
  opar = par(mfrow = vec.par)
  on.exit(par(opar))
  for (i in 1:length(variables)) {
    if (is.null(col)) coli = i + 1
    else coli = col
    if (!is.null(ylabs)) {
      ylab = ylabs[i]
    } else {
      ylab = NULL
    }
    boxplot(dataset$data[variables[i], samples] ~ dataset$metadata[, meta.var],
            horizontal = horizontal, las = 2, main = names.short[i],
            col = coli, cex.axis = cex.axis, ylab = ylab, ...)
  }
}

plotvar_twofactor = function(dataset, variable, meta.var1, meta.var2, colour = "darkblue", title = "",
                             xlabel = NULL, ylabel = NULL)
{
  
  df = data.frame(dataset$data[variable, ], dataset$metadata[, meta.var1],
                  dataset$metadata[, meta.var2])
  if (is.numeric(variable)) n1 = rownames(dataset$data)[variable]
  else n1 = variable
  if (is.numeric(meta.var1)) n2 = colnames(dataset$metadata)[meta.var1]
  else n2 = meta.var1
  if (is.numeric(meta.var2)) n3 = colnames(dataset$metadata)[meta.var2]
  else n3 = meta.var2
  colnames(df) = c("name1", "name2", "name3")
  
  g = ggplot2::ggplot(data = df)
  g = g + ggplot2::aes_string("name2", "name1")
  g = g + ggplot2::geom_boxplot(fill = colour)
  g = g + ggplot2::facet_grid(. ~ name3)
  if (is.null(xlabel)) g = g + ggplot2::xlab(n2)
  else g = g + ggplot2::xlab(xlabel)
  if (is.null(ylabel)) g = g + ggplot2::ylab(n1)
  else g = g + ggplot2::ylab(ylabel)
  if (title != "") g = g + ggplot2::ggtitle(title)
  g
}

##############################SPECTRA PLOTS###############################

"plot_spectra_simple" = function(dataset, samples = NULL, variable.bounds = NULL, xlab = NULL,
                                 ylab = NULL, lty = 1, lwd = 1, col = 1, reverse.x = FALSE, ...) {
  
  if (is.null(xlab)) xlab = get_x_label(dataset)
  if (is.null(ylab)) ylab = get_value_label(dataset)
  
  if (!is.null(dataset$labels$x) && !is.expression(dataset$labels$x) && dataset$labels$x == "mz/rt") {
    if (is.null(variable.bounds)) {
      variables = 1:length(get_x_values_as_text(dataset))
      vars = as.numeric(gsub("/.*", "", get_x_values_as_text(dataset)))
    } else {
      x.vars = as.numeric(gsub("/.*", "", get_x_values_as_text(dataset)))
      variables = which(x.vars > variable.bounds[1] & x.vars < variable.bounds[2])
      vars = x.vars[variables]
    }
  } else {
    if (is.null(variable.bounds)) {
      variables = rownames(dataset$data)
      vars = variables
    }
    else {
      x.vars = get_x_values_as_num(dataset)
      variables = rownames(dataset$data)[x.vars > variable.bounds[1] & x.vars < variable.bounds[2]]
      vars = variables
    }
  }
  if (reverse.x) xlim = c(max(as.numeric(vars)), min(as.numeric(vars)))
  else xlim = range(as.numeric(vars))
  
  if (is.null(samples)) {
    samples = colnames(dataset$data)
  }
  
  matplot(vars, dataset$data[variables, samples, drop = FALSE], type = "l", lty = lty, col = col,
          xlab = xlab, ylab = ylab, xlim = xlim, ...)
}

"plot_spectra" = function(dataset, column.class, func = NULL, samples = NULL,
                          variable.bounds = NULL, xlab = NULL, ylab = NULL, lty = 1,
                          legend.place = "topright", cex = 0.8, reverse.x = FALSE, ...) {
  
  if (is.null(xlab)) xlab = get_x_label(dataset)
  if (is.null(ylab)) ylab = get_value_label(dataset)
  
  if (!is.null(dataset$labels$x) && !is.expression(dataset$labels$x) && dataset$labels$x == "mz/rt") {
    if (is.null(variable.bounds)) {
      variables = 1:length(get_x_values_as_text(dataset))
      vars = as.numeric(gsub("/.*", "", get_x_values_as_text(dataset)))
    } else {
      x.vars = as.numeric(gsub("/.*", "", get_x_values_as_text(dataset)))
      variables = which(x.vars > variable.bounds[1] & x.vars < variable.bounds[2])
      vars = x.vars[variables]
    }
  } else {
    if (is.null(variable.bounds)) {
      variables = rownames(dataset$data)
      vars = variables
    }
    else {
      x.vars = get_x_values_as_num(dataset)
      variables = rownames(dataset$data)[x.vars > variable.bounds[1] & x.vars < variable.bounds[2]]
      vars = variables
    }
  }
  
  if (reverse.x) xlim = c(max(as.numeric(vars)), min(as.numeric(vars)))
  else xlim = range(as.numeric(vars))
  
  if (is.null(samples)) {
    samples = colnames(dataset$data)
    metadata = dataset$metadata[, column.class]
  }
  else {
    metadata = factor(dataset$metadata[samples, column.class])
  }
  if (is.null(func)) {
    matplot(vars, dataset$data[variables, samples], type = "l", lty = lty, col = as.integer(metadata),
            xlab = xlab, ylab = ylab, xlim = xlim, ...)
    if (legend.place != "none")
      legend(legend.place, levels(metadata), cex = cex, fill = sort(as.integer(factor(levels(metadata)))))
  }
  else {
    aggregate.result = aggregate(t(dataset$data[variables, samples]), by = list(metadata), func)
    matplot(vars, t(aggregate.result[-1]), type = "l", lty = 1, col = as.integer(aggregate.result[, 1]),
            xlab = xlab, ylab = ylab, xlim = xlim, ...)
    if (legend.place != "none")
      legend(legend.place, levels(metadata), cex = cex, fill = sort(as.integer(factor(levels(metadata)))))
  }
}

#################################
############2D PLOT##############
#################################

snr_spectra <- function(spec){
  threshold <- mean(base::Filter(isPositive, spec)) / (sd(base::Filter(isPositive, spec[which(apply(spec, 1, cv) > 15), which(apply(spec, 2, cv) > 15)])))
  threshold
}

snr_all <- function(spect, metadata = NULL) {
  tresholds <- unlist(lapply(spect$data, snr_spectra))
  samples <- names(spect$data)
  if (!is.null(metadata)) {
    meta_data <- as.factor(spect$metadata[samples, metadata])
    res <- data.frame(SNR = tresholds, Meta = meta_data, row.names = samples)
    res <- res[order(-res$SNR, res$Meta), ]
  } else {
    res <- data.frame(SNR = tresholds, Samples = samples)
    res <- res[order(-res$SNR), ]
  }
  res
}

plot_2d_spectra <- function(specmine_2d_dataset, title_spectra = "", meta = NULL, spec_samples = NULL) {
  clrs <- c("Greys", "YlGnBu", "Greens", "YlOrRd", "Bluered", "RdBu", "Reds", "Blues", "Picnic", "Rainbow", "Portland", "Jet", "Hot", "Blackbody", "Earth", "Electric", "Viridis", "Cividis")
  
  if (!is.null(spec_samples)) {
    if (typeof(spec_samples) == "character") {
      samples <- names(specmine_2d_dataset$data)[which(names(specmine_2d_dataset$data) == spec_samples)]
    } else if (typeof(spec_samples) == "double") {
      samples <- names(specmine_2d_dataset$data)[spec_samples]
    }
  } else {
    if (!is.null(meta)) {
      snr_spectra <- snr_all(spect = specmine_2d_dataset, metadata = meta)
      samples <- rownames(snr_spectra[stats::ave(-snr_spectra$SNR, snr_spectra$Meta, FUN = rank) <= 2, ])
    } else {
      snr_spectra <- snr_all(spect = specmine_2d_dataset)
      samples <- rownames(snr_spectra)[c(1:2, nrow(snr_spectra) - 1, nrow(snr_spectra))]
    }
  }
  
  butons <- list()
  for (i in (1:(length(samples) + 1))) {
    temp <- list()
    temp$method <- "restyle"
    if (i == length(samples) + 1) {
      temp$args <- list("visible", c(rep(TRUE, length.out = length(samples))))
      temp$label <- "All samples"
    } else if (i != length(samples) + 1) {
      temp$args <- list("visible", c(rep(FALSE, length.out = length(samples))))
      temp$args[[2]][i] <- TRUE
      temp$label <- samples[i]
    }
    butons[[i]] <- temp
  }
  
  p1 <- plotly::plot_ly(showscale = FALSE, hovertemplate = paste(
    "F1: %{x:.0f}<br>",
    "F2: %{y:.0f}<br>",
    "Intensity: %{z:.0f}"
  ))
  
  state <- specmine_2d_dataset$metadata[samples[1], meta]
  index <- 2
  for (sample in samples) {
    if (!is.null(meta)) {
      new_state <- specmine_2d_dataset$metadata[sample, meta]
      if (new_state != state) {
        state <- new_state
        index <- index + 1
      }
      colorstate <- clrs[index]
      p1 <- plotly::add_surface(p1, z = specmine_2d_dataset$data[[sample]], name = paste(state, sample, sep = ": "), visible = FALSE, colorscale = colorstate, showlegend = TRUE, legendgroup = state)
    } else {
      colorstate <- clrs[index]
      p1 <- plotly::add_surface(p1, z = specmine_2d_dataset$data[[sample]], name = sample, visible = FALSE, colorscale = colorstate, showlegend = TRUE)
      index <- index + 1
    }
  }
  
  initial_val <- as.numeric(rownames(specmine_2d_dataset$data[[samples[1]]]))[1]
  final_val <- as.numeric(rownames(specmine_2d_dataset$data[[samples[1]]]))[nrow(specmine_2d_dataset$data[[samples[1]]])]
  if (initial_val < 0) {
    initial_val <- 0
  }
  if (final_val < 0) {
    final_val <- 0
  }
  
  initial_val2 <- as.numeric(colnames(specmine_2d_dataset$data[[samples[1]]]))[1]
  if (initial_val2 < 0) {
    initial_val2 <- 0
  }
  final_val2 <- as.numeric(colnames(specmine_2d_dataset$data[[samples[1]]]))[ncol(specmine_2d_dataset$data[[samples[1]]])]
  if (final_val2 < 0) {
    final_val2 <- 0
  }
  
  seq1 <- seq(initial_val, final_val, length.out = 5)
  seq2 <- seq(initial_val2, final_val2, length.out = 5)
  
  seq_dim1 <- seq(0, nrow(specmine_2d_dataset$data[[samples[1]]]), length.out = 5)
  seq_dim2 <- seq(0, ncol(specmine_2d_dataset$data[[samples[1]]]), length.out = 5)
  
  xaxis_1 <- list(title = "F2", tickvals = seq_dim2, ticktext = seq2)
  yaxis_2 <- list(title = "F1", tickvals = seq_dim1, ticktext = seq1)
  
  p1 <- plotly::layout(
    p1,
    title = title_spectra,
    scene = list(xaxis = xaxis_1, yaxis = yaxis_2, zaxis = list(title = "Intensity")),
    updatemenus = list(list(type = "dropdown", y = 0.8, buttons = butons)),
    legend = list(title = list(text = paste0("<b>", meta)))
  )
  p1
}

##################################
#####PLOT OF MS AND NMR PEAKS#####
##################################

plot_peaks = function(dataset, column.class, samples = NULL,
                      variable.bounds = NULL, xlab = NULL, ylab = NULL,
                      legend.place = "topright", cex = 0.8, reverse.x = FALSE, p.size = 0.5, ...) {
  
  if (is.null(xlab)) xlab = get_x_label(dataset)
  if (is.null(ylab)) ylab = get_value_label(dataset)
  
  if (is.null(variable.bounds)) {
    variables = rownames(dataset$data)
    vars = variables
  }
  else {
    x.vars = get_x_values_as_num(dataset)
    variables = rownames(dataset$data)[x.vars > variable.bounds[1] & x.vars < variable.bounds[2]]
    vars = variables
  }
  
  if (reverse.x) xlim = c(max(as.numeric(vars)), min(as.numeric(vars)))
  else xlim = range(as.numeric(vars))
  
  if (is.null(samples)) {
    samples = colnames(dataset$data)
    metadata = dataset$metadata[, column.class]
  }
  else {
    metadata = factor(dataset$metadata[samples, column.class])
  }
  
  matplot(vars, dataset$data[variables, samples], col = as.integer(metadata),
          xlab = xlab, ylab = ylab, xlim = xlim, type = "p", pch = 16, cex = p.size, ...)
  
  if (legend.place != "none")
    legend(legend.place, levels(metadata), cex = cex, fill = sort(as.integer(factor(levels(metadata)))))
}

##########################################################################################################
## Multiplot from ggplot2

multiplot <- function(plots, plotlist = NULL, file, cols = 1, layout = NULL) {
  plots <- c(plots, plotlist)
  numPlots = length(plots)
  
  if (is.null(layout)) {
    layout <- matrix(seq(1, cols * ceiling(numPlots / cols)),
                     ncol = cols, nrow = ceiling(numPlots / cols))
  }
  
  if (numPlots == 1) {
    plot(plots[[1]])
  } else {
    grid::grid.newpage()
    grid::pushViewport(grid::viewport(layout = grid::grid.layout(nrow(layout), ncol(layout))))
    
    for (i in 1:numPlots) {
      matchidx <- as.data.frame(which(layout == i, arr.ind = TRUE))
      grob_i <- ggplot2::ggplotGrob(plots[[i]])
      grid::grid.draw(
        grob_i,
        vp = grid::viewport(
          layout.pos.row = matchidx$row,
          layout.pos.col = matchidx$col
        )
      )
    }
  }
  
  invisible(plots)
}