From 55461be30a658baeea7c6b2c52862c8fe4e2727f Mon Sep 17 00:00:00 2001 From: Tony Wu Date: Mon, 31 Aug 2026 11:54:51 -0400 Subject: [PATCH 1/3] feat(turnover): Incorporate turnover confidence and classification --- NAMESPACE | 3 + R/MSstatsShiny.R | 1 + R/module-qc-ui.R | 6 +- R/module-statmodel-server.R | 30 +- R/statmodel-server-download-code.R | 36 +- R/statmodel-server-results-table.R | 52 +++ R/statmodel-server-turnover-confidence.R | 171 ++++++++++ R/statmodel-ui-results.R | 5 +- .../test-statmodel-turnover-confidence.R | 319 ++++++++++++++++++ 9 files changed, 617 insertions(+), 6 deletions(-) create mode 100644 R/statmodel-server-turnover-confidence.R create mode 100644 tests/testthat/test-statmodel-turnover-confidence.R diff --git a/NAMESPACE b/NAMESPACE index a41c43d1..763f8e1d 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -51,8 +51,11 @@ importFrom(MSstatsPTM,SkylinetoMSstatsPTMFormat) importFrom(MSstatsPTM,SpectronauttoMSstatsPTMFormat) importFrom(MSstatsPTM,dataProcessPlotsPTM) importFrom(MSstatsPTM,groupComparisonPlotsPTM) +importFrom(MSstatsResponse,calculateConfidence) importFrom(MSstatsResponse,calculatePeptideWeights) +importFrom(MSstatsResponse,calculateQCScore) importFrom(MSstatsResponse,calculateTurnoverRatios) +importFrom(MSstatsResponse,classifyTurnoverProteins) importFrom(MSstatsResponse,doseResponseFit) importFrom(MSstatsResponse,futureExperimentSimulation) importFrom(MSstatsResponse,plot_tpr_power_curve) diff --git a/R/MSstatsShiny.R b/R/MSstatsShiny.R index 9d0ff37a..b22b042b 100644 --- a/R/MSstatsShiny.R +++ b/R/MSstatsShiny.R @@ -39,6 +39,7 @@ #' @importFrom MSstatsPTM dataProcessPlotsPTM groupComparisonPlotsPTM MaxQtoMSstatsPTMFormat PDtoMSstatsPTMFormat FragPipetoMSstatsPTMFormat SkylinetoMSstatsPTMFormat MetamorpheusToMSstatsPTMFormat SpectronauttoMSstatsPTMFormat #' @importFrom MSstatsBioNet exportNetworkToHTML deleteEdgeFromNetwork #' @importFrom MSstatsResponse futureExperimentSimulation run_tpr_simulation plot_tpr_power_curve calculateTurnoverRatios calculatePeptideWeights +#' @importFrom MSstatsResponse calculateQCScore calculateConfidence classifyTurnoverProteins #' @importFrom utils capture.output head packageVersion write.csv #' @importFrom stats aggregate #' @importFrom methods is diff --git a/R/module-qc-ui.R b/R/module-qc-ui.R index c9700983..63c3a0e2 100644 --- a/R/module-qc-ui.R +++ b/R/module-qc-ui.R @@ -125,7 +125,11 @@ qcUI <- function(id) { icon("question-circle", lib = "font-awesome"), div("Compute per-peptide quality weights (coverage, \ intensity, monotonicity, validity) and add them as \ - extra columns to the Turnover Ratios table.", + extra columns to the Turnover Ratios table. The \ + weights also down-weight low-quality peptides in \ + the curve fit and unlock the per-protein \ + confidence score and turnover classification on \ + the statistical modeling page.", class = "icon-tooltip")), value = FALSE ) diff --git a/R/module-statmodel-server.R b/R/module-statmodel-server.R index 10575469..87e380a9 100644 --- a/R/module-statmodel-server.R +++ b/R/module-statmodel-server.R @@ -377,7 +377,34 @@ statmodelServer = function(id, parent_session, loadpage_input, qc_input, ratio_response = FALSE, precalculated_ratios = TRUE ) - list(ComparisonResult = response_results) + + # Per-protein confidence + category/tier, only when the user asked + # for per-peptide weights on the data-processing page. + classification <- NULL + if (turnover_confidence_applies(dia_prepared, increasing)) { + show_modal_spinner(text = "Scoring per-protein confidence...") + classification <- tryCatch( + classify_turnover_fit(dia_prepared, response_results, + preprocess_data()$FeatureLevelData), + error = function(e) { + showNotification( + paste0("Turnover curves were fitted, but per-protein ", + "confidence scores were not calculated: ", + conditionMessage(e)), + type = "warning", duration = 10) + NULL + }, + finally = remove_modal_spinner() + ) + response_results <- merge_turnover_confidence(response_results, + classification) + } else if (turnover_weights_present(dia_prepared)) { + showNotification(turnover_confidence_direction_message(), + type = "message", duration = 10) + } + + list(ComparisonResult = response_results, + TurnoverClassification = classification) } else if (app_template() == TEMPLATES$chemoproteomics) { meta <- condition_metadata() req(!is.null(meta) && "DoseVal" %in% colnames(meta)) @@ -467,6 +494,7 @@ statmodelServer = function(id, parent_session, loadpage_input, qc_input, # Results rendering render_results_table(output, session, data_comparison, SignificantProteins, app_template = app_template) + render_turnover_confidence_table(output, session, data_comparison, app_template = app_template) render_ptm_results_tables(output, session, data_comparison, SignificantProteins) # Download handlers diff --git a/R/statmodel-server-download-code.R b/R/statmodel-server-download-code.R index eaf2f60c..64be7ea9 100644 --- a/R/statmodel-server-download-code.R +++ b/R/statmodel-server-download-code.R @@ -216,8 +216,10 @@ format_r_double <- function(value) { #' Mirrors the app's turnover flow: calculateTurnoverRatios (with the resolved #' tracer constants), an optional calculatePeptideWeights step when #' "Assign feature weights" is checked, then a weighted doseResponseFit and -#' visualizeResponseProtein. Kept as a plain string builder so it can be -#' unit-tested without a running session. +#' visualizeResponseProtein. When weights were requested and the fit is in the +#' synthesis direction, the calculateQCScore / calculateConfidence / +#' classifyTurnoverProteins scoring chain is emitted too. Kept as a plain +#' string builder so it can be unit-tested without a running session. #' #' @param qc_input The QC module input list (the assign_feature_weights checkbox #' lives here). @@ -347,8 +349,36 @@ build_turnover_analysis_code <- function(qc_input, comp_mat, increasing, " increasing = ", increasing, ",\n", " transform_dose = FALSE,\n", " ratio_response = FALSE,\n", - " precalculated_ratios = TRUE\n)\n", + " precalculated_ratios = TRUE\n)\n" + ) + + if (weighting && isTRUE(increasing)) { + code <- paste0( + code, + "\n# Score per-protein confidence and classify turnover behavior. Needs the\n", + "# per-peptide weights above, and is defined for the synthesis direction only.\n", + "# calculateConfidence() and classifyTurnoverProteins() read 'Protein' and\n", + "# 'H_frac'; predictIC50() (called inside the classifier) reads the lower-case\n", + "# fit columns, so the same frame carries both namings.\n", + "classification_input = prepared_data\n", + "classification_input$Protein = as.character(classification_input$protein)\n", + "classification_input$H_frac = classification_input$response\n", + "qc_scores = calculateQCScore(summarized$FeatureLevelData)\n", + "confidence_scores = calculateConfidence(\n", + " weights_df = classification_input,\n", + " fit_df = response_results,\n", + " qc_df = qc_scores,\n", + " feature_data = summarized$FeatureLevelData)\n", + "turnover_classification = classifyTurnoverProteins(\n", + " weights_df = classification_input,\n", + " fit_df = response_results,\n", + " qc_df = qc_scores,\n", + " conf_df = confidence_scores)\n" + ) + } + code <- paste0( + code, "\n# Visualize a single protein's turnover curve\n", "visualizeResponseProtein(\n", " data = prepared_data,\n", diff --git a/R/statmodel-server-results-table.R b/R/statmodel-server-results-table.R index 84876960..4790594c 100644 --- a/R/statmodel-server-results-table.R +++ b/R/statmodel-server-results-table.R @@ -26,6 +26,58 @@ render_results_table = function(output, session, data_comparison, SignificantPro output$number = renderText({ nrow(SignificantProteins()) }) } +#' Register the turnover confidence / classification results panel. +#' +#' The panel only appears for the protein-turnover template when the analysis +#' produced a classification, i.e. when per-peptide weights were calculated on +#' the data-processing page. Unlike the main results table this one has a row +#' per protein in the data, including proteins that produced no fit. +#' @noRd +render_turnover_confidence_table = function(output, session, data_comparison, + app_template = reactive(TEMPLATES$default)) { + ns = session$ns + + classification = reactive({ + if (!isTRUE(app_template() == TEMPLATES$protein_turnover)) { + return(NULL) + } + data_comparison()$TurnoverClassification + }) + + output$turnover_confidence_results = renderUI({ + classified = classification() + if (is.null(classified) || NROW(classified) == 0) { + return(NULL) + } + tagList( + tags$br(), + h2("Turnover confidence and classification"), + h5("One row per protein. ", tags$code("confidence"), " combines the ", + "per-peptide weights, the fit residuals, the light-channel QC score ", + "and a shrinkage factor on heavy-peptide count. ", + tags$code("category"), " describes turnover behavior (fit, ", + "medium_lived, long_lived, fast, no_heavy) and ", tags$code("tier"), + " ranks scoring confidence (HIGH / MEDIUM / LOW). Proteins with no ", + "fit have NA scores."), + tags$br(), + dataTableOutput(ns("turnover_confidence_table")), + downloadButton(ns("download_turnover_confidence"), + "Download confidence scores") + ) + }) + + output$turnover_confidence_table = renderDataTable({ + req(classification()) + }, options = list(scrollX = TRUE)) + + output$download_turnover_confidence = downloadHandler( + filename = function() paste0("Turnover_Confidence-", Sys.Date(), ".csv"), + content = function(file) { + write.csv(classification(), file, row.names = FALSE) + } + ) +} + render_ptm_results_tables = function(output, session, data_comparison, SignificantProteins) { ns = session$ns diff --git a/R/statmodel-server-turnover-confidence.R b/R/statmodel-server-turnover-confidence.R new file mode 100644 index 00000000..f66d678a --- /dev/null +++ b/R/statmodel-server-turnover-confidence.R @@ -0,0 +1,171 @@ +# Per-protein confidence scoring and turnover classification +# (MSstatsResponse::calculateQCScore -> calculateConfidence -> +# classifyTurnoverProteins). +# +# This step only runs when the user asked for per-peptide weights on the +# data-processing page: calculateConfidence averages the `weight` column, so +# without it there is nothing to score. Kept as plain functions so the whole +# chain can be unit-tested without a running session. + +#' Per-protein score columns lifted from the classification table onto the +#' turnover fit result, in display order. +#' @noRd +TURNOVER_CONFIDENCE_COLUMNS <- c("mean_weight", "n_obs", "qc_score", + "n_heavy_peptides", "confidence", + "max_h_frac", "category", "tier") + +#' Feature-level columns calculateQCScore / calculateConfidence read. +#' @noRd +TURNOVER_CONFIDENCE_FEATURE_COLUMNS <- c("PROTEIN", "PEPTIDE", "LABEL", + "INTENSITY", "GROUP") + +#' Are per-peptide weights available on a prepared dose-response frame? +#' +#' @param prepared Output of `prepare_turnover_for_dose_response()`. +#' @return TRUE when the frame has rows and carries a `weight` column. +#' @noRd +turnover_weights_present <- function(prepared) { + !is.null(prepared) && NROW(prepared) > 0 && + "weight" %in% colnames(prepared) +} + +#' Does the confidence / classification step apply to this fit? +#' +#' Two conditions. Weights must be present (they are the input +#' calculateConfidence averages), and the fit must be in the synthesis +#' direction: classifyTurnoverProteins scores `H_frac` against increasing +#' IC50 targets, so a degradation (`L_frac`) fit cannot be classified. +#' +#' @param prepared Output of `prepare_turnover_for_dose_response()`. +#' @param increasing Logical. The fit's trend direction. +#' @return TRUE when `classify_turnover_fit()` can be called. +#' @noRd +turnover_confidence_applies <- function(prepared, increasing) { + turnover_weights_present(prepared) && isTRUE(increasing) +} + +#' Reshape a prepared dose-response frame into classifyTurnoverProteins input. +#' +#' MSstatsResponse wants both naming conventions in one frame: predictIC50() +#' reads the lower-case fit columns (`protein`, `drug`, `dose`, `response`), +#' while calculateConfidence() and classifyTurnoverProteins() read `Protein` +#' and `H_frac`. Deriving the latter from the former (rather than re-deriving +#' them from the raw ratios) guarantees the scores describe exactly the rows +#' that were fitted. +#' +#' @param prepared Output of `prepare_turnover_for_dose_response()`, which must +#' carry a `weight` column. +#' @return `prepared` with `Protein` and `H_frac` columns added. +#' @noRd +prepare_turnover_for_classification <- function(prepared) { + prepared <- as.data.frame(prepared, stringsAsFactors = FALSE) + prepared$Protein <- as.character(prepared$protein) + prepared$H_frac <- prepared$response + prepared +} + +#' Normalize feature-level data for the QC-score / heavy-peptide counts. +#' +#' dataProcess() returns PROTEIN and PEPTIDE as factors; coercing them keeps +#' the joined `Protein` column a character vector so the score columns match +#' the fit result by value rather than by factor level. +#' +#' @param feature_data `preprocess_data()$FeatureLevelData`. +#' @return A data frame with character protein / peptide identifiers. +#' @noRd +prepare_feature_data_for_qc_score <- function(feature_data) { + feature_data <- as.data.frame(feature_data, stringsAsFactors = FALSE) + missing <- setdiff(TURNOVER_CONFIDENCE_FEATURE_COLUMNS, colnames(feature_data)) + if (length(missing) > 0) { + stop("the feature-level data is missing required column(s): ", + paste(missing, collapse = ", "), + ". Re-run protein summarization on the data-processing page.", + call. = FALSE) + } + for (col in c("PROTEIN", "PEPTIDE", "LABEL")) { + feature_data[[col]] <- as.character(feature_data[[col]]) + } + feature_data +} + +#' Score and classify a turnover fit. +#' +#' Runs the MSstatsResponse chain: light-channel QC score per protein, then +#' the combined confidence score (peptide weights x fit residuals x QC score x +#' heavy-peptide shrinkage), then the biological category / confidence tier. +#' +#' The result has one row per protein in the feature-level data, which is a +#' superset of the fitted proteins: proteins that produced no fit come back +#' with NA scores and category `no_heavy`. +#' +#' @param prepared Output of `prepare_turnover_for_dose_response()` with weights. +#' @param fit Output of `doseResponseFit()`. +#' @param feature_data `preprocess_data()$FeatureLevelData`. +#' @param k_shrinkage Numeric. Bayesian shrinkage constant on the heavy-peptide +#' count, passed to `calculateConfidence()`. +#' @return A data frame of per-protein QC, confidence, category and tier. +#' @noRd +classify_turnover_fit <- function(prepared, fit, feature_data, + k_shrinkage = 2) { + weights_df <- prepare_turnover_for_classification(prepared) + features <- prepare_feature_data_for_qc_score(feature_data) + + qc_scores <- calculateQCScore(features) + confidence_scores <- calculateConfidence( + weights_df = weights_df, + fit_df = fit, + qc_df = qc_scores, + feature_data = features, + k_shrinkage = k_shrinkage + ) + classification <- classifyTurnoverProteins( + weights_df = weights_df, + fit_df = fit, + qc_df = qc_scores, + conf_df = confidence_scores + ) + + as.data.frame(classification, stringsAsFactors = FALSE) +} + +#' Attach the per-protein score columns to the turnover fit result. +#' +#' Row order and row count of `fit` are preserved, so the significance +#' filtering and downloads that already read `ComparisonResult` are unaffected; +#' the classification's extra (unfitted) proteins stay in the classification +#' table only. +#' +#' @param fit Output of `doseResponseFit()`. +#' @param classification Output of `classify_turnover_fit()`. +#' @return `fit` with the confidence / category / tier columns appended. +#' @noRd +merge_turnover_confidence <- function(fit, classification) { + if (is.null(classification) || NROW(classification) == 0 || + is.null(fit) || NROW(fit) == 0) { + return(fit) + } + score_cols <- setdiff( + intersect(TURNOVER_CONFIDENCE_COLUMNS, colnames(classification)), + colnames(fit)) + if (length(score_cols) == 0) { + return(fit) + } + + matched <- match(as.character(fit$Protein), + as.character(classification$Protein)) + for (col in score_cols) { + fit[[col]] <- classification[[col]][matched] + } + fit +} + +#' The notification shown when weights were calculated but the fit is a +#' degradation fit, so no classification is possible. +#' @noRd +turnover_confidence_direction_message <- function() { + paste0("Per-protein confidence scores and turnover categories were not ", + "calculated: they are defined for the synthesis direction (heavy ", + "fraction, increasing over time) only. Check \"Synthesis ", + "(heavy-isotope incorporation, increasing)\" and calculate again to ", + "score this fit.") +} diff --git a/R/statmodel-ui-results.R b/R/statmodel-ui-results.R index 6ab1e8c6..32d19e51 100644 --- a/R/statmodel-ui-results.R +++ b/R/statmodel-ui-results.R @@ -36,7 +36,10 @@ create_results_tables <- function(ns) { ), conditionalPanel( condition = "input['loadpage-BIO']!=='PTM'", - uiOutput(ns("table_results")) + uiOutput(ns("table_results")), + # Protein-turnover only, and only when per-peptide weights were + # calculated; empty otherwise. + uiOutput(ns("turnover_confidence_results")) ) ) } \ No newline at end of file diff --git a/tests/testthat/test-statmodel-turnover-confidence.R b/tests/testthat/test-statmodel-turnover-confidence.R new file mode 100644 index 00000000..97c33361 --- /dev/null +++ b/tests/testthat/test-statmodel-turnover-confidence.R @@ -0,0 +1,319 @@ +context("Turnover confidence scoring and classification") + +# ============================================================================ +# Fixtures +# +# `prepared` mirrors prepare_turnover_for_dose_response() output; `features` +# mirrors dataProcess()$FeatureLevelData (factor PROTEIN / PEPTIDE, as +# dataProcess returns them); `fit` mirrors doseResponseFit() output. +# ============================================================================ + +make_prepared <- function(with_weight = TRUE) { + prepared <- data.frame( + protein = rep(c("ProtA", "ProtB"), each = 6), + drug = "time", + dose = rep(c(0, 2, 6), times = 4), + response = c(0.0, 0.3, 0.6, 0.0, 0.35, 0.65, + 0.0, 0.1, 0.2, 0.0, 0.15, 0.25), + BaseSequence = rep(c("PEP1", "PEP2", "PEP1", "PEP2"), each = 3), + stringsAsFactors = FALSE + ) + if (with_weight) { + prepared$weight <- rep(c(0.8, 0.4), each = 6) + } + prepared +} + +make_features <- function() { + features <- expand.grid( + PROTEIN = c("ProtA", "ProtB", "ProtC"), + PEPTIDE = c("PEP1", "PEP2"), + GROUP = c(0, 2, 6), + LABEL = c("H", "L"), + stringsAsFactors = FALSE + ) + # ProtC never incorporated label: light channel only. + features <- features[!(features$PROTEIN == "ProtC" & features$LABEL == "H"), ] + features$RUN <- paste0("run_", features$GROUP) + features$INTENSITY <- seq_len(nrow(features)) * 1000 + features$PROTEIN <- factor(features$PROTEIN) + features$PEPTIDE <- factor(features$PEPTIDE) + features +} + +make_fit <- function() { + data.frame( + Protein = c("ProtA", "ProtB"), + drug = "time", + direction = "increasing", + SSE_Full = c(0.01, 0.05), + SSE_Null = c(0.5, 0.4), + F_statistic = c(20, 8), + P_value = c(0.001, 0.02), + log2FC = c(1.5, 0.8), + adj.pvalue = c(0.002, 0.02), + stringsAsFactors = FALSE + ) +} + +# ============================================================================ +# turnover_weights_present / turnover_confidence_applies +# ============================================================================ + +test_that("turnover_weights_present detects the weight column", { + expect_true(MSstatsShiny:::turnover_weights_present(make_prepared())) + expect_false(MSstatsShiny:::turnover_weights_present( + make_prepared(with_weight = FALSE))) +}) + +test_that("turnover_weights_present is FALSE for NULL or empty input", { + expect_false(MSstatsShiny:::turnover_weights_present(NULL)) + expect_false(MSstatsShiny:::turnover_weights_present( + make_prepared()[0, , drop = FALSE])) +}) + +test_that("turnover_confidence_applies requires weights and the synthesis direction", { + prepared <- make_prepared() + + expect_true(MSstatsShiny:::turnover_confidence_applies(prepared, TRUE)) + expect_false(MSstatsShiny:::turnover_confidence_applies(prepared, FALSE), + info = "classification is only defined for H_frac / increasing fits") + expect_false(MSstatsShiny:::turnover_confidence_applies( + make_prepared(with_weight = FALSE), TRUE), + info = "confidence averages the weight column, so weights are required") +}) + +# ============================================================================ +# prepare_turnover_for_classification +# ============================================================================ + +test_that("prepare_turnover_for_classification adds the columns MSstatsResponse reads", { + result <- MSstatsShiny:::prepare_turnover_for_classification(make_prepared()) + + expect_true(all(c("Protein", "H_frac") %in% colnames(result)), + info = "calculateConfidence / classifyTurnoverProteins read Protein and H_frac") + expect_true(all(c("protein", "drug", "dose", "response") %in% colnames(result)), + info = "predictIC50 still needs the lower-case fit columns") + expect_type(result$Protein, "character") + expect_equal(result$Protein, result$protein) + expect_equal(result$H_frac, result$response, + info = "H_frac must be the exact column that was fitted") +}) + +test_that("prepare_turnover_for_classification coerces a factor protein column", { + prepared <- make_prepared() + prepared$protein <- factor(prepared$protein) + + result <- MSstatsShiny:::prepare_turnover_for_classification(prepared) + + expect_type(result$Protein, "character") + expect_equal(unique(result$Protein), c("ProtA", "ProtB")) +}) + +# ============================================================================ +# prepare_feature_data_for_qc_score +# ============================================================================ + +test_that("prepare_feature_data_for_qc_score makes identifiers character", { + result <- MSstatsShiny:::prepare_feature_data_for_qc_score(make_features()) + + expect_type(result$PROTEIN, "character") + expect_type(result$PEPTIDE, "character") + expect_type(result$LABEL, "character") +}) + +test_that("prepare_feature_data_for_qc_score names the missing columns", { + features <- make_features() + features$LABEL <- NULL + features$INTENSITY <- NULL + + expect_error( + MSstatsShiny:::prepare_feature_data_for_qc_score(features), + "LABEL, INTENSITY") +}) + +# ============================================================================ +# classify_turnover_fit +# ============================================================================ + +test_that("classify_turnover_fit returns a confidence and tier per protein", { + result <- MSstatsShiny:::classify_turnover_fit( + make_prepared(), make_fit(), make_features()) + + expect_true(all(c("Protein", "qc_score", "mean_weight", "n_heavy_peptides", + "confidence", "max_h_frac", "category", "tier") %in% + colnames(result))) + expect_setequal(result$Protein, c("ProtA", "ProtB", "ProtC")) + expect_true(all(result$confidence >= 0 & result$confidence <= 1, + na.rm = TRUE), + info = "confidence is a score in [0, 1]") + expect_true(all(result$tier %in% c("HIGH", "MEDIUM", "LOW"))) +}) + +test_that("classify_turnover_fit keeps proteins that produced no fit", { + result <- MSstatsShiny:::classify_turnover_fit( + make_prepared(), make_fit(), make_features()) + + prot_c <- result[result$Protein == "ProtC", ] + expect_equal(nrow(prot_c), 1L) + expect_true(is.na(prot_c$confidence), + info = "a protein with no fit cannot be scored") + expect_equal(prot_c$category, "no_heavy") +}) + +test_that("classify_turnover_fit weights higher-quality peptides into confidence", { + prepared <- make_prepared() + low <- prepared + low$weight <- low$weight / 4 + + high_conf <- MSstatsShiny:::classify_turnover_fit( + prepared, make_fit(), make_features()) + low_conf <- MSstatsShiny:::classify_turnover_fit( + low, make_fit(), make_features()) + + fitted <- c("ProtA", "ProtB") + expect_true(all( + low_conf$confidence[match(fitted, low_conf$Protein)] < + high_conf$confidence[match(fitted, high_conf$Protein)]), + info = "down-weighting every peptide must lower every protein's confidence") +}) + +test_that("classify_turnover_fit surfaces a missing feature-level column", { + features <- make_features() + features$PEPTIDE <- NULL + + expect_error( + MSstatsShiny:::classify_turnover_fit(make_prepared(), make_fit(), features), + "PEPTIDE") +}) + +# ============================================================================ +# merge_turnover_confidence +# ============================================================================ + +test_that("merge_turnover_confidence appends scores without changing the fit rows", { + fit <- make_fit() + classification <- MSstatsShiny:::classify_turnover_fit( + make_prepared(), fit, make_features()) + + result <- MSstatsShiny:::merge_turnover_confidence(fit, classification) + + expect_equal(nrow(result), nrow(fit), + info = "the unfitted proteins stay out of the fit result") + expect_equal(result$Protein, fit$Protein, + info = "row order must survive so downstream filtering is unaffected") + expect_equal(result$adj.pvalue, fit$adj.pvalue) + expect_true(all(c("confidence", "category", "tier") %in% colnames(result))) + expect_equal(result$confidence, + classification$confidence[match(fit$Protein, + classification$Protein)]) +}) + +test_that("merge_turnover_confidence matches by protein, not by row position", { + fit <- make_fit() + classification <- data.frame( + Protein = c("ProtB", "ProtA"), + confidence = c(0.2, 0.9), + category = c("medium_lived", "fit"), + tier = c("LOW", "HIGH"), + stringsAsFactors = FALSE + ) + + result <- MSstatsShiny:::merge_turnover_confidence(fit, classification) + + expect_equal(result$confidence[result$Protein == "ProtA"], 0.9) + expect_equal(result$tier[result$Protein == "ProtB"], "LOW") +}) + +test_that("merge_turnover_confidence returns the fit unchanged when there is no classification", { + fit <- make_fit() + + expect_equal(MSstatsShiny:::merge_turnover_confidence(fit, NULL), fit) + expect_equal( + MSstatsShiny:::merge_turnover_confidence(fit, data.frame()), fit) +}) + +test_that("merge_turnover_confidence does not overwrite existing fit columns", { + fit <- make_fit() + fit$confidence <- c(99, 99) + classification <- data.frame( + Protein = c("ProtA", "ProtB"), + confidence = c(0.1, 0.2), + tier = c("HIGH", "LOW"), + stringsAsFactors = FALSE + ) + + result <- MSstatsShiny:::merge_turnover_confidence(fit, classification) + + expect_equal(result$confidence, c(99, 99)) + expect_equal(result$tier, c("HIGH", "LOW")) +}) + +# ============================================================================ +# Downloadable analysis code +# ============================================================================ + +test_that("weighted synthesis script emits the confidence and classification chain", { + code <- MSstatsShiny:::build_turnover_analysis_code( + qc_input = stats::setNames(list(TRUE), NAMESPACE_QC$assign_feature_weights), + comp_mat = data.frame(GROUP = c("0hr", "2hr"), TimeVal = c(0, 2), + stringsAsFactors = FALSE), + increasing = TRUE, + tracer_constants = list(values = c("0hr" = 1, "2hr" = 1), + source = CONSTANTS_QC$tracer_source_none, + file = NULL) + ) + + expect_true(grepl("calculateQCScore(summarized$FeatureLevelData)", code, + fixed = TRUE)) + expect_true(grepl("confidence_scores = calculateConfidence(", code, + fixed = TRUE)) + expect_true(grepl("turnover_classification = classifyTurnoverProteins(", code, + fixed = TRUE)) + expect_true(grepl("classification_input$H_frac", code, fixed = TRUE), + info = "the classifier needs the H_frac alias on the fitted frame") + expect_lt(regexpr("calculateConfidence", code, fixed = TRUE), + regexpr("classifyTurnoverProteins", code, fixed = TRUE)) + expect_lt(regexpr("response_results = doseResponseFit", code, fixed = TRUE), + regexpr("calculateConfidence", code, fixed = TRUE)) +}) + +test_that("unweighted script omits the confidence and classification chain", { + code <- MSstatsShiny:::build_turnover_analysis_code( + qc_input = stats::setNames(list(FALSE), NAMESPACE_QC$assign_feature_weights), + comp_mat = data.frame(GROUP = c("0hr", "2hr"), TimeVal = c(0, 2), + stringsAsFactors = FALSE), + increasing = TRUE, + tracer_constants = list(values = c("0hr" = 1, "2hr" = 1), + source = CONSTANTS_QC$tracer_source_none, + file = NULL) + ) + + expect_false(grepl("calculateConfidence", code, fixed = TRUE)) + expect_false(grepl("classifyTurnoverProteins", code, fixed = TRUE)) + expect_false(grepl("calculateQCScore", code, fixed = TRUE)) +}) + +test_that("degradation script omits the classification chain even when weighted", { + code <- MSstatsShiny:::build_turnover_analysis_code( + qc_input = stats::setNames(list(TRUE), NAMESPACE_QC$assign_feature_weights), + comp_mat = data.frame(GROUP = c("0hr", "2hr"), TimeVal = c(0, 2), + stringsAsFactors = FALSE), + increasing = FALSE, + tracer_constants = list(values = c("0hr" = 1, "2hr" = 1), + source = CONSTANTS_QC$tracer_source_none, + file = NULL) + ) + + expect_true(grepl("calculatePeptideWeights", code, fixed = TRUE), + info = "weights still apply to a degradation fit") + expect_false(grepl("classifyTurnoverProteins", code, fixed = TRUE), + info = "classification is H_frac / increasing only") +}) + +test_that("turnover_confidence_direction_message names the checkbox to re-check", { + message_text <- MSstatsShiny:::turnover_confidence_direction_message() + + expect_true(grepl("Synthesis", message_text, fixed = TRUE)) + expect_true(grepl("were not calculated", message_text, fixed = TRUE)) +}) From 94898c1a7b28a18df663be01172fb4751b4f7d1d Mon Sep 17 00:00:00 2001 From: Tony Wu Date: Mon, 31 Aug 2026 12:06:39 -0400 Subject: [PATCH 2/3] remove unnecessary comments --- R/module-qc-ui.R | 7 ++-- R/module-statmodel-server.R | 2 -- R/statmodel-server-download-code.R | 12 +++---- R/statmodel-ui-results.R | 2 -- .../test-statmodel-turnover-confidence.R | 33 ------------------- 5 files changed, 7 insertions(+), 49 deletions(-) diff --git a/R/module-qc-ui.R b/R/module-qc-ui.R index 63c3a0e2..20df8b75 100644 --- a/R/module-qc-ui.R +++ b/R/module-qc-ui.R @@ -126,10 +126,9 @@ qcUI <- function(id) { div("Compute per-peptide quality weights (coverage, \ intensity, monotonicity, validity) and add them as \ extra columns to the Turnover Ratios table. The \ - weights also down-weight low-quality peptides in \ - the curve fit and unlock the per-protein \ - confidence score and turnover classification on \ - the statistical modeling page.", + weights down-weight low-quality peptides in \ + the curve fit and enable per-protein \ + confidence score and turnover classification.", class = "icon-tooltip")), value = FALSE ) diff --git a/R/module-statmodel-server.R b/R/module-statmodel-server.R index 87e380a9..f327a402 100644 --- a/R/module-statmodel-server.R +++ b/R/module-statmodel-server.R @@ -378,8 +378,6 @@ statmodelServer = function(id, parent_session, loadpage_input, qc_input, precalculated_ratios = TRUE ) - # Per-protein confidence + category/tier, only when the user asked - # for per-peptide weights on the data-processing page. classification <- NULL if (turnover_confidence_applies(dia_prepared, increasing)) { show_modal_spinner(text = "Scoring per-protein confidence...") diff --git a/R/statmodel-server-download-code.R b/R/statmodel-server-download-code.R index 64be7ea9..f5d2e20b 100644 --- a/R/statmodel-server-download-code.R +++ b/R/statmodel-server-download-code.R @@ -216,9 +216,9 @@ format_r_double <- function(value) { #' Mirrors the app's turnover flow: calculateTurnoverRatios (with the resolved #' tracer constants), an optional calculatePeptideWeights step when #' "Assign feature weights" is checked, then a weighted doseResponseFit and -#' visualizeResponseProtein. When weights were requested and the fit is in the -#' synthesis direction, the calculateQCScore / calculateConfidence / -#' classifyTurnoverProteins scoring chain is emitted too. Kept as a plain +#' visualizeResponseProtein. When weights are requested and the fit is in the +#' synthesis direction, calculateQCScore / calculateConfidence / +#' classifyTurnoverProteins is applied. Kept as a plain #' string builder so it can be unit-tested without a running session. #' #' @param qc_input The QC module input list (the assign_feature_weights checkbox @@ -355,11 +355,7 @@ build_turnover_analysis_code <- function(qc_input, comp_mat, increasing, if (weighting && isTRUE(increasing)) { code <- paste0( code, - "\n# Score per-protein confidence and classify turnover behavior. Needs the\n", - "# per-peptide weights above, and is defined for the synthesis direction only.\n", - "# calculateConfidence() and classifyTurnoverProteins() read 'Protein' and\n", - "# 'H_frac'; predictIC50() (called inside the classifier) reads the lower-case\n", - "# fit columns, so the same frame carries both namings.\n", + "\n", "classification_input = prepared_data\n", "classification_input$Protein = as.character(classification_input$protein)\n", "classification_input$H_frac = classification_input$response\n", diff --git a/R/statmodel-ui-results.R b/R/statmodel-ui-results.R index 32d19e51..aa3993eb 100644 --- a/R/statmodel-ui-results.R +++ b/R/statmodel-ui-results.R @@ -37,8 +37,6 @@ create_results_tables <- function(ns) { conditionalPanel( condition = "input['loadpage-BIO']!=='PTM'", uiOutput(ns("table_results")), - # Protein-turnover only, and only when per-peptide weights were - # calculated; empty otherwise. uiOutput(ns("turnover_confidence_results")) ) ) diff --git a/tests/testthat/test-statmodel-turnover-confidence.R b/tests/testthat/test-statmodel-turnover-confidence.R index 97c33361..1d677394 100644 --- a/tests/testthat/test-statmodel-turnover-confidence.R +++ b/tests/testthat/test-statmodel-turnover-confidence.R @@ -1,13 +1,5 @@ context("Turnover confidence scoring and classification") -# ============================================================================ -# Fixtures -# -# `prepared` mirrors prepare_turnover_for_dose_response() output; `features` -# mirrors dataProcess()$FeatureLevelData (factor PROTEIN / PEPTIDE, as -# dataProcess returns them); `fit` mirrors doseResponseFit() output. -# ============================================================================ - make_prepared <- function(with_weight = TRUE) { prepared <- data.frame( protein = rep(c("ProtA", "ProtB"), each = 6), @@ -32,7 +24,6 @@ make_features <- function() { LABEL = c("H", "L"), stringsAsFactors = FALSE ) - # ProtC never incorporated label: light channel only. features <- features[!(features$PROTEIN == "ProtC" & features$LABEL == "H"), ] features$RUN <- paste0("run_", features$GROUP) features$INTENSITY <- seq_len(nrow(features)) * 1000 @@ -56,10 +47,6 @@ make_fit <- function() { ) } -# ============================================================================ -# turnover_weights_present / turnover_confidence_applies -# ============================================================================ - test_that("turnover_weights_present detects the weight column", { expect_true(MSstatsShiny:::turnover_weights_present(make_prepared())) expect_false(MSstatsShiny:::turnover_weights_present( @@ -83,10 +70,6 @@ test_that("turnover_confidence_applies requires weights and the synthesis direct info = "confidence averages the weight column, so weights are required") }) -# ============================================================================ -# prepare_turnover_for_classification -# ============================================================================ - test_that("prepare_turnover_for_classification adds the columns MSstatsResponse reads", { result <- MSstatsShiny:::prepare_turnover_for_classification(make_prepared()) @@ -110,10 +93,6 @@ test_that("prepare_turnover_for_classification coerces a factor protein column", expect_equal(unique(result$Protein), c("ProtA", "ProtB")) }) -# ============================================================================ -# prepare_feature_data_for_qc_score -# ============================================================================ - test_that("prepare_feature_data_for_qc_score makes identifiers character", { result <- MSstatsShiny:::prepare_feature_data_for_qc_score(make_features()) @@ -132,10 +111,6 @@ test_that("prepare_feature_data_for_qc_score names the missing columns", { "LABEL, INTENSITY") }) -# ============================================================================ -# classify_turnover_fit -# ============================================================================ - test_that("classify_turnover_fit returns a confidence and tier per protein", { result <- MSstatsShiny:::classify_turnover_fit( make_prepared(), make_fit(), make_features()) @@ -187,10 +162,6 @@ test_that("classify_turnover_fit surfaces a missing feature-level column", { "PEPTIDE") }) -# ============================================================================ -# merge_turnover_confidence -# ============================================================================ - test_that("merge_turnover_confidence appends scores without changing the fit rows", { fit <- make_fit() classification <- MSstatsShiny:::classify_turnover_fit( @@ -249,10 +220,6 @@ test_that("merge_turnover_confidence does not overwrite existing fit columns", { expect_equal(result$tier, c("HIGH", "LOW")) }) -# ============================================================================ -# Downloadable analysis code -# ============================================================================ - test_that("weighted synthesis script emits the confidence and classification chain", { code <- MSstatsShiny:::build_turnover_analysis_code( qc_input = stats::setNames(list(TRUE), NAMESPACE_QC$assign_feature_weights), From 84495e5e51bac38d201b1169415d2c9f63d0bac9 Mon Sep 17 00:00:00 2001 From: Tony Wu Date: Mon, 31 Aug 2026 12:28:32 -0400 Subject: [PATCH 3/3] adjust turnover classification code --- R/statmodel-server-turnover-confidence.R | 53 +++++------------------- 1 file changed, 11 insertions(+), 42 deletions(-) diff --git a/R/statmodel-server-turnover-confidence.R b/R/statmodel-server-turnover-confidence.R index f66d678a..0a883eb3 100644 --- a/R/statmodel-server-turnover-confidence.R +++ b/R/statmodel-server-turnover-confidence.R @@ -1,12 +1,3 @@ -# Per-protein confidence scoring and turnover classification -# (MSstatsResponse::calculateQCScore -> calculateConfidence -> -# classifyTurnoverProteins). -# -# This step only runs when the user asked for per-peptide weights on the -# data-processing page: calculateConfidence averages the `weight` column, so -# without it there is nothing to score. Kept as plain functions so the whole -# chain can be unit-tested without a running session. - #' Per-protein score columns lifted from the classification table onto the #' turnover fit result, in display order. #' @noRd @@ -14,12 +5,12 @@ TURNOVER_CONFIDENCE_COLUMNS <- c("mean_weight", "n_obs", "qc_score", "n_heavy_peptides", "confidence", "max_h_frac", "category", "tier") -#' Feature-level columns calculateQCScore / calculateConfidence read. +#' Feature-level columns calculateQCScore and calculateConfidence read. #' @noRd TURNOVER_CONFIDENCE_FEATURE_COLUMNS <- c("PROTEIN", "PEPTIDE", "LABEL", "INTENSITY", "GROUP") -#' Are per-peptide weights available on a prepared dose-response frame? +#' Checks if per-peptide weights are available in a data frame #' #' @param prepared Output of `prepare_turnover_for_dose_response()`. #' @return TRUE when the frame has rows and carries a `weight` column. @@ -29,12 +20,10 @@ turnover_weights_present <- function(prepared) { "weight" %in% colnames(prepared) } -#' Does the confidence / classification step apply to this fit? +#' Checks if confidence scoring and classification can be applied #' -#' Two conditions. Weights must be present (they are the input -#' calculateConfidence averages), and the fit must be in the synthesis -#' direction: classifyTurnoverProteins scores `H_frac` against increasing -#' IC50 targets, so a degradation (`L_frac`) fit cannot be classified. +#' It can be applied if weights are present and the fit is w.r.t. the synthesis +#' direction. At the moment, classifyTurnoverProteins only handles `H_frac` #' #' @param prepared Output of `prepare_turnover_for_dose_response()`. #' @param increasing Logical. The fit's trend direction. @@ -46,13 +35,6 @@ turnover_confidence_applies <- function(prepared, increasing) { #' Reshape a prepared dose-response frame into classifyTurnoverProteins input. #' -#' MSstatsResponse wants both naming conventions in one frame: predictIC50() -#' reads the lower-case fit columns (`protein`, `drug`, `dose`, `response`), -#' while calculateConfidence() and classifyTurnoverProteins() read `Protein` -#' and `H_frac`. Deriving the latter from the former (rather than re-deriving -#' them from the raw ratios) guarantees the scores describe exactly the rows -#' that were fitted. -#' #' @param prepared Output of `prepare_turnover_for_dose_response()`, which must #' carry a `weight` column. #' @return `prepared` with `Protein` and `H_frac` columns added. @@ -64,11 +46,9 @@ prepare_turnover_for_classification <- function(prepared) { prepared } -#' Normalize feature-level data for the QC-score / heavy-peptide counts. +#' Prepare feature-level data for the QC-score / heavy-peptide counts. #' -#' dataProcess() returns PROTEIN and PEPTIDE as factors; coercing them keeps -#' the joined `Protein` column a character vector so the score columns match -#' the fit result by value rather than by factor level. +#' Specifically, turn PROTEIN, PEPTIDE, and LABEL columns into character columns #' #' @param feature_data `preprocess_data()$FeatureLevelData`. #' @return A data frame with character protein / peptide identifiers. @@ -88,15 +68,8 @@ prepare_feature_data_for_qc_score <- function(feature_data) { feature_data } -#' Score and classify a turnover fit. -#' -#' Runs the MSstatsResponse chain: light-channel QC score per protein, then -#' the combined confidence score (peptide weights x fit residuals x QC score x -#' heavy-peptide shrinkage), then the biological category / confidence tier. -#' -#' The result has one row per protein in the feature-level data, which is a -#' superset of the fitted proteins: proteins that produced no fit come back -#' with NA scores and category `no_heavy`. +#' Score and classify a turnover fit as long-lived vs short-lived and +#' high-quality or low-quality w.r.t. quality scores. #' #' @param prepared Output of `prepare_turnover_for_dose_response()` with weights. #' @param fit Output of `doseResponseFit()`. @@ -128,12 +101,8 @@ classify_turnover_fit <- function(prepared, fit, feature_data, as.data.frame(classification, stringsAsFactors = FALSE) } -#' Attach the per-protein score columns to the turnover fit result. -#' -#' Row order and row count of `fit` are preserved, so the significance -#' filtering and downloads that already read `ComparisonResult` are unaffected; -#' the classification's extra (unfitted) proteins stay in the classification -#' table only. +#' Attach per-protein confidence score columns to the turnover fit statistical +#' result. #' #' @param fit Output of `doseResponseFit()`. #' @param classification Output of `classify_turnover_fit()`.