R/visualization_functions.R

In cgplab/RadAR: Radiomic Analysis with R

#' Draw a clustered heatmap of the radiomic dataset
#'
#' This function uses \code{\link{pheatmap}} to draw a heatmap of feature values (rows) across samples (columns). The rows and columns of the
#' heatmap are ordered based on previous required hierarchical clustering generated by \code{\link{do_hierarchical_clustering}}.
#' Dendrograms for columns and rows are also reported with optional annotation tracks.
#'
#' @param rdr A RadAR object (class \code{\link{SummarizedExperiment}}).
#' @param show_feature_names (logical) Should feature names be displayed in the heatmap?
#' @param show_sample_names (logical) Should sample names be displayed in the heatmap?
#' @param annotation_tracks (character) Tracks to report in the annotation section of the heatmap.
#' Available annotation tracks can be obtained by \code{\link{names(colData(rdr))}}.
#' @param annotation_data (data.frame) A data frame specifying tracks to report in the annotation section. For
#' more details refer to \code{\link{?pheatmap::pheatmap}}.
#' @param annotation_colors (list) List for specifying annotation_col track colors manually.For
#' more details and examples refer to \code{\link{?pheatmap::pheatmap}}.
#' @param ncolors (integer) Number of colors of the color palette.
#' @param save_plot_to (character) If specified, filename where the heatmap will be saved to.
#' @param width_plot (numeric) output file width in inches.
#' @param height_plot (numeric) output file height in inches.
#'
#' @return none
#' @author Matteo Benelli (\email{matteo.benelli@uslcentro.toscana.it})
#' @export
#'
#' @examples
plot_heatmap_hcl <- function(rdr = NULL,
                             show_feature_names = F,
                             show_sample_names = F,
                             annotation_tracks = NULL,
                             annotation_data = NULL,
                             annotation_colors = NULL,
                             ncolors = 50,
                             save_plot_to = NULL,
                             width_plot = 20,
                             height_plot = 15
)
{

  assertthat::assert_that(length(rdr) > 0, msg = "[RadAR] Error: rdr object required")
  assertthat::assert_that(all(annotation_tracks %in% colnames(colData(rdr))), msg = "[RadAR] Error: annotation_tracks contains invalid track names")
#  assertthat::assert_that(which_data %in% c("original", "scaled", "norm"), msg = "[RadAR] Error: Invalid data type")

  if (length(annotation_data) > 0) {
    assertthat::assert_that(assertthat::are_equal(rownames(annotation_data), colnames(rdr)),
                            msg = "[RadAR] Error: Invalid rownames for annotation_data")
  }

  palette_heatmap <- rev(colorspace::sequential_hcl(n = ncolors,
                                                    palette = "Blues3",
                                                    power = 2))
  which_data <- metadata(rdr)$which_data_hcl
  hcl <- metadata(rdr)$hcl

  assertthat::assert_that(length(hcl) > 0,
                          msg = "[RadAR] Error: hierarchical clustering has not been yet computed")

  if (which_data == "normal") {
    data <- assays(rdr)$values
  }
  if (which_data == "scaled") {
    data <- assays(rdr)$scaled_values
    assertthat::assert_that(!is.null(data), msg = "[RadAR] Error: Feature values have not been yet scaled")
  }
  if (which_data == "normalized") {
    data <- assays(rdr)$norm_values
    assertthat::assert_that(!is.null(data), msg = "[RadAR] Error: Feature values have not been yet normalized")
  }


  if (length(annotation_data) == 0 & !is.null(annotation_tracks)) {
    annotation_data <- data.frame(colData(rdr)[, annotation_tracks])
    rownames(annotation_data) <- colnames(rdr)
  }
  colnames(data) <- colnames(rdr)
  pheatmap::pheatmap(data, cluster_rows = hcl$hcl_row, cluster_cols = hcl$hcl_col,
                     color = palette_heatmap,
                     show_rownames = show_feature_names,
                     show_colnames = show_sample_names,
                     annotation_col = annotation_data,
                     annotation_colors = annotation_colors,
                     border_color = NA,
                     scale = ifelse(which_data == "normal", "row", "none"),
                     filename = ifelse(length(save_plot_to)>0, save_plot_to, NA),
                     width = ifelse(length(save_plot_to)>0, width_plot, NA),
                     height = ifelse(length(save_plot_to)>0, height_plot, NA)
  )
}

