R/delta_process.R

In biospi/deamass: Statistical modelling for mass spectrometry omics

#' @import data.table
#' @include generics.R
#' @include seaMass_delta.R
setMethod("process", "seaMass_delta", function(object, chain, job.id) {
  ctrl <- control(object)
  if (ctrl@version != as.character(packageVersion("seaMass")))
    stop(paste0("version mismatch - '", name(object), "' was prepared with seaMass v", ctrl@version, " but is running on v", packageVersion("seaMass")))

  cat(paste0("[", Sys.time(), "]  DELTA-PROCESS name=", name(object)," chain=", chain, "\n"))

  ellipsis <- ctrl@ellipsis
  ellipsis$object <- object
  ellipsis$chains <- chain

  # group dea
  if (ctrl@model != "" && !all(is.na(assay_design(object, as.data.table = T)$Condition))) {
    do.call(paste("dea", ctrl@model, sep = "_"), ellipsis)
  }

  if (increment_completed(filepath(object), "process", job.id) == ctrl@nchain) {
    cat(paste0("[", Sys.time(), "]  DELTA-OUTPUT name=", name(object), "\n"))

    # summarise group de and perform fdr correction
    if (file.exists(file.path(filepath(object), "group.quants.de.index.fst"))) {
      if(ctrl@fdr.model != "") {
        ellipsis <- ctrl@ellipsis
        ellipsis$object <- object
        do.call(paste("fdr", ctrl@fdr.model, sep = "_"), ellipsis)
      } else {
        cat(paste0("[", Sys.time(), "]   getting group quants differential expression summaries...\n"))
        group_quants_de(object, summary = T, as.data.table = T)
      }
    }

    # markdown folder
    report.index <- list()
    root <- file.path(filepath(object), "markdown", "study")
    dir.create(root, recursive = T)
    group <- control(root(object))@group[1]

    # calculate plot limits
    if (ctrl@plots == T) {
      cat(paste0("[", Sys.time(), "]   calculating plot limits...\n"))
      lims <- list()
      if ("group.quants.de" %in% ctrl@plot) lims$group.quants.de <- limits_dists(group_quants_de(object, summary = T, as.data.table = T))
      saveRDS(lims, file.path(filepath(object), "limits.rds"))
    }

    # write out and plot group fdr
    if (file.exists(file.path(filepath(object), "group.quants.fdr.fst"))) {
      cat(paste0("[", Sys.time(), "]   writing group quants differential expression output...\n"))

      DTs.fdr <- fst::read.fst(file.path(filepath(object), "group.quants.fdr.fst"), as.data.table = T)
      DTs.fdr <- split(DTs.fdr, drop = T, by = "Batch")
      for (batch in names(DTs.fdr)) {
        cat(paste0("[", Sys.time(), "]    batch=", batch, "...\n"))

        name <- ifelse(name(object) == "default", "", paste0("__", name(object)))

        # write
        fwrite(DTs.fdr[[batch]], file.path(dirname(filepath(object)), "csv", paste0(tolower(group), "_fdr__", gsub("\\.", "_", batch), name, ".csv")))

        # plot
        if ("group.quants.de.batch" %in% ctrl@plot) {
          #cat(paste0("[", Sys.time(), "]     generating group quants differential expression summary...\n"))
          #group_quants_de(object, summary = T, as.data.table = T)

          cat(paste0("[", Sys.time(), "]     generating group quants differential expression plot...\n"))
          text <- paste0(group, " differential expression for '", gsub("\\.", "' effect, comparison '", batch), "'", name)
          report.index$assay.stdevs <- data.table(
            section = "Study-level", section.order = 0, item = text, item.order = 75000,
            item.href = generate_markdown(
              object,
              plot_group_quants_fdr(object),
              root, paste0("seamass_delta__", name(object), "__group_fdr__", gsub("\\.", "_", batch)),
              text
            )
          )
        }
      }
    }

    # zip
    render_markdown(object, root)

    # save index
    fst::write.fst(rbindlist(report.index), file.path(filepath(object), "report", "study.report.fst"))

    increment_completed(filepath(object))
  }

  return(invisible(NULL))
})



#if ("de.group.quants" %in% ctrl@plot) {
#  cat(paste0("[", Sys.time(), "]   plotting standardised group deviations dea...\n"))
#  DTs <- de_group_quants(object, as.data.table = T)
#  plot_de_group_quants(object, DT, limits_dists(DT, include.zero = T), file = file.path(dirname(filepath(object)), "output", basename(filepath(object)), "log2_de_group_quants.pdf"))
#  rm(DT)
#}

