From 70b2a363e5abf349bb7b020adca04b9082b5547b Mon Sep 17 00:00:00 2001 From: Mia Pras-Raves Date: Tue, 14 Jul 2026 14:16:27 +0200 Subject: [PATCH 01/10] added code for ChemBL (DrugDB) metabolites in DIMS/HMDBparts_main.R --- DIMS/HMDBparts_main.R | 17 +++++++++++++++++ 1 file changed, 17 insertions(+) diff --git a/DIMS/HMDBparts_main.R b/DIMS/HMDBparts_main.R index 1a377eb..72de7ce 100644 --- a/DIMS/HMDBparts_main.R +++ b/DIMS/HMDBparts_main.R @@ -39,4 +39,21 @@ for (scanmode in scanmodes) { save(outlist_part, file = paste0(scanmode, "_hmdb_main.", part_index, ".RData")) start_index = start_index + 1000 } + + # add hmdb parts for drugs in ChEMBL database + CHEMBL_add <- HMDB_mzrange[grep("CHEMBL", rownames(HMDB_mzrange), fixed = TRUE), ] + # remove adducts + CHEMBL_main <- CHEMBL_add[-grep("_", rownames(CHEMBL_add), fixed = TRUE), ] + # sort on m/z value + CHEMBL_main <- CHEMBL_main[order(CHEMBL_main[, column_label]), ] + + # generate hmdb parts of 1000 lines each for ChEMBL entries + nr_parts_chembl <- ceiling(nrow(CHEMBL_main) / 1000) + start_index <- 1 + for (part_index in 1:nr_parts_chembl) { + end_index <- min((start_index + 999), nrow(CHEMBL_main)) + outlist_part <- CHEMBL_main[start_index:end_index, ] + save(outlist_part, file = paste0(scanmode, "_hmdb_main.", part_index + nr_parts, ".RData")) + start_index = start_index + 1000 + } } From 667f0bdae3e638a5d8405b53bdae4cc0eb8e7bf8 Mon Sep 17 00:00:00 2001 From: Mia Pras-Raves Date: Tue, 14 Jul 2026 14:18:13 +0200 Subject: [PATCH 02/10] added code for ChemBL (DrugDB) metabolites in DIMS GenerateExcel step --- DIMS/GenerateExcel.R | 96 +++++++++++++++++++++----- DIMS/export/generate_excel_functions.R | 6 +- 2 files changed, 82 insertions(+), 20 deletions(-) diff --git a/DIMS/GenerateExcel.R b/DIMS/GenerateExcel.R index f33c27d..dd35c8a 100644 --- a/DIMS/GenerateExcel.R +++ b/DIMS/GenerateExcel.R @@ -27,8 +27,6 @@ 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 @@ -39,6 +37,77 @@ load(hmdb_rlvc_file) # load outlist object load("AdductSums_combined.RData") +# Get columns with control intensities +control_intensity_cols <- get_intensities_cols(outlist, control_label) +control_col_idx <- control_intensity_cols$col_idx +control_intensities <- control_intensity_cols$df_intensities + +# Get columns with patient intensities +patient_intensity_cols <- get_intensities_cols(outlist, case_label) +patient_col_idx <- patient_intensity_cols$col_idx +patient_columns <- colnames(patient_intensity_cols$df_intensities) + +intensity_col_ids <- c(control_col_idx, patient_col_idx) + +# Filter out drugs and other metabolites with CHEMBL annotation +outlist_drugdb <- outlist[grep("CHEMBL", outlist$HMDB_ID_all), ] +# Add CHEMBL_code column with all the CHEMBL ID and sort on it +outlist_drugdb <- cbind(outlist_drugdb, "CHEMBL_code" = rownames(outlist_drugdb)) +outlist_drugdb <- outlist_drugdb[order(outlist_drugdb[, "CHEMBL_code"]), ] + +# Create excel for drugs +sheetname <- "Drugs" +wb_intensities_drugs <- openxlsx::createWorkbook("AdductSums") +openxlsx::addWorksheet(wb_intensities_drugs, sheetname) + +dir.create("plots_drugs", showWarnings = FALSE) +# add a column for plots +outlist_drugdb <- cbind(plots = NA, outlist_drugdb) +# if there are any intensities of 0 left, set them to NA for stats +outlist_drugdb[, intensity_col_ids][outlist_drugdb[, intensity_col_ids] == 0] <- NA +outlist_drugdb$HMDB_key <- rownames(outlist_drugdb) +# get intensity columns +intensities_df <- outlist_drugdb[, c(ncol(outlist_drugdb), intensity_col_ids)] + +for (row_index in seq_len(nrow(intensities_df))) { + # get CHEMBL ID + chembl_id <- intensities_df %>% + slice(row_index) %>% + pull(HMDB_key) + + # Transform dataframe to long format + intensities_df_long <- intensities_df_to_long_format(intensities_df, row_index) + + # set plot width to 40 times the number of samples + plot_width <- length(unique(intensities_df_long$Samples)) * 40 + col_width <- plot_width * 2 + + start_row_index <- row_index + 1 + save_plot_to_excel_workbook( + wb_intensities_drugs, + sheetname, + intensities_df_long, + "plots_drugs/plot_", + chembl_id, + plot_width, + col_width, + start_row_index + ) +} +wb_intensities <- set_row_height_col_width_wb( + wb_intensities_drugs, + sheetname, + nrow(outlist), + ncol(outlist), + col_width, + plots_present = TRUE +) + +# write Excel file +openxlsx::writeData(wb_intensities_drugs, sheet = 1, outlist, startCol = 1) +openxlsx::saveWorkbook(wb_intensities_drugs, paste0("Drugs_", project, ".xlsx"), overwrite = TRUE) +rm(wb_intensities_drugs) +unlink("plots", recursive = TRUE) # Filter for biological relevance peaks_in_list <- which(rownames(outlist) %in% rlvnc$HMDB_key) @@ -65,27 +134,18 @@ 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 # add a column for plots outlist <- cbind(plots = NA, outlist) + control_col_idx <- control_col_idx + 1 + intensity_col_ids <- intensity_col_ids + 1 + # three columns will be added for mean, stdev and number of controls; Z-scores start at ncol + 4 startcol <- ncol(outlist) + 4 - # Get columns with control intensities - control_intensity_cols <- get_intensities_cols(outlist, control_label) - control_col_idx <- control_intensity_cols$col_idx - control_intensities <- control_intensity_cols$df_intensities - - # Get columns with patient intensities - patient_intensity_cols <- get_intensities_cols(outlist, case_label) - patient_col_idx <- patient_intensity_cols$col_idx - patient_columns <- colnames(patient_intensity_cols$df_intensities) - - intensity_col_ids <- c(control_col_idx, patient_col_idx) - # if there are any intensities of 0 left, set them to NA for stats outlist[, intensity_col_ids][outlist[, intensity_col_ids] == 0] <- NA @@ -99,7 +159,7 @@ if (z_score == 1) { ) # 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, startcol) # output metabolites filtered on relevance save_to_rdata_and_txt(outlist, "AdductSums_filtered_Zscores") @@ -213,7 +273,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 +297,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/export/generate_excel_functions.R b/DIMS/export/generate_excel_functions.R index 9cf5936..b8384d1 100644 --- a/DIMS/export/generate_excel_functions.R +++ b/DIMS/export/generate_excel_functions.R @@ -9,6 +9,8 @@ get_intensities_cols <- function(outlist, label) { #' 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) + # remove Z-score columns + col_idx <- col_idx[!grepl("_Zscore", colnames(outlist)[col_idx], fixed = TRUE)] df_intensities <- as.data.frame(outlist[, col_idx]) colnames(df_intensities) <- colnames(outlist)[col_idx] return(list(col_idx = col_idx, df_intensities = df_intensities)) @@ -33,8 +35,8 @@ calculate_zscores <- function(outlist, zscore_type, control_cols, stat_filter, i 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)) + outlist$avg_ctrls <- apply(outlist[, control_cols], 1, function(x) mean(as.numeric(x), na.rm = TRUE)) + outlist$sd_ctrls <- apply(outlist[, 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))) { From 6337c25d24936d67d1305ad3e2f57a81267ec410 Mon Sep 17 00:00:00 2001 From: Mia Pras-Raves Date: Fri, 24 Jul 2026 10:17:37 +0200 Subject: [PATCH 03/10] added output for DrugDB in DIMS/GenerateExcel.nf --- DIMS/GenerateExcel.nf | 1 + 1 file changed, 1 insertion(+) diff --git a/DIMS/GenerateExcel.nf b/DIMS/GenerateExcel.nf index 552a8ee..1d6a52f 100644 --- a/DIMS/GenerateExcel.nf +++ b/DIMS/GenerateExcel.nf @@ -15,6 +15,7 @@ process GenerateExcel { tuple path("AdductSums_filtered_Zscores.txt"), path("AdductSums_filtered_robustZ.txt"), path("AdductSums_filtered_outliersremovedZ.txt"), optional: true path("${analysis_id}.xlsx"), emit: project_excel path("Helix_${analysis_id}.xlsx"), optional: true + path("Drugs_${analysis_id}.xlsx"), optional: true script: """ From 46bb4a23e8f9916201c396338841bd12f8c1b749 Mon Sep 17 00:00:00 2001 From: Mia Pras-Raves Date: Fri, 24 Jul 2026 10:21:53 +0200 Subject: [PATCH 04/10] refactored DIMS/GenerateExcel.R for extra DrugDB output Excel file --- DIMS/GenerateExcel.R | 327 ++++++++++++------------------------------- 1 file changed, 89 insertions(+), 238 deletions(-) diff --git a/DIMS/GenerateExcel.R b/DIMS/GenerateExcel.R index dd35c8a..985a437 100644 --- a/DIMS/GenerateExcel.R +++ b/DIMS/GenerateExcel.R @@ -26,7 +26,6 @@ plot <- TRUE export <- TRUE control_label <- "C" case_label <- "P" - # percentage of outliers to remove from calculation of robust scaler perc <- 5 # Z-score for removing outliers with grubbs test @@ -35,268 +34,120 @@ outlier_threshold <- 2 # load HMDB rlvnc table load(hmdb_rlvc_file) -# load outlist object +# load adduct sums object (outlist) load("AdductSums_combined.RData") -# Get columns with control intensities -control_intensity_cols <- get_intensities_cols(outlist, control_label) -control_col_idx <- control_intensity_cols$col_idx -control_intensities <- control_intensity_cols$df_intensities - -# Get columns with patient intensities -patient_intensity_cols <- get_intensities_cols(outlist, case_label) -patient_col_idx <- patient_intensity_cols$col_idx -patient_columns <- colnames(patient_intensity_cols$df_intensities) - -intensity_col_ids <- c(control_col_idx, patient_col_idx) +## DrugDB output # Filter out drugs and other metabolites with CHEMBL annotation outlist_drugdb <- outlist[grep("CHEMBL", outlist$HMDB_ID_all), ] -# Add CHEMBL_code column with all the CHEMBL ID and sort on it -outlist_drugdb <- cbind(outlist_drugdb, "CHEMBL_code" = rownames(outlist_drugdb)) -outlist_drugdb <- outlist_drugdb[order(outlist_drugdb[, "CHEMBL_code"]), ] - -# Create excel for drugs -sheetname <- "Drugs" -wb_intensities_drugs <- openxlsx::createWorkbook("AdductSums") -openxlsx::addWorksheet(wb_intensities_drugs, sheetname) - -dir.create("plots_drugs", showWarnings = FALSE) -# add a column for plots -outlist_drugdb <- cbind(plots = NA, outlist_drugdb) -# if there are any intensities of 0 left, set them to NA for stats -outlist_drugdb[, intensity_col_ids][outlist_drugdb[, intensity_col_ids] == 0] <- NA -outlist_drugdb$HMDB_key <- rownames(outlist_drugdb) -# get intensity columns -intensities_df <- outlist_drugdb[, c(ncol(outlist_drugdb), intensity_col_ids)] - -for (row_index in seq_len(nrow(intensities_df))) { - # get CHEMBL ID - chembl_id <- intensities_df %>% - slice(row_index) %>% - pull(HMDB_key) - - # Transform dataframe to long format - intensities_df_long <- intensities_df_to_long_format(intensities_df, row_index) - - # set plot width to 40 times the number of samples - plot_width <- length(unique(intensities_df_long$Samples)) * 40 - col_width <- plot_width * 2 - - start_row_index <- row_index + 1 - save_plot_to_excel_workbook( - wb_intensities_drugs, - sheetname, - intensities_df_long, - "plots_drugs/plot_", - chembl_id, - plot_width, - col_width, - start_row_index - ) +if (nrow(outlist_drugdb) > 0) { + # Add HMDB_code column with all the CHEMBL ID and an empty description column + outlist_drugdb <- cbind(outlist_drugdb, HMDB_code = rownames(outlist_drugdb), descr = NA) + # sort on CHEMBL ID + outlist_drugdb <- outlist_drugdb[order(outlist_drugdb[, "HMDB_code"]), ] + + if (z_score == 1) { + # calculate Z-scores with outliers removed + outlist_drugdb_zscores <- calculate_zscores( + outlist_drugdb, "_OutlierRemovedZscore", outlier_threshold, + control_label, case_label + ) + # get indices for intensity columns based on control_label and case_label + intensity_col_ids <- get_intensity_col_index(outlist_drugdb_zscores, control_label, case_label) + # Create Excel for drugs + create_excel_output(outlist_drugdb_zscores, intensity_col_ids, "Drugs_", z_score, project) + } else if (z_score == 0) { + # get indices for intensity columns; all columns except those containing HMDB in the name + intensity_col_ids <- 1:ncol(outlist_drugdb) + intensity_col_ids <- intensity_col_ids[-grep("HMDB", colnames(outlist_drugdb))] + # Create Excel for drugs + create_excel_output(outlist_drugdb, intensity_col_ids, "Drugs_", z_score, project) + } } -wb_intensities <- set_row_height_col_width_wb( - wb_intensities_drugs, - sheetname, - nrow(outlist), - ncol(outlist), - col_width, - plots_present = TRUE -) - -# write Excel file -openxlsx::writeData(wb_intensities_drugs, sheet = 1, outlist, startCol = 1) -openxlsx::saveWorkbook(wb_intensities_drugs, paste0("Drugs_", project, ".xlsx"), overwrite = TRUE) -rm(wb_intensities_drugs) -unlink("plots", recursive = TRUE) +## Filtered metabolite output # Filter for biological relevance peaks_in_list <- which(rownames(outlist) %in% rlvnc$HMDB_key) outlist_subset <- outlist[peaks_in_list, ] +outlist_subset$HMDB_code <- rownames(outlist_subset) outlist_subset$HMDB_key <- rownames(outlist_subset) -outlist <- outlist_subset %>% +outlist_filtered <- outlist_subset %>% left_join(rlvnc %>% rename(sec_HMDB_ID_rlvnc = sec_HMDB_ID), by = "HMDB_key") -rownames(outlist) <- outlist$HMDB_key - +rownames(outlist_filtered) <- outlist_filtered$HMDB_key # filter out all irrelevant HMDBs -outlist <- outlist %>% +outlist_filtered <- outlist_filtered %>% tibble::rownames_to_column("rowname") %>% filter(grepl("relevant|Onbekend|Internal", relevance)) %>% tibble::column_to_rownames("rowname") +# sort on HMDB_key +outlist_filtered <- peakgroup_list[order(outlist_filtered[, "HMDB_key"]), ] -# Add HMDB_code column with all the HMDB ID and sort on it -outlist <- cbind(outlist, "HMDB_code" = rownames(outlist)) -outlist <- outlist[order(outlist[, "HMDB_code"]), ] - -# Create excel -sheetname <- "AllPeakGroups" -wb_intensities_zscores <- openxlsx::createWorkbook("SinglePatient") -openxlsx::addWorksheet(wb_intensities_zscores, sheetname) - -# Add Z-scores and create plots if (z_score == 1) { - 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 - # add a column for plots - outlist <- cbind(plots = NA, outlist) - control_col_idx <- control_col_idx + 1 - intensity_col_ids <- intensity_col_ids + 1 - - # three columns will be added for mean, stdev and number of controls; Z-scores start at ncol + 4 - startcol <- ncol(outlist) + 4 - - # if there are any intensities of 0 left, set them to NA for stats - 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) - - # 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 + # calculate Z-scores with outliers removed + outlist_filtered_zscores <- calculate_zscores( + outlist_filtered, "_OutlierRemovedZscore", outlier_threshold, + control_label, case_label ) - - # calculate Z-scores - outlist <- calculate_zscores(outlist, "_Zscore", control_col_idx, NULL, intensity_col_ids, startcol) - - # output metabolites filtered on relevance - save_to_rdata_and_txt(outlist, "AdductSums_filtered_Zscores") + colnames(outlist_filtered_zscores) <- gsub("_OutlierRemovedZscore", "_Zscore", colnames(outlist_filtered_zscores)) + # get indices for intensity columns + intensity_col_ids <- get_intensity_col_index(outlist_filtered_zscores, control_label, case_label) + # Create Excel for biologically relevant metabolites + create_excel_output(outlist_filtered_zscores, intensity_col_ids, "", z_score, project) + # save outlist for GenerateQC step + save(outlist_filtered_zscores, file = "outlist.RData") + # output filtered metabolites after removal of outliers + save_to_rdata_and_txt(outlist_filtered_zscores, "AdductSums_filtered_outliersremovedZ") + # calculate robust Z-scores + outlist_robust_zscore <- calculate_zscores(peakgroup_list, "_RobustZscore", control_col_idx, perc, intensity_col_ids, startcol) # output filtered metabolites with robust scaled Zscores save_to_rdata_and_txt(outlist_robust_zscore, "AdductSums_filtered_robustZ") - # output filtered metabolites after removal of outliers - save_to_rdata_and_txt(outlist_nooutliers, "AdductSums_filtered_outliersremovedZ") - - # use outlier-removed outlist for generating Excel file - outlist <- outlist_nooutliers - colnames(outlist) <- gsub("_OutlierRemovedZscore", "_Zscore", colnames(outlist)) - - # save outlist for GenerateQC step - save(outlist, file = "outlist.RData") + # calculate Z-scores without outlier removal + outlist <- calculate_zscores(peakgroup_list, "_Zscore", control_col_idx, NULL, intensity_col_ids, startcol) + # output metabolites filtered on relevance + save_to_rdata_and_txt(outlist, "AdductSums_filtered_Zscores") +} else if (z_score == 0) { + create_excel_output(outlist_filtered, intensity_col_ids, "", z_score, project) +} - # get Helix IDs for extra Excel file - metabolite_files <- list.files( - path = paste(path_metabolite_groups, "Diagnostics", sep = "/"), - pattern = "*.txt", full.names = FALSE, recursive = FALSE +## Helix output +# get Helix IDs for extra Excel file +metabolite_files <- list.files( + path = paste(path_metabolite_groups, "Diagnostics", sep = "/"), + pattern = "*.txt", full.names = FALSE, recursive = FALSE +) +metab_df_helix <- NULL +for (file_index in seq_along(metabolite_files)) { + infile <- metabolite_files[file_index] + metab_list <- read.table(paste(path_metabolite_groups, "Diagnostics", infile, sep = "/"), + sep = "\t", header = TRUE, quote = "" ) - metab_df_helix <- NULL - for (file_index in seq_along(metabolite_files)) { - infile <- metabolite_files[file_index] - metab_list <- read.table(paste(path_metabolite_groups, "Diagnostics", infile, sep = "/"), - sep = "\t", header = TRUE, quote = "" - ) - metab_df_helix <- rbind(metab_df_helix, metab_list) - } - # get Helix metabolites and unique HMDB IDs and remove ratio HMDBs containing A or L - metab_df_helix <- metab_df_helix %>% - filter(Helix == "ja") %>% - select(c(HMDB_code, HMDB_name)) %>% - rename(H_Name = HMDB_name) - metab_list_helix <- unique(metab_df_helix$HMDB_code) - metab_list_helix <- grep("[AL]", metab_list_helix, value = TRUE, invert = TRUE) + metab_df_helix <- rbind(metab_df_helix, metab_list) +} +# get Helix metabolites and unique HMDB IDs and remove ratio HMDBs containing A or L +metab_df_helix <- metab_df_helix %>% + filter(Helix == "ja") %>% + select(c(HMDB_code, HMDB_name)) %>% + rename(H_Name = HMDB_name) +metab_list_helix <- unique(metab_df_helix$HMDB_code) +metab_list_helix <- grep("[AL]", metab_list_helix, value = TRUE, invert = TRUE) - outlist_helix <- outlist %>% +if (z_score == 1) { + # get intensities for Helix metabolites from dataset + outlist_helix <- outlist_filtered_zscores %>% filter(HMDB_key %in% metab_list_helix) %>% left_join(., metab_df_helix, by = join_by(HMDB_code == HMDB_code)) %>% select( - -c(HMDB_key, sec_HMDB_ID_rlvnc, name, relevance, descr, origin, fluids, tissue, disease, pathway), - -all_of(control_col_idx), -all_of(patient_col_idx) - ) %>% - relocate(c(HMDB_code, H_Name, avg_ctrls, sd_ctrls), .after = plots) %>% - relocate(c(HMDB_name, HMDB_name_all, HMDB_ID_all, sec_HMDB_ID), .after = last_col()) %>% - rename(Name = H_Name) - - # Get intensity columns for controls and patients - intensities_df <- outlist %>% select(HMDB_key, matches("^C|^P[0-9]"), -ends_with("_Zscore")) - - for (row_index in seq_len(nrow(intensities_df))) { - # get HMDB ID - hmdb_id <- intensities_df %>% - slice(row_index) %>% - pull(HMDB_key) - - # Transform dataframe to long format - intensities_df_long <- intensities_df_to_long_format(intensities_df, row_index) - - # set plot width to 40 times the number of samples - plot_width <- length(unique(intensities_df_long$Samples)) * 40 - col_width <- plot_width * 2 - - if (hmdb_id %in% metab_list_helix) { - # Make separate plot for Helix Excel containing all samples - - start_row_index <- row_index + 1 - save_plot_to_excel_workbook( - wb_helix_zscores, - sheetname, - intensities_df_long, - "plots/plot_helix_", - hmdb_id, - plot_width, - col_width, - row_helix - ) - row_helix <- row_helix + 1 - } - - # Remove postive controls and SST mix samples, (e.g. P1001, P1002, P1003, P1005) - intensities_df_long <- intensities_df_long %>% filter(!grepl("^P[0-9]{4}$", Samples)) - - start_row_index <- row_index + 1 - save_plot_to_excel_workbook( - wb_intensities_zscores, - sheetname, - intensities_df_long, - "plots/plot_", - hmdb_id, - plot_width, - col_width, - start_row_index - ) - } - wb_intensities <- set_row_height_col_width_wb( - wb_intensities_zscores, - sheetname, - nrow(outlist), - ncol(outlist), - col_width, - plots_present = TRUE - ) - - wb_helix_intensities <- set_row_height_col_width_wb( - wb_helix_zscores, - sheetname, - nrow(outlist_helix), - ncol(outlist_helix), - col_width, - plots_present = TRUE - ) - openxlsx::writeData(wb_helix_intensities, sheet = 1, outlist_helix, startCol = 1) - openxlsx::saveWorkbook(wb_helix_intensities, paste0("/Helix_", project, ".xlsx"), overwrite = TRUE) - rm(wb_helix_intensities) - - # reorder outlist for Excel file - outlist <- outlist %>% - relocate(c(HMDB_code, HMDB_name_all, descr, avg_ctrls, sd_ctrls), .after = plots) %>% - relocate(all_of(grep("_Zscore", colnames(outlist))), .after = sd_ctrls) %>% - relocate(all_of(c(colnames(control_intensities), patient_columns)), .after = last_col()) -} else { - save(outlist, file = "outlist.RData") - wb_intensities <- set_row_height_col_width_wb( - wb_intensities_zscores, - sheetname, - nrow(outlist), - ncol(outlist), - plot_width = NULL, - plots_present = FALSE - ) - outlist <- outlist %>% - relocate(c(HMDB_name, HMDB_name_all, HMDB_code, HMDB_ID_all)) + -c(HMDB_key, sec_HMDB_ID_rlvnc, name, relevance, descr, origin, fluids, tissue, disease, pathway, monositopic_mass, molecular_formula) #, + ) +} else if (z_score == 0) { + # get intensities for Helix metabolites from dataset + outlist_helix <- outlist_filtered %>% + filter(HMDB_key %in% metab_list_helix) %>% + left_join(., metab_df_helix, by = join_by(HMDB_code == HMDB_code)) %>% + select( + -c(HMDB_key, sec_HMDB_ID_rlvnc, name, relevance, descr, origin, fluids, tissue, disease, pathway, monositopic_mass, molecular_formula) #, + ) } +# Create Excel for Helix +create_excel_output(outlist_helix, intensity_col_ids, "Helix_", z_score, project) -# write Excel file -openxlsx::writeData(wb_intensities_zscores, sheet = 1, outlist, startCol = 1) -openxlsx::saveWorkbook(wb_intensities_zscores, paste0(project, ".xlsx"), overwrite = TRUE) -rm(wb_intensities_zscores) -unlink("plots", recursive = TRUE) From 655b1266842e54f4eb8a2e11a63f7fb8275cb47c Mon Sep 17 00:00:00 2001 From: Mia Pras-Raves Date: Fri, 24 Jul 2026 10:26:38 +0200 Subject: [PATCH 05/10] modified function calculate_zscores, added function create_excel_output in DIMS/export/generate_excel_functions.R --- DIMS/export/generate_excel_functions.R | 218 ++++++++++++++++++++----- 1 file changed, 177 insertions(+), 41 deletions(-) diff --git a/DIMS/export/generate_excel_functions.R b/DIMS/export/generate_excel_functions.R index b8384d1..fc10863 100644 --- a/DIMS/export/generate_excel_functions.R +++ b/DIMS/export/generate_excel_functions.R @@ -16,62 +16,67 @@ get_intensities_cols <- function(outlist, label) { 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 mean and sd - outlist$avg_ctrls <- 0 - outlist$sd_ctrls <- 0 - outlist$nr_ctrls <- length(control_cols) - +#' 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: Method for excluding controls (string) +#' @param stat_filter: Either percentage or outlier threshold used for excluding controls (integer) +#' +#' @returns: peakgroup_list_zscores: same dataframe as the input with added Z-score columns (matrix) +calculate_zscores <- function(peakgroup_list, zscore_type, stat_filter, control_label = "C", case_label = "P") { + # Initialize + peakgroup_list$avg_ctrls <- 0 + peakgroup_list$sd_ctrls <- 0 + peakgroup_list$nr_ctrls <- length(control_col_idx) + + # Get columns with intensities + control_col_idx <- grep(control_label, colnames(peakgroup_list), fixed = TRUE) + control_column_names <- colnames(peakgroup_list)[control_col_idx] + patient_col_idx <- grep(case_label, colnames(peakgroup_list), fixed = TRUE) + patient_column_names <- colnames(peakgroup_list)[patient_col_idx] + intensity_col_ids <- c(control_col_idx, patient_col_idx) + + # calculate mean and standard deviation of controls if (zscore_type == "_Zscore") { - # Calculate mean and sd with all controls - outlist$avg_ctrls <- apply(outlist[, control_cols], 1, function(x) mean(as.numeric(x), na.rm = TRUE)) - outlist$sd_ctrls <- apply(outlist[, control_cols], 1, function(x) sd(as.numeric(x), na.rm = TRUE)) + # using all controls + peakgroup_list$avg_ctrls <- apply(peakgroup_list[, control_col_idx], 1, function(x) mean(as.numeric(x), na.rm = TRUE)) + peakgroup_list$sd_ctrls <- apply(peakgroup_list[, control_col_idx], 1, function(x) sd(as.numeric(x), na.rm = TRUE)) } else { - if (length(control_cols) > 3) { - for (metabolite_index in seq_len(nrow(outlist))) { + if (length(control_col_idx) >= 3) { + 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], - control_cols, stat_filter + # remove outlier controls by using robust scaler + peakgroup_list$avg_ctrls[metabolite_index] <- mean(robust_scaler( + peakgroup_list[metabolite_index, control_col_idx], + control_col_idx, stat_filter )) - outlist$sd_ctrls[metabolite_index] <- sd(robust_scaler( - outlist[metabolite_index, control_cols], - control_cols, stat_filter + peakgroup_list$sd_ctrls[metabolite_index] <- sd(robust_scaler( + peakgroup_list[metabolite_index, control_col_idx], + control_col_idx, stat_filter )) } else { - # Calculate mean, sd and number of remaining controls, remove outlier controls by using grubbs test + # 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_col_idx]), 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 + peakgroup_list_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) - - return(outlist) + colnames(peakgroup_list_zscores) <- paste0(colnames(peakgroup_list)[intensity_col_ids], zscore_type) + colnames(peakgroup_list_zscores) <- gsub("_OutlierRemovedZscore", "_Zscore", colnames(peakgroup_list_zscores)) + peakgroup_list_zscores <- cbind(peakgroup_list, peakgroup_list_zscores) + + return(peakgroup_list_zscores) } robust_scaler <- function(control_intensities, control_col_ids, perc = 5) { @@ -242,3 +247,134 @@ save_plot_to_excel_workbook <- function(excel_workbook, return(excel_workbook) } + +#' Get the indices for intensity columns in a matrix +#' +#' @param peakgroup_list: Dataframe with intensities and possibly Z-scores for all samples (matrix) +#' @param control_label: part of name of all control samples (string) +#' @param case_label: part of name of all patient samples (string) +#' +#' @returns intensity_indices: indices of columns with intensities (vector of integers) +get_intensity_col_index <- function(peakgroup_list, control_label = "C", case_label = "P") { + # remove Zscore columns first + if (any(grep("_Zscore", colnames(peakgroup_list)))) { + peakgroup_list <- peakgroup_list[, -grep("_Zscore", colnames(peakgroup_list))] + } + # get indices for controls + control_col_ids <- grep(control_label, colnames(peakgroup_list), fixed = TRUE) + # get indices for patients + patient_col_ids <- grep(case_label, colnames(peakgroup_list), fixed = TRUE) + # combine + intensity_indices <- c(control_col_ids, patient_col_ids) + + return(intensity_indices) +} + +#' Create an Excel workbook +#' +#' @param peakgroup_list: Dataframe with intensities for all samples (matrix) +#' @param intensity_col_ids: Indices for intensity columns (vector of integers) +#' @param prefix: Prefix to insert before file name (string) +#' @param z_score: Boolean number indicating whether Z-scores should be calculated (integer) +#' @param project: Name of dataset (string) +create_excel_output <- function(peakgroup_list, intensity_col_ids, prefix = "", z_score = 0, project) { + # set up + if (prefix == "Drugs_") { + plotdir <- "plots_drugs" + sheetname <- "AllDrugs" + wb_intensities <- openxlsx::createWorkbook("PeakGroups") + } else if (prefix == "Helix_") { + plotdir <- "plots_helix" + sheetname <- "HelixSelection" + wb_intensities <- openxlsx::createWorkbook("HelixOverview") + } else { + plotdir <- "plots" + sheetname <- "BiologicallyRelevant" + wb_intensities <- openxlsx::createWorkbook("AdductSums") + } + + openxlsx::addWorksheet(wb_intensities, sheetname) + + # Add Z-scores and create plots + if (z_score == 1) { + dir.create(plotdir, showWarnings = FALSE) + row_helix <- 2 # start on row 2 because of header + # get intensity columns + intensities_df <- as.data.frame(peakgroup_list[ , intensity_col_ids]) + sample_names <- colnames(intensities_df) + intensities_df$HMDB_key <- rownames(peakgroup_list) + # add a column for plots + peakgroup_list <- cbind(plots = NA, peakgroup_list) + + for (row_index in seq_len(nrow(intensities_df))) { + # get HMDB ID + hmdb_id <- intensities_df %>% + slice(row_index) %>% + pull(HMDB_key) + + # Transform dataframe to long format + intensities_df_long <- intensities_df_to_long_format(intensities_df, row_index) + + # set plot width to 40 times the number of samples + plot_width <- length(unique(intensities_df_long$Samples)) * 40 + col_width <- plot_width * 2 + + # Remove postive controls and SST mix samples, (e.g. P1001, P1002, P1003, P1005) + intensities_df_long <- intensities_df_long %>% filter(!grepl("^P[0-9]{4}$", Samples)) + + start_row_index <- row_index + 1 + save_plot_to_excel_workbook( + wb_intensities, + sheetname, + intensities_df_long, + paste0(plotdir, "/plot_"), + hmdb_id, + plot_width, + col_width, + start_row_index + ) + } + wb_intensities <- set_row_height_col_width_wb( + wb_intensities, + sheetname, + nrow(peakgroup_list), + ncol(peakgroup_list), + col_width, + plots_present = TRUE + ) + + # reorder outlist for Excel file + peakgroup_list <- peakgroup_list %>% + relocate(c(HMDB_code, HMDB_name_all, avg_ctrls, sd_ctrls), .after = plots) %>% + relocate(all_of(grep("_Zscore", colnames(peakgroup_list))), .after = sd_ctrls) %>% + relocate(all_of(sample_names), .after = last_col()) + if (prefix == "Helix_") { + colnames(peakgroup_list) <- gsub("H_Name", "Name", colnames(peakgroup_list)) + # remove intensity columns + peakgroup_list <- peakgroup_list %>% select(-all_of(sample_names)) + } + } else if (z_score == 0) { + save(peakgroup_list, file = "outlist.RData") + if (!any(grepl("HMDB_code", colnames(peakgroup_list)))) { + peakgroup_list$HMDB_code <- rownames(peakgroup_list) + } + wb_intensities <- set_row_height_col_width_wb( + wb_intensities, + sheetname, + nrow(peakgroup_list), + ncol(peakgroup_list), + plot_width = NULL, + plots_present = FALSE + ) + peakgroup_list <- peakgroup_list %>% + relocate(c(HMDB_name, HMDB_name_all, HMDB_code, HMDB_ID_all, sec_HMDB_ID)) + colnames(peakgroup_list) <- gsub("H_Name", "Name", colnames(peakgroup_list)) + } + + # write Excel file + openxlsx::writeData(wb_intensities, sheet = 1, peakgroup_list, startCol = 1) + openxlsx::saveWorkbook(wb_intensities, paste0(prefix, project, ".xlsx"), overwrite = TRUE) + rm(wb_intensities) + unlink(plotdir, recursive = TRUE) +} + From 77e6c22282bdf07bb53c2233146c6172e712b1a0 Mon Sep 17 00:00:00 2001 From: Mia Pras-Raves Date: Mon, 27 Jul 2026 11:22:52 +0200 Subject: [PATCH 06/10] refactored DIMS/GenerateExcel for extra output DrugDB --- DIMS/GenerateExcel.R | 7 ++++--- DIMS/GenerateExcel.nf | 1 + 2 files changed, 5 insertions(+), 3 deletions(-) diff --git a/DIMS/GenerateExcel.R b/DIMS/GenerateExcel.R index 985a437..f745fd0 100644 --- a/DIMS/GenerateExcel.R +++ b/DIMS/GenerateExcel.R @@ -80,7 +80,7 @@ outlist_filtered <- outlist_filtered %>% filter(grepl("relevant|Onbekend|Internal", relevance)) %>% tibble::column_to_rownames("rowname") # sort on HMDB_key -outlist_filtered <- peakgroup_list[order(outlist_filtered[, "HMDB_key"]), ] +outlist_filtered <- outlist_filtered[order(outlist_filtered[, "HMDB_key"]), ] if (z_score == 1) { # calculate Z-scores with outliers removed @@ -91,6 +91,7 @@ if (z_score == 1) { colnames(outlist_filtered_zscores) <- gsub("_OutlierRemovedZscore", "_Zscore", colnames(outlist_filtered_zscores)) # get indices for intensity columns intensity_col_ids <- get_intensity_col_index(outlist_filtered_zscores, control_label, case_label) + control_col_idx <- get_intensity_col_index(outlist_filtered_zscores, control_label, "none") # Create Excel for biologically relevant metabolites create_excel_output(outlist_filtered_zscores, intensity_col_ids, "", z_score, project) # save outlist for GenerateQC step @@ -98,11 +99,11 @@ if (z_score == 1) { # output filtered metabolites after removal of outliers save_to_rdata_and_txt(outlist_filtered_zscores, "AdductSums_filtered_outliersremovedZ") # calculate robust Z-scores - outlist_robust_zscore <- calculate_zscores(peakgroup_list, "_RobustZscore", control_col_idx, perc, intensity_col_ids, startcol) + outlist_robust_zscore <- calculate_zscores(outlist_filtered, "_RobustZscore", perc, control_label, case_label) # output filtered metabolites with robust scaled Zscores save_to_rdata_and_txt(outlist_robust_zscore, "AdductSums_filtered_robustZ") # calculate Z-scores without outlier removal - outlist <- calculate_zscores(peakgroup_list, "_Zscore", control_col_idx, NULL, intensity_col_ids, startcol) + outlist <- calculate_zscores(outlist_filtered, "_Zscore", NULL, control_label, case_label) # output metabolites filtered on relevance save_to_rdata_and_txt(outlist, "AdductSums_filtered_Zscores") } else if (z_score == 0) { diff --git a/DIMS/GenerateExcel.nf b/DIMS/GenerateExcel.nf index 1d6a52f..af27a7e 100644 --- a/DIMS/GenerateExcel.nf +++ b/DIMS/GenerateExcel.nf @@ -16,6 +16,7 @@ process GenerateExcel { path("${analysis_id}.xlsx"), emit: project_excel path("Helix_${analysis_id}.xlsx"), optional: true path("Drugs_${analysis_id}.xlsx"), optional: true + path("Drugs_in_dataset.RData"), emit: outlist_drugs, optional: true script: """ From aea03a9338af659044b0ff69cae629a314978e54 Mon Sep 17 00:00:00 2001 From: Mia Pras-Raves Date: Mon, 27 Jul 2026 11:23:29 +0200 Subject: [PATCH 07/10] extra functions for refactored DIMS/GenerateExcel for extra output DrugDB --- DIMS/export/generate_excel_functions.R | 6 ++---- 1 file changed, 2 insertions(+), 4 deletions(-) diff --git a/DIMS/export/generate_excel_functions.R b/DIMS/export/generate_excel_functions.R index fc10863..4facef9 100644 --- a/DIMS/export/generate_excel_functions.R +++ b/DIMS/export/generate_excel_functions.R @@ -27,14 +27,12 @@ calculate_zscores <- function(peakgroup_list, zscore_type, stat_filter, control_ # Initialize peakgroup_list$avg_ctrls <- 0 peakgroup_list$sd_ctrls <- 0 - peakgroup_list$nr_ctrls <- length(control_col_idx) - # Get columns with intensities + # Get columns indices with intensities control_col_idx <- grep(control_label, colnames(peakgroup_list), fixed = TRUE) - control_column_names <- colnames(peakgroup_list)[control_col_idx] patient_col_idx <- grep(case_label, colnames(peakgroup_list), fixed = TRUE) - patient_column_names <- colnames(peakgroup_list)[patient_col_idx] intensity_col_ids <- c(control_col_idx, patient_col_idx) + peakgroup_list$nr_ctrls <- length(control_col_idx) # calculate mean and standard deviation of controls if (zscore_type == "_Zscore") { From c7fa2313842fbbe53799cf7dce036fb8058d56e8 Mon Sep 17 00:00:00 2001 From: Mia Pras-Raves Date: Mon, 27 Jul 2026 11:24:16 +0200 Subject: [PATCH 08/10] modifications to DIMS/GenerateViolinPlots for DrugDB output --- DIMS/GenerateViolinPlots.R | 13 ++++++++++++- DIMS/GenerateViolinPlots.nf | 1 + 2 files changed, 13 insertions(+), 1 deletion(-) diff --git a/DIMS/GenerateViolinPlots.R b/DIMS/GenerateViolinPlots.R index 9d1f3e7..980edf0 100644 --- a/DIMS/GenerateViolinPlots.R +++ b/DIMS/GenerateViolinPlots.R @@ -17,12 +17,17 @@ path_metabolite_groups <- cmd_args[3] file_ratios_metabolites <- cmd_args[4] file_expected_biomarkers_iem <- cmd_args[5] file_explanation <- cmd_args[6] +file_drugs_in_dataset <- cmd_args[7] # 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) + +# load dataframe with intensities and Z-scores for drug metabolites +intensities_zscore_drugs_df <- get(load("Drugs_in_dataset.RData")) + # 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) @@ -50,6 +55,7 @@ number_of_metabolites <- list( control_ids <- get_colnames_by_prefix(intensities_zscore_df, "C") patient_ids <- get_colnames_by_prefix(intensities_zscore_df, "P") +patient_ids_drugs <- get_colnames_by_prefix(intensities_zscore_drugs_df, "P") all_sample_ids <- c(control_ids, patient_ids) number_of_samples <- list( controls = length(control_ids), @@ -68,11 +74,16 @@ zscore_patients_df <- intensities_zscore_ratios_df %>% 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")) +zscore_pat_drugs_df <- intensities_zscore_drugs_df %>% + select(HMDB_code, HMDB_name, any_of(paste0(patient_ids_drugs, "_Zscore"))) %>% + rename_with(~ str_remove(.x, "_Zscore"), .cols = contains("_Zscore")) + #### Make violin plots ##### make_and_save_violin_plot_pdfs( zscore_patients_df, zscore_controls_df, + zscore_pat_drugs_df, path_metabolite_groups, nr_plots_perpage, number_of_samples, diff --git a/DIMS/GenerateViolinPlots.nf b/DIMS/GenerateViolinPlots.nf index ec65a2e..24f16fe 100755 --- a/DIMS/GenerateViolinPlots.nf +++ b/DIMS/GenerateViolinPlots.nf @@ -7,6 +7,7 @@ process GenerateViolinPlots { input: path(outlist_zscores) val(analysis_id) + path(outlist_drugs) output: path('Diagnostics/*.pdf'), emit: diag_plot_files, optional: true From 9beceb74e9671aad7de006e5f78b97fff18e6d8e Mon Sep 17 00:00:00 2001 From: Mia Pras-Raves Date: Mon, 27 Jul 2026 11:25:12 +0200 Subject: [PATCH 09/10] modifications to DIMS/export/generate_violin_plots_functions.R for DrugDB output --- DIMS/export/generate_violin_plots_functions.R | 66 ++++++++++++++++--- 1 file changed, 58 insertions(+), 8 deletions(-) diff --git a/DIMS/export/generate_violin_plots_functions.R b/DIMS/export/generate_violin_plots_functions.R index 468bbde..1f53ee0 100644 --- a/DIMS/export/generate_violin_plots_functions.R +++ b/DIMS/export/generate_violin_plots_functions.R @@ -27,6 +27,11 @@ prepare_intensities_zscore_df <- function(intensities_zscore_df) { #' @returns sample_colnames: a vector of column names all containing the prefix. get_colnames_by_prefix <- function(dataframe, prefix) { sample_colnames <- grep(paste0("^", prefix), colnames(dataframe), value = TRUE) + # remove Z-score columns from intensity_col_names + if (any(grepl("_Zscore", sample_colnames))) { + sample_colnames <- sample_colnames[-grep("_Zscore", sample_colnames)] + } + return(sample_colnames) } @@ -119,6 +124,7 @@ calculate_zscore_ratios <- function(metabolites_ratios_df, intensities_zscores_d #' #' @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 zscore_pat_drugs_df: dataframe with Z-scores for all patient samples for drug metabolites #' @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 @@ -129,6 +135,7 @@ calculate_zscore_ratios <- function(metabolites_ratios_df, intensities_zscores_d make_and_save_violin_plot_pdfs <- function( zscore_patients_df, zscore_controls_df, + zscore_pat_drugs_df, path_metabolite_groups, nr_plots_perpage, number_of_samples, @@ -175,6 +182,13 @@ make_and_save_violin_plot_pdfs <- function( if (grepl("Diagnost", pdf_dir)) { # make list of metabolites that exceed alarm values for this patient top_metabs_patient <- get_top_metabolites_df(patient_id, dims_helix_table) + # for drugs: make list of top highest and lowest Z-scores for this patient + top_drugs_patient <- prepare_toplist( + patient_id, + zscore_pat_drugs_df, + number_of_metabolites$highest, + number_of_metabolites$lowest + ) } else { # make list of top highest and lowest Z-scores for this patient top_metabs_patient <- prepare_toplist( @@ -183,9 +197,16 @@ make_and_save_violin_plot_pdfs <- function( number_of_metabolites$highest, number_of_metabolites$lowest ) + # for drugs: make list of top highest and lowest Z-scores for this patient + top_drugs_patient <- prepare_toplist( + patient_id, + zscore_pat_drugs_df, + number_of_metabolites$highest, + number_of_metabolites$lowest + ) } # generate normal violin plots - create_pdf_violin_plots(pdf_dir, patient_id, metab_perpage, top_metabs_patient, explanation_violin_plot) + create_pdf_violin_plots(pdf_dir, patient_id, metab_perpage, top_metabs_patient, top_drugs_patient, explanation_violin_plot) } } } @@ -597,8 +618,9 @@ prepare_toplist <- function(patient_id, zscore_patients, num_of_highest_metaboli #' @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 top_drugs_pt: dataframe with increased and decreased drug 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) { +create_pdf_violin_plots <- function(pdf_dir, patient_id, metab_perpage, top_metab_pt, top_drugs_patient, explanation) { # set parameters for plots plot_height <- 9.6 plot_width <- 6 @@ -684,13 +706,41 @@ create_pdf_violin_plots <- function(pdf_dir, patient_id, metab_perpage, top_meta suppressWarnings(print(ggplot_object)) } + # put table of drugs into PDF file, if not empty + if (!is.null(dim(top_drugs_patient))) { + max_rows_per_page <- 35 + total_rows <- nrow(top_drugs_patient) + number_of_pages <- ceiling(total_rows / max_rows_per_page) + + # 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")) + ) + + 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_drugs_patient[start_row:end_row, ] + + table_grob <- tableGrob(page_data, theme = table_theme, rows = NULL) + + grid.arrange( + table_grob, + top = paste0("Top deviating drug metabolites for patient: ", patient_id) + ) + } + } + # 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) + if (grepl("Diagnost", pdf_dir)) { + 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) + } } } From 4f6c4321f91a80039c490a37328e2c833c3d2ba5 Mon Sep 17 00:00:00 2001 From: Mia Pras-Raves Date: Mon, 27 Jul 2026 11:25:59 +0200 Subject: [PATCH 10/10] modified name of input object in DIMS/GenerateQCOutput.R --- DIMS/GenerateQCOutput.R | 16 ++++++++-------- 1 file changed, 8 insertions(+), 8 deletions(-) diff --git a/DIMS/GenerateQCOutput.R b/DIMS/GenerateQCOutput.R index d22c337..269d82e 100644 --- a/DIMS/GenerateQCOutput.R +++ b/DIMS/GenerateQCOutput.R @@ -22,7 +22,7 @@ source(paste0(export_scripts_dir, "generate_qc_output_functions.R")) # load init files load(init_file) -# load outlist from GenerateExcel +# load outlist_filtered_zscores from GenerateExcel load("outlist.RData") # load combined adducts for each scanmodus load("AdductSums_positive.RData") @@ -39,10 +39,10 @@ control_label <- "C" #### CHECK NUMBER OF CONTROLS #### file_name <- "Check_number_of_controls.txt" min_num_controls <- 25 -check_number_of_controls(outlist, min_num_controls, file_name) +check_number_of_controls(outlist_filtered_zscores, min_num_controls, file_name) #### INTERNAL STANDARDS #### -is_list <- outlist[grep("Internal standard", outlist[, "relevance"], fixed = TRUE), ] +is_list <- outlist_filtered_zscores[grep("Internal standard", outlist_filtered_zscores[, "relevance"], fixed = TRUE), ] is_codes <- rownames(is_list) # check if there is data present for all the samples that the pipeline started with @@ -268,7 +268,7 @@ save_internal_standard_plot( # 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) -column_list <- colnames(outlist) +column_list <- colnames(outlist_filtered_zscores) patterns <- c("^(P1002\\.)[[:digit:]]+_", "^(P1003\\.)[[:digit:]]+_", "^(P1005\\.)[[:digit:]]+_") positive_controls_index <- grepl(pattern = paste(patterns, collapse = "|"), column_list) positive_control_list <- column_list[positive_controls_index] @@ -304,21 +304,21 @@ if (length(positive_control_list) > 0) { pa_sample_name <- positive_control_list[grepl(pos_ctrl_samplename, positive_control_list)] pa_codes <- c("HMDB0000824", "HMDB0000725", "HMDB0000123") pa_names <- c("Propionylcarnitine", "Propionylglycine", "Glycine") - pa_data <- get_pos_ctrl_data(outlist, pa_sample_name, pa_codes, pa_names) + pa_data <- get_pos_ctrl_data(outlist_filtered_zscores, pa_sample_name, pa_codes, pa_names) positive_control <- rbind(positive_control, pa_data) } if (any(grepl("^P1003", pos_ctrl))) { pku_sample_name <- positive_control_list[grepl(pos_ctrl_samplename, positive_control_list)] pku_codes <- c("HMDB0000159") pku_names <- c("L-Phenylalanine") - pku_data <- get_pos_ctrl_data(outlist, pku_sample_name, pku_codes, pku_names) + pku_data <- get_pos_ctrl_data(outlist_filtered_zscores, pku_sample_name, pku_codes, pku_names) positive_control <- rbind(positive_control, pku_data) } if (any(grepl("^P1005", pos_ctrl))) { lpi_sample_name <- positive_control_list[grepl(pos_ctrl_samplename, positive_control_list)] lpi_codes <- c("HMDB0000904", "HMDB0000641", "HMDB0000182") lpi_names <- c("Citrulline", "L-Glutamine", "L-Lysine") - lpi_data <- get_pos_ctrl_data(outlist, lpi_sample_name, lpi_codes, lpi_names) + lpi_data <- get_pos_ctrl_data(outlist_filtered_zscores, lpi_sample_name, lpi_codes, lpi_names) positive_control <- rbind(positive_control, lpi_data) } } @@ -354,7 +354,7 @@ is_pos_intensities <- get_is_intensities(outlist_tot_pos, is_codes = is_codes) # 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_metabolites_df <- outlist_filtered_zscores %>% filter(HMDB_code %in% sst_components$HMDB_ID) sst_sample_column_index <- grep("P1001", colnames(sst_metabolites_df)) # Check if SST mix sample(s) are present