feat(protein turnover): Create initial protein turnover functions for ratio calculation

DESCRIPTION

@@ -32,6 +32,8 @@ Imports:
     stats,
     parallel,
     data.table,
+    tidyr,
+    stringr,
     plotly
 Suggests:
     MSstats,

NAMESPACE

@@ -1,12 +1,14 @@
 # Generated by roxygen2: do not edit by hand
 
 export(MSstatsPrepareDoseResponseFit)
+export(calculateTurnoverRatios)
 export(convertGroupToNumericDose)
 export(doseResponseFit)
 export(futureExperimentSimulation)
 export(plot_tpr_power_curve)
 export(predictIC50)
 export(predictIC50Parallel)
+export(selectTopNPeptides)
 export(run_tpr_simulation)
 export(visualizeResponseProtein)
 import(dplyr)
@@ -19,8 +21,12 @@ importFrom(dplyr,filter)
 importFrom(dplyr,group_by)
 importFrom(dplyr,if_else)
 importFrom(dplyr,mutate)
+importFrom(dplyr,pull)
+importFrom(dplyr,rename)
 importFrom(dplyr,select)
+importFrom(dplyr,slice_head)
 importFrom(dplyr,summarise)
+importFrom(dplyr,ungroup)
 importFrom(ggplot2,aes)
 importFrom(ggplot2,annotate)
 importFrom(ggplot2,element_blank)
@@ -54,3 +60,7 @@ importFrom(stats,pf)
 importFrom(stats,quantile)
 importFrom(stats,rlnorm)
 importFrom(stats,rnorm)
+importFrom(stringr,str_detect)
+importFrom(stringr,str_extract)
+importFrom(stringr,str_remove_all)
+importFrom(tidyr,pivot_wider)

R/Fit_Isotonic_Regression.R

@@ -5,9 +5,14 @@
 #'
 #' @param data Protein-level data, formatted with MSstatsPreparedoseResponseFit().
 #' @param weights Optional numeric vector of weights. Defaults to equal weights.
-#' @param increasing Logical. If TRUE, fits a non-decreasing model. If FALSE, fits non-increasing.
+#' @param increasing Logical or character. If TRUE, fits a non-decreasing model. If FALSE, fits non-increasing.
+#'   If "both", fits both increasing and decreasing models and selects the one with better fit (lower SSE).
+#'   Recommended for exploratory analysis, but ~2x slower. Default is FALSE.
 #' @param transform_dose Logical. If TRUE, applies log10(dose + 1) transformation. Default is TRUE.
 #' @param ratio_response Logical. If TRUE, converts log2 abundance to ratios relative to DMSO. Default is FALSE.
+#' @param precalculated_ratios Logical. If TRUE, assumes the response column contains pre-calculated ratios
+#'   (not log2 values) and skips internal ratio calculation. This is useful for protein turnover data
+#'   or other datasets where ratios are computed externally. Default is FALSE.
 #'
 #' @return A data frame with protein-wise F-test results and BH-adjusted p-values.
 #'
@@ -43,15 +48,33 @@
 #'   ratio_response = FALSE
 #' )
 #'
-#' # View results
-#' print(interaction_results)
+#' # Example 2: Fit both directions and select better
+#' interaction_both <- doseResponseFit(
+#'   data = example_data,
+#'   increasing = "both",
+#'   transform_dose = TRUE,
+#'   ratio_response = FALSE
+#' )
 #'
-#' # Check significant interactions
-#' significant <- interaction_results[interaction_results$adj.pvalue < 0.05, ]
-#' print(paste("Found", nrow(significant), "significant interactions"))
+#' # Check which direction was selected for each protein
+#' table(interaction_both$direction)
+#'
+#' # Example 3: Using pre-calculated ratios (e.g., from protein turnover)
+#' # Assume you have pre-calculated ratio data
+#' turnover_data <- example_data
+#' turnover_data$response <- 2^turnover_data$response /
+#'                           mean(2^turnover_data$response[turnover_data$dose == 0])
+#'
+#' interaction_results_precalc <- doseResponseFit(
+#'   data = turnover_data,
+#'   increasing = FALSE,
+#'   transform_dose = TRUE,
+#'   ratio_response = FALSE,
+#'   precalculated_ratios = TRUE
+#' )
 #'
 #' \dontrun{
-#' # Example 2: Full dataset analysis
+#' # Example 4: Full dataset analysis
 #' full_results <- doseResponseFit(
 #'   data = prepared_data,
 #'   increasing = FALSE,
