refactor(createTurnoverRatios): Speed up code to use data.table

.Rbuildignore

@@ -1,2 +1,4 @@
 ^.*\.Rproj$
 ^\.Rproj\.user$
+^\.positai$
+^\.claude$

.gitignore

@@ -3,3 +3,4 @@
 .RData
 .Ruserdata
 *.Rproj
+.positai

DESCRIPTION

@@ -35,7 +35,6 @@ Imports:
     stats,
     parallel,
     data.table,
-    tidyr,
     stringr,
     plotly,
     utils

NAMESPACE

@@ -17,7 +17,14 @@ export(visualizeResponseProtein)
 import(dplyr)
 importFrom(BiocParallel,bplapply)
 importFrom(BiocParallel,bpparam)
+importFrom(data.table,":=")
+importFrom(data.table,as.data.table)
+importFrom(data.table,copy)
+importFrom(data.table,data.table)
+importFrom(data.table,dcast)
+importFrom(data.table,is.data.table)
 importFrom(data.table,rbindlist)
+importFrom(data.table,setnames)
 importFrom(dplyr,across)
 importFrom(dplyr,all_of)
 importFrom(dplyr,any_of)
@@ -35,7 +42,6 @@ importFrom(dplyr,n_distinct)
 importFrom(dplyr,pull)
 importFrom(dplyr,rename)
 importFrom(dplyr,select)
-importFrom(dplyr,slice_head)
 importFrom(dplyr,summarise)
 importFrom(dplyr,ungroup)
 importFrom(ggplot2,aes)
@@ -78,6 +84,5 @@ importFrom(stats,rnorm)
 importFrom(stringr,str_detect)
 importFrom(stringr,str_extract)
 importFrom(stringr,str_remove_all)
-importFrom(tidyr,pivot_wider)
 importFrom(utils,setTxtProgressBar)
 importFrom(utils,txtProgressBar)

R/Utils.R

@@ -14,6 +14,10 @@ utils::globalVariables(c(
   # Variables from plotHitRateMSstatsResponse
   "adj.pvalue",
 
+  # Variables from calculateTurnoverRatios (data.table NSE)
+  ".SD", "Protein", "BaseSequence", "Label", "TimeVal", "Run",
+  "Heavy", "Light", "Total", "H_frac", "L_frac", "tracer_factor",
+
   # Add any other NSE variables used in your package
   "drug", "protein", "dose", "x", "y", "y_pred"
 ))

R/protein_turnover_ratio_helper.R

@@ -63,8 +63,7 @@
 #' }
 #'
 #' @export
