added DEqual option to use custom deg_tstats, other plot options

R/DEqual.R

@@ -10,9 +10,12 @@
 #' t-statistics from degradation time adjusting for the six brain regions from
 #' degradation experiment data used for determining `rse_tx`.
 #'
-#' @param DE a `data.frame()` with one column containing the t-statistics from
-#' Differential Expression, typically generated with `limma::topTable()`.
-#' The `rownames(DE)` should be transcript GENCODE IDs.
+#' @param DE a `data.frame()` with a column "t" containing the t-statistics
+#' from Differential Expression, typically generated with `limma::topTable()`.
+#' `rownames(DE)` must have transcript Ensembl/Gencode IDs.
+#' @param deg_tstats an optional`data.frame()` with a column "t" containing
+#' t-statistics resulted from a degradation experiment. Default is the
+#' internal `qsvaR::degradation_tstats` from the package authors.
 #'
 #' @return a `ggplot` object of the DE t-statistic vs
 #' the DE statistic from degradation
@@ -36,48 +39,67 @@
 #'
 #' ## Create the DEqual plot
 #' DEqual(random_de)
-DEqual <- function(DE) {
+DEqual <- function(DE, deg_tstats = qsvaR::degradation_tstats, show.legend = TRUE,
+                        show.cor = c('caption','corner-top','corner-bottom','none')) {
     ## For R CMD check
     DE_t <- degradation_t <- NULL
-
+    show.cor=rlang::arg_match(show.cor)
     ## Check input
-    # stopifnot("t" %in% colnames(DE))
-    # stopifnot(!is.null(rownames(DE)))
-
-    # Check if input is a dataframe
     if (!is.data.frame(DE)) {
       stop("The input to DEqual is not a dataframe.", call. = FALSE)
       }
-
     # Check if 't' is in the column names of DE
     if (!("t" %in% colnames(DE))) {
         stop("'t' is not a column in 'DE'.", call. = FALSE)
       }
-
     # Check if DE has non-null row names
     if (is.null(rownames(DE))) {
         stop("Row names of 'DE' are NULL.", call. = FALSE)
-      }
+    }
+    if (!missing(deg_tstats)) {
+        if (!is.data.frame(deg_tstats)) {
+            stop("The 'deg_tstats' parameter is not a dataframe.", call. = FALSE)
+        }
+        if (!("t" %in% colnames(deg_tstats))) {
+            stop("'t' is not a column in 'deg_tstats'.", call. = FALSE)
+        }
+        if (is.null(rownames(deg_tstats))) {
+            stop("Row names of 'deg_tstats' are NULL.", call. = FALSE)
+        }
+        if (!all(grepl("^ENST\\d+", rownames(deg_tstats)))) {
+            stop("The row names of 'deg_tstats' must be ENSEMBL or Gencode IDs (ENST...)", call. = FALSE)
+        }
+    }
 
     ## Locate common transcripts
-    deg_tstats = qsvaR::degradation_tstats
     whichTx <- which_tx_names(rownames(DE),rownames(deg_tstats))
-    #rownames(deg_tstats) = check_tx_names(rownames(DE),rownames(qsvaR::degradation_tstats),'rownames(DE)','qsvaR::degradation_tstats')
-    #common = intersect(rownames(deg_tstats), rownames(DE))
     common = qsvaR::normalize_tx_names(rownames(DE)[whichTx])
     stopifnot(length(common) > 0)
     rownames(deg_tstats) <- qsvaR::normalize_tx_names(rownames(deg_tstats))
     ## Create dataframe with common transcripts
     common_data <- data.frame(
-        degradation_t = deg_tstats[common, ],
-        DE_t = DE[whichTx, ]
+        degradation_t = deg_tstats[common, "t"],
+        DE_t = DE[whichTx, "t"]
     )
+    cor_val <- signif(cor(common_data[, 1], common_data[, 2]), 2))
     p <- ggplot(common_data, aes(x = DE_t, y = degradation_t)) +
         xlab("DE t-statistic") +
         ylab("Degradation t-statistic") +
-        geom_bin2d(bins = 70) +
+        geom_bin2d(bins = 70, show.legend = FALSE) +
         scale_fill_continuous(type = "viridis") +
-        theme_bw() +
-        labs(caption = paste0("correlation: ", signif(cor(common_data[, 1], common_data[, 2]), 2)))
+        theme_bw()
+        # labs(caption = paste0("correlation: ", cor_val)
+    if (show.cor != 'none') {
+      switch(show.cor,
+        'caption' = {
+          p <- p + labs(caption = paste0("correlation: ", cor_val))
+        },
+        'corner-top' = {
+          p <- p + annotate("text", x = max(DE_t), y = max(degradation_t), label = paste0("correlation: ", cor_val), hjust = 1, vjust = 1)
+        },
+        'corner-bottom' = {
+          p <- p + annotate("text", x = max(DE_t), y = min(degradation_t), label = paste0("correlation: ", cor_val), hjust = 1, vjust = 0)
+        }
+    }
     return(p)
 }