From cd95388b37e4d9b7b9c7dcb1be3f5b9b1acf0a6b Mon Sep 17 00:00:00 2001 From: Tony Wu Date: Mon, 31 Aug 2026 20:13:38 -0400 Subject: [PATCH 1/3] refactor(createTurnoverRatios): Speed up code to use data.table --- .Rbuildignore | 2 + .gitignore | 1 + DESCRIPTION | 1 - NAMESPACE | 9 +++- R/Utils.R | 4 ++ R/protein_turnover_ratio_helper.R | 77 +++++++++++++--------------- man/calculatePeptideWeights.Rd | 84 ------------------------------- man/kendall_monotonicity.Rd | 42 ++++++++++++++++ 8 files changed, 91 insertions(+), 129 deletions(-) delete mode 100644 man/calculatePeptideWeights.Rd diff --git a/.Rbuildignore b/.Rbuildignore index 91114bf..ea84996 100644 --- a/.Rbuildignore +++ b/.Rbuildignore @@ -1,2 +1,4 @@ ^.*\.Rproj$ ^\.Rproj\.user$ +^\.positai$ +^\.claude$ diff --git a/.gitignore b/.gitignore index f4f606b..209e83c 100644 --- a/.gitignore +++ b/.gitignore @@ -3,3 +3,4 @@ .RData .Ruserdata *.Rproj +.positai diff --git a/DESCRIPTION b/DESCRIPTION index b56850d..d7e65f5 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -35,7 +35,6 @@ Imports: stats, parallel, data.table, - tidyr, stringr, plotly, utils diff --git a/NAMESPACE b/NAMESPACE index 4ae76f3..a782199 100644 --- a/NAMESPACE +++ b/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) diff --git a/R/Utils.R b/R/Utils.R index a35cb5a..a2dae9d 100644 --- a/R/Utils.R +++ b/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" )) diff --git a/R/protein_turnover_ratio_helper.R b/R/protein_turnover_ratio_helper.R index 5a60576..ae2ea96 100644 --- a/R/protein_turnover_ratio_helper.R +++ b/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 diff --git a/man/calculatePeptideWeights.Rd b/man/calculatePeptideWeights.Rd deleted file mode 100644 index ea3fbba..0000000 --- a/man/calculatePeptideWeights.Rd +++ /dev/null @@ -1,84 +0,0 @@ -% Generated by roxygen2: do not edit by hand -% Please edit documentation in R/protein_turnover_ratio_helper.R -\name{calculatePeptideWeights} -\alias{calculatePeptideWeights} -\title{Calculate quality-based weights for peptide measurements} -\usage{ -calculatePeptideWeights( - data, - protein_col = "Protein", - peptide_col = "BaseSequence", - time_col = "TimeVal", - response_col = "H_frac", - light_intensity_col = "Light", - validity_threshold = 1.3, - top_n_peptides = NULL -) -} -\arguments{ -\item{data}{Data frame with peptide-level measurements (output from calculateTurnoverRatios)} - -\item{protein_col}{Character. Column containing protein identifiers. Default = "Protein"} - -\item{peptide_col}{Character. Column containing peptide identifiers. Default = "BaseSequence"} - -\item{time_col}{Character. Column containing timepoint values. Default = "TimeVal"} - -\item{response_col}{Character. Column containing response values (e.g., H_frac). Default = "H_frac"} - -\item{light_intensity_col}{Character. Column containing light channel intensity. Default = "Light"} - -\item{validity_threshold}{Numeric. Maximum allowed response value. Default = 1.3} - -\item{top_n_peptides}{Numeric. Top n peptides based on median light channel intensity. If NULL, no filtering (all get weight 1).} -} -\value{ -Input data frame with added columns: -\itemize{ -\item n_obs: Total number of observations for this peptide -\item k_obs: Number of non-zero timepoints where this peptide is detected -\item coverage_per_peptide: k_obs / n (per-peptide detection proportion) -\item p_protein: Protein-level mean detection rate across all its peptides -\item coverage_score: P(X <= k_obs | n, p_protein) — binomial CDF coverage score -\item light_intensity_score: 1 (no filter) or binary top-N indicator (per protein) -\item monotonicity_score: Kendall correlation (time vs response), floored at 0 -\item validity_flag: 0 if any invalid values, 1 otherwise -\item weight: Combined quality weight (product of all components) -} -} -\description{ -Calculates weights based on coverage, signal intensity, monotonicity, and data validity. -Designed for protein turnover data but applicable to any dose/time-response data. -} -\details{ -Coverage is scored via a binomial CDF: for each peptide, P(X <= k | n, p) where k is -the number of non-zero timepoints detected, n is the total non-zero timepoints in the -experiment, and p is the protein-level mean detection rate across all its peptides. -This penalizes peptides with unusually low coverage relative to the protein's norm. -} -\examples{ -\dontrun{ -# Calculate ratios first -ratios <- calculateTurnoverRatios(feature_data) - -# Add quality weights -ratios_weighted <- calculatePeptideWeights(ratios) - -# Inspect coverage diagnostics -ratios_weighted \%>\% - group_by(Protein, BaseSequence) \%>\% - slice(1) \%>\% - select(Protein, BaseSequence, k_obs, p_protein, coverage_score, monotonicity_score, weight) - -# Use with doseResponseFit -result <- doseResponseFit( - data = ratios_weighted, - weights = ratios_weighted$weight, - increasing = TRUE -) - -# Use stricter validity threshold -ratios_weighted_strict <- calculatePeptideWeights(ratios, validity_threshold = 1.0) -} - -} diff --git a/man/kendall_monotonicity.Rd b/man/kendall_monotonicity.Rd index 6ec1f61..0f65e80 100644 --- a/man/kendall_monotonicity.Rd +++ b/man/kendall_monotonicity.Rd @@ -4,14 +4,47 @@ \alias{kendall_monotonicity} \title{Kendall monotonicity score, robust to sparse or all-missing data} \usage{ +kendall_monotonicity(time, response) + kendall_monotonicity(time, response) } \arguments{ \item{time}{Numeric vector of timepoints.} \item{response}{Numeric vector of responses (same length as \code{time}).} + +\item{data}{Data frame with peptide-level measurements (output from calculateTurnoverRatios)} + +\item{protein_col}{Character. Column containing protein identifiers. Default = "Protein"} + +\item{peptide_col}{Character. Column containing peptide identifiers. Default = "BaseSequence"} + +\item{time_col}{Character. Column containing timepoint values. Default = "TimeVal"} + +\item{response_col}{Character. Column containing response values (e.g., H_frac). Default = "H_frac"} + +\item{light_intensity_col}{Character. Column containing light channel intensity. Default = "Light"} + +\item{validity_threshold}{Numeric. Maximum allowed response value. Default = 1.3} + +\item{top_n_peptides}{Numeric. Top n peptides based on median light channel intensity. If NULL, no filtering (all get weight 1).} } \value{ +A single numeric monotonicity score in [0, 1]. + +Input data frame with added columns: +\itemize{ +\item n_obs: Total number of observations for this peptide +\item k_obs: Number of non-zero timepoints where this peptide is detected +\item coverage_per_peptide: k_obs / n (per-peptide detection proportion) +\item p_protein: Protein-level mean detection rate across all its peptides +\item coverage_score: P(X <= k_obs | n, p_protein) — binomial CDF coverage score +\item light_intensity_score: 1 (no filter) or binary top-N indicator (per protein) +\item monotonicity_score: Kendall correlation (time vs response), floored at 0 +\item validity_flag: 0 if any invalid values, 1 otherwise +\item weight: Combined quality weight (product of all components) +} + A single numeric monotonicity score in [0, 1]. } \description{ @@ -20,5 +53,14 @@ with fewer than two finite (time, response) pairs as non-monotonic (score 0) instead of letting \code{cor(use = "complete.obs")} error on zero complete pairs. A zero-variance group (>= 2 points but constant) yields \code{NA} from \code{cor()}, which is also mapped to 0. + +Calculates weights based on coverage, signal intensity, monotonicity, and data validity. +Designed for protein turnover data but applicable to any dose/time-response data. +} +\details{ +Coverage is scored via a binomial CDF: for each peptide, P(X <= k | n, p) where k is +the number of non-zero timepoints detected, n is the total non-zero timepoints in the +experiment, and p is the protein-level mean detection rate across all its peptides. +This penalizes peptides with unusually low coverage relative to the protein's norm. } \keyword{internal} From b913771802e259579025a0383a0743181a2cb1a0 Mon Sep 17 00:00:00 2001 From: Tony Wu Date: Mon, 31 Aug 2026 20:59:54 -0400 Subject: [PATCH 2/3] fix docs, add unit tests --- R/protein_turnover_ratio_helper.R | 24 --- man/calculatePeptideWeights.Rd | 84 ++++++++ man/kendall_monotonicity.Rd | 42 ---- .../test-protein_turnover_ratio_helper.R | 192 ++++++++++++++++++ 4 files changed, 276 insertions(+), 66 deletions(-) create mode 100644 man/calculatePeptideWeights.Rd create mode 100644 tests/testthat/test-protein_turnover_ratio_helper.R diff --git a/R/protein_turnover_ratio_helper.R b/R/protein_turnover_ratio_helper.R index ae2ea96..871b651 100644 --- a/R/protein_turnover_ratio_helper.R +++ b/R/protein_turnover_ratio_helper.R @@ -259,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 diff --git a/man/calculatePeptideWeights.Rd b/man/calculatePeptideWeights.Rd new file mode 100644 index 0000000..ea3fbba --- /dev/null +++ b/man/calculatePeptideWeights.Rd @@ -0,0 +1,84 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/protein_turnover_ratio_helper.R +\name{calculatePeptideWeights} +\alias{calculatePeptideWeights} +\title{Calculate quality-based weights for peptide measurements} +\usage{ +calculatePeptideWeights( + data, + protein_col = "Protein", + peptide_col = "BaseSequence", + time_col = "TimeVal", + response_col = "H_frac", + light_intensity_col = "Light", + validity_threshold = 1.3, + top_n_peptides = NULL +) +} +\arguments{ +\item{data}{Data frame with peptide-level measurements (output from calculateTurnoverRatios)} + +\item{protein_col}{Character. Column containing protein identifiers. Default = "Protein"} + +\item{peptide_col}{Character. Column containing peptide identifiers. Default = "BaseSequence"} + +\item{time_col}{Character. Column containing timepoint values. Default = "TimeVal"} + +\item{response_col}{Character. Column containing response values (e.g., H_frac). Default = "H_frac"} + +\item{light_intensity_col}{Character. Column containing light channel intensity. Default = "Light"} + +\item{validity_threshold}{Numeric. Maximum allowed response value. Default = 1.3} + +\item{top_n_peptides}{Numeric. Top n peptides based on median light channel intensity. If NULL, no filtering (all get weight 1).} +} +\value{ +Input data frame with added columns: +\itemize{ +\item n_obs: Total number of observations for this peptide +\item k_obs: Number of non-zero timepoints where this peptide is detected +\item coverage_per_peptide: k_obs / n (per-peptide detection proportion) +\item p_protein: Protein-level mean detection rate across all its peptides +\item coverage_score: P(X <= k_obs | n, p_protein) — binomial CDF coverage score +\item light_intensity_score: 1 (no filter) or binary top-N indicator (per protein) +\item monotonicity_score: Kendall correlation (time vs response), floored at 0 +\item validity_flag: 0 if any invalid values, 1 otherwise +\item weight: Combined quality weight (product of all components) +} +} +\description{ +Calculates weights based on coverage, signal intensity, monotonicity, and data validity. +Designed for protein turnover data but applicable to any dose/time-response data. +} +\details{ +Coverage is scored via a binomial CDF: for each peptide, P(X <= k | n, p) where k is +the number of non-zero timepoints detected, n is the total non-zero timepoints in the +experiment, and p is the protein-level mean detection rate across all its peptides. +This penalizes peptides with unusually low coverage relative to the protein's norm. +} +\examples{ +\dontrun{ +# Calculate ratios first +ratios <- calculateTurnoverRatios(feature_data) + +# Add quality weights +ratios_weighted <- calculatePeptideWeights(ratios) + +# Inspect coverage diagnostics +ratios_weighted \%>\% + group_by(Protein, BaseSequence) \%>\% + slice(1) \%>\% + select(Protein, BaseSequence, k_obs, p_protein, coverage_score, monotonicity_score, weight) + +# Use with doseResponseFit +result <- doseResponseFit( + data = ratios_weighted, + weights = ratios_weighted$weight, + increasing = TRUE +) + +# Use stricter validity threshold +ratios_weighted_strict <- calculatePeptideWeights(ratios, validity_threshold = 1.0) +} + +} diff --git a/man/kendall_monotonicity.Rd b/man/kendall_monotonicity.Rd index 0f65e80..6ec1f61 100644 --- a/man/kendall_monotonicity.Rd +++ b/man/kendall_monotonicity.Rd @@ -4,47 +4,14 @@ \alias{kendall_monotonicity} \title{Kendall monotonicity score, robust to sparse or all-missing data} \usage{ -kendall_monotonicity(time, response) - kendall_monotonicity(time, response) } \arguments{ \item{time}{Numeric vector of timepoints.} \item{response}{Numeric vector of responses (same length as \code{time}).} - -\item{data}{Data frame with peptide-level measurements (output from calculateTurnoverRatios)} - -\item{protein_col}{Character. Column containing protein identifiers. Default = "Protein"} - -\item{peptide_col}{Character. Column containing peptide identifiers. Default = "BaseSequence"} - -\item{time_col}{Character. Column containing timepoint values. Default = "TimeVal"} - -\item{response_col}{Character. Column containing response values (e.g., H_frac). Default = "H_frac"} - -\item{light_intensity_col}{Character. Column containing light channel intensity. Default = "Light"} - -\item{validity_threshold}{Numeric. Maximum allowed response value. Default = 1.3} - -\item{top_n_peptides}{Numeric. Top n peptides based on median light channel intensity. If NULL, no filtering (all get weight 1).} } \value{ -A single numeric monotonicity score in [0, 1]. - -Input data frame with added columns: -\itemize{ -\item n_obs: Total number of observations for this peptide -\item k_obs: Number of non-zero timepoints where this peptide is detected -\item coverage_per_peptide: k_obs / n (per-peptide detection proportion) -\item p_protein: Protein-level mean detection rate across all its peptides -\item coverage_score: P(X <= k_obs | n, p_protein) — binomial CDF coverage score -\item light_intensity_score: 1 (no filter) or binary top-N indicator (per protein) -\item monotonicity_score: Kendall correlation (time vs response), floored at 0 -\item validity_flag: 0 if any invalid values, 1 otherwise -\item weight: Combined quality weight (product of all components) -} - A single numeric monotonicity score in [0, 1]. } \description{ @@ -53,14 +20,5 @@ with fewer than two finite (time, response) pairs as non-monotonic (score 0) instead of letting \code{cor(use = "complete.obs")} error on zero complete pairs. A zero-variance group (>= 2 points but constant) yields \code{NA} from \code{cor()}, which is also mapped to 0. - -Calculates weights based on coverage, signal intensity, monotonicity, and data validity. -Designed for protein turnover data but applicable to any dose/time-response data. -} -\details{ -Coverage is scored via a binomial CDF: for each peptide, P(X <= k | n, p) where k is -the number of non-zero timepoints detected, n is the total non-zero timepoints in the -experiment, and p is the protein-level mean detection rate across all its peptides. -This penalizes peptides with unusually low coverage relative to the protein's norm. } \keyword{internal} diff --git a/tests/testthat/test-protein_turnover_ratio_helper.R b/tests/testthat/test-protein_turnover_ratio_helper.R new file mode 100644 index 0000000..acf2ad5 --- /dev/null +++ b/tests/testthat/test-protein_turnover_ratio_helper.R @@ -0,0 +1,192 @@ +# Tests for calculateTurnoverRatios and its helpers. +# +# calculateTurnoverRatios was refactored from dplyr/tidyr to data.table; these +# tests pin the numeric contract (aggregation, pivot, H/L fractions, tracer +# normalization) so a future rewrite of the internals is verifiable. + +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 + ) + # P1 carries the two modified peptides, P2 only the third + 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() + # Two charge states per (protein, peptide, time, run, label) -> one row out + 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)) + + # max is the default: verify one cell directly + 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)) + # constants are re-keyed through parse_timepoint, so order of names is irrelevant + 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") + # protein-level shape: peptide column is the protein column (MSstatsShiny path) + 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", { + # Single run: with replicate runs calculatePeptideWeights counts rows rather + # than distinct timepoints in k_obs, which drives coverage_per_peptide above 1 + # and yields NaN weights out of pbinom(). That predates the data.table + # refactor and is not exercised here. + 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)) +}) From 7d508823433509f127897ab06a9a8d3aaf66995f Mon Sep 17 00:00:00 2001 From: Tony Wu Date: Mon, 31 Aug 2026 21:01:35 -0400 Subject: [PATCH 3/3] remove unnecessary comments --- tests/testthat/test-protein_turnover_ratio_helper.R | 13 ------------- 1 file changed, 13 deletions(-) diff --git a/tests/testthat/test-protein_turnover_ratio_helper.R b/tests/testthat/test-protein_turnover_ratio_helper.R index acf2ad5..55b7147 100644 --- a/tests/testthat/test-protein_turnover_ratio_helper.R +++ b/tests/testthat/test-protein_turnover_ratio_helper.R @@ -1,8 +1,4 @@ # Tests for calculateTurnoverRatios and its helpers. -# -# calculateTurnoverRatios was refactored from dplyr/tidyr to data.table; these -# tests pin the numeric contract (aggregation, pivot, H/L fractions, tracer -# normalization) so a future rewrite of the internals is verifiable. make_feature_data <- function() { grid <- expand.grid( @@ -14,7 +10,6 @@ make_feature_data <- function() { CHARGE = c(2, 3), stringsAsFactors = FALSE ) - # P1 carries the two modified peptides, P2 only the third grid <- grid[!(grid$PROTEIN == "P2" & grid$PEPTIDE != "THIRDPEP"), ] grid <- grid[!(grid$PROTEIN == "P1" & grid$PEPTIDE == "THIRDPEP"), ] grid$INTENSITY <- seq_len(nrow(grid)) * 100 @@ -63,7 +58,6 @@ test_that("fractions sum to one and Total is the H + L sum", { test_that("duplicate features are collapsed by agg_function", { feat <- make_feature_data() - # Two charge states per (protein, peptide, time, run, label) -> one row out n_expected <- nrow(unique(feat[, c("PROTEIN", "PEPTIDE", "GROUP", "RUN")])) res_max <- calculateTurnoverRatios(feat) expect_equal(nrow(res_max), n_expected) @@ -71,7 +65,6 @@ test_that("duplicate features are collapsed by agg_function", { res_min <- calculateTurnoverRatios(feat, agg_function = min) expect_true(all(res_max$Heavy >= res_min$Heavy)) - # max is the default: verify one cell directly 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" & @@ -128,7 +121,6 @@ test_that("tracer normalization rescales H_frac and floors L_frac at zero", { tracer_constants = tc) expect_true("tracer_factor" %in% names(norm)) - # constants are re-keyed through parse_timepoint, so order of names is irrelevant 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)) @@ -147,7 +139,6 @@ test_that("custom column names and heavy/light labels are honored", { names(feat)[names(feat) == "PROTEIN"] <- "Protein" names(feat)[names(feat) == "GROUP"] <- "Condition" feat$LABEL <- ifelse(feat$LABEL == "H", "Heavy", "Light") - # protein-level shape: peptide column is the protein column (MSstatsShiny path) res <- calculateTurnoverRatios( feat, time_col = "Condition", peptide_col = "Protein", protein_col = "Protein", heavy_label = "Heavy", light_label = "Light" @@ -178,10 +169,6 @@ test_that("kendall_monotonicity floors at zero and tolerates sparse input", { }) test_that("calculatePeptideWeights consumes calculateTurnoverRatios output", { - # Single run: with replicate runs calculatePeptideWeights counts rows rather - # than distinct timepoints in k_obs, which drives coverage_per_peptide above 1 - # and yields NaN weights out of pbinom(). That predates the data.table - # refactor and is not exercised here. feat <- make_feature_data() res <- calculateTurnoverRatios(feat[feat$RUN == "R1", ]) weighted <- calculatePeptideWeights(res)