isotope correction updates

R/correct-isotope.R

@@ -1,7 +1,7 @@
-#' Manual correction for interference contributed by another feature
+#' Manual interference correction
 #'
 #' @description
-#' The interference is substracted using following formula:
+#' The interference is subtracted using following formula:
 #' \deqn{Value_{Corrected} = Value_{Feature} - Factor_{Contribution} * Value_{Interfering Feature}}
 #'
 #' @param data MidarExperiment object
@@ -15,20 +15,24 @@
 
 #  Example:  mexp <- correct_interference_manually(mexp, "feature_intensity", "PC 32:0 | SM 36:1 M+3", "SM 36:1", 0.0106924, "PC 32:0")
 
-correct_interference_manually <- function(data = NULL, variable, feature, interfering_feature, relative_contribution, updated_feature_id = NULL) {
+correct_interference_manually <- function(data = NULL, variable, feature, interfering_feature, relative_contribution, updated_feature_id = NA) {
 
   check_data(data)
-
   variable_var <- rlang::ensym(variable)
 
-  updated_feature_id <- ifelse(is.null(updated_feature_id) | is.na(updated_feature_id), "", updated_feature_id)
+  updated_feature_id <- ifelse(is.null(updated_feature_id) | is.na(updated_feature_id), NA_character_, updated_feature_id)
 
   # Validate input
-  if (!feature %in% data@annot_features$feature_id) cli::cli_abort("Feature is not present in the dataset")
-  if (!interfering_feature %in% data@annot_features$feature_id) cli::cli_abort("Interfering feature is not present in the dataset")
-  if (!variable %in% names(data@dataset)) stop(glue::glue("Variable `{variable` is not defined in the dataset"))
-  if (relative_contribution < 0 | relative_contribution >= 1) cli::cli_abort("`relative_contribution` must be between 0 and 1")
-  if (updated_feature_id %in% data@annot_features$feature_id) cli::cli_abort("Mew fFeature name must not present already be present in the dataset")
+  if (!feature %in% data@annot_features$feature_id) cli::cli_abort(col_red("Feature is not present in the dataset"))
+  if (!interfering_feature %in% data@annot_features$feature_id) cli::cli_abort(col_red("Interfering feature is not present in the dataset"))
+  if (!variable %in% names(data@dataset)) cli::cli_abort(col_red("Variable `{variable}` is not defined in the dataset"))
+  if (is.numeric(relative_contribution) & relative_contribution < 0 ) cli::cli_abort(col_red("col_red(`relative_contribution` must be between larger than 0"))
+  if (!is.na(updated_feature_id) & updated_feature_id %in% data@annot_features$feature_id) cli::cli_abort(col_red("New feature id must not present in the dataset"))
+
+  if (!"interference_corrected" %in% names(data@dataset)) {
+    data@dataset <- data@dataset |>
+      mutate(interference_corrected = FALSE)
+  }
 
   # Correction
   data@dataset <- data@dataset |>
@@ -41,8 +45,10 @@ correct_interference_manually <- function(data = NULL, variable, feature, interf
       interference_corrected = if_else(.data$feature_id == feature, TRUE, .data$interference_corrected)
     )
 
-  data@dataset <- data@dataset |>
-    mutate(feature_id = if_else(.data$feature_id == feature & nchar(stringr::str_squish(.data$updated_feature_id)) > 0, updated_feature_id, .data$feature_id))
+  if(!is.na(updated_feature_id)) {
+    data@dataset <- data@dataset |>
+      mutate(feature_id = if_else(.data$feature_id == feature & .data$interference_corrected, updated_feature_id, .data$feature_id))
+  }
 
   data
 }

tests/testthat/test-correct-isotope.R

@@ -1,6 +1,6 @@
-mexp <- readRDS(file = testthat::test_path("testdata/MHQuant_demo.rds"))
-mexp_proc <- mexp
+mexp_orig <- readRDS(file = testthat::test_path("testdata/MHQuant_demo.rds"))
 
+mexp <- mexp_orig
 
 test_that("correct_interferences corrects overlapping interferences", {
 
@@ -47,3 +47,60 @@ test_that("correct_interferences corrects overlapping interferences", {
                7256.0500353124)
 
 })
+
+
+
+
+
+test_that("correct_interferences corrects overlapping interferences", {
+
+  # d18:2 is interfering with d18:1, which in turn is interfering with d18:0
+  # the code does not correct for M+4 isotope interference
+
+
+  mexp2 <- mexp_orig
+
+
+  expect_equal(mexp2@dataset |>
+                 filter(feature_id == "S1P d18:1 [M>60]",
+                        analysis_id == "008_LTR_LTR01") |>
+                 pull(feature_intensity),
+               85299)
+
+
+  mexp2 <- correct_interference_manually(mexp2,
+                                         variable = "feature_intensity",
+                                         feature = "S1P d18:1 [M>60]",
+                                         interfering_feature = "S1P d18:2 [M>60]",
+                                         relative_contribution = 0.0315385)
+
+
+
+  expect_equal(mexp2@dataset |>
+                 filter(feature_id == "S1P d18:1 [M>60]",
+                        analysis_id == "008_LTR_LTR01") |>
+                 pull(feature_intensity),
+               84304.7172490)
+
+  # test renaming of feature and additional correction
+  mexp2 <- correct_interference_manually(mexp2,
+                                         variable = "feature_intensity",
+                                         feature = "S1P d18:0 [M>60]",
+                                         interfering_feature = "S1P d18:1 [M>60]",
+                                         relative_contribution = 0.0315872, updated_feature_id = "S1P d18:0 [M>60] corrected")
+
+  expect_equal(mexp2@dataset |>
+                 filter(feature_id == "S1P d18:0 [M>60] corrected",
+                        analysis_id == "008_LTR_LTR01") |>
+                 pull(feature_intensity),
+               7256.0500353124)
+
+
+  expect_error(  mexp2 <- correct_interference_manually(mexp2,
+                                                        variable = "feature_intensity",
+                                                        feature = "S1P d18:0 [M>601]",
+                                                        interfering_feature = "S1P d18:1 [M>60]",
+                                                        relative_contribution = 0.0315872, updated_feature_id = "S1P d18:0 [M>60] corrected"),
+                 "Feature is not present in the dataset, please adjust `feature` argument"
+  )
+})

tests/testthat/test-metadata-import.R

@@ -216,6 +216,8 @@ test_that("Check import of inconsitent metadata", {
                "`is_quantifier` is inconsistently defined, i.e., not for one or more features. Please")
   expect_error(mexp <- midar:::import_metadata_features(mexp, path = path, sheet = "Features_Inconsist_val"),
                "`valid_feature` is inconsistently defined, i.e., not for one or more features. Please")
+  expect_error(mexp <- midar:::import_metadata_features(mexp, path = path, sheet = "Features_inval_interf"))
+  expect_error(mexp <- midar:::import_metadata_features(mexp, path = path, sheet = "Features_miss_interf"))
 })
 
 test_that("Replacing specific undefined metadata", {