@@ -67,7 +90,8 @@
 doseResponseFit = function(data, weights = NULL,
                            increasing = FALSE,
                            transform_dose = TRUE,
-                           ratio_response = FALSE) {
+                           ratio_response = FALSE,
+                           precalculated_ratios = FALSE) {
 
   drug_list = unique(data$drug[data$drug != "DMSO"])
 
@@ -81,6 +105,7 @@ doseResponseFit = function(data, weights = NULL,
       F_statistic = numeric(),
       P_value = numeric(),
       adj.pvalue = numeric(),
+      direction = character(),
       stringsAsFactors = FALSE
     ))
   }
@@ -104,7 +129,6 @@ doseResponseFit = function(data, weights = NULL,
           if (is.null(weights)) {
             w = rep(1, length(y))
           } else if (length(weights) == nrow(data)) {
-            # If weights provided for full dataset, subset them
             protein_drug_idx = which(data$protein == protein_list[i] &
                                        data$drug %in% c("DMSO", drug_type))
             w = weights[protein_drug_idx]
@@ -115,21 +139,64 @@ doseResponseFit = function(data, weights = NULL,
             w = rep(1, length(y))
           }
 
-          #w = if (is.null(weights)) rep(1, length(y)) else weights
-
-          fit = fitIsotonicRegression(
-            x = x, y = y, w = w,
-            increasing = increasing,
-            transform_x = transform_dose,
-            ratio_y = ratio_response,
-            test_significance = TRUE
-          )
-
-          results_temp = fit$f_test
-          results_temp$protein = protein_list[i]
-          results_temp$drug = drug_type
-          results_temp = results_temp[, c("protein", "drug", setdiff(names(results_temp), c("protein", "drug")))]
-          results_list[[i]] = results_temp
+          # Handle "both" option for increasing parameter
+          if (is.character(increasing) && increasing == "both") {
+            # Fit both directions - keep both results for FDR correction
+            fit_dec = fitIsotonicRegression(
+              x = x, y = y, w = w,
+              increasing = FALSE,
+              transform_x = transform_dose,
+              ratio_y = ratio_response,
+              precalculated_ratios = precalculated_ratios,
+              test_significance = TRUE
+            )
+
+            fit_inc = fitIsotonicRegression(
+              x = x, y = y, w = w,
+              increasing = TRUE,
+              transform_x = transform_dose,
+              ratio_y = ratio_response,
+              precalculated_ratios = precalculated_ratios,
+              test_significance = TRUE
+            )
+
+            # Store both results temporarily for FDR correction
+            results_dec = fit_dec$f_test
+            results_dec$protein = protein_list[i]
+            results_dec$drug = drug_type
+            results_dec$direction = "decreasing"
+            results_dec$F_statistic_alt = fit_inc$f_test$F_statistic  # Store alternative F for later comparison
+
+            results_inc = fit_inc$f_test
+            results_inc$protein = protein_list[i]
+            results_inc$drug = drug_type
+            results_inc$direction = "increasing"
+            results_inc$F_statistic_alt = fit_dec$f_test$F_statistic
+
+            # Add both to results list
+            results_list[[length(results_list) + 1]] = results_dec
+            results_list[[length(results_list) + 1]] = results_inc
+
+          } else {
+            # Original behavior - single direction
+            fit = fitIsotonicRegression(
+              x = x, y = y, w = w,
+              increasing = increasing,
+              transform_x = transform_dose,
+              ratio_y = ratio_response,
+              precalculated_ratios = precalculated_ratios,
+              test_significance = TRUE
+            )
+            chosen_direction = if(is.logical(increasing) && increasing) "increasing" else "decreasing"
+
+            results_temp = fit$f_test
+            results_temp$protein = protein_list[i]
+            results_temp$drug = drug_type
+            results_temp$direction = chosen_direction
+            results_temp = results_temp[, c("protein", "drug", "direction",
+                                            setdiff(names(results_temp), c("protein", "drug", "direction")))]
+            results_list[[i]] = results_temp
+          }
         })
       }, error = function(e) {
         warning(paste("Error for drug:", drug_type, "protein:", protein_list[i], ":", conditionMessage(e)))
@@ -139,7 +206,62 @@ doseResponseFit = function(data, weights = NULL,
     # Combine and adjust p-values for this drug
     if (length(results_list) > 0) {
       results_df = do.call(rbind, results_list)
-      results_df$adj.pvalue = p.adjust(results_df$P_value, method = "BH")
+
+      # If using "both", we need to handle FDR correction and selection
+      if (is.character(increasing) && increasing == "both") {
+        # Separate FDR correction for increasing and decreasing
+        inc_idx = which(results_df$direction == "increasing")
+        dec_idx = which(results_df$direction == "decreasing")
+
+        # Initialize adj.pvalue column
+        results_df$adj.pvalue = NA
+
+        # FDR correction within each direction
+        if (length(inc_idx) > 0) {
+          results_df$adj.pvalue[inc_idx] = p.adjust(results_df$P_value[inc_idx], method = "BH")
+        }
+        if (length(dec_idx) > 0) {
+          results_df$adj.pvalue[dec_idx] = p.adjust(results_df$P_value[dec_idx], method = "BH")
+        }
+
+        # Now select the better direction for each protein
+        proteins_unique = unique(results_df$protein)
+        final_results = list()
+
+        for (prot in proteins_unique) {
+          prot_rows = results_df[results_df$protein == prot, ]
+
+          if (nrow(prot_rows) == 2) {
+            # Compare F-statistics to select better fit
+            dec_row = prot_rows[prot_rows$direction == "decreasing", ]
+            inc_row = prot_rows[prot_rows$direction == "increasing", ]
+
+            if (dec_row$F_statistic >= inc_row$F_statistic) {
+              selected_row = dec_row
+            } else {
+              selected_row = inc_row
+            }
+
+            # Remove the F_statistic_alt column
+            selected_row$F_statistic_alt = NULL
+            final_results[[length(final_results) + 1]] = selected_row
+          } else {
+            # Only one direction (shouldn't happen, but handle it)
+            prot_rows$F_statistic_alt = NULL
+            final_results[[length(final_results) + 1]] = prot_rows
+          }
+        }
+
+        results_df = do.call(rbind, final_results)
+
+      } else {
+        # Standard FDR correction when direction is fixed
+        results_df$adj.pvalue = p.adjust(results_df$P_value, method = "BH")
+      }
+
+      # Reorder columns to put direction after drug
+      results_df = results_df[, c("protein", "drug", "direction",
+                                  setdiff(names(results_df), c("protein", "drug", "direction")))]
       all_results[[drug_type]] = results_df
     }
   }
@@ -154,11 +276,13 @@ doseResponseFit = function(data, weights = NULL,
       F_statistic = numeric(),
       P_value = numeric(),
       adj.pvalue = numeric(),
+      direction = character(),
       stringsAsFactors = FALSE
     ))
   }
 
   final_results = do.call(rbind, all_results)
