Merge branch 'main' of github.com:pinheiroGroup/Kinbiont.jl

Author: aivuk

Project.toml

@@ -1,7 +1,7 @@
 name = "Kinbiont"
 uuid = "3b7c519f-be5e-4159-a208-6486ce96bb37"
 authors = ["Fabrizio Angaroni", "Edgar Zanella Alvarenga"]
-version = "1.0.2"
+version = "1.0.10"
 
 [deps]
 AbstractTrees = "1520ce14-60c1-5f80-bbc7-55ef81b5835c"

README.md

@@ -1,10 +1,13 @@
-<p align="center">
-  <img src="https://kinbiont.fuzue.org/assets/logo.png">
-</p>
-
 
 Ecological and evolutionary processes of microbes are characterized by observables like growth rates and biomass yield, inferred from kinetics experiments. 
 Across conditions, these observables map response patterns such as antibiotic growth inhibition and yield dependence on substrate.
 But how do we extract ecological and evolutionary insights from massive datasets of time-resolved microbial data? Here we introduce Kinbiont — an ecosystem of numerical methods combining state-of-the-art solvers for ordinary differential equations, non-linear optimization, signal processing, and interpretable machine learning algorithms.
 Kinbiont provides a comprehensive, model-based analysis pipeline, covering all aspects of microbial kinetics data, from preprocessing to result interpretation. 
