R/gs_mds.R

In GeneTonic: Enjoy Analyzing And Integrating The Results From Differential Expression Analysis And Functional Enrichment Analysis

#' Multi Dimensional Scaling plot for gene sets
#'
#' Multi Dimensional Scaling plot for gene sets, extracted from a `res_enrich`
#' object
#'
#' @param res_enrich A `data.frame` object, storing the result of the functional
#' enrichment analysis. See more in the main function, [GeneTonic()], to check the
#' formatting requirements (a minimal set of columns should be present).
#' @param res_de A `DESeqResults` object.
#' @param annotation_obj A `data.frame` object with the feature annotation
#' information, with at least two columns, `gene_id` and `gene_name`.
#' @param n_gs Integer value, corresponding to the maximal number of gene sets to
#' be included (from the top ranked ones). Defaults to the number of rows of
#' `res_enrich`
#' @param gs_ids Character vector, containing a subset of `gs_id` as they are
#' available in `res_enrich`. Lists the gene sets to be included, additionally to
#' the ones specified via `n_gs`. Defaults to NULL.
#' @param similarity_measure Character, currently defaults to `kappa_matrix`, to
#' specify how to compute the similarity measure between gene sets
#' @param mds_k Integer value, number of dimensions to compute in the multi
#' dimensional scaling procedure
#' @param mds_labels Integer, defines the number of labels to be plotted on top
#' of the scatter plot for the provided gene sets.
#' @param mds_colorby Character specifying the column of `res_enrich` to be used
#' for coloring the plotted gene sets. Defaults sensibly to `z_score`.
#' @param gs_labels Character vector, containing a subset of `gs_id` as they are
#' available in `res_enrich`. Lists the gene sets to be labeled.
#' @param plot_title Character string, used as title for the plot. If left `NULL`,
#' it defaults to a general description of the plot and of the DE contrast
#' @param return_data Logical, whether the function should just return the
#' data.frame of the MDS coordinates, related to the original `res_enrich`
#' object. Defaults to FALSE.
#'
#' @return A `ggplot` object
#'
#' @seealso [create_kappa_matrix()] is used to calculate the similarity between
#' gene sets
#'
#' @export
#'
#' @examples
#'
#' library("macrophage")
#' library("DESeq2")
#' library("org.Hs.eg.db")
#' library("AnnotationDbi")
#'
#' # dds object
#' data("gse", package = "macrophage")
#' dds_macrophage <- DESeqDataSet(gse, design = ~line + condition)
#' rownames(dds_macrophage) <- substr(rownames(dds_macrophage), 1, 15)
#' dds_macrophage <- estimateSizeFactors(dds_macrophage)
#'
#' # annotation object
#' anno_df <- data.frame(
#'   gene_id = rownames(dds_macrophage),
#'   gene_name = mapIds(org.Hs.eg.db,
#'                      keys = rownames(dds_macrophage),
#'                      column = "SYMBOL",
#'                      keytype = "ENSEMBL"),
#'   stringsAsFactors = FALSE,
#'   row.names = rownames(dds_macrophage)
#' )
#'
#' # res object
#' data(res_de_macrophage, package = "GeneTonic")
#' res_de <- res_macrophage_IFNg_vs_naive
#'
#' # res_enrich object
#' data(res_enrich_macrophage, package = "GeneTonic")
#' res_enrich <- shake_topGOtableResult(topgoDE_macrophage_IFNg_vs_naive)
#' res_enrich <- get_aggrscores(res_enrich, res_de, anno_df)
#'
#' gs_mds(res_enrich,
#'        res_de,
#'        anno_df,
#'        n_gs = 200,
#'        mds_labels = 10)
gs_mds <- function(res_enrich,
                   res_de,
                   annotation_obj,
                   n_gs = nrow(res_enrich),
                   gs_ids = NULL,
                   similarity_measure = "kappa_matrix",
                   mds_k = 2,
                   mds_labels = 0,
                   mds_colorby = "z_score",
                   gs_labels = NULL,
                   plot_title = NULL,
                   return_data = FALSE) { # or aggr_score

  similarity_measure <- match.arg(similarity_measure,
                                  c("kappa_matrix", "overlap_matrix"))

  # require res_enrich to have aggregated scores and so
  if (!("z_score" %in% colnames(res_enrich))) {
    res_enrich <- get_aggrscores(res_enrich,
                                 res_de,
                                 annotation_obj = annotation_obj)
  }

  n_gs <- min(n_gs, nrow(res_enrich))

  gs_to_use <- unique(
    c(
      res_enrich$gs_id[seq_len(n_gs)],  # the ones from the top
      gs_ids[gs_ids %in% res_enrich$gs_id]  # the ones specified from the custom list
    )
  )


  if (similarity_measure == "kappa_matrix") {
    my_simmat <- create_kappa_matrix(res_enrich, n_gs = n_gs, gs_ids = gs_ids)

  } else if (similarity_measure == "overlap_matrix") {
    my_simmat <- create_jaccard_matrix(res_enrich, n_gs = n_gs, gs_ids = gs_ids, return_sym = TRUE)
  } # else ...

  # subset here, internally
  res_enrich <- res_enrich[gs_to_use, ]
  mds_labels <- min(mds_labels, nrow(res_enrich))

  mysets <- res_enrich[["gs_id"]]
  mysets_names <- res_enrich[["gs_description"]]


  mds_gs <- cmdscale((1 - my_simmat), eig = TRUE, k = mds_k)

  mds_gs_df <- data.frame(
    dim1 = mds_gs$points[, 1],
    dim2 = mds_gs$points[, 2],
    gs_id = mysets,
    gs_name = mysets_names,
    gs_DEcount = res_enrich$DE_count,
    gs_colby = res_enrich[[mds_colorby]],
    gs_text = paste0(mysets, ": ", mysets_names),
    stringsAsFactors = FALSE
  )

  if (return_data) {
    return(mds_gs_df)
  }

  max_z <- max(abs(range(mds_gs_df$gs_colby)))
  limit <- max_z * c(-1, 1)

  this_contrast <- (sub(".*p-value: (.*)", "\\1", mcols(res_de, use.names = TRUE)["pvalue", "description"]))

  p <- ggplot(mds_gs_df, aes_string(x = "dim1",
                                    y = "dim2",
                                    text = "gs_text")) +
    geom_point(aes_string(color = "gs_colby",
                          size = "gs_DEcount")) +
    scale_color_gradient2(limit = limit,
                          low = muted("deepskyblue"),
                          mid = "lightyellow",
                          high = muted("firebrick")) +
    labs(x = "Dimension 1",
         y = "Dimension 2",
         col = mds_colorby,
         size = "Geneset \nmembers"
    ) +
    theme_bw()

  if (mds_labels > 0) {
    label_these <- mds_gs_df$gs_id[seq_len(mds_labels)]
  } else {
    label_these <- NULL
  }

  if (!is.null(gs_labels)) {
    if (!all(gs_labels %in% res_enrich$gs_id)) {
      warning("Not all specified geneset ids were found in the `res_enrich` object")
    }
    label_those <- mds_gs_df[mds_gs_df$gs_id %in% gs_labels, "gs_id"]
  } else {
    label_those <- NULL
  }

  df_gs_labels <- mds_gs_df[mds_gs_df$gs_id %in% unique(c(label_these, label_those)), ]

  p <- p + geom_label_repel(
    aes_string(label = "gs_name"), data = df_gs_labels, size = 3)

  # handling the title
  if (is.null(plot_title)) {
    p <- p + ggtitle(paste0("Geneset MDS plot - ", this_contrast))
  } else {
    p <- p + ggtitle(plot_title)
  }

  return(p)
}