R/admin.R

In crickbabs/RNASeq-DESeq: DESeq2 Wrapper For Pre-Specified Analyses

#' R6 class representing the set of parameters used in an analysis
#' 
#' @export
ParamList <- R6::R6Class("ParamList",
                    private = list(
                      params=list(),
                      descriptions=list(alpha="The threshold of statistical significance is {deparse(alpha)}",
                                        lfcThreshold="The threshold for absolute effect size is {}",
                                        title="The name of the project is '{}'",
                                        top_n_variable="Only use {} genes for unsupervised clustering",
                                        script="The project was generated by '{}'",
                                        showCategory= "Only show top {} enriched categories in plots",
                                        seed="A random seed of {} is used to ensure reproduciblity",
                                        filterFun="{if (is.null(filterFun)) 'Default' else deparse(filterFun)} independent filtering",
                                        baseMeanMin="Discard transcripts with few average counts per sample than {}"
                                        ),
                      defaults=list()
                    ),
                    public = list(
                      #' @description
                      #' Create a new set of parameters.
                      #' @param defaults Named list of default values.  Names are the parameters, and the values will be their default.
                      #' @return An object that will store all future values of analysis parameters.
                      initialize = function(defaults=list()) {
                        private$defaults <- defaults
                      },
                      #' @description
                      #' Set the value of a parameter
                      #' @param id The name of the parameter to be set.
                      #' @param value The value the parameter should taken henceforth; if missing, it will take the default value.
                      #' @param description A string describing what the purpose of the parameter is.
                      #' @param div Logical, whether to mention in the markdown report what the value has been set to.
                      set = function(id, value, description="", div=TRUE) {
                        if (missing(value)) {
                          if (id %in% names(private$defaults)) {
                            value <- private$defaults[[id]]
                          } else {
                            eg <- readLines(system.file("templates/example.spec", package="DESdemonA"))
                            eg <- eg[grep(paste0("^\\s+", id),eg)]
                            if (length(eg)>0) {
                              stop("Attempt to set '", id, "' but no value or default provided.\n", "You may need to update your spec file, as new settings have been introduced.  Maybe lines like:\n", paste(eg, collapse="\n"))
                            } else {
                              stop("Attempt to set '", id, "' but no value or default provided.\n")
                            }
                          }
                        }
                        if (is.null(value)) {
                          private$params[id] <- list(NULL)
                        } else {
                          private$params[[id]]=value
                        }
                        if (description=="") {
                          if (!(id %in% names(private$descriptions))) { # provide tautological definition if 
                            description=paste0("The value of ", id, " is {", id, "}")
                          } else {
                            description=private$descriptions[id]}
                        } 
                        description <- sub("\\{\\}", paste0("{", id, "}"), description)
                        private$descriptions[[id]]=description
                        if (div & isTRUE(getOption('knitr.in.progress'))) {
#                          cat(knitr::knit_child(text=knitr::knit_expand(text=c("```{block, type='rparam'}", self$describe(id), "```")), quiet=TRUE))
                          cat("\n\n<div class=\"rparam\">", self$describe(id), "</div>\n\n")
                        }
                        invisible(self$get(id))
                      },
                      #' @description
                      #' Get the value that the parameter is currently set to.
                      #' @param id Name of the value you want to access.
                      get = function(id) {
                        if (!id %in% names(private$params)) {
                          stop(id, " has not yet been initialized")
                        }
                        ret <- private$params[[id]]
                        if (is.call(ret)) eval(ret) else ret
                      },
                      #' @description
                      #' Turn the mutable object into a list
                      publish = function() {
                        lapply(private$params, eval)
                      },
                      #' @description
                      #' Get a text description of what the setting is, and what value it currently takes.
                      #' @param id Name of the value you want to access.
                      describe = function(id) {
                        if (missing(id)) {
                          map(private$descriptions[names(private$params)],  function(d) glue::glue_data(.x =private$params, d))
                        } else {
                          glue::glue_data(private$params,
                                          private$descriptions[[id]])
                        }
                      }
                    )
                    )