-For the stable documentation please consult [Documentation](https://kinbiont.fuzue.org/)
+
+
+
+For stable documentation please consult [Documentation](https://kinbiont.fuzue.org/). 
+
+Pre-print at  [Biorxiv](https://www.biorxiv.org/content/10.1101/2024.09.09.611847v1) .
+
+Data and scripts to reproduce the paper results at [Kinbiont utilities](https://github.com/pinheiroGroup/Kinbiont_utilities)

src/Fit_one_file_functions.jl

@@ -43,7 +43,7 @@ Fits a logarithmic-linear model to data from a .csv file. The function assumes t
 - `pt_avg=7`:  Number of points used in the rolling average smoothing.
 - `pt_smoothing_derivative=7`: Number of points for evaluating specific growth rate. If less than 2, uses interpolation; otherwise, a sliding window approach is used.
 - `pt_min_size_of_win=7`: Minimum size of the exponential windows in terms of the number of smoothed points.
-- `type_of_win="maximum"`: Method for selecting the exponential phase window. Options are `"maximum"` or `"global_thr"`.
+- `type_of_win="maximum"`: Method for selecting the exponential phase window. Options are "maximum"` or `"global_thr"` "max_with_min_OD".
 - `threshold_of_exp=0.9`: Threshold in quantile to define the exponential windows, between 0 and 1.
 - `do_blank_subtraction="avg_blank"`: Method for blank subtraction. Options include `"NO"`, `"avg_subtraction"`, and `"time_avg"`.
 - `blank_value=0.0`: Average value of the blank, used only if `do_blank_subtraction` is not `"NO"`.
@@ -143,8 +143,11 @@ function fit_one_file_Log_Lin(
 
     # excluding blank data and discarded wells
     if length(list_of_blank) > 0
+
         names_of_cols = filter!(e -> !(e in list_of_blank), names_of_cols)
+
     end
+    
     if length(list_of_discarded) > 0
 
         names_of_cols = filter!(e -> !(e in list_of_discarded), names_of_cols)
@@ -1662,7 +1665,7 @@ function segment_gr_analysis_file(
     if write_res == true
 
         CSV.write(
-            string(path_to_results,results, label_exp, "_results.csv"),
+            string(path_to_results,label_exp, "_results.csv"),
             Tables.table(Matrix(results)),
         )
 

src/Fit_one_well_functions.jl

@@ -34,7 +34,7 @@ Fits a logarithmic-linear model to data from a .csv file. The function assumes t
 - `pt_avg=7`:  Number of points used in the rolling average smoothing.
 - `pt_smoothing_derivative=7`: Number of points for evaluating specific growth rate. If less than 2, uses interpolation; otherwise, a sliding window approach is used.
 - `pt_min_size_of_win=7`: Minimum size of the exponential windows in terms of the number of smoothed points.
-- `type_of_win="maximum"`: Method for selecting the exponential phase window. Options are `"maximum"` or `"global_thr"`.
+- `type_of_win="maximum"`: Method for selecting the exponential phase window. Options are `"maximum"`, `"global_thr"` or `"max_with_min_OD"`.
 - `threshold_of_exp=0.9`: Threshold in quantile to define the exponential windows, between 0 and 1.
 - `do_blank_subtraction="avg_blank"`: Method for blank subtraction. Options include `"NO"`, `"avg_subtraction"`, and `"time_avg"`.
 - `blank_value=0.0`: Average value of the blank, used only if `do_blank_subtraction` is not `"NO"`.
@@ -85,7 +85,7 @@ function fitting_one_well_Log_Lin(
     method_multiple_scattering_correction="interpolation",
     calibration_OD_curve="NA",
     thr_lowess=0.05,
-    start_exp_win_thr=0.05
+    start_exp_win_thr=0.01,
     )
     if multiple_scattering_correction == true
 
@@ -120,6 +120,10 @@ function fitting_one_well_Log_Lin(
     t_start = 0.0
 
     if type_of_win == "maximum"
+
+
+
+
         for yy = 1:(index_of_max-2)
             if specific_gr[index_of_max-yy] >= lb_of_distib
                 t_start = copy(specific_gr_times[index_of_max-yy])
@@ -187,17 +191,75 @@ function fitting_one_well_Log_Lin(
 
     end
 
+
+    if type_of_win == "max_with_min_OD"
+
+        index_od_over_thr = findfirst(data_smooted[2, :] .> start_exp_win_thr)
+
+        if isnothing(index_od_over_thr) == false && length(specific_gr[index_od_over_thr:end]) > 2
+
+            lb_of_distib = quantile(specific_gr[index_od_over_thr:end], threshold_of_exp)
+
+            index_of_max = argmax(specific_gr[index_od_over_thr:end])[1] + index_od_over_thr - 1
+
+
+
+            for yy = 1:(index_of_max-2)
+                if specific_gr[index_of_max-yy] >= lb_of_distib
+                    t_start = copy(specific_gr_times[index_of_max-yy])
+                else
+                    if specific_gr[(index_of_max-yy-1)] < lb_of_distib
+                        break
+                    end
+
+                end
+            end
+
+            # index of start
+            index_of_t_start = findfirst(x -> x > t_start, data_smooted[1, :])[1]
+
+            if index_od_over_thr[1] > index_of_t_start
+                    index_of_t_start = index_od_over_thr[1]
+                   t_start = specific_gr_times[index_of_t_start]
+            end
+
+            # searching t_end of the exp phase
+            t_end = specific_gr_times[end]
+
+            for yy = index_of_max:(eachindex(specific_gr)[end]-1)
+                if specific_gr[yy] >= lb_of_distib
+                    t_end = copy(specific_gr_times[yy])
+                else
+                    if specific_gr[(yy+1)] < lb_of_distib
+                        break
+                    end
+                end
+            end
+            index_of_t_end = findfirst(x -> x > t_end, data_smooted[1, :])[1]
+
+
+        else
+            # the conditions are not satisfied i fix t_end and  t_start to be discarded later and have all results to missing
+            t_end = 100
+            t_start = 200
+            index_of_t_start = 1
+            index_of_t_end = index_of_t_start + 2 * pt_min_size_of_win
+
+        end
+    end
+
+
     # checking the minimum size of the window before fitting
     if (index_of_t_end - index_of_t_start) < pt_min_size_of_win
         index_of_t_start = convert(Int, index_of_max - floor(pt_min_size_of_win / 2))
         index_of_t_end = convert(Int, index_of_max + floor(pt_min_size_of_win / 2))
 
         if index_of_t_start < 1
-            index_of_t_start =1 
+            index_of_t_start =1
         end
 
         if index_of_t_end > length(data_smooted[1, :])
-            index_of_t_end = length(data_smooted[1, :]) 
+            index_of_t_end = length(data_smooted[1, :])
         end
     end
 
@@ -213,7 +275,7 @@ function fitting_one_well_Log_Lin(
 
         sigma_a = sigma_b = r = zeros(N)
         Theoretical_fitting = coeff_1 .+ data_to_fit_times .* coeff_2
-        
+
         Cantrell_errors = sqrt(sum((data_to_fit_values - coeff_2 * data_to_fit_times .- coeff_1) .^ 2) / (N - 2))  # goodness of fit
         sigma_b = sqrt(1 / sum((data_to_fit_times .- mean_x) .^ 2))
         sigma_a = Cantrell_errors * sqrt(1 / N + mean_x^2 * sigma_b^2)
@@ -246,6 +308,11 @@ function fitting_one_well_Log_Lin(
 
     else
 
+
+        data_to_fit_times  = missing
+        Theoretical_fitting = missing
+        confidence_band = missing
+
         results_lin_log_fit = [
             label_exp,
             name_well,
@@ -266,11 +333,10 @@ function fitting_one_well_Log_Lin(
     end
 
 
-    Kinbiont_res_one_well_log_lin = ("Log-lin", results_lin_log_fit, hcat(data_to_fit_times, data_to_fit_values), data_smooted, confidence_band)
+    Kinbiont_res_one_well_log_lin = ("Log-lin", results_lin_log_fit, hcat(data_to_fit_times, Theoretical_fitting), data_smooted, confidence_band)
 
     return Kinbiont_res_one_well_log_lin
-end
-
+  end
 
 """
     fitting_one_well_ODE_constrained(
@@ -746,7 +812,7 @@ function ODE_Model_selection(
 
         temp_start_param = param_array[mm]
 
-   
+
 
         # generating IC
         # setting initial conditions
@@ -1123,7 +1189,7 @@ end
     opt_params...
     )
 
-This function fits an Ordinary Differential Equation (ODE) model to segmented time-series data. Users provide fixed change points, and the function  models to each segment defined by these points. 
+This function fits an Ordinary Differential Equation (ODE) model to segmented time-series data. Users provide fixed change points, and the function  models to each segment defined by these points.
 
 # Arguments:
 
@@ -1455,7 +1521,7 @@ function segmentation_ODE(
     pt_avg=1, # number of the point to generate intial condition
     smoothing=true, # the smoothing is done or not?
     path_to_results="NA",
-    win_size=14, #  
+    win_size=14, #
     pt_smooth_derivative=7,
     beta_smoothing_ms=2.0,
     multiple_scattering_correction=false, # if true uses the given calibration curve to fix the data
@@ -1775,7 +1841,7 @@ function segment_gr_analysis(
     thr_lowess=0.05, # keyword argument of lowees smoothing
     type_of_detection="slinding_win",
     type_of_curve="original",
-    win_size=14, #  
+    win_size=14, #
     n_bins=40,
     method_peaks_detection="peaks_prominence"
 )
@@ -1870,4 +1936,4 @@ export ODE_Model_selection
 export one_well_morris_sensitivity
 export selection_ODE_fixed_intervals
 export segmentation_ODE
-export segment_gr_analysis
+export segment_gr_analysis

src/NL_fit_one_well.jl

@@ -358,7 +358,7 @@ function fit_NL_model_with_sensitivity(data::Matrix{Float64}, # dataset first ro
 
     best_res_param = fin_param[:, index_best+1]
 
-    best_fitted_model = model_function(best_res_param[3:(end-3)], data[1, :])
+    best_fitted_model = model_function(best_res_param[4:(end-3)], data[1, :])
 
     if write_res == true
         mkpath(path_to_results)
@@ -588,7 +588,7 @@ function fit_NL_model_bootstrap(data::Matrix{Float64}, # dataset first row times
     index_best = findmin(fin_param[end, 2:end])[2]
 
     best_res_param = fin_param[:, index_best+1]
-    best_fitted_model = model_function(best_res_param[3:(end-3)], data[1, :])
+    best_fitted_model = model_function(best_res_param[4:(end-3)], data[1, :])
 
     if write_res == true
         mkpath(path_to_results)
@@ -807,7 +807,7 @@ function NL_error_blanks(data::Matrix{Float64}, # dataset first row times second
     index_best = findmin(fin_param[end, 2:end])[2]
 
     best_res_param = fin_param[:, index_best+1]
-    best_fitted_model = model_function(best_res_param[3:(end-3)], data[1, :])
+    best_fitted_model = model_function(best_res_param[4:(end-3)], data[1, :])
 
     if write_res == true
         mkpath(path_to_results)