FDR cutoff alpha not applied in rnaseq.Rmd

[Mira0507]

The alpha value is configured in workflows/rnaseq/downstream/config.yaml, as shown below:

  # By default DESeq2 assumes alpha of 0.1 when filtering low-count genes. This
  # is also used to select significant genes in the various helper functions.
  alpha: 0.1

This alpha is intended to be added to metadata(res)$alpha in the make_results function, enabling the extraction of this value in the my_summary function, as shown below:

make_results <- function(dds_name, label, dds_list=NULL, ...){

  if (!is.null(dds_list)){
    dds <- dds_list[[dds_name]]
  } else {
    # Get the dds object from the global dds_list.
    dds <- get_dds(dds_name)
  }

  # Save a copy of the arbitrary arguments
  dots <- list(...)         # <===== arguments added, including alpha!!! 

# Lines skipped...
# Lines skipped...

  # Call results() with the subset of dots that it accepts.
  results_dots <- lcdbwf:::match_from_dots(dots, DESeq2::results)     
  res <- do.call(DESeq2::results, results_dots)        # <===== res built!!! 

  # Adjust log2FoldChange for LRT test
  if (!is.null(dots$test) && dots$test == 'LRT') {
    res$log2FoldChange <- 0           
    warning("All log2FoldChange values in the DESeq2 results object have been set to 0. See https://github.com/lcdb/lcdb-wf/blob/master/docs/rnaseq-rmd.rst?plain=1#L269.")
  }

  # Checks for LRT test and non-NULL type
  if (!is.null(dots$type) && !is.null(dots$test) && dots$test == 'LRT' && !test_detected) {
    stop("You cannot pass a non-NULL or missing type to make_results with test == 'LRT'. For LRT, LFC values are set to 0 and should not be passed to lfcShrink. Use type == NULL in make_results for LRT DDS objects.")
  } else if (!is.null(dots$type) && !is.null(dots$test) && dots$test == 'LRT' && test_detected) {
    stop("You cannot pass a non-NULL or missing type to make_results with an LRT dds object. For LRT, LFC values are set to 0 and should not be passed to lfcShrink. Use type == NULL in make_results for LRT DDS objects.")
  }

  # Call lfcShrink if applicable
  if (!is.null(dots$type) && dots$test != 'LRT') {
    dots[['res']] <- res
    dots[['dds']] <- dds

    lfcShrink_dots <- lcdbwf:::match_from_dots(dots, DESeq2::lfcShrink)
    res <- do.call(DESeq2::lfcShrink, lfcShrink_dots)

    S4Vectors::metadata(res)$type <- dots$type          # <======= metadata of `res` updated!!!
  } else {
    S4Vectors::metadata(res)$type <- NULL
  }

  return(
    list(
      res=res,
      dds=dds_name,
      label=label
    )
  )

my_summary <- function(res, dds, name, ...){
  dds_label <- dds
  if (class(dds) != 'character'){
    stop("expecting dds to be a string")
  }
  dds <- lcdbwf:::get_dds(dds_label)
  alpha <- metadata(res)$alpha         # <===== alpha extracted from `res` obj

# Lines skipped... 

Turned out that this is not a bug, but it would be helpful to explicitly add the alpha argument in the rnaseq.Rmd, as shown below:

# results_01 chunk:
contr_01_main <- lcdbwf:::make_results(
  dds_name='main',
  label='LFC>0',
  contrast=treatment - control,
  type='ashr',
  alpha=config$main$alpha         # <==== updated!!!
)