Feature/DIMS_adductsum_perclass

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 = ""))

DIMS/preprocessing/sum_intensities_adducts.R

@@ -1,21 +1,22 @@
-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
   hmdb_in_peaklist <- strsplit(peakgroup_list$all_hmdb_ids, ";")
   # 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)
 }
 

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"))
 })