Feature/DIMS_v3.5_harmonize_Rscripts

.github/workflows/dims_lint.yml

@@ -13,7 +13,7 @@ jobs:
     runs-on: ubuntu-latest
 
     container:
-      image: docker://umcugenbioinf/dims:1.3
+      image: ghcr.io/umcugenetics/dims:v1.4.0
 
     defaults:
       run:

.github/workflows/dims_test.yml

@@ -12,7 +12,7 @@ jobs:
     runs-on: ubuntu-latest
 
     container:
-      image: docker://umcugenbioinf/dims:1.3
+      image: ghcr.io/umcugenetics/dims:v1.4.0
 
     defaults:
         run:

DIMS/AssignToBins.R

@@ -9,24 +9,22 @@ breaks_filepath <- cmd_args[2]
 trim_parameters_filepath <- cmd_args[3]
 resol <- as.numeric(cmd_args[4])
 
-# load breaks_file: contains breaks_fwhm, breaks_fwhm_avg,
+# Initialize
+options(digits = 16)
+dims_thresh <- 100
+pos_results <- NULL
+neg_results <- NULL
+
+# load breaks_file: contains breaks_fwhm & breaks_fwhm_avg
 load(breaks_filepath)
+pos_bins <- rep(0, length(breaks_fwhm) - 1)
+neg_bins <- pos_bins
 # load trim parameters file: contains trim_left_neg, trim_left_pos, trim_right_neg & trim_right_pos
 load(trim_parameters_filepath)
 
-# get sample name
+# get name of the technical replicate
 techrep_name <- sub("\\..*$", "", basename(mzml_filepath))
 
-options(digits = 16)
-
-# Initialize
-pos_results <- NULL
-neg_results <- NULL
-bins <- rep(0, length(breaks_fwhm) - 1)
-pos_bins <- bins
-neg_bins <- bins
-dims_thresh <- 100
-
 # read in the data for 1 sample
 raw_data <- suppressMessages(xcms::xcmsRaw(mzml_filepath))
 
