Add human readable log to DESeq2 pipelines

michael-kotliar · michael-kotliar · commit 062f5901bdd3 · 2025-08-27T15:17:29.000-04:00
diff --git a/tools/deseq-lrt-step-1.cwl b/tools/deseq-lrt-step-1.cwl
@@ -169,6 +169,11 @@ outputs:
     outputBinding:
       glob: "*_all_contrasts.tsv"
 
+  human_log:
+    type: File?
+    outputBinding:
+      glob: "error_report.txt"
+
   stdout_log:
     type: stdout
 
@@ -177,4 +182,4 @@ outputs:
 
 baseCommand: [run_deseq_lrt_step_1.sh]
 stdout: error_msg.txt
-stderr: error_report.txt
+stderr: deseq_lrt_step_1_stderr.log
diff --git a/tools/deseq-lrt-step-2.cwl b/tools/deseq-lrt-step-2.cwl
@@ -125,6 +125,11 @@ outputs:
     outputBinding:
       glob: "*_diff_expr.tsv"
 
+  human_log:
+    type: File?
+    outputBinding:
+      glob: "error_report.txt"
+
   stdout_log:
     type: stdout
 
@@ -133,4 +138,4 @@ outputs:
 
 baseCommand: [run_deseq_lrt_step_2.sh]
 stdout: error_msg.txt
-stderr: error_report.txt
+stderr: deseq_lrt_step_2_stderr.log
diff --git a/tools/dockerfiles/scripts/diff_tools/modules/io.R b/tools/dockerfiles/scripts/diff_tools/modules/io.R
@@ -215,51 +215,17 @@ export_morpheus_html_heatmap <- function(gct_location, rootname, color_scheme=NU
     )
 }
 
-load_metadata <- function(location, aliases, batch=NULL){
+load_metadata <- function(location){
     metadata <- utils::read.table(
         location,
         sep=get_file_type(location),
         header=TRUE,
         stringsAsFactors=FALSE,
         row.names=1,
         quote="",
-        colClasses="character"                                               # to prevent loading numbers as numerical types
-    ) %>% dplyr::mutate_at(base::colnames(.), forcats::fct_inorder)                       # assigns levels bases on the first appearance in the metadata
-    base::rownames(metadata) <- adjust_names(base::rownames(metadata))                   # to make it correspond to the aliases
-
-    if (base::anyDuplicated(base::rownames(metadata)) > 0) {
-        base::print(
-            base::paste(
-                "Exiting: the first columns in",
-                "the metadata shouldn't include",
-                "any duplicates."
-            )
-        )
-        base::quit(save="no", status=1, runLast=FALSE)
-    }
-
-    if (!generics::setequal(base::rownames(metadata), aliases)){                        # we already cheched that neither metadata nor aliases have duplicates
-        base::print(
-            base::paste(
-                "Exiting: the sample names in the",
-                "metadata should correposnd to the",
-                "sample names provided in the aliases",
-                "parameter."
-            )
-        )
-        base::quit(save="no", status=1, runLast=FALSE)
-    }
-
-    if (!is.null(batch) && !(batch %in% base::colnames(metadata))){
-        base::print(
-            base::paste(
-                "Exiting: the batch parameter",
-                "should be used as one of the",
-                "samples metadata column."
-            )
-        )
-        base::quit(save="no", status=1, runLast=FALSE)
-    }
+        colClasses="character"                                                          # to prevent loading numbers as numerical types
+    ) %>% dplyr::mutate_at(base::colnames(.), forcats::fct_inorder)                     # assigns levels bases on the first appearance in the metadata
+    base::rownames(metadata) <- adjust_names(base::rownames(metadata))                  # to make it correspond to the aliases
 
     for (i in 1:length(base::colnames(metadata))){
         base::print(
@@ -269,7 +235,6 @@ load_metadata <- function(location, aliases, batch=NULL){
             )
         )
     }
-
     return (metadata)
 }
 
