diff --git a/DIMS/SumAdducts.R b/DIMS/SumAdducts.R index 489ae17..523e2a9 100644 --- a/DIMS/SumAdducts.R +++ b/DIMS/SumAdducts.R @@ -11,12 +11,8 @@ source(paste0(scripts_dir, "sum_intensities_adducts.R")) if (grepl("positive_hmdb", hmdbpart_main_file)) { scanmode <- "positive" - # for the adduct sum: include adducts M+Na (1) and M+K (2) - adducts <- c(1, 2) } else if (grepl("negative_hmdb", hmdbpart_main_file)) { scanmode <- "negative" - # for the adduct sum: include adduct M+Cl (1) - adducts <- c(1) } # load input files @@ -28,5 +24,5 @@ hmdb_main_part <- get(load(hmdbpart_main_file)) batch_number <- strsplit(basename(hmdbpart_main_file), ".", fixed = TRUE)[[1]][2] # sum adducts and save output -adductsum <- sum_intensities_adducts(peakgroup_list, hmdb_main_part, adducts, z_score) +adductsum <- sum_intensities_adducts_perclass(peakgroup_list, hmdb_main_part, z_score) save(adductsum, file = paste(scanmode, "_", batch_number, "_SummedAdducts.RData", sep = "")) diff --git a/DIMS/preprocessing/sum_intensities_adducts.R b/DIMS/preprocessing/sum_intensities_adducts.R index bb5c3ab..5f1a44f 100644 --- a/DIMS/preprocessing/sum_intensities_adducts.R +++ b/DIMS/preprocessing/sum_intensities_adducts.R @@ -1,13 +1,14 @@ -sum_intensities_adducts <- function(peakgroup_list, hmdb_part, adducts, z_score) { +sum_intensities_adducts_perclass <- function(peakgroup_list, hmdb_part, 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"]) + hmdb_part_info <- cbind(HMDB_id = rownames(hmdb_part), + CompoundName = hmdb_part[, "CompoundName"], + adducts = hmdb_part[, "adducts"]) # create overview of row indices for each metabolite_adduct combination in peaklist # split the all_hmdb_ids column into list with each id as a value @@ -15,7 +16,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)), @@ -49,12 +50,14 @@ sum_intensities_adducts <- function(peakgroup_list, hmdb_part, adducts, z_score) for (hmdb_index in 1:nrow(hmdb_part_info)) { compound <- hmdb_part_info[hmdb_index, "HMDB_id"] + adducts <- hmdb_part_info[hmdb_index, "adducts"] + adducts <- as.numeric(strsplit(adducts, ",")[[1]]) compound_plus_adducts <- 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_id %in% compound_plus_adducts) metab_indices <- as.numeric(hmdb_in_peaklist_rownr$row_id[metab_row]) - + # find intensities and sum them ints <- peakgroup_list[metab_indices, int_cols] total <- apply(ints, 2, sum) @@ -66,7 +69,7 @@ sum_intensities_adducts <- function(peakgroup_list, hmdb_part, adducts, z_score) names_long <- c(names_long, hmdb_part_info[hmdb_index, "CompoundName"]) } } - + if (!is.null(adductsum)) { rownames(adductsum) <- names adductsum <- cbind(adductsum, "HMDB_name" = names_long) @@ -75,7 +78,7 @@ sum_intensities_adducts <- function(peakgroup_list, hmdb_part, adducts, z_score) hmdb_info <- hmdb_part[names, cols_hmdb_info, drop = FALSE] adductsum <- cbind(adductsum, hmdb_info) } - + return(adductsum) } diff --git a/DIMS/tests/testthat/test_sum_intensities_adducts.R b/DIMS/tests/testthat/test_sum_intensities_adducts.R index fcea89f..3c6fb79 100644 --- a/DIMS/tests/testthat/test_sum_intensities_adducts.R +++ b/DIMS/tests/testthat/test_sum_intensities_adducts.R @@ -1,9 +1,11 @@ # unit tests for SumAdducts -# function: sum_intensities_adducts +# function: sum_intensities_adducts_perclass source("../../preprocessing/sum_intensities_adducts.R") -# test sum_intensities_adducts -testthat::test_that("adduct sums are correctly calculated", { +suppressPackageStartupMessages(library("dplyr")) + +# test sum_intensities_adducts_perclass +testthat::test_that("adduct sums per class are correctly calculated", { # create peakgroup_list to test on in diagnostics setting test_peakgroup_list <- data.frame(matrix(NA, nrow = 4, ncol = 23)) colnames(test_peakgroup_list) <- c("mzmed.pgrp", "nrsamples", "ppmdev", "assi_HMDB", "all_hmdb_names", @@ -20,20 +22,21 @@ testthat::test_that("adduct sums are correctly calculated", { test_peakgroup_list[ , "assi_HMDB"], sep=";") test_peakgroup_list[ , grep("C", colnames(test_peakgroup_list))] <- 1000*(1:16) test_peakgroup_list[ , grep("P", colnames(test_peakgroup_list))] <- 10000*(1:16) - + # create HMDB part object - test_hmdb_main_part <- matrix(NA, nrow = 2, ncol = 8) + test_hmdb_main_part <- matrix(NA, nrow = 2, ncol = 9) colnames(test_hmdb_main_part) <- c("HMDB_ID_all", "sec_HMDB_ID", "CompoundName", "HMDB_name_all", - "Composition", "MNeutral", "MNeg", "Mpos") + "Composition", "MNeutral", "MNeg", "Mpos", "adducts") rownames(test_hmdb_main_part) <- c("HMDB1234567", "HMDB7654321") test_hmdb_main_part[, 1] <- c("HMDB1234567;HMDB0000567", "HMDB7654321;HMDB0000321") test_hmdb_main_part[, "CompoundName"] <- c("metab_1", "metab_2") - - expect_type(sum_intensities_adducts(test_peakgroup_list, test_hmdb_main_part, c(1, 2), 1), "character") - expect_equal(rownames(sum_intensities_adducts(test_peakgroup_list, test_hmdb_main_part, c(1, 2), 1))[1], "HMDB1234567", TRUE) - expect_equal(colnames(sum_intensities_adducts(test_peakgroup_list, test_hmdb_main_part, c(1, 2), 1)), - c("C101.1", "C102.1", "P2.1", "P3.1", "HMDB_name", "HMDB_ID_all", "sec_HMDB_ID", "HMDB_name_all"), TRUE) - expect_equal(as.vector(sum_intensities_adducts(test_peakgroup_list, test_hmdb_main_part, c(1, 2), 1)[1, c(1:4)]), - c("6000", "18000", "60000", "180000"), TRUE) + test_hmdb_main_part[, "adducts"] <- c("1,2", "1,7") + + expect_type(sum_intensities_adducts_perclass(test_peakgroup_list, test_hmdb_main_part, 1), "character") + expect_equal(rownames(sum_intensities_adducts_perclass(test_peakgroup_list, test_hmdb_main_part, 1))[1], "HMDB1234567") + expect_equal(colnames(sum_intensities_adducts_perclass(test_peakgroup_list, test_hmdb_main_part, 1)), + c("C101.1", "C102.1", "P2.1", "P3.1", "HMDB_name", "HMDB_ID_all", "sec_HMDB_ID", "HMDB_name_all")) + expect_equal(as.vector(sum_intensities_adducts_perclass(test_peakgroup_list, test_hmdb_main_part, 1)[1, c(1:4)]), + c("6000", "18000", "60000", "180000")) })