Merge pull request #118 from UMCUGenetics/revert-95-feature/refactor_DIMS_GenerateBreaks

Revert "Feature/refactor_DIMS_GenerateBreaks"

Author: mraves2

DIMS/GenerateBreaks.R

@@ -5,19 +5,56 @@ suppressPackageStartupMessages(library("xcms"))
 cmd_args <- commandArgs(trailingOnly = TRUE)
 
 filepath <- cmd_args[1]
-trim <- as.numeric(cmd_args[2])
-resol <- as.numeric(cmd_args[3])
+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
+breaks_fwhm <- NULL
+breaks_fwhm_avg <- NULL
+bins <- NULL
 
 # read in mzML file
 raw_data <- suppressMessages(xcms::xcmsRaw(filepath))
 
-# get trim parameters and save them to file
-get_trim_parameters(raw_data@scantime, raw_data@polarity, trim)
+# Get time values for positive and negative scans
+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
 
-# create breaks of bins for intensities. Bin size is a function of fwhm which is a function of m/z
-get_breaks_for_bins(raw_data$mzrange, resol)
+# trim (remove) scans at the start and end for negative
+trim_left_neg  <- round(neg_times[length(neg_times) * trim])
+trim_right_neg <- round(neg_times[length(neg_times) * (1 - trim)])
 
-# Determine maximum m/z and save to file
+# Mass range m/z
+low_mz  <- raw_data@mzrange[1]
 high_mz <- raw_data@mzrange[2]
-save(high_mz, file = "highest_mz.RData")
 
+# determine number of segments (bins)
+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]
+  resol_mz <- resol * (1 / sqrt(2) ^ (log2(start_segment / 200)))
+  fwhm_segment <- start_segment / resol_mz
+  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
+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")

DIMS/GenerateBreaks.nf

@@ -7,13 +7,14 @@ process GenerateBreaks {
     input:
        tuple(val(file_id), path(mzML_file))
 
+
     output:
        path('breaks.fwhm.RData'), emit: breaks
        path('trim_params.RData'), emit: trim_params
        path('highest_mz.RData'), emit: highest_mz
 
     script:
         """
-        Rscript ${baseDir}/CustomModules/DIMS/GenerateBreaks.R $mzML_file $params.trim $params.resolution 
+        Rscript ${baseDir}/CustomModules/DIMS/GenerateBreaks.R $mzML_file ./ $params.trim $params.resolution 
         """
 }