update 1.7

ChangeLog

@@ -0,0 +1,7 @@
+changes in version 1.7 (2022-11-22)
+	* add results_apeglm to apply apeglm shrinkage on all contrasts
+	* add print_genes to print scaled values from plot_genes
+	* fix order of contrasts in results_all and check_contrasts
+
+changes in version 1.6 (2021-09-01)
+	* adding leadingN to fgsea_all output with number of leading edge genes

DESCRIPTION

@@ -1,6 +1,6 @@
 Package: hciR
 Title: RNA-seq workflows at HCI
-Version: 1.6
+Version: 1.7
 Author: Chris Stubben
 Maintainer: Chris Stubben <[email protected]>
 Imports:
@@ -24,7 +24,8 @@ Imports:
     ggplot2,
     fgsea
 Suggests:
-    clusterProfiler
+    clusterProfiler,
+    dendextend
 License: GPL-3
 Encoding: UTF-8
 Description: This package simplifies the R code in a differential expression analysis for markdown reports. The package includes functions to run DESeq2 using sample and count tibbles as input, get annotated DESeq results for all pairwise comparisons and create interactive plots and other visualizations.

NAMESPACE

@@ -29,6 +29,7 @@ export(plot_interactions)
 export(plot_ma)
 export(plot_pca)
 export(plot_volcano)
+export(print_genes)
 export(r_log)
 export(read_RSEM)
 export(read_RnaSeqMetrics)
@@ -42,6 +43,7 @@ export(read_ipa)
 export(read_sample_files)
 export(results2)
 export(results_all)
+export(results_apeglm)
 export(sort_counts)
 export(sort_samples)
 export(summary_deseq)

R/check_contrasts.R

@@ -34,8 +34,7 @@ check_contrasts <- function(trt, vs="all"){
 	   # DESeq2 will create factor for character strings in alphabetical order
       n <- sort(unique(trt))
    }
