diff --git a/.github/workflows/dims_lint.yml b/.github/workflows/dims_lint.yml
index 7a068e0d..7797c996 100644
--- a/.github/workflows/dims_lint.yml
+++ b/.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:
diff --git a/.github/workflows/dims_test.yml b/.github/workflows/dims_test.yml
index e8c4078c..791a750e 100644
--- a/.github/workflows/dims_test.yml
+++ b/.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:
diff --git a/DIMS/AssignToBins.R b/DIMS/AssignToBins.R
index 8b31af7e..a076c534 100644
--- a/DIMS/AssignToBins.R
+++ b/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,6 +45,7 @@ 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"])
@@ -54,18 +53,12 @@ mean_neg <- weighted.mean(tic_intensity_neg[ , "tic_intensity"], tic_intensity_n
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
diff --git a/DIMS/AssignToBins.nf b/DIMS/AssignToBins.nf
index d3bc79ae..6e28a617 100644
--- a/DIMS/AssignToBins.nf
+++ b/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:
diff --git a/DIMS/AveragePeaks.R b/DIMS/AveragePeaks.R
index 7114e3c4..b2a37041 100644
--- a/DIMS/AveragePeaks.R
+++ b/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"))
diff --git a/DIMS/AveragePeaks.nf b/DIMS/AveragePeaks.nf
index f50bd874..d1d69127 100644
--- a/DIMS/AveragePeaks.nf
+++ b/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:
diff --git a/DIMS/CollectAveraged.R b/DIMS/CollectAveraged.R
index e6466d9e..8463b6ff 100755
--- a/DIMS/CollectAveraged.R
+++ b/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")
diff --git a/DIMS/CollectAveraged.nf b/DIMS/CollectAveraged.nf
index fc65bf21..b3d34fba 100644
--- a/DIMS/CollectAveraged.nf
+++ b/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:
diff --git a/DIMS/CollectFilled.R b/DIMS/CollectFilled.R
index 4cd25fbc..6f982490 100644
--- a/DIMS/CollectFilled.R
+++ b/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"))
}
diff --git a/DIMS/CollectFilled.nf b/DIMS/CollectFilled.nf
index b0025286..8e2dafee 100644
--- a/DIMS/CollectFilled.nf
+++ b/DIMS/CollectFilled.nf
@@ -1,7 +1,7 @@
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:
@@ -9,11 +9,11 @@ process CollectFilled {
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
"""
}
diff --git a/DIMS/CollectSumAdducts.R b/DIMS/CollectSumAdducts.R
index 28b5bf0d..81c1a935 100755
--- a/DIMS/CollectSumAdducts.R
+++ b/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
diff --git a/DIMS/EvaluateTics.R b/DIMS/EvaluateTics.R
index 521430e3..f113e5f9 100644
--- a/DIMS/EvaluateTics.R
+++ b/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")
-
diff --git a/DIMS/EvaluateTics.nf b/DIMS/EvaluateTics.nf
index 2cd8bb57..a8b1538d 100644
--- a/DIMS/EvaluateTics.nf
+++ b/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:
diff --git a/DIMS/FillMissing.R b/DIMS/FillMissing.R
index a523bab6..75ff8a7a 100755
--- a/DIMS/FillMissing.R
+++ b/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)
+
diff --git a/DIMS/FillMissing.nf b/DIMS/FillMissing.nf
index 50227026..2478442c 100644
--- a/DIMS/FillMissing.nf
+++ b/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:
diff --git a/DIMS/GenerateBreaks.R b/DIMS/GenerateBreaks.R
index 007d5ab7..8461d4c3 100644
--- a/DIMS/GenerateBreaks.R
+++ b/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,16 +32,20 @@ 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)
@@ -54,7 +53,7 @@ for (i in 1:nr_segments) {
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")
diff --git a/DIMS/GenerateBreaks.nf b/DIMS/GenerateBreaks.nf
index c486010d..e9e3452f 100644
--- a/DIMS/GenerateBreaks.nf
+++ b/DIMS/GenerateBreaks.nf
@@ -1,7 +1,7 @@
process GenerateBreaks {
tag "DIMS GenerateBreaks"
label 'GenerateBreaks'
- container = 'docker://umcugenbioinf/dims:1.3'
+ container = 'ghcr.io/umcugenetics/dims:v1.4.0'
shell = ['/bin/bash', '-euo', 'pipefail']
input:
@@ -15,6 +15,6 @@ process GenerateBreaks {
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
"""
}
diff --git a/DIMS/GenerateExcel.R b/DIMS/GenerateExcel.R
index f33c27de..910f1e86 100644
--- a/DIMS/GenerateExcel.R
+++ b/DIMS/GenerateExcel.R
@@ -18,17 +18,12 @@ path_metabolite_groups <- cmd_args[5]
# load in function scripts
source(paste0(export_scripts_dir, "generate_excel_functions.R"))
-# set the number of digits for floats
+# Initialize
options(digits = 16)
-
-# Initialise
plot <- TRUE
export <- TRUE
control_label <- "C"
case_label <- "P"
-
-# setting outdir to export files to the working directory
-outdir <- "./"
# percentage of outliers to remove from calculation of robust scaler
perc <- 5
# Z-score for removing outliers with grubbs test
@@ -65,7 +60,7 @@ openxlsx::addWorksheet(wb_intensities_zscores, sheetname)
# Add Z-scores and create plots
if (z_score == 1) {
- dir.create(paste0(outdir, "/plots"), showWarnings = FALSE)
+ dir.create("./plots", showWarnings = FALSE)
wb_helix_zscores <- openxlsx::createWorkbook("SinglePatient")
openxlsx::addWorksheet(wb_helix_zscores, sheetname)
row_helix <- 2 # start on row 2 because of header
@@ -90,16 +85,15 @@ if (z_score == 1) {
outlist[, intensity_col_ids][outlist[, intensity_col_ids] == 0] <- NA
# calculate robust Z-scores
- outlist_robust_zscore <- calculate_zscores(outlist, "_RobustZscore", control_col_idx, perc, intensity_col_ids, startcol)
+ outlist_robust_zscore <- calculate_zscores(outlist, "_RobustZscore", control_col_idx, perc, intensity_col_ids)
# calculate Z-scores after removal of outliers in Control samples with grubbs test
outlist_nooutliers <- calculate_zscores(
- outlist, "_OutlierRemovedZscore", control_col_idx, outlier_threshold,
- intensity_col_ids, startcol
+ outlist, "_OutlierRemovedZscore", control_col_idx, outlier_threshold, intensity_col_ids
)
# calculate Z-scores
- outlist <- calculate_zscores(outlist, "_Zscore", control_intensities, NULL, intensity_col_ids, startcol)
+ outlist <- calculate_zscores(outlist, "_Zscore", control_col_idx, NULL, intensity_col_ids)
# output metabolites filtered on relevance
save_to_rdata_and_txt(outlist, "AdductSums_filtered_Zscores")
@@ -213,7 +207,7 @@ if (z_score == 1) {
plots_present = TRUE
)
openxlsx::writeData(wb_helix_intensities, sheet = 1, outlist_helix, startCol = 1)
- openxlsx::saveWorkbook(wb_helix_intensities, paste0(outdir, "/Helix_", project, ".xlsx"), overwrite = TRUE)
+ openxlsx::saveWorkbook(wb_helix_intensities, paste0("Helix_", project, ".xlsx"), overwrite = TRUE)
rm(wb_helix_intensities)
# reorder outlist for Excel file
@@ -237,6 +231,6 @@ if (z_score == 1) {
# write Excel file
openxlsx::writeData(wb_intensities_zscores, sheet = 1, outlist, startCol = 1)
-openxlsx::saveWorkbook(wb_intensities_zscores, paste0(outdir, "/", project, ".xlsx"), overwrite = TRUE)
+openxlsx::saveWorkbook(wb_intensities_zscores, paste0(project, ".xlsx"), overwrite = TRUE)
rm(wb_intensities_zscores)
unlink("plots", recursive = TRUE)
diff --git a/DIMS/GenerateExcel.nf b/DIMS/GenerateExcel.nf
index 552a8ee9..a4badbe6 100644
--- a/DIMS/GenerateExcel.nf
+++ b/DIMS/GenerateExcel.nf
@@ -1,7 +1,7 @@
process GenerateExcel {
tag "DIMS GenerateExcel"
label 'GenerateExcel'
- container = 'docker://umcugenbioinf/dims:1.3'
+ container = 'ghcr.io/umcugenetics/dims:v1.4.0'
shell = ['/bin/bash', '-euo', 'pipefail']
input:
diff --git a/DIMS/GenerateQCOutput.R b/DIMS/GenerateQCOutput.R
index cb3d865f..5d52d6ee 100644
--- a/DIMS/GenerateQCOutput.R
+++ b/DIMS/GenerateQCOutput.R
@@ -3,9 +3,6 @@ library("reshape2")
library("openxlsx")
suppressMessages(library("dplyr"))
-# set the number of digits for floats
-options(digits = 16)
-
# define parameters
cmd_args <- commandArgs(trailingOnly = TRUE)
@@ -15,11 +12,17 @@ dims_matrix <- cmd_args[3]
sst_components_file <- cmd_args[4]
export_scripts_dir <- cmd_args[5]
-outdir <- "./"
-
# load in function scripts
source(paste0(export_scripts_dir, "generate_qc_output_functions.R"))
+# Initialize
+options(digits = 16)
+control_label <- "C"
+# get current date
+rundate <- Sys.Date()
+# create a directory for plots
+dir.create("plots", showWarnings = FALSE)
+
# load init files
load(init_file)
# load outlist from GenerateExcel
@@ -28,22 +31,14 @@ load("outlist.RData")
load("AdductSums_positive.RData")
load("AdductSums_negative.RData")
-# get current date
-rundate <- Sys.Date()
-
-# create a directory for plots in project directory
-dir.create(paste0(outdir, "/plots"), showWarnings = FALSE)
-
-control_label <- "C"
-
-#### CHECK NUMBER OF CONTROLS ####
+# Check number of controls
if (any(grepl("nr_ctrls", colnames(outlist)))) {
file_name <- "Check_number_of_controls.txt"
min_num_controls <- 25
check_number_of_controls(outlist, min_num_controls, file_name)
}
-#### INTERNAL STANDARDS ####
+# INTERNAL STANDARDS
is_list <- outlist[grep("Internal standard", outlist[, "relevance"], fixed = TRUE), ]
is_codes <- rownames(is_list)
@@ -53,7 +48,7 @@ if (length(sample_names_nodata) == 0) {
sample_names_nodata <- "none"
}
write.table(sample_names_nodata,
- file = paste(outdir, "sample_names_nodata.txt", sep = "/"),
+ file = "sample_names_nodata.txt",
row.names = FALSE, col.names = FALSE, quote = FALSE
)
if (!is.null(sample_names_nodata)) {
@@ -75,7 +70,7 @@ is_pos <- get_internal_standards(is_list, "pos", outlist_tot_pos, dims_matrix, r
is_neg <- get_internal_standards(is_list, "neg", outlist_tot_neg, dims_matrix, rundate, project)
# Save results
-save(is_pos, is_neg, is_summed, file = paste0(outdir, "/", project, "_IS_results.RData"))
+save(is_pos, is_neg, is_summed, file = paste0(project, "_IS_results.RData"))
# number of samples, for plotting length and width
sample_count <- length(repl_pattern)
@@ -92,15 +87,15 @@ plot_width <- 9 + 0.35 * sample_count
plot_height <- plot_width / 2.5
save_internal_standard_plot(
- is_neg, "barplot", "Interne Standaard (Neg)", outdir,
+ is_neg, "barplot", "Interne Standaard (Neg)",
"IS_bar_all_neg", plot_width, plot_height
)
save_internal_standard_plot(
- is_pos, "barplot", "Interne Standaard (Pos)", outdir,
+ is_pos, "barplot", "Interne Standaard (Pos)",
"IS_bar_all_pos", plot_width, plot_height
)
save_internal_standard_plot(
- is_summed, "barplot", "Interne Standaard (Summed)", outdir,
+ is_summed, "barplot", "Interne Standaard (Summed)",
"IS_bar_all_sum", plot_width, plot_height
)
@@ -110,15 +105,15 @@ plot_height <- plot_width / 2.5
save_internal_standard_plot(
is_neg, "lineplot", "Interne Standaard (Neg)",
- outdir, "IS_line_all_neg", plot_width, plot_height
+ "IS_line_all_neg", plot_width, plot_height
)
save_internal_standard_plot(
is_pos, "lineplot", "Interne Standaard (Pos)",
- outdir, "IS_line_all_pos", plot_width, plot_height
+ "IS_line_all_pos", plot_width, plot_height
)
save_internal_standard_plot(
is_summed, "lineplot", "Interne Standaard (Sum)",
- outdir, "IS_line_all_sum", plot_width, plot_height
+ "IS_line_all_sum", plot_width, plot_height
)
## bar plots with a selection of IS
@@ -202,7 +197,7 @@ if (dims_matrix == "Plasma") {
} else if (dims_matrix == "DBS") {
is_below_threshold_neg <- find_is_below_threshold(is_neg_selection_subset, threshold_is_dbs_neg, is_neg_selection, "neg")
is_below_threshold_pos <- find_is_below_threshold(is_pos_selection_subset, threshold_is_dbs_pos, is_pos_selection, "pos")
- is_below_threshold_sum <- find_is_below_threshold(is_sum_selection_subset, threshold_is_dbs_sum, is_neg_selection, "sum")
+ is_below_threshold_sum <- find_is_below_threshold(is_sum_selection_subset, threshold_is_dbs_sum, is_sum_selection, "sum")
is_below_threshold <- rbind(is_below_threshold_pos, is_below_threshold_neg, is_below_threshold_sum)
} else {
# generate empty table
@@ -223,50 +218,50 @@ if (nrow(is_below_threshold) > 0) {
# bar plot either with or without minimal intensity lines
if (add_min_intens_lines) {
save_internal_standard_plot(
- is_neg_selection_subset, "barplot", "Interne Standaard (Neg)", outdir,
+ is_neg_selection_subset, "barplot", "Interne Standaard (Neg)",
"IS_bar_select_neg", plot_width, plot_height, hline_data_neg
)
save_internal_standard_plot(
- is_pos_selection_subset, "barplot", "Interne Standaard (Pos)", outdir,
+ is_pos_selection_subset, "barplot", "Interne Standaard (Pos)",
"IS_bar_select_pos", plot_width, plot_height, hline_data_pos
)
save_internal_standard_plot(
- is_sum_selection_subset, "barplot", "Interne Standaard (Sum)", outdir,
+ is_sum_selection_subset, "barplot", "Interne Standaard (Sum)",
"IS_bar_select_sum", plot_width, plot_height, hline_data_sum
)
} else {
save_internal_standard_plot(
- is_neg_selection_subset, "barplot", "Interne Standaard (Neg)", outdir,
+ is_neg_selection_subset, "barplot", "Interne Standaard (Neg)",
"IS_bar_select_neg", plot_width, plot_height
)
save_internal_standard_plot(
- is_pos_selection_subset, "barplot", "Interne Standaard (Pos)", outdir,
+ is_pos_selection_subset, "barplot", "Interne Standaard (Pos)",
"IS_bar_select_pos", plot_width, plot_height
)
save_internal_standard_plot(
- is_sum_selection_subset, "barplot", "Interne Standaard (Sum)", outdir,
+ is_sum_selection_subset, "barplot", "Interne Standaard (Sum)",
"IS_bar_select_sum", plot_width, plot_height
)
}
-## line plots with a selection of IS
+# line plots with a selection of IS
plot_width <- 8 + 0.2 * sample_count
plot_height <- plot_width / 2.0
save_internal_standard_plot(
- is_neg_selection_subset, "lineplot", "Interne Standaard (Neg)", outdir,
+ is_neg_selection_subset, "lineplot", "Interne Standaard (Neg)",
"IS_line_select_neg", plot_width, plot_height
)
save_internal_standard_plot(
- is_pos_selection_subset, "lineplot", "Interne Standaard (Pos)", outdir,
+ is_pos_selection_subset, "lineplot", "Interne Standaard (Pos)",
"IS_line_select_pos", plot_width, plot_height
)
save_internal_standard_plot(
- is_sum_selection_subset, "lineplot", "Interne Standaard (Sum)", outdir,
+ is_sum_selection_subset, "lineplot", "Interne Standaard (Sum)",
"IS_line_select_sum", plot_width, plot_height
)
-### POSITIVE CONTROLS CHECK
+# POSITIVE CONTROLS
# these positive controls need to be in the samplesheet, in order to make the positive_control.RData file
# Positive control samples all have the format P1002.x, P1003.x and P1005.x (where x is a number)
@@ -293,7 +288,7 @@ if (length(pos_contr_warning) == 0) {
pos_contr_warning <- "No positive controls found"
}
write.table(pos_contr_warning,
- file = paste(outdir, "positive_controls_warning.txt", sep = "/"),
+ file = "positive_controls_warning.txt",
row.names = FALSE, col.names = FALSE, quote = FALSE
)
@@ -332,17 +327,15 @@ if (length(positive_control_list) > 0) {
positive_control$Project <- project
# Save results
- save(positive_control, file = paste0(outdir, "/", project, "_positive_control.RData"))
+ save(positive_control, file = paste0(project, "_positive_control.RData"))
# round the Z-scores to 2 digits
positive_control$Zscore <- round_df(positive_control$Zscore, 2)
write.xlsx(positive_control,
- file = paste0(outdir, "/", project, "_positive_control.xlsx"),
+ file = paste0(project, "_positive_control.xlsx"),
sheetName = "Sheet1", col.names = TRUE, row.names = TRUE, append = FALSE
)
}
-### SST components output ####
-
# Internal standards lists, calculate coefficients of variation
if ("plots" %in% colnames(is_list)) {
intensity_col_ids <- 2:(which(colnames(is_list) == "HMDB_name") - 1)
@@ -354,7 +347,7 @@ is_list_intensities <- get_is_intensities(is_list, int_cols = intensity_col_ids)
is_neg_intensities <- get_is_intensities(outlist_tot_neg, is_codes = is_codes)
is_pos_intensities <- get_is_intensities(outlist_tot_pos, is_codes = is_codes)
-# SST components
+# SST COMPONENTS
sst_components <- read.csv(sst_components_file, header = TRUE, sep = "\t")
sst_metabolites_df <- outlist %>% filter(HMDB_code %in% sst_components$HMDB_ID)
sst_sample_column_index <- grep("P1001", colnames(sst_metabolites_df))
@@ -401,7 +394,7 @@ setColWidths(wb, 3, cols = 1, widths = 24)
addWorksheet(wb, "SST components")
openxlsx::writeData(wb, sheet = 4, sst_intensities_df)
setColWidths(wb, 4, cols = 1:3, widths = 24)
-xlsx_name <- paste0(outdir, "/", project, "_IS_SST.xlsx")
+xlsx_name <- paste0(project, "_IS_SST.xlsx")
openxlsx::saveWorkbook(wb, xlsx_name, overwrite = TRUE)
rm(wb)
@@ -410,19 +403,19 @@ if (sum(grepl("P1001", colnames(sst_intensities_df))) > 0) {
zscore_column <- grep("_Zscore", colnames(sst_intensities_df))[1]
sst_intensities_df_qc <- sst_intensities_df[sst_intensities_df[, zscore_column] < 2, ]
sst_intensities_df_qc <- select(sst_intensities_df_qc, -c("CV_controls"))
- write.table(sst_intensities_df_qc, file = paste(outdir, "sst_qc.txt", sep = "/"), row.names = FALSE, sep = "\t")
+ write.table(sst_intensities_df_qc, file = "sst_qc.txt", row.names = FALSE, sep = "\t")
} else {
- write.table("no SST sample present", file = paste(outdir, "sst_qc.txt", sep = "/"), row.names = FALSE, col.names = FALSE)
+ write.table("no SST sample present", file = "sst_qc.txt", row.names = FALSE, col.names = FALSE)
}
-### MISSING M/Z CHECK
-# check the outlist_identified_(negative/positive).RData files for missing m/z values and save to file
-# Load the outlist_identified files + remove the loaded files
-load(paste0(outdir, "/outlist_identified_negative.RData"))
-mzmed_pgrp_ident_neg <- outlist_ident$mzmed.pgrp
-load(paste0(outdir, "/outlist_identified_positive.RData"))
-mzmed_pgrp_ident_pos <- outlist_ident$mzmed.pgrp
-rm(outlist_ident)
+# MISSING M/Z CHECK
+# check the peakgroup_list_identified_(negative/positive).RData files for missing m/z values and save to file
+# Load the peakgroup_list_identified files + remove the loaded files
+load("./peakgroup_list_identified_negative.RData")
+mzmed_pgrp_ident_neg <- peakgroup_list_ident$mzmed.pgrp
+load("./peakgroup_list_identified_positive.RData")
+mzmed_pgrp_ident_pos <- peakgroup_list_ident$mzmed.pgrp
+rm(peakgroup_list_ident)
# Check for missing mz values, if present returned with vector of missing mz values
mz_missing_neg <- check_missing_mz(mzmed_pgrp_ident_neg, "Negative")
@@ -430,6 +423,6 @@ mz_missing_pos <- check_missing_mz(mzmed_pgrp_ident_pos, "Positive")
# Write both scanmodes to missing_mz_warning file
lapply(c(mz_missing_neg, mz_missing_pos), write,
- file = paste0(outdir, "/missing_mz_warning.txt"),
+ file = "./missing_mz_warning.txt",
append = TRUE, ncolumns = 1000
)
diff --git a/DIMS/GenerateQCOutput.nf b/DIMS/GenerateQCOutput.nf
index b39d0584..0c9244fc 100644
--- a/DIMS/GenerateQCOutput.nf
+++ b/DIMS/GenerateQCOutput.nf
@@ -1,7 +1,7 @@
process GenerateQCOutput {
tag "DIMS GenerateQCOutput"
label 'GenerateQCOutput'
- container = 'docker://umcugenbioinf/dims:1.3'
+ container = 'ghcr.io/umcugenetics/dims:v1.4.0'
shell = ['/bin/bash', '-euo', 'pipefail']
input:
diff --git a/DIMS/GenerateViolinPlots.R b/DIMS/GenerateViolinPlots.R
index 9d1f3e70..d57910b6 100644
--- a/DIMS/GenerateViolinPlots.R
+++ b/DIMS/GenerateViolinPlots.R
@@ -1,12 +1,11 @@
-# load packages
+# load required packages
suppressPackageStartupMessages(library("dplyr"))
+suppressPackageStartupMessages(library("gridExtra"))
library(reshape2)
library(openxlsx)
library(ggplot2)
-suppressPackageStartupMessages(library("gridExtra"))
library(stringr)
-options(digits = 16)
# define parameters
cmd_args <- commandArgs(trailingOnly = TRUE)
@@ -18,11 +17,31 @@ file_ratios_metabolites <- cmd_args[4]
file_expected_biomarkers_iem <- cmd_args[5]
file_explanation <- cmd_args[6]
+# Initialize
+options(digits = 16)
+iem_variables <- list(
+ top_number_iem_diseases = 5,
+ threshold_iem = 5
+)
+# number of diseases that score highest in dIEM algorithm to plot
+top_number_iem_diseases <- 5
+# probability score cut-off for plotting the top diseases
+threshold_iem <- 5
+nr_plots_perpage <- 20
+zscore_cutoff <- 5
+protocol_name <- "DIMS_PL_DIAG"
+number_of_metabolites <- list(
+ highest = 20,
+ lowest = 10
+)
+
# load functions
source(paste0(export_scripts_dir, "generate_violin_plots_functions.R"))
+
# load dataframe with intensities and Z-scores for all samples
intensities_zscore_df <- get(load("outlist.RData"))
rm(outlist)
+
# read input files
metabolites_ratios_df <- read.csv(file_ratios_metabolites, sep = ";", stringsAsFactors = FALSE)
expected_biomarkers_df <- read.csv(file_expected_biomarkers_iem, sep = ";", stringsAsFactors = FALSE)
@@ -33,21 +52,6 @@ expected_biomarkers_df <- expected_biomarkers_df %>%
)
explanation_violin_plot <- readLines(file_explanation)
-# Set global variables
-iem_variables <- list(
- top_number_iem_diseases = 5,
- threshold_iem = 5
-)
-top_number_iem_diseases <- 5 # number of diseases that score highest in algorithm to plot
-threshold_iem <- 5 # probability score cut-off for plotting the top diseases
-nr_plots_perpage <- 20 # number of violin plots per page in PDF
-zscore_cutoff <- 5
-protocol_name <- "DIMS_PL_DIAG"
-number_of_metabolites <- list(
- highest = 20,
- lowest = 10
-)
-
control_ids <- get_colnames_by_prefix(intensities_zscore_df, "C")
patient_ids <- get_colnames_by_prefix(intensities_zscore_df, "P")
all_sample_ids <- c(control_ids, patient_ids)
@@ -69,7 +73,7 @@ zscore_controls_df <- intensities_zscore_ratios_df %>%
select(HMDB_code, HMDB_name, any_of(paste0(control_ids, "_Zscore"))) %>%
rename_with(~ str_remove(.x, "_Zscore"), .cols = contains("_Zscore"))
-#### Make violin plots #####
+# Make violin plots
make_and_save_violin_plot_pdfs(
zscore_patients_df,
zscore_controls_df,
@@ -82,12 +86,12 @@ make_and_save_violin_plot_pdfs(
number_of_metabolites
)
-#### Run the IEM algorithm #########
+# Run the IEM algorithm
diem_probability_score <- run_diem_algorithm(expected_biomarkers_df, zscore_patients_df, patient_ids)
save_prob_scores_to_excel(diem_probability_score, run_name)
-#### Generate dIEM plots #########
+# Generate dIEM plots
patient_no_iem <- make_and_save_diem_plots(
diem_probability_score,
patient_ids,
diff --git a/DIMS/GenerateViolinPlots.nf b/DIMS/GenerateViolinPlots.nf
index ec65a2e9..4dbe55f4 100755
--- a/DIMS/GenerateViolinPlots.nf
+++ b/DIMS/GenerateViolinPlots.nf
@@ -1,7 +1,7 @@
process GenerateViolinPlots {
tag "DIMS GenerateViolinPlots"
label 'GenerateViolinPlots'
- container = 'docker://umcugenbioinf/dims:1.3'
+ container = 'ghcr.io/umcugenetics/dims:v1.4.0'
shell = ['/bin/bash', '-euo', 'pipefail']
input:
diff --git a/DIMS/HMDBparts.R b/DIMS/HMDBparts.R
index 8c2234f4..346e0a58 100755
--- a/DIMS/HMDBparts.R
+++ b/DIMS/HMDBparts.R
@@ -5,8 +5,9 @@ db_file <- cmd_args[1]
breaks_file <- cmd_args[2]
standard_run <- cmd_args[3]
-# load file with binning breaks
+# load breaks_file: contains breaks_fwhm & breaks_fwhm_avg
load(breaks_file)
+# determine minimum and maximum m/z values
min_mz <- round(breaks_fwhm[1])
max_mz <- round(breaks_fwhm[length(breaks_fwhm)])
diff --git a/DIMS/HMDBparts.nf b/DIMS/HMDBparts.nf
index 760b28d0..f254e8f9 100644
--- a/DIMS/HMDBparts.nf
+++ b/DIMS/HMDBparts.nf
@@ -1,7 +1,7 @@
process HMDBparts {
tag "DIMS HMDBparts"
label 'HMDBparts'
- container = 'docker://umcugenbioinf/dims:1.3'
+ container = 'ghcr.io/umcugenetics/dims:v1.4.0'
shell = ['/bin/bash', '-euo', 'pipefail']
input:
diff --git a/DIMS/HMDBparts_main.R b/DIMS/HMDBparts_main.R
index 1a377eb6..43b1f131 100644
--- a/DIMS/HMDBparts_main.R
+++ b/DIMS/HMDBparts_main.R
@@ -4,7 +4,12 @@ cmd_args <- commandArgs(trailingOnly = TRUE)
db_file <- cmd_args[1]
breaks_file <- cmd_args[2]
+# Initialize
+options(digits = 16)
+
+# load HMDB database file
load(db_file)
+# load breaks_file: contains breaks_fwhm & breaks_fwhm_avg
load(breaks_file)
# get minimum and maximum m/z in dataset
diff --git a/DIMS/HMDBparts_main.nf b/DIMS/HMDBparts_main.nf
index b38bac08..0b49dfa9 100644
--- a/DIMS/HMDBparts_main.nf
+++ b/DIMS/HMDBparts_main.nf
@@ -1,7 +1,7 @@
process HMDBparts_main {
tag "DIMS HMDBparts_main"
label 'HMDBparts_main'
- container = 'docker://umcugenbioinf/dims:1.3'
+ container = 'ghcr.io/umcugenetics/dims:v1.4.0'
shell = ['/bin/bash', '-euo', 'pipefail']
input:
diff --git a/DIMS/MakeInit.nf b/DIMS/MakeInit.nf
index 7aae0e46..51932633 100644
--- a/DIMS/MakeInit.nf
+++ b/DIMS/MakeInit.nf
@@ -1,7 +1,7 @@
process MakeInit {
tag "DIMS MakeInit"
label 'MakeInit'
- container = 'docker://umcugenbioinf/dims:1.3'
+ container = 'ghcr.io/umcugenetics/dims:v1.4.0'
shell = ['/bin/bash', '-euo', 'pipefail']
input:
diff --git a/DIMS/PeakFinding.R b/DIMS/PeakFinding.R
index 74f5d31e..f3500b49 100644
--- a/DIMS/PeakFinding.R
+++ b/DIMS/PeakFinding.R
@@ -1,3 +1,4 @@
+# loead required packages
library(dplyr)
# define parameters
@@ -6,12 +7,15 @@ cmd_args <- commandArgs(trailingOnly = TRUE)
replicate_rdatafile <- cmd_args[1]
resol <- as.numeric(cmd_args[2])
preprocessing_scripts_dir <- cmd_args[3]
-# use fixed theshold between noise and signal for peak
-peak_thresh <- 2000
# source functions script
source(paste0(preprocessing_scripts_dir, "peak_finding_functions.R"))
+# Initialize
+options(digits = 16)
+# theshold between noise and signal for peak
+peak_thresh <- 2000
+
# Load output of AssignToBins (peak_list) for a technical replicate
load(replicate_rdatafile)
techrepl_name <- colnames(peak_list$pos)[1]
@@ -19,9 +23,6 @@ techrepl_name <- colnames(peak_list$pos)[1]
# load list of technical replicates per sample that passed threshold filter
techreps_passed <- read.table("replicates_per_sample.txt", sep=",")
-# Initialize
-options(digits = 16)
-
# do peak finding
scanmodes <- c("positive", "negative")
for (scanmode in scanmodes) {
diff --git a/DIMS/PeakFinding.nf b/DIMS/PeakFinding.nf
index 1d02e505..a3565239 100644
--- a/DIMS/PeakFinding.nf
+++ b/DIMS/PeakFinding.nf
@@ -1,7 +1,7 @@
process PeakFinding {
tag "DIMS PeakFinding ${rdata_file}"
label 'PeakFinding'
- container = 'docker://umcugenbioinf/dims:1.3'
+ container = 'ghcr.io/umcugenetics/dims:v1.4.0'
shell = ['/bin/bash', '-euo', 'pipefail']
input:
diff --git a/DIMS/PeakGrouping.R b/DIMS/PeakGrouping.R
index e0a5c828..ee44ee31 100644
--- a/DIMS/PeakGrouping.R
+++ b/DIMS/PeakGrouping.R
@@ -1,9 +1,9 @@
-# define parameters
-cmd_args <- commandArgs(trailingOnly = TRUE)
-
# load required packages
library("dplyr")
+# define parameters
+cmd_args <- commandArgs(trailingOnly = TRUE)
+
hmdb_part_file <- cmd_args[1]
preprocessing_scripts_dir <- cmd_args[2]
ppm <- as.numeric(cmd_args[3])
@@ -11,6 +11,7 @@ ppm <- as.numeric(cmd_args[3])
# load in function scripts
source(paste0(preprocessing_scripts_dir, "peak_grouping_functions.R"))
+# Initialize
options(digits = 16)
# load part of the HMDB
@@ -42,8 +43,6 @@ outlist_df$height.pkt <- as.numeric(outlist_df$height.pkt)
rm(outlist_total)
sample_names <- unique(outlist_df$samplenr)
-## peak grouping
-peakgrouplist <- NULL
# limit the peaklist to the m/z range in the HMDB part, with ppm tolerance
minmz_hmdbpart <- min(hmdb_add_iso[, column_label])
maxmz_hmdbpart <- max(hmdb_add_iso[, column_label])
@@ -54,10 +53,10 @@ outlist_mzrange <- outlist_df[outlist_df$mzmed.pkt > (minmz_hmdbpart - mz_tolera
outlist_sorted <- outlist_mzrange %>% dplyr::arrange(desc(height.pkt))
# find peak groups
-ints_sorted <- find_peak_groups(outlist_sorted, mz_tolerance, sample_names)
+peakgrouplist <- find_peak_groups(outlist_sorted, mz_tolerance, sample_names)
# do annotation
-peakgrouplist_identified <- annotate_peak_groups(ints_sorted, hmdb_add_iso, column_label, mz_tolerance)
+peakgrouplist_identified <- annotate_peak_groups(peakgrouplist, hmdb_add_iso, column_label, mz_tolerance)
# write output to file
save(peakgrouplist_identified, file = paste0(batch_number, "_", scanmode, "_identified.RData"))
diff --git a/DIMS/PeakGrouping.nf b/DIMS/PeakGrouping.nf
index 6cc4ddb0..a69bf006 100644
--- a/DIMS/PeakGrouping.nf
+++ b/DIMS/PeakGrouping.nf
@@ -1,7 +1,7 @@
process PeakGrouping {
tag "DIMS PeakGrouping ${hmdbpart_file}"
label 'PeakGrouping'
- container = 'docker://umcugenbioinf/dims:1.3'
+ container = 'ghcr.io/umcugenetics/dims:v1.4.0'
shell = ['/bin/bash', '-euo', 'pipefail']
input:
diff --git a/DIMS/SumAdducts.R b/DIMS/SumAdducts.R
index 489ae179..c0dbfd6a 100644
--- a/DIMS/SumAdducts.R
+++ b/DIMS/SumAdducts.R
@@ -1,4 +1,6 @@
+# load required libraries
library("dplyr")
+
# define parameters
cmd_args <- commandArgs(trailingOnly = TRUE)
@@ -20,7 +22,7 @@ if (grepl("positive_hmdb", hmdbpart_main_file)) {
}
# load input files
-collect_file <- paste0("outlist_identified_", scanmode, ".RData")
+collect_file <- paste0("peakgroup_list_identified_", scanmode, ".RData")
peakgroup_list <- get(load(collect_file))
hmdb_main_part <- get(load(hmdbpart_main_file))
diff --git a/DIMS/SumAdducts.nf b/DIMS/SumAdducts.nf
index 4b3f9650..f089307d 100644
--- a/DIMS/SumAdducts.nf
+++ b/DIMS/SumAdducts.nf
@@ -1,7 +1,7 @@
process SumAdducts {
tag "DIMS SumAdducts ${hmdbpart_main_file}"
label 'SumAdducts'
- container = 'docker://umcugenbioinf/dims:1.3'
+ container = 'ghcr.io/umcugenetics/dims:v1.4.0'
shell = ['/bin/bash', '-euo', 'pipefail']
input:
diff --git a/DIMS/Utils/RawFiles.nf b/DIMS/Utils/RawFiles.nf
deleted file mode 100644
index 8cc0efee..00000000
--- a/DIMS/Utils/RawFiles.nf
+++ /dev/null
@@ -1,12 +0,0 @@
-def extractRawfilesFromDir(dir) {
- // Original code from: https://github.com/SciLifeLab/Sarek - MIT License - Copyright (c) 2016 SciLifeLab
- dir = dir.tokenize().collect{"$it/*.raw"}
- Channel
- .fromPath(dir, type:'file')
- .ifEmpty { error "No raw files found in ${dir}." }
- .map { rawfiles_path ->
- def file_id = rawfiles_path.getSimpleName()
- [file_id, rawfiles_path]
- }
-}
-
diff --git a/DIMS/Utils/add_lab_id_and_onderzoeksnummer.R b/DIMS/Utils/add_lab_id_and_onderzoeksnummer.R
deleted file mode 100644
index f2ff13ed..00000000
--- a/DIMS/Utils/add_lab_id_and_onderzoeksnummer.R
+++ /dev/null
@@ -1,16 +0,0 @@
-add_lab_id_and_onderzoeksnummer <- function(df_metabs_helix) {
- #' Adding labnummer and Onderzoeksnummer to a dataframe
- #'
- #' @param df_metabs_helix: dataframe with patient data to be uploaded to Helix
- #'
- #' @return: dataframe with added labnummer and Onderzoeksnummer columns
-
- # Split patient number into labnummer and Onderzoeksnummer
- for (row in 1:nrow(df_metabs_helix)) {
- df_metabs_helix[row, "labnummer"] <- gsub("^P|\\.[0-9]*", "", df_metabs_helix[row, "Patient"])
- labnummer_split <- strsplit(as.character(df_metabs_helix[row, "labnummer"]), "M")[[1]]
- df_metabs_helix[row, "Onderzoeksnummer"] <- paste0("MB", labnummer_split[1], "/", labnummer_split[2])
- }
-
- return(df_metabs_helix)
-}
diff --git a/DIMS/Utils/atomic_info.R b/DIMS/Utils/atomic_info.R
deleted file mode 100644
index 94d0e83a..00000000
--- a/DIMS/Utils/atomic_info.R
+++ /dev/null
@@ -1,92 +0,0 @@
-## adapted from globalAssignments.HPC.R
-# relative abundancies from theoretical mass and composition
-# not a function, but a large amount of objects that will become available in memory
-# refactor: include in get_element_info
-# check if all these elements are necessary. Only use snake_make for Hmass -> hydrogen_mass.
-# move source library to higher level
-options(digits = 16)
-
-suppressPackageStartupMessages(library(lattice))
-
-# The following list was copied from Rdisop elements.R and corrected for C, H, O, Cl, S according to NIST
-Ba <- list(name = "Ba", mass = 130, isotope = list(mass = c(-0.093718, 0, -0.094958, 0, -0.095514, -0.094335, -0.095447, -0.094188, -0.094768), abundance = c(0.00106, 0, 0.00101, 0, 0.02417, 0.06592, 0.07854, 0.1123, 0.717)))
-Br <- list(name = "Br", mass = 79, isotope = list(mass = c(-0.0816639, 0, -0.083711), abundance = c(0.5069, 0, 0.4931)))
-C <- list(name = "C", mass = 12, isotope = list(mass = c(0, 0.003354838, 0.003241989), abundance = c(0.9893, 0.0107, 0)))
-Ca <- list(name = "Ca", mass = 40, isotope = list(mass = c(-0.0374094, 0, -0.0413824, -0.0412338, -0.0445194, 0, -0.046311, 0, -0.047467), abundance = c(0.96941, 0, 0.00647, 0.00135, 0.02086, 0, 4e-05, 0, 0.00187)))
-Cl <- list(name = "Cl", mass = 35, isotope = list(mass = c(-0.03114732, 0, -0.03409741), abundance = c(0.7576, 0, 0.2424)))
-Cr <- list(name = "Cr", mass = 50, isotope = list(mass = c(-0.0539536, 0, -0.0594902, -0.0593487, -0.0611175), abundance = c(0.04345, 0, 0.83789, 0.09501, 0.02365)))
-Cu <- list(name = "Cu", mass = 63, isotope = list(mass = c(-0.0704011, 0, -0.0722071), abundance = c(0.6917, 0, 0.3083)))
-F <- list(name = "F", mass = 19, isotope = list(mass = c(-0.00159678), abundance = c(1)))
-Fe <- list(name = "Fe", mass = 54, isotope = list(mass = c(-0.0603873, 0, -0.0650607, -0.0646042, -0.0667227), abundance = c(0.058, 0, 0.9172, 0.022, 0.0028)))
-H <- list(name = "H", mass = 1, isotope = list(mass = c(0.00782503207, 0.014101778, 0.01604928), abundance = c(0.999885, 0.000115, 0)))
-Hg <- list(name = "Hg", mass = 196, isotope = list(mass = c(-0.034193, 0, -0.033257, -0.031746, -0.0317, -0.029723, -0.029383, 0, -0.026533), abundance = c(0.0015, 0, 0.0997, 0.1687, 0.231, 0.1318, 0.2986, 0, 0.0687)))
-I <- list(name = "I", mass = 127, isotope = list(mass = c(-0.095527), abundance = c(1)))
-K <- list(name = "K", mass = 39, isotope = list(mass = c(-0.0362926, -0.0360008, -0.0381746), abundance = c(0.932581, 0.000117, 0.067302)))
-Li <- list(name = "Li", mass = 6, isotope = list(mass = c(0.0151214, 0.016003), abundance = c(0.075, 0.925)))
-Mg <- list(name = "Mg", mass = 24, isotope = list(mass = c(-0.0149577, -0.0141626, -0.0174063), abundance = c(0.7899, 0.1, 0.1101)))
-Mn <- list(name = "Mn", mass = 55, isotope = list(mass = c(-0.0619529), abundance = c(1)))
-N <- list(name = "N", mass = 14, isotope = list(mass = c(0.003074002, 0.00010897), abundance = c(0.99634, 0.00366)))
-Na <- list(name = "Na", mass = 23, isotope = list(mass = c(-0.0102323), abundance = c(1)))
-Ni <- list(name = "Ni", mass = 58, isotope = list(mass = c(-0.0646538, 0, -0.0692116, -0.0689421, -0.0716539, 0, -0.0720321), abundance = c(0.68077, 0, 0.26223, 0.0114, 0.03634, 0, 0.00926)))
-O <- list(name = "O", mass = 16, isotope = list(mass = c(-0.00508538044, -0.0008683, -0.0008397), abundance = c(0.99757, 0.000381, 0.00205)))
-P <- list(name = "P", mass = 31, isotope = list(mass = c(-0.026238), abundance = c(1)))
-S <- list(name = "S", mass = 32, isotope = list(mass = c(-0.027929, -0.02854124, -0.0321331, 0, -0.03291924), abundance = c(0.9499, 0.0075, 0.0425, 0, 1e-04)))
-Se <- list(name = "Se", mass = 74, isotope = list(mass = c(-0.0775254, 0, -0.080788, -0.0800875, -0.0826924, 0, -0.0834804, 0, -0.0833022), abundance = c(0.0089, 0, 0.0936, 0.0763, 0.2378, 0, 0.4961, 0, 0.0873)))
-Si <- list(name = "Si", mass = 28, isotope = list(mass = c(-0.0230729, -0.0235051, -0.0262293), abundance = c(0.9223, 0.0467, 0.031)))
-Sn <- list(name = "Sn", mass = 112, isotope = list(mass = c(-0.095174, 0, -0.097216, -0.096652, -0.098253, -0.097044, -0.098391, -0.09669, -0.0978009, 0, -0.0965596, 0, -0.0947257), abundance = c(0.0097, 0, 0.0065, 0.0034, 0.1453, 0.0768, 0.2423, 0.0859, 0.3259, 0, 0.0463, 0, 0.0579)))
-Zn <- list(name = "Zn", mass = 64, isotope = list(mass = c(-0.0708552, 0, -0.0739653, -0.0728709, -0.0751541, 0, -0.074675), abundance = c(0.486, 0, 0.279, 0.041, 0.188, 0, 0.006)))
-
-# The following list is for our own use
-NH4 <- list(name = "NH4", mass = 18, isotope = list(mass = c(0.03437, 0.03141, -0.95935)), abundance = c(0.995, 0.004, 0.001)) # SISweb: 18.03437 100 19.03141 0.4 19.04065 0.1
-Ac <- list(name = "Ac", mass = 60, isotope = list(mass = c(0.02114, 0.02450, 0.02538)), abundance = c(0.975, 0.021, 0.004)) # SISweb: 60.02114 100 61.02450 2.2 62.02538 0.4
-NaCl <- list(name = "NaCl", mass = 58, isotope = list(mass = c(-0.04137, 0, -0.04433)), abundance = c(0.755, 0, 0.245)) # SISweb: 57.95862 100 59.95567 32.4
-NaCl2 <- list(name = "NaCl2", mass = 116, isotope = list(mass = c(-0.08274, 0, -0.08866)), abundance = c(0.755, 0, 0.245)) # SISweb: 57.95862 100 59.95567 32.4
-NaCl3 <- list(name = "NaCl3", mass = 174, isotope = list(mass = c(-0.12411, 0, -0.13299)), abundance = c(0.755, 0, 0.245)) # SISweb: 57.95862 100 59.95567 32.4
-NaCl4 <- list(name = "NaCl4", mass = 232, isotope = list(mass = c(-0.16548, 0, -0.17732)), abundance = c(0.755, 0, 0.245)) # SISweb: 57.95862 100 59.95567 32.4
-NaCl5 <- list(name = "NaCl5", mass = 290, isotope = list(mass = c(-0.20685, 0, -0.22165)), abundance = c(0.755, 0, 0.245)) # SISweb: 57.95862 100 59.95567 32.4
-For <- list(name = "For", mass = 45, isotope = list(mass = c(-0.00233, 0.00103)), abundance = c(0.989, 0.011)) # SISweb: 46.00549 100 47.00885 1.1 (47.0097 0.1) 48.00973 0.4
-Na2 <- list(name = "2Na-H", mass = 46, isotope = list(mass = c(-1.0282896)), abundance = c(1)) # SISweb for Na2: 45.97954 100 # minus 1 H !
-Met <- list(name = "CH3OH", mass = 32, isotope = list(mass = c(1.034045, 1.037405)), abundance = c(0.989, 0.011)) # SISweb: 32.02622 100 33.02958 1.1 33.0325 0.1 34.03046 0.2
-CH3OH <- list(name = "CH3OH", mass = 32, isotope = list(mass = c(1.034045, 1.037405)), abundance = c(0.989, 0.011)) # SISweb: 32.02622 100 33.02958 1.1 33.0325 0.1 34.03046 0.2
-Na3 <- list(name = "3Na-2H", mass = 69, isotope = list(mass = c(-2.0463469)), abundance = c(1)) # SISweb for Na2: 45.97954 100 # minus 1 H !
-KCl <- list(name = "KCl", mass = 74, isotope = list(mass = c(-0.06744, 0.92961, -0.06744, 0.92772)), abundance = c(0.7047, 0.2283, 0.0507, 0.0162)) # SISweb: 73.93256 100 75.92961 32.4 75.93067 7.2 77.92772 2.3
-H2PO4 <- list(name = "H2PO4", mass = 97, isotope = list(mass = c(-0.03091)), abundance = c(1))
-HSO4 <- list(name = "HSO4", mass = 97, isotope = list(mass = c(-0.04042, 0, -0.04462)), abundance = c(0.96, 0, 0.04))
-Met2 <- list(name = "Met2", mass = 64, isotope = list(mass = c(1.060265, 1.013405)), abundance = c(0.978, 0.022))
-Met3 <- list(name = "Met3", mass = 96, isotope = list(mass = c(1.086485, 1.089845)), abundance = c(0.969, 0.031))
-Met4 <- list(name = "Met4", mass = 128, isotope = list(mass = c(1.112705, 1.116065)), abundance = c(0.959, 0.041))
-Met5 <- list(name = "Met5", mass = 160, isotope = list(mass = c(1.20935, 1.142285)), abundance = c(0.949, 0.051))
-NaminH <- list(name = "Na-H", mass = 21, isotope = list(mass = c(-0.02571416)), abundance = c(1))
-KminH <- list(name = "K-H", mass = 37, isotope = list(mass = c(-0.05194, 0.94617)), abundance = c(0.9328, 0.0672))
-H2O <- list(name = "H2O", mass = -19, isotope = list(mass = c(-0.01894358)), abundance = c(1))
-NaK <- list(name = "NaK-H", mass = 61, isotope = list(mass = c(-0.054345, 0.943765)), abundance = c(0.9328, 0.0672))
-min2H <- list(name = "min2H", mass = -2, isotope = list(mass = c(-0.0151014)), abundance = c(1))
-plus2H <- list(name = "plus2H", mass = 2, isotope = list(mass = c(0.0151014)), abundance = c(1))
-plus2Na <- list(name = "plus2Na", mass = 46, isotope = list(mass = c(-0.02046), abundance = c(1)))
-plusNaH <- list(name = "plusNaH", mass = 24, isotope = list(mass = c(-0.00295588), abundance = c(1)))
-plusKH <- list(name = "plusKH", mass = 40, isotope = list(mass = c(-0.029008, 0.969101)), abundance = c(0.9328, 0.0672))
-plusHNH <- list(name = "plusHNH", mass = 19, isotope = list(mass = c(0.04164642)), abundance = c(1))
-min3H <- list(name = "min3H", mass = -3, isotope = list(mass = c(-0.02182926)), abundance = c(1))
-plus3H <- list(name = "plus3H", mass = 3, isotope = list(mass = c(0.02182926)), abundance = c(1))
-plus3Na <- list(name = "plus3Na", mass = 68, isotope = list(mass = c(0.96931)), abundance = c(1))
-plus2NaH <- list(name = "plus2NaH", mass = 47, isotope = list(mass = c(0.2985712)), abundance = c(1))
-plusNa2H <- list(name = "plusNa2H", mass = 25, isotope = list(mass = c(0.00432284)), abundance = c(1))
-
-# add first isotope of Cl because of high abundance
-Cl37 <- list(name = "Cl37", mass = 37, isotope = list(mass = c(-0.03409741)), abundance = c(1))
-
-all_elements <- list(Ba, Br, C, Ca, Cl, Cr, Cu, F, Fe, H, Hg, I, K, Li, Mg, Mn, N, Na, Ni, O, P, S, Se, Si, Sn, Zn)
-all_adducts <- list(Ba, Br, Ca, Cl, Cl37, Cr, Cu, Fe, Hg, I, K, Li, Mg, Mn, Na, Ni, Se, Si, Sn, Zn,
- NH4, Ac, NaCl, For, Na2, CH3OH, NaCl2, NaCl3, NaCl4, NaCl5, Na3, KCl, H2PO4, HSO4,
- Met2, Met3, Met4, Met5, NaminH, KminH, H2O, NaK, min2H, plus2H, plus2Na, plusNaH,
- plusKH, min3H, plus3H, plusHNH, plus3Na, plus2NaH, plusNa2H)
-
-atoms_inuse <- c("P", "O", "N", "C", "H", "S", "Cl", "D", "13C", "34S", "18O", "37Cl")
-atomic_weights <- c(30.97376163, 15.99491463, 14.0030740052, 12.0000, 1.0078250321, 31.9720707, 34.968852721, 2.0141017778, 13.0033548378, 33.96786690, 17.9991610, 36.96590259)
-electron <- 0.00054858
-
-hydrogen_mass <- Hmass <- H$mass + H$isotope$mass[1]
-Dmass <- H$mass + 1 + H$isotope$mass[2]
-Tmass <- H$mass + 2 + H$isotope$mass[3]
-C13mass <- C$mass + 1 + C$isotope$mass[2]
-N15mass <- N$mass + 1 + N$isotope$mass[2]
diff --git a/DIMS/Utils/calculate_zscores.R b/DIMS/Utils/calculate_zscores.R
deleted file mode 100644
index 12365e41..00000000
--- a/DIMS/Utils/calculate_zscores.R
+++ /dev/null
@@ -1,63 +0,0 @@
-## adapted from statistics_z.R
-# refactor: change column names from avg.ctrls to avg_ctrls, sd.ctrls to sd_ctrls
-# check logic of parameter adducts
-calculate_zscores <- function(peakgroup_list, adducts) {
- #' Calculate Z-scores for peak groups based on average and standard deviation of controls
- #'
- #' @param peakgroup_list: Peak group list (matrix)
- #' @param sort_col: Column to sort on (string)
- #' @param adducts: Parameter indicating whether there are adducts in the list (boolean)
- #'
- #' @return peakgroup_list_dedup: de-duplicated peak group list (matrix)
-
- case_label <- "P"
- control_label <- "C"
- # get index for new column names
- startcol <- ncol(peakgroup_list) + 3
-
- # calculate mean and standard deviation for Control group
- ctrl_cols <- grep(control_label, colnames(peakgroup_list), fixed = TRUE)
- case_cols <- grep(case_label, colnames(peakgroup_list), fixed = TRUE)
- int_cols <- c(ctrl_cols, case_cols)
- # set al zeros to NA
- peakgroup_list[, int_cols][peakgroup_list[, int_cols] == 0] <- NA
- ctrl_ints <- peakgroup_list[, ctrl_cols, drop = FALSE]
- peakgroup_list$avg.ctrls <- apply(ctrl_ints, 1, function(x) mean(as.numeric(x), na.rm = TRUE))
- peakgroup_list$sd.ctrls <- apply(ctrl_ints, 1, function(x) sd(as.numeric(x), na.rm = TRUE))
-
- # set new column names and calculate Z-scores
- colnames_zscores <- NULL
- for (col_index in int_cols) {
- col_name <- colnames(peakgroup_list)[col_index]
- colnames_zscores <- c(colnames_zscores, paste0(col_name, "_Zscore"))
- zscores_1col <- (as.numeric(as.vector(unlist(peakgroup_list[, col_index]))) -
- peakgroup_list$avg.ctrls) / peakgroup_list$sd.ctrls
- peakgroup_list <- cbind(peakgroup_list, zscores_1col)
- }
-
- # apply new column names to columns at end plus avg and sd columns
- colnames(peakgroup_list)[startcol:ncol(peakgroup_list)] <- colnames_zscores
-
- # add ppm deviation column
- zscore_cols <- grep("Zscore", colnames(peakgroup_list), fixed = TRUE)
- if (!adducts) {
- if ((dim(peakgroup_list[, zscore_cols])[2] + 6) != (startcol - 1)) {
- ppmdev <- array(1:nrow(peakgroup_list), dim = c(nrow(peakgroup_list)))
- # calculate ppm deviation
- for (i in 1:nrow(peakgroup_list)) {
- if (!is.na(peakgroup_list$theormz_HMDB[i]) &&
- !is.null(peakgroup_list$theormz_HMDB[i]) &&
- (peakgroup_list$theormz_HMDB[i] != "")) {
- ppmdev[i] <- 10^6 * (as.numeric(as.vector(peakgroup_list$mzmed.pgrp[i])) -
- as.numeric(as.vector(peakgroup_list$theormz_HMDB[i]))) /
- as.numeric(as.vector(peakgroup_list$theormz_HMDB[i]))
- } else {
- ppmdev[i] <- NA
- }
- }
- peakgroup_list <- cbind(peakgroup_list[, 1:4], ppmdev = ppmdev, peakgroup_list[, 5:ncol(peakgroup_list)])
- }
- }
-
- return(peakgroup_list)
-}
diff --git a/DIMS/Utils/check_overlap.R b/DIMS/Utils/check_overlap.R
deleted file mode 100644
index 66492e86..00000000
--- a/DIMS/Utils/check_overlap.R
+++ /dev/null
@@ -1,24 +0,0 @@
-## adapted from checkOverlap.R
-check_overlap <- function(range1, range2) {
- #' Modify range1 and range2 in case of overlap
- #'
- #' @param range1: Vector of m/z values for first peak (float)
- #' @param range2: Vector of m/z values for second peak (float)
- #'
- #' @return new_ranges: list of two ranges (list)
-
- # Check for overlap
- if (length(intersect(range1, range2)) == 2) {
- if (length(range1) >= length(range2)) {
- range1 <- range1[-length(range1)]
- } else {
- range2 <- range2[-1]
- }
- } else if (length(intersect(range1, range2)) == 3) {
- range1 <- range1[-length(range1)]
- range2 <- range2[-1]
- }
- new_ranges <- list("range1" = range1, "range2" = range2)
- return(new_ranges)
-}
-
diff --git a/DIMS/Utils/check_same_samplename.R b/DIMS/Utils/check_same_samplename.R
deleted file mode 100644
index 6e808410..00000000
--- a/DIMS/Utils/check_same_samplename.R
+++ /dev/null
@@ -1,10 +0,0 @@
-check_same_samplename <- function(int_col_name, zscore_col_name) {
- #' A check to see if intensity and Z-score columns match
- #'
- #' @param int_col_name: name of an intensity column (string)
- #' @param zscore_col_name: name of a Z-score column (string)
- #'
- #' @return: match or mismatch of the columns (boolean)
-
- paste0(int_col_name, "_Zscore") == zscore_col_name
-}
diff --git a/DIMS/Utils/create_violin_plots.R b/DIMS/Utils/create_violin_plots.R
deleted file mode 100644
index be4f1f6d..00000000
--- a/DIMS/Utils/create_violin_plots.R
+++ /dev/null
@@ -1,122 +0,0 @@
-# remove parameter default value
-# add explanation variable to parameters
-create_violin_plots <- function(pdf_dir, pt_name, metab_perpage, top_metab_pt = NULL) {
- #' Create violin plots for each patient
- #'
- #' @param pdf_dir: location where to save the pdf file (string)
- #' @param pt_name: patient code (string)
- #' @param metab_perpage: list of dataframe with a dataframe for each page in the violinplot pdf (list)
- #' @param top_metab_pt: dataframe with increased and decrease metabolites for this patient (dataframe)
-
- # set parameters for plots
- plot_height <- 9.6
- plot_width <- 6
- fontsize <- 1
- circlesize <- 0.8
- colors_4plot <- c("#22E4AC", "#00B0F0", "#504FFF", "#A704FD", "#F36265", "#DA0641")
- # green blue blue/purple purple orange red
-
- # patient plots, create the PDF device
- pt_name_sub <- pt_name
- suffix <- ""
- if (grepl("Diagnostics", pdf_dir) & is_diagnostic_patient(pt_name)) {
- prefix <- "MB"
- suffix <- "_DIMS_PL_DIAG"
- # substitute P and M in P2020M00001 into right format for Helix
- pt_name_sub <- gsub("[PM]", "", pt_name)
- pt_name_sub <- gsub("\\..*", "", pt_name_sub)
- } else if (grepl("Diagnostics", pdf_dir)) {
- prefix <- "Dx_"
- } else if (grepl("IEM", pdf_dir)) {
- prefix <- "IEM_"
- } else {
- prefix <- "R_"
- }
-
- pdf(paste0(pdf_dir, "/", prefix, pt_name_sub, suffix, ".pdf"),
- onefile = TRUE,
- width = plot_width,
- height = plot_height)
-
- # page headers:
- page_headers <- names(metab_perpage)
-
- # put table into PDF file, if not empty
- if (!is.null(dim(top_metab_pt))) {
- plot.new()
- # get the names and numbers in the table aligned
- table_theme <- ttheme_default(core = list(fg_params = list(hjust = 0, x = 0.05, fontsize = 6)),
- colhead = list(fg_params = list(fontsize = 8, fontface = "bold")))
- grid.table(top_metab_pt, theme = table_theme, rows = NULL)
- # g <- tableGrob(top_metab_pt)
- # grid.draw(g)
- text(x = 0.45, y = 1.02, paste0("Top deviating metabolites for patient: ", pt_name), font = 1, cex = 1)
- }
-
- # violin plots
- for (page_index in 1:length(metab_perpage)) {
- # extract list of metabolites to plot on a page
- metab_list_2plot <- metab_perpage[[page_index]]
- # extract original data for patient of interest (pt_name) before cut-offs
- pt_list_2plot_orig <- metab_list_2plot[which(metab_list_2plot$variable == pt_name), ]
- # cut off Z-scores higher than 20 or lower than -5 (for nicer plots)
- metab_list_2plot$value[metab_list_2plot$value > 20] <- 20
- metab_list_2plot$value[metab_list_2plot$value < -5] <- -5
- # extract data for patient of interest (pt_name)
- pt_list_2plot <- metab_list_2plot[which(metab_list_2plot$variable == pt_name), ]
- # restore original Z-score before cut-off, for showing Z-scores in PDF
- pt_list_2plot$value_orig <- pt_list_2plot_orig$value
- # remove patient of interest (pt_name) from list; violins will be made up of controls and other patients
- metab_list_2plot <- metab_list_2plot[-which(metab_list_2plot$variable == pt_name), ]
- # subtitle per page
- sub_perpage <- gsub("_", " ", page_headers[page_index])
- # for IEM plots, put subtitle on two lines
- sub_perpage <- gsub("probability", "\nprobability", sub_perpage)
- # add size parameter for showing Z-score of patient per metabolite
- z_size <- rep(3, nrow(pt_list_2plot))
- # set size to 0 if row is empty
- z_size[is.na(pt_list_2plot$value)] <- 0
-
- # draw violin plot. shape=22 gives square for patient of interest
- ggplot_object <- ggplot(metab_list_2plot, aes(x = value, y = HMDB_name)) +
- theme(axis.text.y = element_text(size = rel(fontsize)), plot.caption = element_text(size = rel(fontsize))) +
- xlim(-5, 20) +
- geom_violin(scale = "width") +
- geom_point(data = pt_list_2plot, aes(color = value), size = 3.5 * circlesize, shape = 22, fill = "white") +
- scale_fill_gradientn(
- colors = colors_4plot, values = NULL, space = "Lab", na.value = "grey50", guide = "colourbar",
- aesthetics = "colour"
- ) +
- # add Z-score value for patient of interest at x=16
- geom_text(
- data = pt_list_2plot, aes(16, label = paste0("Z=", round(value_orig, 2))), hjust = "left", vjust = +0.2,
- size = z_size
- ) +
- # add labels. Use font Courier to get all the plots in the same location.
- labs(x = "Z-scores", y = "Metabolites", subtitle = sub_perpage, color = "z-score") +
- theme(axis.text.y = element_text(family = "Courier", size = 6)) +
- # do not show legend
- theme(legend.position = "none") +
- # add title
- ggtitle(label = paste0("Results for patient ", pt_name)) +
- # add vertical lines
- geom_vline(xintercept = 2, col = "grey", lwd = 0.5, lty = 2) +
- geom_vline(xintercept = -2, col = "grey", lwd = 0.5, lty = 2)
-
- suppressWarnings(print(ggplot_object))
-
- }
-
- # add explanation of violin plots, version number etc.
- plot(NA, xlim = c(0, 5), ylim = c(0, 5), bty = "n", xaxt = "n", yaxt = "n", xlab = "", ylab = "")
- if (length(explanation) > 0) {
- text(0.2, 5, explanation[1], pos = 4, cex = 0.8)
- for (line_index in 2:length(explanation)) {
- text_y_position <- 5 - (line_index * 0.2)
- text(-0.2, text_y_position, explanation[line_index], pos = 4, cex = 0.5)
- }
- }
-
- # close the PDF file
- dev.off()
-}
diff --git a/DIMS/Utils/do_peakfinding.R b/DIMS/Utils/do_peakfinding.R
deleted file mode 100644
index fff6afe8..00000000
--- a/DIMS/Utils/do_peakfinding.R
+++ /dev/null
@@ -1,55 +0,0 @@
-## adapted from findpeaks.Gauss.HPC.R
-# NB: this function will be taken up into PeakFinding.R
-# variables with fixed values will be removed from function parameters
-# int_factor, scale, outdir, plot, thresh, width, height
-do_peakfinding <- function(sample_avgtechrepl, int_factor, scale, resol, outdir, scanmode, plot, thresh, width, height) {
- #' start peak finding
- #'
- #' @param sample_avgtechrepl: Dataframe with binned intensities averaged over technical replicates for a sample
- #' @param int_factor: Value used to calculate area under Gaussian curve (integer)
- #' @param scale: Initial value used to estimate scaling parameter (integer)
- #' @param resol: Value for resolution (integer)
- #' @param outdir: Path for output directory (string)
- #' @param scanmode: Scan mode, positive or negative (string)
- #' @param plot: Parameter indicating whether plots should be made (boolean)
- #' @param thresh: Value for noise level threshold (integer)
- #' @param width: Value for width of plot (integer)
- #' @param height: Value for height of plot (integer)
- #'
- #' @return save output to file
-
- sample_name <- colnames(sample_avgtechrepl)[1]
-
- # turn dataframe with intensities into a named list
- ints_fullrange <- as.vector(sample_avgtechrepl)
- names(ints_fullrange) <- rownames(sample_avgtechrepl)
-
- # initialise list to store results for all peaks
- allpeaks_values <- list("mean" = NULL, "area" = NULL, "nr" = NULL,
- "min" = NULL, "max" = NULL, "qual" = NULL, "spikes" = 0)
-
- # look for m/z range for all peaks
- allpeaks_values <- search_mzrange(ints_fullrange, allpeaks_values, int_factor, scale, resol,
- outdir, sample_name, scanmode,
- plot, width, height, thresh)
-
- # turn the list into a dataframe
- outlist_persample <- NULL
- outlist_persample <- cbind("samplenr" = allpeaks_values$nr,
- "mzmed.pkt" = allpeaks_values$mean,
- "fq" = allpeaks_values$qual,
- "mzmin.pkt" = allpeaks_values$min,
- "mzmax.pkt" = allpeaks_values$max,
- "height.pkt" = allpeaks_values$area)
-
- # remove peaks with height = 0
- outlist_persample <- outlist_persample[outlist_persample[, "height.pkt"] != 0, ]
-
- # save output to file
- save(outlist_persample, file = paste0(sample_name, "_", scanmode, ".RData"))
-
- # generate text output to log file on number of spikes for this sample
- # spikes are peaks that are too narrow, e.g. 1 data point
- cat(paste("There were", allpeaks_values$spikes, "spikes"))
-}
-
diff --git a/DIMS/Utils/estimate_area.R b/DIMS/Utils/estimate_area.R
deleted file mode 100644
index 9f9042ea..00000000
--- a/DIMS/Utils/estimate_area.R
+++ /dev/null
@@ -1,23 +0,0 @@
-estimate_area <- function(mass_max, resol, scale, sigma, int_factor) {
- #' Estimate area of Gaussian curve
- #'
- #' @param mass_max: Value for m/z at maximum intensity of a peak (float)
- #' @param resol: Value for resolution (integer)
- #' @param scale: Value for peak width (float)
- #' @param sigma: Value for standard deviation (float)
- #' @param int_factor: Value used to calculate area under Gaussian curve (integer)
- #'
- #' @return area_curve: Value for area under the Gaussian curve (float)
-
- # generate a mass_vector with equally spaced m/z values
- fwhm <- get_fwhm(mass_max, resol)
- mz_min <- mass_max - 2 * fwhm
- mz_max <- mass_max + 2 * fwhm
- mz_range <- mz_max - mz_min
- mass_vector2 <- seq(mz_min, mz_max, length = 1000)
-
- # estimate area under the curve
- area_curve <- sum(scale * dnorm(mass_vector2, mass_max, sigma)) / 100
-
- return(area_curve)
-}
diff --git a/DIMS/Utils/fit_gaussian.R b/DIMS/Utils/fit_gaussian.R
deleted file mode 100644
index 6be96600..00000000
--- a/DIMS/Utils/fit_gaussian.R
+++ /dev/null
@@ -1,319 +0,0 @@
-## adapted from fitGaussian.R
-# variables with fixed values will be removed from function parameters
-# scale, outdir, plot, width, height
-# max_index doesn't need to be passed to this function, can be determined here.
-# remove plot sections (commented out)
-# several functions need to be loaded before this function can run
-fit_gaussian <- function(mass_vector2, mass_vector, int_vector, max_index, scale, resol,
- outdir, force, use_bounds, plot, scanmode,
- int_factor, width, height) {
- #' Fit 1, 2, 3 or 4 Gaussian peaks in small region of m/z
- #'
- #' @param mass_vector2: Vector of equally spaced m/z values (float)
- #' @param mass_vector: Vector of m/z values for a region of interest (float)
- #' @param int_vector: Value used to calculate area under Gaussian curve (integer)
- #' @param max_index: Index in int_vector with the highest intensity (integer)
- #' @param scale: Initial value used to estimate scaling parameter (integer)
- #' @param resol: Value for resolution (integer)
- #' @param outdir: Path for output directory (string)
- #' @param force: Number of local maxima in int_vector (integer)
- #' @param use_bounds: Boolean to indicate whether boundaries are to be used
- #' @param plot: Parameter indicating whether plots should be made (boolean)
- #' @param scanmode: Scan mode, positive or negative (string)
- #' @param int_factor: Value used to calculate area under Gaussian curve (integer)
- #' @param width: Value for width of plot (integer)
- #' @param height: Value for height of plot (integer)
- #'
- #' @return roi_value_list: list of fit values for region of interest (list)
-
- # Initialise
- peak_mean <- NULL
- peak_area <- NULL
- peak_qual <- NULL
- peak_min <- NULL
- peak_max <- NULL
- fit_quality1 <- 0.15
- fit_quality <- 0.2
-
- # One local maximum:
- if (force == 1) {
- # determine fit values for 1 Gaussian peak (mean, scale, sigma, qual)
- fit_values <- fit_1peak(mass_vector2, mass_vector, int_vector, max_index, scale, resol,
- plot, fit_quality1, use_bounds)
- # set initial value for scale factor
- scale <- 2
- # test if the mean is outside the m/z range
- if (fit_values$mean[1] < mass_vector[1] || fit_values$mean[1] > mass_vector[length(mass_vector)]) {
- # run this function again with fixed boundaries
- return(fit_gaussian(mass_vector2, mass_vector, int_vector, max_index, scale, resol,
- outdir, force = 1, use_bounds = TRUE, plot, scanmode, int_factor, width, height))
- } else {
- # test if the fit is bad
- if (fit_values$qual > fit_quality1) {
- # Try to fit two curves; find two local maxima
- new_index <- which(diff(sign(diff(int_vector))) == -2) + 1
- # test if there are two indices in new_index
- if (length(new_index) != 2) {
- new_index <- round(length(mass_vector) / 3)
- new_index <- c(new_index, 2 * new_index)
- }
- # run this function again with two local maxima
- return(fit_gaussian(mass_vector2, mass_vector, int_vector, new_index,
- scale, resol, outdir, force = 2, use_bounds = FALSE,
- plot, scanmode, int_factor, width, height))
- # good fit
- } else {
- peak_mean <- c(peak_mean, fit_values$mean)
- peak_area <- c(peak_area, estimate_area(fit_values$mean, resol, fit_values$scale,
- fit_values$sigma, int_factor))
- peak_qual <- fit_values$qual
- peak_min <- mass_vector[1]
- peak_max <- mass_vector[length(mass_vector)]
- }
- }
-
- #### Two local maxima; need at least 6 data points for this ####
- } else if (force == 2 && (length(mass_vector) > 6)) {
- # determine fit values for 2 Gaussian peaks (mean, scale, sigma, qual)
- fit_values <- fit_2peaks(mass_vector2, mass_vector, int_vector, max_index, scale, resol,
- use_bounds, plot, fit_quality, int_factor)
- # test if one of the means is outside the m/z range
- if (fit_values$mean[1] < mass_vector[1] || fit_values$mean[1] > mass_vector[length(mass_vector)] ||
- fit_values$mean[2] < mass_vector[1] || fit_values$mean[2] > mass_vector[length(mass_vector)]) {
- # check if fit quality is bad
- if (fit_values$qual > fit_quality) {
- # run this function again with fixed boundaries
- return(fit_gaussian(mass_vector2, mass_vector, int_vector, max_index, scale, resol,
- outdir, force = 2, use_bounds = TRUE,
- plot, scanmode, int_factor, width, height))
- } else {
- # check which mean is outside range and remove it from the list of means
- # NB: peak_mean and other variables have not been given values from 2-peak fit yet!
- for (i in 1:length(fit_values$mean)){
- if (fit_values$mean[i] < mass_vector[1] || fit_values$mean[i] > mass_vector[length(mass_vector)]) {
- peak_mean <- c(peak_mean, -i)
- peak_area <- c(peak_area, -i)
- } else {
- peak_mean <- c(peak_mean, fit_values$mean[i])
- peak_area <- c(peak_area, fit_values$area[i])
- }
- }
- peak_qual <- fit_values$qual
- peak_min <- mass_vector[1]
- peak_max <- mass_vector[length(mass_vector)]
- }
- # if all means are within range
- } else {
- # check for bad fit
- if (fit_values$qual > fit_quality) {
- # Try to fit three curves; find three local maxima
- new_index <- which(diff(sign(diff(int_vector))) == -2) + 1
- # test if there are three indices in new_index
- if (length(new_index) != 3) {
- new_index <- round(length(mass_vector) / 4)
- new_index <- c(new_index, 2 * new_index, 3 * new_index)
- }
- # run this function again with three local maxima
- return(fit_gaussian(mass_vector2, mass_vector, int_vector, new_index,
- scale, resol, outdir, force = 3, use_bounds = FALSE,
- plot, scanmode, int_factor, width, height))
- # good fit, all means are within m/z range
- } else {
- # check if means are within 3 ppm and sum if so
- tmp <- fit_values$qual
- nr_means_new <- -1
- nr_means <- length(fit_values$mean)
- while (nr_means != nr_means_new) {
- nr_means <- length(fit_values$mean)
- fit_values <- within_ppm(fit_values$mean, fit_values$scale, fit_values$sigma, fit_values$area,
- mass_vector2, mass_vector, ppm = 4, resol, plot)
- nr_means_new <- length(fit_values$mean)
- }
- # restore original quality score
- fit_values$qual <- tmp
-
- for (i in 1:length(fit_values$mean)){
- peak_mean <- c(peak_mean, fit_values$mean[i])
- peak_area <- c(peak_area, fit_values$area[i])
- }
- peak_qual <- fit_values$qual
- peak_min <- mass_vector[1]
- peak_max <- mass_vector[length(mass_vector)]
- }
- }
-
- # Three local maxima; need at least 6 data points for this
- } else if (force == 3 && (length(mass_vector) > 6)) {
- # determine fit values for 3 Gaussian peaks (mean, scale, sigma, qual)
- fit_values <- fit_3peaks(mass_vector2, mass_vector, int_vector, max_index, scale, resol,
- use_bounds, plot, fit_quality, int_factor)
- # test if one of the means is outside the m/z range
- if (fit_values$mean[1] < mass_vector[1] || fit_values$mean[1] > mass_vector[length(mass_vector)] ||
- fit_values$mean[2] < mass_vector[1] || fit_values$mean[2] > mass_vector[length(mass_vector)] ||
- fit_values$mean[3] < mass_vector[1] || fit_values$mean[3] > mass_vector[length(mass_vector)]) {
- # check if fit quality is bad
- if (fit_values$qual > fit_quality) {
- # run this function again with fixed boundaries
- return(fit_gaussian(mass_vector2, mass_vector, int_vector, max_index, scale, resol,
- outdir, force, use_bounds = TRUE,
- plot, scanmode, int_factor, width, height))
- } else {
- # check which mean is outside range and remove it from the list of means
- # NB: peak_mean and other variables have not been given values from 2-peak fit yet!
- for (i in 1:length(fit_values$mean)) {
- if (fit_values$mean[i] < mass_vector[1] || fit_values$mean[i] > mass_vector[length(mass_vector)]) {
- peak_mean <- c(peak_mean, -i)
- peak_area <- c(peak_area, -i)
- } else {
- peak_mean <- c(peak_mean, fit_values$mean[i])
- peak_area <- c(peak_area, fit_values$area[i])
- }
- }
- peak_qual <- fit_values$qual
- peak_min <- mass_vector[1]
- peak_max <- mass_vector[length(mass_vector)]
- }
- # if all means are within range
- } else {
- # check for bad fit
- if (fit_values$qual > fit_quality) {
- # Try to fit four curves; find four local maxima
- new_index <- which(diff(sign(diff(int_vector))) == -2) + 1
- # test if there are four indices in new_index
- if (length(new_index) != 4) {
- new_index <- round(length(mass_vector) / 5)
- new_index <- c(new_index, 2 * new_index, 3 * new_index, 4 * new_index)
- }
- # run this function again with four local maxima
- return(fit_gaussian(mass_vector2, mass_vector, int_vector, new_index, scale, resol,
- outdir, force = 4, use_bounds = FALSE, plot, scanmode,
- int_factor, width, height))
- # good fit, all means are within m/z range
- } else {
- # check if means are within 4 ppm and sum if so
- tmp <- fit_values$qual
- nr_means_new <- -1
- nr_means <- length(fit_values$mean)
- while (nr_means != nr_means_new) {
- nr_means <- length(fit_values$mean)
- fit_values <- within_ppm(fit_values$mean, fit_values$scale, fit_values$sigma, fit_values$area,
- mass_vector2, mass_vector, ppm = 4, resol, plot)
- nr_means_new <- length(fit_values$mean)
- }
- # restore original quality score
- fit_values$qual <- tmp
-
- for (i in 1:length(fit_values$mean)){
- peak_mean <- c(peak_mean, fit_values$mean[i])
- peak_area <- c(peak_area, fit_values$area[i])
- }
- peak_qual <- fit_values$qual
- peak_min <- mass_vector[1]
- peak_max <- mass_vector[length(mass_vector)]
-
- }
- }
-
- #### Four local maxima; need at least 6 data points for this ####
- } else if (force == 4 && (length(mass_vector) > 6)) {
- # determine fit values for 4 Gaussian peaks (mean, scale, sigma, qual)
- fit_values <- fit_4peaks(mass_vector2, mass_vector, int_vector, max_index, scale, resol,
- use_bounds, plot, fit_quality, int_factor)
- # test if one of the means is outside the m/z range
- if (fit_values$mean[1] < mass_vector[1] || fit_values$mean[1] > mass_vector[length(mass_vector)] ||
- fit_values$mean[2] < mass_vector[1] || fit_values$mean[2] > mass_vector[length(mass_vector)] ||
- fit_values$mean[3] < mass_vector[1] || fit_values$mean[3] > mass_vector[length(mass_vector)] ||
- fit_values$mean[4] < mass_vector[1] || fit_values$mean[4] > mass_vector[length(mass_vector)]) {
- # check if quality of fit is bad
- if (fit_values$qual > fit_quality) {
- # run this function again with fixed boundaries
- return(fit_gaussian(mass_vector2, mass_vector, int_vector, max_index, scale, resol,
- outdir, force, use_bounds = TRUE,
- plot, scanmode, int_factor, width, height))
-
- } else {
- # check which mean is outside range and remove it from the list of means
- # NB: peak_mean and other variables have not been given values from 2-peak fit yet!
- for (i in 1:length(fit_values$mean)) {
- if (fit_values$mean[i] < mass_vector[1] | fit_values$mean[i] > mass_vector[length(mass_vector)]) {
- peak_mean <- c(peak_mean, -i)
- peak_area <- c(peak_area, -i)
- } else {
- peak_mean <- c(peak_mean, fit_values$mean[i])
- peak_area <- c(peak_area, fit_values$area[i])
- }
- }
- peak_qual <- fit_values$qual
- peak_min <- mass_vector[1]
- peak_max <- mass_vector[length(mass_vector)]
- }
- # if all means are within range
- } else {
- # check for bad fit
- if (fit_values$qual > fit_quality) {
- # Try to fit 1 curve, force = 5
- return(fit_gaussian(mass_vector2, mass_vector, int_vector, max_index, scale, resol,
- outdir, force = 5, use_bounds = FALSE,
- plot, scanmode, int_factor, width, height))
- # good fit, all means are within m/z range
- } else {
- # check if means are within 4 ppm and sum if so
- tmp <- fit_values$qual
- nr_means_new <- -1
- nr_means <- length(fit_values$mean)
- while (nr_means != nr_means_new) {
- nr_means <- length(fit_values$mean)
- fit_values <- within_ppm(fit_values$mean, fit_values$scale, fit_values$sigma, fit_values$area,
- mass_vector2, mass_vector, ppm = 4, resol, plot)
- nr_means_new <- length(fit_values$mean)
- }
- # restore original quality score
- fit_values$qual <- tmp
-
- for (i in 1:length(fit_values$mean)){
- peak_mean <- c(peak_mean, fit_values$mean[i])
- peak_area <- c(peak_area, fit_values$area[i])
- }
- peak_qual <- fit_values$qual
- peak_min <- mass_vector[1]
- peak_max <- mass_vector[length(mass_vector)]
-
- }
- }
-
- # More than four local maxima; fit 1 peak.
- } else {
- scale <- 2
- fit_quality1 <- 0.40
- use_bounds <- TRUE
- max_index <- which(int_vector == max(int_vector))
- fit_values <- fit_1peak(mass_vector2, mass_vector, int_vector, max_index, scale, resol,
- plot, fit_quality1, use_bounds)
- # check for bad fit
- if (fit_values$qual > fit_quality1) {
- # remove
- if (plot) dev.off()
- # get fit values from fit_optim
- fit_values <- fit_optim(mass_vector, int_vector, resol, plot, scanmode, int_factor, width, height)
- peak_mean <- c(peak_mean, fit_values$mean)
- peak_area <- c(peak_area, fit_values$area)
- peak_min <- fit_values$min
- peak_max <- fit_values$max
- peak_qual <- 0
- } else {
- peak_mean <- c(peak_mean, fit_values$mean)
- peak_area <- c(peak_area, estimate_area(fit_values$mean, resol, fit_values$scale, fit_values$sigma, int_factor))
- peak_qual <- fit_values$qual
- peak_min <- mass_vector[1]
- peak_max <- mass_vector[length(mass_vector)]
- }
- }
-
- # put all values for this region of interest into a list
- roi_value_list <- list("mean" = peak_mean,
- "area" = peak_area,
- "qual" = peak_qual,
- "min" = peak_min,
- "max" = peak_max)
- return(roi_value_list)
-}
diff --git a/DIMS/Utils/fit_gaussians.R b/DIMS/Utils/fit_gaussians.R
deleted file mode 100644
index 05bff860..00000000
--- a/DIMS/Utils/fit_gaussians.R
+++ /dev/null
@@ -1,196 +0,0 @@
-# Gaussian fit functions
-## adapted from fitG.R, fit2G.R, fit3G.R and fit4G.R (combined)
-fit_1gaussian <- function(mass_vector, int_vector, sigma, query_mass, scale, use_bounds) {
- #' Fit a Gaussian curve for a peak with given parameters
- #'
- #' @param mass_vector: Vector of masses (float)
- #' @param int_vector: Vector of intensities (float)
- #' @param sigma: Value for width of the peak (float)
- #' @param query_mass: Value for mass at center of peak (float)
- #' @param scale: Value for scaling intensities (float)
- #' @param use_bounds: Boolean to indicate whether boundaries are to be used
- #'
- #' @return opt_fit: list of parameters and values describing the optimal fit
-
- # define optimization function for optim based on normal distribution
- opt_f <- function(params) {
- d <- params[2] * dnorm(mass_vector, mean = params[1], sd = sigma)
- sum((d - int_vector) ^ 2)
- }
- if (use_bounds) {
- # determine lower and upper boundaries
- lower <- c(mass_vector[1], 0, mass_vector[1], 0)
- upper <- c(mass_vector[length(mass_vector)], Inf, mass_vector[length(mass_vector)], Inf)
- # get optimal value for fitted Gaussian curve
- opt_fit <- optim(c(as.numeric(query_mass), as.numeric(scale)),
- opt_f, control = list(maxit = 10000), method = "L-BFGS-B",
- lower = lower, upper = upper)
- } else {
- opt_fit <- optim(c(as.numeric(query_mass), as.numeric(scale)),
- opt_f, control = list(maxit = 10000))
- }
- return(opt_fit)
-}
-
-
-fit_2gaussians <- function(mass_vector, int_vector, sigma1, sigma2,
- query_mass1, scale1,
- query_mass2, scale2, use_bounds) {
- #' Fit two Gaussian curves for a peak with given parameters
- #'
- #' @param mass_vector: Vector of masses (float)
- #' @param int_vector: Vector of intensities (float)
- #' @param sigma1: Value for width of the first peak (float)
- #' @param sigma2: Value for width of the second peak (float)
- #' @param query_mass1: Value for mass at center of first peak (float)
- #' @param scale1: Value for scaling intensities for first peak (float)
- #' @param query_mass2: Value for mass at center of second peak (float)
- #' @param scale2: Value for scaling intensities for second peak (float)
- #' @param use_bounds: Boolean to indicate whether boundaries are to be used
- #'
- #' @return opt_fit: list of parameters and values describing the optimal fit
-
- # define optimization function for optim based on normal distribution
- opt_f <- function(params) {
- d <- params[2] * dnorm(mass_vector, mean = params[1], sd = sigma1) +
- params[4] * dnorm(mass_vector, mean = params[3], sd = sigma2)
- sum((d - int_vector) ^ 2)
- }
-
- if (use_bounds) {
- # determine lower and upper boundaries
- lower <- c(mass_vector[1], 0, mass_vector[1], 0)
- upper <- c(mass_vector[length(mass_vector)], Inf, mass_vector[length(mass_vector)], Inf)
- # get optimal value for 2 fitted Gaussian curves
- if (is.null(query_mass2) && is.null(scale2) && is.null(sigma2)) {
- sigma2 <- sigma1
- opt_fit <- optim(c(as.numeric(query_mass1), as.numeric(scale1),
- as.numeric(query_mass1), as.numeric(scale1)),
- opt_f, control = list(maxit = 10000),
- method = "L-BFGS-B", lower = lower, upper = upper)
- } else {
- opt_fit <- optim(c(as.numeric(query_mass1), as.numeric(scale1),
- as.numeric(query_mass2), as.numeric(scale2)),
- opt_f, control = list(maxit = 10000),
- method = "L-BFGS-B", lower = lower, upper = upper)
- }
- } else {
- if (is.null(query_mass2) && is.null(scale2) && is.null(sigma2)) {
- sigma2 <- sigma1
- opt_fit <- optim(c(as.numeric(query_mass1), as.numeric(scale1),
- as.numeric(query_mass1), as.numeric(scale1)),
- opt_f, control = list(maxit = 10000))
- } else {
- opt_fit <- optim(c(as.numeric(query_mass1), as.numeric(scale1),
- as.numeric(query_mass2), as.numeric(scale2)),
- opt_f, control = list(maxit = 10000))
- }
- }
- return(opt_fit)
-}
-
-
-fit_3gaussians <- function(mass_vector, int_vector, sigma1, sigma2, sigma3,
- query_mass1, scale1,
- query_mass2, scale2,
- query_mass3, scale3, use_bounds) {
- #' Fit three Gaussian curves for a peak with given parameters
- #'
- #' @param mass_vector: Vector of masses (float)
- #' @param int_vector: Vector of intensities (float)
- #' @param sigma1: Value for width of the first peak (float)
- #' @param sigma2: Value for width of the second peak (float)
- #' @param sigma3: Value for width of the third peak (float)
- #' @param query_mass1: Value for mass at center of first peak (float)
- #' @param scale1: Value for scaling intensities for first peak (float)
- #' @param query_mass2: Value for mass at center of second peak (float)
- #' @param scale2: Value for scaling intensities for second peak (float)
- #' @param query_mass3: Value for mass at center of third peak (float)
- #' @param scale3: Value for scaling intensities for third peak (float)
- #' @param use_bounds: Boolean to indicate whether boundaries are to be used
- #'
- #' @return opt_fit: list of parameters and values describing the optimal fit
-
- # define optimization function for optim based on normal distribution
- opt_f <- function(params) {
- d <- params[2] * dnorm(mass_vector, mean = params[1], sd = sigma1) +
- params[4] * dnorm(mass_vector, mean = params[3], sd = sigma2) +
- params[6] * dnorm(mass_vector, mean = params[5], sd = sigma3)
- sum((d - int_vector) ^ 2)
- }
-
- if (use_bounds) {
- # determine lower and upper boundaries
- lower <- c(mass_vector[1], 0, mass_vector[1], 0, mass_vector[1], 0)
- upper <- c(mass_vector[length(mass_vector)], Inf, mass_vector[length(mass_vector)], Inf,
- mass_vector[length(mass_vector)], Inf)
- # get optimal value for 3 fitted Gaussian curves
- opt_fit <- optim(c(query_mass1, scale1,
- query_mass2, scale2,
- query_mass3, scale3),
- opt_f, control = list(maxit = 10000),
- method = "L-BFGS-B", lower = lower, upper = upper)
- } else {
- opt_fit <- optim(c(query_mass1, scale1,
- query_mass2, scale2,
- query_mass3, scale3),
- opt_f, control = list(maxit = 10000))
- }
- return(opt_fit)
-}
-
-fit_4gaussians <- function(mass_vector, int_vector, sigma1, sigma2, sigma3, sigma4,
- query_mass1, scale1,
- query_mass2, scale2,
- query_mass3, scale3,
- query_mass4, scale4, use_bounds) {
- #' Fit four Gaussian curves for a peak with given parameters
- #'
- #' @param mass_vector: Vector of masses (float)
- #' @param int_vector: Vector of intensities (float)
- #' @param sigma1: Value for width of the first peak (float)
- #' @param sigma2: Value for width of the second peak (float)
- #' @param sigma3: Value for width of the third peak (float)
- #' @param sigma4: Value for width of the fourth peak (float)
- #' @param query_mass1: Value for mass at center of first peak (float)
- #' @param scale1: Value for scaling intensities for first peak (float)
- #' @param query_mass2: Value for mass at center of second peak (float)
- #' @param scale2: Value for scaling intensities for second peak (float)
- #' @param query_mass3: Value for mass at center of third peak (float)
- #' @param scale3: Value for scaling intensities for third peak (float)
- #' @param query_mass4: Value for mass at center of fourth peak (float)
- #' @param scale4: Value for scaling intensities for fourth peak (float)
- #' @param use_bounds: Boolean to indicate whether boundaries are to be used
- #'
- #' @return opt_fit: list of parameters and values describing the optimal fit
-
- # define optimization function for optim based on normal distribution
- opt_f <- function(params) {
- d <- params[2] * dnorm(mass_vector, mean = params[1], sd = sigma1) +
- params[4] * dnorm(mass_vector, mean = params[3], sd = sigma2) +
- params[6] * dnorm(mass_vector, mean = params[5], sd = sigma3) +
- params[8] * dnorm(mass_vector, mean = params[7], sd = sigma4)
- sum((d - int_vector) ^ 2)
- }
-
- if (use_bounds) {
- # determine lower and upper boundaries
- lower <- c(mass_vector[1], 0, mass_vector[1], 0, mass_vector[1], 0, mass_vector[1], 0)
- upper <- c(mass_vector[length(mass_vector)], Inf, mass_vector[length(mass_vector)], Inf,
- mass_vector[length(mass_vector)], Inf, mass_vector[length(mass_vector)], Inf)
- # get optimal value for 3 fitted Gaussian curves
- opt_fit <- optim(c(query_mass1, scale1,
- query_mass2, scale2,
- query_mass3, scale3,
- query_mass4, scale4),
- opt_f, control = list(maxit = 10000),
- method = "L-BFGS-B", lower = lower, upper = upper)
- } else {
- opt_fit <- optim(c(query_mass1, scale1,
- query_mass2, scale2,
- query_mass3, scale3,
- query_mass4, scale4),
- opt_f, control = list(maxit = 10000))
- }
- return(opt_fit)
-}
diff --git a/DIMS/Utils/fit_init.R b/DIMS/Utils/fit_init.R
deleted file mode 100644
index ecae3375..00000000
--- a/DIMS/Utils/fit_init.R
+++ /dev/null
@@ -1,47 +0,0 @@
-## adapted from fitGaussianInit.R
-# variables with fixed values will be removed from function parameters
-# scale, outdir, plot, width, height
-# mz_index, start_index, end_index, sample_name not used.
-# fit_gaussian should be defined before this function is called.
-fit_init <- function(mass_vector, int_vector, int_factor, scale, resol,
- outdir, sample_name, scanmode, plot, width, height,
- mz_index, start_index, end_index) {
- #' Determine initial fit of Gaussian curve to small region of m/z
- #'
- #' @param mass_vector: Vector of m/z values for a region of interest (float)
- #' @param int_vector: Vector of intensities for a region of interest (float)
- #' @param int_factor: Value used to calculate area under Gaussian curve (integer)
- #' @param scale: Initial value used to estimate scaling parameter (integer)
- #' @param resol: Value for resolution (integer)
- #' @param outdir: Path for output directory (string)
- #' @param sample_name: Sample name (string)
- #' @param scanmode: Scan mode, positive or negative (string)
- #' @param plot: Parameter indicating whether plots should be made (boolean)
- #' @param width: Value for width of plot (integer)
- #' @param height: Value for height of plot (integer)
- #' @param mz_index: Index of m/z value with non-zero intensity (integer)
- #' @param start_index: Index of start of m/z range in mass_vector (integer)
- #' @param end_index: Index of end of m/z range in mass_vector (integer)
- #'
- #' @return roi_value_list: list of fit values for region of interest (list)
-
- # define mass_diff as difference between last and first value of mass_vector
- mass_diff <- mass_vector[length(mass_vector)] - mass_vector[1]
- # generate a second mass_vector with equally spaced m/z values
- mass_vector2 <- seq(mass_vector[1], mass_vector[length(mass_vector)],
- length = mass_diff * int_factor)
-
- # Find the index in int_vector with the highest intensity
- max_index <- which(int_vector == max(int_vector))
- roi_values <- fit_gaussian(mass_vector2, mass_vector, int_vector, max_index,
- scale, resol, outdir, force = length(max_index),
- use_bounds = FALSE, plot, scanmode, int_factor, width, height)
- # put all values for this region of interest into a list
- roi_value_list <- list("mean" = roi_values$mean,
- "area" = roi_values$area,
- "qual" = roi_values$qual,
- "min" = roi_values$min,
- "max" = roi_values$max)
- return(roi_value_list)
-}
-
diff --git a/DIMS/Utils/fit_optim.R b/DIMS/Utils/fit_optim.R
deleted file mode 100644
index 404cc7a7..00000000
--- a/DIMS/Utils/fit_optim.R
+++ /dev/null
@@ -1,47 +0,0 @@
-## adapted from generateGaussian.R
-# variables with fixed values will be removed from function parameters
-# plot, width, height
-# fit_gaussian should be defined before this function is called.
-fit_optim <- function(mass_vector, int_vector, resol,
- plot, scanmode, int_factor, width, height) {
- #' Determine optimized fit of Gaussian curve to small region of m/z
- #'
- #' @param mass_vector: Vector of m/z values for a region of interest (float)
- #' @param int_vector: Vector of intensities for a region of interest (float)
- #' @param resol: Value for resolution (integer)
- #' @param plot: Parameter indicating whether plots should be made (boolean)
- #' @param scanmode: Scan mode, positive or negative (string)
- #' @param int_factor: Value used to calculate area under Gaussian curve (integer)
- #' @param width: Value for width of plot (integer)
- #' @param height: Value for height of plot (integer)
- #'
- #' @return roi_value_list: list of fit values for region of interest (list)
-
- factor <- 1.5
- # Find the index in int_vector with the highest intensity
- max_index <- which(int_vector == max(int_vector))[1]
- mass_max <- mass_vector[max_index]
- int_max <- int_vector[max_index]
- # get peak width
- fwhm <- get_fwhm(mass_max, resol)
- # simplify the peak shape: represent it by a triangle
- mass_max_simple <- c(mass_max - factor * fwhm, mass_max, mass_max + factor * fwhm)
- int_max_simple <- c(0, int_max, 0)
-
- # define mass_diff as difference between last and first value of mass_max_simple
- mass_diff <- mass_max_simple[length(mass_max_simple)] - mass_max_simple[1]
- # generate a second mass_vector with equally spaced m/z values
- mass_vector2 <- seq(mass_max_simple[1], mass_max_simple[length(mass_max_simple)],
- length = mass_diff * int_factor)
- sigma <- get_stdev(mass_vector2, int_max_simple)
- scale <- optimize_gaussfit(mass_vector2, int_max_simple, sigma, mass_max)
-
- # get an estimate of the area under the peak
- area <- estimate_area(mass_max, resol, scale, sigma, int_factor)
- # put all values for this region of interest into a list
- roi_value_list <- list("mean" = mass_max,
- "area" = area,
- "min" = mass_vector2[1],
- "max" = mass_vector2[length(mass_vector2)])
- return(roi_value_list)
-}
diff --git a/DIMS/Utils/fit_peaks.R b/DIMS/Utils/fit_peaks.R
deleted file mode 100644
index 02d54086..00000000
--- a/DIMS/Utils/fit_peaks.R
+++ /dev/null
@@ -1,381 +0,0 @@
-## adapted from fit1Peak.R, fit2peaks.R, fit3peaks.R and fit4peaks.R (combined)
-# variables with fixed values will be removed from function parameters
-# plot, int_factor
-fit_1peak <- function(mass_vector2, mass_vector, int_vector, max_index, scale, resol, plot, fit_quality, use_bounds) {
- #' Fit 1 Gaussian peak in small region of m/z
- #'
- #' @param mass_vector2: Vector of equally spaced m/z values (float)
- #' @param mass_vector: Vector of m/z values for a region of interest (float)
- #' @param int_vector: Value used to calculate area under Gaussian curve (integer)
- #' @param max_index: Index in int_vector with the highest intensity (integer)
- #' @param scale: Initial value used to estimate scaling parameter (integer)
- #' @param resol: Value for resolution (integer)
- #' @param plot: Parameter indicating whether plots should be made (boolean)
- #' @param fit_quality: Value indicating quality of fit of Gaussian curve (float)
- #' @param use_bounds: Boolean to indicate whether boundaries are to be used
- #'
- #' @return roi_value_list: list of fit values for region of interest (list)
-
- if (length(int_vector) < 3) {
- message("Range too small, no fit possible!")
- } else {
- if ((length(int_vector) == 4)) {
- # fit 1 peak
- mu <- weighted.mean(mass_vector, int_vector)
- sigma <- get_stdev(mass_vector, int_vector)
- fitted_peak <- fit_1gaussian(mass_vector, int_vector, sigma, mu, scale, use_bounds)
- } else {
- # set range vector
- if ((length(mass_vector) - length(max_index)) < 2) {
- range1 <- c((length(mass_vector) - 4) : length(mass_vector))
- } else if (length(max_index) < 2) {
- range1 <- c(1:5)
- } else {
- range1 <- c(max_index[1] - 2, max_index[1] - 1, max_index[1], max_index[1] + 1, max_index[1] + 2)
- }
- if (range1[1] == 0) range1 <- range1[-1]
- # remove NA
- if (length(which(is.na(int_vector[range1]))) != 0) {
- range1 <- range1[-which(is.na(int_vector[range1]))]
- }
- # fit 1 peak
- mu <- weighted.mean(mass_vector[range1], int_vector[range1])
- sigma <- get_stdev(mass_vector[range1], int_vector[range1])
- fitted_peak <- fit_1gaussian(mass_vector, int_vector, sigma, mu, scale, use_bounds)
- }
-
- p1 <- fitted_peak$par
-
- # get new value for fit quality and scale
- fq_new <- get_fit_quality(mass_vector, int_vector, p1[1], p1[1], resol, p1[2], sigma)$fq_new
- scale_new <- 1.2 * scale
-
- # bad fit
- if (fq_new > fit_quality) {
- # optimize scaling factor
- fq <- 0
- scale <- 0
- if (sum(int_vector) > sum(p1[2] * dnorm(mass_vector, p1[1], sigma))) {
- while ((round(fq, digits = 3) != round(fq_new, digits = 3)) && (scale_new < 10000)) {
- fq <- fq_new
- scale <- scale_new
- # fit 1 peak
- fitted_peak <- fit_1gaussian(mass_vector, int_vector, sigma, mu, scale, use_bounds)
- p1 <- fitted_peak$par
- # get new value for fit quality and scale
- fq_new <- get_fit_quality(mass_vector, int_vector, p1[1], p1[1], resol, p1[2], sigma)$fq_new
- scale_new <- 1.2 * scale
- }
- } else {
- while ((round(fq, digits = 3) != round(fq_new, digits = 3)) && (scale_new < 10000)) {
- fq <- fq_new
- scale <- scale_new
- # fit 1 peak
- fitted_peak <- fit_1gaussian(mass_vector, int_vector, sigma, mu, scale, use_bounds)
- p1 <- fitted_peak$par
- # get new value for fit quality and scale
- fq_new <- get_fit_quality(mass_vector, int_vector, p1[1], p1[1], resol, p1[2], sigma)$fq_new
- scale_new <- 0.8 * scale
- }
- }
- # use optimized scale factor to fit 1 peak
- if (fq < fq_new) {
- fitted_peak <- fit_1gaussian(mass_vector, int_vector, sigma, mu, scale, use_bounds)
- p1 <- fitted_peak$par
- fq_new <- fq
- }
- }
- }
-
- roi_value_list <- list("mean" = p1[1], "scale" = p1[2], "sigma" = sigma, "qual" = fq_new)
- return(roi_value_list)
-}
-
-fit_2peaks <- function(mass_vector2, mass_vector, int_vector, max_index, scale, resol, use_bounds = FALSE,
- plot = FALSE, fit_quality, int_factor) {
- #' Fit 2 Gaussian peaks in small region of m/z
- #'
- #' @param mass_vector2: Vector of equally spaced m/z values (float)
- #' @param mass_vector: Vector of m/z values for a region of interest (float)
- #' @param int_vector: Value used to calculate area under Gaussian curve (integer)
- #' @param max_index: Index in int_vector with the highest intensity (integer)
- #' @param scale: Initial value used to estimate scaling parameter (integer)
- #' @param resol: Value for resolution (integer)
- #' @param plot: Parameter indicating whether plots should be made (boolean)
- #' @param fit_quality: Value indicating quality of fit of Gaussian curve (float)
- #' @param use_bounds: Boolean to indicate whether boundaries are to be used
- #' @param int_factor: Value used to calculate area under Gaussian curve (integer)
- #'
- #' @return roi_value_list: list of fit values for region of interest (list)
-
- peak_mean <- NULL
- peak_area <- NULL
- peak_scale <- NULL
- peak_sigma <- NULL
-
- # set range vectors for 2 peaks
- range1 <- c(max_index[1] - 2, max_index[1] - 1, max_index[1], max_index[1] + 1, max_index[1] + 2)
- if (range1[1] == 0) range1 <- range1[-1]
- range2 <- c(max_index[2] - 2, max_index[2] - 1, max_index[2], max_index[2] + 1, max_index[2] + 2)
- if (length(mass_vector) < range2[length(range2)]) range2 <- range2[-length(range2)]
- range1 <- check_overlap(range1, range2)[[1]]
- range2 <- check_overlap(range1, range2)[[2]]
- # check for negative or 0
- remove <- which(range1 < 1)
- if (length(remove) > 0) range1 <- range1[-remove]
- remove <- which(range2 < 1)
- if (length(remove) > 0) range2 <- range2[-remove]
- # remove NA
- if (length(which(is.na(int_vector[range1]))) != 0) range1 <- range1[-which(is.na(int_vector[range1]))]
- if (length(which(is.na(int_vector[range2]))) != 0) range2 <- range2[-which(is.na(int_vector[range2]))]
-
- # fit 2 peaks, first separately, then together
- mu1 <- weighted.mean(mass_vector[range1], int_vector[range1])
- sigma1 <- get_stdev(mass_vector[range1], int_vector[range1])
- fitted_peak <- fit_1gaussian(mass_vector[range1], int_vector[range1], sigma1, mu1, scale, use_bounds)
- p1 <- fitted_peak$par
- # second peak
- mu2 <- weighted.mean(mass_vector[range2], int_vector[range2])
- sigma2 <- get_stdev(mass_vector[range2], int_vector[range2])
- fitted_peak <- fit_1gaussian(mass_vector[range2], int_vector[range2], sigma2, mu2, scale, use_bounds)
- p2 <- fitted_peak$par
- # combined
- fitted_2peaks <- fit_2gaussians(mass_vector, int_vector, sigma1, sigma2, p1[1], p1[2], p2[1], p2[2], use_bounds)
- pc <- fitted_2peaks$par
-
- # get fit quality
- if (is.null(sigma2)) sigma2 <- sigma1
- sum_fit <- (pc[2] * dnorm(mass_vector, pc[1], sigma1)) +
- (pc[4] * dnorm(mass_vector, pc[3], sigma2))
- fq <- get_fit_quality(mass_vector, int_vector, sort(c(pc[1], pc[3]))[1], sort(c(pc[1], pc[3]))[2],
- resol, sum_fit = sum_fit)$fq_new
-
- # get parameter values
- area1 <- estimate_area(pc[1], resol, pc[2], sigma1, int_factor)
- area2 <- estimate_area(pc[3], resol, pc[4], sigma2, int_factor)
- peak_area <- c(peak_area, area1)
- peak_area <- c(peak_area, area2)
- peak_mean <- c(peak_mean, pc[1])
- peak_mean <- c(peak_mean, pc[3])
- peak_scale <- c(peak_scale, pc[2])
- peak_scale <- c(peak_scale, pc[4])
- peak_sigma <- c(peak_sigma, sigma1)
- peak_sigma <- c(peak_sigma, sigma2)
-
- roi_value_list <- list("mean" = peak_mean, "scale" = peak_scale, "sigma" = peak_sigma, "area" = peak_area, "qual" = fq)
- return(roi_value_list)
-}
-
-fit_3peaks <- function(mass_vector2, mass_vector, int_vector, max_index, scale, resol, use_bounds = FALSE,
- plot = FALSE, fit_quality, int_factor) {
- #' Fit 3 Gaussian peaks in small region of m/z
- #'
- #' @param mass_vector2: Vector of equally spaced m/z values (float)
- #' @param mass_vector: Vector of m/z values for a region of interest (float)
- #' @param int_vector: Value used to calculate area under Gaussian curve (integer)
- #' @param max_index: Index in int_vector with the highest intensity (integer)
- #' @param scale: Initial value used to estimate scaling parameter (integer)
- #' @param resol: Value for resolution (integer)
- #' @param plot: Parameter indicating whether plots should be made (boolean)
- #' @param fit_quality: Value indicating quality of fit of Gaussian curve (float)
- #' @param use_bounds: Boolean to indicate whether boundaries are to be used
- #' @param int_factor: Value used to calculate area under Gaussian curve (integer)
- #'
- #' @return roi_value_list: list of fit values for region of interest (list)
-
- peak_mean <- NULL
- peak_area <- NULL
- peak_scale <- NULL
- peak_sigma <- NULL
-
- # set range vectors for 3 peaks
- range1 <- c(max_index[1] - 2, max_index[1] - 1, max_index[1], max_index[1] + 1, max_index[1] + 2)
- range2 <- c(max_index[2] - 2, max_index[2] - 1, max_index[2], max_index[2] + 1, max_index[2] + 2)
- range3 <- c(max_index[3] - 2, max_index[3] - 1, max_index[3], max_index[3] + 1, max_index[3] + 2)
- remove <- which(range1 < 1)
- if (length(remove) > 0) {
- range1 <- range1[-remove]
- }
- remove <- which(range2 < 1)
- if (length(remove) > 0) {
- range2 <- range2[-remove]
- }
- if (length(mass_vector) < range3[length(range3)]) range3 <- range3[-length(range3)]
- range1 <- check_overlap(range1, range2)[[1]]
- range2 <- check_overlap(range1, range2)[[2]]
- range2 <- check_overlap(range2, range3)[[1]]
- range3 <- check_overlap(range2, range3)[[2]]
- # check for negative or 0
- remove <- which(range1 < 1)
- if (length(remove) > 0) range1 <- range1[-remove]
- remove <- which(range2 < 1)
- if (length(remove) > 0) range2 <- range2[-remove]
- remove <- which(range3 < 1)
- if (length(remove) > 0) range3 <- range3[-remove]
- # remove NA
- if (length(which(is.na(int_vector[range1]))) != 0) range1 <- range1[-which(is.na(int_vector[range1]))]
- if (length(which(is.na(int_vector[range2]))) != 0) range2 <- range2[-which(is.na(int_vector[range2]))]
- if (length(which(is.na(int_vector[range3]))) != 0) range3 <- range3[-which(is.na(int_vector[range3]))]
-
- # fit 3 peaks, first separately, then together
- mu1 <- weighted.mean(mass_vector[range1], int_vector[range1])
- sigma1 <- get_stdev(mass_vector[range1], int_vector[range1])
- fitted_peak <- fit_1gaussian(mass_vector[range1], int_vector[range1], sigma1, mu1, scale, use_bounds)
- p1 <- fitted_peak$par
- # second peak
- mu2 <- weighted.mean(mass_vector[range2], int_vector[range2])
- sigma2 <- get_stdev(mass_vector[range2], int_vector[range2])
- fitted_peak <- fit_1gaussian(mass_vector[range2], int_vector[range2], sigma2, mu2, scale, use_bounds)
- p2 <- fitted_peak$par
- # third peak
- mu3 <- weighted.mean(mass_vector[range3], int_vector[range3])
- sigma3 <- get_stdev(mass_vector[range3], int_vector[range3])
- fitted_peak <- fit_1gaussian(mass_vector[range3], int_vector[range3], sigma3, mu3, scale, use_bounds)
- p3 <- fitted_peak$par
- # combined
- fitted_3peaks <- fit_3gaussians(mass_vector, int_vector, sigma1, sigma2, sigma3,
- p1[1], p1[2], p2[1], p2[2], p3[1], p3[2], use_bounds)
- pc <- fitted_3peaks$par
-
- # get fit quality
- sum_fit = (pc[2] * dnorm(mass_vector, pc[1], sigma1)) +
- (pc[4] * dnorm(mass_vector, pc[3], sigma2)) +
- (pc[6] * dnorm(mass_vector, pc[5], sigma3))
- fq <- get_fit_quality(mass_vector, int_vector, sort(c(pc[1], pc[3], pc[5]))[1], sort(c(pc[1], pc[3], pc[5]))[3],
- resol, sum_fit = sum_fit)$fq_new
-
- # get parameter values
- area1 <- estimate_area(pc[1], resol, pc[2], sigma1, int_factor)
- area2 <- estimate_area(pc[3], resol, pc[4], sigma2, int_factor)
- area3 <- estimate_area(pc[5], resol, pc[6], sigma3, int_factor)
- peak_area <- c(peak_area, area1)
- peak_area <- c(peak_area, area2)
- peak_area <- c(peak_area, area3)
- peak_mean <- c(peak_mean, pc[1])
- peak_mean <- c(peak_mean, pc[3])
- peak_mean <- c(peak_mean, pc[5])
- peak_scale <- c(peak_scale, pc[2])
- peak_scale <- c(peak_scale, pc[4])
- peak_scale <- c(peak_scale, pc[6])
- peak_sigma <- c(peak_sigma, sigma1)
- peak_sigma <- c(peak_sigma, sigma2)
- peak_sigma <- c(peak_sigma, sigma3)
-
- roi_value_list <- list("mean" = peak_mean, "scale" = peak_scale, "sigma" = peak_sigma, "area" = peak_area, "qual" = fq)
- return(roi_value_list)
-}
-
-fit_4peaks <- function(mass_vector2, mass_vector, int_vector, max_index, scale, resol, use_bounds = FALSE,
- plot = FALSE, fit_quality, int_factor) {
- #' Fit 4 Gaussian peaks in small region of m/z
- #'
- #' @param mass_vector2: Vector of equally spaced m/z values (float)
- #' @param mass_vector: Vector of m/z values for a region of interest (float)
- #' @param int_vector: Value used to calculate area under Gaussian curve (integer)
- #' @param max_index: Index in int_vector with the highest intensity (integer)
- #' @param scale: Initial value used to estimate scaling parameter (integer)
- #' @param resol: Value for resolution (integer)
- #' @param plot: Parameter indicating whether plots should be made (boolean)
- #' @param fit_quality: Value indicating quality of fit of Gaussian curve (float)
- #' @param use_bounds: Boolean to indicate whether boundaries are to be used
- #' @param int_factor: Value used to calculate area under Gaussian curve (integer)
- #'
- #' @return roi_value_list: list of fit values for region of interest (list)
-
- peak_mean <- NULL
- peak_area <- NULL
- peak_scale <- NULL
- peak_sigma <- NULL
-
- # set range vectors for 4 peaks
- range1 <- c(max_index[1] - 2, max_index[1] - 1, max_index[1], max_index[1] + 1, max_index[1] + 2)
- range2 <- c(max_index[2] - 2, max_index[2] - 1, max_index[2], max_index[2] + 1, max_index[2] + 2)
- range3 <- c(max_index[3] - 2, max_index[3] - 1, max_index[3], max_index[3] + 1, max_index[3] + 2)
- range4 <- c(max_index[4] - 2, max_index[4] - 1, max_index[4], max_index[4] + 1, max_index[4] + 2)
- if (range1[1] == 0) range1 <- range1[-1]
- if (length(mass_vector) < range4[length(range4)]) range4 <- range4[-length(range4)]
- range1 <- check_overlap(range1, range2)[[1]]
- range2 <- check_overlap(range1, range2)[[2]]
- range2 <- check_overlap(range2, range3)[[1]]
- range3 <- check_overlap(range2, range3)[[2]]
- range3 <- check_overlap(range3, range4)[[1]]
- range4 <- check_overlap(range3, range4)[[2]]
- remove <- which(range4 > length(mass_vector))
- if (length(remove) > 0) {
- range4 <- range4[-remove]
- }
- # check for negative or 0
- remove <- which(range1 < 1)
- if (length(remove) > 0) range1 <- range1[-remove]
- remove <- which(range2 < 1)
- if (length(remove) > 0) range2 <- range2[-remove]
- remove <- which(range3 < 1)
- if (length(remove) > 0) range3 <- range3[-remove]
- remove <- which(range4 < 1)
- if (length(remove) > 0) range4 <- range4[-remove]
- # remove NA
- if (length(which(is.na(int_vector[range1]))) != 0) range1 <- range1[-which(is.na(int_vector[range1]))]
- if (length(which(is.na(int_vector[range2]))) != 0) range2 <- range2[-which(is.na(int_vector[range2]))]
- if (length(which(is.na(int_vector[range3]))) != 0) range3 <- range3[-which(is.na(int_vector[range3]))]
- if (length(which(is.na(int_vector[range4]))) != 0) range4 <- range4[-which(is.na(int_vector[range4]))]
-
- # fit 4 peaks, first separately, then together
- mu1 <- weighted.mean(mass_vector[range1], int_vector[range1])
- sigma1 <- get_stdev(mass_vector[range1], int_vector[range1])
- fitted_peak <- fit_1gaussian(mass_vector[range1], int_vector[range1], sigma1, mu1, scale, use_bounds)
- p1 <- fitted_peak$par
- # second peak
- mu2 <- weighted.mean(mass_vector[range2], int_vector[range2])
- sigma2 <- get_stdev(mass_vector[range2], int_vector[range2])
- fitted_peak <- fit_1gaussian(mass_vector[range2], int_vector[range2], sigma2, mu2, scale, use_bounds)
- p2 <- fitted_peak$par
- # third peak
- mu3 <- weighted.mean(mass_vector[range3], int_vector[range3])
- sigma3 <- get_stdev(mass_vector[range3], int_vector[range3])
- fitted_peak <- fit_1gaussian(mass_vector[range3], int_vector[range3], sigma3, mu3, scale, use_bounds)
- p3 <- fitted_peak$par
- # fourth peak
- mu4 <- weighted.mean(mass_vector[range4], int_vector[range4])
- sigma4 <- get_stdev(mass_vector[range4], int_vector[range4])
- fitted_peak <- fit_1gaussian(mass_vector[range4], int_vector[range4], sigma4, mu4, scale, use_bounds)
- p4 <- fitted_peak$par
- # combined
- fitted_4peaks <- fit_4gaussians(mass_vector, int_vector, sigma1, sigma2, sigma3, sigma3,
- p1[1], p1[2], p2[1], p2[2], p3[1], p3[2], p4[1], p4[2], use_bounds)
- pc <- fitted_4peaks$par
-
- # get fit quality
- sum_fit <- (pc[2] * dnorm(mass_vector, pc[1], sigma1)) +
- (pc[4] * dnorm(mass_vector, pc[3], sigma2)) +
- (pc[6] * dnorm(mass_vector, pc[5], sigma3)) +
- (pc[8] * dnorm(mass_vector, pc[7], sigma3))
- fq <- get_fit_quality(mass_vector, int_vector,
- sort(c(pc[1], pc[3], pc[5], pc[7]))[1], sort(c(pc[1], pc[3], pc[5], pc[7]))[4],
- resol, sum_fit = sum_fit)$fq_new
-
- # get parameter values
- area1 <- estimate_area(pc[1], resol, pc[2], sigma1, int_factor)
- area2 <- estimate_area(pc[3], resol, pc[4], sigma2, int_factor)
- area3 <- estimate_area(pc[5], resol, pc[6], sigma3, int_factor)
- area4 <- estimate_area(pc[7], resol, pc[8], sigma4, int_factor)
- peak_area <- c(peak_area, area1)
- peak_area <- c(peak_area, area2)
- peak_area <- c(peak_area, area3)
- peak_area <- c(peak_area, area4)
- peak_mean <- c(peak_mean, pc[1])
- peak_mean <- c(peak_mean, pc[3])
- peak_mean <- c(peak_mean, pc[5])
- peak_mean <- c(peak_mean, pc[7])
- peak_scale <- c(peak_scale, pc[2])
- peak_scale <- c(peak_scale, pc[4])
- peak_scale <- c(peak_scale, pc[6])
- peak_scale <- c(peak_scale, pc[8])
- peak_sigma <- c(peak_sigma, sigma1)
- peak_sigma <- c(peak_sigma, sigma2)
- peak_sigma <- c(peak_sigma, sigma3)
- peak_sigma <- c(peak_sigma, sigma4)
-
- roi_value_list <- list("mean" = peak_mean, "scale" = peak_scale, "sigma" = peak_sigma, "area" = peak_area, "qual" = fq)
- return(roi_value_list)
-}
-
diff --git a/DIMS/Utils/get_element_info.R b/DIMS/Utils/get_element_info.R
deleted file mode 100644
index a619f6a0..00000000
--- a/DIMS/Utils/get_element_info.R
+++ /dev/null
@@ -1,24 +0,0 @@
-## adapted from elementInfo.R, which is adapted from Rdisop function .getElement
-# refactor: check where library is initialised: library(Rdisop)
-get_element_info <- function(name, elements = NULL) {
- #' Get info on m/z and isotopes for all chemical elements
- #'
- #' @param name: Name of adduct, e.g. Na (string)
- #' @param elements: List of all adducts to take into account (list of strings)
- #'
- #' @return element_info: peak group list with filled-in intensities (matrix)
-
- # get information on all elements
- if (!is.list(elements) || length(elements) == 0 ) {
- elements <- initializePSE()
- }
- # extract information for a particular adduct
- if (name == "CH3OH+H") {
- # regular_expr should be exact match for name, except for methanol
- regular_expr <- "^CH3OH\\+H$"
- } else {
- regular_expr <- paste0("^", name, "$")
- }
- element_info <- elements[[grep(regular_expr, sapply(elements, function(x) { x$name }))]]
- return(element_info)
-}
\ No newline at end of file
diff --git a/DIMS/Utils/get_fit_quality.R b/DIMS/Utils/get_fit_quality.R
deleted file mode 100644
index 23b1947d..00000000
--- a/DIMS/Utils/get_fit_quality.R
+++ /dev/null
@@ -1,34 +0,0 @@
-## adapted from getFitQuality.R
-# parameter not used: mu_last
-get_fit_quality <- function(mass_vector, int_vector, mu_first, mu_last, resol, scale = NULL, sigma = NULL, sum_fit = NULL) {
- #' Fit 1 Gaussian peak in small region of m/z
- #'
- #' @param mass_vector: Vector of m/z values for a region of interest (float)
- #' @param int_vector: Value used to calculate area under Gaussian curve (integer)
- #' @param mu_first: Value for first peak (float)
- #' @param scale: Initial value used to estimate scaling parameter (integer)
- #' @param resol: Value for resolution (integer)
- #' @param sum_fit: Value indicating quality of fit of Gaussian curve (float)
- #'
- #' @return list_params: list of parameters indicating quality of fit (list)
- if (is.null(sum_fit)) {
- mass_vector_int <- mass_vector
- int_vector_int <- int_vector
- # get new fit quality
- fq_new <- mean(abs((scale * dnorm(mass_vector_int, mu_first, sigma)) - int_vector_int) /
- rep((max(scale * dnorm(mass_vector_int, mu_first, sigma)) / 2), length(mass_vector_int)))
- } else {
- sum_fit_int <- sum_fit
- int_vector_int <- int_vector
- mass_vector_int <- mass_vector
- # get new fit quality
- fq_new <- mean(abs(sum_fit_int - int_vector_int) / rep(max(sum_fit_int) /2, length(sum_fit_int)))
- }
-
- # Prevent division by 0
- if (is.nan(fq_new)) fq_new <- 1
-
- list_params <- list("fq_new" = fq_new, "x_int" = mass_vector_int, "y_int" = int_vector_int)
- return(list_params)
-}
-
diff --git a/DIMS/Utils/get_fwhm.R b/DIMS/Utils/get_fwhm.R
deleted file mode 100644
index 75439216..00000000
--- a/DIMS/Utils/get_fwhm.R
+++ /dev/null
@@ -1,21 +0,0 @@
-## adapted from getFwhm.R
-get_fwhm <- function(query_mass, resol) {
- #' Calculate fwhm (full width at half maximum intensity) for a peak
- #'
- #' @param query_mass: Value for mass (float)
- #' @param resol: Value for resolution (integer)
- #'
- #' @return fwhm: Value for full width at half maximum (float)
-
- # set aberrant values of query_mass to zero
- if (is.nan(query_mass)) query_mass <- 0
- if (is.na(query_mass)) query_mass <- 0
- if (is.null(query_mass)) query_mass <- 0
- if (query_mass < 0) query_mass <- 0
- # calculate resolution at given m/z value
- resol_mz <- resol * (1 / sqrt(2) ^ (log2(query_mass / 200)))
- # calculate full width at half maximum
- fwhm <- query_mass / resol_mz
- return(fwhm)
-}
-
diff --git a/DIMS/Utils/get_patient_data_to_helix.R b/DIMS/Utils/get_patient_data_to_helix.R
deleted file mode 100644
index 1eeaf7a6..00000000
--- a/DIMS/Utils/get_patient_data_to_helix.R
+++ /dev/null
@@ -1,39 +0,0 @@
-get_patient_data_to_helix <- function(metab_interest_sorted, metab_list_all) {
- #' Get patient data to be uploaded to Helix
- #'
- #' @param metab_interest_sorted: list of dataframes with metabolite Z-scores for each sample/patient (list)
- #' @param metab_list_all: list of tables with metabolites for Helix and violin plots (list)
- #'
- #' @return: dataframe with patient data with only metabolites for Helix and violin plots
- #' with Helix name, high/low Z-score cutoffs
-
- # Combine Z-scores of metab groups together
- df_all_metabs_zscores <- bind_rows(metab_interest_sorted)
- # Change columnnames
- colnames(df_all_metabs_zscores) <- c("HMDB_name", "Patient", "Z_score")
- # Change Patient column to character instead of factor
- df_all_metabs_zscores$Patient <- as.character(df_all_metabs_zscores$Patient)
-
- # Delete whitespaces HMDB_name
- df_all_metabs_zscores$HMDB_name <- str_trim(df_all_metabs_zscores$HMDB_name, "right")
-
- # Split HMDB_name column on "nitine;" for match dims_helix_table
- df_all_metabs_zscores$HMDB_name_split <- str_split_fixed(df_all_metabs_zscores$HMDB_name, "nitine;", 2)[, 1]
-
- # Combine stofgroepen
- dims_helix_table <- bind_rows(metab_list_all)
- # Filter table for metabolites for Helix
- dims_helix_table <- dims_helix_table %>% filter(Helix == "ja")
- # Split HMDB_name column on "nitine;" for match df_all_metabs_zscores
- dims_helix_table$HMDB_name_split <- str_split_fixed(dims_helix_table$HMDB_name, "nitine;", 2)[, 1]
- dims_helix_table <- dims_helix_table %>% select(HMDB_name_split, Helix_naam, high_zscore, low_zscore)
-
- # Filter DIMS results for metabolites for Helix
- df_metabs_helix <- df_all_metabs_zscores %>% filter(HMDB_name_split %in% dims_helix_table$HMDB_name_split)
- # Combine dims_helix_table and df_metabs_helix, adding Helix codes etc.
- df_metabs_helix <- df_metabs_helix %>% left_join(dims_helix_table, by = join_by(HMDB_name_split))
-
- df_metabs_helix <- df_metabs_helix %>% select(HMDB_name, Patient, Z_score, Helix_naam, high_zscore, low_zscore)
-
- return(df_metabs_helix)
-}
diff --git a/DIMS/Utils/get_stdev.R b/DIMS/Utils/get_stdev.R
deleted file mode 100644
index a3851876..00000000
--- a/DIMS/Utils/get_stdev.R
+++ /dev/null
@@ -1,22 +0,0 @@
-## adapted from getSD.R
-get_stdev <- function(mass_vector, int_vector, resol = 140000) {
- #' Calculate standard deviation to determine width of a peak
- #'
- #' @param mass_vector: Vector of 3 mass values (float)
- #' @param int_vector: Vector of 3 intensities (float)
- #' @param resol: Value for resolution (integer)
- #'
- #' @return stdev: Value for standard deviation
- # find maximum intensity in vector
- max_index <- which(int_vector == max(int_vector))
- # find corresponding mass at maximum intensity
- max_mass <- mass_vector[max_index]
- # calculate resolution at given m/z value
- resol_mz <- resol * (1 / sqrt(2) ^ (log2(max_mass / 200)))
- # calculate full width at half maximum
- fwhm <- max_mass / resol_mz
- # calculate standard deviation
- stdev <- (fwhm / 2) * 0.85
- return(stdev)
-}
-
diff --git a/DIMS/Utils/identify_noisepeaks.R b/DIMS/Utils/identify_noisepeaks.R
deleted file mode 100644
index cf79279f..00000000
--- a/DIMS/Utils/identify_noisepeaks.R
+++ /dev/null
@@ -1,103 +0,0 @@
-## adapted from ident.hires.noise.HPC
-# refactor: remove variables slope, incpt, ppm_iso_fixed
-# combine with function get_element_info
-# modified identify function to also look for adducts and their isotopes
-identify_noisepeaks <- function(peakgroup_list, all_adducts, scanmode = "Negative", look4 = c("Cl", "Ac"),
- noise_mz = NULL, resol = 140000, slope = 0, incpt = 0, ppm_fixed, ppm_iso_fixed) {
- #' Replace intensities that are zero with random value
- #'
- #' @param peakgroup_list: Peak group list (matrix)
- #' @param all_adducts: List of adducts to take into account (list of strings)
- #' @param scanmode: Scan mode, positive or negative (string)
- #' @param look4: List of adducts to look for (list of strings)
- #' @param noise_mz: All known noise peaks (matrix)
- #' @param resol: Value for resolution (integer)
- #' @param slope: Value for slope for mass correction (float)
- #' @param incpt: Value for intercept for mass correction (float)
- #' @param ppm_fixed: Value for distance between two values of mass (integer)
- #' @param ppm_iso_fixed: Value for distance between two values of mass for isotope peaks (integer)
- #'
- #' @return final_outlist: peak group list with filled-in intensities (matrix)
-
- options(stringsAsFactors = FALSE)
- metlin <- assi <- iso <- rep("", nrow(peakgroup_list))
- theormz <- nisos <- expint <- conf <- rep(0, nrow(peakgroup_list))
-
- # add adducts to identification list
- if (scanmode == "Positive") {
- adduct_scanmode <- "+"
- } else {
- adduct_scanmode <- "-"
- }
- # make a copy of noise_mz
- noise_mz_orig <- noise_mz
-
- # loop over type of adduct
- for (adduct_index in 1:length(look4)) {
- noise_mz_adduct <- noise_mz_orig
- noise_mz_adduct[, "CompoundName"] <- as.character(noise_mz_orig[, "CompoundName"])
-
- if (look4[adduct_index] == "H2O") {
- add2label <- paste0("[M-", look4[adduct_index], "]", adduct_scanmode)
- } else {
- add2label <- paste0("[M+", look4[adduct_index], "]", adduct_scanmode)
- }
-
- noise_mz_adduct[, "CompoundName"] <- paste0(noise_mz_adduct[, "CompoundName"], add2label)
- adduct_info <- get_element_info(look4[adduct_index], all_adducts)
- if (scanmode == "Positive") {
- adduct_mass <- adduct_info$mass[1] + adduct_info$isotope$mass[1] - hydrogen_mass
- } else {
- adduct_mass <- adduct_info$mass[1] + adduct_info$isotope$mass[1] + hydrogen_mass
- }
-
- # loop over compounds in database
- for (compound_index in 1:nrow(noise_mz_adduct)) {
- # construct information for compound + adduct:
- if (scanmode == "Positive") {
- noise_mz_adduct[compound_index, "Mpos"] <- as.numeric(noise_mz_adduct[compound_index, "Mpos"]) + adduct_mass
- noise_mz_adduct[compound_index, "MNeg"] <- 0
- } else {
- noise_mz_adduct[compound_index, "Mpos"] <- 0
- noise_mz_adduct[compound_index, "MNeg"] <- as.numeric(noise_mz_adduct[compound_index, "MNeg"]) + adduct_mass
- }
- }
- noise_mz <- rbind(noise_mz, noise_mz_adduct)
- }
-
- if (scanmode == "Positive") {
- theor_mcol <- as.numeric(noise_mz[, "Mpos"])
- } else {
- theor_mcol <- as.numeric(noise_mz[, "MNeg"])
- }
-
- # get mz information from peakgroup_list
- mcol <- peakgroup_list[, "mzmed.pgrp"]
- # if column with average intensities is missing, calculate it:
- if (!("avg.int" %in% colnames(peakgroup_list))) {
- mzmaxcol <- which(colnames(peakgroup_list) == "mzmax.pgrp")
- endcol <- ncol(peakgroup_list)
- peakgroup_list[, "avg.int"] <- apply(peakgroup_list[, (mzmaxcol + 1):(endcol)], 1, mean)
- }
-
- # do indentification using own database:
- for (row_index in 1:nrow(noise_mz)) {
- theor_mz <- theor_mcol[row_index]
-
- # set tolerance for mz accuracy of main peak
- mtol <- theor_mz * ppm_fixed / 1000000
- # find main peak
- selp <- which(mcol > (theor_mz - mtol) & mcol < (theor_mz + mtol))
- # if there is more than one candidate peak for main, select best one based on mz_diff
- if (length(selp) > 1) {
- selp <- selp[abs(mcol[selp] - theor_mz) == min(abs(mcol[selp] - theor_mz))]
- }
- if (length(selp) == 1) {
- assi[selp] <- paste(assi[selp], as.character(noise_mz[row_index, "CompoundName"]), sep = ";")
- theormz[selp] <- theor_mz
- }
- }
-
- final_outlist <- cbind(peakgroup_list, assi, theormz, conf, nisos, iso, expint, metlin)
- return(final_outlist)
-}
diff --git a/DIMS/Utils/is_diagnostic_patient.R b/DIMS/Utils/is_diagnostic_patient.R
deleted file mode 100644
index 58208383..00000000
--- a/DIMS/Utils/is_diagnostic_patient.R
+++ /dev/null
@@ -1,11 +0,0 @@
-is_diagnostic_patient <- function(patient_column) {
- #' Check for Diagnostics patients with correct patient number (e.g. starting with "P2024M")
- #'
- #' @param patient_column: a column from dataframe with IDs (character vector)
- #'
- #' @return: a logical vector with TRUE or FALSE for each element (vector)
-
- diagnostic_patients <- grepl("^P[0-9]{4}M", patient_column)
-
- return(diagnostic_patients)
-}
diff --git a/DIMS/Utils/merge_duplicate_rows.R b/DIMS/Utils/merge_duplicate_rows.R
deleted file mode 100644
index 25afd312..00000000
--- a/DIMS/Utils/merge_duplicate_rows.R
+++ /dev/null
@@ -1,58 +0,0 @@
-## adapted from mergeDuplicatedRows.R
-merge_duplicate_rows <- function(peakgroup_list) {
- #' Merge identification info for peak groups with the same mass
- #'
- #' @param peakgroup_list: Peak group list (matrix)
- #'
- #' @return peakgroup_list_dedup: de-duplicated peak group list (matrix)
-
- collapse <- function(column_label, peakgroup_list, index_dup) {
- #' Collapse identification info for peak groups with the same mass
- #'
- #' @param column_label: Name of column in peakgroup_list (string)
- #' @param peakgroup_list: Peak group list (matrix)
- #' @param index_dup: Index of duplicate peak group (integer)
- #'
- #' @return collapsed_items: Semicolon-separated list of info (string)
- # get the item(s) that need to be collapsed
- list_items <- as.vector(peakgroup_list[index_dup, column_label])
- # remove NA
- if (length(which(is.na(list_items))) > 0) list_items <- list_items[-which(is.na(list_items))]
- collapsed_items <- paste(list_items, collapse = ";")
- return(collapsed_items)
- }
-
- options(digits = 16)
- collect <- NULL
- remove <- NULL
-
- # check for peak groups with identical mass
- index_dup <- which(duplicated(peakgroup_list[, "mzmed.pgrp"]))
-
- while (length(index_dup) > 0) {
- # get the index for the peak group which is double
- peaklist_index <- which(peakgroup_list[, "mzmed.pgrp"] == peakgroup_list[index_dup[1], "mzmed.pgrp"])
- single_peakgroup <- peakgroup_list[peaklist_index[1], , drop = FALSE]
-
- # use function collapse to concatenate info
- single_peakgroup[, "assi_HMDB"] <- collapse("assi_HMDB", peakgroup_list, peaklist_index)
- single_peakgroup[, "iso_HMDB"] <- collapse("iso_HMDB", peakgroup_list, peaklist_index)
- single_peakgroup[, "HMDB_code"] <- collapse("HMDB_code", peakgroup_list, peaklist_index)
- single_peakgroup[, "all_hmdb_ids"] <- collapse("all_hmdb_ids", peakgroup_list, peaklist_index)
- single_peakgroup[, "sec_hmdb_ids"] <- collapse("sec_hmdb_ids", peakgroup_list, peaklist_index)
- if (single_peakgroup[, "sec_hmdb_ids"] == ";") single_peakgroup[, "sec_hmdb_ids"] < NA
-
- # keep track of deduplicated entries
- collect <- rbind(collect, single_peakgroup)
- remove <- c(remove, peaklist_index)
-
- # remove current entry from index
- index_dup <- index_dup[-which(peakgroup_list[index_dup, "mzmed.pgrp"] == peakgroup_list[index_dup[1], "mzmed.pgrp"])]
- }
-
- # remove duplicate entries
- if (!is.null(remove)) peakgroup_list <- peakgroup_list[-remove, ]
- # append deduplicated entries
- peakgroup_list_dedup <- rbind(peakgroup_list, collect)
- return(peakgroup_list_dedup)
-}
diff --git a/DIMS/Utils/optimize_gaussfit.R b/DIMS/Utils/optimize_gaussfit.R
deleted file mode 100644
index 2d7f95da..00000000
--- a/DIMS/Utils/optimize_gaussfit.R
+++ /dev/null
@@ -1,23 +0,0 @@
-## adapted from optimizeGauss.R
-optimize_gaussfit <- function(mass_vector, int_vector, sigma, mass_max) {
- #' Optimize fit of Gaussian curve to small region of m/z
- #'
- #' @param mass_vector: Vector of m/z values for a region of interest (float)
- #' @param int_vector: Vector of intensities for a region of interest (float)
- #' @param sigma: Value for standard deviation (float)
- #' @param mass_max: Value for mass at center of peak (float)
- #'
- #' @return opt_fit: list of fit values for region of interest (list)
-
- # define optimization function for optim based on normal distribution
- opt_f <- function(p, mass_vector, int_vector, sigma, mass_max) {
- curve <- p * dnorm(mass_vector, mass_max, sigma)
- return((max(curve) - max(int_vector))^2)
- }
-
- # get optimal value for fitted Gaussian curve
- opt_fit <- optimize(opt_f, c(0, 100000), tol = 0.0001, mass_vector, int_vector, sigma, mass_max)
-
- return(opt_fit$minimum)
-}
-
diff --git a/DIMS/Utils/output_helix.R b/DIMS/Utils/output_helix.R
deleted file mode 100644
index 09fa1a9b..00000000
--- a/DIMS/Utils/output_helix.R
+++ /dev/null
@@ -1,39 +0,0 @@
-output_for_helix <- function(protocol_name, df_metabs_helix) {
- #' Get the output dataframe for Helix
- #'
- #' @param protocol_name: protocol name (string)
- #' @param df_metabs_helix: dataframe with metabolite Z-scores for patients (dataframe)
- #'
- #' @return: dataframe with patient metabolite Z-scores in correct format for Helix
-
- # Remove positive controls
- df_metabs_helix <- df_metabs_helix %>% filter(is_diagnostic_patient(Patient))
-
- # Add 'Vial' column, each patient has unique ID
- df_metabs_helix <- df_metabs_helix %>%
- group_by(Patient) %>%
- mutate(Vial = cur_group_id()) %>%
- ungroup()
-
- # Split patient number into labnummer and Onderzoeksnummer
- df_metabs_helix <- add_lab_id_and_onderzoeksnummer(df_metabs_helix)
-
- # Add column with protocol name
- df_metabs_helix$Protocol <- protocol_name
-
- # Change name Z_score and Helix_naam columns to Amount and Name
- change_columns <- c(Amount = "Z_score", Name = "Helix_naam")
- df_metabs_helix <- df_metabs_helix %>% rename(all_of(change_columns))
-
- # Select only necessary columns and set them in correct order
- df_metabs_helix <- df_metabs_helix %>%
- select(c(Vial, labnummer, Onderzoeksnummer, Protocol, Name, Amount))
-
- # Remove duplicate patient-metabolite combinations ("leucine + isoleucine + allo-isoleucin_Z-score" is added 3 times)
- df_metabs_helix <- df_metabs_helix %>%
- group_by(Onderzoeksnummer, Name) %>%
- distinct() %>%
- ungroup()
-
- return(df_metabs_helix)
-}
diff --git a/DIMS/Utils/prepare_alarmvalues.R b/DIMS/Utils/prepare_alarmvalues.R
deleted file mode 100644
index 94ffbfd2..00000000
--- a/DIMS/Utils/prepare_alarmvalues.R
+++ /dev/null
@@ -1,59 +0,0 @@
-prepare_alarmvalues <- function(pt_name, dims_helix_table) {
- #' Create a dataframe with all metabolites that exceed the min and max Z-score cutoffs
- #'
- #' @param pt_name: patient code (string)
- #' @param dims_helix_table: dataframe with metabolite Z-scores for each patient and Helix info (dataframe)
- #'
- #' @return: dataframe with metabolites that exceed the min and max Z-score cutoffs for the selected patient
-
- # extract data for patient of interest (pt_name)
- pt_metabs_helix <- dims_helix_table %>% filter(Patient == pt_name)
- pt_metabs_helix$Z_score <- round(pt_metabs_helix$Z_score, 2)
-
- # Make empty dataframes for metabolites above or below alarmvalues
- pt_list_high <- data.frame(HMDB_name = character(), Z_score = numeric())
- pt_list_low <- data.frame(HMDB_name = character(), Z_score = numeric())
-
- # Loop over individual metabolites
- for (metab in unique(pt_metabs_helix$HMDB_name)){
- # Get data for individual metabolite
- pt_metab <- pt_metabs_helix %>% filter(HMDB_name == metab)
- # print(pt_metab)
-
- # Check if zscore is positive of negative
- if (pt_metab$Z_score > 0) {
- # Get specific alarmvalue for metabolite
- high_zscore_cutoff_metab <- pt_metabs_helix %>% filter(HMDB_name == metab) %>% pull(high_zscore)
-
- # If zscore is above the alarmvalue, add to pt_list_high table
- if (pt_metab$Z_score > high_zscore_cutoff_metab) {
- pt_metab_high <- pt_metab %>% select(HMDB_name, Z_score)
- pt_list_high <- rbind(pt_list_high, pt_metab_high)
- }
- } else {
- # Get specific alarmvalue for metabolite
- low_zscore_cutoff_metab <- pt_metabs_helix %>% filter(HMDB_name == metab) %>% pull(low_zscore)
-
- # If zscore is below the alarmvalue, add to pt_list_low table
- if (pt_metab$Z_score < low_zscore_cutoff_metab) {
- pt_metab_low <- pt_metab %>% select(HMDB_name, Z_score)
- pt_list_low <- rbind(pt_list_low, pt_metab_low)
- }
- }
- }
-
- # sort tables on zscore
- pt_list_high <- pt_list_high %>% arrange(desc(Z_score))
- pt_list_low <- pt_list_low %>% arrange(Z_score)
- # add lines for increased, decreased
- extra_line1 <- c("Increased", "")
- extra_line2 <- c("Decreased", "")
- # combine the two lists
- top_metab_pt <- rbind(extra_line1, pt_list_high, extra_line2, pt_list_low)
- # remove row names
- rownames(top_metab_pt) <- NULL
- # change column names for display
- colnames(top_metab_pt) <- c("Metabolite", "Z-score")
-
- return(top_metab_pt)
-}
diff --git a/DIMS/Utils/prepare_data.R b/DIMS/Utils/prepare_data.R
deleted file mode 100644
index ea3efe5b..00000000
--- a/DIMS/Utils/prepare_data.R
+++ /dev/null
@@ -1,51 +0,0 @@
-# unused variables will be removed: metab_list_alarm
-prepare_data <- function(metab_list_all, zscore_patients_local) {
- #' Combine patient Z-scores with metabolite info
- #'
- #' @param metab_list_all: list of dataframes with metabolite information for different stofgroepen (list)
- #' @param zscore_patients_local: dataframe with metabolite Z-scores for all patient
- #'
- #' @return: list of dataframes for each stofgroep with data for each metabolite and patient/control per row
-
- # remove "_Zscore" from column (patient) names
- colnames(zscore_patients_local) <- gsub("_Zscore", "", colnames(zscore_patients_local))
- # put data into pages, max 20 violin plots per page in PDF
- metab_interest_sorted <- list()
- metab_category <- c()
- for (metab_class_index in 1:length(metab_list_all)) {
- metab_class <- names(metab_list_all)[metab_class_index]
- metab_list <- metab_list_all[[metab_class_index]]
- if (ncol(metab_list) > 2) {
- # third column are the alarm values, so reduce the data frame to 2 columns and save list
- metab_list_alarm <- metab_list
- metab_list <- metab_list[, c(1, 2)]
- }
- # make sure that all HMDB_names have 45 characters
- for (metab_index in 1:length(metab_list$HMDB_name)) {
- if (is.character(metab_list$HMDB_name[metab_index])) {
- hmdb_name_separated <- strsplit(metab_list$HMDB_name[metab_index], "")[[1]]
- } else {
- hmdb_name_separated <- "strspliterror"
- }
- if (length(hmdb_name_separated) <= 45) {
- hmdb_name_separated <- c(hmdb_name_separated, rep(" ", 45 - length(hmdb_name_separated)))
- } else {
- hmdb_name_separated <- c(hmdb_name_separated[1:42], "...")
- }
- metab_list$HMDB_name[metab_index] <- paste0(hmdb_name_separated, collapse = "")
- }
- # find metabolites and ratios in data frame zscore_patients_local
- metab_interest <- inner_join(metab_list, zscore_patients_local[-2], by = "HMDB_code")
- # remove column "HMDB_code"
- metab_interest <- metab_interest[, -which(colnames(metab_interest) == "HMDB_code")]
- # put the data frame in long format
- metab_interest_melt <- reshape2::melt(metab_interest, id.vars = "HMDB_name")
- # sort on metabolite names (HMDB_name)
- sort_order <- order(metab_interest_melt$HMDB_name)
- metab_interest_sorted[[metab_class_index]] <- metab_interest_melt[sort_order, ]
- metab_category <- c(metab_category, metab_class)
- }
- names(metab_interest_sorted) <- metab_category
-
- return(metab_interest_sorted)
-}
diff --git a/DIMS/Utils/prepare_data_perpage.R b/DIMS/Utils/prepare_data_perpage.R
deleted file mode 100644
index 61d2e15c..00000000
--- a/DIMS/Utils/prepare_data_perpage.R
+++ /dev/null
@@ -1,58 +0,0 @@
-# remove default variable values
-prepare_data_perpage <- function(metab_interest_sorted, metab_interest_contr, nr_plots_perpage, nr_pat = 20, nr_contr = 30) {
- #' Combine patient and control data for each page of the violinplot pdf
- #'
- #' @param metab_interest_sorted: list of dataframes with data for each metabolite and patient (list)
- #' @param metab_interest_contr: list of dataframes with data for each metabolite and control (list)
- #' @param nr_plots_perpage: number of plots per page in the violinplot pdf (integer)
- #' @param nr_pat: number of patients (integer)
- #' @param nr_contr: number of controls (integer)
- #'
- #' @return: list of dataframes with metabolite Z-scores for each patient and control,
- #' the length of list is the number of pages for the violinplot pdf (list)
-
- total_nr_pages <- 0
- metab_perpage <- list()
- metab_category <- c()
- for (metab_class_index in 1:length(metab_interest_sorted)) {
- # split list into pages, each page containing max nr_plots_perpage (20) compounds
- metab_interest_perclass <- metab_interest_sorted[[metab_class_index]]
- metab_class <- names(metab_interest_sorted)[metab_class_index]
- # add controls
- metab_interest_contr_perclass <- metab_interest_contr[[metab_class_index]]
- # number of pages for this class
- nr_pages <- ceiling(length(unique(metab_interest_perclass$HMDB_name)) / nr_plots_perpage)
- for (page_nr in 1:nr_pages) {
- total_nr_pages <- total_nr_pages + 1
- select_rows_start <- (nr_pat * nr_plots_perpage * (page_nr - 1)) + 1
- select_rows_end <- nr_pat * nr_plots_perpage * page_nr
- metab_onepage_pat <- metab_interest_perclass[select_rows_start:select_rows_end, ]
- # same for controls
- select_rows_start_contr <- (nr_contr * nr_plots_perpage * (page_nr - 1)) + 1
- select_rows_end_contr <- nr_contr * nr_plots_perpage * page_nr
- metab_onepage_pcontr <- metab_interest_contr_perclass[select_rows_start_contr:select_rows_end_contr, ]
- # add controls
- metab_onepage <- rbind(metab_onepage_pat, metab_onepage_pcontr)
- # if a page has fewer than nr_plots_perpage plots, fill page with empty plots
- na_rows <- which(is.na(metab_onepage$HMDB_name))
- if (length(na_rows) > 0) {
- # repeat the patient and control variables
- metab_onepage$variable[na_rows] <- metab_onepage$variable[1:(nr_pat + nr_contr)]
- # for HMDB name, substitute a number of spaces
- for (row_nr in na_rows) {
- metab_onepage$HMDB_name[row_nr] <- paste0(rep("_", ceiling(row_nr / (nr_pat + nr_contr))), collapse = "")
- }
- metab_onepage$HMDB_name <- gsub("_", " ", metab_onepage$HMDB_name)
- # leave the values at NA
- }
- # put data for one page into object with data for all pages
- metab_perpage[[total_nr_pages]] <- metab_onepage
- # create list of page headers
- metab_category <- c(metab_category, paste(metab_class, page_nr, sep = "_"))
- }
- }
- # add page headers to list
- names(metab_perpage) <- metab_category
-
- return(metab_perpage)
-}
diff --git a/DIMS/Utils/prepare_toplist.R b/DIMS/Utils/prepare_toplist.R
deleted file mode 100644
index cc616c40..00000000
--- a/DIMS/Utils/prepare_toplist.R
+++ /dev/null
@@ -1,35 +0,0 @@
-prepare_toplist <- function(pt_name, zscore_patients_copy) {
- #' Create a dataframe with the top 20 highest and top 10 lowest metabolites per patient
- #'
- #' @param pt_name: patient code (string)
- #' @param zscore_patients_copy: dataframe with metabolite Z-scores per patient (dataframe)
- #'
- #' @return: dataframe with 30 metabolites and Z-scores (dataframe)
-
- # set parameters for table
- top_highest <- 20
- top_lowest <- 10
-
- # extract data for patient of interest (pt_name)
- pt_list <- zscore_patients_copy[, c(1, 2, which(colnames(zscore_patients_copy) == pt_name))]
- # sort metabolites on Z-scores for this patient
- pt_list_sort <- sort(pt_list[, 3], index.return = TRUE)
- # determine top highest and lowest Z-scores for this patient
- pt_list_sort <- sort(pt_list[, 3], index.return = TRUE)
- pt_list_low <- pt_list[pt_list_sort$ix[1:top_lowest], ]
- pt_list_high <- pt_list[pt_list_sort$ix[length(pt_list_sort$ix):(length(pt_list_sort$ix) - top_highest + 1)], ]
- # round off Z-scores
- pt_list_low[, 3] <- round(as.numeric(pt_list_low[, 3]), 2)
- pt_list_high[, 3] <- round(as.numeric(pt_list_high[, 3]), 2)
- # add lines for increased, decreased
- extra_line1 <- c("Increased", "", "")
- extra_line2 <- c("Decreased", "", "")
- top_metab_pt <- rbind(extra_line1, pt_list_high, extra_line2, pt_list_low)
- # remove row names
- rownames(top_metab_pt) <- NULL
-
- # change column names for display
- colnames(top_metab_pt) <- c("HMDB_ID", "Metabolite", "Z-score")
-
- return(top_metab_pt)
-}
diff --git a/DIMS/Utils/replace_zeros.R b/DIMS/Utils/replace_zeros.R
deleted file mode 100644
index 0ef06a3c..00000000
--- a/DIMS/Utils/replace_zeros.R
+++ /dev/null
@@ -1,64 +0,0 @@
-## adapted from replaceZeros.R
-# this function does two things: replace zeros with random value and identify noise peaks
-# refactor: split into two functions
-# remove parameters outdir and thresh
-# make hard-coded path to file with noise peaks into variable
-# remove variables outdir, thresh
-replace_zeros <- function(peakgroup_list, repl_pattern, scanmode, resol, outdir, thresh, ppm) {
- #' Replace intensities that are zero with random value
- #'
- #' @param peakgroup_list: Peak group list (matrix)
- #' @param repl_pattern: Replication pattern (list of strings)
- #' @param scanmode: Scan mode, positive or negative (string)
- #' @param resol: Value for resolution (integer)
- #' @param outdir: Path for output directory (string)
- #' @param thresh: Value for threshold (integer)
- #' @param ppm: Value for distance between two values of mass (integer)
- #'
- #' @return final_outlist: peak group list with filled-in intensities (matrix)
-
- # replace zeros
- if (!is.null(peakgroup_list)) {
- for (sample_index in 1:length(names(repl_pattern))) {
- sample_peaks <- peakgroup_list[, names(repl_pattern)[sample_index]]
- zero_intensity <- which(sample_peaks <= 0)
- if (!length(zero_intensity)) {
- next
- }
- for (zero_index in 1:length(zero_intensity)) {
- area <- fit_optim(peakgroup_list[zero_intensity[zero_index], "mzmed.pgrp"], thresh,
- resol, FALSE, scanmode, int_factor = 1 * 10^5, 1, 1)$area
- peakgroup_list[zero_intensity[zero_index], names(repl_pattern)[sample_index]] <- rnorm(n = 1, mean = area,
- sd = 0.25 * area)
- }
- }
-
- # Add column with average intensity
- peakgroup_list <- cbind(peakgroup_list, "avg.int" = apply(peakgroup_list[, 7:(ncol(peakgroup_list) - 4)], 1, mean))
-
- if (scanmode == "negative") {
- label <- "MNeg"
- label2 <- "Negative"
- # look for adducts in negative mode
- look4_adducts <- c("Cl", "Cl37", "For", "NaCl", "KCl", "H2PO4", "HSO4", "Na-H", "K-H", "H2O", "I")
- } else {
- label <- "Mpos"
- label2 <- "Positive"
- # look for adducts in positive mode
- look4_adducts <- c("Na", "K", "NaCl", "NH4", "2Na-H", "CH3OH", "KCl", "NaK-H")
- }
-
- # Identify noise peaks
- noise_mz <- read.table(file = "/hpc/dbg_mz/tools/db/TheoreticalMZ_NegPos_incNaCl.txt",
- sep = "\t", header = TRUE, quote = "")
- noise_mz <- noise_mz[(noise_mz[, label] != 0), 1:4]
- outlist_withnoise <- identify_noisepeaks(peakgroup_list, all_adducts, scanmode = label2,
- noise_mz, look4 = look4_adducts, resol = resol,
- slope = 0, incpt = 0, ppm_fixed = ppm, ppm_iso_fixed = ppm)
- noise_info <- outlist_withnoise[, c("assi", "theormz")]
- colnames(noise_info) <- c("assi_noise", "theormz_noise")
-
- final_outlist <- cbind(peakgroup_list, noise_info)
- return(final_outlist)
- }
-}
diff --git a/DIMS/Utils/search_mzrange.R b/DIMS/Utils/search_mzrange.R
deleted file mode 100644
index dcdc42a2..00000000
--- a/DIMS/Utils/search_mzrange.R
+++ /dev/null
@@ -1,180 +0,0 @@
-## adapted from searchMZRange.R
-# variables with fixed values will be removed from function parameters
-# int_factor, scale, outdir, plot, thresh, width, height
-# allpeaks_values should be generated here, not passed on from do_peakfinding
-search_mzrange <- function(ints_fullrange, allpeaks_values, int_factor, scale, resol,
- outdir, sample_name, scanmode, plot, width, height, thresh) {
- #' Divide the full m/z range into regions of interest with min, max and mean m/z
- #'
- #' @param ints_fullrange: Named list of intensities (float)
- #' @param allpeaks_values: Empty list to store results for all peaks
- #' @param int_factor: Value used to calculate area under Gaussian curve (integer)
- #' @param scale: Initial value used to estimate scaling parameter (integer)
- #' @param resol: Value for resolution (integer)
- #' @param outdir: Path for output directory (string)
- #' @param sample_name: Sample name (string)
- #' @param scanmode: Scan mode, positive or negative (string)
- #' @param plot: Parameter indicating whether plots should be made (boolean)
- #' @param width: Value for width of plot (integer)
- #' @param height: Value for height of plot (integer)
- #' @param thresh: Value for noise level threshold (integer)
- #'
- #' @return allpeaks_values: list of m/z regions of interest
-
- # find indices where intensity is not equal to zero
- nonzero_indices <- as.vector(which(ints_fullrange != 0))
-
- # bad infusion. These should have been taken out in AverageTechReplicates
- if (length(nonzero_indices) == 0) return(allpeaks_values)
-
- # initialize
- end_index <- NULL
- start_index <- nonzero_indices[1]
- # maximum length of region of interest
- max_roi_length <- 15
-
- # find regions of interest
- for (mz_index in 1:length(nonzero_indices)) {
- # check whether mz_index is smaller than length(nonzero_indices).
- # only false if mz_index == length(nonzero_indixes).
- # second check is true at the end of a peak.
- if (mz_index < length(nonzero_indices) && (nonzero_indices[mz_index + 1] - nonzero_indices[mz_index]) > 1) {
- end_index <- nonzero_indices[mz_index]
- # get m/z values and intensities for this region of interest
- mass_vector <- as.numeric(names(ints_fullrange)[c(start_index:end_index)])
- int_vector <- as.vector(ints_fullrange[c(start_index:end_index)])
- # check whether the vector of intensities is not empty.
- if (length(int_vector) != 0) {
- # check if intensity is above threshold or the maximum intensity is NaN
- if (max(int_vector) < thresh || is.nan(max(int_vector))) {
- # go to next region of interest
- start_index <- nonzero_indices[mz_index + 1]
- next
- }
- # check if there are more intensities than maximum for region of interest
- if (length(int_vector) > max_roi_length) {
- # trim lowest intensities to zero
- int_vector[which(int_vector < min(int_vector) * 1.1)] <- 0
- # split the range into multiple sub ranges
- sub_range <- int_vector
- names(sub_range) <- mass_vector
- allpeaks_values <- search_mzrange(sub_range, allpeaks_values, int_factor,
- scale, resol, outdir, sample_name, scanmode,
- plot, width, height, thresh)
- # A proper peak needs to have at least 3 intensities above threshold
- } else if (length(int_vector) > 3) {
- # check if the sum of intensities is above zero. Why is this necessary?
- if (sum(int_vector) == 0) next
- # get initial fit values
- roi_values <- fit_init(mass_vector, int_vector, int_factor, scale, resol,
- outdir, sample_name, scanmode, plot, width, height,
- mz_index, start_index, end_index)
- print(roi_values)
- if (roi_values$qual[1] == 1) {
- # get optimized fit values
- roi_values <- fit_optim(mass_vector, int_vector, resol, plot,
- scanmode, int_factor, width, height)
- # add region of interest to list of all peaks
- allpeaks_values$mean <- c(allpeaks_values$mean, roi_values$mean)
- allpeaks_values$area <- c(allpeaks_values$area, roi_values$area)
- allpeaks_values$nr <- c(allpeaks_values$nr, sample_name)
- allpeaks_values$min <- c(allpeaks_values$min, roi_values$min)
- allpeaks_values$max <- c(allpeaks_values$max, roi_values$max)
- allpeaks_values$qual <- c(allpeaks_values$qual, 0)
- allpeaks_values$spikes <- allpeaks_values$spikes + 1
-
- } else {
- for (j in 1:length(roi_values$mean)){
- allpeaks_values$mean <- c(allpeaks_values$mean, roi_values$mean[j])
- allpeaks_values$area <- c(allpeaks_values$area, roi_values$area[j])
- allpeaks_values$nr <- c(allpeaks_values$nr, sample_name)
- allpeaks_values$min <- c(allpeaks_values$min, roi_values$min[1])
- allpeaks_values$max <- c(allpeaks_values$max, roi_values$max[1])
- allpeaks_values$qual <- c(allpeaks_values$qual, roi_values$qual[1])
- }
- }
-
- } else {
-
- roi_values <- fit_optim(mass_vector, int_vector, resol,
- plot, scanmode, int_factor, width, height)
- allpeaks_values$mean <- c(allpeaks_values$mean, roi_values$mean)
- allpeaks_values$area <- c(allpeaks_values$area, roi_values$area)
- allpeaks_values$nr <- c(allpeaks_values$nr, sample_name)
- allpeaks_values$min <- c(allpeaks_values$min, roi_values$min)
- allpeaks_values$max <- c(allpeaks_values$max, roi_values$max)
- allpeaks_values$qual <- c(allpeaks_values$qual, 0)
- allpeaks_values$spikes <- allpeaks_values$spikes + 1
- }
- }
- start_index <- nonzero_indices[mz_index + 1]
- }
- }
-
- # last little range
- end_index <- nonzero_indices[length(nonzero_indices)]
- mass_vector <- as.numeric(names(ints_fullrange)[c(start_index:end_index)])
- int_vector <- as.vector(ints_fullrange[c(start_index:end_index)])
-
- if (length(int_vector) != 0) {
- # check if intensity above threshold
- if (max(int_vector) < thresh || is.nan(max(int_vector))) {
- # do nothing
- } else {
- # check if there are more intensities than maximum for region of interest
- if (length(int_vector) > max_roi_length) {
- # trim lowest intensities to zero
- int_vector[which(int_vector < min(int_vector) * 1.1)] <- 0
- # split the range into multiple sub ranges
- sub_range <- int_vector
- names(sub_range) <- mass_vector
-
- allpeaks_values <- search_mzrange(sub_range, allpeaks_values, int_factor, scale, resol,
- outdir, sample_name, scanmode,
- plot, width, height, thresh)
-
- } else if (length(int_vector) > 3) {
- # Check only zeros
- if (sum(int_vector) == 0) next
-
- roi_values <- fit_init(mass_vector, int_vector, int_factor, scale, resol,
- outdir, sample_name, scanmode, plot, width, height,
- mz_index, start_index, end_index)
- if (roi_values$qual[1] == 1) {
- roi_values <- fit_optim(mass_vector, int_vector, resol,
- plot, scanmode, int_factor, width, height)
-
- allpeaks_values$mean <- c(allpeaks_values$mean, roi_values$mean)
- allpeaks_values$area <- c(allpeaks_values$area, roi_values$area)
- allpeaks_values$nr <- c(allpeaks_values$nr, sample_name)
- allpeaks_values$min <- c(allpeaks_values$min, roi_values$min)
- allpeaks_values$max <- c(allpeaks_values$max, roi_values$max)
- allpeaks_values$qual <- c(allpeaks_values$qual, 0)
- allpeaks_values$spikes <- allpeaks_values$spikes + 1
-
- } else {
- for (j in 1:length(roi_values$mean)){
- allpeaks_values$mean <- c(allpeaks_values$mean, roi_values$mean[j])
- allpeaks_values$area <- c(allpeaks_values$area, roi_values$area[j])
- allpeaks_values$nr <- c(allpeaks_values$nr, sample_name)
- allpeaks_values$min <- c(allpeaks_values$min, roi_values$min[1])
- allpeaks_values$max <- c(allpeaks_values$max, roi_values$max[1])
- allpeaks_values$qual <- c(allpeaks_values$qual, roi_values$qual[1])
- }
- }
- } else {
- roi_values <- fit_optim(mass_vector, int_vector, resol,
- plot, scanmode, int_factor, width, height)
- allpeaks_values$mean <- c(allpeaks_values$mean, roi_values$mean)
- allpeaks_values$area <- c(allpeaks_values$area, roi_values$area)
- allpeaks_values$nr <- c(allpeaks_values$nr, sample_name)
- allpeaks_values$min <- c(allpeaks_values$min, roi_values$min)
- allpeaks_values$max <- c(allpeaks_values$max, roi_values$max)
- allpeaks_values$qual <- c(allpeaks_values$qual, 0)
- allpeaks_values$spikes <- allpeaks_values$spikes + 1
- }
- }
- }
- return(allpeaks_values)
-}
-
diff --git a/DIMS/Utils/sum_curves.R b/DIMS/Utils/sum_curves.R
deleted file mode 100644
index 542026b6..00000000
--- a/DIMS/Utils/sum_curves.R
+++ /dev/null
@@ -1,35 +0,0 @@
-## adapted from sumCurves.R
-# variables with fixed values will be removed from function parameters
-# plot
-# parameter half_max not used
-sum_curves <- function(mean1, mean2, scale1, scale2, sigma1, sigma2, mass_vector2, mass_vector, resol, plot) {
- #' Sum two curves
- #'
- #' @param mean1: Value for mean m/z of first peak (float)
- #' @param mean2: Value for mean m/z of second peak (float)
- #' @param scale1: Initial value used to estimate scaling parameter for first peak (integer)
- #' @param scale2: Initial value used to estimate scaling parameter for second peak (integer)
- #' @param sigma1: Value for standard deviation for first peak (float)
- #' @param sigma2: Value for standard deviation for second peak (float)
- #' @param mass_vector2: Vector of equally spaced m/z values (float)
- #' @param mass_vector: Vector of m/z values for a region of interest (float)
- #' @param resol: Value for resolution (integer)
- #' @param plot: Parameter indicating whether plots should be made (boolean)
- #'
- #' @return list_params: list of parameters indicating quality of fit (list)
-
- sum_fit <- (scale1 * dnorm(mass_vector2, mean1, sigma1)) + (scale2 * dnorm(mass_vector2, mean2, sigma2))
-
- mean1_plus2 <- weighted.mean(c(mean1, mean2), c(max(scale1 * dnorm(mass_vector2, mean1, sigma1)),
- max(scale2 * dnorm(mass_vector2, mean2, sigma2))))
-
- # get new values for parameters
- fwhm <- get_fwhm(mean1_plus2, resol)
- area <- max(sum_fit)
- scale <- scale1 + scale2
- sigma <- (fwhm / 2) * 0.85
-
- list_params <- list("mean" = mean1_plus2, "area" = area, "scale" = scale, "sigma" = sigma)
- return(list_params)
-}
-
diff --git a/DIMS/Utils/sum_intensities_adducts.R b/DIMS/Utils/sum_intensities_adducts.R
deleted file mode 100644
index 7b4734c9..00000000
--- a/DIMS/Utils/sum_intensities_adducts.R
+++ /dev/null
@@ -1,76 +0,0 @@
-sum_intensities_adducts <- function(peakgroup_list, hmdb_part, adducts, z_score) {
- #' Sum intensities for different adducts of the same metabolite
- #'
- #' @param peakgroup_list: Peak group list (matrix)
- #' @param hmdb_part: Matrix of metabolites , part of the HMDB (matrix)
- #' @param adducts: Vector of adducts (vector of integers)
- #' @param z_score: Value indicating whether Z-scores have been calculated (integer)
- #'
- #' @return adductsum: peak group list with summed intensities (matrix)
- hmdb_codes <- rownames(hmdb_part)
- hmdb_names <- hmdb_part[, 1]
-
- # create overview of row indices for each metabolite_adduct combination in peaklist
- hmdb_in_peaklist <- peakgroup_list$HMDB_code
- # avoid rows with only "" in HMDB_code column
- hmdb_in_peaklist[which(hmdb_in_peaklist == "")] <- ";"
- hmdb_in_peaklist_rownr <- c()
- for (row_nr in 1:length(hmdb_in_peaklist)) {
- hmdb_split <- strsplit(hmdb_in_peaklist[row_nr], ";")[[1]]
- hmdb_rownr <- cbind(hmdb_split, row_nr)
- hmdb_in_peaklist_rownr <- rbind(hmdb_in_peaklist_rownr, hmdb_rownr)
- }
- hmdb_in_peaklist_rownr <- as.data.frame(hmdb_in_peaklist_rownr)
- # remove NA, if any
- if (sum(is.na(hmdb_in_peaklist_rownr$hmdb_split)) > 0 ) {
- hmdb_in_peaklist_rownr <- hmdb_in_peaklist_rownr[-which(is.na(hmdb_in_peaklist_rownr$hmdb_split)), ]
- }
-
- # find intensity columns in peakgroup_list
- if (z_score == 1) {
- int_cols_C <- grep("C", colnames(peakgroup_list)[1:which(colnames(peakgroup_list) == "avg.ctrls")])
- int_cols_P <- grep("P", colnames(peakgroup_list)[1:which(colnames(peakgroup_list) == "avg.ctrls")])
- int_cols <- c(int_cols_C, int_cols_P)
- } else {
- int_cols_start <- which(colnames(peakgroup_list) == "nrsamples") + 1
- int_cols_end <- which(colnames(peakgroup_list) == "assi_HMDB") - 1
- int_cols <- c(int_cols_start:int_cols_end)
- }
-
- # initialize
- names <- NULL
- adductsum <- NULL
- names_long <- NULL
-
- # find adducts of each metabolite and sum the intensities
- if (length(hmdb_codes) == 0) {
- return(adductsum)
- }
-
- for (hmdb_index in 1:length(hmdb_codes)) {
- compound <- hmdb_codes[hmdb_index]
- compound_plus <- c(compound, paste(compound, adducts, sep = "_"))
-
- # find indices of rows in peakgroup_list that contain compound plus adducts
- metab_row <- which(hmdb_in_peaklist_rownr$hmdb_split %in% compound_plus)
- metab_indices <- as.numeric(hmdb_in_peaklist_rownr$row_nr[metab_row])
-
- # find intensities and sum them
- ints <- peakgroup_list[metab_indices, int_cols]
- total <- apply(ints, 2, sum)
-
- # add to adductsum
- if (sum(total) != 0) {
- names <- c(names, compound)
- adductsum <- rbind(adductsum, total)
- names_long <- c(names_long, hmdb_names[hmdb_index])
- }
- }
-
- if (!is.null(adductsum)) {
- rownames(adductsum) <- names
- adductsum <- cbind(adductsum, "HMDB_name" = names_long)
- }
-
- return(adductsum)
-}
diff --git a/DIMS/Utils/within_ppm.R b/DIMS/Utils/within_ppm.R
deleted file mode 100644
index abdb0d44..00000000
--- a/DIMS/Utils/within_ppm.R
+++ /dev/null
@@ -1,64 +0,0 @@
-## adapted from isWithinXppm.R
-# variables with fixed values will be removed from function parameters
-# plot
-within_ppm <- function(mean, scale, sigma, area, mass_vector2, mass_vector, ppm = 4, resol, plot) {
- #' Test whether two mass ranges are within ppm distance of each other
- #'
- #' @param mean: Value for mean m/z (float)
- #' @param scale: Initial value used to estimate scaling parameter (integer)
- #' @param sigma: Value for standard deviation (float)
- #' @param area: Value for area under the curve (float)
- #' @param mass_vector2: Vector of equally spaced m/z values (float)
- #' @param mass_vector: Vector of m/z values for a region of interest (float)
- #' @param ppm: Value for distance between two values of mass (integer)
- #' @param resol: Value for resolution (integer)
- #' @param plot: Parameter indicating whether plots should be made (boolean)
- #'
- #' @return list_params: list of parameters indicating quality of fit (list)
-
- # sort
- index <- order(mean)
- mean <- mean[index]
- scale <- scale[index]
- sigma <- sigma[index]
- area <- area[index]
-
- summed <- NULL
- remove <- NULL
-
- if (length(mean) > 1) {
- for (i in 2:length(mean)) {
- if ((abs(mean[i - 1] - mean[i]) / mean[i - 1]) * 10^6 < ppm) {
-
- # avoid double occurance in sum
- if ((i - 1) %in% summed) next
-
- result_values <- sum_curves(mean[i - 1], mean[i], scale[i - 1], scale[i], sigma[i - 1], sigma[i],
- mass_vector2, mass_vector, resol, plot)
- summed <- c(summed, i - 1, i)
- if (is.nan(result_values$mean)) result_values$mean <- 0
- mean[i - 1] <- result_values$mean
- mean[i] <- result_values$mean
- area[i - 1] <- result_values$area
- area[i] <- result_values$area
- scale[i - 1] <- result_values$scale
- scale[i] <- result_values$scale
- sigma[i - 1] <- result_values$sigma
- sigma[i] <- result_values$sigma
-
- remove <- c(remove, i)
- }
- }
- }
-
- if (length(remove) != 0) {
- mean <- mean[-c(remove)]
- area <- area[-c(remove)]
- scale <- scale[-c(remove)]
- sigma <- sigma[-c(remove)]
- }
-
- list_params <- list("mean" = mean, "area" = area, "scale" = scale, "sigma" = sigma, "qual" = NULL)
- return(list_params)
-}
-
diff --git a/DIMS/export/generate_excel_functions.R b/DIMS/export/generate_excel_functions.R
index 9cf5936e..2e6481dd 100644
--- a/DIMS/export/generate_excel_functions.R
+++ b/DIMS/export/generate_excel_functions.R
@@ -1,132 +1,141 @@
# Functions for GenerateExcel
-get_intensities_cols <- function(outlist, label) {
- #' Get the indices of the control columns and a dataframe of intensities of the controls
- #'
- #' @param outlist: dataframe with intensities for all samples
- #' @param label: string used by grep to get the correct columns
- #'
- #' @returns: list with 2 items:
- #' col_idx: vector with indices of the control columns
- #' df_intensities: dataframe with the intensities of the controls
- col_idx <- grep(label, colnames(outlist), fixed = TRUE)
- df_intensities <- as.data.frame(outlist[, col_idx])
- colnames(df_intensities) <- colnames(outlist)[col_idx]
+
+#' Get the indices of the control columns and a dataframe of intensities of the controls
+#'
+#' @param peakgroup_list: Dataframe with intensities for all samples (matrix)
+#' @param label: Label used by grep to get the correct column (string)
+#'
+#' @returns: List with 2 items:
+#' col_idx: vector with indices of the control columns
+#' df_intensities: dataframe with the intensities of the controls
+get_intensities_cols <- function(peakgroup_list, label) {
+ # get the column indices
+ col_idx <- grep(label, colnames(peakgroup_list), fixed = TRUE)
+ # get intensity columns
+ df_intensities <- as.data.frame(peakgroup_list[, col_idx])
+ colnames(df_intensities) <- colnames(peakgroup_list)[col_idx]
+
return(list(col_idx = col_idx, df_intensities = df_intensities))
}
-calculate_zscores <- function(outlist, zscore_type, control_cols, stat_filter, intensity_col_ids, startcol) {
- #' Calculate the Z-scores with different methods for excluding controls
- #'
- #' @param outlist: dataframe with intensities for all samples
- #' @param zscore_type: string with method for excluding controls
- #' @param control_cols: vector with indices of the control columns
- #' @param stat_filter: integer used for excluding controls, either percentage or outlier threshold
- #' @param intensity_col_ids: vector with indices of the samples for which to calculate Z-scores
- #' @param startcol: integer of the column from where to add the Z-score columns
- #'
- #' @returns: outlist: same dataframe as the input with added Z-score columns
+#' Calculate the Z-scores with different methods for excluding outliers in controls
+#'
+#' @param peakgroup_list: dataframe with intensities for all samples (matrix)
+#' @param zscore_type: string with method for excluding controls
+#' @param control_cols: vector with indices of the control columns
+#' @param stat_filter: integer used for excluding controls, either percentage or outlier threshold
+#' @param intensity_col_ids: vector with indices of the samples for which to calculate Z-scores
+#'
+#' @returns: peakgroup_list: same dataframe as the input with added Z-score columns (matrix)
+calculate_zscores <- function(peakgroup_list, zscore_type, control_cols, stat_filter, intensity_col_ids) {
# Calculate mean and sd
- outlist$avg_ctrls <- 0
- outlist$sd_ctrls <- 0
- outlist$nr_ctrls <- length(control_cols)
+ peakgroup_list$avg_ctrls <- 0
+ peakgroup_list$sd_ctrls <- 0
+ peakgroup_list$nr_ctrls <- length(control_cols)
if (zscore_type == "_Zscore") {
# Calculate mean and sd with all controls
- outlist$avg_ctrls <- apply(control_cols, 1, function(x) mean(as.numeric(x), na.rm = TRUE))
- outlist$sd_ctrls <- apply(control_cols, 1, function(x) sd(as.numeric(x), na.rm = TRUE))
+ peakgroup_list$avg_ctrls <- apply(peakgroup_list[, control_cols], 1, function(x) mean(as.numeric(x), na.rm = TRUE))
+ peakgroup_list$sd_ctrls <- apply(peakgroup_list[, control_cols], 1, function(x) sd(as.numeric(x), na.rm = TRUE))
} else {
if (length(control_cols) > 3) {
- for (metabolite_index in seq_len(nrow(outlist))) {
+ for (metabolite_index in seq_len(nrow(peakgroup_list))) {
if (zscore_type == "_RobustZscore") {
# Calculate mean and sd, remove outlier controls by using robust scaler
- outlist$avg_ctrls[metabolite_index] <- mean(robust_scaler(
- outlist[metabolite_index, control_cols],
+ peakgroup_list$avg_ctrls[metabolite_index] <- mean(robust_scaler(
+ peakgroup_list[metabolite_index, control_cols],
control_cols, stat_filter
))
- outlist$sd_ctrls[metabolite_index] <- sd(robust_scaler(
- outlist[metabolite_index, control_cols],
+ peakgroup_list$sd_ctrls[metabolite_index] <- sd(robust_scaler(
+ peakgroup_list[metabolite_index, control_cols],
control_cols, stat_filter
))
} else {
# Calculate mean, sd and number of remaining controls, remove outlier controls by using grubbs test
intensities_without_outliers <- remove_outliers_grubbs(
- as.numeric(outlist[metabolite_index, control_cols]),
+ as.numeric(peakgroup_list[metabolite_index, control_cols]),
stat_filter
)
- outlist$avg_ctrls[metabolite_index] <- mean(intensities_without_outliers)
- outlist$sd_ctrls[metabolite_index] <- sd(intensities_without_outliers)
- outlist$nr_ctrls[metabolite_index] <- length(intensities_without_outliers)
+ peakgroup_list$avg_ctrls[metabolite_index] <- mean(intensities_without_outliers)
+ peakgroup_list$sd_ctrls[metabolite_index] <- sd(intensities_without_outliers)
+ peakgroup_list$nr_ctrls[metabolite_index] <- length(intensities_without_outliers)
}
}
}
}
# Calculate Z-scores
- outlist_zscores <- apply(outlist[, intensity_col_ids, drop = FALSE], 2, function(col) {
- (as.numeric(col) - outlist$avg_ctrls) / outlist$sd_ctrls
+ outlist_zscores <- apply(peakgroup_list[, intensity_col_ids, drop = FALSE], 2, function(col) {
+ (as.numeric(col) - peakgroup_list$avg_ctrls) / peakgroup_list$sd_ctrls
})
- outlist <- cbind(outlist, outlist_zscores)
- colnames(outlist)[startcol:ncol(outlist)] <- paste0(colnames(outlist)[intensity_col_ids], zscore_type)
+ colnames(outlist_zscores) <- paste0(colnames(peakgroup_list)[intensity_col_ids], zscore_type)
+ peakgroup_list <- cbind(peakgroup_list, outlist_zscores)
- return(outlist)
+ return(peakgroup_list)
}
+#' Robust scaler: remove outlier values in controls
+#'
+#' @param control_intensities: Intensities for control samples (vector of float)
+#' @param control_col_ids: Column names for control samples (vector of string)
+#' @param perc: Percentage of outliers which will be removed from controls (float)
+#'
+#' @returns trimmed_control_intensities: Intensities trimmed for outliers (vector of float)
robust_scaler <- function(control_intensities, control_col_ids, perc = 5) {
- #' Calculate robust scaler: Z-score based on controls without outliers
- #'
- #' @param control_intensities: Matrix with intensities for control samples
- #' @param control_col_ids: Vector with column names for control samples
- #' @param perc: Percentage of outliers which will be removed from controls (float)
- #'
- #' @return trimmed_control_intensities: Intensities trimmed for outliers
+ # determine how many values will be removed, based on percentage
nr_to_remove <- ceiling(length(control_col_ids) * perc / 100)
+ # sort intensities
sorted_control_intensities <- sort(as.numeric(control_intensities))
+ # remove highest and lowest intensities
start_index <- nr_to_remove + 1
end_index <- length(sorted_control_intensities) - nr_to_remove
trimmed_control_intensities <- sorted_control_intensities[start_index:end_index]
+
return(trimmed_control_intensities)
}
+#' Remove outlier values using Grubb's test
+#'
+#' @param control_intensities: Intensities for control samples (vector of float)
+#' @param outlier_threshold: Threshold for outliers to be removed from controls (float)
+#'
+#' @returns trimmed_control_intensities: Intensities trimmed for outliers (vector of float)
remove_outliers_grubbs <- function(control_intensities, outlier_threshold = 2) {
- #' Remove outliers per metabolite according to Grubb's test
- #'
- #' @param control_intensities: Vector with intensities for control samples
- #' @param outlier_threshold: Threshold for outliers which will be removed from controls (float)
- #'
- #' @return trimmed_control_intensities: Intensities trimmed for outliers
+ # calculate Z-scores for all intensities
mean_permetabolite <- mean(as.numeric(control_intensities))
stdev_permetabolite <- sd(as.numeric(control_intensities))
zscores_permetabolite <- (control_intensities - mean_permetabolite) / stdev_permetabolite
- # remove intensities with a zscore_permetabolite greater than outlier_threshold
+ # remove intensities with a Z-score greater than outlier_threshold
if (sum(zscores_permetabolite > outlier_threshold) > 0) {
trimmed_control_intensities <- control_intensities[-which(zscores_permetabolite > outlier_threshold)]
} else {
trimmed_control_intensities <- control_intensities
}
+
return(trimmed_control_intensities)
}
+#' Save a dataframe to RData and txt
+#'
+#' @param df: Dataframe (matrix)
+#' @param file_name: File name (string)
save_to_rdata_and_txt <- function(df, file_name) {
- #' Save a dataframe to RData and txt
- #'
- #' @param df: dataframe
- #' @param file_name: string with the file name
save(df, file = paste0(file_name, ".RData"))
write.table(df, file = paste0(file_name, ".txt"), sep = "\t", row.names = FALSE)
}
+#' Set the row height and column width of the Excel
+#'
+#' @param wb: An openxlsx workbook (S4 object)
+#' @param sheetname: Name of the workbook sheet (string)
+#' @param num_rows_df: Number of rows in dataframe (integer)
+#' @param num_col_df: Number of columns in dataframe (integer)
+#' @param plot_width: Width of the plots to be added (integer)
+#' @param plots_present: Parameter that indicates whether plots are added to the workbook (boolean)
+#'
+#' @returns wb: Workbook object with changed row heights and column widths
set_row_height_col_width_wb <- function(wb, sheetname, num_rows_df, num_cols_df, plot_width, plots_present) {
- #' Change the row height and column width of the Excel
- #'
- #' @param wb: an openxlsx workbook (S4 object)
- #' @param sheetname: name of the workbook sheet (string)
- #' @param num_rows_df: number of rows in dataframe (int)
- #' @param num_col_df: number of columns in dataframe (int)
- #' @param plot_width: width of the plots to be added (int)
- #' @param plots_present: boolean if plots are added to the workbook (boolean)
- #'
- #' @returns wb: a workbook object with changed row heights and column widths
if (plots_present) {
openxlsx::setColWidths(wb, sheetname, cols = 1, widths = plot_width / 20)
openxlsx::setRowHeights(wb, sheetname, rows = c(seq(2, num_rows_df + 1)), heights = 560 / 4)
@@ -144,9 +153,10 @@ set_row_height_col_width_wb <- function(wb, sheetname, num_rows_df, num_cols_df,
#' pivot to long format, arrange Samples nummerically, change Sample names, get group size and
#' set Intensities to numeric.
#'
-#' @param intensities_df: a dataframe with HMDB_key column and intensities for all samples
+#' @param intensities_df: Dataframe with HMDB_key column and intensities (matrix)
+#' @param row_index: Index of row (integer)
#'
-#' @returns intensities_df_long: a dataframe with on each row a sample and their intensity
+#' @returns intensities_df_long: Dataframe with on each row a sample and its intensity (matrix)
intensities_df_to_long_format <- function(intensities_df, row_index) {
intensities_df_long <- intensities_df %>%
slice(row_index) %>%
@@ -170,8 +180,8 @@ intensities_df_to_long_format <- function(intensities_df, row_index) {
#' Create a plot of intensities of samples for Excel
#' Use boxplot if group size is above 2, otherwise use a dash/line
#'
-#' @param intensities_df_long: a dataframe with on each row a sample and their intensity
-#' @param hmdb_id: HMDB ID of the selected metabolite
+#' @param intensities_df_long: Dataframe with on each row a sample and its intensity (matrix)
+#' @param hmdb_id: HMDB ID of the selected metabolite (string)
#'
#' @returns boxplot_object: ggplot2 object containing the plot of intensities
create_boxplot <- function(intensities_df_long, hmdb_id) {
@@ -198,18 +208,17 @@ create_boxplot <- function(intensities_df_long, hmdb_id) {
return(boxplot_object)
}
-#' Make and save a boxplot of intensities to an Excel workbook
-#'
-#' For the Helix Excel the positive controls and SST mix samples are removed.
+#' Make and save a boxplot of intensities in png format and insert into an Excel workbook
+#' For the Helix Excel, the positive controls and SST mix samples are removed.
#'
-#' @param excel_workbook: an openxlsx Workbook object
-#' @param sheetname: a string containing the sheetname where the plots are to be placed
-#' @param intensities_df: a dataframe containing intensities for controls and patients of a specific HMDB ID
-#' @param file_path: a string containing the filepath for the png
-#' @param hmdb_id: a string containing the HMDB ID that the intensities_df contains data for
-#' @param plot_width: an integer containing the plot width for the png
-#' @param col_width: an integer containing the width of the column that has the plots
-#' @param start_row_index: an integer containing the index of the row where the plot has to be placed
+#' @param excel_workbook: An openxlsx Workbook object (workbook object)
+#' @param sheetname: Name of the sheet where the plots are to be placed (string)
+#' @param intensities_df: Dataframe containing intensities for controls and patients of a specific HMDB ID (matrix)
+#' @param file_path: Filepath for the png (string)
+#' @param hmdb_id: HMDB ID corresponding to intensities_df (string)
+#' @param plot_width: Plot width for the png (integer)
+#' @param col_width: Width of the column that has the plots (integer)
+#' @param start_row_index: Index of the row where the plot has to be placed (integer)
save_plot_to_excel_workbook <- function(excel_workbook,
sheetname,
intensities_df,
diff --git a/DIMS/export/generate_qc_output_functions.R b/DIMS/export/generate_qc_output_functions.R
index 9f672adb..759e9a01 100644
--- a/DIMS/export/generate_qc_output_functions.R
+++ b/DIMS/export/generate_qc_output_functions.R
@@ -1,32 +1,33 @@
-check_number_of_controls <- function(outlist, min_num_controls, file_name) {
- #' Check te number of controls for all metabolites and report findings in a txt file
- #'
- #' @param outlist: Dataframe with intensities and Z-scores for all samples and controls
- #' @param min_num_controls: Integer that is the minimum number of controls
- #' @param file_name: String containing the file name
-
- outlist_under_ctrls <- outlist %>%
+# functions for checking the quality of the dataset (internal standards, positive controls and SST sample)
+
+#' Check the number of control samples for each metabolite and report findings in a txt file
+#'
+#' @param peakgroup_list: Dataframe with intensities and Z-scores for all samples (matrix)
+#' @param min_num_controls: Minimum number of controls (integer)
+#' @param file_name: Output file name (string)
+check_number_of_controls <- function(peakgroup_list, min_num_controls, file_name) {
+ peakgroups_below_nr_control <- peakgroup_list %>%
filter(nr_ctrls < min_num_controls) %>%
select(HMDB_name, HMDB_code)
- if (nrow(outlist_under_ctrls) == 0) {
+ if (nrow(peakgroups_below_nr_control) == 0) {
writeLines(paste0("All metabolites have", min_num_controls, "or more controls."), file_name)
} else {
writeLines(paste0("These metabolites have less than", min_num_controls, "controls. \n"), file_name)
- write.table(outlist_under_ctrls, file_name, append = TRUE, row.names = FALSE, col.names = TRUE, sep = "\t")
+ write.table(peakgroups_below_nr_control, file_name, append = TRUE, row.names = FALSE, col.names = TRUE, sep = "\t")
}
}
+#' Get the internal standards data
+#'
+#' @param internal_stand_df: dataframe with all data for all the internal standards (matrix)
+#' @param scanmode: positive, negative or summed (string)
+#' @param is_subset_filter: filter/threshold for outlier control removal (float)
+#' @param dims_matrix: matrix used, e.g. Plasma, Research, etc. (string)
+#' @param rundate: date of pipeline run (Date object)
+#' @param project: project name (string)
+#'
+#' @returns internal_stand: dataframe with the intensities of the internal standards for all samples (matrix)
get_internal_standards <- function(internal_stand_df, scanmode, is_subset_filter, dims_matrix, rundate, project) {
- #' Get the internal standards data
- #'
- #' @param internal_stand_df: dataframe with all data for all the internal standards
- #' @param scanmode: positive, negative or summed (string)
- #' @param is_subset_filter: filter/threshold for outlier control removal (float)
- #' @param dims_matrix: matrix used, e.g. Plasma, Research, etc. (string)
- #' @param rundate: date of pipeline run (Date object)
- #' @param project: project name (string)
- #'
- #' @returns internal_stand: dataframe with the intensity of the internal standards for all samples
if (scanmode == "summed") {
internal_stand <- internal_stand_df[c(names(is_subset_filter), "HMDB_code")]
internal_stand$HMDB_name <- internal_stand_df$name
@@ -49,26 +50,23 @@ get_internal_standards <- function(internal_stand_df, scanmode, is_subset_filter
return(internal_stand)
}
+#' Generate and save internal standard plot
+#'
+#' @param plot_data: Dataframe with the data to be plotted (matrix)
+#' @param plot_type: Type of plot (string)
+#' @param plot_title: Title for the plot (string)
+#' @param file_name: Name of the file (string)
+#' @param plot_width: Width of the plot (int)
+#' @param plot_height: Height of the plot (int)
+#' @param hline_data: Dataframe with values for the minimal intensity line (matrix)
save_internal_standard_plot <- function(
plot_data,
plot_type,
plot_title,
- outdir,
file_name,
plot_width,
plot_height,
hline_data = NULL) {
- #' Generate and save internal standard plot
- #'
- #' @param plot_data: dataframe with the data to be plotted
- #' @param plot_type: type of plot (string)
- #' @param plot_title: title for the plot (string)
- #' @param outdir: directory where the plot needs to be saved (string)
- #' @param file_name: name of the file (string)
- #' @param plot_width: width of the plot (int)
- #' @param plot_height: height of the plot (int)
- #' @param hline_data: values for the minimal intensity line (dataframe)
- #'
# Check if plot_data contains data
if (nrow(plot_data) == 0) {
@@ -116,21 +114,21 @@ save_internal_standard_plot <- function(
)
}
- ggplot2::ggsave(paste0(outdir, "/plots/", file_name, ".png"),
+ ggplot2::ggsave(paste0("./plots/", file_name, ".png"),
plot = plot, height = plot_height, width = plot_width, units = "in"
)
}
-get_pos_ctrl_data <- function(outlist, sample_name, hmdb_codes, hmdb_names) {
- #' Get the positive control data
- #'
- #' @param outlist: dataframe with intensities for all samples
- #' @param sample_name: positive control sample name(s) (string)
- #' @param hmdb_codes: HMDB codes of the positive control metabolites (vector)
- #' @param hmdb_names: HMDB names of the positive control metabolites (vector)
- #'
- #' @returns: pos_ctrl_data: dataframe with intensities and Z-scores of the positive control metabolites
- pos_ctrl_data <- outlist[hmdb_codes, c("HMDB_code", "name", sample_name)]
+#' Get the positive control data
+#'
+#' @param peakgroup_list: Dataframe with intensities for all samples (matrix)
+#' @param sample_name: positive control sample name(s) (string)
+#' @param hmdb_codes: HMDB codes of the positive control metabolites (vector)
+#' @param hmdb_names: HMDB names of the positive control metabolites (vector)
+#'
+#' @returns: pos_ctrl_data: Dataframe with intensities and Z-scores of the positive control metabolites (matrix)
+get_pos_ctrl_data <- function(peakgroup_list, sample_name, hmdb_codes, hmdb_names) {
+ pos_ctrl_data <- peakgroup_list[hmdb_codes, c("HMDB_code", "name", sample_name)]
pos_ctrl_data <- reshape2::melt(pos_ctrl_data, id.vars = c("HMDB_code", "name"))
colnames(pos_ctrl_data) <- c("HMDB_code", "HMDB_name", "Sample", "Zscore")
pos_ctrl_data$HMDB_name <- hmdb_names
@@ -140,26 +138,27 @@ get_pos_ctrl_data <- function(outlist, sample_name, hmdb_codes, hmdb_names) {
return(pos_ctrl_data)
}
+#' Round numbers to a set number of digits for numeric values
+#'
+#' @param df: Dataframe containing numeric values (matrix)
+#' @param digits: Number of digits to round off to (integer)
+#'
+#' @return df: Dataframe with rounded numbers (matrix)
round_df <- function(df, digits) {
- #' Round numbers to a set number of digits for numeric values
- #'
- #' @param df: Dataframe containing numeric values
- #' @param digits: Number of digits to round off to (integer)
- #'
- #' @return df: Dataframe with rounded numbers
numeric_columns <- sapply(df, mode) == "numeric"
df[numeric_columns] <- round(df[numeric_columns], digits)
+
return(df)
}
+#' Get internal standard intensities
+#'
+#' @param is_data: Dataframe with intensities of internal standards (matrix)
+#' @param int_cols: Indices of internal standard columns (vector of integers)
+#' @param is_codes: Internal standard codes (vector of strings)
+#'
+#' @returns: is_intensities: Dataframe with intensities of internal standards (matrix)
get_is_intensities <- function(is_data, int_cols = NULL, is_codes = NULL) {
- #' Get internal standard intensities
- #'
- #' @param is_data: dataframe with intensities of internal standards
- #' @param int_cols: default = NULL, if present indices of internal standard columns
- #' @param is_codes: default = NULL, if present internal standard codes
- #'
- #' @returns: is_intensities: dataframe with intensities of internal standards
if (is.null(is_codes)) {
is_intensities <- is_data[, int_cols]
} else {
@@ -168,15 +167,16 @@ get_is_intensities <- function(is_data, int_cols = NULL, is_codes = NULL) {
}
is_intensities <- calc_coefficient_of_variation(is_intensities)
is_intensities <- cbind(IS_name = is_data$HMDB_name, is_intensities)
+
return(is_intensities)
}
+#' Calculate coefficent of variation (cv) based on standard deviation (sd) and mean
+#'
+#' @param intensity_list: Intensities (matrix)
+#'
+#' @return intensity_list_with_cv: Intensities and cv, mean, sd (matrix)
calc_coefficient_of_variation <- function(intensity_list) {
- #' Calculate coefficent of variation (cv) based on standard deviation (sd) and mean
- #'
- #' @param intensity_list: Matrix with intensities
- #'
- #' @return intensity_list_with_cv: Matrix with intensities and cv, mean, sd
intensity_list <- as.data.frame(
apply(intensity_list, 2, function(x) round(as.numeric(x), 0)),
row.names = rownames(intensity_list)
@@ -190,24 +190,23 @@ calc_coefficient_of_variation <- function(intensity_list) {
sd = sd_allsamples,
intensity_list
)
+
return(intensity_list_with_cv)
}
+#' Check if all m/z values are present
+#'
+#' @param mzmed_pgrp_ident: All m/z values for a specific scanmode (vector of float)
+#' @param scanmode: Scan mode, positive or negative (string)
+#'
+#' @return results_mz_missing: Either the missing mz values or a message that no values are missing (string)
check_missing_mz <- function(mzmed_pgrp_ident, scanmode) {
- #' Check if all m/z values are present
- #'
- #' @param mzmed_pgrp_ident: Vector of all m/z values for a specific scanmode
- #' @param scanmode: String with the scanmode, positive or negative
- #'
- #' @return results_mz_missing: String with either the missing mz values or
- #' message that no values are missing.
-
# retrieve all unique m/z values in whole numbers and check if all are available
mzmed_pgrp_ident <- unique(round(mzmed_pgrp_ident, digits = 0))
# m/z range for a standard run = 70-600
mz_range <- seq(70, 599, by = 1)
mz_missing <- setdiff(mz_range, mzmed_pgrp_ident)
- # check if m/z are missing and make an .txt file with information
+ # check if m/z are missing and make a txt file with information
mz_missing_group <- cumsum(c(1, diff(mz_missing) != 1))
if (length(mz_missing_group) > 1) {
results_mz_missing <- c(paste0("Missing m/z values ", scanmode, " mode"))
@@ -218,28 +217,32 @@ check_missing_mz <- function(mzmed_pgrp_ident, scanmode) {
return(results_mz_missing)
}
+#' Create a list of all internal standards with intensity below a threshold value
+#'
+#' @param is_selection_subset: Intensities for each internal standard in each sample (matrix)
+#' @param thresholds: Threshold values for a given scan mode and matrix (vector of integers)
+#' @param is_names: Names of internal standards for a given scan mode (vector of strings)
+#' @param scanmode: Scan mode to include in output (string)
+#'
+#' @return is_below_threshold: All samples for which internal standard intensity is below threshold (matrix)
find_is_below_threshold <- function(is_selection_subset, thresholds, is_names, scanmode) {
- #' Create a list of all internal standards with intensity below a threshold value
- #'
- #' @param is_selection_subset: Matrix with intensities for each internal standard in each sample
- #' @param thresholds: Threshold values for a given scan mode and matrix
- #' @param is_names: Array of names of internal standards for a given scan mode
- #' @param scanmode: string indicating scan mode to include in output
- #'
- #' @return is_below_threshold: Matrix listing all samples for which internal standard intensity is below threshold
-
- # initialize; get the headers of the matrix
- is_below_threshold <- is_selection_subset[0, ]
# for every line, check if intensity is below the appropriate threshold
+ below_threshold_index <- c()
for (line_index in seq_len(nrow(is_selection_subset))) {
is_selected <- is_selection_subset$HMDB_name[line_index]
thresh_selected <- thresholds[which(is_names == is_selected)]
if (is_selection_subset$Intensity[line_index] < thresh_selected) {
- is_below_threshold <- rbind(is_below_threshold, is_selection_subset[line_index, ])
+ below_threshold_index <- c(below_threshold_index, line_index)
}
}
+ is_below_threshold <- is_selection_subset[below_threshold_index, ]
# add information on scan mode
- is_below_threshold <- cbind(is_below_threshold, scanmode = rep(scanmode, nrow(is_below_threshold)))
+ if (nrow(is_below_threshold) > 0) {
+ is_below_threshold <- cbind(is_below_threshold, scanmode = scanmode)
+ } else {
+ is_below_threshold$scanmode_char <- character(0)
+ }
+
return(is_below_threshold)
}
diff --git a/DIMS/export/generate_violin_plots_functions.R b/DIMS/export/generate_violin_plots_functions.R
index c314ce72..32c0b33d 100644
--- a/DIMS/export/generate_violin_plots_functions.R
+++ b/DIMS/export/generate_violin_plots_functions.R
@@ -1,12 +1,15 @@
+# functions for creating violin plot PDF output
+
#' Preparing the intensities and Z-score dataframe.
#' Certain columns are removed, the HMDB_code and HMDB_name column are moved forward,
#' the avg_ctrls and sd_ctrls columns are renamed and the column type of all columns containing numbers
#' is changed to numeric.
#'
-#' @param intensities_zscore_df: dataframe with intensities, Z-scores and metabolite information for all samples
+#' @param intensities_zscore_df: Dataframe with intensities, Z-scores and metabolite information for all
+#' samples (matrix)
#'
-#' @returns intensities_zscore_df: a dataframe containing intensities, Z-scores, HMDB IDs, HMDB names and
-#' the mean and average of all controls
+#' @return intensities_zscore_df: Dataframe containing intensities, Z-scores, HMDB IDs, HMDB names and
+#' renamed columns for mean and average of all controls
prepare_intensities_zscore_df <- function(intensities_zscore_df) {
intensities_zscore_df <- intensities_zscore_df %>%
select(-c(
@@ -16,41 +19,44 @@ prepare_intensities_zscore_df <- function(intensities_zscore_df) {
relocate(c(HMDB_code, HMDB_name)) %>%
rename(mean_controls = avg_ctrls, sd_controls = sd_ctrls) %>%
mutate(across(!c(HMDB_name, HMDB_code), as.numeric))
+
return(intensities_zscore_df)
}
-#' Get all column names containing a specific prefix.
+#' Get all column names containing a specific prefix
#'
-#' @param dataframe: dataframe containing multiple columns with Z-scores
-#' @param prefix: a string of a prefix to be searched in the column names, e.g. "P" or "C".
+#' @param dataframe: Dataframe containing multiple columns with Z-scores (matrix)
+#' @param prefix: Prefix to be searched in the column names, e.g. "P" or "C" (string)
#'
-#' @returns sample_colnames: a vector of column names all containing the prefix.
+#' @return sample_colnames: Column names all containing the prefix (vector of strings)
get_colnames_by_prefix <- function(dataframe, prefix) {
sample_colnames <- grep(paste0("^", prefix), colnames(dataframe), value = TRUE)
+
return(sample_colnames)
}
#' Remove the suffix from a vector of names
#'
-#' @param vector_names: vector containing names with or without a suffix
-#' @param suffix: string containing the suffix to be removed
+#' @param vector_names: Names with or without a suffix (vector of strings)
+#' @param suffix: Suffix to be removed (string)
#'
-#' @returns names_no_suffix: a vector of unique names without the suffix
+#' @return names_no_suffix: Unique names without the suffix (vector of strings)
remove_suffix_from_items <- function(vector_names, suffix) {
names_no_suffix <- unique(gsub("_Zscore", "", vector_names))
+
return(names_no_suffix)
}
#' Add Zscores for multiple ratios to the dataframe
#'
-#' @param outlist: dataframe containing intensities and Z-scores for all controls and patients
-#' @param metabolites_ratios_df: dataframe containing numerators and denominators for all ratios
-#' @param all_sample_ids: vector of sample IDS, controls and patients
+#' @param peakgroup_list: Dataframe containing intensities and Z-scores for all controls and patients (matrix)
+#' @param metabolites_ratios_df: Dataframe containing numerators and denominators for all ratios (matrix)
+#' @param all_sample_ids: Sample names for controls and patients (vector of strings)
#'
-#' @returns intensities_zscore_ratios_df: dataframe containing intensities and Z-scores for all controls and patients
-#' for all metabolites and ratios
-add_zscores_ratios_to_df <- function(outlist, metabolites_ratios_df, all_sample_ids) {
- intensities_zscores_df <- prepare_intensities_zscore_df(outlist)
+#' @return intensities_zscore_ratios_df: Dataframe containing intensities and Z-scores for all controls and patients
+#' for all metabolites and ratios (matrix)
+add_zscores_ratios_to_df <- function(peakgroup_list, metabolites_ratios_df, all_sample_ids) {
+ intensities_zscores_df <- prepare_intensities_zscore_df(peakgroup_list)
# calculate Z-scores for the ratios
zscore_ratios_df <- calculate_zscore_ratios(metabolites_ratios_df, intensities_zscores_df, all_sample_ids)
@@ -61,11 +67,11 @@ add_zscores_ratios_to_df <- function(outlist, metabolites_ratios_df, all_sample_
#' Calculate Z-scores for ratios
#'
-#' @param metabolites_ratios_df: dataframe containing numerators and denominators for all ratios
-#' @param intensities_zscores_df: dataframe containing intensities and Z-scores for all controls and patients
-#' @param intensity_col_names: vector of sample IDS, controls and patients
+#' @param metabolites_ratios_df: Dataframe containing numerators and denominators for all ratios (matrix)
+#' @param intensities_zscores_df: Dataframe containing intensities and Z-scores for all samples (matrix)
+#' @param intensity_col_names: Sample names (vector of strings)
#'
-#' @returns zscore_ratios_df: dataframe containing Z-scores for all ratios for all samples
+#' @return zscore_ratios_df: Dataframe containing Z-scores for all ratios for all samples (matrix)
calculate_zscore_ratios <- function(metabolites_ratios_df, intensities_zscores_df, intensity_col_names) {
# remove Z-score columns from intensity_col_names
if (any(grepl("_Zscore", intensity_col_names))) {
@@ -102,7 +108,7 @@ calculate_zscore_ratios <- function(metabolites_ratios_df, intensities_zscores_d
intensity_col_names
)
# calculate the intensity ratio for each sample
- zscore_ratios_df[row_index, intensity_cols_index] <- log2(numerator_intensities / denominator_intensities)
+ zscore_ratios_df[row_index, intensity_cols_index] <- numerator_intensities / denominator_intensities
}
control_intensities_cols_index <- grep("^C[^_]*$", colnames(intensities_zscores_df), perl = TRUE)
@@ -123,15 +129,15 @@ calculate_zscore_ratios <- function(metabolites_ratios_df, intensities_zscores_d
#' Make and save violin plots for each patient in a PDF
#'
-#' @param zscore_patients_df: dataframe with Z-scores for all patient samples
-#' @param zscore_controls_df: dataframe with Z-scores for all control samples
-#' @param path_metabolite_groups: string containing the path for the metabolite groups directories
-#' @param nr_plots_perpage: integer containing the number of metabolites on a plot per page
-#' @param number_of_samples: list containing the number of patient and control samples
-#' @param run_name: string containing the run name
-#' @param protocol_name: string containing the protocol name
-#' @param explanation_violin_plot: vector of strings containing the explanation of the violin plots
-#' @param number_of_metabolites: list containing the number of metabolites for the top and lowest table
+#' @param zscore_patients_df: Dataframe with Z-scores for all patient samples (matrix)
+#' @param zscore_controls_df: Dataframe with Z-scores for all control samples (matrix)
+#' @param path_metabolite_groups: Path for the metabolite groups directories (string)
+#' @param nr_plots_perpage: Number of metabolites on a plot per page (integer)
+#' @param number_of_samples: Number of patient and control samples (list of integers)
+#' @param run_name: Run name (string)
+#' @param protocol_name: Protocol name (string)
+#' @param explanation_violin_plot: Explanation of the violin plots (vector of strings)
+#' @param number_of_metabolites: Number of metabolites for the top and lowest table (list of integers)
make_and_save_violin_plot_pdfs <- function(
zscore_patients_df,
zscore_controls_df,
@@ -170,6 +176,8 @@ make_and_save_violin_plot_pdfs <- function(
if (any(is_diagnostic_patients(dims_helix_table$Sample))) {
# transform dataframe for Helix output
output_helix <- transform_metab_df_to_helix_df(protocol_name, dims_helix_table)
+ # round the Z-scores to 1 decimal
+ output_helix$Amount <- round(output_helix$Amount, 1)
# save the DIMS Helix dataframe
path_helixfile <- paste0("./output_Helix_", run_name, ".csv")
write.csv(output_helix, path_helixfile, quote = FALSE, row.names = FALSE)
@@ -196,11 +204,11 @@ make_and_save_violin_plot_pdfs <- function(
}
}
-#' Get a list with dataframes for all off the metabolite group in a directory
+#' Get a list with dataframes for all of the metabolite group in a directory
#'
-#' @param dir_with_subdirs: directory containing txt files with metabolites per group (string)
+#' @param dir_with_subdirs: Directory containing txt files with metabolites per group (string)
#'
-#' @returns list_of_dataframes: list with dataframes with info on metabolites (list of dataframes)
+#' @return list_of_dataframes: Dataframes with info on metabolites (list of dataframes)
get_list_dataframes_from_dir <- function(dir_with_subdirs) {
# get a list of all metabolite files
txt_files_paths <- list.files(dir_with_subdirs, pattern = "*.txt", recursive = FALSE, full.names = TRUE)
@@ -213,11 +221,11 @@ get_list_dataframes_from_dir <- function(dir_with_subdirs) {
#' Merge patient Z-scores with metabolite info
#'
-#' @param list_df_metabolite_groups: list of dataframes with metabolite information for different metabolite classes (list)
-#' @param zscore_df: dataframe with metabolite Z-scores for all patient
+#' @param list_df_metabolite_groups: Dataframes with metabolite information for different metabolite classes (list)
+#' @param zscore_df: Dataframe with metabolite Z-scores for all patient (matrix)
#'
-#' @return list_dfs_metabs_info_zscores: list of dataframes for each metabolite class
-#' containing info and zscores for all samples
+#' @return list_dfs_metabs_info_zscores: Dataframes for each metabolite class
+#' containing info and zscores for all samples (list)
merge_metabolite_info_zscores <- function(list_df_metabolite_groups, zscore_df) {
# remove HMDB_name column and "_Zscore" from column (patient) names
zscore_df <- zscore_df %>%
@@ -252,16 +260,16 @@ merge_metabolite_info_zscores <- function(list_df_metabolite_groups, zscore_df)
return(list_dfs_metabs_info_zscores)
}
-#' Combine patient and control data for each page of the violinplot pdf
+#' Combine patient and control data for each page of the violin plot pdf
#'
-#' @param metab_interest_patients: list of dataframes with data for each metabolite and patient (list)
-#' @param metab_interest_controls: list of dataframes with data for each metabolite and control (list)
-#' @param number_of_plots_per_page: number of plots per page in the violinplot pdf (integer)
-#' @param number_of_patients: number of patients (integer)
-#' @param number_of_controls: number of controls (integer)
+#' @param metab_interest_patients: Dataframes with data for each metabolite and patient (list)
+#' @param metab_interest_controls: Dataframes with data for each metabolite and control (list)
+#' @param number_of_plots_per_page: Number of plots per page in the violinplot pdf (integer)
+#' @param number_of_patients: Number of patients (integer)
+#' @param number_of_controls: Number of controls (integer)
#'
-#' @return list_metabolite_df_per_page: list of dataframes with metabolite Z-scores for each patient and control,
-#' the length of list is the number of pages for the violinplot pdf (list)
+#' @return list_metabolite_df_per_page: Dataframes with metabolite Z-scores for each patient and control,
+#' the length of list is the number of pages for the violin plot pdf (list)
get_data_per_metabolite_class <- function(
metab_interest_patients,
metab_interest_controls,
@@ -306,7 +314,7 @@ get_data_per_metabolite_class <- function(
#' @param list_metabolite_classes: list of tables with metabolites for Helix and violin plots (list)
#'
#' @return df_zscores_to_helix: dataframe with patient data with only metabolites for Helix and violin plots
-#' with Helix name, high/low Z-score cutoffs
+#' with Helix name, high/low Z-score cutoffs
prepare_helix_patient_data <- function(list_dfs_metab_classes_zscores, list_metabolite_classes) {
# Combine Z-scores of metab groups together
metabolite_zscore_dataframe <- bind_rows(list_dfs_metab_classes_zscores)
@@ -340,13 +348,13 @@ prepare_helix_patient_data <- function(list_dfs_metab_classes_zscores, list_meta
#' Getting the intensities for calculating ratio Z-scores
#' Retrieving a vector of intensities for a particular fraction side of the ratios for all samples.
#'
-#' @param ratios_metabs_df: dataframe with HMDB codes for the ratios (dataframe)
-#' @param row_index: index of the row in the ratios_metabs_df (integer)
-#' @param intensities_zscore_df: dataframe with intensities for each sample (dataframe)
-#' @param fraction_side: either numerator or denominator, which side of the fraction (string)
-#' @param intensity_cols: names of the columns that contain the intensities (string)
+#' @param ratios_metabs_df: Dataframe with HMDB codes for the ratios (matrix)
+#' @param row_index: Index of the row in the ratios_metabs_df (integer)
+#' @param intensities_zscore_df: Dataframe with intensities for each sample (matrix)
+#' @param fraction_side: Either numerator or denominator, which side of the fraction (string)
+#' @param intensity_cols: Names of the columns that contain the intensities (vector of strings)
#'
-#' @returns fraction_side_intensity: a vector of intensities (vector of integers)
+#' @return fraction_side_intensity: Intensities (vector of floats)
get_intensities_fraction_side <- function(ratios_metabs_df, row_index, intensities_zscore_df, fraction_side, intensity_cols) {
# get the HMDB ID(s) for the given fraction side
fraction_side_hmdb_ids <- ratios_metabs_df[row_index, fraction_side]
@@ -370,31 +378,33 @@ get_intensities_fraction_side <- function(ratios_metabs_df, row_index, intensiti
}
# vector of intensities for all samples
fraction_side_intensity <- as.numeric(fraction_side_intensity)
+
return(fraction_side_intensity)
}
#' Get the sample IDs for columns that have Z-score and intensities
#'
-#' @param colnames_zscore_cols: vector of sample IDs from the dataframe containing Z-scores (vector of strings)
-#' @param colnames_intensity_cols: vector of sample IDs form the dataframe containing intensities (vector of strings)
+#' @param colnames_zscore_cols: Column names from the dataframe containing Z-scores (vector of strings)
+#' @param colnames_intensity_cols: Sample names from the dataframe containing intensities (vector of strings)
#'
-#' @returns colnames_intersect: vector of sample IDs that are in both input vectors, ending on "_Zscore" (vector of strings)
+#' @return colnames_intersect: Sample names that are in both input vectors, ending on "_Zscore" (vector of strings)
get_sample_ids_with_zscores <- function(colnames_zscore_cols, colnames_intensity_cols) {
colnames_intersect <- intersect(
paste0(colnames_intensity_cols, "_Zscore"),
grep("_Zscore", colnames_zscore_cols, value = TRUE)
)
+
return(colnames_intersect)
}
#' Pad or truncate HMDB names to a fixed width
#' Add spaces or remove HMDB name characters till the length of the name equals the 'width'
#'
-#' @param metabolite_info_df: A dataframe containing a column `HMDB_name` (character).
-#' @param width: Integer target width for the display names. Default is 45.
-#' @param pad_character: Single character used for padding. Default is a space `" "`.
+#' @param metabolite_info_df: A dataframe containing columns `HMDB_code` and `HMDB_name` (matrix)
+#' @param width: Target width for the display names. Default is 45 (integer)
+#' @param pad_character: Single character used for padding. Default is a space `" "` (string)
#'
-#' @return metabolite_info_df: A dataframe where the HMDB names are transformed
+#' @return metabolite_info_df: Dataframe where the HMDB names are transformed (matrix)
pad_truncate_hmdb_names <- function(metabolite_info_df, width, pad_character) {
# Change the HMDB_name column so all names have 45 characters
# remove characters if name is longer and add "..."
@@ -409,15 +419,15 @@ pad_truncate_hmdb_names <- function(metabolite_info_df, width, pad_character) {
return(metabolite_info_df)
}
-#' Get a list of dataframes for each chunk
+#' Get a list of dataframes for each chunk of list of metabolites
#' For each chunk, get a dataframe containing the metabolites in that chunk and add it to the list
#'
-#' @param metabolites_in_chunks: list of vectors, each containing metabolites
-#' @param metabolite_class_patients_df: dataframe of Z-scores for all patient
-#' @param metabolite_class_controls_df: dataframe of Z-scores for all control
-#' @param number_of_plots_per_page: integer containing the number of metabolites per plot per page
+#' @param metabolites_in_chunks: List of metabolites (list of vectors of strings)
+#' @param metabolite_class_patients_df: Dataframe of Z-scores for all patients (matrix)
+#' @param metabolite_class_controls_df: Dataframe of Z-scores for all controls (matrix)
+#' @param number_of_plots_per_page: Number of metabolites per plot per page (integer)
#'
-#' @returns page_plot_data_list: a list of dataframes containing Z-scores
+#' @return patients_controls_df_chunk: List of dataframes containing Z-scores (list of matrices)
get_list_page_plot_data <- function(
metabolites_in_chunks,
metabolite_class_patients_df,
@@ -443,10 +453,10 @@ get_list_page_plot_data <- function(
#' Create the order of metabolites and add empty strings if the number of metabolites is lower than
#' the number of plots per page.
#'
-#' @param number_of_plots_per_page: integer containing the number of metabolites per plot per page
-#' @param metabolite_names_chunk: list of vectors, each containing metabolites
+#' @param number_of_plots_per_page: Number of metabolites per plot per page (integer)
+#' @param metabolite_names_chunk: List of metabolites (list of vectors of strings)
#'
-#' @returns metabolite_order: a vector containing all metabolites and possibly empty strings
+#' @return metabolite_order: All metabolites and possibly empty strings (vector of strings)
make_metabolite_order <- function(number_of_plots_per_page, metabolite_names_chunk) {
# Add empty dummy's to extend the number of metabs to the nr_plots_perpage
number_of_plots_missing <- number_of_plots_per_page - length(metabolite_names_chunk)
@@ -456,14 +466,16 @@ make_metabolite_order <- function(number_of_plots_per_page, metabolite_names_chu
} else {
metabolite_order <- metabolite_names_chunk
}
+
return(metabolite_order)
}
-#' Check for Diagnostics patients with correct patient number (e.g. starting with "P2024M")
+#' Check sample name for Diagnostics patients for vector of sample names
+#' (e.g. starting with "PYYYYM")
#'
-#' @param patient_column: a column from dataframe with IDs (character vector)
+#' @param patient_column: Column names from dataframe (vector of strings)
#'
-#' @return: a logical vector with TRUE or FALSE for each element (vector)
+#' @return diagnostic_patients: Vector with TRUE or FALSE for each element (vector of booleans)
is_diagnostic_patients <- function(patient_column) {
diagnostic_patients <- grepl("^P[0-9]{4}M", patient_column)
@@ -472,10 +484,10 @@ is_diagnostic_patients <- function(patient_column) {
#' Get the output dataframe for Helix
#'
-#' @param protocol_name: protocol name (string)
-#' @param df_metabs_helix: dataframe with metabolite Z-scores for patients (dataframe)
+#' @param protocol_name: Protocol name (string)
+#' @param df_metabs_helix: Dataframe with metabolite Z-scores for patients (matrix)
#'
-#' @return: dataframe with patient metabolite Z-scores in correct format for Helix
+#' @return df_metabs_helix: Same dataframe in correct format for Helix (matrix)
transform_metab_df_to_helix_df <- function(protocol_name, df_metabs_helix) {
# Remove positive controls
df_metabs_helix <- df_metabs_helix %>% filter(is_diagnostic_patients(Sample))
@@ -500,7 +512,7 @@ transform_metab_df_to_helix_df <- function(protocol_name, df_metabs_helix) {
df_metabs_helix <- df_metabs_helix %>%
select(c(Vial, labnummer, Onderzoeksnummer, Protocol, Name, Amount))
- # Remove duplicate patient-metabolite combinations ("leucine + isoleucine + allo-isoleucin_Z-score" is added 3 times)
+ # Remove duplicate patient-metabolite combinations
df_metabs_helix <- df_metabs_helix %>%
group_by(Onderzoeksnummer, Name) %>%
distinct() %>%
@@ -511,9 +523,9 @@ transform_metab_df_to_helix_df <- function(protocol_name, df_metabs_helix) {
#' Adding labnummer and Onderzoeksnummer to a dataframe
#'
-#' @param df_metabs_helix: dataframe with patient data to be uploaded to Helix
+#' @param df_metabs_helix: Dataframe with patient data to be uploaded to Helix (matrix)
#'
-#' @return: dataframe with added labnummer and Onderzoeksnummer columns
+#' @return df_metabs_helix: Same dataframe with labnummer and Onderzoeksnummer columns (matrix)
add_lab_id_and_onderzoeksnr <- function(df_metabs_helix) {
# Split patient number into labnummer and Onderzoeksnummer
for (row in seq_len(nrow(df_metabs_helix))) {
@@ -521,15 +533,17 @@ add_lab_id_and_onderzoeksnr <- function(df_metabs_helix) {
labnummer_split <- strsplit(as.character(df_metabs_helix[row, "labnummer"]), "M")[[1]]
df_metabs_helix[row, "Onderzoeksnummer"] <- paste0("MB", labnummer_split[1], "/", labnummer_split[2])
}
+
return(df_metabs_helix)
}
#' Create a dataframe with all metabolites that exceed the min and max Z-score cutoffs
#'
-#' @param patient_name: patient code (string)
-#' @param dims_helix_table: dataframe with metabolite Z-scores for each patient and Helix info (dataframe)
+#' @param patient_name: Patient name (string)
+#' @param dims_helix_table: Dataframe with metabolite Z-scores for each patient and Helix info (matrix)
#'
-#' @return: dataframe with metabolites that exceed the min and max Z-score cutoffs for the selected patient
+#' @return top_metab_patient: Dataframe with metabolites that exceed the min and max Z-score cutoffs
+#' for the selected patient (matrix)
get_top_metabolites_df <- function(patient_name, dims_helix_table) {
# extract data for patient of interest (patient_name)
patient_metabs_helix <- dims_helix_table %>%
@@ -565,12 +579,12 @@ get_top_metabolites_df <- function(patient_name, dims_helix_table) {
#' Create a dataframe with the top 20 highest and top 10 lowest metabolites per patient
#'
-#' @param pt_name: patient code (string)
-#' @param zscore_patients: dataframe with metabolite Z-scores per patient (dataframe)
-#' @param top_highest: the number of metabolites with the highest Z-score to display in the table (numeric)
-#' @param top_lowest: the number of metabolites with the lowest Z-score to display in the table (numeric)
+#' @param patient_id: Patient name (string)
+#' @param zscore_patients: Dataframe with metabolite Z-scores per patient (dataframe)
+#' @param num_of_highest_metabolites: Number of metabolites with highest Z-scores to display in the table (numeric)
+#' @param num_of_lowest_metabolites: Number of metabolites with lowest Z-scores to display in the table (numeric)
#'
-#' @return: dataframe with 30 metabolites and Z-scores (dataframe)
+#' @return top_metab_patient: dataframe with 30 metabolites and Z-scores (dataframe)
prepare_toplist <- function(patient_id, zscore_patients, num_of_highest_metabolites, num_of_lowest_metabolites) {
patient_df <- zscore_patients %>%
select(HMDB_code, HMDB_name, !!sym(patient_id)) %>%
@@ -587,24 +601,24 @@ prepare_toplist <- function(patient_id, zscore_patients, num_of_highest_metaboli
# add lines for increased, decreased
extra_line1 <- c("Increased", "", "")
extra_line2 <- c("Decreased", "", "")
- top_metab_pt <- rbind(extra_line1, patient_df_high, extra_line2, patient_df_low)
+ top_metab_patient <- rbind(extra_line1, patient_df_high, extra_line2, patient_df_low)
# remove row names
- rownames(top_metab_pt) <- NULL
+ rownames(top_metab_patient) <- NULL
# change column names for display
- colnames(top_metab_pt) <- c("HMDB_ID", "Metabolite", "Z-score")
+ colnames(top_metab_patient) <- c("HMDB_ID", "Metabolite", "Z-score")
- return(top_metab_pt)
+ return(top_metab_patient)
}
#' Create a pdf with table with metabolites and violin plots
#'
-#' @param pdf_dir: location where to save the pdf file (string)
-#' @param patient_id: patient id (string)
-#' @param metab_perpage: list of dataframes, each dataframe contains data for a page in de pdf (list)
-#' @param top_metab_pt: dataframe with increased and decreased metabolites for this patient (dataframe)
-#' @param explanation: text that explains the violin plots and the pipeline version (string)
-create_pdf_violin_plots <- function(pdf_dir, patient_id, metab_perpage, top_metab_pt, explanation) {
+#' @param pdf_dir: Location for saving the pdf file (string)
+#' @param patient_id: Patient name (string)
+#' @param metab_perpage: List of dataframes, each dataframe contains data for a page in de pdf (list)
+#' @param top_metab_patient: Dataframe with increased and decreased metabolites for this patient (matrix)
+#' @param explanation_violin_plot: Text that explains the violin plots and the pipeline version (vector of strings)
+create_pdf_violin_plots <- function(pdf_dir, patient_id, metab_perpage, top_metab_patient, explanation_violin_plot) {
# set parameters for plots
plot_height <- 9.6
plot_width <- 6
@@ -618,13 +632,13 @@ create_pdf_violin_plots <- function(pdf_dir, patient_id, metab_perpage, top_meta
# patient plots, create the PDF device
patient_id_sub <- patient_id
suffix <- ""
- if (grepl("Diagnostics", pdf_dir) && is_diagnostic_patients(patient_id)) {
+ if (grepl("Diagnost", pdf_dir) && is_diagnostic_patients(patient_id)) {
prefix <- "MB"
suffix <- "_DIMS_PL_DIAG"
# substitute P and M in P2020M00001 into right format for Helix
patient_id_sub <- gsub("[PM]", "", patient_id)
patient_id_sub <- gsub("\\..*", "", patient_id_sub)
- } else if (grepl("Diagnostics", pdf_dir)) {
+ } else if (grepl("Diagnost", pdf_dir)) {
prefix <- "Dx_"
} else if (grepl("IEM", pdf_dir)) {
prefix <- "IEM_"
@@ -642,9 +656,9 @@ create_pdf_violin_plots <- function(pdf_dir, patient_id, metab_perpage, top_meta
page_headers <- names(metab_perpage)
# put table into PDF file, if not empty
- if (!is.null(dim(top_metab_pt))) {
+ if (!is.null(dim(top_metab_patient))) {
max_rows_per_page <- 35
- total_rows <- nrow(top_metab_pt)
+ total_rows <- nrow(top_metab_patient)
number_of_pages <- ceiling(total_rows / max_rows_per_page)
# get the names and numbers in the table aligned
@@ -656,7 +670,7 @@ create_pdf_violin_plots <- function(pdf_dir, patient_id, metab_perpage, top_meta
for (page in seq(number_of_pages)) {
start_row <- (page - 1) * max_rows_per_page + 1
end_row <- min(page * max_rows_per_page, total_rows)
- page_data <- top_metab_pt[start_row:end_row, ]
+ page_data <- top_metab_patient[start_row:end_row, ]
table_grob <- tableGrob(page_data, theme = table_theme, rows = NULL)
@@ -696,13 +710,13 @@ create_pdf_violin_plots <- function(pdf_dir, patient_id, metab_perpage, top_meta
suppressWarnings(print(ggplot_object))
}
- # add explanation of violin plots, version number etc.
+ # add explanation_violin_plot of violin plots, version number etc.
plot(NA, xlim = c(0, 5), ylim = c(0, 5), bty = "n", xaxt = "n", yaxt = "n", xlab = "", ylab = "")
- if (length(explanation) > 0) {
- text(0.2, 5, explanation[1], pos = 4, cex = 0.8)
- for (line_index in 2:length(explanation)) {
+ if (length(explanation_violin_plot) > 0) {
+ text(0.2, 5, explanation_violin_plot[1], pos = 4, cex = 0.8)
+ for (line_index in 2:length(explanation_violin_plot)) {
text_y_position <- 5 - (line_index * 0.2)
- text(-0.2, text_y_position, explanation[line_index], pos = 4, cex = 0.5)
+ text(-0.2, text_y_position, explanation_violin_plot[line_index], pos = 4, cex = 0.5)
}
}
@@ -712,12 +726,12 @@ create_pdf_violin_plots <- function(pdf_dir, patient_id, metab_perpage, top_meta
#' Create violin plots
#'
-#' @param metab_zscores_df: dataframe with Z-scores for all samples (dataframe)
-#' @param patient_zscore_df: dataframe with Z-scores for the specified patient (dataframe)
-#' @param sub_perpage: subtitle of the page (string)
-#' @param patient_id: the patient id of the selected patient (string)
+#' @param metab_zscores_df: Dataframe with Z-scores for all samples (matrix)
+#' @param patient_zscore_df: Dataframe with Z-scores for the specified patient (matrix)
+#' @param sub_perpage: Subtitle of the page (string)
+#' @param patient_id: Patient id of the selected patient (string)
#'
-#' @returns ggpplot_object: a violin plot of metabolites that highlights the selected patient (ggplot object)
+#' @return ggplot_object: Violin plot of metabolites that highlights the selected patient (ggplot object)
create_violin_plot <- function(metab_zscores_df, patient_zscore_df, sub_perpage, patient_id) {
fontsize <- 1
circlesize <- 0.8
@@ -776,11 +790,11 @@ create_violin_plot <- function(metab_zscores_df, patient_zscore_df, sub_perpage,
#' Run the dIEM algorithm (DOI: 10.3390/ijms21030979)
#'
-#' @param expected_biomarkers_df: dataframe with information for HMDB codes about IEMs (dataframe)
-#' @param zscore_patients: dataframe containing Z-scores for patient (dataframe)
-#' @param sample_cols: vector containing column names with intensities and Z-scores for patients (vector)
+#' @param expected_biomarkers_df: Dataframe with information for HMDB codes about IEMs (matrix)
+#' @param zscore_patients_df: Dataframe containing Z-scores for all patients (matrix)
+#' @param sample_cols: Column names with intensities and Z-scores for patients (vector of strings)
#'
-#' @returns probability_score: a dataframe with probability scores for IEMs for each patient (dataframe)
+#' @return probability_score: a dataframe with probability scores for IEMs for each patient (matrix)
run_diem_algorithm <- function(expected_biomarkers_df, zscore_patients_df, sample_cols) {
# Rank the metabolites for each patient individually
ranking_patients <- zscore_patients_df %>%
@@ -852,9 +866,9 @@ run_diem_algorithm <- function(expected_biomarkers_df, zscore_patients_df, sampl
#' Ranking Z-scores for a patient, separate for positive and negative Z-scores
#'
-#' @param zscore_col: vector with Z-scores for a single patient (vector of integers)
+#' @param zscore_col: Z-scores for a single patient (vector of floats)
#'
-#' @returns ranking: a vector of the ranking of the Z-scores (vector of integers)
+#' @return ranking: Ranking of the Z-scores (vector of integers)
rank_patient_zscores <- function(zscore_col) {
# Create ranking column with default NA values
ranking <- rep(NA_real_, length(zscore_col))
@@ -870,10 +884,10 @@ rank_patient_zscores <- function(zscore_col) {
return(ranking)
}
-#' Save the probability score dataframe as an Excel file
+#' Save the probability score dataframe to an Excel file
#'
-#' @param probability_score: a dataframe containing probability scores for each patient (dataframe)
-#' @param run_name: name of the run, for the file name (string)
+#' @param probability_score: Dataframe containing probability scores for each patient (matrix)
+#' @param run_name: Name of the run, to be included in the file name (string)
save_prob_scores_to_excel <- function(probability_score, run_name) {
# Create conditional formatting for output Excel sheet. Colors according to values.
wb <- createWorkbook()
@@ -889,16 +903,18 @@ save_prob_scores_to_excel <- function(probability_score, run_name) {
#' Make and save dIEM plots
#'
-#' @param diem_probability_score: dataframe with dIEM probability scores
-#' @param patient_col_names: vector containing all patient column names
-#' @param expected_biomarkers_df: dataframe with information for HMDB codes about IEMs
-#' @param zscore_patients_df: dataframe containing Z-scores for all patients
-#' @param zscore_controls_df: dataframe containing Z-scores for all controls
-#' @param nr_plots_perpage: integer containing the number of metabolites per page
-#' @param number_of_samples: list containing the number of patients and controls
-#' @param number_of_metabolites: list containing the number of metabolites for the top and lowest table
+#' @param diem_probability_score: Dataframe with dIEM probability scores (matrix)
+#' @param patient_col_names: All patient column names (vector of strings)
+#' @param expected_biomarkers_df: Dataframe with information for HMDB codes about IEMs (matrix)
+#' @param zscore_patients_df: Dataframe containing Z-scores for all patients (matrix)
+#' @param zscore_controls_df: Dataframe containing Z-scores for all controls (matrix)
+#' @param nr_plots_perpage: Number of metabolites per page (integer)
+#' @param number_of_samples: Number of patients and controls (list of integers)
+#' @param number_of_metabolites: Number of metabolites for the top and lowest table (list of integers)
+#' @param iem_variables: Top number of IEMs and threshold for dIEM (list of integers)
+#' @param explanation_violin_plot: Text that explains the violin plots and the pipeline version (vector of strings)
#'
-#' @returns patient_no_iem: vector of patient IDs that have no IEMs
+#' @return patient_no_iem: Patient IDs for which no IEM data is available (vector of strings)
make_and_save_diem_plots <- function(
diem_probability_score,
patient_col_names,
@@ -910,6 +926,7 @@ make_and_save_diem_plots <- function(
number_of_metabolites,
iem_variables,
explanation_violin_plot) {
+ # create output folder
diem_plot_dir <- paste("./dIEM_plots", sep = "/")
dir.create(diem_plot_dir)
@@ -968,11 +985,11 @@ make_and_save_diem_plots <- function(
#' Get the IEM probabilities for a patient for all diseases
#'
-#' @param patient_top_iems_probs: dataframe containing the probability scores for diseases for a patient
-#' @param expected_biomarkers_df: dataframe with information for HMDB codes about IEMs
-#' @param patient_id: string containing the patien ID
+#' @param patient_top_iems_probs: Dataframe containing the probability scores for diseases for a patient (matrix)
+#' @param expected_biomarkers_df: Dataframe with information for HMDB codes about IEMs (matrix)
+#' @param patient_id: Patient ID (string)
#'
-#' @returns list_metabolites_iems: list of dataframes containing the HMDB codes and names for all diseases
+#' @return list_metabolites_iems: Dataframes containing the HMDB codes and names for all diseases (list of matrices)
get_probabilities_top_iems <- function(patient_top_iems_probs, expected_biomarkers_df, patient_id) {
# Get the metabolites for each IEM and their probability
list_metabolites_iems <- list()
@@ -998,8 +1015,8 @@ get_probabilities_top_iems <- function(patient_top_iems_probs, expected_biomarke
#' Save a list of patient IDs to a text file
#'
-#' @param threshold_iem: integer containing the IEM threshold
-#' @param patient_no_iem: vector containing patient IDs
+#' @param threshold_iem: IEM threshold (integer)
+#' @param patient_no_iem: Patient IDs (vector of strings)
save_patient_no_iem <- function(threshold_iem, patient_no_iem) {
patient_no_iem <- c(
paste0(
diff --git a/DIMS/preprocessing/average_peaks_functions.R b/DIMS/preprocessing/average_peaks_functions.R
index beed29bb..b7ad6f6a 100644
--- a/DIMS/preprocessing/average_peaks_functions.R
+++ b/DIMS/preprocessing/average_peaks_functions.R
@@ -1,13 +1,15 @@
+# function for averaging peak intensities for different technical replicates
+
+#' Average the intensity of peaks that occur in different technical replicates of a biological sample
+#'
+#' @param peaklist_allrepl_sorted: Dataframe with peaks sorted on median m/z (float)
+#' @param sample_name: String with sample name (string)
+#'
+#' @return averaged_peaks: matrix of averaged peaks (float)
average_peaks_per_sample <- function(peaklist_allrepl_sorted, sample_name) {
- #' Average the intensity of peaks that occur in different technical replicates of a biological sample
- #'
- #' @param peaklist_allrepl_sorted: Dataframe with peaks sorted on median m/z (float)
- #'
- #' @return averaged_peaks: matrix of averaged peaks (float)
-
# initialize
averaged_peaks <- peaklist_allrepl_sorted[0, ]
- # set ppm as fixed value, not the same ppm as in peak grouping
+ # set ppm as fixed value, not the same ppm as in peak group annotation
ppm_peak <- 2
while (nrow(peaklist_allrepl_sorted) > 1) {
diff --git a/DIMS/preprocessing/collect_filled_functions.R b/DIMS/preprocessing/collect_filled_functions.R
index abcebbee..57ad982c 100644
--- a/DIMS/preprocessing/collect_filled_functions.R
+++ b/DIMS/preprocessing/collect_filled_functions.R
@@ -1,13 +1,13 @@
-# CollectFilled functions
+# functions for collecting peak groups after fill missing
+#' Collapse identification info for peak groups with the same mass
+#'
+#' @param column_label: Name of column in peakgroup_list (string)
+#' @param peakgroup_list: Peak group list (matrix)
+#' @param index_dup: Index of duplicate peak group (integer)
+#'
+#' @return collapsed_items: Semicolon-separated list of info (string)
collapse_information <- function(column_label, peakgroup_list, index_dup) {
- #' Collapse identification info for peak groups with the same mass
- #'
- #' @param column_label: Name of column in peakgroup_list (string)
- #' @param peakgroup_list: Peak group list (matrix)
- #' @param index_dup: Index of duplicate peak group (integer)
- #'
- #' @return collapsed_items: Semicolon-separated list of info (string)
# get the item(s) that need to be collapsed
list_items <- as.vector(peakgroup_list[index_dup, column_label])
# remove NA
@@ -18,18 +18,18 @@ collapse_information <- function(column_label, peakgroup_list, index_dup) {
return(collapsed_items)
}
+#' Merge identification info for peak groups with the same mass
+#'
+#' @param peakgroup_list: Peak group list (matrix)
+#'
+#' @return peakgroup_list_dedup: de-duplicated peak group list (matrix)
merge_duplicate_rows <- function(peakgroup_list) {
- #' Merge identification info for peak groups with the same mass
- #'
- #' @param peakgroup_list: Peak group list (matrix)
- #'
- #' @return peakgroup_list_dedup: de-duplicated peak group list (matrix)
-
- options(digits = 16)
- collect <- NULL
- remove <- NULL
+ # initialize
+ dedup_peakgroups <- NULL
+ remove_peakgroup_indices <- NULL
# check for peak groups with identical mass
+ peakgroup_list[, "mzmed.pgrp"] <- as.character(peakgroup_list[, "mzmed.pgrp"])
index_dup <- which(duplicated(peakgroup_list[, "mzmed.pgrp"]))
while (length(index_dup) > 0) {
@@ -43,35 +43,38 @@ merge_duplicate_rows <- function(peakgroup_list) {
single_peakgroup[, "HMDB_code"] <- collapse_information("HMDB_code", peakgroup_list, peaklist_index)
single_peakgroup[, "all_hmdb_ids"] <- collapse_information("all_hmdb_ids", peakgroup_list, peaklist_index)
single_peakgroup[, "sec_hmdb_ids"] <- collapse_information("sec_hmdb_ids", peakgroup_list, peaklist_index)
- if (single_peakgroup[, "sec_hmdb_ids"] == ";") single_peakgroup[, "sec_hmdb_ids"] < NA
+ if (single_peakgroup[, "sec_hmdb_ids"] == ";") {
+ single_peakgroup[, "sec_hmdb_ids"] <- NA
+ }
# keep track of deduplicated entries
- collect <- rbind(collect, single_peakgroup)
- remove <- c(remove, peaklist_index)
+ dedup_peakgroups <- rbind(dedup_peakgroups, single_peakgroup)
+ remove_peakgroup_indices <- c(remove_peakgroup_indices, peaklist_index)
# remove current entry from index
index_dup <- index_dup[-which(peakgroup_list[index_dup, "mzmed.pgrp"] == peakgroup_list[index_dup[1], "mzmed.pgrp"])]
}
# remove duplicate entries
- if (!is.null(remove)) {
- peakgroup_list <- peakgroup_list[-remove, ]
+ if (!is.null(remove_peakgroup_indices)) {
+ peakgroup_list <- peakgroup_list[-remove_peakgroup_indices, ]
}
# append deduplicated entries
- peakgroup_list_dedup <- rbind(peakgroup_list, collect)
+ peakgroup_list_dedup <- rbind(peakgroup_list, dedup_peakgroups)
+ peakgroup_list[, "mzmed.pgrp"] <- as.numeric(peakgroup_list[, "mzmed.pgrp"])
return(peakgroup_list_dedup)
}
+#' Calculate Z-scores for peak groups based on average and standard deviation of controls
+#'
+#' @param peakgroup_list: Peak group list (matrix)
+#'
+#' @return peakgroup_list_zscores: peak group list with Z-scores (matrix)
calculate_zscores_peakgrouplist <- function(peakgroup_list) {
- #' Calculate Z-scores for peak groups based on average and standard deviation of controls
- #'
- #' @param peakgroup_list: Peak group list (matrix)
- #'
- #' @return peakgroup_list_dedup: de-duplicated peak group list (matrix)
-
+ # define case and control labels (fixed)
case_label <- "P"
control_label <- "C"
- # get index for new column names
+ # get index for new columns
startcol <- ncol(peakgroup_list) + 3
# calculate mean and standard deviation for Control group
ctrl_cols <- grep(control_label, colnames(peakgroup_list), fixed = TRUE)
@@ -84,49 +87,48 @@ calculate_zscores_peakgrouplist <- function(peakgroup_list) {
peakgroup_list$sd.ctrls <- apply(ctrl_ints, 1, function(x) sd(as.numeric(x), na.rm = TRUE))
# set new column names and calculate Z-scores
- colnames_zscores <- NULL
+ peakgroup_list_zscores <-peakgroup_list
+ colnames_zscores <- paste0(colnames(peakgroup_list)[int_cols], "_Zscore")
for (col_index in int_cols) {
- col_name <- colnames(peakgroup_list)[col_index]
- colnames_zscores <- c(colnames_zscores, paste0(col_name, "_Zscore"))
zscores_1col <- (as.numeric(as.vector(unlist(peakgroup_list[, col_index]))) -
peakgroup_list$avg.ctrls) / peakgroup_list$sd.ctrls
- peakgroup_list <- cbind(peakgroup_list, zscores_1col)
+ peakgroup_list_zscores <- cbind(peakgroup_list_zscores, zscores_1col)
}
# apply new column names to columns at end plus avg and sd columns
- colnames(peakgroup_list)[startcol:ncol(peakgroup_list)] <- colnames_zscores
+ colnames(peakgroup_list_zscores)[startcol:ncol(peakgroup_list_zscores)] <- colnames_zscores
- return(peakgroup_list)
+ return(peakgroup_list_zscores)
}
+#' Calculate ppm deviation between observed mass and expected theoretical mass
+#'
+#' @param peakgroup_list: Peak group list (matrix)
+#'
+#' @return peakgroup_list_ppm: peak group list with ppm column (matrix)
calculate_ppm_deviation <- function(peakgroup_list) {
- #' Calculate ppm deviation between observed mass and expected theoretical mass
- #'
- #' @param peakgroup_list: Peak group list (matrix)
- #'
- #' @return peakgroup_list_ppm: peak group list with ppm column (matrix)
-
# make sure values in columns mzmed.pgrp and theormz_HMDB are numeric
peakgroup_list$mzmed.pgrp <- as.numeric(peakgroup_list$mzmed.pgrp)
peakgroup_list$theormz_HMDB <- as.numeric(peakgroup_list$theormz_HMDB)
# calculate ppm deviation
- for (row_index in seq_len(nrow(peakgroup_list))) {
- observed_mz <- peakgroup_list$mzmed.pgrp[row_index]
- theor_mz <- peakgroup_list$theormz_HMDB[row_index]
- peakgroup_list$ppmdev[row_index] <- 10^6 * (observed_mz - theor_mz) / theor_mz
+ peakgroup_list_ppm <- peakgroup_list
+ for (row_index in seq_len(nrow(peakgroup_list_ppm))) {
+ observed_mz <- peakgroup_list_ppm$mzmed.pgrp[row_index]
+ theor_mz <- peakgroup_list_ppm$theormz_HMDB[row_index]
+ peakgroup_list_ppm$ppmdev[row_index] <- 10^6 * (observed_mz - theor_mz) / theor_mz
}
- return(peakgroup_list)
+ return(peakgroup_list_ppm)
}
+#' Put columns in peak group list in correct order
+#'
+#' @param peakgroup_list: Peak group list (matrix)
+#'
+#' @return peakgroup_ordered: peak group list with columns in correct order (matrix)
order_columns_peakgrouplist <- function(peakgroup_list) {
- #' Put columns in peak group list in correct order
- #'
- #' @param peakgroup_list: Peak group list (matrix)
- #'
- #' @return peakgroup_ordered: peak group list with columns in correct order (matrix)
-
+ # retain original column names
original_colnames <- colnames(peakgroup_list)
mass_columns <- c(grep("mzm", original_colnames), grep("nrsamples", original_colnames))
if (any(grepl("avg.int", original_colnames))) {
@@ -140,7 +142,6 @@ order_columns_peakgrouplist <- function(peakgroup_list) {
zscore_columns <- grep("_Zscore", original_colnames)
# create peak group list with columns in correct order
peakgroup_ordered <- peakgroup_list[ , c(mass_columns, descriptive_columns, intensity_columns, control_columns, zscore_columns)]
-
+
return(peakgroup_ordered)
}
-
diff --git a/DIMS/preprocessing/collect_sum_adducts_functions.R b/DIMS/preprocessing/collect_sum_adducts_functions.R
index a65e3f7d..45e30310 100644
--- a/DIMS/preprocessing/collect_sum_adducts_functions.R
+++ b/DIMS/preprocessing/collect_sum_adducts_functions.R
@@ -1,37 +1,38 @@
+# functions for combining sum adducts
+
+#' Combine all AdductSum parts in 1 dataframe
+#'
+#' @param scanmode: string with the scan mode, either "positive" or "negative"
+#'
+#' @returns: adductsums_all: dataframe with all adducts for one scan mode
combine_sum_adduct_parts <- function(scanmode) {
- #' Combine all AdductSum parts in 1 dataframe
- #'
- #' @param scanmode: string with the scanmodus, either positive or negative
- #'
- #' @returns: outlist_tot: dataframe with all adducts for one scanmodus
adductsum_part_files <- list.files("./", pattern = scanmode)
- outlist_tot <- NULL
+ adductsums_all <- NULL
for (i in seq_along(adductsum_part_files)) {
load(adductsum_part_files[i])
- outlist_tot <- rbind(outlist_tot, adductsum)
+ adductsums_all <- rbind(adductsums_all, adductsum)
}
- return(outlist_tot)
+ return(adductsums_all)
}
+#' Combine the scan modes; add intensities if both scan modes are present
+#'
+#' @param outlist_pos_adducts_hmdb: intensities for adducts in the positive scan mode (data frame)
+#' @param outlist_neg_adducts_hmdb: intensities for adducts in the negative scan mode (data frame)
+#'
+#' @returns: outlist: dataframe with intensities for all metabolites present in either or both scanmodes
combine_scanmodes_intensities <- function(outlist_pos_adducts_hmdb, outlist_neg_adducts_hmdb) {
- #' Combine the scanmodes and when present in both scanmodes add intensities
- #'
- #' @param outlist_pos_adducts_hmdb: dataframe with adducts for the positive scanmodus
- #' @param outlist_neg_adducts_hmdb: dataframe with adducts for the positive scanmodus
- #'
- #' @returns: outlist: dataframe with intensities for all metabolites present in either or both scanmodes
-
- # Only continue with patients (columns) that are in both pos and neg, so patients that are in both
+ # Only continue with samples (columns) that are in both pos and neg
samples_both_modes <- intersect(colnames(outlist_neg_adducts_hmdb), colnames(outlist_pos_adducts_hmdb))
outlist_neg_adducts_hmdb <- outlist_neg_adducts_hmdb[, samples_both_modes]
outlist_pos_adducts_hmdb <- outlist_pos_adducts_hmdb[, samples_both_modes]
- # Find indexes of neg hmdb code that are also found in pos and vice versa
+ # Find indices of metabolites in neg that are also found in pos and vice versa
index_neg <- which(rownames(outlist_neg_adducts_hmdb) %in% rownames(outlist_pos_adducts_hmdb))
index_pos <- which(rownames(outlist_pos_adducts_hmdb) %in% rownames(outlist_neg_adducts_hmdb))
- # Get intensities of metabs present in both modes from pos modus
+ # Get intensities of metabs present in both modes from pos mode
outlist_combi_pos_ints <- outlist_pos_adducts_hmdb[rownames(outlist_pos_adducts_hmdb)[index_pos], ] %>% select(-c(HMDB_name, HMDB_name_all, HMDB_ID_all, sec_HMDB_ID))
# Get intensities of metabs present in both modes from neg modus
@@ -40,12 +41,12 @@ combine_scanmodes_intensities <- function(outlist_pos_adducts_hmdb, outlist_neg_
# HMDB info for metabs present in both modes
outlist_combi_info <- outlist_pos_adducts_hmdb[rownames(outlist_pos_adducts_hmdb)[index_pos], ] %>% select(HMDB_name, HMDB_name_all, HMDB_ID_all, sec_HMDB_ID)
- # Combine positive and negative numbers and paste back HMDB column
+ # Sum positive and negative intensities and put back HMDB columns
outlist_combi_ints <- apply(outlist_combi_pos_ints, 2, as.numeric) + apply(outlist_combi_neg_ints, 2, as.numeric)
rownames(outlist_combi_ints) <- rownames(outlist_combi_pos_ints)
outlist_combi <- cbind(outlist_combi_ints, outlist_combi_info)
- # Get outlist with metabs in either pos or neg
+ # Get remaining metabolites which are not present in both scan modes
outlist_pos <- outlist_pos_adducts_hmdb[-index_pos, ]
outlist_neg <- outlist_neg_adducts_hmdb[-index_neg, ]
diff --git a/DIMS/preprocessing/evaluate_tics_functions.R b/DIMS/preprocessing/evaluate_tics_functions.R
index 8cb782e2..08b01d8d 100644
--- a/DIMS/preprocessing/evaluate_tics_functions.R
+++ b/DIMS/preprocessing/evaluate_tics_functions.R
@@ -1,12 +1,13 @@
# EvaluateTics functions
+
+#' Find technical replicates with a total intensity below a threshold
+#'
+#' @param repl_pattern: List of samples with corresponding technical replicates (strings)
+#' @param thresh2remove: Threshold value for acceptance or rejection of total intensity (integer)
+#'
+#' @return remove_tech_reps: Array of rejected technical replicates (strings)
find_bad_replicates <- function(repl_pattern, thresh2remove) {
- #' Find technical replicates with a total intensity below a threshold
- #'
- #' @param repl_pattern: List of samples with corresponding technical replicates (strings)
- #' @param thresh2remove: Threshold value for acceptance or rejection of total intensity (integer)
- #'
- #' @return remove_tech_reps: Array of rejected technical replicates (strings)
-
+ # initialize
remove_pos <- NULL
remove_neg <- NULL
cat("Pklist sum threshold to remove technical replicate:", thresh2remove, "\n")
@@ -49,21 +50,20 @@ find_bad_replicates <- function(repl_pattern, thresh2remove) {
col.names = FALSE,
sep = "\t"
)
-
+
# combine removed technical replicates from pos and neg
remove_tech_reps <- list(pos = remove_pos, neg = remove_neg)
return(remove_tech_reps)
}
+#' Remove technical replicates with insufficient quality from a biological sample
+#'
+#' @param bad_samples: Array of technical replicates of insufficient quality (strings)
+#' @param repl_pattern: List of samples with corresponding technical replicates (strings)
+#' @param nr_replicates: Number of technical replicates per biological sample (integer)
+#'
+#' @return repl_pattern_filtered: list of technical replicates of sufficient quality (strings)
remove_from_repl_pattern <- function(bad_samples, repl_pattern, nr_replicates) {
- #' Remove technical replicates with insufficient quality from a biological sample
- #'
- #' @param bad_samples: Array of technical replicates of insufficient quality (strings)
- #' @param repl_pattern: List of samples with corresponding technical replicates (strings)
- #' @param nr_replicates: Number of technical replicates per biological sample (integer)
- #'
- #' @return repl_pattern_filtered: list of technical replicates of sufficient quality (strings)
-
# collect list of samples to remove from replication pattern
remove_from_group <- NULL
for (sample_nr in 1:length(repl_pattern)){
@@ -89,14 +89,13 @@ remove_from_repl_pattern <- function(bad_samples, repl_pattern, nr_replicates) {
return(repl_pattern_filtered)
}
+#' Create an overview of technical replicates with sufficient quality from a biological sample
+#'
+#' @param repl_pattern_filtered: List of samples with corresponding technical replicates (strings)
+#' @param scanmode: Scan mode "positive" or "negative" (string)
+#'
+#' @return allsamples_techreps_scanmode: Matrix of technical replicates of sufficient quality (strings)
get_overview_tech_reps <- function(repl_pattern_filtered, scanmode) {
- #' Create an overview of technical replicates with sufficient quality from a biological sample
- #'
- #' @param repl_pattern_filtered: List of samples with corresponding technical replicates (strings)
- #' @param scanmode: Scan mode "positive" or "negative" (string)
- #'
- #' @return allsamples_techreps_scanmode: Matrix of technical replicates of sufficient quality (strings)
-
allsamples_techreps_scanmode <- matrix("", ncol = 3, nrow = length(repl_pattern_filtered))
for (sample_nr in 1:length(repl_pattern_filtered)) {
allsamples_techreps_scanmode[sample_nr, 1] <- names(repl_pattern_filtered)[sample_nr]
@@ -105,4 +104,3 @@ get_overview_tech_reps <- function(repl_pattern_filtered, scanmode) {
allsamples_techreps_scanmode[, 3] <- scanmode
return(allsamples_techreps_scanmode)
}
-
diff --git a/DIMS/preprocessing/fill_missing_functions.R b/DIMS/preprocessing/fill_missing_functions.R
index b55bc029..12f0ff1b 100644
--- a/DIMS/preprocessing/fill_missing_functions.R
+++ b/DIMS/preprocessing/fill_missing_functions.R
@@ -1,13 +1,14 @@
-fill_missing_intensities <- function(peakgroup_list, repl_pattern, thresh, disable_randomness = FALSE) {
- #' Replace intensities that are zero with random value
- #'
- #' @param peakgroup_list: Peak groups (matrix)
- #' @param repl_pattern: Replication pattern (list of strings)
- #' @param thresh: Value for threshold between noise and signal (integer)
- #' @param thresh: Variable which indicates whether randomness should be disabled (boolean)
- #'
- #' @return final_outlist: peak groups with filled-in intensities (matrix)
+# function for fill missing (zero) intensities with random noise
+#' Replace intensities that are zero with random value
+#'
+#' @param peakgroup_list: Peak groups with zero intensities (matrix)
+#' @param repl_pattern: Replication pattern (list of strings)
+#' @param thresh: Value for threshold between noise and signal (integer)
+#' @param disable_randomness: Variable which indicates whether randomness should be disabled (boolean)
+#'
+#' @return final_outlist: peak groups with filled-in intensities (matrix)
+fill_missing_intensities <- function(peakgroup_list, repl_pattern, thresh, disable_randomness = FALSE) {
# for unit test, turn off randomness
if (disable_randomness) {
set.seed(123)
@@ -30,8 +31,8 @@ fill_missing_intensities <- function(peakgroup_list, repl_pattern, thresh, disab
# Add column with average intensity; find intensity columns first
int_cols <- which(colnames(peakgroup_list) %in% names(repl_pattern))
- peakgroup_list <- cbind(peakgroup_list, "avg.int" = apply(peakgroup_list[, int_cols], 1, mean))
+ final_outlist <- cbind(peakgroup_list, "avg.int" = apply(peakgroup_list[, int_cols], 1, mean))
- return(peakgroup_list)
+ return(final_outlist)
}
}
diff --git a/DIMS/preprocessing/peak_finding_functions.R b/DIMS/preprocessing/peak_finding_functions.R
index 30d16996..84469b2d 100644
--- a/DIMS/preprocessing/peak_finding_functions.R
+++ b/DIMS/preprocessing/peak_finding_functions.R
@@ -1,15 +1,15 @@
# functions for peak finding
+# NB: add package R.utils to docker image, remove seqToIntervals function from this file
+# replace seqToIntervals to R.utils::seqToIntervals in function search_regions_of_interest
+#' Divide the full m/z range into regions of interest with indices
+#'
+#' @param ints_fullrange: Matrix with m/z values and intensities (float)
+#'
+#' @return regions_of_interest: matrix of m/z regions of interest (integer)
search_regions_of_interest <- function(ints_fullrange) {
- #' Divide the full m/z range into regions of interest with indices
- #'
- #' @param ints_fullrange: Matrix with m/z values and intensities (float)
- #'
- #' @return regions_of_interest: matrix of m/z regions of interest (integer)
-
# find indices where intensity is not equal to zero
nonzero_positions <- as.vector(which(ints_fullrange$int != 0))
-
# find regions of interest (look for consecutive numbers)
regions_of_interest_consec <- seqToIntervals(nonzero_positions)
# add length of regions of interest
@@ -75,24 +75,23 @@ search_regions_of_interest <- function(ints_fullrange) {
# sort on first index
if (nrow(regions_of_interest_final) > 1){
- regions_of_interest_sorted <- regions_of_interest_final %>% as.data.frame %>% dplyr::arrange(from)
+ regions_of_interest <- regions_of_interest_final %>% as.data.frame %>% dplyr::arrange(from)
} else {
- regions_of_interest_sorted <- regions_of_interest_final
+ regions_of_interest <- regions_of_interest_final
}
- return(regions_of_interest_sorted)
+ return(regions_of_interest)
}
+#' Fit Gaussian peak for each region of interest and integrate area under the curve
+#'
+#' @param ints_fullrange: Named list of intensities (float)
+#' @param regions_of_interest: Named list of intensities (float)
+#' @param resol: Value for resolution (integer)
+#' @param peak_thresh: Value for noise level threshold (integer) # NOT USED YET
+#'
+#' @return allpeaks_values: matrix of integrated peaks
integrate_peaks <- function(ints_fullrange, regions_of_interest, resol, peak_thresh) {
- #' Fit Gaussian peak for each region of interest and integrate area under the curve
- #'
- #' @param ints_fullrange: Named list of intensities (float)
- #' @param regions_of_interest: Named list of intensities (float)
- #' @param resol: Value for resolution (integer)
- #' @param peak_thresh: Value for noise level threshold (integer) # NOT USED YET
- #'
- #' @return allpeaks_values: matrix of integrated peaks
-
# initialize dataframe to store results for all peaks
allpeaks_values <- matrix(0, nrow = nrow(regions_of_interest), ncol = 5)
colnames(allpeaks_values) <- c("mzmed.pkt", "fq", "mzmin.pkt", "mzmax.pkt", "height.pkt")
@@ -191,14 +190,13 @@ seqToIntervals <- function(idx) {
return(res)
}
+#' Calculate fwhm (full width at half maximum intensity) for a peak
+#'
+#' @param query_mass: Value for mass (float)
+#' @param resol: Value for resolution (integer)
+#'
+#' @return fwhm: Value for full width at half maximum (float)
get_fwhm <- function(query_mass, resol) {
- #' Calculate fwhm (full width at half maximum intensity) for a peak
- #'
- #' @param query_mass: Value for mass (float)
- #' @param resol: Value for resolution (integer)
- #'
- #' @return fwhm: Value for full width at half maximum (float)
-
# set aberrant values of query_mass to default value of 200
if (is.na(query_mass)) {
query_mass <- 200
@@ -214,25 +212,16 @@ get_fwhm <- function(query_mass, resol) {
return(fwhm)
}
-
-# from https://rdrr.io/cran/rvmethod/src/R/gaussfit.R
-#' Gaussian Function from package rvmethod
+#' Gaussian function based on stats::dnorm
+#' Calculate area under the Gaussian curve
+#'
+#' @param mz_range: Vector of values at which to evaluate the Gaussian (vector)
+#' @param mu: center of the Gaussian peak (float)
+#' @param sigma: Spread of the Gaussian peak (float)
#'
-#' This function returns the unnormalized (height of 1.0) Gaussian curve with a
-#' given center and spread.
-#'
-#' @param x the vector of values at which to evaluate the Gaussian
-#' @param mu the center of the Gaussian
-#' @param sigma the spread of the Gaussian (must be greater than 0)
-#' @return vector of values of the Gaussian
-#' @examples x = seq(-4, 4, length.out = 100)
-#' y = gaussfunc(x, 0, 1)
-#' plot(x, y)
-#'
-#' @import stats
-#'
-#' @export
-gaussfunc <- function(x, mu, sigma) {
- return(exp(-((x - mu) ^ 2) / (2 * (sigma ^ 2))))
-}
+#' @return calc_area: integrated area under the Gaussian curve (float)
+gaussfunc <- function(mz_range, mu, sigma) {
+ calc_area <- stats::dnorm(mz_range, mu, sigma) / stats::dnorm(mu, mu, sigma)
+ return(calc_area)
+}
diff --git a/DIMS/preprocessing/peak_grouping_functions.R b/DIMS/preprocessing/peak_grouping_functions.R
index 93ff96f5..9251f6c7 100644
--- a/DIMS/preprocessing/peak_grouping_functions.R
+++ b/DIMS/preprocessing/peak_grouping_functions.R
@@ -1,13 +1,13 @@
-# functions for PeakGrouping
-find_peak_groups <- function(outlist_sorted, mz_tolerance, sample_names) {
- #' find peaks in all samples with query m/z values and form peak groups
- #'
- #' @param outlist_sorted: matrix of peaks (mzmed, intensity) in all samples
- #' @param mz_tolerance: Value for mass tolerance around query m/z (float)
- #' @param sample_names: vector of sample names (vector of strings)
- #'
- #' @return ints_sorted: matrix of peak groups
+# functions for generating peak groups and annotating them
+#' find peaks in all samples with query m/z values and form peak groups
+#'
+#' @param outlist_sorted: matrix of peaks (mzmed, intensity) in all samples
+#' @param mz_tolerance: Value for mass tolerance around query m/z (float)
+#' @param sample_names: vector of sample names (vector of strings)
+#'
+#' @return ints_sorted: matrix of peak groups
+find_peak_groups <- function(outlist_sorted, mz_tolerance, sample_names) {
# set up object for intensities for all samples
ints_allsamps <- matrix(0, nrow = nrow(outlist_sorted), ncol = 3 + (length(sample_names)))
colnames(ints_allsamps) <- c("mzmed.pgrp", "mzmin.pgrp", "mzmax.pgrp", sample_names)
@@ -63,16 +63,15 @@ find_peak_groups <- function(outlist_sorted, mz_tolerance, sample_names) {
return(ints_sorted)
}
+#' annotate peak groups; assign metabolites (adducts, isotopes) with suitable mass from HMDB
+#'
+#' @param ints_sorted: matrix of peak groups
+#' @param hmdb_add_iso: subset of HMDB (matrix)
+#' @param column_label: column name with appropriate m/z values for scan mode (string)
+#' @param mz_tolerance: Value for mass tolerance around query m/z (float)
+#'
+#' @return peakgrouplist_identified: matrix of peak groups with annotation
annotate_peak_groups <- function(ints_sorted, hmdb_add_iso, column_label, mz_tolerance) {
- #' annotate peak groups; assign metabolites (adducts, isotopes) with suitable mass from HMDB
- #'
- #' @param ints_sorted: matrix of peak groups
- #' @param hmdb_add_iso: subset of HMDB (matrix)
- #' @param column_label: column name with appropriate m/z values for scan mode (string)
- #' @param mz_tolerance: Value for mass tolerance around query m/z (float)
- #'
- #' @return peakgrouplist_identified: matrix of peak groups with annotation
-
# Initialize matrix for annotation
assigned_hmdb <- matrix("", nrow = nrow(ints_sorted), ncol = 7)
colnames(assigned_hmdb) <- c("assi_HMDB", "all_hmdb_names", "iso_HMDB", "HMDB_code",
@@ -80,7 +79,7 @@ annotate_peak_groups <- function(ints_sorted, hmdb_add_iso, column_label, mz_tol
# make sure isotope entries in the hmdb part have no HMDB IDs in the HMDB_ID_all column
hmdb_add_iso[grep("iso", rownames(hmdb_add_iso)), "HMDB_ID_all"] <- ""
-
+
# for each peak group, find all entries in HMDB part with mass within ppm range
for (row_number in 1:nrow(ints_sorted)) {
# initialize to make sure there's no information from the previous peak group
@@ -99,14 +98,12 @@ annotate_peak_groups <- function(ints_sorted, hmdb_add_iso, column_label, mz_tol
reference_mass <- ints_sorted[row_number, "mzmed.pgrp"]
# select indices for all HMDB entries with mass between +/- ppm tolerance
select_from_hmdb <- which(hmdb_add_iso[, column_label] > (reference_mass - mz_tolerance) &
- hmdb_add_iso[, column_label] < (reference_mass + mz_tolerance))
+ hmdb_add_iso[, column_label] < (reference_mass + mz_tolerance))
if (length(select_from_hmdb) > 0) {
# get dataframe of all entries which are selected
select_hmdb_df <- hmdb_add_iso[select_from_hmdb, ]
# separate into main metabolites, metabolites with adducts and isotopes
# main metabolites have no "_" in their name
- # if there are rownames with "_", choose only those without "_"
- # check if any rownames contain a hmdb code without "_"
if (any(!grepl("_", rownames(select_hmdb_df)))) {
grep_noiso_noadduct <- which(!grepl("_", rownames(select_hmdb_df)))
} else {
@@ -161,8 +158,8 @@ annotate_peak_groups <- function(ints_sorted, hmdb_add_iso, column_label, mz_tol
peakgrouplist_identified[, "theormz_HMDB"] <- as.numeric(peakgrouplist_identified[, "theormz_HMDB"])
peakgrouplist_identified[, "mzmed.pgrp"] <- as.numeric(peakgrouplist_identified[, "mzmed.pgrp"])
peakgrouplist_identified[, "ppmdev"] <- 1000000 * (peakgrouplist_identified[, "mzmed.pgrp"] -
- peakgrouplist_identified[, "theormz_HMDB"]) /
- peakgrouplist_identified[, "theormz_HMDB"]
+ peakgrouplist_identified[, "theormz_HMDB"]) /
+ peakgrouplist_identified[, "theormz_HMDB"]
return(peakgrouplist_identified)
}
diff --git a/DIMS/preprocessing/sum_intensities_adducts.R b/DIMS/preprocessing/sum_intensities_adducts.R
index 30c9cdc9..4217b36f 100644
--- a/DIMS/preprocessing/sum_intensities_adducts.R
+++ b/DIMS/preprocessing/sum_intensities_adducts.R
@@ -1,12 +1,14 @@
+# functions for summing intensities from different adducts of the same metabolite
+
+#' Sum intensities for different adducts of the same metabolite
+#'
+#' @param peakgroup_list: Peak group list (matrix)
+#' @param hmdb_part: Matrix of metabolites , part of the HMDB (matrix)
+#' @param adducts: Vector of adducts (vector of integers)
+#' @param z_score: Value indicating whether Z-scores have been calculated (integer)
+#'
+#' @return adductsum: peak group list with summed intensities (matrix)
sum_intensities_adducts <- function(peakgroup_list, hmdb_part, adducts, z_score) {
- #' Sum intensities for different adducts of the same metabolite
- #'
- #' @param peakgroup_list: Peak group list (matrix)
- #' @param hmdb_part: Matrix of metabolites , part of the HMDB (matrix)
- #' @param adducts: Vector of adducts (vector of integers)
- #' @param z_score: Value indicating whether Z-scores have been calculated (integer)
- #'
- #' @return adductsum: peak group list with summed intensities (matrix)
hmdb_part_info <- cbind(HMDB_id = rownames(hmdb_part), CompoundName = hmdb_part[, "CompoundName"])
# create overview of row indices for each metabolite_adduct combination in peaklist
@@ -15,7 +17,7 @@ sum_intensities_adducts <- function(peakgroup_list, hmdb_part, adducts, z_score)
# avoid rows with only "" in HMDB_code column
hmdb_in_peaklist[which(hmdb_in_peaklist == "")] <- ";"
hmdb_in_peaklist_rownr <- c()
-
+
# create dataframe with for each HMDB id a row number
hmdb_in_peaklist_rownr <- data.frame(
row_id = rep(seq_along(hmdb_in_peaklist), lengths(hmdb_in_peaklist)),
@@ -78,4 +80,3 @@ sum_intensities_adducts <- function(peakgroup_list, hmdb_part, adducts, z_score)
return(adductsum)
}
-
diff --git a/DIMS/tests/testthat/_snaps/generate_qc_output/test_barplot.png b/DIMS/tests/testthat/_snaps/generate_qc_output/test_barplot.png
index 9024fa35..b656b83a 100644
Binary files a/DIMS/tests/testthat/_snaps/generate_qc_output/test_barplot.png and b/DIMS/tests/testthat/_snaps/generate_qc_output/test_barplot.png differ
diff --git a/DIMS/tests/testthat/_snaps/generate_violin_plots.md b/DIMS/tests/testthat/_snaps/generate_violin_plots.md
index 598ec3fd..83d6efc1 100644
--- a/DIMS/tests/testthat/_snaps/generate_violin_plots.md
+++ b/DIMS/tests/testthat/_snaps/generate_violin_plots.md
@@ -3,19 +3,19 @@
Code
content_pdf_violinplots
Output
- [1] "Top deviating metabolites for patient: P2025M1\n Metabolite Z.score\n Increased\n metab1 2.45\n Decreased\n metab11 −1.51\n"
- [2] " Results for patient P2025M1\n test acyl carnitines\n metab1 Z=0.31\nMetabolites\n metab3 Z=2.34\n −5 0 5 10 15 20\n Z−scores\n"
- [3] " Results for patient P2025M1\n test crea gua\n metab4 Z=0.84\nMetabolites\n metab11 Z=−0.46\n −5 0 5 10 15 20\n Z−scores\n"
- [4] " Unit test Generate Violin Plots\nUnit test Generate Violin Plots\n"
+ [1] "Top deviating metabolites for patient: P2025M1\n\n\n\n\n Metabolite Z.score\n Increased\n\n metab1 2.45\n\n Decreased\n\n metab11 −1.51\n"
+ [2] " Results for patient P2025M1\n test acyl carnitines\n\n\n\n\n metab1 Z=0.31\nMetabolites\n\n\n\n\n metab3 Z=2.34\n\n\n\n\n −5 0 5 10 15 20\n Z−scores\n"
+ [3] " Results for patient P2025M1\n test crea gua\n\n\n\n\n metab4 Z=0.84\nMetabolites\n\n\n\n\n metab11 Z=−0.46\n\n\n\n\n −5 0 5 10 15 20\n Z−scores\n"
+ [4] " Unit test Generate Violin Plots\nUnit test Generate Violin Plots\n"
# save_prob_scores_to_excel: Saving the probability score dataframe as an Excel file
- Disease P2025M1 P2025M2 P2025M3 P2025M4
- 1 Disease A 10.900 -10.90000000000000 49.90000000000000 -49.9
- 2 Disease B 0.953 0.00000000000000 2.29000000000000 0.0
- 3 Disease C 12.100 0.00000000000000 0.00000000000000 12.1
- 4 Disease D 0.000 -12.50000000000000 0.00000000000000 18.2
- 5 Disease E 44.300 0.00000000000000 0.00000000000000 28.1
- 6 Disease F 0.000 -77.40000000000001 -77.40000000000001 0.0
- 7 Disease G -38.700 38.70000000000000 38.70000000000000 -38.7
+ Disease P2025M1 P2025M2 P2025M3 P2025M4
+ 1 Disease A 10.900 -10.9 49.90 -49.9
+ 2 Disease B 0.953 0.0 2.29 0.0
+ 3 Disease C 12.100 0.0 0.00 12.1
+ 4 Disease D 0.000 -12.5 0.00 18.2
+ 5 Disease E 44.300 0.0 0.00 28.1
+ 6 Disease F 0.000 -77.4 -77.40 0.0
+ 7 Disease G -38.700 38.7 38.70 -38.7
diff --git a/DIMS/tests/testthat/_snaps/generate_violin_plots/violin-plot-p2025m1.svg b/DIMS/tests/testthat/_snaps/generate_violin_plots/violin-plot-p2025m1.svg
index 89edec18..4fdc0eee 100644
--- a/DIMS/tests/testthat/_snaps/generate_violin_plots/violin-plot-p2025m1.svg
+++ b/DIMS/tests/testthat/_snaps/generate_violin_plots/violin-plot-p2025m1.svg
@@ -40,8 +40,8 @@
-
-
+
+
Z=2.34
diff --git a/DIMS/tests/testthat/test_collect_filled.R b/DIMS/tests/testthat/test_collect_filled.R
index 0b711963..72953f04 100644
--- a/DIMS/tests/testthat/test_collect_filled.R
+++ b/DIMS/tests/testthat/test_collect_filled.R
@@ -20,10 +20,9 @@ testthat::test_that("Duplicate rows in a peak group list are correctly merged",
test_peakgroup_list_dup <- test_peakgroup_list[c(1, 2, 2, 3), ]
# after merging duplicate rows, the test peak group list should have 3 rows
- expect_equal(nrow(merge_duplicate_rows(test_peakgroup_list_dup)), 3, TRUE, tolerance = 0.001)
+ expect_equal(nrow(merge_duplicate_rows(test_peakgroup_list_dup)), 3, tolerance = 0.001)
expect_equal(merge_duplicate_rows(test_peakgroup_list_dup)[3, "all_hmdb_ids"],
- paste(test_peakgroup_list_dup[2, "all_hmdb_ids"], test_peakgroup_list_dup[3, "all_hmdb_ids"], sep = ";"),
- TRUE)
+ paste(test_peakgroup_list_dup[2, "all_hmdb_ids"], test_peakgroup_list_dup[3, "all_hmdb_ids"], sep = ";"))
})
testthat::test_that("Z-scores are correctly calculated in CollectFilled", {
@@ -35,11 +34,11 @@ testthat::test_that("Z-scores are correctly calculated in CollectFilled", {
test_peakgroup_list_noz <- test_peakgroup_list_noz[ , -grep("_Zscore", colnames(test_peakgroup_list_noz))]
# after calculate_zscores_peakgrouplist, there should be 4 columns with _Zscore in the name
- expect_equal(length(grep("_Zscore", colnames(calculate_zscores_peakgrouplist(test_peakgroup_list_noz)))), 4, TRUE, tolerance = 0.001)
+ expect_equal(length(grep("_Zscore", colnames(calculate_zscores_peakgrouplist(test_peakgroup_list_noz)))), 4, tolerance = 0.001)
# after calculate_zscores_peakgrouplist, the 4 columns with _Zscore in the name should be filled non-zero
- expect_equal(calculate_zscores_peakgrouplist(test_peakgroup_list_noz)$C101.1_Zscore[1], -0.7071, TRUE, tolerance = 0.00001)
- expect_equal(calculate_zscores_peakgrouplist(test_peakgroup_list_noz)$P2.1_Zscore[4], 12.0208, TRUE, tolerance = 0.00001)
+ expect_equal(calculate_zscores_peakgrouplist(test_peakgroup_list_noz)$C101.1_Zscore[1], -0.7071, tolerance = 0.00001)
+ expect_equal(calculate_zscores_peakgrouplist(test_peakgroup_list_noz)$P2.1_Zscore[4], 12.0208, tolerance = 0.00001)
})
diff --git a/DIMS/tests/testthat/test_generate_excel.R b/DIMS/tests/testthat/test_generate_excel.R
index e4a7ff38..8e56e8b8 100644
--- a/DIMS/tests/testthat/test_generate_excel.R
+++ b/DIMS/tests/testthat/test_generate_excel.R
@@ -34,40 +34,40 @@ testthat::test_that("get_intensities_cols: Get indices of columns and dataframe
testthat::test_that("calculate_zscores: Calculating Z-scores using different methods for excluding controls", {
test_outlist <- read.delim(test_path("fixtures", "test_outlist.txt"))
+ test_outlist <- as.data.frame(test_outlist)
control_intensities <- read.delim(test_path("fixtures", "test_control_intensities.txt"))
control_col_idx <- c(2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13)
intensity_col_ids <- c(2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15)
- startcol <- ncol(test_outlist) + 4
perc <- 5
outlier_threshold <- 2
- expect_type(calculate_zscores(test_outlist, "_Zscore", control_intensities, NULL, intensity_col_ids, startcol), "list")
- expect_identical(colnames(calculate_zscores(test_outlist, "_Zscore", control_intensities, NULL, intensity_col_ids, startcol)),
+ expect_type(calculate_zscores(test_outlist, "_Zscore", control_col_idx, NULL, intensity_col_ids), "list")
+ expect_identical(colnames(calculate_zscores(test_outlist, "_Zscore", control_col_idx, NULL, intensity_col_ids)),
c("plots", "C101.1", "C102.1", "C103.1", "C104.1", "C105.1", "C106.1", "C107.1", "C108.1", "C109.1", "C110.1",
"C111.1", "C112.1", "P2.1", "P3.1", "HMDB_name", "HMDB_name_all", "HMDB_ID_all", "sec_HMDB_ID",
"HMDB_key", "sec_HMDB_ID_rlvnc", "name", "relevance", "descr", "origin", "fluids", "tissue", "disease",
"pathway", "HMDB_code", "avg_ctrls", "sd_ctrls", "nr_ctrls", "C101.1_Zscore", "C102.1_Zscore", "C103.1_Zscore",
"C104.1_Zscore", "C105.1_Zscore", "C106.1_Zscore", "C107.1_Zscore", "C108.1_Zscore", "C109.1_Zscore",
"C110.1_Zscore", "C111.1_Zscore", "C112.1_Zscore", "P2.1_Zscore", "P3.1_Zscore"))
- expect_equal(round(calculate_zscores(test_outlist, "_Zscore", control_intensities, NULL, intensity_col_ids, startcol)$avg_ctrls, 3),
+ expect_equal(round(calculate_zscores(test_outlist, "_Zscore", control_col_idx, NULL, intensity_col_ids)$avg_ctrls, 3),
c(16129.167, 1150.0, 1231.250, 4015.833), tolerance = 0.001)
- expect_equal(calculate_zscores(test_outlist, "_Zscore", control_intensities, NULL, intensity_col_ids, startcol)$P2.1_Zscore,
+ expect_equal(calculate_zscores(test_outlist, "_Zscore", control_col_idx, NULL, intensity_col_ids)$P2.1_Zscore,
c(-0.2544103, 32.4586955, 13.6066674, 0.4037668), tolerance = 0.001)
- expect_type(calculate_zscores(test_outlist, "_RobustZscore", control_col_idx, perc, intensity_col_ids, startcol), "list")
- expect_identical(colnames(calculate_zscores(test_outlist, "_RobustZscore", control_col_idx, perc, intensity_col_ids, startcol))[34:47],
+ expect_type(calculate_zscores(test_outlist, "_RobustZscore", control_col_idx, perc, intensity_col_ids), "list")
+ expect_identical(colnames(calculate_zscores(test_outlist, "_RobustZscore", control_col_idx, perc, intensity_col_ids))[34:47],
c("C101.1_RobustZscore", "C102.1_RobustZscore",
"C103.1_RobustZscore", "C104.1_RobustZscore", "C105.1_RobustZscore", "C106.1_RobustZscore",
"C107.1_RobustZscore", "C108.1_RobustZscore", "C109.1_RobustZscore", "C110.1_RobustZscore",
"C111.1_RobustZscore", "C112.1_RobustZscore", "P2.1_RobustZscore", "P3.1_RobustZscore"))
- expect_equal(calculate_zscores(test_outlist, "_RobustZscore", control_col_idx, perc, intensity_col_ids, startcol)$avg_ctrls,
+ expect_equal(calculate_zscores(test_outlist, "_RobustZscore", control_col_idx, perc, intensity_col_ids)$avg_ctrls,
c(1255.0, 1110.0, 1227.5, 2811.5), tolerance = 0.001)
- expect_equal(calculate_zscores(test_outlist, "_RobustZscore", control_col_idx, perc, intensity_col_ids, startcol)$P2.1_RobustZscore,
+ expect_equal(calculate_zscores(test_outlist, "_RobustZscore", control_col_idx, perc, intensity_col_ids)$P2.1_RobustZscore,
c(9.1511750, 46.9804468, 16.8039663, 0.8565111), tolerance = 0.001)
- expect_type(calculate_zscores(test_outlist, "_OutlierRemovedZscore", control_col_idx, outlier_threshold, intensity_col_ids, startcol), "list")
- expect_identical(colnames(calculate_zscores(test_outlist, "_OutlierRemovedZscore", control_col_idx, outlier_threshold, intensity_col_ids, startcol))[34:47],
+ expect_type(calculate_zscores(test_outlist, "_OutlierRemovedZscore", control_col_idx, outlier_threshold, intensity_col_ids), "list")
+ expect_identical(colnames(calculate_zscores(test_outlist, "_OutlierRemovedZscore", control_col_idx, outlier_threshold, intensity_col_ids))[34:47],
c("C101.1_OutlierRemovedZscore",
"C102.1_OutlierRemovedZscore", "C103.1_OutlierRemovedZscore", "C104.1_OutlierRemovedZscore",
"C105.1_OutlierRemovedZscore", "C106.1_OutlierRemovedZscore", "C107.1_OutlierRemovedZscore",
@@ -75,11 +75,11 @@ testthat::test_that("calculate_zscores: Calculating Z-scores using different met
"C111.1_OutlierRemovedZscore", "C112.1_OutlierRemovedZscore", "P2.1_OutlierRemovedZscore",
"P3.1_OutlierRemovedZscore")
)
- expect_equal(calculate_zscores(test_outlist, "_OutlierRemovedZscore", control_col_idx, outlier_threshold, intensity_col_ids, startcol)$avg_ctrls,
+ expect_equal(calculate_zscores(test_outlist, "_OutlierRemovedZscore", control_col_idx, outlier_threshold, intensity_col_ids)$avg_ctrls,
c(1231.818, 1077.273, 1231.250, 2649.091), tolerance = 0.001)
- expect_equal(calculate_zscores(test_outlist, "_OutlierRemovedZscore", control_col_idx, outlier_threshold, intensity_col_ids, startcol)$nr_ctrls,
+ expect_equal(calculate_zscores(test_outlist, "_OutlierRemovedZscore", control_col_idx, outlier_threshold, intensity_col_ids)$nr_ctrls,
c(11, 11, 12, 11))
- expect_equal(calculate_zscores(test_outlist, "_OutlierRemovedZscore", control_col_idx, outlier_threshold, intensity_col_ids, startcol)$P2.1_OutlierRemovedZscore,
+ expect_equal(calculate_zscores(test_outlist, "_OutlierRemovedZscore", control_col_idx, outlier_threshold, intensity_col_ids)$P2.1_OutlierRemovedZscore,
c(8.9955723, 44.9136860, 13.6066674, 0.9345077), tolerance = 0.001)
})
diff --git a/DIMS/tests/testthat/test_generate_qc_output.R b/DIMS/tests/testthat/test_generate_qc_output.R
index 07be0178..d4c23dde 100644
--- a/DIMS/tests/testthat/test_generate_qc_output.R
+++ b/DIMS/tests/testthat/test_generate_qc_output.R
@@ -60,7 +60,7 @@ testthat::test_that("Save internal standard plots", {
expect_silent(
save_internal_standard_plot(
- test_plot_data, "barplot", "Test barplot", temp_dir,
+ test_plot_data, "barplot", "Test barplot",
file_name_barplot, 6, 4, test_hline_data
)
)
@@ -71,7 +71,7 @@ testthat::test_that("Save internal standard plots", {
out_file_lineplot <- file.path(temp_dir, "plots", paste0(file_name_lineplot, ".png"))
expect_silent(
save_internal_standard_plot(
- test_plot_data, "lineplot", "Test lineplot", temp_dir,
+ test_plot_data, "lineplot", "Test lineplot",
file_name_lineplot, 6, 4
)
)
@@ -80,7 +80,7 @@ testthat::test_that("Save internal standard plots", {
test_plot_data <- test_plot_data[-c(1, 2, 3, 4, 5, 6), ]
expect_identical(save_internal_standard_plot(
- test_plot_data, "barplot", "Test barplot", temp_dir,
+ test_plot_data, "barplot", "Test barplot",
file_name_barplot, 6, 4, test_hline_data
), NULL)
@@ -88,7 +88,7 @@ testthat::test_that("Save internal standard plots", {
expect_silent(
save_internal_standard_plot(
- test_plot_data, "barplot", "Test barplot", temp_dir,
+ test_plot_data, "barplot", "Test barplot",
file_name_barplot_select, 6, 4, test_hline_data
)
)
diff --git a/DIMS/tests/testthat/test_generate_violin_plots.R b/DIMS/tests/testthat/test_generate_violin_plots.R
index d5fe3a8e..45b18356 100644
--- a/DIMS/tests/testthat/test_generate_violin_plots.R
+++ b/DIMS/tests/testthat/test_generate_violin_plots.R
@@ -676,7 +676,7 @@ testthat::test_that("add_zscores_ratios_to_df: Add Zscores for multiple ratios t
test_metabolites_ratios_df,
test_all_sample_ids
)$C101.1,
- c(1000, 1200, 1300, 1400, 1500, 1600, -0.2630344, -0.1069152, 11.8533096)
+ c(1000, 1200, 1300, 1400, 1500, 1600, 0.83333, 0.9285714, 3700.00000), tolerance = 0.0001
)
expect_equal(
add_zscores_ratios_to_df(
@@ -684,7 +684,7 @@ testthat::test_that("add_zscores_ratios_to_df: Add Zscores for multiple ratios t
test_metabolites_ratios_df,
test_all_sample_ids
)$C101.1_Zscore,
- c(0.45, 1.67, -1.86, 0.58, 2.47, -0.56, -0.5899371, 0.4858991, -0.4552026),
+ c(0.45, 1.67, -1.86, 0.58, 2.47, -0.56, -0.4574, 0.5552, -0.4486),
tolerance = 0.0001
)
})
@@ -711,11 +711,11 @@ testthat::test_that("calculate_zscore_ratios: Calculate Zscores for ratios", {
)
expect_equal(
calculate_zscore_ratios(test_metabolites_ratios_df, test_outlist_df, test_all_sample_ids)$C101.1,
- c(-0.2630344, -0.1069152, 11.8533096)
+ c(0.83333, 0.9285714, 3700.00000), tolerance = 0.0001
)
expect_equal(
calculate_zscore_ratios(test_metabolites_ratios_df, test_outlist_df, test_all_sample_ids)$C101.1_Zscore,
- c(-0.5899371, 0.4858991, -0.4552026),
+ c(-0.4574, 0.5552, -0.4486),
tolerance = 0.0001
)
})