Feat peptide weighting

NAMESPACE

@@ -51,6 +51,7 @@ importFrom(MSstatsPTM,SkylinetoMSstatsPTMFormat)
 importFrom(MSstatsPTM,SpectronauttoMSstatsPTMFormat)
 importFrom(MSstatsPTM,dataProcessPlotsPTM)
 importFrom(MSstatsPTM,groupComparisonPlotsPTM)
+importFrom(MSstatsResponse,calculatePeptideWeights)
 importFrom(MSstatsResponse,calculateTurnoverRatios)
 importFrom(MSstatsResponse,doseResponseFit)
 importFrom(MSstatsResponse,futureExperimentSimulation)

R/MSstatsShiny.R

@@ -38,7 +38,7 @@
 #' @importFrom MSstats MSstatsQualityMetricsPlot
 #' @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
+#' @importFrom MSstatsResponse futureExperimentSimulation run_tpr_simulation plot_tpr_power_curve calculateTurnoverRatios calculatePeptideWeights
 #' @importFrom utils capture.output head packageVersion write.csv
 #' @importFrom stats aggregate
 #' @importFrom methods is

R/constants.R

@@ -146,6 +146,7 @@ NAMESPACE_QC = list(
   reference_norm = "reference_norm",
   remove_norm_channel = "remove_norm_channel",
   features_used = "features_used",
+  assign_feature_weights = "assign_feature_weights",
   n_feat = "n_feat",
   cens_int = "censInt",
   null1 = "null1",
@@ -175,6 +176,8 @@ NAMESPACE_QC = list(
   standards_type_section = "standards_type_section",
   reference_norm_panel = "reference_norm_panel",
   lf_options_panel = "lf_options_panel",
+  feature_subset_panel = "feature_subset_panel",
+  feature_weights_panel = "feature_weights_panel",
   features_topn_panel = "features_topn_panel",
   censoring_section = "censoring_section",
   mbi_panel = "mbi_panel",

R/module-qc-ui.R

@@ -102,21 +102,33 @@ qcUI <- function(id) {
         # Feature subset, missing-value handling, imputation and summary method: label-free branch
         shinyjs::hidden(div(
           id = ns(NAMESPACE_QC$lf_options_panel),
-
-          # features
-
-          #h4("Used features"),
-          radioButtons(ns("features_used"),
-                       label = h4("Feature subset",class = "icon-wrapper",icon("question-circle", lib = "font-awesome"),
-                                  div("What features to use in \
-                                   summarization. All features or a subset of \
-                                   features can be used.", class = "icon-tooltip")),
-                       c("Use all features" = "all", "Use top N features" = "topN",
-                         "Remove uninformative features & outliers" = "highQuality")),
-          #),
+          div(
+            id = ns(NAMESPACE_QC$feature_subset_panel),
+            radioButtons(ns("features_used"),
+                         label = h4("Feature subset",class = "icon-wrapper",icon("question-circle", lib = "font-awesome"),
+                                    div("What features to use in \
+                                     summarization. All features or a subset of \
+                                     features can be used.", class = "icon-tooltip")),
+                         c("Use all features" = "all", "Use top N features" = "topN",
+                           "Remove uninformative features & outliers" = "highQuality")),
+            shinyjs::hidden(div(
+              id = ns(NAMESPACE_QC$features_topn_panel),
+              uiOutput(ns("features"))
+            ))
+          ),
           shinyjs::hidden(div(
-            id = ns(NAMESPACE_QC$features_topn_panel),
-            uiOutput(ns("features"))
+            id = ns(NAMESPACE_QC$feature_weights_panel),
+            h4("Feature Weighting"),
+            checkboxInput(
+              ns(NAMESPACE_QC$assign_feature_weights),
+              label = tags$div("Assign feature weights", class = "icon-wrapper",
+                               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.",
+                                   class = "icon-tooltip")),
+              value = FALSE
+            )
           )),
           #uiOutput("features"),
           tags$hr(),

R/module-statmodel-server.R

@@ -362,8 +362,10 @@ statmodelServer = function(id, parent_session, loadpage_input, qc_input,
             ratios <- turnover_ratios()
             increasing <- isTRUE(input[[NAMESPACE_STATMODEL$modeling_response_curve_increasing_trend]])
             dia_prepared <- prepare_turnover_for_dose_response(ratios, increasing = increasing)
+            turnover_weights <- if ("weight" %in% colnames(dia_prepared)) dia_prepared$weight else NULL
             response_results <- doseResponseFit(
               data      = dia_prepared,
+              weights   = turnover_weights,
               increasing    = increasing,
               transform_dose = FALSE,
               ratio_response = FALSE,
@@ -476,6 +478,7 @@ statmodelServer = function(id, parent_session, loadpage_input, qc_input,
             req(turnover_ratios())
             increasing <- isTRUE(input[[NAMESPACE_STATMODEL$modeling_response_curve_increasing_trend]])
             dia_prepared <- prepare_turnover_for_dose_response(turnover_ratios(), add_zero_timepoint = TRUE, increasing = increasing)
+            turnover_weights <- if ("weight" %in% colnames(dia_prepared)) dia_prepared$weight else NULL
           } else {
             meta <- condition_metadata()
             req(!is.null(meta) && "DoseVal" %in% colnames(meta))
@@ -496,6 +499,8 @@ statmodelServer = function(id, parent_session, loadpage_input, qc_input,
                 data = dia_prepared,
                 protein_name = input[[NAMESPACE_STATMODEL$visualization_which_protein]],
                 drug_name = "time",
+                weights = turnover_weights,
+                show_weights = !is.null(turnover_weights),
                 ratio_response = FALSE,
                 show_ic50 = TRUE,
                 add_ci = FALSE,

R/qc-server-sidebar.R

@@ -55,6 +55,19 @@ qc_show_log_section <- function(dda_dia, bio, template) {
   qc_show_lf(dda_dia, bio) && !isTRUE(template == TEMPLATES$protein_turnover)
 }
 
+#' Show the feature-subset radio: label-free branch, but hidden for the
+#' protein-turnover template (which always summarizes on all features).
+#' @noRd
+qc_show_feature_subset <- function(dda_dia, bio, template) {
+  qc_show_lf(dda_dia, bio) && !isTRUE(template == TEMPLATES$protein_turnover)
+}
+
+#' Show the feature-weighting checkbox: protein-turnover template only.
+#' @noRd
+qc_show_feature_weights <- function(template) {
+  isTRUE(template == TEMPLATES$protein_turnover)
+}
+
 #' Show the profile-plot options (summary toggle, feature legend).
 #' @noRd
 qc_show_profileplot_options <- function(type1) {
@@ -175,6 +188,21 @@ register_qc_visibility_observers <- function(input, session, loadpage_input, app
     )
   })
 
+  observe({
+    shinyjs::toggle(
+      NAMESPACE_QC$feature_subset_panel,
+      condition = qc_show_feature_subset(
+        loadpage_input()$DDA_DIA, loadpage_input()$BIO, app_template()
+      )
+    )
+  })
+  observe({
+    shinyjs::toggle(
+      NAMESPACE_QC$feature_weights_panel,
+      condition = qc_show_feature_weights(app_template())
+    )
+  })
+
   # Plot-tab option panels.
   observe({
     shinyjs::toggle(

R/qc-server-turnover.R

@@ -91,6 +91,16 @@ register_qc_turnover <- function(input, output, session, app_template, get_data,
     turnover_ratios()
   }, ignoreInit = TRUE)
 
+  turnover_ratios_display <- reactive({
+    ratios <- turnover_ratios()
+    req(ratios)
+    if (isTRUE(input[[NAMESPACE_QC$assign_feature_weights]]) && nrow(ratios) > 0) {
+      calculatePeptideWeights(ratios)
+    } else {
+      ratios
+    }
+  })
+
   output$turnover_ratios_panel <- renderUI({
     req(!is.null(app_template) && !is.null(app_template()) &&
           app_template() == TEMPLATES$protein_turnover)
@@ -106,24 +116,24 @@ register_qc_turnover <- function(input, output, session, app_template, get_data,
   })
 
   output$turnover_ratios_table_ui <- renderUI({
-    req(turnover_ratios())
+    req(turnover_ratios_display())
     ns <- session$ns
     enable("download_turnover_ratios")
     dataTableOutput(ns("turnover_ratios_table"))
   })
 
   output$turnover_ratios_table <- renderDataTable({
-    turnover_ratios()
-  })
+    turnover_ratios_display()
+  }, options = list(scrollX = TRUE))
 
   output$download_turnover_ratios <- downloadHandler(
     filename = function() {
       paste0("Turnover_Ratios-", Sys.Date(), ".csv")
     },
     content = function(file) {
-      write.csv(turnover_ratios(), file, row.names = FALSE)
+      write.csv(turnover_ratios_display(), file, row.names = FALSE)
     }
   )
 
-  turnover_ratios
+  turnover_ratios_display
 }

R/statmodel-server-comparisons.R

@@ -366,7 +366,7 @@ build_all_pair_contrast = function(input, condition_list, contrast, comp_list, r
 #' @param ratios Data frame from calculateTurnoverRatios (Protein, GROUP,
 #'   H_frac and L_frac columns required).
 #' @param add_zero_timepoint If TRUE, inserts a synthetic dose=0 / response=0
-#'   row for any (protein[, BaseSequence]) group missing one.
+#'   row for any peptide group missing one.
 #' @param increasing If TRUE (default), models synthesis from H_frac
 #'   (increasing trend over time). If FALSE, models degradation from L_frac
 #'   (decreasing trend over time).
@@ -383,6 +383,9 @@ prepare_turnover_for_dose_response <- function(ratios, add_zero_timepoint = FALS
   if ("BaseSequence" %in% colnames(result)) {
     keep_cols <- c(keep_cols, "BaseSequence")
   }
+  if ("weight" %in% colnames(result)) {
+    keep_cols <- c(keep_cols, "weight")
+  }
   result <- result[, keep_cols]
 
   if (add_zero_timepoint) {
@@ -403,6 +406,9 @@ prepare_turnover_for_dose_response <- function(ratios, add_zero_timepoint = FALS
       # Synthesis (H_frac): heavy fraction is 0 at t=0 (no incorporation yet).
       # Degradation (L_frac): light fraction is 1 at t=0 (pre-existing pool intact).
       zero_rows$response <- if (isTRUE(increasing)) 0 else 1
+      if ("weight" %in% keep_cols) {
+        zero_rows$weight <- 1
+      }
       result <- rbind(zero_rows[, keep_cols], result)
     }
   }

R/statmodel-server-download-code.R

@@ -13,6 +13,11 @@ generate_analysis_code = function(qc_input, loadpage_input, comp_mat, input, app
     codes = paste(codes, "\n# Set up dose response analysis\n", sep = "")
     codes = paste(codes, "library(MSstatsResponse)\n", sep = "")
 
+    if (isTRUE(app_template == TEMPLATES$protein_turnover)) {
+      codes = paste(codes, build_turnover_analysis_code(qc_input, comp_mat, increasing), sep = "")
+      return(codes)
+    }
+
     # Serialize the contrast matrix as a data frame
     codes = paste(codes, "group_metadata = data.frame(\n", sep = "")
     codes = paste(codes, "  GROUP = c(\"", paste(comp_mat$GROUP, collapse = "\",\""), "\"),\n", sep = "")
@@ -148,4 +153,121 @@ generate_analysis_code = function(qc_input, loadpage_input, comp_mat, input, app
   }
 
   return(codes)
+}
+
+
+#' Generate reproducible code for the protein-turnover dose-response pipeline.
+#'
+#' Mirrors the app's turnover flow: calculateTurnoverRatios (with the
+#' user-entered 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.
+#'
+#' @param qc_input The QC module input list (tracer_* numerics and
+#'   assign_feature_weights checkbox live here).
+#' @param comp_mat The turnover contrast matrix (GROUP + TimeVal columns); its
+#'   GROUP column supplies the condition names the tracer constants are keyed by.
+#' @param increasing Logical passed through to the fit / visualization.
+#' @return A character scalar of R code, appended after `library(MSstatsResponse)`.
+#' @noRd
+build_turnover_analysis_code <- function(qc_input, comp_mat, increasing) {
+  conditions <- as.character(comp_mat$GROUP)
+  weighting  <- isTRUE(qc_input[[NAMESPACE_QC$assign_feature_weights]])
+
+  # Serialize the tracer constants keyed by original condition name, matching
+  # the app (calculateTurnoverRatios parses these names into timepoints).
+  tracer_vals <- vapply(conditions, function(cond) {
+    val <- qc_input[[paste0("tracer_", make.names(cond))]]
+    if (is.null(val)) 1.0 else as.numeric(val)
+  }, numeric(1))
+  tracer_pairs <- paste0("  \"", conditions, "\" = ", tracer_vals, collapse = ",\n")
+
+  code <- paste0(
+    "\n# Tracer constants entered per condition on the data-processing page\n",
+    "tracer_constants = c(\n", tracer_pairs, "\n)\n",
+
+    "\n# Calculate turnover (Heavy/Light) ratios. Use protein-level data when any\n",
+    "# condition has replicate runs; otherwise fall back to feature-level data.\n",
+    "pld = summarized$ProteinLevelData\n",
+    "samples_per_condition = tapply(pld$RUN, pld$GROUP, function(x) length(unique(x)))\n",
+    "if (any(samples_per_condition > 1, na.rm = TRUE)) {\n",
+    "  turnover_ratios = calculateTurnoverRatios(\n",
+    "    summarized$ProteinLevelData,\n",
+    "    channel_col = \"LABEL\", heavy_label = \"H\", light_label = \"L\",\n",
+    "    time_col = \"GROUP\", peptide_col = \"Protein\", protein_col = \"Protein\",\n",
+    "    intensity_col = \"LogIntensities\", run_col = \"RUN\",\n",
+    "    agg_function = max, normalize_tracer = TRUE, tracer_constants = tracer_constants)\n",
+    "} else {\n",
+    "  turnover_ratios = calculateTurnoverRatios(\n",
+    "    summarized$FeatureLevelData,\n",
+    "    channel_col = \"LABEL\", heavy_label = \"H\", light_label = \"L\",\n",
+    "    time_col = \"GROUP\", peptide_col = \"PEPTIDE\", protein_col = \"PROTEIN\",\n",
+    "    intensity_col = \"INTENSITY\", run_col = \"RUN\",\n",
+    "    agg_function = max, normalize_tracer = TRUE, tracer_constants = tracer_constants)\n",
+    "}\n"
+  )
+
+  if (weighting) {
+    code <- paste0(
+      code,
+      "\n# Assign per-peptide quality weights (coverage, intensity, monotonicity,\n",
+      "# validity). Adds a 'weight' column used to down-weight low-quality peptides.\n",
+      "turnover_ratios = calculatePeptideWeights(turnover_ratios)\n"
+    )
+  }
+
+  # Column mapping matches prepare_turnover_for_dose_response(). Real turnover
+  # designs include a 0hr baseline, so no synthetic t=0 anchor is needed here.
+  frac_col <- if (isTRUE(increasing)) "H_frac" else "L_frac"
+  target_cols <- if (weighting) {
+    "c(\"protein\", \"drug\", \"dose\", \"response\", \"BaseSequence\", \"weight\")"
+  } else {
+    "c(\"protein\", \"drug\", \"dose\", \"response\", \"BaseSequence\")"
+  }
+
+  code <- paste0(
+    code,
+    "\n# Map columns to the dose-response format\n",
+    "frac_col = \"", frac_col, "\"\n",
+    "prepared_data = turnover_ratios[!is.na(turnover_ratios[[frac_col]]), ]\n",
+    "prepared_data$protein  = as.character(prepared_data$Protein)\n",
+    "prepared_data$drug     = \"time\"\n",
+    "prepared_data$dose     = as.numeric(prepared_data$TimeVal)\n",
+    "prepared_data$response = prepared_data[[frac_col]]\n",
+    "keep_cols = intersect(", target_cols, ", colnames(prepared_data))\n",
+    "prepared_data = prepared_data[, keep_cols, drop = FALSE]\n"
+  )
+
+  weights_arg <- if (weighting) "  weights = prepared_data$weight,\n" else ""
+
+  code <- paste0(
+    code,
+    "\n# Fit turnover time-course curves\n",
+    "response_results = doseResponseFit(\n",
+    "  data = prepared_data,\n",
+    weights_arg,
+    "  increasing = ", increasing, ",\n",
+    "  transform_dose = FALSE,\n",
+    "  ratio_response = FALSE,\n",
+    "  precalculated_ratios = TRUE\n)\n",
+
+    "\n# Visualize a single protein's turnover curve\n",
+    "visualizeResponseProtein(\n",
+    "  data = prepared_data,\n",
+    "  protein_name = \"Enter protein name here\",\n",
+    "  drug_name = \"time\",\n",
+    if (weighting) "  weights = prepared_data$weight,\n  show_weights = TRUE,\n" else "",
+    "  ratio_response = FALSE,\n",
+    "  show_ic50 = TRUE,\n",
+    "  add_ci = FALSE,\n",
+    "  transform_dose = FALSE,\n",
+    "  n_samples = 1000,\n",
+    "  increasing = ", increasing, ",\n",
+    "  precalculated_ratios = TRUE,\n",
+    "  color_by = \"BaseSequence\",\n",
+    "  target_response = 0.5\n)\n"
+  )
+
+  code
 }

R/statmodel-server-visualization.R

@@ -238,6 +238,7 @@ create_download_plot_handler <- function(output, input, contrast, preprocess_dat
               }
               increasing <- isTRUE(input[[NAMESPACE_STATMODEL$modeling_response_curve_increasing_trend]])
               dia_prepared <- prepare_turnover_for_dose_response(ratios, add_zero_timepoint = TRUE, increasing = increasing)
+              turnover_weights <- if ("weight" %in% colnames(dia_prepared)) dia_prepared$weight else NULL
             } else {
               if (isTRUE(app_template() == TEMPLATES$chemoproteomics)) {
                 meta <- tryCatch(condition_metadata(), error = function(e) NULL)
@@ -260,6 +261,8 @@ create_download_plot_handler <- function(output, input, contrast, preprocess_dat
                 data = dia_prepared,
                 protein_name = input[[NAMESPACE_STATMODEL$visualization_which_protein]],
                 drug_name = "time",
+                weights = turnover_weights,
+                show_weights = !is.null(turnover_weights),
                 ratio_response = FALSE,
                 show_ic50 = TRUE,
                 add_ci = FALSE,