#' Draw correlation matrix of radiomic features
#'
#' This function draws a correlation matrix of the radiomic dataset by \code{\link{corrplot}} or \code{\link{pheatmap}}.
#' Correlation matrix can be visualized by clustered correlation plot or heatmap.
#'
#' @param rdr A RadAR object (class \code{\link{SummarizedExperiment}}).
#' @param view_as (character) How to visualize correlation matrix? It can be "corrplot" or "heatmap".
#' @param method_correlation (character) Which method use to calculate correlation. It can be one of the following:
#' "pearson", "kendall", "spearman".
#' @param method_hcl (character) Which method use for clustering matrix columns and rows. Available methods are
#' "ward.D", "ward.D2", "single", "complete","average" , "mcquitty","median", "centroid".
#' @param which_data (character) Which data use for the computation of the correlation coefficients. It can be one of the following: "normal",
#' "scaled", "normalized".
#' @param cex_labels (numeric) Text label cex.
#' @param srt_labels (numeric) Text label string rotation in degrees.
#' @param ncolors (numeric) Number of colors to generate the palette.
#' @param save_plot_to (character) If specified, filename where the heatmap will be saved to.
#' @param width_plot (numeric) output file width in inches.
#' @param height_plot (numeric) output file height in inches.
#' @param show_names (logical) Should feature names (labels) be visualized?
#' @param show_dendrogram (logical) Should dendogram (for columns and rows) be drawn (valid for heatmap visualzation)?
#'
#' @return none
#' @author Matteo Benelli (\email{matteo.benelli@uslcentro.toscana.it})
#' @export
#'
#' @examples
plot_correlation_matrix <- function(rdr = NULL,
                                    view_as=c("heatmap", "corrplot"),
                                    method_correlation = "pearson",
                                    method_hcl = "ward.D",
                                    which_data = "normal",
                                    cex_labels = .5,
                                    srt_labels = 45,
                                    ncolors = 50,
                                    save_plot_to=NULL,
                                    width_plot = 20,
                                    height_plot = 15,
                                    show_names = T,
                                    show_dendrogram = T

)
{

  if (assertthat::are_equal( which_data, c("normal", "scaled", "normalized"))) {
    assertthat::assert_that(1<0, msg = "[RadAR] which_data needs to be specified")
  }

  if (length(view_as) > 1) {
    assertthat::assert_that(1<0, msg = "[RadAR] view_as needs to be specified")
  }
  methods_hcl <- c("ward.D",
                   "ward.D2",
                   "single",
                   "complete",
                   "average" ,
                   "mcquitty",
                   "median",
                   "centroid")
  assertthat::assert_that(length(rdr) > 0, msg = "[RadAR] Error: rdr object required")
  assertthat::assert_that(view_as %in% c("heatmap", "corrplot"), msg = "[RadAR] Error: Invalid visualization setting")
  assertthat::assert_that(method_correlation %in% c("pearson", "kendall", "spearman"), msg = "[RadAR] Error: Invalid correlation method setting")
  assertthat::assert_that(which_data %in% c("normal", "scaled", "normalized"), msg = "[RadAR] Error: Invalid data type")
  assertthat::assert_that(cex_labels > 0, msg = "[RadAR] IError: nvalid cex value for labels")
  assertthat::assert_that(srt_labels > 0, msg = "[RadAR] Error: Invalid srt value for labels")

  palette_cor <- colorspace::diverge_hcl(n = ncolors,
                                         palette = "Blue-Red3",
                                         power = 2)

  if (which_data == "normal") {
    data <- assays(rdr)$values
  }
  if (which_data == "scaled") {
    data <- assays(rdr)$scaled_values
    assertthat::assert_that(!is.null(data), msg = "[RadAR] Error: Feature values have not been yet scaled")
  }
  if (which_data == "normalized") {
    assertthat::assert_that(!is.null(data), msg = "[RadAR] Error: Feature values have not been yet normalized")
  }
  res <- cor(t(data), method = method_correlation, use = "na")
  if (view_as == "corrplot") {
    if (length(save_plot_to) > 0) {
      pdf (file = save_plot_to, width = width_plot, height = height_plot)
    }
    corrplot::corrplot(res, type = "upper", order = "hclust",
                       tl.col = "black",
                       tl.srt = srt_labels,
                       tl.cex = cex_labels,
                       tl.pos = ifelse(show_names, "lt", "n"),
                       col = palette_cor,
                       method = "color",
                       hclust.method = method_hcl
                       )
    if (length(save_plot_to) > 0) {
      dev.off()
    }

  }
  if (view_as == "heatmap") {
    pheatmap::pheatmap(res, symm = T,
                       col = palette_cor,
                       border_color=NA,
                       symbreaks = T,
                       clustering_distance_rows = "euclidean",
                       clustering_distance_cols = "euclidean",
                       clustering_method = method_hcl,
                       fontsize_row = cex_labels*12,
                       fontsize_col = cex_labels*12,
                       filename = ifelse(length(save_plot_to)>0, save_plot_to, NA),
                       width = ifelse(length(save_plot_to)>0, width_plot, NA),
                       height = ifelse(length(save_plot_to)>0, height_plot, NA),
                       show_rownames = show_names,
                       show_colnames = show_names,
                       treeheight_row = ifelse(show_dendrogram, 50, 0),
                       treeheight_col = ifelse(show_dendrogram, 50, 0)

    )
  }
}