-   n <- rev(n)
-   contrast <- utils::combn(n, 2)
+   contrast <- utils::combn(n, 2)[2:1,, drop=FALSE]
    if( vs == "combined"){
       ## if two columns are combined into a single trt group, compare within first group
       n1 <- apply(contrast, 2, function(x) length(unique( gsub("[ _-].+", "", x)))==1)

R/print_genes.R

@@ -0,0 +1,33 @@
+#' Print the scaled values from plot_genes and optionally add branches from cutree
+#'
+#' @param x a tibble from \code{\link{top_counts}}
+#' @param cut number of branches to cut
+#'
+#' @return tibble with id, branch and rlog values used by \code{\link{plot_genes}}
+#'
+#' @author Chris Stubben
+#'
+#' @examples
+#' x <- top_counts(pasilla$results, pasilla$rlog)
+#' plot_genes(x, c("condition", "type"), scale="row", annotation_names_col=FALSE)
+#' print_genes(x, cut=5)
+#' @export
+
+print_genes <- function(x, cut){
+  # scale by rows
+  y <- scale(t(as_matrix(x)))
+  x1 <- as_tbl(t(y))
+  attr(x1, "colData") <- attr(x, "colData")
+  p <- plot_genes(x1, scale="none", silent=TRUE)
+  ## reorder rows and columns in matrix to match heatmap
+  y2 <- y[p$tree_col$order, p$tree_row$order]
+  y2 <- round(y2,3)
+  ## convert to tibble (easier to print and save)
+  y2 <- as_tbl(t(y2))
+  if(!missing(cut)){
+     ## cut column dendrogram
+     n <- dendextend::cutree(p$tree_row, cut, order=FALSE)
+     y2 <- tibble::add_column(y2, branch=n, .after="id")
+  }
+  y2
+}

R/results_all.R

@@ -55,8 +55,7 @@ results_all <- function( object, biomart,  vs="all", subset, relevel, alpha = 0.
     }else{
       n <- relevel
    }
-   n <- rev(n)
-   contrast <- utils::combn(n, 2)
+   contrast <- utils::combn(n, 2)[2:1,, drop=FALSE]
    if( vs == "combined"){
       ## if two columns are combined into a single trt group, compare within first group
       n1 <- apply(contrast, 2, function(x) length(unique( gsub("[ _-].+", "", x)))==1)

R/results_apeglm.R

@@ -0,0 +1,95 @@
+#' Run all results from a DESeq analysis using apeglm shrinkage
+#'
+#' @param object a DESeqDataSet
+#' @param biomart annotations from \code{read_biomart} with column 1 matching row names in results
+#' @param alpha the significance cutoff for the adjusted p-value cutoff (FDR)
+#' @param add_columns a vector of biomart columns to add to result table, default
+#'        gene_name, biotype, chromosome, description and human_homolog if present
+#' @param trt Compare groups within trt group, default is first term in the design formula
+#' @param simplify return a tibble if only 1 contrast present
+#' @param \dots additional options passed to \code{results}
+#'
+#' @return A list of tibbles for each contrast
+#'
+#' @author Chris Stubben
+#'
+#' @examples
+#' \dontrun{
+#' library(hciRdata)
+#' res <- results_apeglm(pasilla$dds, fly98)
+#' res
+#' }
+#' @export
+
+results_apeglm <- function( object, biomart, alpha = 0.05, add_columns, trt, simplify=TRUE,  ...){
+   message("Using adjusted p-value < ", alpha)
+
+   if(missing(trt)){
+      n <- as.character( DESeq2::design(object))
+      ## [1] "~"   "trt"
+      # drop intercept if "0 + trt"?
+      n[2] <- gsub("^0 \\+ ", "", n[2])
+      ## multiple terms in design formula  ~ trt + mouse
+      if(grepl(" + ", n[2], fixed=TRUE)){
+         n[2] <- gsub(" \\+.*", "", n[2])
+      }
+      trt <- n[2]
+   }
+   ## sample groups
+   grps <- levels(object[[trt]])
+   n <- length(grps)
+   contrast <- utils::combn(grps, 2)[2:1,, drop=FALSE]
+  # contrast names
+   vs <- apply(contrast, 2, paste, collapse = " vs ")
+   if(length(vs) == 0) stop("No contrasts found")
+   ## padded for message
+   vs1 <- sprintf(paste0("%-", max(nchar(vs))+2, "s"), paste0(vs, ":") )
+
+   ## resultNames for contrasts
+   resNames <- paste0( trt, "_", gsub(" ", "_", vs))
+   ## number contrasts for each reference level 1 1 1 2 2 3
+   n2 <- rep(1:(n-1), (n-1):1)
+
+   res <- vector("list", length(vs))
+   names(res) <- vs
+   if(missing(add_columns)){
+        add_columns <- c("gene_name", "biotype", "chromosome", "description")
+        if(!missing(biomart)){
+           if("human_homolog" %in% names(biomart)) add_columns <- c(add_columns, "human_homolog")
+        }
+   }
+   message("Adding apeglm's shrunken estimates to log2FoldChange, saving unshrunken in MLE_log2FC")
+   for(i in 1:(n-1)){
+      if(i > 1){
+       #  message("Setting ", grps[i], " as reference level")
+        ## create new levels...
+        object[[trt]] <- stats::relevel(object[[trt]], ref = grps[i])
+        object <- suppressMessages( DESeq2::nbinomWaldTest(object))
+      }
+      for(j in which(n2 %in% i)){
+        res1 <- DESeq2::results(object, name=resNames[j] , alpha = alpha, ...)
+          mle_log2FC <- res1$log2FoldChange
+          res1 <-  DESeq2::lfcShrink(object, resNames[j], res=res1, type = "apeglm", quiet=TRUE)
+          res1$MLE_log2FC <- mle_log2FC
+
+          ft <- S4Vectors::metadata(res1)$filterThreshold
+          x <- suppressMessages( summary_deseq(res1) )
+          message(j, ". ", vs1[j], x[1,2], " up and ", x[2,2], " down regulated" )
+          if(!missing(biomart)){
+              # suppress messages like 70 rows in results are missing from biomart table and print once
+              res1 <- suppressMessages( annotate_results( res1, biomart, add_columns) )
+          }else{
+             ## gene_name or id?
+             res1 <- tibble::as_tibble(tibble::rownames_to_column(data.frame(res1), var="id"))
+          }
+          attr(res1, "contrast") <- vs[j]
+          attr(res1, "alpha") <- alpha
+          attr(res1, "filterThreshold") <- ft
+          res[[j]] <- res1
+      }
+   }
+   n1 <-is.na(res[[1]][[3]])
+   if(any(n1)) message("Note: ", sum(n1), " rows in results are missing from biomart table")
+   if(simplify & length(res) == 1)  res <- res[[1]]
+   res
+}

R/summary_tibble.R

@@ -23,8 +23,8 @@ summary_tibble <- function(x){
       unique = sapply(x, function(y) length( stats::na.omit( unique(y)))),
          NAs = sapply(x, function(y) sum(is.na(y) ) ),
         ### suppressWarnings to avoid :  min(c(NA, NA), na.rm=TRUE)
-        min  = suppressWarnings( apply(x, 2, min, na.rm=TRUE )),
-        max  = suppressWarnings( apply(x, 2, max, na.rm=TRUE )),
+        min  = trimws(suppressWarnings( apply(x, 2, min, na.rm=TRUE ))),
+        max  = trimws(suppressWarnings( apply(x, 2, max, na.rm=TRUE ))),
         ## will be slow with many rows...
        top3 = sapply(x, top3)
    )

R/write_gsea_rnk.R

@@ -66,6 +66,8 @@ write_gsea_rnk <- function(res, write=TRUE, protein_coding = TRUE, na_pvalue = F
           x <- dplyr::filter(y, gene_name != "") %>%
                 dplyr::select(gene_name, log2FoldChange)
       }
+      ## remove NA fold change.
+      x <- dplyr::filter(x, !is.na(log2FoldChange))
       ## remove duplicates
       x <- dplyr::arrange(x, gene_name, dplyr::desc( abs(log2FoldChange)))
       n <- duplicated(x$gene_name)