# plot fdr
#plot_fdr(DTs.fdr[[name]])
#ggplot2::ggsave(file.path(dirname(filepath(object)), "output", basename(filepath(object)), paste0("log2_fdr_group_quants.", gsub("\\.", "_", name), ".pdf")), width = 8, height = 8)
# plot volcano
#plot_volcano(DTs.fdr[[name]])
#ggplot2::ggsave(file.path(dirname(filepath(object)), "output", basename(filepath(object)), paste0("log2_fdr_group_quants_volcano.", gsub("\\.", "_", name), ".pdf")), width = 8, height = 8)
# plot fc
#plot_volcano(DTs.fdr[[name]], stdev.col = "s", x.col = "m", y.col = "s")
#ggplot2::ggsave(file.path(dirname(filepath(object)), "output", basename(filepath(object)), paste0("log2_fdr_group_quants_fc.", gsub("\\.", "_", name), ".pdf")), width = 8, height = 8)

# if ("fdr.group.quants" %in% ctrl@plot) {
#   DTs <- list(
#     DTs.fdr[[name]],
#     de_group_quants(object, unique(DTs.fdr[[name]][, .(Effect, Contrast, Baseline)]), as.data.table = T)
#   )
#   plot_fdr_group_quants(object, DTs, limits_dists(DTs, include.zero = T), file = file.path(dirname(filepath(object)), "output", basename(filepath(object)), paste0("log2_fdr_group_quants_dists.", gsub("\\.", "_", name), ".pdf")))
#   rm(DTs)
# }


#if ("de.component.deviations" %in% ctrl@plot) {
#cat(paste0("[", Sys.time(), "]    plotting component deviations differential expression...\n"))
#parallel_lapply(groups(object, as.data.table = T)[, Group], function(item, object) {
#  item <- substr(item, 0, 60)
#  plot_de_component_deviations(object, item, file = file.path(dirname(filepath(object)), "output", basename(filepath(object)), "log2_de_component_deviations", paste0(item, ".pdf")))
#}, nthread = ctrl@nthread)
#}

# component deviation dea
# if (ctrl@component.deviations == T && control(object@fit)@component.model == "independent") {
#   ellipsis$type <- "component.deviations"
#   do.call(paste("dea", ctrl@model, sep = "_"), ellipsis)
# }

# if (ctrl@component.deviations == T && control(object@fit)@component.model == "independent") {
#   # summarise component deviation de and perform fdr correction
#   if (file.exists(file.path(filepath(object), "component.deviations.index.fst"))) {
#     if(ctrl@fdr.model != "") {
#       ellipsis <- ctrl@ellipsis
#       ellipsis$object <- object
#       ellipsis$type <- "component.deviations"
#       do.call(paste("fdr", ctrl@fdr.model, sep = "_"), ellipsis)
#     } else {
#       cat(paste0("[", Sys.time(), "]   getting component deviations differential expression summaries...\n"))
#       component_deviations_de(object, summary = T, as.data.table = T)
#     }
#   }
#
#   # write out component deviations fdr
#   if (file.exists(file.path(filepath(object), "fdr.component.deviations.fst"))) {
#     DTs.fdr <- split(fdr_component_deviations(object, as.data.table = T), drop = T, by = "Batch")
#     for (name in names(DTs.fdr)) {
#       # save
#       fwrite(DTs.fdr[[name]], file.path(dirname(filepath(object)), "output", basename(filepath(object)), paste0("log2_fdr_component_deviations.", gsub("\\.", "_", name), ".csv")))
#       # plot fdr
#       plot_fdr(DTs.fdr[[name]])
#       ggplot2::ggsave(file.path(dirname(filepath(object)), "output", basename(filepath(object)), paste0("log2_fdr_component_deviations.", gsub("\\.", "_", name), ".pdf")), width = 8, height = 8)
#       # plot volcano
#       plot_volcano(DTs.fdr[[name]])
#       ggplot2::ggsave(file.path(dirname(filepath(object)), "output", basename(filepath(object)), paste0("log2_fdr_component_deviations_volcano.", gsub("\\.", "_", name), ".pdf")), width = 8, height = 8)
#       # plot fc
#       plot_volcano(DTs.fdr[[name]], stdev.col = "s", x.col = "m", y.col = "s")
#       ggplot2::ggsave(file.path(dirname(filepath(object)), "output", basename(filepath(object)), paste0("log2_fdr_component_deviations_fc.", gsub("\\.", "_", name), ".pdf")), width = 8, height = 8)
#     }
#   }