From 4d0dc2197d93693b3b4d646a5a3ce75c02d84544 Mon Sep 17 00:00:00 2001 From: Tony Wu Date: Mon, 20 Jul 2026 15:43:05 -0400 Subject: [PATCH 1/2] fix(turnover): handle scenarios where L_frac is NA for all timepoints --- R/protein_turnover_ratio_helper.R | 31 ++++++++++++++++++++++++++----- 1 file changed, 26 insertions(+), 5 deletions(-) diff --git a/R/protein_turnover_ratio_helper.R b/R/protein_turnover_ratio_helper.R index 8f910c4..a33f3e6 100644 --- a/R/protein_turnover_ratio_helper.R +++ b/R/protein_turnover_ratio_helper.R @@ -242,6 +242,30 @@ parse_timepoint <- function(time_strings) { #' 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 @@ -297,11 +321,8 @@ calculatePeptideWeights <- function( ) %>% group_by(across(all_of(c(protein_col, peptide_col)))) %>% mutate( - monotonicity_score = pmax(0, - cor(.data[[time_col]], .data[[response_col]], - method = "kendall", use = "complete.obs") - ), - monotonicity_score = if_else(is.na(monotonicity_score), 0, monotonicity_score) + monotonicity_score = kendall_monotonicity(.data[[time_col]], + .data[[response_col]]) ) %>% ungroup() %>% mutate( From bd83ab3fd46c26c2783f7e4f9af95760810c3fb4 Mon Sep 17 00:00:00 2001 From: Tony Wu Date: Mon, 20 Jul 2026 15:48:25 -0400 Subject: [PATCH 2/2] edit docs --- DESCRIPTION | 2 +- NAMESPACE | 11 ++++ R/ExperimentalDesignSimulation.R | 3 +- R/TPR_Power_Curve.R | 6 +-- R/protein_turnover_ratio_helper.R | 48 ++++++++--------- man/DIA_MSstats_Normalized.Rd | 2 +- man/calculateConfidence.Rd | 90 +++++++++++++++++++++++++++++++ man/calculatePeptideWeights.Rd | 21 +++++--- man/calculateQCScore.Rd | 61 +++++++++++++++++++++ man/classifyTurnoverProteins.Rd | 77 ++++++++++++++++++++++++++ man/kendall_monotonicity.Rd | 24 +++++++++ 11 files changed, 307 insertions(+), 38 deletions(-) create mode 100644 man/calculateConfidence.Rd create mode 100644 man/calculateQCScore.Rd create mode 100644 man/classifyTurnoverProteins.Rd create mode 100644 man/kendall_monotonicity.Rd diff --git a/DESCRIPTION b/DESCRIPTION index d29519e..8250789 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -28,7 +28,6 @@ License: Artistic-2.0 Encoding: UTF-8 Depends: R (>= 4.5.0) LazyData: false -RoxygenNote: 7.3.3 Imports: BiocParallel, ggplot2, @@ -54,3 +53,4 @@ Suggests: VignetteBuilder: knitr Roxygen: list(markdown = TRUE) biocViews: Proteomics, MassSpectrometry, StatisticalMethod, Software, Regression +Config/roxygen2/version: 8.0.0 diff --git a/NAMESPACE b/NAMESPACE index de45ddb..4ae76f3 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -1,8 +1,11 @@ # Generated by roxygen2: do not edit by hand export(MSstatsPrepareDoseResponseFit) +export(calculateConfidence) export(calculatePeptideWeights) +export(calculateQCScore) export(calculateTurnoverRatios) +export(classifyTurnoverProteins) export(convertGroupToNumericDose) export(doseResponseFit) export(futureExperimentSimulation) @@ -17,12 +20,18 @@ importFrom(BiocParallel,bpparam) importFrom(data.table,rbindlist) importFrom(dplyr,across) importFrom(dplyr,all_of) +importFrom(dplyr,any_of) importFrom(dplyr,arrange) +importFrom(dplyr,case_when) +importFrom(dplyr,coalesce) +importFrom(dplyr,dense_rank) importFrom(dplyr,distinct) importFrom(dplyr,filter) importFrom(dplyr,group_by) importFrom(dplyr,if_else) +importFrom(dplyr,left_join) importFrom(dplyr,mutate) +importFrom(dplyr,n_distinct) importFrom(dplyr,pull) importFrom(dplyr,rename) importFrom(dplyr,select) @@ -59,7 +68,9 @@ importFrom(plotly,ggplotly) importFrom(plotly,layout) importFrom(stats,approx) importFrom(stats,cor) +importFrom(stats,median) importFrom(stats,p.adjust) +importFrom(stats,pbinom) importFrom(stats,pf) importFrom(stats,quantile) importFrom(stats,rlnorm) diff --git a/R/ExperimentalDesignSimulation.R b/R/ExperimentalDesignSimulation.R index 1b08c95..890fbf8 100644 --- a/R/ExperimentalDesignSimulation.R +++ b/R/ExperimentalDesignSimulation.R @@ -494,8 +494,7 @@ simulateChemoProteinLevelNonParametric = function(N_proteins = 3000, #' @param concentration_count Number of concentrations in simulation #' #' @return A list containing the plot and plot data -#' @importFrom ggplot2 ggplot aes geom_bar geom_text labs scale_fill_manual -#' scale_y_continuous theme_classic theme element_text +#' @importFrom ggplot2 ggplot aes geom_bar geom_text labs scale_fill_manual scale_y_continuous theme_classic theme element_text #' @import dplyr plotHitRateMSstatsResponse = function(results, rep_count, concentration_count) { diff --git a/R/TPR_Power_Curve.R b/R/TPR_Power_Curve.R index 8f066e6..3de844d 100644 --- a/R/TPR_Power_Curve.R +++ b/R/TPR_Power_Curve.R @@ -217,8 +217,7 @@ run_tpr_simulation <- function(rep_range, concentrations, dose_range, #' @param show_legend Logical. Whether to display the legend. #' #' @return A ggplot object. -#' @importFrom ggplot2 ggplot aes geom_line geom_point scale_x_continuous -#' scale_y_continuous scale_color_manual labs theme_bw theme element_text +#' @importFrom ggplot2 ggplot aes geom_line geom_point scale_x_continuous scale_y_continuous scale_color_manual labs theme_bw theme element_text #' @noRd .make_tpr_panel <- function(data, k_grid, show_legend = FALSE) { rep_levels <- sort(unique(data$N_rep)) @@ -288,8 +287,7 @@ run_tpr_simulation <- function(rep_range, concentrations, dose_range, #' plot_tpr_power_curve(results) #' } #' -#' @importFrom ggplot2 ggplot aes geom_line geom_point scale_x_continuous -#' scale_y_continuous scale_color_manual labs theme_bw theme element_text +#' @importFrom ggplot2 ggplot aes geom_line geom_point scale_x_continuous scale_y_continuous scale_color_manual labs theme_bw theme element_text #' @importFrom plotly ggplotly layout #' @export plot_tpr_power_curve <- function(simulation_results, static = FALSE) { diff --git a/R/protein_turnover_ratio_helper.R b/R/protein_turnover_ratio_helper.R index a33f3e6..0419159 100644 --- a/R/protein_turnover_ratio_helper.R +++ b/R/protein_turnover_ratio_helper.R @@ -187,6 +187,30 @@ parse_timepoint <- function(time_strings) { return(hours) } +#' 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) +} + #' Calculate quality-based weights for peptide measurements #' @@ -242,30 +266,6 @@ parse_timepoint <- function(time_strings) { #' 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/DIA_MSstats_Normalized.Rd b/man/DIA_MSstats_Normalized.Rd index 24ef2c5..a559c0f 100644 --- a/man/DIA_MSstats_Normalized.Rd +++ b/man/DIA_MSstats_Normalized.Rd @@ -9,7 +9,7 @@ A data frame with protein-level abundance values and associated MSstats metadata column names. } \usage{ -DIA_MSstats_Normalized +data(DIA_MSstats_Normalized) } \description{ This dataset contains normalized protein-level data from a DIA-MS diff --git a/man/calculateConfidence.Rd b/man/calculateConfidence.Rd new file mode 100644 index 0000000..4eda356 --- /dev/null +++ b/man/calculateConfidence.Rd @@ -0,0 +1,90 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/protein_turnover_ratio_helper.R +\name{calculateConfidence} +\alias{calculateConfidence} +\title{Calculate per-protein confidence score for turnover fits} +\usage{ +calculateConfidence( + weights_df, + fit_df, + qc_df, + feature_data, + protein_col = "Protein", + weight_col = "weight", + sse_col = "SSE_Full", + qc_protein_col = "PROTEIN", + qc_score_col = "qc_score", + feature_protein_col = "PROTEIN", + feature_peptide_col = "PEPTIDE", + feature_label_col = "LABEL", + heavy_label = "H", + k_shrinkage = 2 +) +} +\arguments{ +\item{weights_df}{Output of \code{calculatePeptideWeights()}.} + +\item{fit_df}{Output of \code{doseResponseFit()} (must contain an SSE column).} + +\item{qc_df}{Output of \code{calculateQCScore()}.} + +\item{feature_data}{Raw feature-level data (used to count heavy peptides per protein).} + +\item{protein_col}{Character. Column in weights_df / fit_df identifying proteins. Default = "Protein"} + +\item{weight_col}{Character. Column in weights_df containing per-observation weight. Default = "weight"} + +\item{sse_col}{Character. Column in fit_df containing SSE. Default = "SSE_Full"} + +\item{qc_protein_col}{Character. Protein column in qc_df. Default = "PROTEIN"} + +\item{qc_score_col}{Character. QC score column in qc_df. Default = "qc_score"} + +\item{feature_protein_col}{Character. Protein column in feature_data. Default = "PROTEIN"} + +\item{feature_peptide_col}{Character. Peptide column in feature_data. Default = "PEPTIDE"} + +\item{feature_label_col}{Character. Label column in feature_data. Default = "LABEL"} + +\item{heavy_label}{Character. Value indicating heavy channel. Default = "H"} + +\item{k_shrinkage}{Numeric. Bayesian shrinkage constant for the peptide-count factor. +Larger values penalize low-peptide proteins more strongly. Default = 2} +} +\value{ +The input fit_df with additional columns: +\itemize{ +\item mean_weight: average per-observation weight +\item n_obs: number of observations used in the fit +\item qc_score: joined from qc_df +\item n_heavy_peptides: count from feature_data +\item pep_factor: n_heavy / (n_heavy + k_shrinkage) +\item confidence: combined score in [0, 1] +} +} +\description{ +Combines peptide-quality weights, fit residuals, light-channel QC, and +a Bayesian shrinkage factor on heavy-peptide count into a single +per-protein confidence score in [0, 1]. Higher values indicate the +dose-response fit is supported by clean, complete, abundant data. +} +\details{ +Formula: + +confidence = mean_weight * 1/(1 + SSE_Full) * qc_score * n_heavy / (n_heavy + k_shrinkage) + +The \code{n_heavy / (n_heavy + k)} term is a Bayesian / Laplace smoothing +factor that penalizes thin peptide support: with \code{k_shrinkage = 2}, a +1-peptide protein is capped at 1/3 of its otherwise-achievable score. +} +\examples{ +\dontrun{ +qc <- calculateQCScore(df_feat) +weights <- calculatePeptideWeights(ratios) +fit <- doseResponseFit(ratios, increasing = TRUE, precalculated_ratios = TRUE) +conf <- calculateConfidence(weights, fit, qc, df_feat) + +conf \%>\% arrange(desc(confidence)) \%>\% head() +} + +} diff --git a/man/calculatePeptideWeights.Rd b/man/calculatePeptideWeights.Rd index f0c8f67..ea3fbba 100644 --- a/man/calculatePeptideWeights.Rd +++ b/man/calculatePeptideWeights.Rd @@ -35,10 +35,13 @@ calculatePeptideWeights( \value{ Input data frame with added columns: \itemize{ -\item n_obs: Number of observations for this peptide -\item coverage_score: Proportion of timepoints observed -\item light_intensity_score: Normalized median light intensity (per protein) -\item monotonicity_score: Kendall correlation (time vs response), 0 if decreasing +\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) } @@ -47,6 +50,12 @@ Input data frame with added columns: 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 @@ -55,11 +64,11 @@ ratios <- calculateTurnoverRatios(feature_data) # Add quality weights ratios_weighted <- calculatePeptideWeights(ratios) -# Inspect weights +# Inspect coverage diagnostics ratios_weighted \%>\% group_by(Protein, BaseSequence) \%>\% slice(1) \%>\% - select(Protein, BaseSequence, coverage_score, monotonicity_score, weight) + select(Protein, BaseSequence, k_obs, p_protein, coverage_score, monotonicity_score, weight) # Use with doseResponseFit result <- doseResponseFit( diff --git a/man/calculateQCScore.Rd b/man/calculateQCScore.Rd new file mode 100644 index 0000000..50f2ca6 --- /dev/null +++ b/man/calculateQCScore.Rd @@ -0,0 +1,61 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/protein_turnover_ratio_helper.R +\name{calculateQCScore} +\alias{calculateQCScore} +\title{Calculate per-protein QC score from light-channel coverage} +\usage{ +calculateQCScore( + feature_data, + protein_col = "PROTEIN", + peptide_col = "PEPTIDE", + label_col = "LABEL", + intensity_col = "INTENSITY", + time_col = "GROUP", + light_label = "L" +) +} +\arguments{ +\item{feature_data}{Data frame from MSstats dataProcess()$FeatureLevelData.} + +\item{protein_col}{Character. Column containing protein identifiers. Default = "PROTEIN"} + +\item{peptide_col}{Character. Column containing peptide sequences. Default = "PEPTIDE"} + +\item{label_col}{Character. Column containing Heavy/Light labels. Default = "LABEL"} + +\item{intensity_col}{Character. Column with intensity values. Default = "INTENSITY"} + +\item{time_col}{Character. Column containing timepoint information. Default = "GROUP"} + +\item{light_label}{Character. Value in label_col indicating light channel. Default = "L"} +} +\value{ +Data frame with one row per protein containing: +\itemize{ +\item protein identifier (column name per \code{protein_col}) +\item n_light_peptides: distinct peptides observed in light channel +\item observed_cells: distinct (peptide x timepoint) cells observed +\item n_max_possible: n_light_peptides x n_distinct_timepoints +\item qc_score: observed_cells / n_max_possible, capped at 1 +} +} +\description{ +Computes the fraction of expected (peptide x timepoint) light-channel +observations actually detected for each protein. This is a pure +measurement-quality metric -- it ignores the heavy channel entirely, so it +remains meaningful even for proteins where no heavy incorporation occurred +(e.g., very long-lived proteins, or peptides without heavy-label residues). +} +\details{ +For each protein with n distinct light peptides observed across T total +timepoints in the experiment, the maximum possible (peptide x timepoint) +light measurements is n * T. The qc_score is the fraction of those that +were actually observed, capped at 1. +} +\examples{ +\dontrun{ +qc <- calculateQCScore(quant_data$FeatureLevelData) +qc \%>\% arrange(desc(qc_score)) \%>\% head() +} + +} diff --git a/man/classifyTurnoverProteins.Rd b/man/classifyTurnoverProteins.Rd new file mode 100644 index 0000000..26fe58a --- /dev/null +++ b/man/classifyTurnoverProteins.Rd @@ -0,0 +1,77 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/protein_turnover_ratio_helper.R +\name{classifyTurnoverProteins} +\alias{classifyTurnoverProteins} +\title{Classify proteins into turnover categories and confidence tiers} +\usage{ +classifyTurnoverProteins( + weights_df, + fit_df, + qc_df, + conf_df, + high_quantile = 0.85, + low_quantile = 0.25, + min_obs = 3, + target_short = 0.21, + target_long = 0.5 +) +} +\arguments{ +\item{weights_df}{Output of \code{calculatePeptideWeights()}.} + +\item{fit_df}{Output of \code{doseResponseFit()}.} + +\item{qc_df}{Output of \code{calculateQCScore()}.} + +\item{conf_df}{Output of \code{calculateConfidence()}.} + +\item{high_quantile}{Numeric. Upper percentile cutoff for HIGH tier. Default = 0.85 (top 15\%).} + +\item{low_quantile}{Numeric. Lower percentile cutoff for LOW tier. Default = 0.25 (bottom 25\%).} + +\item{min_obs}{Numeric. Minimum observations required for HIGH tier. Default = 3.} + +\item{target_short}{Numeric. Lower response target passed to predictIC50() for lifetime +classification. Default = 0.21.} + +\item{target_long}{Numeric. Upper response target passed to predictIC50(). Default = 0.50.} +} +\value{ +Data frame with one row per protein containing input columns plus: +\itemize{ +\item max_h_frac: per-protein maximum H_frac +\item category: one of \code{fit}, \code{medium_lived}, \code{long_lived}, \code{fast}, \code{no_heavy} +\item tier: one of \code{HIGH}, \code{MEDIUM}, \code{LOW} +} +} +\description{ +Assigns each protein a biological \code{category} describing its turnover +behavior and a \code{tier} (HIGH / MEDIUM / LOW) summarizing scoring confidence. +Combines QC + confidence + IC50 predictions in one call. +} +\details{ +Categories: +\itemize{ +\item \code{fit}: IC50 reached at the long-target response (default 0.50) +\item \code{medium_lived}: reached the short target (0.21) but not the long target +\item \code{long_lived}: failed both targets and max H_frac stayed below 0.5 +\item \code{fast}: failed both targets but max H_frac exceeds 0.5 (IC50 below observed range) +\item \code{no_heavy}: no fit was possible (no paired heavy peptides) +} + +Tiers use percentile cutoffs computed from the input data. Proteins with a +fit are tiered on \code{confidence}; \code{no_heavy} proteins are tiered on \code{qc_score}. +HIGH tier additionally requires a minimum number of observations. +} +\examples{ +\dontrun{ +qc <- calculateQCScore(df_feat) +weights <- calculatePeptideWeights(ratios) +fit <- doseResponseFit(ratios, increasing = TRUE, precalculated_ratios = TRUE) +conf <- calculateConfidence(weights, fit, qc, df_feat, k_shrinkage = 2) + +final <- classifyTurnoverProteins(weights, fit, qc, conf) +final \%>\% count(category, tier) +} + +} diff --git a/man/kendall_monotonicity.Rd b/man/kendall_monotonicity.Rd new file mode 100644 index 0000000..6ec1f61 --- /dev/null +++ b/man/kendall_monotonicity.Rd @@ -0,0 +1,24 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/protein_turnover_ratio_helper.R +\name{kendall_monotonicity} +\alias{kendall_monotonicity} +\title{Kendall monotonicity score, robust to sparse or all-missing data} +\usage{ +kendall_monotonicity(time, response) +} +\arguments{ +\item{time}{Numeric vector of timepoints.} + +\item{response}{Numeric vector of responses (same length as \code{time}).} +} +\value{ +A single numeric monotonicity score in [0, 1]. +} +\description{ +Computes \verb{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 \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. +} +\keyword{internal}