Release/v2.7.0 to main

.github/workflows/dims_test.yml

@@ -23,7 +23,7 @@ jobs:
         uses: actions/checkout@v4
 
       - name: Install dependencies
-        run: Rscript -e "install.packages(c('testthat', 'withr'))"
+        run: Rscript -e "install.packages(c('testthat', 'withr', 'vdiffr', 'pdftools'))"
 
       - name: Run tests
         run: Rscript tests/testthat.R

.lintr

@@ -1,6 +1,7 @@
 linters: linters_with_defaults(
     line_length_linter(127),
     object_usage_linter = NULL,
-    return_linter = NULL
+    return_linter = NULL,
+    pipe_consistency_linter = lintr::pipe_consistency_linter("auto")
   )
 encoding: "UTF-8"

DIMS/AssignToBins.R

@@ -6,11 +6,13 @@ cmd_args <- commandArgs(trailingOnly = TRUE)
 
 mzml_filepath <- cmd_args[1]
 breaks_filepath <- cmd_args[2]
-resol <- as.numeric(cmd_args[3])
+trim_parameters_filepath <- cmd_args[3]
+resol <- as.numeric(cmd_args[4])
 
 # load breaks_file: contains breaks_fwhm, breaks_fwhm_avg,
-# trim_left_neg, trim_left_pos, trim_right_neg & trim_right_pos
 load(breaks_filepath)
+# load trim parameters file: contains trim_left_neg, trim_left_pos, trim_right_neg & trim_right_pos
+load(trim_parameters_filepath)
 
 # get sample name
 techrep_name <- sub("\\..*$", "", basename(mzml_filepath))

DIMS/AssignToBins.nf

@@ -5,15 +5,15 @@ process AssignToBins {
     shell = ['/bin/bash', '-euo', 'pipefail']
 
     input:
-       tuple(val(file_id), path(mzML_file), path(breaks_file))
+       tuple(val(file_id), path(mzML_file), path(breaks_file), path(trim_params_file))
 
     output:
        path("${file_id}.RData"), emit: rdata_file
        path("${file_id}_TIC.txt"), emit: tic_txt_file
 
     script:
         """
-        Rscript ${baseDir}/CustomModules/DIMS/AssignToBins.R $mzML_file $breaks_file $params.resolution
+        Rscript ${baseDir}/CustomModules/DIMS/AssignToBins.R $mzML_file $breaks_file $trim_params_file $params.resolution
         """
 }
 

DIMS/AveragePeaks.R

@@ -0,0 +1,37 @@
+library(dplyr)
+
+# define parameters
+cmd_args <- commandArgs(trailingOnly = TRUE)
+
+sample_name <- cmd_args[1]
+techreps <- cmd_args[2]
+scanmode <- cmd_args[3]
+preprocessing_scripts_dir <- cmd_args[4]
+tech_reps <- strsplit(techreps, ";")[[1]]
+
+# load in function scripts
+source(paste0(preprocessing_scripts_dir, "average_peaks_functions.R"))
+
+# Initialize per sample
+peaklist_allrepl <- NULL
+nr_repl_persample <- 0
+averaged_peaks <- matrix(0, nrow = 0, ncol = 6) 
+colnames(averaged_peaks) <- c("samplenr", "mzmed.pkt", "fq", "mzmin.pkt", "mzmax.pkt", "height.pkt")
+
+# load RData files of technical replicates belonging to biological sample
+for (file_nr in 1:length(tech_reps)) {
+  tech_repl_file <- paste0(tech_reps[file_nr], "_", scanmode, ".RData")
+  tech_repl <- get(load(tech_repl_file))
+  # combine data for all technical replicates
+  peaklist_allrepl <- rbind(peaklist_allrepl, tech_repl)
+}
+# sort on mass
+peaklist_allrepl_df <- as.data.frame(peaklist_allrepl)
+peaklist_allrepl_df$mzmed.pkt <- as.numeric(peaklist_allrepl_df$mzmed.pkt) 
+peaklist_allrepl_df$height.pkt <- as.numeric(peaklist_allrepl_df$height.pkt) 
+peaklist_allrepl_sorted <- peaklist_allrepl_df %>% arrange(mzmed.pkt)
+
+# average over technical replicates
+averaged_peaks <- average_peaks_per_sample(peaklist_allrepl_sorted, sample_name)
+save(averaged_peaks, file = paste0("AvgPeaks_", sample_name, "_", scanmode, ".RData"))
+

DIMS/AveragePeaks.nf

@@ -0,0 +1,18 @@
+process AveragePeaks {
+    tag "DIMS AveragePeaks"
+    label 'AveragePeaks'
+    container = 'docker://umcugenbioinf/dims:1.3'
+    shell = ['/bin/bash', '-euo', 'pipefail']
+
+    input:
+       path(rdata_files)
+       tuple val(sample_id), val(tech_reps), val(scanmode)
+
+    output:
+       path 'AvgPeaks_*.RData'
+
+    script:
+        """
+        Rscript ${baseDir}/CustomModules/DIMS/AveragePeaks.R $sample_id $tech_reps $scanmode $params.preprocessing_scripts_dir
+        """
+}

DIMS/CollectAveraged.R

@@ -0,0 +1,18 @@
+# define parameters
+cmd_args <- commandArgs(trailingOnly = TRUE)
+
+scripts_dir <- cmd_args[1]
+
+# for each scan mode, collect all averaged peak lists per biological sample
+scanmodes <- c("positive", "negative")
+for (scanmode in scanmodes) {
+  # get list of files
+  filled_files <- list.files("./", full.names = TRUE, pattern = paste0(scanmode, ".RData"))
+  # load files and combine into one object
+  outlist_total <- NULL
+  for (file_nr in 1:length(filled_files)) {
+    peaklist_averaged <- get(load(filled_files[file_nr]))
+    outlist_total <- rbind(outlist_total, peaklist_averaged)
+  }
+  save(outlist_total, file = paste0("AvgPeaks_", scanmode, ".RData"))
+}