+  rownames(final_results) = NULL
   return(final_results)
 }
 
@@ -173,6 +297,8 @@ doseResponseFit = function(data, weights = NULL,
 #' @param increasing Logical. If TRUE, fits a non-decreasing model. If FALSE, fits non-increasing.
 #' @param transform_x Logical. If TRUE, applies log10(x + 1) transformation. Default is TRUE.
 #' @param ratio_y Logical. If TRUE, converts log2 abundance to ratios relative to DMSO. Default is FALSE.
+#' @param precalculated_ratios Logical. If TRUE, assumes y contains pre-calculated ratios and skips
+#'   internal ratio calculation. Default is FALSE.
 #' @param test_significance Logical. If TRUE, performs F-test to assess significance.
 #'
 #' @return A list representing the isotonic regression model (class = "isotonic_model").
@@ -182,6 +308,7 @@ fitIsotonicRegression = function(x, y, w = rep(1, length(y)),
                                  increasing = FALSE,
                                  transform_x = TRUE,
                                  ratio_y = FALSE,
+                                 precalculated_ratios = FALSE,
                                  test_significance = FALSE) {
   stopifnot(length(x) == length(y), length(y) == length(w))
 
@@ -191,11 +318,17 @@ fitIsotonicRegression = function(x, y, w = rep(1, length(y)),
     x = log10(x + 1)
   }
 
-  if (ratio_y) {
+  # Handle ratio calculation
+  if (precalculated_ratios) {
+    # y already contains ratios - use as-is
+    # No transformation needed
+  } else if (ratio_y) {
+    # Original behavior: convert log2 to ratios
     y = 2^y
-    baseline = mean(y[x == 0], na.rm = TRUE)
+    baseline = mean(y[x_org == 0], na.rm = TRUE)
     y = y / baseline
   }
+  # Otherwise y stays as-is (log2 values)
 
   order_idx = order(x)
   x_sorted = x[order_idx]
@@ -263,16 +396,15 @@ fitIsotonicRegression = function(x, y, w = rep(1, length(y)),
     y_pred_new = stats::approx(x, y_final, xout = x, rule = 2)$y
   }
 
-  #y_pred_new = stats::approx(x, y_final, xout = x, rule = 2)$y
-
   result = list(
     x = x,
     y_pred = y_pred_new,
     original_x = x_org[order_idx],
     original_y = y[order_idx],
     increasing = increasing,
     transform_x = transform_x,
-    ratio_y = ratio_y
+    ratio_y = ratio_y,
+    precalculated_ratios = precalculated_ratios
   )
 
   if (test_significance) {