Add classify turnover proteins

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)

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

R/module-qc-ui.R

@@ -125,7 +125,10 @@ 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 down-weight low-quality peptides in \
+                                the curve fit and enable per-protein \
+                                confidence score and turnover classification.",
                                    class = "icon-tooltip")),
               value = FALSE
             )

R/module-statmodel-server.R

@@ -377,7 +377,32 @@ statmodelServer = function(id, parent_session, loadpage_input, qc_input,
               ratio_response = FALSE,
               precalculated_ratios = TRUE
             )
-            list(ComparisonResult = response_results)
+
+            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 +492,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

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 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
 #'   lives here).
@@ -347,8 +349,32 @@ 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",
+      "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",

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
 

R/statmodel-server-turnover-confidence.R

@@ -0,0 +1,140 @@
+#' 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 and calculateConfidence read.
+#' @noRd
+TURNOVER_CONFIDENCE_FEATURE_COLUMNS <- c("PROTEIN", "PEPTIDE", "LABEL",
+                                         "INTENSITY", "GROUP")
+
+#' 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.
+#' @noRd
+turnover_weights_present <- function(prepared) {
+  !is.null(prepared) && NROW(prepared) > 0 &&
+    "weight" %in% colnames(prepared)
+}
+
+#' Checks if confidence scoring and classification can be applied
+#'
+#' 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.
+#' @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.
+#'
+#' @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
+}
+
+#' Prepare feature-level data for the QC-score / heavy-peptide counts.
+#'
+#' 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.
+#' @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 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()`.
+#' @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 per-protein confidence score columns to the turnover fit statistical 
+#' result.
+#'
+#' @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.")
+}

R/statmodel-ui-results.R

@@ -36,7 +36,8 @@ create_results_tables <- function(ns) {
     ),
     conditionalPanel(
       condition = "input['loadpage-BIO']!=='PTM'",
-      uiOutput(ns("table_results"))
+      uiOutput(ns("table_results")),
+      uiOutput(ns("turnover_confidence_results"))
     )
   )
 }