diff --git a/.DS_Store b/.DS_Store
deleted file mode 100644
index 259c1f8..0000000
Binary files a/.DS_Store and /dev/null differ
diff --git a/DESCRIPTION b/DESCRIPTION
index e3e40c1..db8921e 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -32,6 +32,8 @@ Imports:
     stats,
     parallel,
     data.table,
+    tidyr,
+    stringr,
     plotly
 Suggests:
     MSstats,
diff --git a/NAMESPACE b/NAMESPACE
index 24d4bd3..d1f1db2 100644
--- a/NAMESPACE
+++ b/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)
diff --git a/R/Fit_Isotonic_Regression.R b/R/Fit_Isotonic_Regression.R
index 48724f5..59bfe87 100644
--- a/R/Fit_Isotonic_Regression.R
+++ b/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,8 +396,6 @@ 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,
@@ -272,7 +403,8 @@ fitIsotonicRegression = function(x, y, w = rep(1, length(y)),
     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) {
diff --git a/R/IC50_prediction.R b/R/IC50_prediction.R
index 2a4859e..06b25a1 100644
--- a/R/IC50_prediction.R
+++ b/R/IC50_prediction.R
@@ -4,6 +4,7 @@
                            increasing,
                            transform_dose,
                            ratio_response,
+                           precalculated_ratios,
                            bootstrap,
                            prot,
                            drug_type,
@@ -22,6 +23,7 @@
       IC50 = NA,
       IC50_lower_bound = NA,
       IC50_upper_bound = NA,
+      direction = NA,
       stringsAsFactors = FALSE
     ))
   }
@@ -35,6 +37,7 @@
       IC50 = NA,
       IC50_lower_bound = NA,
       IC50_upper_bound = NA,
+      direction = NA,
       stringsAsFactors = FALSE
     ))
   }
@@ -43,8 +46,13 @@
   x = x[order_idx]
   y = y[order_idx]
 
-  y_unlog = 2^y
-  baseline = mean(y_unlog[x == 0], na.rm = TRUE)
+  # Calculate baseline for validation
+  if (precalculated_ratios) {
+    baseline = mean(y[x == 0], na.rm = TRUE)
+  } else {
+    y_unlog = 2^y
+    baseline = mean(y_unlog[x == 0], na.rm = TRUE)
+  }
 
   # Check if baseline is valid
   if (is.na(baseline) || !is.finite(baseline)) {
@@ -55,14 +63,15 @@
       IC50 = NA,
       IC50_lower_bound = NA,
       IC50_upper_bound = NA,
+      direction = NA,
       stringsAsFactors = FALSE
     ))
   }
 