diff --git a/tools/dockerfiles/scripts/diff_tools/modules/logger.R b/tools/dockerfiles/scripts/diff_tools/modules/logger.R
@@ -0,0 +1,27 @@
+import("futile.logger", attach=FALSE)
+
+export(
+    "setup",
+    "info"
+)
+
+
+setup <- function(location, name="hlog", header=NULL, level=futile.logger::INFO, format="~m"){
+    futile.logger::flog.logger(
+        name,
+        level,
+        appender=futile.logger::appender.file(location),
+        layout=futile.logger::layout.format(format)
+    )
+    if (!is.null(header)){
+        futile.logger:::.log_level(
+            header, level=level, name=name,
+            capture=FALSE, logger=NULL
+        )
+    }
+}
+
+info <- function(message, name="hlog", skip_stdout=FALSE){
+    futile.logger::flog.info(message, name=name) 
+    if (!skip_stdout) {print(message)}
+}
diff --git a/tools/dockerfiles/scripts/diff_tools/run_deseq_lrt_step_1.R b/tools/dockerfiles/scripts/diff_tools/run_deseq_lrt_step_1.R
@@ -15,12 +15,13 @@ options(error=function(){traceback(3); quit(save="no", status=1, runLast=FALSE)}
 HERE <- (function() {return (dirname(sub("--file=", "", commandArgs(trailingOnly=FALSE)[grep("--file=", commandArgs(trailingOnly=FALSE))])))})()
 suppressMessages(io <- modules::use(file.path(HERE, "modules/io.R")))
 suppressMessages(analyses <- modules::use(file.path(HERE, "modules/analyses.R")))
+suppressMessages(logger <- modules::use(file.path(HERE, "modules/logger.R")))
 
 
 assert_args <- function(args){
 
     if (length(args$expression) != length(args$aliases)) {
-        print(
+        logger$info(
             paste(
                 "Exiting: the --expression and --aliases",
                 "parameters should have the same length."
@@ -30,7 +31,7 @@ assert_args <- function(args){
     }
 
     if (is.null(args$batch) != is.null(args$correction)){
-        print(
+        logger$info(
             paste(
                 "Exiting: for batch correction both",
                 "--batch and --correction parameters",
@@ -44,7 +45,7 @@ assert_args <- function(args){
         grepl("\\*", args$design) ||
         grepl("\\*", args$reduced)
     ){
-        print(
+        logger$info(
             paste(
                 "Exiting: both of the --design and",
                 "--reduced formulas should be povided",
@@ -60,7 +61,7 @@ assert_args <- function(args){
         !is.null(args$batch) &&
         !(args$batch %in% unique(all.vars(args$design)))
     ){
-        print(
+        logger$info(
             paste(
                 "Exiting: the --batch parameter should",
                 "also be used in the --design formula."
@@ -71,7 +72,7 @@ assert_args <- function(args){
 
     args$aliases <- io$adjust_names(args$aliases)
     if (anyDuplicated(args$aliases) > 0) {
-        print(
+        logger$info(
             paste(
                 "Exiting: the sample names provided",
                 "in the --aliases parameter shouldn't",
@@ -276,7 +277,12 @@ get_args <- function() {
         type="integer",
         default=1
     )
-    args <- assert_args(parser$parse_args(commandArgs(trailingOnly=TRUE)))
+    args <- parser$parse_args(commandArgs(trailingOnly=TRUE))
+    logger$setup(
+        file.path(dirname(ifelse(args$output == "", "./", args$output)), "error_report.txt"),
+        header="DESeq2 LRT Step 1 (run_deseq_lrt_step_1.R)"
+    )
+    args <- assert_args(args)
     return(args)
 }
 
@@ -289,13 +295,43 @@ print(paste("Setting parallelizations to", args$cpus, "cores."))
 register(MulticoreParam(args$cpus))
 
 print("Loading metadata.")
-metadata <- io$load_metadata(
-    location=args$metadata,
-    aliases=args$aliases,
-    batch=args$batch
-)
+metadata <- io$load_metadata(args$metadata)
 print(metadata)
 
+if (anyDuplicated(rownames(metadata)) > 0) {
+    logger$info(
+        paste(
+            "Exiting: the first columns in",
+            "the metadata shouldn't include",
+            "any duplicates."
+        )
+    )
+    quit(save="no", status=1, runLast=FALSE)
+}
+
+if (!setequal(rownames(metadata), args$aliases)){                               # we already cheched that neither metadata nor aliases have duplicates
+    logger$info(
+        paste(
+            "Exiting: the sample names in the",
+            "metadata should correspond to the",
+            "sample names provided in the aliases",
+            "parameter."
+        )
+    )
+    quit(save="no", status=1, runLast=FALSE)
+}
+
+if (!is.null(args$batch) && !(args$batch %in% colnames(metadata))){
+    logger$info(
+        paste(
+            "Exiting: the batch parameter",
+            "should be used as one of the",
+            "samples metadata column."
+        )
+    )
+    quit(save="no", status=1, runLast=FALSE)
+}
+
 print("Loading gene expression.")
 expression_data <- io$load_expression_data(
     locations=args$expression,
diff --git a/tools/dockerfiles/scripts/diff_tools/run_deseq_lrt_step_2.R b/tools/dockerfiles/scripts/diff_tools/run_deseq_lrt_step_2.R
@@ -13,6 +13,7 @@ HERE <- (function() {return (dirname(sub("--file=", "", commandArgs(trailingOnly
 suppressMessages(io <- modules::use(file.path(HERE, "modules/io.R")))
 suppressMessages(analyses <- modules::use(file.path(HERE, "modules/analyses.R")))
 suppressMessages(graphics <- modules::use(file.path(HERE, "modules/graphics.R")))
+suppressMessages(logger <- modules::use(file.path(HERE, "modules/logger.R")))
 
 
 assert_args <- function(args){
@@ -138,7 +139,12 @@ get_args <- function() {
         type="integer",
         default=1
     )
-    args <- assert_args(parser$parse_args(commandArgs(trailingOnly=TRUE)))
+    args <- parser$parse_args(commandArgs(trailingOnly=TRUE))
+    logger$setup(
+        file.path(dirname(ifelse(args$output == "", "./", args$output)), "error_report.txt"),
+        header="DESeq2 LRT Step 2 (run_deseq_lrt_step_2.R)"
+    )
+    args <- assert_args(args)
     return(args)
 }
 
@@ -184,7 +190,7 @@ selected_contrasts <- all_contrasts[intersect(args$target, names(all_contrasts))
 print(selected_contrasts)
 
 if (length(selected_contrasts) == 0) {
-    print(
+    logger$info(
         paste(
             "Exiting: no target contrasts found."
         )
@@ -259,7 +265,7 @@ collected_deseq_results <- collected_deseq_results %>%
                            rownames_to_column(var="Feature")                                 # we need to to join with cluster information
 print(head(collected_deseq_results))
 if(!any(grepl("_padj$", colnames(collected_deseq_results)))){                                # to check that we have at least one contrast successfully calculated
-    print(
+    logger$info(
         paste(
             "Exiting: neither of the selected",
             "contrasts were successfully calculated."
diff --git a/workflows/deseq-lrt-step-1.cwl b/workflows/deseq-lrt-step-1.cwl
@@ -321,6 +321,14 @@ outputs:
       All contrasts produced by DESeq2 Wald tests.
       RDS format.
 
+  human_log:
+    type: File?
+    outputSource: deseq_lrt_step_1/human_log
+    label: "Human readable error log"
+    doc: |
+      Human readable error log
+      from the deseq_lrt_step_1 step.
+
   stdout_log:
     type: File
     outputSource: deseq_lrt_step_1/stdout_log
@@ -371,5 +379,6 @@ steps:
     - diff_expr_tsv
     - all_contrasts_rds
     - all_contrasts_tsv
+    - human_log
     - stdout_log
     - stderr_log
diff --git a/workflows/deseq-lrt-step-2.cwl b/workflows/deseq-lrt-step-2.cwl
@@ -231,6 +231,14 @@ outputs:
       Morpheus compatible heatmap of normalized read counts.
       GCT format.
 
+  human_log:
+    type: File?
+    outputSource: deseq_lrt_step_2/human_log
+    label: "Human readable error log"
+    doc: |
+      Human readable error log
+      from the deseq_lrt_step_2 step.
+
   stdout_log:
     type: File
     outputSource: deseq_lrt_step_2/stdout_log
@@ -268,5 +276,6 @@ steps:
     - read_counts_gct
     - read_counts_html
     - diff_expr_tsv
+    - human_log
     - stdout_log
     - stderr_log