##' Populate a folder with DESdemonA starter scripts
##'
##' To start a DESdemonA-based project, we provide a sample set
##' of scripts to point you in the right direction. Firstly there is
##' an '00_init.r' file that will create a universal DESeqDataSet object
##' from your quantified counts file.  You may need to edit this to link
##' it to where your quantification pipeline stores its results, and to
##' ensure that the full set of metadata is inserted into the colData.
##'
##' There is an example '.spec' file - rename and use this as a basis
##' for your statistical analysis plan, or if you have an existing one,
##' delete the example one and copy the existing one into the folder instead.
##' 
##' The main analysis is run via "01_analyse.r" - you should render this
##' via rmarkdown.  It will look for every '.spec' file in the current
##' directory.
##'
##' There will also be a "02_further_steps.r" script at some point. This
##' will give concrete examples of how you might want to extract results
##' for further programmatic use, to build upon the html report that
##' rendering the "01_analyse.r" will provide.
##'
##' There's also a DESCRIPTION file, so that it is easy to turn your
##' analysis into a re-distributable R package.
##'
##' The recommended usage is, at the start of project development,
##' to simply call 'DESdemonA::get_started()' in the relevant directory,
##' as the defaults path and files are sufficient - it will refuse
##' to overwrite existing files, so is safe in that sense.
##' 
##' @title Initiate a DESdemonA project
##' @param files Which files to retreive from the DESdemonA project
##' @param path Where to copy the files to
##' @return 
##' @author Gavin Kelly
#' 
#' @export

get_started <- function(files = dir(system.file("templates",package="DESdemonA")),
                path=".",
                yml="",
                overwrite=FALSE,
                file_col="filename",
                name_col="name",
                ...
                ) {
  args <- list(...)
  if (yml!="") {
    yml_args <- read_yml(yml)
    ind <- setdiff(yml_args, args)
    args[ind] <- yml_args[ind]
  }
  defaults <- list(
    nfcore="results",
    metadata=system.file("extdata/metadata.xlsx", package="babsrnaseq"),
    file_col=deparse(substitute(file_col)),
    name_col=deparse(substitute(name_col)),
    counts=quote(file.path(nfcore, "star_rsem")),
    org_package="",
    project=basename(getwd()),
    author=getOption("usethis.full_name")
  )
  ind <- setdiff(names(defaults), names(args))
  args <- c(args, defaults[ind])
  args <- lapply(args, eval, args)
  pre_exist <- file.exists(file.path(path, files))
  if (any(pre_exist) && (!overwrite)) {
    stop(paste(file.path(path, files), collapse=", "), " already exist. Remove or rename them")
  }
  for (fname in files) {
    if (file.exists(fname) && overwrite) {
      unlink(fname)
    }
    usethis::use_template(fname, save_as=fname, data=args, package="DESdemonA")
  }
}


##' Run DESdemonA Report on existing counts object
##'
##' This will generate a standard report on the DESeq2 data object you
##' provide it. It will store data objects in the `data` directory, so
##' that will need to be created, as will the results folder.
##'
##' 
##' @param dds The DESeqDataSet object that you want to run the report
##'   on. It needs the basic set of `colData` columns that are used in
##'   the analysis plan. `colnames(dds)` will be used as labels in
##'   plot, etc. In addition, if its `mcols` has columns of the that
##'   are set to `entrez` and/or `symbol`, these will get added to the
##'   report.  `metadata(dds)$org <- "org.Mm.eg.db"` will ensure that
##'   the correct annotation libraries are used, but this is only strictly
##'   necessary for functional annotation (which also requires those additional
##'   columns in the `mcols` property.
##' 
##' @param spec_file The Analaysis Plan
##' @param results Directory in which to store excel results
##' @param output_file The name of the html report.
##' @param title HTML Title of the document
##' @param autor The name of the author to appear on the report
##' @return 
##' @author Gavin Kelly
#' 
#' @export

run_report <- function(dds, spec_file, results="results", output_file="analyse.html", title="RNASeq Analysis", author=Sys.info["user"]) {
  count_source=deparse1(substitute(dds))
  repeat{
    fname <- paste0(tempfile("analyse", tmpdir="."), ".r")
    if (!file.exists(fname)) break
  }
  file.copy(system.file("templates/01_analyse.r", package="DESdemonA"), fname)
  rmarkdown::render(fname,
                    output_file=output_file,
                    params=list(res_dir=results,
                                spec_file=spec_file,
                                count_source=dds,
                                param_call=list(count_source=count_source)
                                )
                    )
  unlink(fname)
}


read_yml <- function(file) {
  lines <- readLines(file)
  lines <- lines[grepl("  .*:", lines)]
  field <- gsub("  (*[^:]*):.*", "\\1", lines)
  value <- gsub("  *[^:]*: *", "", lines)
  setNames(as.list(value), field)
}