R/chemoDivPlot.R

In chemodiv: Analysing Chemodiversity of Phytochemical Data

#' Plot chemodiversity
#'
#' Function to conveniently create basic plots of the different types of
#' chemodiversity measurements calculated by functions in the package.
#' This function exists to provide an easy way to make basic chemodiversity
#' plots. As functions in the package output data in standard formats,
#' customized plots are easily created with \code{\link[ggplot2]{ggplot2}}.
#'
#' The function can create four different types of plots,
#' (using \code{\link[ggplot2]{ggplot2}}) depending on which input data
#' is supplied:
#' \itemize{
#' \item Function argument \code{compDisMat}. A compound dissimilarity matrix
#' will be plotted as a dendrogram visualizing how
#' structurally/biosynthetically similar different compounds are to each other.
#' \item Function argument \code{divData}. Diversity/evenness values will
#' be plotted as a boxplot.
#' \item Function argument \code{divProfData}. A diversity profile, plotting
#' (Functional) Hill diversity at different values of q will be plotted
#' as a line plot.
#' \item Function argument \code{sampDisMat}. A sample dissimilarity matrix
#' will be plotted as an NMDS plot.
#' \item Function argument \code{groupData}. Grouping data (e.g. population,
#' species etc.) may be supplied, to plot each group in different
#' boxes/lines/colours.
#' }
#' Note that this function can take any combination of the four arguments
#' as input, and argument names should always be specified to ensure
#' each dataset is correctly plotted. If including the function
#' argument \code{sampDisMat}, a Nonmetric Multidimensional Scaling (NMDS)
#' will be performed, which may take time for larger datasets.
#'
#' @param compDisMat Compound dissimilarity matrix, generated by
#' the \code{\link{compDis}} function. Note that only a single matrix
#' should be supplied, and not the whole list.
#' @param divData Diversity/evenness data frame,
#' generated by the \code{\link{calcDiv}} function. This data frame can
#' contain a single or multiple columns with diversity/evenness measures.
#' @param divProfData Diversity profile, generated by
#' the \code{\link{calcDivProf}} function. Note that the whole list
#' outputted by the \code{\link{calcDivProf}} function should be supplied.
#' @param sampDisMat Sample dissimilarity matrix, generated by
#' the \code{\link{sampDis}} function. This can be either the list of one or
#' both matrices outputted by the function, or a single matrix directly.
#' @param groupData Grouping data. Should be either a vector or a data frame
#' with a single column.
#'
#' @return The specified chemodiversity plots.
#'
#' @export
#'
#' @examples
#' minimalDiv <- calcDiv(minimalSampData, minimalCompDis, type = "FuncHillDiv")
#' groups <- c("A", "A", "B", "B")
#' chemoDivPlot(divData = minimalDiv, groupData = groups)
#'
#' data(alpinaCompDis)
#' data(alpinaSampDis)
#' data(alpinaPopData)
#' alpinaDiv <- calcDiv(sampleData = alpinaSampData, compDisMat = alpinaCompDis,
#' type = "FuncHillDiv")
#' alpinaDivProf <- calcDivProf(sampleData = alpinaSampData,
#' compDisMat = alpinaCompDis, type = "FuncHillDiv",
#' qMin = 0, qMax = 2, step = 0.2)
#' chemoDivPlot(compDisMat = alpinaCompDis, divData = alpinaDiv,
#' divProfData = alpinaDivProf, sampDisMat = alpinaSampDis,
#' groupData = alpinaPopData)
chemoDivPlot <- function(compDisMat = NULL,
                         divData = NULL,
                         divProfData = NULL,
                         sampDisMat = NULL,
                         groupData = NULL) {

  allPlots <- list()

  if (is.data.frame(groupData)) {
    groupData <- as.vector(groupData[,1])
  }

  if (!is.null(compDisMat)) { # Compound dendrogram

    compDisMatClust <- stats::hclust(stats::as.dist(compDisMat),
                                     method = "average")
    compDisMatClustDend <- stats::as.dendrogram(compDisMatClust)
    compDisMatClustDendData <- ggdendro::dendro_data(compDisMatClustDend)

    compDisMatTreePlot <- ggplot2::ggplot() +
      ggplot2::geom_segment(data = compDisMatClustDendData$segments,
                            ggplot2::aes(x = .data$x,
                                         y = .data$y,
                                         xend = .data$xend,
                                         yend = .data$yend)) +
      ggplot2::geom_text(data = compDisMatClustDendData$labels,
                         ggplot2::aes(x = .data$x,
                                      y = .data$y,
                                      label = .data$label),
                         hjust = -0.1, angle = 0) +
      ggplot2::scale_y_reverse(limits = c(1, -0.5),
                               breaks = c(1, 0.75, 0.5, 0.25, 0)) +
      ggplot2::ylab("Dissimilarity") +
      ggplot2::ggtitle("") +
      ggplot2::theme(axis.title.y = ggplot2::element_blank(),
                     axis.text.y = ggplot2::element_blank(),
                     axis.ticks.y = ggplot2::element_blank(),
                     panel.grid.major = ggplot2::element_blank(),
                     panel.grid.minor = ggplot2::element_blank(),
                     panel.border = ggplot2::element_blank(),
                     panel.background = ggplot2::element_blank(),
                     axis.line.x = ggplot2::element_line(color = "black",
                                                         size = 0.5)) +
      ggplot2::coord_flip()

    allPlots[["compDisMatTreePlot"]] <- compDisMatTreePlot
  }

  if (!is.null(divData)) { # Boxplot of diversity/evenness

    divDatadf <- as.data.frame(divData)

    if (is.null(groupData)) {
      message("No grouping data provided.")
      groupData <- rep("NoGroup", nrow(divDatadf))
    }
    divDatadf$Group <- groupData

    for (i in 1:(ncol(divDatadf)-1)) {
      allPlots[[paste0("divPlot", colnames(divDatadf)[i])]] <- local({
        i <- i
        currentCol <- colnames(divDatadf)[i]

        divPlot <- ggplot2::ggplot(data = divDatadf,
                                   ggplot2::aes(x = .data$Group,
                                                y = .data[[currentCol]],
                                                fill = .data$Group)) +
          ggplot2::geom_boxplot(outlier.shape = NA) +
          ggplot2::geom_jitter(height = 0, width = 0.1, shape = 21) +
          ggplot2::ylab(currentCol) +
          ggplot2::theme(text = ggplot2::element_text(size = 15),
                         legend.position = "none")
      })
    }
  }

  if (!is.null(divProfData)) { # Diversity profile

    if (is.null(groupData)) {
      message("No grouping data provided.")
      groupData <- rep("NoGroup", nrow(divProfData$divProf))
    }
    divProf <- divProfData$divProf

    # Get mean data in order
    divProfMean1 <- stats::aggregate(divProf,
                                     by = list(Group = groupData), mean,
                                     na.rm = TRUE)
    divProfMean2 <- as.data.frame(t(divProfMean1[, 2:ncol(divProfMean1)]))
    colnames(divProfMean2) <- divProfMean1$Group
    qAll <- seq(from = divProfData$qMin,
                to = divProfData$qMax,
                by = divProfData$step)
    divProfMean2$q <- qAll
    divHillLong <- tidyr::pivot_longer(divProfMean2,
                                       1:(ncol(divProfMean2) - 1),
                                       names_to = "Group",
                                       values_to = "Diversity")

    # Get individual sample data in order
    divProfInd <- as.data.frame(t(divProf))
    divProfInd$q <- qAll
    divHillLongInd <- tidyr::pivot_longer(divProfInd,
                                          1:(ncol(divProfInd) - 1),
                                          names_to = "Individual",
                                          values_to = "Diversity")
    divHillLongInd$Group <- rep(groupData, length(unique(divProfInd$q)))

    divProfPlot <- ggplot2::ggplot() +
      ggplot2::geom_line(data = divHillLongInd,
                         ggplot2::aes(x = .data$q,
                                      y = .data$Diversity,
                                      group = .data$Individual,
                                      color = .data$Group),
                         linewidth = 0.5, alpha = 0.15) +
      ggplot2::geom_line(data = divHillLong,
                         ggplot2::aes(x = .data$q,
                                      y = .data$Diversity,
                                      color = .data$Group),
                         linewidth = 2) +
      ggplot2::xlab("Diversity order (q)") +
      ggplot2::ylab(divProfData$type) +
      ggplot2::theme(text = ggplot2::element_text(size = 15))

    allPlots[["divProfPlot"]] <- divProfPlot
  }

  if (!is.null(sampDisMat)) { # NMDS

    if (is.null(groupData)) {
      message("No grouping data provided.")

      if (is.matrix(sampDisMat)) { # If input is a single matrix
        groupData <- rep("NoGroup", nrow(sampDisMat))

      } else { # If input is a list with multiple matrices
        groupData <- rep("NoGroup", nrow(sampDisMat[[1]]))
      }
    }

    # If input is a single matrix
    if (is.matrix(sampDisMat)) {

      # Suppressing iteration output for vegan::metaMDS, this requires
      # capture.output() rather than the more common suppressMessage()
      utils::capture.output(NMDS <- vegan::metaMDS(sampDisMat,
                                                   autotransform = FALSE))
      NMDSCoords <- as.data.frame(NMDS$points)
      NMDSCoords$Group <- groupData
      NMDSPlot <- ggplot2::ggplot(data = NMDSCoords,
                                  ggplot2::aes(x = .data$MDS1,
                                               y = .data$MDS2,
                                               color = .data$Group)) +
        ggplot2::geom_point(size = 4, alpha = 0.5) +
        ggplot2::theme(text = ggplot2::element_text(size = 15))
      allPlots[["NMDSPlot"]] <- NMDSPlot

    } else { # If input is a list with a BrayCurtis and a GenUniFrac matrix

      utils::capture.output(BCNMDS <- vegan::metaMDS(sampDisMat$BrayCurtis,
                                                     autotransform = FALSE))
      BCNMDSCoords <- as.data.frame(BCNMDS$points)
      BCNMDSCoords$Group <- groupData
      BCNMDSPlot <- ggplot2::ggplot(data = BCNMDSCoords,
                                    ggplot2::aes(x = .data$MDS1,
                                                 y = .data$MDS2,
                                                 color = .data$Group)) +
        ggplot2::geom_point(size = 4, alpha = 0.5) +
        ggplot2::theme(text = ggplot2::element_text(size = 15)) +
        ggplot2::ggtitle("Bray-Curtis NMDS")
      allPlots[["BCNMDSPlot"]] <- BCNMDSPlot


      utils::capture.output(GUNMDS <- vegan::metaMDS(sampDisMat$GenUniFrac,
                                                     autotransform = FALSE))
      GUNMDSCoords <- as.data.frame(GUNMDS$points)
      GUNMDSCoords$Group <- groupData
      GUNMDSPlot <- ggplot2::ggplot(data = GUNMDSCoords,
                                    ggplot2::aes(x = .data$MDS1,
                                                 y = .data$MDS2,
                                                 color = .data$Group)) +
        ggplot2::geom_point(size = 4, alpha = 0.5) +
        ggplot2::theme(text = ggplot2::element_text(size = 15)) +
        ggplot2::ggtitle("Generalized UniFrac NMDS")
      allPlots[["GUNMDSPlot"]] <- GUNMDSPlot
    }
  }
  return(gridExtra::grid.arrange(grobs = allPlots,
                                 ncol = ceiling(sqrt(length(allPlots)))))
}