-#' @importFrom dplyr filter mutate group_by summarise arrange slice_head pull ungroup select rename
-#' @importFrom tidyr pivot_wider
+#' @importFrom data.table as.data.table copy data.table dcast is.data.table setnames :=
 #' @importFrom stringr str_remove_all str_extract str_detect
 calculateTurnoverRatios <- function(
     feature_data,
@@ -89,58 +88,55 @@ calculateTurnoverRatios <- function(
     stop("Missing required columns: ", paste(missing_cols, collapse = ", "))
   }
 
-  # Process all proteins
-  df <- feature_data %>%
-    mutate(
-      Protein = .data[[protein_col]],
-      BaseSequence = str_remove_all(.data[[peptide_col]], "\\[.*?\\]"),
-      Label = .data[[channel_col]],
-      TimeVal = parse_timepoint(.data[[time_col]]),
-      Intensity = .data[[intensity_col]],
-      Run = .data[[run_col]]
-    ) %>%
-    filter(!is.na(TimeVal)) %>%
-    filter(Label %in% c(heavy_label, light_label))
+  # Process all proteins. Copy so the caller's data.table is never modified
+  # by the := assignments below.
+  df <- if (is.data.table(feature_data)) {
+    copy(feature_data)
+  } else {
+    as.data.table(feature_data)
+  }
+
+  df[, c("Protein", "BaseSequence", "Label",
+         "TimeVal", "Intensity", "Run") := list(
+    df[[protein_col]],
+    str_remove_all(df[[peptide_col]], "\\[.*?\\]"),
+    df[[channel_col]],
+    parse_timepoint(df[[time_col]]),
+    df[[intensity_col]],
+    df[[run_col]]
+  )]
+  df <- df[!is.na(TimeVal) & Label %in% c(heavy_label, light_label)]
 
   if (nrow(df) == 0) {
     warning("No data found matching Heavy/Light labels")
-    return(data.frame())
+    return(data.table())
   }
 
-  # Apply optional peptide selector per protein
+  # Apply optional peptide selector per protein (.SD excludes Protein)
   if (!is.null(peptide_selector)) {
-    df <- df %>%
-      group_by(Protein) %>%
-      group_modify(~ peptide_selector(.x)) %>%
-      ungroup()
+    df <- df[, as.data.table(peptide_selector(.SD)), by = Protein]
   }
 
   # Aggregate duplicates (multiple features/transitions for same peptide)
   # This can happen when there are multiple charge states or modifications
-  df <- df %>%
-    group_by(Protein, BaseSequence, TimeVal, Run, Label) %>%
-    summarise(Intensity = agg_function(Intensity, na.rm = TRUE), .groups = "drop")
+  df <- df[, list(Intensity = agg_function(Intensity, na.rm = TRUE)),
+           keyby = list(Protein, BaseSequence, TimeVal, Run, Label)]
 
   # Pivot to wide format (keep replicates separate)
-  df_wide <- df %>%
-    select(Protein, BaseSequence, TimeVal, Run, Label, Intensity) %>%
-    pivot_wider(names_from = Label, values_from = Intensity)
+  df_wide <- dcast(df, Protein + BaseSequence + TimeVal + Run ~ Label,
+                   value.var = "Intensity")
 
   # Rename heavy/light columns if they're not "Heavy" and "Light"
   if (heavy_label != "Heavy") {
-    df_wide <- df_wide %>% rename(Heavy = all_of(heavy_label))
+    setnames(df_wide, heavy_label, "Heavy")
   }
   if (light_label != "Light") {
-    df_wide <- df_wide %>% rename(Light = all_of(light_label))
+    setnames(df_wide, light_label, "Light")
   }
 
-  df_wide <- df_wide %>%
-    filter(!is.na(Heavy) & !is.na(Light)) %>%
-    mutate(
-      Total = Heavy + Light,
-      H_frac = Heavy / Total,
-      L_frac = Light / Total
-    )
+  df_wide <- df_wide[!is.na(Heavy) & !is.na(Light)]
+  df_wide[, Total := Heavy + Light]
+  df_wide[, c("H_frac", "L_frac") := list(Heavy / Total, Light / Total)]
 
   # Optional tracer normalization
   if (normalize_tracer) {
@@ -150,15 +146,12 @@ calculateTurnoverRatios <- function(
 
     names(tracer_constants) = as.character(parse_timepoint(names(tracer_constants)))
 
-    df_wide <- df_wide %>%
-      mutate(
-        tracer_factor = tracer_constants[as.character(TimeVal)],
-        H_frac = H_frac / tracer_factor,
-        L_frac =  pmax(1 - H_frac, 0)
-      )
+    df_wide[, tracer_factor := unname(tracer_constants[as.character(TimeVal)])]
+    df_wide[, H_frac := H_frac / tracer_factor]
+    df_wide[, L_frac := pmax(1 - H_frac, 0)]
   }
 
-  df_wide
+  df_wide[]
 }
 
 #' Parse timepoint strings to numeric hours
@@ -266,30 +259,6 @@ kendall_monotonicity <- function(time, response) {
 #' ratios_weighted_strict <- calculatePeptideWeights(ratios, validity_threshold = 1.0)
 #' }
 #'
-#' Kendall monotonicity score, robust to sparse or all-missing data
-#'
-#' Computes `max(0, Kendall's tau)` between time and response, treating a group
-#' with fewer than two finite (time, response) pairs as non-monotonic (score 0)
-#' instead of letting `cor(use = "complete.obs")` error on zero complete pairs.
-#' A zero-variance group (>= 2 points but constant) yields `NA` from `cor()`,
-#' which is also mapped to 0.
-#'
-#' @param time Numeric vector of timepoints.
-#' @param response Numeric vector of responses (same length as `time`).
-#'
-#' @return A single numeric monotonicity score in \[0, 1\].
-#'
-#' @keywords internal
-#' @importFrom stats cor
-kendall_monotonicity <- function(time, response) {
-  ok <- is.finite(time) & is.finite(response)
-  if (sum(ok) < 2) {
-    return(0)
-  }
-  score <- suppressWarnings(cor(time[ok], response[ok], method = "kendall"))
-  if (is.na(score)) 0 else max(0, score)
-}
-
 #' @export
 #' @importFrom dplyr group_by mutate ungroup across all_of distinct summarise left_join if_else dense_rank
 #' @importFrom stats cor pbinom median

tests/testthat/test-protein_turnover_ratio_helper.R

@@ -0,0 +1,179 @@
+# Tests for calculateTurnoverRatios and its helpers.
+
+make_feature_data <- function() {
+  grid <- expand.grid(
+    PROTEIN = c("P1", "P2"),
+    PEPTIDE = c("PEP[+80]TIDE", "SEQVENCE", "THIRDPEP"),
+    GROUP = c("0hr", "1hr", "12hrs", "168hrs"),
+    RUN = c("R1", "R2"),
+    LABEL = c("H", "L"),
+    CHARGE = c(2, 3),
+    stringsAsFactors = FALSE
+  )
+  grid <- grid[!(grid$PROTEIN == "P2" & grid$PEPTIDE != "THIRDPEP"), ]
+  grid <- grid[!(grid$PROTEIN == "P1" & grid$PEPTIDE == "THIRDPEP"), ]
+  grid$INTENSITY <- seq_len(nrow(grid)) * 100
+  rownames(grid) <- NULL
+  grid
+}
+
+test_that("parse_timepoint converts hours, days and weeks to hours", {
+  expect_equal(parse_timepoint(c("0hr", "1hr", "12hrs", "168hrs")),
+               c(0, 1, 12, 168))
+  expect_equal(parse_timepoint(c("2d", "1day")), c(48, 24))
+  expect_equal(parse_timepoint(c("1w", "2weeks")), c(168, 336))
+  expect_true(is.na(parse_timepoint("control")))
+})
+
+test_that("calculateTurnoverRatios errors on missing required columns", {
+  feat <- make_feature_data()
+  expect_error(
+    calculateTurnoverRatios(feat[, setdiff(names(feat), "RUN")]),
+    "Missing required columns: RUN"
+  )
+})
+
+test_that("calculateTurnoverRatios returns the documented columns", {
+  res <- calculateTurnoverRatios(make_feature_data())
+  expect_true(data.table::is.data.table(res))
+  expect_equal(names(res),
+               c("Protein", "BaseSequence", "TimeVal", "Run",
+                 "Heavy", "Light", "Total", "H_frac", "L_frac"))
+  expect_gt(nrow(res), 0)
+})
+
+test_that("modifications are stripped from the peptide sequence", {
+  res <- calculateTurnoverRatios(make_feature_data())
+  expect_setequal(unique(res$BaseSequence),
+                  c("PEPTIDE", "SEQVENCE", "THIRDPEP"))
+  expect_false(any(grepl("[", res$BaseSequence, fixed = TRUE)))
+})
+
+test_that("fractions sum to one and Total is the H + L sum", {
+  res <- calculateTurnoverRatios(make_feature_data())
+  expect_equal(res$Total, res$Heavy + res$Light)
+  expect_equal(res$H_frac + res$L_frac, rep(1, nrow(res)))
+  expect_equal(res$H_frac, res$Heavy / res$Total)
+})
+
+test_that("duplicate features are collapsed by agg_function", {
+  feat <- make_feature_data()
+  n_expected <- nrow(unique(feat[, c("PROTEIN", "PEPTIDE", "GROUP", "RUN")]))
+  res_max <- calculateTurnoverRatios(feat)
+  expect_equal(nrow(res_max), n_expected)
+
+  res_min <- calculateTurnoverRatios(feat, agg_function = min)
+  expect_true(all(res_max$Heavy >= res_min$Heavy))
+
+  cell <- feat[feat$PROTEIN == "P1" & feat$PEPTIDE == "SEQVENCE" &
+                 feat$GROUP == "1hr" & feat$RUN == "R1" & feat$LABEL == "H", ]
+  target <- res_max[res_max$Protein == "P1" & res_max$BaseSequence == "SEQVENCE" &
+                      res_max$TimeVal == 1 & res_max$Run == "R1", ]
+  expect_equal(target$Heavy, max(cell$INTENSITY))
+})
+
+test_that("unparseable timepoints and non-H/L labels are dropped", {
+  feat <- make_feature_data()
+  noise <- rbind(
+    transform(feat[1:4, ], GROUP = "control"),
+    transform(feat[1:4, ], LABEL = "M")
+  )
+  res_clean <- calculateTurnoverRatios(feat)
+  res_noisy <- calculateTurnoverRatios(rbind(feat, noise))
+  expect_equal(as.data.frame(res_noisy), as.data.frame(res_clean))
+})
+
+test_that("rows without both channels are dropped", {
+  feat <- make_feature_data()
+  drop <- feat$PROTEIN == "P1" & feat$PEPTIDE == "SEQVENCE" &
+    feat$GROUP == "1hr" & feat$RUN == "R2" & feat$LABEL == "L"
+  res <- calculateTurnoverRatios(feat[!drop, ])
+  orphan <- res[res$Protein == "P1" & res$BaseSequence == "SEQVENCE" &
+                  res$TimeVal == 1 & res$Run == "R2", ]
+  expect_equal(nrow(orphan), 0)
+  expect_false(anyNA(res$Heavy))
+  expect_false(anyNA(res$Light))
+})
+
+test_that("no matching Heavy/Light labels warns and returns no rows", {
+  feat <- transform(make_feature_data(), LABEL = "X")
+  expect_warning(res <- calculateTurnoverRatios(feat),
+                 "No data found matching Heavy/Light labels")
+  expect_equal(nrow(res), 0)
+})
+
+test_that("peptide_selector is applied within each protein", {
+  feat <- make_feature_data()
+  keep_first <- function(df) {
+    df[df$BaseSequence == sort(unique(df$BaseSequence))[1], ]
+  }
+  res <- calculateTurnoverRatios(feat, peptide_selector = keep_first)
+  by_protein <- tapply(res$BaseSequence, res$Protein, function(x) unique(x))
+  expect_equal(unname(by_protein[["P1"]]), "PEPTIDE")
+  expect_equal(unname(by_protein[["P2"]]), "THIRDPEP")
+})
+
+test_that("tracer normalization rescales H_frac and floors L_frac at zero", {
+  feat <- make_feature_data()
+  tc <- c("0hr" = 1.0, "1hr" = 0.95, "12hrs" = 0.85, "168hrs" = 0.75)
+  plain <- calculateTurnoverRatios(feat)
+  norm <- calculateTurnoverRatios(feat, normalize_tracer = TRUE,
+                                  tracer_constants = tc)
+
+  expect_true("tracer_factor" %in% names(norm))
+  expect_equal(norm$tracer_factor, unname(tc[match(norm$TimeVal, c(0, 1, 12, 168))]))
+  expect_equal(norm$H_frac, plain$H_frac / norm$tracer_factor)
+  expect_equal(norm$L_frac, pmax(1 - norm$H_frac, 0))
+  expect_true(all(norm$L_frac >= 0))
+})
+
+test_that("tracer normalization requires constants", {
+  expect_error(
+    calculateTurnoverRatios(make_feature_data(), normalize_tracer = TRUE),
+    "tracer_constants was not provided"
+  )
+})
+
+test_that("custom column names and heavy/light labels are honored", {
+  feat <- make_feature_data()
+  names(feat)[names(feat) == "PROTEIN"] <- "Protein"
+  names(feat)[names(feat) == "GROUP"] <- "Condition"
+  feat$LABEL <- ifelse(feat$LABEL == "H", "Heavy", "Light")
+  res <- calculateTurnoverRatios(
+    feat, time_col = "Condition", peptide_col = "Protein",
+    protein_col = "Protein", heavy_label = "Heavy", light_label = "Light"
+  )
+  expect_true(all(c("Heavy", "Light") %in% names(res)))
+  expect_setequal(unique(res$BaseSequence), c("P1", "P2"))
+})
+
+test_that("a data.table input is not modified in place", {
+  feat <- data.table::as.data.table(make_feature_data())
+  before <- data.table::copy(feat)
+  calculateTurnoverRatios(feat)
+  expect_equal(feat, before)
+})
+
+test_that("data.frame and data.table inputs give the same result", {
+  feat <- make_feature_data()
+  expect_equal(calculateTurnoverRatios(feat),
+               calculateTurnoverRatios(data.table::as.data.table(feat)))
+})
+
+test_that("kendall_monotonicity floors at zero and tolerates sparse input", {
+  expect_equal(kendall_monotonicity(c(1, 2, 3), c(1, 2, 3)), 1)
+  expect_equal(kendall_monotonicity(c(1, 2, 3), c(3, 2, 1)), 0)
+  expect_equal(kendall_monotonicity(1, 1), 0)
+  expect_equal(kendall_monotonicity(c(1, 2), c(NA, NA)), 0)
+  expect_equal(kendall_monotonicity(c(1, 2, 3), c(5, 5, 5)), 0)
+})
+
+test_that("calculatePeptideWeights consumes calculateTurnoverRatios output", {
+  feat <- make_feature_data()
+  res <- calculateTurnoverRatios(feat[feat$RUN == "R1", ])
+  weighted <- calculatePeptideWeights(res)
+  expect_true(all(c("coverage_score", "monotonicity_score", "validity_flag",
+                    "weight") %in% names(weighted)))
+  expect_equal(nrow(weighted), nrow(res))
+  expect_true(all(weighted$weight >= 0 & weighted$weight <= 1))
+})