# draw boxplot for single features  -----------------------------------------------------------
#' Draw boxplot + stripchart of selected feature(s)
#'
#' This function generates boxplots + stripcharts of selected feature(s). Data are stratified by
#' user defined conditions.
#'
#' @param rdr A RadAR object (class \code{\link{SummarizedExperiment}}).
#' @param feature_names (character) Which feature(s) should be plotted. Available feature names are
#' obtained by \code{rowData(rdr)$feature_name}.
#' @param conditions (numeric, character) Vector of labels specifying conditions for each sample in rdr.
#' Should have the same length of \code{ncol(rdr)}.
#' @param which_data (character) Which data use for plot. It can be one of the following: "normal", "scaled", "normalized".
#' @param max_range (logical) Should plot(s) be ranged between min and max values?
#' @param user_range (numeric) A vector of length 2 reporting min and max ylim values.
#'
#' @return none
#' @author Matteo Benelli (\email{matteo.benelli@uslcentro.toscana.it})
#' @export
#'
#' @examples
plot_features <- function(rdr = NULL,
                          feature_names= NULL,
                          conditions = NULL,
                          which_data = "scaled",
                          max_range = F,
                          user_range = NULL
)

{
  assertthat::assert_that(length(rdr) > 0, msg = "[RadAR] Error: rdr object required")
  assertthat::assert_that(length(feature_names) > 0 & length(feature_names) <= 12,
                          msg = "[RadAR] Error: feature_names should be specified. Maximum of 12 features admitted")
  assertthat::assert_that(length(conditions) == ncol(rdr), msg = "[RadAR] Error: conditions need to be specified")
  assertthat::assert_that(which_data %in% c("normal", "scaled", "normalized"), msg = "[RadAR] Error: Invalid data type")
  assertthat::assert_that(all(feature_names %in% rownames(rdr)), msg = "[RadAR] Error: Invalid feature names")
  n_features_to_plot <- length(feature_names)
  if (length(user_range) > 0) {
    assertthat::assert_that(length(user_range) ==2, msg = "[RadAR] Error: user_range should be a vector of length 2")
  }

  if (which_data == "normal") {
    data <- assays(rdr)$values
  }
  if (which_data == "scaled") {
    data <- assays(rdr)$scaled_values
    assertthat::assert_that(!is.null(data), msg = "[RadAR] Error: Feature values have not been yet scaled")
  }
  if (which_data == "normalized") {
    data <- assays(rdr)$norm_values
    assertthat::assert_that(!is.null(data), msg = "[RadAR] Error: Feature values have not been yet normalized")
  }

  par (mfrow = n2mfrow(n_features_to_plot))
  for (i in 1: length(feature_names)) {
    range_feature <- range(data[feature_names[i], ], na.rm = T)
    if (max_range & which_data != "normal") {
      range_feature <- c(0, 1)
    }
    if (length(user_range) == 2) {
      range_feature <- as.numeric(user_range)
    }

    boxplot(data[feature_names[i], ] ~ conditions, varwidth = T, frame.plot = F, las = 1,
            pch = 19, outline = F, ylim = range_feature, main = feature_names[i],
            xlab = "", ylab = "Feature value")
    stripchart(data[feature_names[i], ] ~ conditions, pch = 19, col = "grey60", add = T, vertical = T,
               method = "jitter")
    grid()
  }
  par (mfrow = c(1, 1))

}