DIMS/CollectAveraged.nf

@@ -0,0 +1,17 @@
+process CollectAveraged {
+    tag "DIMS CollectAveraged"
+    label 'CollectAveraged'
+    container = 'docker://umcugenbioinf/dims:1.3'
+    shell = ['/bin/bash', '-euo', 'pipefail']
+
+    input:
+       path(averaged_files)
+
+    output:
+       path('AvgPeaks*.RData'), emit: averaged_peaks
+
+    script:
+        """
+        Rscript ${baseDir}/CustomModules/DIMS/CollectAveraged.R
+        """
+}

DIMS/CollectFilled.R

@@ -1,14 +1,11 @@
-## adapted from 10-collectSamplesFilled.R
-
 # define parameters
 cmd_args <- commandArgs(trailingOnly = TRUE)
 
-scripts_dir <- cmd_args[1]
+preprocessing_scripts_dir <- cmd_args[1]
 ppm <- as.numeric(cmd_args[2])
 z_score <- as.numeric(cmd_args[3])
 
-source(paste0(scripts_dir, "merge_duplicate_rows.R"))
-source(paste0(scripts_dir, "calculate_zscores.R"))
+source(paste0(preprocessing_scripts_dir, "collect_filled_functions.R"))
 
 # for each scan mode, collect all filled peak group lists
 scanmodes <- c("positive", "negative")
@@ -17,7 +14,7 @@ for (scanmode in scanmodes) {
   filled_files <- list.files("./", full.names = TRUE, pattern = paste0(scanmode, "_identified_filled"))
   # load files and combine into one object
   outlist_total <- NULL
-  for (file_nr in 1:length(filled_files)) {
+  for (file_nr in seq_along(filled_files)) {
     peakgrouplist_filled <- get(load(filled_files[file_nr]))
     outlist_total <- rbind(outlist_total, peakgrouplist_filled)
   }
@@ -30,37 +27,17 @@ for (scanmode in scanmodes) {
   repl_pattern <- get(load(pattern_file))
   # calculate Z-scores
   if (z_score == 1) {
-    outlist_stats <- calculate_zscores(outlist_total, adducts = FALSE)
-    nr_removed_samples <- length(which(repl_pattern[] == "character(0)"))
-    order_index_int <- order(colnames(outlist_stats)[8:(length(repl_pattern) - nr_removed_samples + 7)])
-    outlist_stats_more <- cbind(
-      outlist_stats[, 1:7],
-      outlist_stats[, (length(repl_pattern) - nr_removed_samples + 8):(length(repl_pattern) - nr_removed_samples + 8 + 6)],
-      outlist_stats[, 8:(length(repl_pattern) - nr_removed_samples + 7)][order_index_int],
-      outlist_stats[, (length(repl_pattern) - nr_removed_samples + 5 + 10):ncol(outlist_stats)]
-    )
-    # sort Z-score columns and append to peak group list
-    tmp_index <- grep("_Zscore", colnames(outlist_stats_more), fixed = TRUE)
-    tmp_index_order <- order(colnames(outlist_stats_more[, tmp_index]))
-    tmp <- outlist_stats_more[, tmp_index[tmp_index_order]]
-    outlist_stats_more <- outlist_stats_more[, -tmp_index]
-    outlist_stats_more <- cbind(outlist_stats_more, tmp)
-    outlist_total <- outlist_stats_more
+    outlist_stats <- calculate_zscores_peakgrouplist(outlist_total)
   }
-
-  # make a copy of the outlist
-  outlist_ident <- outlist_total
-  # take care of NAs in theormz_noise
-  outlist_ident$theormz_noise[which(is.na(outlist_ident$theormz_noise))] <- 0
-  outlist_ident$theormz_noise <- as.numeric(outlist_ident$theormz_noise)
-  outlist_ident$theormz_noise[which(is.na(outlist_ident$theormz_noise))] <- 0
-  outlist_ident$theormz_noise <- as.numeric(outlist_ident$theormz_noise)
+  # calculate ppm deviation
+  outlist_withppm <- calculate_ppm_deviation(outlist_stats)
+  #  put columns in correct order
+  outlist_ident <- order_columns_peakgrouplist(outlist_withppm)
 
-  # Extra output in Excel-readable format:
-  remove_columns <- c("fq.best", "fq.worst", "mzmin.pgrp", "mzmax.pgrp")
-  remove_colindex <- which(colnames(outlist_ident) %in% remove_columns)
-  outlist_ident <- outlist_ident[, -remove_colindex]
+  # generate output in Excel-readable format:
+  remove_columns <- c("mzmin.pgrp", "mzmax.pgrp")
+  outlist_ident <- outlist_ident[, -which(colnames(outlist_ident) %in% remove_columns)]
   write.table(outlist_ident, file = paste0("outlist_identified_", scanmode, ".txt"), sep = "\t", row.names = FALSE)
-  # output in RData format
+  # export output in RData format
   save(outlist_ident, file = paste0("outlist_identified_", scanmode, ".RData"))
 }

DIMS/CollectFilled.nf

@@ -14,6 +14,6 @@ process CollectFilled {
 
     script:
         """
-        Rscript ${baseDir}/CustomModules/DIMS/CollectFilled.R $params.scripts_dir $params.ppm $params.zscore
+        Rscript ${baseDir}/CustomModules/DIMS/CollectFilled.R $params.preprocessing_scripts_dir $params.ppm $params.zscore
         """
 }