@@ -47,25 +45,20 @@ tic_intensity_pos <- tic_intensity_persample[tic_intensity_persample[ , "retenti
                                              tic_intensity_persample[ , "retention_time"] < max(pos_times_trimmed), ]
 tic_intensity_neg <- tic_intensity_persample[tic_intensity_persample[ , "retention_time"] > min(neg_times_trimmed) &
                                              tic_intensity_persample[ , "retention_time"] < max(neg_times_trimmed), ]
+
 # calculate weighted mean of intensities for pos and neg separately
 mean_pos <- weighted.mean(tic_intensity_pos[ , "tic_intensity"], tic_intensity_pos[ , "tic_intensity"])
 mean_neg <- weighted.mean(tic_intensity_neg[ , "tic_intensity"], tic_intensity_neg[ , "tic_intensity"])
 # intensity per scan should be at least 80% of weighted mean
 dims_thresh_pos <- 0.8 * mean_pos
 dims_thresh_neg <- 0.8 * mean_neg
 
-# Generate an index with which to select values for each mode
-#pos_index <- which(raw_data_matrix[, "time"] %in% pos_times)
-#neg_index <- which(raw_data_matrix[, "time"] %in% neg_times)
 # select only data from scans which pass the dims_thresh_pos and *_neg filter
 pos_times_pass <- tic_intensity_pos[which(tic_intensity_pos[ , "tic_intensity"] > dims_thresh_pos), "retention_time"]
 neg_times_pass <- tic_intensity_neg[which(tic_intensity_neg[ , "tic_intensity"] > dims_thresh_neg), "retention_time"]
-# Generate an index with which to select values for each mode
-pos_index <- which(raw_data_matrix[, "time"] %in% pos_times_pass)
-neg_index <- which(raw_data_matrix[, "time"] %in% neg_times_pass)
 # Separate each mode into its own matrix
-pos_raw_data_matrix <- raw_data_matrix[pos_index, ]
-neg_raw_data_matrix <- raw_data_matrix[neg_index, ]
+pos_raw_data_matrix <- raw_data_matrix[which(raw_data_matrix[, "time"] %in% pos_times_pass), ]
+neg_raw_data_matrix <- raw_data_matrix[which(raw_data_matrix[, "time"] %in% neg_times_pass), ]
 
 # Get index for binning intensity values
 bin_indices_pos <- cut(
@@ -83,8 +76,7 @@ bin_indices_neg <- cut(
   labels = FALSE
 )
 
-# Get the list of intensity values for each bin, and add the
-# intensity values which are in the same bin
+# Get the list of intensity values for each bin and sum the intensities
 if (nrow(pos_raw_data_matrix) > 0) {
   # set NA in intensities to zero
   pos_raw_data_matrix[is.na(pos_raw_data_matrix[, "intensity"]), "intensity"] <- 0

DIMS/AssignToBins.nf

@@ -1,7 +1,7 @@
 process AssignToBins {
     tag "DIMS AssignToBins ${file_id}"
     label 'AssignToBins'
-    container = 'docker://umcugenbioinf/dims:1.3'
+    container = 'ghcr.io/umcugenetics/dims:v1.4.0'
     shell = ['/bin/bash', '-euo', 'pipefail']
 
     input:

DIMS/AveragePeaks.R

@@ -1,3 +1,4 @@
+# load required packages
 library(dplyr)
 
 # define parameters
@@ -12,9 +13,9 @@ tech_reps <- strsplit(techreps, ";")[[1]]
 # load in function scripts
 source(paste0(preprocessing_scripts_dir, "average_peaks_functions.R"))
 
-# Initialize per sample
+# Initialize
+options(digits = 16)
 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")
 
@@ -25,13 +26,15 @@ for (file_nr in 1:length(tech_reps)) {
   # combine data for all technical replicates
   peaklist_allrepl <- rbind(peaklist_allrepl, tech_repl)
 }
-# sort on mass
+# make sure mass and intensity columns are numeric
 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) 
+# sort on mass
 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

@@ -1,7 +1,7 @@
 process AveragePeaks {
     tag "DIMS AveragePeaks"
     label 'AveragePeaks'
-    container = 'docker://umcugenbioinf/dims:1.3'
+    container = 'ghcr.io/umcugenetics/dims:v1.4.0'
     shell = ['/bin/bash', '-euo', 'pipefail']
 
     input:

DIMS/CollectAveraged.R

@@ -1,7 +1,5 @@
-# define parameters
-cmd_args <- commandArgs(trailingOnly = TRUE)
-
-scripts_dir <- cmd_args[1]
+# Initialize
+options(digits = 16)
 
 # for each scan mode, collect all averaged peak lists per biological sample
 scanmodes <- c("positive", "negative")

DIMS/CollectAveraged.nf

@@ -1,7 +1,7 @@
 process CollectAveraged {
     tag "DIMS CollectAveraged"
     label 'CollectAveraged'
-    container = 'docker://umcugenbioinf/dims:1.3'
+    container = 'ghcr.io/umcugenetics/dims:v1.4.0'
     shell = ['/bin/bash', '-euo', 'pipefail']
 
     input:

DIMS/CollectFilled.R

@@ -2,44 +2,46 @@
 cmd_args <- commandArgs(trailingOnly = TRUE)
 
 preprocessing_scripts_dir <- cmd_args[1]
-ppm <- as.numeric(cmd_args[2])
-z_score <- as.numeric(cmd_args[3])
+z_score <- as.numeric(cmd_args[2])
 
 source(paste0(preprocessing_scripts_dir, "collect_filled_functions.R"))
 
+# Initialize
+options(digits = 16)
+
 # for each scan mode, collect all filled peak group lists
 scanmodes <- c("positive", "negative")
 for (scanmode in scanmodes) {
   # get list of files
   filled_files <- list.files("./", full.names = TRUE, pattern = paste0(scanmode, "_identified_filled"))
   # load files and combine into one object
-  outlist_total <- NULL
+  peakgroup_list_total <- NULL
   for (file_nr in seq_along(filled_files)) {
     peakgrouplist_filled <- get(load(filled_files[file_nr]))
-    outlist_total <- rbind(outlist_total, peakgrouplist_filled)
+    peakgroup_list_total <- rbind(peakgroup_list_total, peakgrouplist_filled)
   }
   # remove duplicates; peak groups with exactly the same m/z
-  outlist_total <- merge_duplicate_rows(outlist_total)
+  peakgroup_list_total <- merge_duplicate_rows(peakgroup_list_total)
   # sort on mass
-  outlist_total <- outlist_total[order(outlist_total[, "mzmed.pgrp"]), ]
+  peakgroup_list_total <- peakgroup_list_total[order(peakgroup_list_total[, "mzmed.pgrp"]), ]
   # load replication pattern
   pattern_file <- paste0(scanmode, "_repl_pattern.RData")
   repl_pattern <- get(load(pattern_file))
   # calculate Z-scores
   if (z_score == 1) {
-    outlist_stats <- calculate_zscores_peakgrouplist(outlist_total)
+    peakgroup_list_stats <- calculate_zscores_peakgrouplist(peakgroup_list_total)
   } else {
-    outlist_stats <- outlist_total
+    peakgroup_list_stats <- peakgroup_list_total
   }
   # calculate ppm deviation
-  outlist_withppm <- calculate_ppm_deviation(outlist_stats)
+  peakgroup_list_withppm <- calculate_ppm_deviation(peakgroup_list_stats)
   #  put columns in correct order
-  outlist_ident <- order_columns_peakgrouplist(outlist_withppm)
+  peakgroup_list_ident <- order_columns_peakgrouplist(peakgroup_list_withppm)
 
   # 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)
+  peakgroup_list_ident <- peakgroup_list_ident[, -which(colnames(peakgroup_list_ident) %in% remove_columns)]
+  write.table(peakgroup_list_ident, file = paste0("peakgroup_list_identified_", scanmode, ".txt"), sep = "\t", row.names = FALSE)
   # export output in RData format
-  save(outlist_ident, file = paste0("outlist_identified_", scanmode, ".RData"))
+  save(peakgroup_list_ident, file = paste0("peakgroup_list_identified_", scanmode, ".RData"))
 }

DIMS/CollectFilled.nf

@@ -1,19 +1,19 @@
 process CollectFilled {
     tag "DIMS CollectFilled"
     label 'CollectFilled'
-    container = 'docker://umcugenbioinf/dims:1.3'
+    container = 'ghcr.io/umcugenetics/dims:v1.4.0'
     shell = ['/bin/bash', '-euo', 'pipefail']
 
     input:
        path(filled_files)
        each path(replication_pattern)
 
     output:
-       path('outlist*.txt')
-       path('outlist*.RData'), emit: filled_pgrlist
+       path('peakgroup_list*.txt')
+       path('peakgroup_list*.RData'), emit: filled_pgrlist
 
     script:
         """
-        Rscript ${baseDir}/CustomModules/DIMS/CollectFilled.R $params.preprocessing_scripts_dir $params.ppm $params.zscore
+        Rscript ${baseDir}/CustomModules/DIMS/CollectFilled.R $params.preprocessing_scripts_dir $params.zscore
         """
 }

DIMS/CollectSumAdducts.R

@@ -1,5 +1,4 @@
-## Combining all AdductSums part files for each scanmode and
-# combine intensities if present in both scanmodes
+# load required packages
 suppressMessages(library("dplyr"))
 
 # define parameters

DIMS/EvaluateTics.R

@@ -1,4 +1,4 @@
-# load packages
+# load required packages
 library("ggplot2")
 library("gridExtra")
 
@@ -18,6 +18,9 @@ preprocessing_scripts_dir <- cmd_args[7]
 # load functions
 source(paste0(preprocessing_scripts_dir, "evaluate_tics_functions.R"))
 
+# Initialize
+options(digits = 16)
+
 # load init_file: contains repl_pattern
 load(init_file)
 
@@ -44,7 +47,7 @@ print(remove_tech_reps)
 remove_neg <- remove_tech_reps$neg
 repl_pattern_filtered <- remove_from_repl_pattern(remove_neg, repl_pattern, nr_replicates)
 save(repl_pattern_filtered, file = "negative_repl_pattern.RData")
-# get an overview of suitable technical replicates for both negative mode
+# get an overview of suitable technical replicates for negative scan mode
 allsamples_techreps_neg <- get_overview_tech_reps(repl_pattern_filtered, "negative")
 
 # positive scan mode
@@ -63,9 +66,7 @@ write.table(allsamples_techreps_both_scanmodes,
             sep = ","
 )
 
-
-## generate TIC plots
-# get all txt files
+# generate TIC plots using all tic files
 tic_files <- list.files("./", full.names = TRUE, pattern = "*TIC.txt")
 all_samps <- sub("_TIC\\..*$", "", basename(tic_files))
 
@@ -141,4 +142,3 @@ tic_plot_pdf <- marrangeGrob(
 ggsave(filename = paste0(run_name, "_TICplots.pdf"),
        tic_plot_pdf, width = 21, height = 29.7, units = "cm")
 
-

DIMS/EvaluateTics.nf

@@ -1,7 +1,7 @@
 process EvaluateTics {
     tag "DIMS EvaluateTics"
     label 'EvaluateTics'
-    container = 'docker://umcugenbioinf/dims:1.3'
+    container = 'ghcr.io/umcugenetics/dims:v1.4.0'
     shell = ['/bin/bash', '-euo', 'pipefail']
 
     input:

DIMS/FillMissing.R

@@ -8,6 +8,9 @@ thresh <- as.numeric(cmd_args[3])
 # load in function scripts
 source(paste0(preprocessing_scripts_dir, "fill_missing_functions.R"))
 
+# Initialize
+options(digits = 16)
+
 # determine scan mode
 if (grepl("_pos", peakgrouplist_file)) {
   scanmode <- "positive"
@@ -30,3 +33,4 @@ outputfile_name <- gsub(".RData", "_filled.RData", peakgrouplist_file)
 
 # save output
 save(peakgrouplist_filled, file = outputfile_name)
+

DIMS/FillMissing.nf

@@ -1,7 +1,7 @@
 process FillMissing {
     tag "DIMS FillMissing ${peakgrouplist_file}"
     label 'FillMissing'
-    container = 'docker://umcugenbioinf/dims:1.3'
+    container = 'ghcr.io/umcugenetics/dims:v1.4.0'
     shell = ['/bin/bash', '-euo', 'pipefail']
 
     input:

DIMS/GenerateBreaks.R

@@ -5,18 +5,13 @@ suppressPackageStartupMessages(library("xcms"))
 cmd_args <- commandArgs(trailingOnly = TRUE)
 
 filepath <- cmd_args[1]
-outdir <- cmd_args[2]
-trim <- as.numeric(cmd_args[3])
-resol <- as.numeric(cmd_args[4])
-
-# initialize
-trim_left_pos <- NULL
-trim_right_pos <- NULL
-trim_left_neg <- NULL
-trim_right_neg <- NULL
+trim <- as.numeric(cmd_args[2])
+resol <- as.numeric(cmd_args[3])
+
+# Initialize
+options(digits = 16)
 breaks_fwhm <- NULL
 breaks_fwhm_avg <- NULL
-bins <- NULL
 
 # read in mzML file
 raw_data <- suppressMessages(xcms::xcmsRaw(filepath))
@@ -26,8 +21,8 @@ pos_times <- raw_data@scantime[raw_data@polarity == "positive"]
 neg_times <- raw_data@scantime[raw_data@polarity == "negative"]
 
 # trim (remove) scans at the start and end for positive
-trim_left_pos  <- round(pos_times[length(pos_times) * (trim * 1.5)]) # 15% aan het begin
-trim_right_pos <- round(pos_times[length(pos_times) * (1 - (trim * 0.5))]) # 5% aan het eind
+trim_left_pos  <- round(pos_times[length(pos_times) * (trim * 1.5)]) # 15% at the start
+trim_right_pos <- round(pos_times[length(pos_times) * (1 - (trim * 0.5))]) # 5% at the end
 
 # trim (remove) scans at the start and end for negative
 trim_left_neg  <- round(neg_times[length(neg_times) * trim])
@@ -37,24 +32,28 @@ trim_right_neg <- round(neg_times[length(neg_times) * (1 - trim)])
 low_mz  <- raw_data@mzrange[1]
 high_mz <- raw_data@mzrange[2]
 
-# determine number of segments (bins)
+# determine number of segments
 nr_segments <- 2 * (high_mz - low_mz)
 segment <- seq(from = low_mz, to = high_mz, length.out = nr_segments + 1)
 
-# determine start and end of each bin.
-for (i in 1:nr_segments) {
-  start_segment <- segment[i]
-  end_segment <- segment[i+1]
+# create bins for each segment
+for (segment_index in 1:nr_segments) {
+  # determine start and end of each bin.
+  start_segment <- segment[segment_index]
+  end_segment <- segment[segment_index + 1]
+  # determine resolution for this mz range
   resol_mz <- resol * (1 / sqrt(2) ^ (log2(start_segment / 200)))
+  # determine full width at half maximum of peaks in this segment
   fwhm_segment <- start_segment / resol_mz
+  # determine boundaries of bins
   breaks_fwhm <- c(breaks_fwhm, seq(from = (start_segment + fwhm_segment), to = end_segment, by = 0.2 * fwhm_segment))
   # average the m/z instead of start value
   range <- seq(from = (start_segment + fwhm_segment), to = end_segment, by = 0.2 * fwhm_segment)
   delta_mz <- range[2] - range[1]
   breaks_fwhm_avg <- c(breaks_fwhm_avg, range + 0.5 * delta_mz)
 }
 
-# generate output file
+# generate output files
 save(breaks_fwhm, breaks_fwhm_avg, file = "breaks.fwhm.RData")
 save(trim_left_pos, trim_right_pos, trim_left_neg, trim_right_neg, file = "trim_params.RData")
 save(high_mz, file = "highest_mz.RData")