-
-  y_ratio = y_unlog / baseline
-
-  if (!ratio_response){
+  # Adjust target response based on data type
+  if (precalculated_ratios) {
+    adjusted_target_response = target_response
+  } else if (!ratio_response) {
     dmso_mean = mean(y[x == 0], na.rm = TRUE)
     adjusted_target_response = dmso_mean + log2(target_response)
   } else {
@@ -71,17 +80,70 @@
 
   result = NULL
   ic50_est = NA
+  chosen_direction = NA
 
   tryCatch({
-    fit_try = fitIsotonicRegression(x, y,
-                                    increasing = increasing,
-                                    transform_x = transform_dose,
-                                    ratio_y = ratio_response,
-                                    test_significance = FALSE)
-    ic50_est = PredictIC50(fit_try, target_response = adjusted_target_response)
+    # Handle "both" option
+    if (is.character(increasing) && increasing == "both") {
+      # Fit both directions
+      fit_dec = fitIsotonicRegression(x, y,
+                                      increasing = FALSE,
+                                      transform_x = transform_dose,
+                                      ratio_y = ratio_response,
+                                      precalculated_ratios = precalculated_ratios,
+                                      test_significance = FALSE)
+
+      fit_inc = fitIsotonicRegression(x, y,
+                                      increasing = TRUE,
+                                      transform_x = transform_dose,
+                                      ratio_y = ratio_response,
+                                      precalculated_ratios = precalculated_ratios,
+                                      test_significance = FALSE)
+
+      # Try to predict IC50 for both
+      ic50_dec = PredictIC50(fit_dec, target_response = adjusted_target_response)
+      ic50_inc = PredictIC50(fit_inc, target_response = adjusted_target_response)
+
+      # Select the one that successfully predicts IC50
+      if (!is.na(ic50_dec) && is.na(ic50_inc)) {
+        ic50_est = ic50_dec
+        chosen_direction = "decreasing"
+        fit_try = fit_dec
+        increasing_final = FALSE
+      } else if (is.na(ic50_dec) && !is.na(ic50_inc)) {
+        ic50_est = ic50_inc
+        chosen_direction = "increasing"
+        fit_try = fit_inc
+        increasing_final = TRUE
+      } else if (!is.na(ic50_dec) && !is.na(ic50_inc)) {
+        # Both work - choose based on better fit (could use SSE if stored)
+        # Default to decreasing for now
+        ic50_est = ic50_dec
+        chosen_direction = "decreasing"
+        fit_try = fit_dec
+        increasing_final = FALSE
+      } else {
+        ic50_est = NA
+        chosen_direction = NA
+        increasing_final = FALSE
+      }
+    } else {
+      # Original single-direction behavior
+      fit_try = fitIsotonicRegression(x, y,
+                                      increasing = increasing,
+                                      transform_x = transform_dose,
+                                      ratio_y = ratio_response,
+                                      precalculated_ratios = precalculated_ratios,
+                                      test_significance = FALSE)
+      ic50_est = PredictIC50(fit_try, target_response = adjusted_target_response)
+      chosen_direction = if(is.logical(increasing) && increasing) "increasing" else "decreasing"
+      increasing_final = increasing
+    }
+
     ic50_est = ifelse(is.na(ic50_est), NA, ic50_est)
   }, error = function(e) {
     ic50_est <- NA
+    chosen_direction <- NA
   })
 
   if (is.na(ic50_est)) {
@@ -90,21 +152,30 @@
       drug = drug_type,
       IC50 = NA,
       IC50_lower_bound = NA,
-      IC50_upper_bound = NA
+      IC50_upper_bound = NA,
+      direction = chosen_direction
     )
   } else {
-        if (transform_dose) {
+    if (transform_dose) {
       ic50 = -log10(ic50_est)
     } else {
       ic50 = ic50_est
     }
 
     if (bootstrap) {
-      if(ratio_response){
+      if (precalculated_ratios) {
+        bootstrap_res = bootstrapIC50Precalculated(
+          dose = x, response = y,
+          n_samples = n_samples, alpha = alpha,
+          increasing = increasing_final,
+          target_response = adjusted_target_response,
+          transform_x = transform_dose
+        )
+      } else if (ratio_response) {
         bootstrap_res = bootstrapIC50(
           dose = x, response = y,
           n_samples = n_samples, alpha = alpha,
-          increasing = increasing,
+          increasing = increasing_final,
           target_response = target_response,
           transform_x = transform_dose
         )
@@ -112,7 +183,7 @@
         bootstrap_res = bootstrapIC50LogScale(
           x = x, y = y,
           n_samples = n_samples, alpha = alpha,
-          increasing = increasing,
+          increasing = increasing_final,
           target_response = target_response,
           transform_x = transform_dose
         )
@@ -128,7 +199,8 @@
       drug = drug_type,
       IC50 = ic50,
       IC50_lower_bound = lower,
-      IC50_upper_bound = upper
+      IC50_upper_bound = upper,
+      direction = chosen_direction
     )
   }
 
@@ -246,7 +318,6 @@ bootstrapIC50 = function(dose, response, n_samples = 1000, alpha = 0.10,
 #' @param alpha Significance level for confidence interval (default = 0.05).
 #' @param increasing Logical. Fit non-decreasing if TRUE.
 #' @param target_response Numeric value for response level (default = 0.5).
-#'
 #' @param transform_x Logical. If TRUE, applies log10(x + 1) transformation. Default = TRUE.
 #' @return List with mean IC50, CI bounds, and transformed estimates.
 #' @importFrom stats quantile
@@ -313,3 +384,71 @@ bootstrapIC50LogScale = function(x, y, n_samples = 1000, alpha = 0.05,
     ci_upper_transform = ci_upper_transform
   ))
 }
+
+#' Bootstrap IC50 with pre-calculated ratios
+#'
+#' @param dose Numeric vector of dose values.
+#' @param response Numeric vector of response values (pre-calculated ratios).
+#' @param n_samples Number of bootstrap samples. Default = 1000.
+#' @param alpha Significance level for confidence interval. Default = 0.10.
+#' @param increasing Logical. Fit non-decreasing if TRUE.
+#' @param target_response Numeric value for response level.
+#' @param transform_x Logical. If TRUE, applies log10(x + 1) transformation. Default = TRUE.
+#'
+#' @return List with mean IC50, CI bounds, and transformed estimates.
+#' @importFrom stats quantile
+#' @importFrom dplyr filter select mutate group_by summarise arrange distinct
+bootstrapIC50Precalculated = function(dose, response, n_samples = 1000, alpha = 0.10,
+                                      increasing = FALSE, target_response = 0.5,
+                                      transform_x = TRUE) {
+  ic50_vals = numeric(n_samples)
+  df = data.frame(dose = dose, response = response)
+
+  for (i in seq_len(n_samples)) {
+    boot_df = df %>%
+      dplyr::group_by(dose) %>%
+      dplyr::sample_frac(size = 1, replace = TRUE) %>%
+      dplyr::ungroup()
+
+    x_sample = boot_df$dose
+    y_sample = boot_df$response
+
+    tryCatch({
+      fit_sample = fitIsotonicRegression(x_sample, y_sample,
+                                         increasing = increasing,
+                                         transform_x = transform_x,
+                                         ratio_y = FALSE,
+                                         precalculated_ratios = TRUE,
+                                         test_significance = FALSE)
+      ic50_est = PredictIC50(fit_sample, target_response = target_response)
+      ic50_vals[i] = ifelse(is.na(ic50_est), NA, ic50_est)
+    }, error = function(e) {
+      ic50_vals[i] = NA
+    })
+  }
+
+  ic50_values_clean = ic50_vals[!is.na(ic50_vals) & ic50_vals != 0]
+  ci_lower = stats::quantile(ic50_values_clean, alpha / 2, na.rm = TRUE)
+  ci_upper = stats::quantile(ic50_values_clean, 1 - alpha / 2, na.rm = TRUE)
+  mean_ic50 = mean(ic50_values_clean, na.rm = TRUE)
+
+  if (transform_x) {
+    mean_ic50_transform = -log10(mean_ic50)
+    ci_lower_transform = -log10(ci_lower)
+    ci_upper_transform = -log10(ci_upper)
+  } else {
+    mean_ic50_transform = mean_ic50
+    ci_lower_transform = ci_lower
+    ci_upper_transform = ci_upper
+  }
+
+  return(list(
+    ic50_values = ic50_vals,
+    mean_ic50 = mean_ic50,
+    ci_lower = ci_lower,
+    ci_upper = ci_upper,
+    mean_ic50_transform = mean_ic50_transform,
+    ci_lower_transform = ci_lower_transform,
+    ci_upper_transform = ci_upper_transform
+  ))
+}
diff --git a/R/PredictIC50.R b/R/PredictIC50.R
index 8b9dc85..2c48b1e 100644
--- a/R/PredictIC50.R
+++ b/R/PredictIC50.R
@@ -3,17 +3,20 @@
 #' @param data A data frame with columns: protein, drug, dose, response.
 #' @param n_samples Number of bootstrap samples. Default = 1000.
 #' @param alpha Confidence level. Default = 0.10.
-#' @param increasing Logical. If TRUE, fit a non-decreasing trend. Default = FALSE.
+#' @param increasing Logical or character. If TRUE, fit a non-decreasing trend. If FALSE, fit non-increasing.
+#'   If "both", fits both directions and selects the better fit. Default = FALSE.
 #' @param transform_dose Logical. If TRUE, applies log10(dose + 1) transformation. Default = TRUE.
 #' @param ratio_response Logical. If TRUE, use ratio response; else use log2 scale. Default = TRUE.
+#' @param precalculated_ratios Logical. If TRUE, assumes response column contains pre-calculated ratios.
+#'   Default is FALSE.
 #' @param bootstrap Logical. If FALSE, skip confidence interval bootstrap estimation and only return IC50. Default = TRUE.
-#'@param BPPARAM A \code{BiocParallelParam} for parallel processing. The recommended usage
-#'   is to \emph{register} a backend once (e.g., \code{register(MulticoreParam(workers=4))} on
-#'   Linux/macOS or \code{register(SnowParam(workers=4, type="SOCK"))} on Windows) and pass
-#'   \code{BPPARAM = bpparam()}. Default \code{bpparam()}.
-#' @param target_response Numeric, the response fraction (e.g., 0.5, 0.25, 0.75). Default = 0.5.
+#' @param BPPARAM A \code{BiocParallelParam} for parallel processing. Default \code{bpparam()}.
+#' @param target_response Numeric, the response fraction (e.g., 0.5, 0.25, 0.75, 1.5).
+#'   For decreasing curves: use values < 1 (e.g., 0.5 for IC50).
+#'   For increasing curves: use values > 1 (e.g., 1.5 for 50% increase).
+#'   Default = 0.5 (auto-adjusts to 1.5 for increasing curves).
 #'
-#' @return A data frame with columns: protein, drug, IC50, IC50_lower_bound, IC50_upper_bound.
+#' @return A data frame with columns: protein, drug, IC50, IC50_lower_bound, IC50_upper_bound, direction.
 #'
 #' @examples
 #' # Load example data
@@ -47,8 +50,16 @@
 #' )
 #' print(ic50_quick)
 #'
+#' # Example 2: Fit both increasing and decreasing curves
+#' ic50_both <- predictIC50(
+#'   data = example_data,
+#'   increasing = "both",
+#'   bootstrap = FALSE
+#' )
+#' print(ic50_both)
+#'
 #' \dontrun{
-#' # Example 2: Full IC50 estimation with bootstrap confidence intervals
+#' # Example 3: Full IC50 estimation with bootstrap confidence intervals
 #' ic50_results <- predictIC50(
 #'   data = prepared_data,
 #'   n_samples = 1000,
@@ -57,7 +68,7 @@
 #'   bootstrap = TRUE
 #' )
 #'
-#' # Example 3: Parallel processing for large datasets
+#' # Example 4: Parallel processing for large datasets
 #' library(BiocParallel)
 #' ic50_parallel <- predictIC50(
 #'   data = prepared_data,
@@ -66,12 +77,31 @@
 #'   bootstrap = TRUE
 #' )
 #'
-#' # Example 4: IC50 at different response levels (IC25, IC75)
+#' # Example 5: IC50 at different response levels (IC25, IC75)
 #' ic25_results <- predictIC50(
 #'   data = prepared_data,
 #'   target_response = 0.25,
 #'   bootstrap = TRUE
 #' )
+#'
+#' # Example 6: Increasing curves with EC50 at 1.5x baseline
+#' ec50_results <- predictIC50(
+#'   data = prepared_data,
+#'   increasing = TRUE,
+#'   target_response = 1.5,
+#'   bootstrap = TRUE
+#' )
+#'
+#' # Example 7: Using pre-calculated ratios
+#' turnover_data <- prepared_data
+#' turnover_data$response <- 2^turnover_data$response /
+#'                           mean(2^turnover_data$response[turnover_data$dose == 0])
+#'
+#' ic50_precalc <- predictIC50(
+#'   data = turnover_data,
+#'   precalculated_ratios = TRUE,
+#'   bootstrap = TRUE
+#' )
 #' }
 #'
 #' @export
@@ -84,9 +114,28 @@ predictIC50 = function(data,
                        increasing = FALSE,
                        transform_dose = TRUE,
                        ratio_response = TRUE,
+                       precalculated_ratios = FALSE,
                        bootstrap = TRUE,
                        BPPARAM = bpparam(),
-                       target_response = 0.5) {
+                       target_response = NULL) {
+
+  # Handle target_response
+  if (is.null(target_response)) {
+    if (precalculated_ratios) {
+      stop("When using precalculated_ratios = TRUE, you must explicitly specify target_response.\n",
+           "For decreasing responses (e.g., inhibition): use values < 1 (e.g., 0.5 for IC50)\n",
+           "For increasing responses (e.g., activation): use values > 1 (e.g., 1.5 for EC50)")
+    } else {
+      # Set defaults for package-calculated ratios
+      if (is.logical(increasing) && increasing == TRUE) {
+        target_response = 1.5
+        message("Using target_response = 1.5 for increasing curves (EC50 at 50% increase)")
+      } else {
+        target_response = 0.5
+        message("Using target_response = 0.5 for decreasing curves (IC50 at 50% reduction)")
+      }
+    }
+  }
 
   # Create list of protein-drug combinations to process
   loop_list = data %>%
@@ -106,6 +155,7 @@ predictIC50 = function(data,
       increasing = increasing,
       transform_dose = transform_dose,
       ratio_response = ratio_response,
+      precalculated_ratios = precalculated_ratios,
       bootstrap = bootstrap,
       prot = temp[[2]],
       drug_type = temp[[1]],
diff --git a/R/Visualize_Isotonic_Fit.R b/R/Visualize_Isotonic_Fit.R
index 05b790c..d0376e6 100644
--- a/R/Visualize_Isotonic_Fit.R
+++ b/R/Visualize_Isotonic_Fit.R
@@ -3,14 +3,20 @@
 #' @param data Protein-level dataset (e.g., output of MSstatsPrepareDoseResponseFit).
 #' @param protein_name Character. Protein name to plot.
 #' @param drug_name Character. Drug name to plot.
-#' @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 and selects better fit.
 #' @param ratio_response Logical. If TRUE, compute IC50 on ratio scale; if FALSE, use log2 intensities.
+#' @param precalculated_ratios Logical. If TRUE, assumes response contains pre-calculated ratios. Default is FALSE.
 #' @param transform_dose Logical. If TRUE, applies log10(dose + 1). Default is TRUE.
 #' @param show_ic50 Logical. If TRUE, adds vertical line and annotation for IC50.
+#' @param target_response Numeric. Target response level for IC50/EC50 calculation. Default = 0.5.
 #' @param add_ci Logical. Include IC50 95% confidence interval bands if TRUE. Default is FALSE.
 #' @param n_samples Number of bootstrap samples if including confidence intervals. Default is 1000.
-#' @param alpha Alpha level for confidence intervals. Default is 0.05.
-#' @param y_lab Character. Label for the y-axis. Default is "Ratio Response".
+#' @param alpha Alpha level for confidence intervals. Default is 0.10.
+#' @param color_by Character. Column name to color points by (e.g., "replicate", "peptide"). Default is NULL.
+#' @param x_lab Character. Custom label for the x-axis. If NULL, uses default based on transform_dose.
+#' @param y_lab Character. Custom label for the y-axis. If NULL, uses default based on ratio_response.
+#' @param title Character. Custom plot title. If NULL, auto-generates from drug and protein names.
 #'
 #' @return A ggplot object.
 #'
@@ -45,22 +51,36 @@
 #' )
 #' print(plot1)
 #'
-#' # Example 2: Add IC50 annotation
+#' # Example 2: Add IC50 annotation with custom labels
 #' plot2 <- visualizeResponseProtein(
 #'   data = prepared_data,
 #'   protein_name = "PROTEIN_A",
 #'   drug_name = "Drug1",
 #'   ratio_response = TRUE,
 #'   show_ic50 = TRUE,
-#'   add_ci = FALSE
+#'   x_lab = "Drug Concentration (M)",
+#'   y_lab = "Relative Protein Abundance",
+#'   title = "Drug Response for PROTEIN_A"
 #' )
 #' print(plot2)
 #'
+#' # Example 3: Color points by replicate
+#' plot3 <- visualizeResponseProtein(
+#'   data = prepared_data,
+#'   protein_name = "PROTEIN_A",
+#'   drug_name = "Drug1",
+#'   ratio_response = TRUE,
+#'   show_ic50 = TRUE,
+#'   color_by = "replicate"
+#' )
+#' print(plot3)
+#'
 #' @export
 visualizeResponseProtein = function(data,
                                     protein_name,
                                     drug_name,
                                     ratio_response = TRUE,
+                                    precalculated_ratios = FALSE,
                                     transform_dose = TRUE,
                                     show_ic50 = TRUE,
                                     target_response = 0.5,
@@ -68,7 +88,10 @@ visualizeResponseProtein = function(data,
                                     n_samples = 1000,
                                     alpha = 0.10,
                                     increasing = FALSE,
-                                    y_lab = 'Ratio Response') {
+                                    color_by = NULL,
+                                    x_lab = NULL,
+                                    y_lab = NULL,
+                                    title = NULL) {
   # Subset to single protein + drug (+ DMSO)
   data_single = data %>%
     dplyr::filter(protein == protein_name & drug %in% c("DMSO", drug_name))
@@ -80,26 +103,59 @@ visualizeResponseProtein = function(data,
   x = x_org[order_idx]
   y = y_org[order_idx]
 
-  # Fit isotonic model using renamed function
-  fit = fitIsotonicRegression(
-    x = x, y = y,
-    increasing = increasing,
-    transform_x = transform_dose,
-    ratio_y = ratio_response,
-    test_significance = FALSE
-  )
+  # Handle "both" option by fitting and selecting better model
+  if (is.character(increasing) && increasing == "both") {
+    fit_dec = fitIsotonicRegression(
+      x = x, y = y,
+      increasing = FALSE,
+      transform_x = transform_dose,
+      ratio_y = ratio_response,
+      precalculated_ratios = precalculated_ratios,
+      test_significance = TRUE
+    )
+
+    fit_inc = fitIsotonicRegression(
+      x = x, y = y,
+      increasing = TRUE,
+      transform_x = transform_dose,
+      ratio_y = ratio_response,
+      precalculated_ratios = precalculated_ratios,
+      test_significance = TRUE
+    )
+
+    # Select better fit
+    if (fit_dec$f_test$SSE_Full <= fit_inc$f_test$SSE_Full) {
+      fit = fit_dec
+      increasing_final = FALSE
+    } else {
+      fit = fit_inc
+      increasing_final = TRUE
+    }
+  } else {
+    # Fit isotonic model
+    fit = fitIsotonicRegression(
+      x = x, y = y,
+      increasing = increasing,
+      transform_x = transform_dose,
+      ratio_y = ratio_response,
+      precalculated_ratios = precalculated_ratios,
+      test_significance = FALSE
+    )
+    increasing_final = increasing
+  }
 
   # Optional: bootstrap CIs
   ci_bounds = NULL
   if (add_ci) {
     suppressWarnings({
-    ic50_est = predictIC50(data_single,
-                           ratio_response = ratio_response,
-                           transform_dose = transform_dose,
-                           increasing = increasing,
-                           target_response = target_response,
-                           n_samples = n_samples,
-                           alpha = alpha)
+      ic50_est = predictIC50(data_single,
+                             ratio_response = ratio_response,
+                             precalculated_ratios = precalculated_ratios,
+                             transform_dose = transform_dose,
+                             increasing = increasing_final,
+                             target_response = target_response,
+                             n_samples = n_samples,
+                             alpha = alpha)
     })
     ci_bounds = list(
       ci_lower = ic50_est$IC50_lower_bound,
@@ -110,13 +166,18 @@ visualizeResponseProtein = function(data,
   p = plotIsotonic(
     fit = fit,
     ratio = ratio_response,
+    precalculated_ratios = precalculated_ratios,
     show_ic50 = show_ic50,
     target_response = target_response,
     ci = ci_bounds,
     drug_name = drug_name,
     protein_name = protein_name,
+    x_lab = x_lab,
     y_lab = y_lab,
-    transform_x = transform_dose
+    title = title,
+    transform_x = transform_dose,
+    color_by = color_by,
+    original_data = data_single
   )
 
   return(p)
@@ -125,17 +186,22 @@ visualizeResponseProtein = function(data,
 #' Plot Isotonic Regression Model
 #'
 #' @param fit A model object returned by fitIsotonicRegression().
-#' @param ratio Logical. If TRUE, shows plot on the ratio scale relative to DMSO (i.e. 0-1 scale). Default is FALSE.
+#' @param ratio Logical. If TRUE, shows plot on the ratio scale relative to DMSO. Default is FALSE.
+#' @param precalculated_ratios Logical. If TRUE, response values are pre-calculated ratios. Default is FALSE.
 #' @param show_ic50 Logical. If TRUE, adds vertical line and annotation for IC50.
+#' @param target_response Numeric. Target response level for IC50/EC50. Default = 0.5.
 #' @param drug_name Drug name for plotting data.
 #' @param protein_name Protein name for plot.
-#' @param x_lab Label for x-axis.
-#' @param y_lab Label for y-axis.
-#' @param title Title for the plot.
-#' @param ci Logical. Include IC50 95% confidence interval bands if TRUE. Default is FALSE.
+#' @param x_lab Character. Label for x-axis. If NULL, uses default based on transform_x.
+#' @param y_lab Character. Label for y-axis. If NULL, uses default based on ratio/precalculated_ratios.
+#' @param title Character. Title for the plot. If NULL, auto-generates.
+#' @param ci List with ci_lower and ci_upper. Include IC50 confidence interval bands if provided. Default is NULL.
 #' @param legend Logical. Show legend if TRUE.
 #' @param theme_style ggplot2 theme name to apply (default = "classic").
 #' @param original_label Logical. If TRUE, replace x-axis tick labels with original dose labels.
+#' @param transform_x Logical. Whether doses were log-transformed. Default = TRUE.
+#' @param color_by Character. Column name to color points by. Default is NULL.
+#' @param original_data Data frame containing original data with grouping variables.
 #'
 #' @return A ggplot object.
 #'
@@ -145,45 +211,49 @@ visualizeResponseProtein = function(data,
 #' @importFrom ggplot2 theme element_rect element_blank element_text
 plotIsotonic = function(fit,
                         ratio = TRUE,
+                        precalculated_ratios = FALSE,
                         show_ic50 = FALSE,
-                        target_response = target_response,
+                        target_response = 0.5,
                         drug_name = NULL,
                         protein_name = NULL,
-                        x_lab = expression(Log[10] ~ "[drug (M)]"),
-                        y_lab = "Log2 Intensity",
+                        x_lab = NULL,
+                        y_lab = NULL,
                         title = NULL,
                         ci = NULL,
                         legend = FALSE,
                         theme_style = "classic",
                         original_label = FALSE,
-                        transform_x = TRUE) {
+                        transform_x = TRUE,
+                        color_by = NULL,
+                        original_data = NULL) {
 
   # Construct title if not provided
-  if (is.null(title) && !is.null(drug_name) && !is.null(protein_name)) {
-    title = paste0(drug_name, ": ", protein_name, " Response")
-  } else if (is.null(title)) {
-    title = "Isotonic Regression Fit"
+  if (is.null(title)) {
+    if (!is.null(drug_name) && !is.null(protein_name)) {
+      title = paste0(drug_name, ": ", protein_name, " Response")
+    } else {
+      title = "Isotonic Regression Fit"
+    }
   }
 
-  if (ratio) {
-    y_lab = "Ratio relative to DMSO"
-  } else {
-    y_lab = expression(Log[2] ~ "Intensity")
+  # Set default labels if not provided
+  if (is.null(y_lab)) {
+    if (precalculated_ratios || ratio) {
+      y_lab = "Ratio relative to DMSO"
+    } else {
+      y_lab = expression(Log[2] ~ "Intensity")
+    }
   }
 
-  if (transform_x) {
-    x_lab = expression(Log[10] ~ "[drug (M)]")
-  } else {
-    x_lab = "x"
+  if (is.null(x_lab)) {
+    if (transform_x) {
+      x_lab = expression(Log[10] ~ "[drug (M)]")
+    } else {
+      x_lab = "Dose"
+    }
   }
 
   # Prepare data
-  #fit_df = data.frame(dose = log10(fit$x), y_pred = fit$y_pred)
-  #orig_df = data.frame(x = log10(fit$x), y = fit$original_y)
-
-  # Prepare data - fit$x is in the correct scale based on transform_x
-  fit_df = data.frame(dose = fit$x, y_pred = fit$y_pred)
-  orig_df = data.frame(x = fit$x, y = fit$original_y)
   if (transform_x) {
     fit_df = data.frame(dose = log10(fit$x), y_pred = fit$y_pred)
     orig_df = data.frame(x = log10(fit$x), y = fit$original_y)
@@ -192,16 +262,43 @@ plotIsotonic = function(fit,
     orig_df = data.frame(x = fit$x, y = fit$original_y)
   }
 
-  # Base plot
-  model_fit_plot = ggplot2::ggplot() +
-    ggplot2::geom_point(data = orig_df,
-                        ggplot2::aes(x = x, y = y, color = "Observed Data"), size = 2) +
-    ggplot2::geom_line(data = fit_df,
-                       ggplot2::aes(x = dose, y = y_pred, color = "Isotonic Regression Fit"), size = 1.2) +
-    ggplot2::scale_color_manual(name = NULL,
-                                values = c("Observed Data" = "black",
-                                           "Isotonic Regression Fit" = "orange")) +
-    ggplot2::labs(x = x_lab, y = y_lab, title = title) +
+  # Add color grouping variable if requested
+  if (!is.null(color_by) && !is.null(original_data)) {
+    if (!color_by %in% names(original_data)) {
+      stop(paste("Column", color_by, "not found in data"))
+    }
+    order_idx = order(original_data$dose)
+    orig_df[[color_by]] = original_data[[color_by]][order_idx]
+  }
+
+  # Base plot with conditional coloring
+  if (!is.null(color_by) && !is.null(original_data)) {
+    # Color points by grouping variable
+    model_fit_plot = ggplot2::ggplot() +
+      ggplot2::geom_point(data = orig_df,
+                          ggplot2::aes(x = x, y = y, color = .data[[color_by]]),
+                          size = 2) +
+      ggplot2::geom_line(data = fit_df,
+                         ggplot2::aes(x = dose, y = y_pred),
+                         color = "orange", size = 1.2) +
+      ggplot2::labs(x = x_lab, y = y_lab, title = title, color = color_by)
+  } else {
+    # Default: black points, orange line
+    model_fit_plot = ggplot2::ggplot() +
+      ggplot2::geom_point(data = orig_df,
+                          ggplot2::aes(x = x, y = y, color = "Observed Data"),
+                          size = 2) +
+      ggplot2::geom_line(data = fit_df,
+                         ggplot2::aes(x = dose, y = y_pred, color = "Isotonic Regression Fit"),
+                         size = 1.2) +
+      ggplot2::scale_color_manual(name = NULL,
+                                  values = c("Observed Data" = "black",
+                                             "Isotonic Regression Fit" = "orange")) +
+      ggplot2::labs(x = x_lab, y = y_lab, title = title)
+  }
+
+  # Apply common theme
+  model_fit_plot = model_fit_plot +
     ggplot2::ylim(0, NA) +
     ggplot2::theme_minimal() +
     ggplot2::theme(
@@ -212,34 +309,29 @@ plotIsotonic = function(fit,
       plot.title = ggplot2::element_text(hjust = 0.5, face = "bold", size = 16)
     )
 
-  if (!legend) {
+  # Handle legend display
+  if (!legend && is.null(color_by)) {
     model_fit_plot = model_fit_plot +
       ggplot2::theme(legend.position = "none")
   }
 
   # Optional IC50
   if (show_ic50) {
-    if (ratio) {
-      target_response = target_response
+    if (precalculated_ratios || ratio) {
+      adjusted_target = target_response
     } else {
       if (transform_x) {
-        #dmso_idx = which(orig_df$x == 0)  # log10(0+1) = 0
         dmso_idx = which(is.infinite(orig_df$x) & orig_df$x < 0)
-        target_response = mean(orig_df$y[dmso_idx], na.rm = TRUE) + log2(target_response)
-
       } else {
-        dmso_idx = which(fit$original_x == 0)  # Raw dose = 0
-        target_response = mean(orig_df$y[dmso_idx], na.rm = TRUE) + log2(target_response)
-        }
-      #dmso_idx = which(orig_df$x == -Inf)
-      #target_response = mean(orig_df$y[dmso_idx], na.rm = TRUE) - 1
+        dmso_idx = which(fit$original_x == 0)
+      }
+      adjusted_target = mean(orig_df$y[dmso_idx], na.rm = TRUE) + log2(target_response)
     }
 
-    ic50_pred = PredictIC50(fit, target_response = target_response)
+    ic50_pred = PredictIC50(fit, target_response = adjusted_target)
 
-    # fix approx warning for handling duplicate x values
+    # Handle duplicate x values for approx
     if (length(unique(fit$x)) < length(fit$x)) {
-      # Aggregate y values for duplicate x values before interpolation
       x_unique = unique(fit$x)
       y_aggregated = sapply(x_unique, function(xi) {
         mean(fit$y_pred[fit$x == xi])
@@ -250,27 +342,16 @@ plotIsotonic = function(fit,
         y_ic50 = stats::approx(fit$x, fit$y_pred, xout = ic50_pred)$y
       })
     }
-    #y_ic50 = stats::approx(fit$x, fit$y_pred, xout = ic50_pred)$y
-
 
     if (transform_x) {
-      # ic50_pred is on log10(dose+1) scale, convert to pIC50 for label
-      ic50_pred_transform = -log10(ic50_pred) #pIC50
+      ic50_pred_transform = -log10(ic50_pred)
       ic50_label = paste("pIC50 =", round(ic50_pred_transform, 2))
-      x_pos = log10(ic50_pred)  # Already on correct scale for plotting
+      x_pos = log10(ic50_pred)
     } else {
-      # ic50_pred is on raw scale, use as-is
       ic50_label = paste("IC50 =", round(ic50_pred, 2))
       x_pos = ic50_pred
     }
 
-   # model_fit_plot = model_fit_plot +
-    #  ggplot2::geom_point(ggplot2::aes(x = log10(ic50_pred), y = y_ic50),
-    #                      shape = 23, size = 3.5, fill = "red", color = "black") +
-    #  ggplot2::annotate("text", x = log10(ic50_pred) + 0.35, y = y_ic50,
-     #                   label = paste("pIC50 =", round(ic50_pred_transform, 2)),
-     #                   vjust = -0.5, hjust = 0.5, size = 4, fontface = "bold")
-
     model_fit_plot = model_fit_plot +
       ggplot2::geom_point(ggplot2::aes(x = x_pos, y = y_ic50),
                           shape = 23, size = 3.5, fill = "red", color = "black") +
@@ -278,30 +359,28 @@ plotIsotonic = function(fit,
                         label = ic50_label,
                         vjust = -0.5, hjust = 0.5, size = 4, fontface = "bold")
 
-    # optional display IC50 confidence interval
-    if (show_ic50 && !is.null(ci)) {
+    # Optional display IC50 confidence interval
+    if (show_ic50 && !is.null(ci) && !is.na(ci$ci_lower) && !is.na(ci$ci_upper)) {
       if (transform_x) {
-        # CI bounds are in pIC50, convert to log10(dose+1) scale for plotting
-        ci_lower_dose = 10^(-ci$ci_lower)  # pIC50 to dose
+        ci_lower_dose = 10^(-ci$ci_lower)
         ci_upper_dose = 10^(-ci$ci_upper)
-        ci_lower_x = log10(ci_lower_dose)  # dose to log10(dose+1)
+        ci_lower_x = log10(ci_lower_dose)
         ci_upper_x = log10(ci_upper_dose)
       } else {
-        # CI bounds are already in raw scale
         ci_lower_x = ci$ci_lower
         ci_upper_x = ci$ci_upper
       }
 
       model_fit_plot = model_fit_plot +
         ggplot2::geom_segment(ggplot2::aes(x = ci_lower_x, xend = ci_lower_x,
-                                           y = 0, yend = target_response),
+                                           y = 0, yend = adjusted_target),
                               linetype = "dotted", color = "gray40", linewidth = 0.8) +
         ggplot2::geom_segment(ggplot2::aes(x = ci_upper_x, xend = ci_upper_x,
-                                           y = 0, yend = target_response),
+                                           y = 0, yend = adjusted_target),
                               linetype = "dotted", color = "gray40", linewidth = 0.8) +
         ggplot2::annotate("rect",
                           xmin = ci_lower_x, xmax = ci_upper_x,
-                          ymin = 0, ymax = target_response,
+                          ymin = 0, ymax = adjusted_target,
                           alpha = 0.1, fill = "red")
     }
   }
diff --git a/R/protein_turnover_ratio_helper.R b/R/protein_turnover_ratio_helper.R
new file mode 100644
index 0000000..65041c3
--- /dev/null
+++ b/R/protein_turnover_ratio_helper.R
@@ -0,0 +1,209 @@
+#' Calculate turnover ratios from MSstats FeatureLevelData
+#'
+#' @param feature_data Data frame from MSstats dataProcess()$FeatureLevelData
+#' @param channel_col Character. Name of column containing Heavy/Light labels. Default = "LABEL"
+#' @param heavy_label Character. Value in channel_col indicating heavy channel. Default = "H"
+#' @param light_label Character. Value in channel_col indicating light channel. Default = "L"
+#' @param time_col Character. Column containing timepoint information. Default = "GROUP"
+#' @param peptide_col Character. Column containing peptide sequences. Default = "PEPTIDE"
+#' @param protein_col Character. Column containing protein identifiers. Default = "PROTEIN"
+#' @param intensity_col Character. Column with intensity values. Default = "INTENSITY"
+#' @param run_col Character. Column identifying technical replicates. Default = "RUN"
+#' @param peptide_selector Optional function to select peptides.
+#'   Function should take a data frame (grouped by protein) and return filtered data frame.
+#'   Default = NULL (use all peptides).
+#' @param agg_function Function to aggregate duplicate peptide measurements (multiple charge states,
+#'   transitions, etc.). Default = max (takes highest signal).
+#' @param normalize_tracer Logical. If TRUE, normalize by tracer incorporation. Default = FALSE
+#' @param tracer_constants Named numeric vector. Tracer constants for each timepoint.
+#'   Required if normalize_tracer = TRUE
+#'
+#' @return Data frame with columns: Protein, BaseSequence, TimeVal, Run, Heavy, Light, Total, H_frac, L_frac
+#'
+#' @examples
+#' # Basic usage - all proteins, all peptides
+#' ratios <- calculateTurnoverRatios(
+#'   feature_data = quant_data$FeatureLevelData
+#' )
+#'
+#' # With peptide selection (top 10 by signal per protein)
+#' top10_selector <- function(df) {
+#'   top_peptides <- df %>%
+#'     group_by(BaseSequence) %>%
+#'     summarise(total_signal = sum(Intensity, na.rm = TRUE), .groups = "drop") %>%
+#'     arrange(desc(total_signal)) %>%
+#'     slice_head(n = 10) %>%
+#'     pull(BaseSequence)
+#'
+#'   df %>% filter(BaseSequence %in% top_peptides)
+#' }
+#'
+#' ratios_top10 <- calculateTurnoverRatios(
+#'   feature_data = quant_data$FeatureLevelData,
+#'   peptide_selector = top10_selector
+#' )
+#'
+#' # Use different aggregation function (e.g., median for robustness)
+#' ratios_median <- calculateTurnoverRatios(
+#'   feature_data = quant_data$FeatureLevelData,
+#'   agg_function = median
+#' )
+#'
+#' # With tracer normalization
+#' tracer_consts <- c("0hr" = 1.0, "1hr" = 0.95, "12hrs" = 0.85, "168hrs" = 0.75)
+#' ratios_norm <- calculateTurnoverRatios(
+#'   feature_data = quant_data$FeatureLevelData,
+#'   normalize_tracer = TRUE,
+#'   tracer_constants = tracer_consts
+#' )
+#'
+#' # Filter for specific protein after calculation
+#' protein_ratios <- ratios %>% filter(Protein == "A0A2K5TXF6")
+#'
+#' @export
+#' @importFrom dplyr filter mutate group_by summarise arrange slice_head pull ungroup select rename
+#' @importFrom tidyr pivot_wider
+#' @importFrom stringr str_remove_all str_extract str_detect
+calculateTurnoverRatios <- function(
+    feature_data,
+    channel_col = "LABEL",
+    heavy_label = "H",
+    light_label = "L",
+    time_col = "GROUP",
+    peptide_col = "PEPTIDE",
+    protein_col = "PROTEIN",
+    intensity_col = "INTENSITY",
+    run_col = "RUN",
+    peptide_selector = NULL,
+    agg_function = max,
+    normalize_tracer = FALSE,
+    tracer_constants = NULL
+) {
+
+  # Check required columns
+  required_cols <- c(protein_col, channel_col, time_col, peptide_col, intensity_col, run_col)
+  missing_cols <- setdiff(required_cols, colnames(feature_data))
+  if (length(missing_cols) > 0) {
+    stop("Missing required columns: ", paste(missing_cols, collapse = ", "))
+  }
+
+  # Process all proteins
+  df <- feature_data %>%
+    mutate(
+      Protein = .data[[protein_col]],
+      BaseSequence = str_remove_all(.data[[peptide_col]], "\\[.*?\\]"),
+      Label = .data[[channel_col]],
+      TimeVal = parse_timepoint(.data[[time_col]]),
+      Intensity = .data[[intensity_col]],
+      Run = .data[[run_col]]
+    ) %>%
+    filter(!is.na(TimeVal)) %>%
+    filter(Label %in% c(heavy_label, light_label))
+
+  if (nrow(df) == 0) {
+    warning("No data found matching Heavy/Light labels")
+    return(data.frame())
+  }
+
+  # Apply optional peptide selector per protein
+  if (!is.null(peptide_selector)) {
+    df <- df %>%
+      group_by(Protein) %>%
+      group_modify(~ peptide_selector(.x)) %>%
+      ungroup()
+  }
+
+  # Aggregate duplicates (multiple features/transitions for same peptide)
+  # This can happen when there are multiple charge states or modifications
+  df <- df %>%
+    group_by(Protein, BaseSequence, TimeVal, Run, Label) %>%
+    summarise(Intensity = agg_function(Intensity, na.rm = TRUE), .groups = "drop")
+
+  # Pivot to wide format (keep replicates separate)
+  df_wide <- df %>%
+    select(Protein, BaseSequence, TimeVal, Run, Label, Intensity) %>%
+    pivot_wider(names_from = Label, values_from = Intensity)
+
+  # Rename heavy/light columns if they're not "Heavy" and "Light"
+  if (heavy_label != "Heavy") {
+    df_wide <- df_wide %>% rename(Heavy = all_of(heavy_label))
+  }
+  if (light_label != "Light") {
+    df_wide <- df_wide %>% rename(Light = all_of(light_label))
+  }
+
+  df_wide <- df_wide %>%
+    filter(!is.na(Heavy) & !is.na(Light)) %>%
+    mutate(
+      Total = Heavy + Light,
+      H_frac = Heavy / Total,
+      L_frac = Light / Total
+    )
+
+  # Optional tracer normalization
+  if (normalize_tracer) {
+    if (is.null(tracer_constants)) {
+      stop("normalize_tracer = TRUE but tracer_constants was not provided")
+    }
+
+    df_wide <- df_wide %>%
+      mutate(
+        tracer_factor = tracer_constants[as.character(TimeVal)],
+        H_frac = H_frac / tracer_factor
+      )
+  }
+
+  df_wide
+}
+
+#' Parse timepoint strings to numeric hours
+#'
+#' @param time_strings Character vector of timepoint labels
+#'
+#' @return Numeric vector of hours
+#'
+#' @examples
+#' parse_timepoint(c("0hr", "1hr", "12hrs", "168hrs"))
+#' # Returns: 0 1 12 168
+#'
+#' @keywords internal
+parse_timepoint <- function(time_strings) {
+  # Extract numeric part
+  numeric_part <- as.numeric(str_extract(time_strings, "^[0-9]+"))
+
+  # Handle different units
+  is_days <- str_detect(time_strings, "d|day")
+  is_weeks <- str_detect(time_strings, "w|week")
+
+  hours <- numeric_part
+  hours[is_days] <- numeric_part[is_days] * 24
+  hours[is_weeks] <- numeric_part[is_weeks] * 24 * 7
+
+  return(hours)
+}
+
+#' Example peptide selector: Top N by total signal
+#'
+#' @param df Data frame with BaseSequence and Intensity columns (for one protein)
+#' @param n Number of top peptides to select
+#'
+#' @return Filtered data frame
+#'
+#' @examples
+#' # Use as peptide_selector
+#' ratios <- calculateTurnoverRatios(
+#'   feature_data = quant_data$FeatureLevelData,
+#'   peptide_selector = function(df) selectTopNPeptides(df, n = 10)
+#' )
+#'
+#' @export
+selectTopNPeptides <- function(df, n = 10) {
+  top_peptides <- df %>%
+    group_by(BaseSequence) %>%
+    summarise(total_signal = sum(Intensity, na.rm = TRUE), .groups = "drop") %>%
+    arrange(desc(total_signal)) %>%
+    slice_head(n = n) %>%
+    pull(BaseSequence)
+
+  df %>% filter(BaseSequence %in% top_peptides)
+}
diff --git a/man/DoseResponseFit.Rd b/man/DoseResponseFit.Rd
index 3033264..c5db8a8 100644
--- a/man/DoseResponseFit.Rd
+++ b/man/DoseResponseFit.Rd
@@ -9,7 +9,8 @@ doseResponseFit(
   weights = NULL,
   increasing = FALSE,
   transform_dose = TRUE,
-  ratio_response = FALSE
+  ratio_response = FALSE,
+  precalculated_ratios = FALSE
 )
 }
 \arguments{
@@ -17,11 +18,17 @@ doseResponseFit(
 
 \item{weights}{Optional numeric vector of weights. Defaults to equal weights.}
 
-\item{increasing}{Logical. If TRUE, fits a non-decreasing model. If FALSE, fits non-increasing.}
+\item{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.}
 
 \item{transform_dose}{Logical. If TRUE, applies log10(dose + 1) transformation. Default is TRUE.}
 
 \item{ratio_response}{Logical. If TRUE, converts log2 abundance to ratios relative to DMSO. Default is FALSE.}
+
+\item{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.}
 }
 \value{
 A data frame with protein-wise F-test results and BH-adjusted p-values.
@@ -62,15 +69,34 @@ interaction_results <- doseResponseFit(
   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 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
+)
 
-# Check significant interactions
-significant <- interaction_results[interaction_results$adj.pvalue < 0.05, ]
-print(paste("Found", nrow(significant), "significant interactions"))
 
 \dontrun{
-# Example 2: Full dataset analysis
+# Example 4: Full dataset analysis
 full_results <- doseResponseFit(
   data = prepared_data,
   increasing = FALSE,
diff --git a/man/PredictIC50.Rd b/man/PredictIC50.Rd
index 70b7b8a..3d956f6 100644
--- a/man/PredictIC50.Rd
+++ b/man/PredictIC50.Rd
@@ -11,9 +11,10 @@ predictIC50(
   increasing = FALSE,
   transform_dose = TRUE,
   ratio_response = TRUE,
+  precalculated_ratios = FALSE,
   bootstrap = TRUE,
   BPPARAM = bpparam(),
-  target_response = 0.5
+  target_response = NULL
 )
 }
 \arguments{
@@ -23,23 +24,27 @@ predictIC50(
 
 \item{alpha}{Confidence level. Default = 0.10.}
 
-\item{increasing}{Logical. If TRUE, fit a non-decreasing trend. Default = FALSE.}
+\item{increasing}{Logical or character. If TRUE, fit a non-decreasing trend. If FALSE, fit non-increasing.
+If "both", fits both directions and selects the better fit. Default = FALSE.}
 
 \item{transform_dose}{Logical. If TRUE, applies log10(dose + 1) transformation. Default = TRUE.}
 
 \item{ratio_response}{Logical. If TRUE, use ratio response; else use log2 scale. Default = TRUE.}
 
+\item{precalculated_ratios}{Logical. If TRUE, assumes response column contains pre-calculated ratios.
+Default is FALSE.}
+
 \item{bootstrap}{Logical. If FALSE, skip confidence interval bootstrap estimation and only return IC50. Default = TRUE.}
 
-\item{BPPARAM}{A \code{BiocParallelParam} for parallel processing. The recommended usage
-is to \emph{register} a backend once (e.g., \code{register(MulticoreParam(workers=4))} on
-Linux/macOS or \code{register(SnowParam(workers=4, type="SOCK"))} on Windows) and pass
-\code{BPPARAM = bpparam()}. Default \code{bpparam()}.}
+\item{BPPARAM}{A \code{BiocParallelParam} for parallel processing. Default \code{bpparam()}.}
 
-\item{target_response}{Numeric, the response fraction (e.g., 0.5, 0.25, 0.75). Default = 0.5.}
+\item{target_response}{Numeric, the response fraction (e.g., 0.5, 0.25, 0.75, 1.5).
+For decreasing curves: use values < 1 (e.g., 0.5 for IC50).
+For increasing curves: use values > 1 (e.g., 1.5 for 50\% increase).
+Default = 0.5 (auto-adjusts to 1.5 for increasing curves).}
 }
 \value{
-A data frame with columns: protein, drug, IC50, IC50_lower_bound, IC50_upper_bound.
+A data frame with columns: protein, drug, IC50, IC50_lower_bound, IC50_upper_bound, direction.
 }
 \description{
 Predict IC50 (dose where response = target) for each protein and drug
@@ -76,8 +81,16 @@ ic50_quick <- predictIC50(
 )
 print(ic50_quick)
 
+# Example 2: Fit both increasing and decreasing curves
+ic50_both <- predictIC50(
+  data = example_data,
+  increasing = "both",
+  bootstrap = FALSE
+)
+print(ic50_both)
+
 \dontrun{
-# Example 2: Full IC50 estimation with bootstrap confidence intervals
+# Example 3: Full IC50 estimation with bootstrap confidence intervals
 ic50_results <- predictIC50(
   data = prepared_data,
   n_samples = 1000,
@@ -86,7 +99,7 @@ ic50_results <- predictIC50(
   bootstrap = TRUE
 )
 
-# Example 3: Parallel processing for large datasets
+# Example 4: Parallel processing for large datasets
 library(BiocParallel)
 ic50_parallel <- predictIC50(
   data = prepared_data,
@@ -95,12 +108,31 @@ ic50_parallel <- predictIC50(
   bootstrap = TRUE
 )
 
-# Example 4: IC50 at different response levels (IC25, IC75)
+# Example 5: IC50 at different response levels (IC25, IC75)
 ic25_results <- predictIC50(
   data = prepared_data,
   target_response = 0.25,
   bootstrap = TRUE
 )
+
+# Example 6: Increasing curves with EC50 at 1.5x baseline
+ec50_results <- predictIC50(
+  data = prepared_data,
+  increasing = TRUE,
+  target_response = 1.5,
+  bootstrap = TRUE
+)
+
+# Example 7: Using pre-calculated ratios
+turnover_data <- prepared_data
+turnover_data$response <- 2^turnover_data$response /
+                          mean(2^turnover_data$response[turnover_data$dose == 0])
+
+ic50_precalc <- predictIC50(
+  data = turnover_data,
+  precalculated_ratios = TRUE,
+  bootstrap = TRUE
+)
 }
 
 }
diff --git a/man/VisualizeResponseProtein.Rd b/man/VisualizeResponseProtein.Rd
index 7d89203..c0e6dd3 100644
--- a/man/VisualizeResponseProtein.Rd
+++ b/man/VisualizeResponseProtein.Rd
@@ -9,6 +9,7 @@ visualizeResponseProtein(
   protein_name,
   drug_name,
   ratio_response = TRUE,
+  precalculated_ratios = FALSE,
   transform_dose = TRUE,
   show_ic50 = TRUE,
   target_response = 0.5,
@@ -16,7 +17,10 @@ visualizeResponseProtein(
   n_samples = 1000,
   alpha = 0.1,
   increasing = FALSE,
-  y_lab = "Ratio Response"
+  color_by = NULL,
+  x_lab = NULL,
+  y_lab = NULL,
+  title = NULL
 )
 }
 \arguments{
@@ -28,19 +32,30 @@ visualizeResponseProtein(
 
 \item{ratio_response}{Logical. If TRUE, compute IC50 on ratio scale; if FALSE, use log2 intensities.}
 
+\item{precalculated_ratios}{Logical. If TRUE, assumes response contains pre-calculated ratios. Default is FALSE.}
+
 \item{transform_dose}{Logical. If TRUE, applies log10(dose + 1). Default is TRUE.}
 
 \item{show_ic50}{Logical. If TRUE, adds vertical line and annotation for IC50.}
 
+\item{target_response}{Numeric. Target response level for IC50/EC50 calculation. Default = 0.5.}
+
 \item{add_ci}{Logical. Include IC50 95\% confidence interval bands if TRUE. Default is FALSE.}
 
 \item{n_samples}{Number of bootstrap samples if including confidence intervals. Default is 1000.}
 
-\item{alpha}{Alpha level for confidence intervals. Default is 0.05.}
+\item{alpha}{Alpha level for confidence intervals. Default is 0.10.}
+
+\item{increasing}{Logical or character. If TRUE, fits a non-decreasing model. If FALSE, fits non-increasing.
+If "both", fits both and selects better fit.}
+
+\item{color_by}{Character. Column name to color points by (e.g., "replicate", "peptide"). Default is NULL.}
 
-\item{increasing}{Logical. If TRUE, fits a non-decreasing model. If FALSE, fits non-increasing.}
+\item{x_lab}{Character. Custom label for the x-axis. If NULL, uses default based on transform_dose.}
 
-\item{y_lab}{Character. Label for the y-axis. Default is "Ratio Response".}
+\item{y_lab}{Character. Custom label for the y-axis. If NULL, uses default based on ratio_response.}
+
+\item{title}{Character. Custom plot title. If NULL, auto-generates from drug and protein names.}
 }
 \value{
 A ggplot object.
@@ -79,15 +94,28 @@ plot1 <- visualizeResponseProtein(
 )
 print(plot1)
 
-# Example 2: Add IC50 annotation
+# Example 2: Add IC50 annotation with custom labels
 plot2 <- visualizeResponseProtein(
   data = prepared_data,
   protein_name = "PROTEIN_A",
   drug_name = "Drug1",
   ratio_response = TRUE,
   show_ic50 = TRUE,
-  add_ci = FALSE
+  x_lab = "Drug Concentration (M)",
+  y_lab = "Relative Protein Abundance",
+  title = "Drug Response for PROTEIN_A"
 )
 print(plot2)
 
+# Example 3: Color points by replicate
+plot3 <- visualizeResponseProtein(
+  data = prepared_data,
+  protein_name = "PROTEIN_A",
+  drug_name = "Drug1",
+  ratio_response = TRUE,
+  show_ic50 = TRUE,
+  color_by = "replicate"
+)
+print(plot3)
+
 }
diff --git a/man/bootstrapIC50Precalculated.Rd b/man/bootstrapIC50Precalculated.Rd
new file mode 100644
index 0000000..05f8442
--- /dev/null
+++ b/man/bootstrapIC50Precalculated.Rd
@@ -0,0 +1,37 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/IC50_prediction.R
+\name{bootstrapIC50Precalculated}
+\alias{bootstrapIC50Precalculated}
+\title{Bootstrap IC50 with pre-calculated ratios}
+\usage{
+bootstrapIC50Precalculated(
+  dose,
+  response,
+  n_samples = 1000,
+  alpha = 0.1,
+  increasing = FALSE,
+  target_response = 0.5,
+  transform_x = TRUE
+)
+}
+\arguments{
+\item{dose}{Numeric vector of dose values.}
+
+\item{response}{Numeric vector of response values (pre-calculated ratios).}
+
+\item{n_samples}{Number of bootstrap samples. Default = 1000.}
+
+\item{alpha}{Significance level for confidence interval. Default = 0.10.}
+
+\item{increasing}{Logical. Fit non-decreasing if TRUE.}
+
+\item{target_response}{Numeric value for response level.}
+
+\item{transform_x}{Logical. If TRUE, applies log10(x + 1) transformation. Default = TRUE.}
+}
+\value{
+List with mean IC50, CI bounds, and transformed estimates.
+}
+\description{
+Bootstrap IC50 with pre-calculated ratios
+}
diff --git a/man/calculateTurnoverRatios.Rd b/man/calculateTurnoverRatios.Rd
new file mode 100644
index 0000000..67e47a3
--- /dev/null
+++ b/man/calculateTurnoverRatios.Rd
@@ -0,0 +1,100 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/protein_turnover_ratio_helper.R
+\name{calculateTurnoverRatios}
+\alias{calculateTurnoverRatios}
+\title{Calculate turnover ratios from MSstats FeatureLevelData}
+\usage{
+calculateTurnoverRatios(
+  feature_data,
+  channel_col = "LABEL",
+  heavy_label = "H",
+  light_label = "L",
+  time_col = "GROUP",
+  peptide_col = "PEPTIDE",
+  protein_col = "PROTEIN",
+  intensity_col = "INTENSITY",
+  run_col = "RUN",
+  peptide_selector = NULL,
+  agg_function = max,
+  normalize_tracer = FALSE,
+  tracer_constants = NULL
+)
+}
+\arguments{
+\item{feature_data}{Data frame from MSstats dataProcess()$FeatureLevelData}
+
+\item{channel_col}{Character. Name of column containing Heavy/Light labels. Default = "LABEL"}
+
+\item{heavy_label}{Character. Value in channel_col indicating heavy channel. Default = "H"}
+
+\item{light_label}{Character. Value in channel_col indicating light channel. Default = "L"}
+
+\item{time_col}{Character. Column containing timepoint information. Default = "GROUP"}
+
+\item{peptide_col}{Character. Column containing peptide sequences. Default = "PEPTIDE"}
+
+\item{protein_col}{Character. Column containing protein identifiers. Default = "PROTEIN"}
+
+\item{intensity_col}{Character. Column with intensity values. Default = "INTENSITY"}
+
+\item{run_col}{Character. Column identifying technical replicates. Default = "RUN"}
+
+\item{peptide_selector}{Optional function to select peptides.
+Function should take a data frame (grouped by protein) and return filtered data frame.
+Default = NULL (use all peptides).}
+
+\item{agg_function}{Function to aggregate duplicate peptide measurements (multiple charge states,
+transitions, etc.). Default = max (takes highest signal).}
+
+\item{normalize_tracer}{Logical. If TRUE, normalize by tracer incorporation. Default = FALSE}
+
+\item{tracer_constants}{Named numeric vector. Tracer constants for each timepoint.
+Required if normalize_tracer = TRUE}
+}
+\value{
+Data frame with columns: Protein, BaseSequence, TimeVal, Run, Heavy, Light, Total, H_frac, L_frac
+}
+\description{
+Calculate turnover ratios from MSstats FeatureLevelData
+}
+\examples{
+# Basic usage - all proteins, all peptides
+ratios <- calculateTurnoverRatios(
+  feature_data = quant_data$FeatureLevelData
+)
+
+# With peptide selection (top 10 by signal per protein)
+top10_selector <- function(df) {
+  top_peptides <- df \%>\%
+    group_by(BaseSequence) \%>\%
+    summarise(total_signal = sum(Intensity, na.rm = TRUE), .groups = "drop") \%>\%
+    arrange(desc(total_signal)) \%>\%
+    slice_head(n = 10) \%>\%
+    pull(BaseSequence)
+
+  df \%>\% filter(BaseSequence \%in\% top_peptides)
+}
+
+ratios_top10 <- calculateTurnoverRatios(
+  feature_data = quant_data$FeatureLevelData,
+  peptide_selector = top10_selector
+)
+
+# Use different aggregation function (e.g., median for robustness)
+ratios_median <- calculateTurnoverRatios(
+  feature_data = quant_data$FeatureLevelData,
+  agg_function = median
+)
+
+# With tracer normalization
+tracer_consts <- c("0hr" = 1.0, "1hr" = 0.95, "12hrs" = 0.85, "168hrs" = 0.75)
+ratios_norm <- calculateTurnoverRatios(
+  feature_data = quant_data$FeatureLevelData,
+  normalize_tracer = TRUE,
+  tracer_constants = tracer_consts
+)
+
+# Filter for specific protein after calculation
+protein_ratios <- ratios \%>\% filter(Protein == "A0A2K5TXF6")
+
+}
diff --git a/man/fitIsotonicRegression.Rd b/man/fitIsotonicRegression.Rd
index 6ac6996..594acfe 100644
--- a/man/fitIsotonicRegression.Rd
+++ b/man/fitIsotonicRegression.Rd
@@ -11,6 +11,7 @@ fitIsotonicRegression(
   increasing = FALSE,
   transform_x = TRUE,
   ratio_y = FALSE,
+  precalculated_ratios = FALSE,
   test_significance = FALSE
 )
 }
@@ -27,6 +28,9 @@ fitIsotonicRegression(
 
 \item{ratio_y}{Logical. If TRUE, converts log2 abundance to ratios relative to DMSO. Default is FALSE.}
 
+\item{precalculated_ratios}{Logical. If TRUE, assumes y contains pre-calculated ratios and skips
+internal ratio calculation. Default is FALSE.}
+
 \item{test_significance}{Logical. If TRUE, performs F-test to assess significance.}
 }
 \value{
diff --git a/man/parse_timepoint.Rd b/man/parse_timepoint.Rd
new file mode 100644
index 0000000..aa73b64
--- /dev/null
+++ b/man/parse_timepoint.Rd
@@ -0,0 +1,23 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/protein_turnover_ratio_helper.R
+\name{parse_timepoint}
+\alias{parse_timepoint}
+\title{Parse timepoint strings to numeric hours}
+\usage{
+parse_timepoint(time_strings)
+}
+\arguments{
+\item{time_strings}{Character vector of timepoint labels}
+}
+\value{
+Numeric vector of hours
+}
+\description{
+Parse timepoint strings to numeric hours
+}
+\examples{
+parse_timepoint(c("0hr", "1hr", "12hrs", "168hrs"))
+# Returns: 0 1 12 168
+
+}
+\keyword{internal}
diff --git a/man/plotIsotonic.Rd b/man/plotIsotonic.Rd
index e50d99f..31d62ae 100644
--- a/man/plotIsotonic.Rd
+++ b/man/plotIsotonic.Rd
@@ -7,44 +7,57 @@
 plotIsotonic(
   fit,
   ratio = TRUE,
+  precalculated_ratios = FALSE,
   show_ic50 = FALSE,
-  target_response = target_response,
+  target_response = 0.5,
   drug_name = NULL,
   protein_name = NULL,
-  x_lab = expression(Log[10] ~ "[drug (M)]"),
-  y_lab = "Log2 Intensity",
+  x_lab = NULL,
+  y_lab = NULL,
   title = NULL,
   ci = NULL,
   legend = FALSE,
   theme_style = "classic",
   original_label = FALSE,
-  transform_x = TRUE
+  transform_x = TRUE,
+  color_by = NULL,
+  original_data = NULL
 )
 }
 \arguments{
 \item{fit}{A model object returned by fitIsotonicRegression().}
 
-\item{ratio}{Logical. If TRUE, shows plot on the ratio scale relative to DMSO (i.e. 0-1 scale). Default is FALSE.}
+\item{ratio}{Logical. If TRUE, shows plot on the ratio scale relative to DMSO. Default is FALSE.}
+
+\item{precalculated_ratios}{Logical. If TRUE, response values are pre-calculated ratios. Default is FALSE.}
 
 \item{show_ic50}{Logical. If TRUE, adds vertical line and annotation for IC50.}
 
+\item{target_response}{Numeric. Target response level for IC50/EC50. Default = 0.5.}
+
 \item{drug_name}{Drug name for plotting data.}
 
 \item{protein_name}{Protein name for plot.}
 
-\item{x_lab}{Label for x-axis.}
+\item{x_lab}{Character. Label for x-axis. If NULL, uses default based on transform_x.}
 
-\item{y_lab}{Label for y-axis.}
+\item{y_lab}{Character. Label for y-axis. If NULL, uses default based on ratio/precalculated_ratios.}
 
-\item{title}{Title for the plot.}
+\item{title}{Character. Title for the plot. If NULL, auto-generates.}
 
-\item{ci}{Logical. Include IC50 95\% confidence interval bands if TRUE. Default is FALSE.}
+\item{ci}{List with ci_lower and ci_upper. Include IC50 confidence interval bands if provided. Default is NULL.}
 
 \item{legend}{Logical. Show legend if TRUE.}
 
 \item{theme_style}{ggplot2 theme name to apply (default = "classic").}
 
 \item{original_label}{Logical. If TRUE, replace x-axis tick labels with original dose labels.}
+
+\item{transform_x}{Logical. Whether doses were log-transformed. Default = TRUE.}
+
+\item{color_by}{Character. Column name to color points by. Default is NULL.}
+
+\item{original_data}{Data frame containing original data with grouping variables.}
 }
 \value{
 A ggplot object.
diff --git a/man/selectTopNPeptides.Rd b/man/selectTopNPeptides.Rd
new file mode 100644
index 0000000..6e3bbb4
--- /dev/null
+++ b/man/selectTopNPeptides.Rd
@@ -0,0 +1,27 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/protein_turnover_ratio_helper.R
+\name{selectTopNPeptides}
+\alias{selectTopNPeptides}
+\title{Example peptide selector: Top N by total signal}
+\usage{
+selectTopNPeptides(df, n = 10)
+}
+\arguments{
+\item{df}{Data frame with BaseSequence and Intensity columns (for one protein)}
+
+\item{n}{Number of top peptides to select}
+}
+\value{
+Filtered data frame
+}
+\description{
+Example peptide selector: Top N by total signal
+}
+\examples{
+# Use as peptide_selector
+ratios <- calculateTurnoverRatios(
+  feature_data = quant_data$FeatureLevelData,
+  peptide_selector = function(df) selectTopNPeptides(df, n = 10)
+)
+
+}
diff --git a/vignettes/Protein_turnover_vignettes.Rmd b/vignettes/Protein_turnover_vignettes.Rmd
new file mode 100644
index 0000000..72bcb9f
--- /dev/null
+++ b/vignettes/Protein_turnover_vignettes.Rmd
@@ -0,0 +1,109 @@
+---
+title: "MSstatsResponse - Protein Turnover Vignette "
+output: html_document
+date: "2026-03-11"
+---
+
+```{r setup, include=FALSE}
+knitr::opts_chunk$set(echo = TRUE)
+```
+
+Load pacakges --- 
+```{r}
+library(MSstatsResponse)
+library(tidyverse)
+
+```
+
+Load most recent data --- 
+```{r}
+quant_data_fraction_impute = readRDS("../../../ProteinTurnover/data/quant_data_fraction_imputed_new.rds")
+
+```
+
+manual reformatting (for now)
+```{r}
+# rename columns --
+df_test = quant_data_fraction_impute$FeatureLevelData %>%
+  mutate(
+    protein_raw   = PROTEIN,
+    PROTEIN       = sub("_leu6$", "", PROTEIN),
+   # timepoint     = as.numeric(str_extract(GROUP, "\\d+")),
+    Channel      = str_detect(protein_raw, "_leu6$")
+  ) %>%
+  filter(
+    is.finite(INTENSITY),
+    !is.na(GROUP)
+  )
+
+df_test$LABEL = ifelse(df_test$Channel == TRUE, "H", "L")
+
+
+tracer_constants <- c(
+  `0` = 1, `1` = 0.286, `4` = 0.294, `12` = 0.261,
+  `24` = 0.304, `48` = 0.266, `96` = 0.263, `168` = 0.313
+)
+
+```
+
+NEW FUNCTION --- 'calculateTurnoverRatios'
+```{r}
+test_ratio_calc = calculateTurnoverRatios(
+    df_test,
+    channel_col = "LABEL", #IsotopeLabel
+    heavy_label = "H",
+    light_label = "L",
+    time_col = "GROUP",
+    peptide_col = "PEPTIDE",
+    protein_col = "PROTEIN",
+    intensity_col = "INTENSITY",
+    run_col = "RUN",
+    peptide_selector = NULL,
+    normalize_tracer = TRUE,
+    tracer_constants = tracer_constants)
+
+```
+
+must rename columns --- 
+```{r}
+#test_ratio_calc$drug = "Synthesis"
+
+#test_ratio_calc = MSstatsPrepareDoseResponseFit(
+ # data = test_ratio_calc,
+ # protein_column = 'Protein',
+ # log_abundance_column = 'H_frac',
+ # dose_column = 'TimeVal',
+ # drug_column = 'drug')
+
+# these column names are needed for all MSstatsResponse functions
+test_ratio_calc$protein = test_ratio_calc$Protein
+test_ratio_calc$drug = "Synthesis"
+test_ratio_calc$dose = test_ratio_calc$TimeVal
+test_ratio_calc$response = test_ratio_calc$H_frac
+
+
+# add color_by option for visualization 
+test_ratio_calc %>% filter (protein == example_proteins[1]) # should have 71 rows 
+
+```
+
+
+```{r}
+# visualize response
+visualizeResponseProtein(
+  test_ratio_calc,
+  protein_name = example_proteins[1],
+  drug_name = "Synthesis",
+  precalculated_ratios = TRUE,
+  ratio_response = FALSE,
+  transform_dose = FALSE,
+  increasing = TRUE,
+  show_ic50 = TRUE,
+  add_ci = FALSE,
+  target_response = 0.5, 
+  y_lab = "relative abundance",x_lab = "time (hrs)",
+  color_by = 'BaseSequence')
+
+
+```
+