Add files via upload

POTATO_TOMATO.py

@@ -0,0 +1,19 @@
+
+from matplotlib.figure import Figure
+
+
+def plot_TOMATO(FD):
+    figure1 = Figure(figsize=(9, 7), dpi=100)
+    subplot1 = figure1.add_subplot(111)
+
+    F = FD[:, 0]
+    PD_nm = FD[:, 1]
+
+    subplot1.set_xlabel("Distance [$\\mu$m]")
+    subplot1.set_ylabel("Force [pN]")
+    subplot1.plot(PD_nm, F, color='gray')
+    subplot1.tick_params('both', direction='in')
+    subplot1.set_ylim([min(F), max(F)])
+    subplot1.set_xlim([min(PD_nm) - 10, max(PD_nm) + 10])
+
+    return figure1

POTATO_config.py

@@ -1,12 +1,14 @@
 """Copyright 2021 Helmholtz-Zentrum für Infektionsforschung GmbH"""
+
 """default values for each data type"""
 
 default_values_HF = {
     'Downsampling rate': '30',
     'Butterworth filter degree': '4',
     'Cut-off frequency': '0.005',
     'Force threshold, pN': '5',
-    'Z-score': '3',
+    'Z-score force': '3',
+    'Z-score distance': '3',
     'Step d': '10',
     'Moving median window size': '800',
     'STD difference threshold': '0.05',
@@ -18,34 +20,36 @@
     'Butterworth filter degree': '2',
     'Cut-off frequency': '0.5',
     'Force threshold, pN': '5',
-    'Z-score': '3',
+    'Z-score force': '3',
+    'Z-score distance': '3',
     'Step d': '3',
     'Moving median window size': '20',
     'STD difference threshold': '0.05',
     'Data frequency, Hz': '1000'
 }
 
 default_values_CSV = {
-    'Downsampling rate': '2',
-    'Butterworth filter degree': '4',
-    'Cut-off frequency': '0.005',
+    'Downsampling rate': '1',
+    'Butterworth filter degree': '1',
+    'Cut-off frequency': '0.01',
     'Force threshold, pN': '5',
-    'Z-score': '3',
+    'Z-score force': '2.5',
+    'Z-score distance': '3',
     'Step d': '10',
-    'Moving median window size': '250',
+    'Moving median window size': '120',
     'STD difference threshold': '0.05',
-    'Data frequency, Hz': '1000'
+    'Data frequency, Hz': '20'
 }
 
 default_values_FIT = {
     'Persistance-Length ds, nm': '40',
     'Persistance-Length ds, upper bound, nm': '80',
     'Persistance-Length ds, lower bound, nm': '12',
     'Persistance-Length ss, nm': '1',
-    'Contour-Length ds, nm': '1256',
+    'Contour-Length ds, nm': '830',
     'Contour-Length ss, nm': '0',
     'Stiffness ds, pN': '500',
-    'Stiffness ds, upper bound, pN': '600',
+    'Stiffness ds, upper bound, pN': '1500',
     'Stiffness ds, lower bound, pN': '400',
     'Stiffness ss, pN': '800',
     'Force offset, pN': '0',

POTATO_find_steps.py

@@ -40,7 +40,7 @@ def moving_median(input_data, column_number, window_size):
 
 
 # sorting the data based on a x times STD threshold (normal distibuted noise vs extreme values from steps)
-def cut_off(input_array, column_number, mm, std, n_of_std):
+def cut_off(input_array, column_number, mm, std, z_score):
     # sort values - inside STD region, above STD region and below STD region
     F_values_inside = []
     PD_values_inside = []
@@ -60,17 +60,18 @@ def cut_off(input_array, column_number, mm, std, n_of_std):
 
     i = 0
     for n in range(0, len(input_array), 1):
-        if input_array[n, column_number] > mm[int(i)] + n_of_std * std:
+        if input_array[n, column_number] > mm[int(i)] + z_score * std:
             F_dt_above.append(input_array[n, 2])
             F_values_above.append(input_array[n, 0])
             PD_values_above.append(input_array[n, 1])
             PD_dt_above.append(input_array[n, 3])
 
-        elif input_array[n, column_number] < mm[int(i)] - n_of_std * std:
+        elif input_array[n, column_number] < mm[int(i)] - z_score * std:
             F_dt_below.append(input_array[n, 2])
             F_values_below.append(input_array[n, 0])
             PD_values_below.append(input_array[n, 1])
             PD_dt_below.append(input_array[n, 3])
+
         else:
             F_dt_inside.append(input_array[n, 2])
             F_values_inside.append(input_array[n, 0])
@@ -88,7 +89,7 @@ def cut_off(input_array, column_number, mm, std, n_of_std):
 
 
 # searching for minima in the force derivation to identify unfolding events
-def find_steps_F(input_settings, filename_i, Force_Distance, der_arr):
+def find_steps_F(input_settings, filename_i, Force_Distance, der_arr, orientation):
     global y_vector_F
     global F_mm2_STD2_positive
     global F_mm2_STD2_negative
@@ -98,7 +99,7 @@ def find_steps_F(input_settings, filename_i, Force_Distance, der_arr):
 
     STD_1 = STD(der_arr, 2)
     F_mm = moving_median(der_arr, 2, input_settings['window_size'])
-    Above, Inside, Below, inside_indices_F = cut_off(der_arr, 2, F_mm, STD_1, input_settings['z-score'])
+    Above, Inside, Below, inside_indices_F = cut_off(der_arr, 2, F_mm, STD_1, input_settings['z-score_f'])
 
     F_mm2_STD2_positive = []
     F_mm2_STD2_negative = []
@@ -107,15 +108,15 @@ def find_steps_F(input_settings, filename_i, Force_Distance, der_arr):
     while abs(STD_1 - STD(Inside, 2)) / STD_1 > input_settings['STD_diff']:
         F_mm = moving_median(Inside, 2, input_settings['window_size'])
         STD_1 = STD(Inside, 2)
-        Above, Inside, Below, inside_indices_F = cut_off(der_arr, 2, F_mm, STD_1, input_settings['z-score'])
+        Above, Inside, Below, inside_indices_F = cut_off(der_arr, 2, F_mm, STD_1, input_settings['z-score_f'])
         n_runs = n_runs + 1
 
     if STD_1 < 0.05:
         STD_1 = 0.05
 
     print('STD is', STD_1)
 
-    Above, Inside, Below, inside_indices_F = cut_off(der_arr, 2, F_mm, STD_1, input_settings['z-score'])
+    Above, Inside, Below, inside_indices_F = cut_off(der_arr, 2, F_mm, STD_1, input_settings['z-score_f'])
     F_mm = moving_median(Inside, 2, input_settings['window_size'])
 
     y_vector_F = []
@@ -129,8 +130,8 @@ def find_steps_F(input_settings, filename_i, Force_Distance, der_arr):
             F_mm.insert(n, F_mm[last])
 
     for i in range(len(F_mm)):
-        F_mm2_STD2_positive.append(F_mm[i] + input_settings['z-score'] * STD_1)
-        F_mm2_STD2_negative.append(F_mm[i] - input_settings['z-score'] * STD_1)
+        F_mm2_STD2_positive.append(F_mm[i] + input_settings['z-score_f'] * STD_1)
+        F_mm2_STD2_negative.append(F_mm[i] - input_settings['z-score_f'] * STD_1)
 
     # find the step points
     # for those steps that cross the STD2 threshold -> find the closest 0 values prior/following to the crossing one
@@ -158,6 +159,7 @@ def find_steps_F(input_settings, filename_i, Force_Distance, der_arr):
         PD_start_F.append(der_arr[i_start, 1])
         dict1 = {
             "filename": filename_i,
+            "orientation": orientation,
             "Derivation of": 'Force',
             'step #': n_steps,
             'F1': der_arr[i_start, 0],
@@ -176,7 +178,7 @@ def find_steps_F(input_settings, filename_i, Force_Distance, der_arr):
 
 
 # searching for maxima in the distance derivation to identify unfolding events
-def find_steps_PD(input_settings, filename_i, Force_Distance, der_arr):
+def find_steps_PD(input_settings, filename_i, Force_Distance, der_arr, orientation):
     global y_vector_PD
     global PD_mm2_STD2_positive
     global PD_mm2_STD2_negative
@@ -187,7 +189,7 @@ def find_steps_PD(input_settings, filename_i, Force_Distance, der_arr):
     STD_1 = STD(der_arr, 3)
     PD_mm = moving_median(der_arr, 3, input_settings['window_size'])
 
-    Above, Inside, Below, inside_indices_PD = cut_off(der_arr, 3, PD_mm, STD_1, input_settings['z-score'])
+    Above, Inside, Below, inside_indices_PD = cut_off(der_arr, 3, PD_mm, STD_1, input_settings['z-score_d'])
 
     PD_mm2_STD2_positive = []
     PD_mm2_STD2_negative = []
@@ -196,15 +198,15 @@ def find_steps_PD(input_settings, filename_i, Force_Distance, der_arr):
     while abs(STD_1 - STD(Inside, 3)) / STD_1 > input_settings['STD_diff']:
         PD_mm = moving_median(Inside, 3, input_settings['window_size'])
         STD_1 = STD(Inside, 3)
-        Above, Inside, Below, inside_indices_PD = cut_off(der_arr, 3, PD_mm, STD_1, input_settings['z-score'])
+        Above, Inside, Below, inside_indices_PD = cut_off(der_arr, 3, PD_mm, STD_1, input_settings['z-score_d'])
         n_runs = n_runs + 1
 
     if STD_1 < 0.05:
         STD_1 = 0.05
 
     print('STD is', STD_1)
 
-    Above, Inside, Below, inside_indices_PD = cut_off(der_arr, 3, PD_mm, STD_1, input_settings['z-score'])
+    Above, Inside, Below, inside_indices_PD = cut_off(der_arr, 3, PD_mm, STD_1, input_settings['z-score_d'])
     PD_mm = moving_median(Inside, 3, input_settings['window_size'])
 
     y_vector_PD = []
@@ -218,13 +220,11 @@ def find_steps_PD(input_settings, filename_i, Force_Distance, der_arr):
             PD_mm.insert(n, PD_mm[last])
 
     for i in range(len(PD_mm)):
-        PD_mm2_STD2_positive.append(PD_mm[i] + input_settings['z-score'] * STD_1)
-        PD_mm2_STD2_negative.append(PD_mm[i] - input_settings['z-score'] * STD_1)
+        PD_mm2_STD2_positive.append(PD_mm[i] + input_settings['z-score_d'] * STD_1)
+        PD_mm2_STD2_negative.append(PD_mm[i] - input_settings['z-score_d'] * STD_1)
     # find the step points
     # for those steps that cross the 3*STD2 threshold -> find the closest 0 values prior/following to the crossing one
 
-    # for local minima
-
     loc_max = argrelextrema(Above[:, 3], np.greater)
 
     n_steps = 1
@@ -249,6 +249,7 @@ def find_steps_PD(input_settings, filename_i, Force_Distance, der_arr):
 
         dict1 = {
             "filename": filename_i,
+            "orientation": orientation,
             "Derivation of": 'Distance',
             'step #': n_steps,
             'F1': der_arr[i_start, 0],

POTATO_fitting.py

@@ -10,18 +10,33 @@
 """define the functions used for fitting"""
 
 
-def fitting_ds(filename_i, input_settings, export_data, input_fitting, i_start, Force_Distance, derivation_array, F_low):
+def fitting_ds(filename_i, input_settings, export_data, input_fitting, i_start, Force_Distance, derivation_array, F_low, TOMATO_param):
     global model_ds, fit_ds
     global ds_fit_dict
     global f_fitting_region_ds, d_fitting_region_ds
     global export_fit_ds
     global fitting_model
 
-    start_step1 = np.where(derivation_array[:, 1] == i_start)
-    start_step1 = start_step1[0][0]
+    if TOMATO_param == 0:
+        start_step1 = np.where(derivation_array[:, 1] == i_start)
+        start_step1 = start_step1[0][0]
 
-    f_fitting_region_ds = Force_Distance[0:start_step1 * input_settings['step_d'] + len(F_low), 0]
-    d_fitting_region_ds = Force_Distance[0:start_step1 * input_settings['step_d'] + len(F_low), 1]
+        f_fitting_region_ds = Force_Distance[0:start_step1 * input_settings['step_d'] + len(F_low), 0]
+        d_fitting_region_ds = Force_Distance[0:start_step1 * input_settings['step_d'] + len(F_low), 1]
+
+    elif TOMATO_param == 1:
+        i_start = Force_Distance[-1, 1]
+        f_fitting_region_ds = Force_Distance[0]
+        d_fitting_region_ds = Force_Distance[1]
+
+    delta_f = f_fitting_region_ds[-1] - f_fitting_region_ds[0]
+    delta_d = d_fitting_region_ds[-1] - d_fitting_region_ds[0]
+    d_f_ratio = delta_d / delta_f
+
+    # downsample the data used for fitting with a dD/dF ratio
+    while len(f_fitting_region_ds) > 100 * d_f_ratio:
+        f_fitting_region_ds = f_fitting_region_ds[::2]
+        d_fitting_region_ds = d_fitting_region_ds[::2]
 
     model_ds = lk.inverted_odijk("ds_part").subtract_independent_offset() + lk.force_offset("ds_part")
 
@@ -63,14 +78,16 @@ def fitting_ds(filename_i, input_settings, export_data, input_fitting, i_start,
     # calculate the integral until the first unfolding step
     # used to calculate the work done by the machine
     distance_integral = np.arange(min(Force_Distance[:, 1]), i_start)
-    ds_integral = model_ds(distance_integral, fit_ds)
+    ds_integral = model_ds(distance_integral, fit_ds.params)
     area_ds = simps(ds_integral)
     print("area_ds = " + str(area_ds))
 
     # export the fitting parameters
     ds_fit_dict = {
         'filename': filename_i,
         'model': 'WLC',
+        'model_ds': model_ds,
+        'fit_model': fit_ds,
         'log_likelihood': fit_qual,
         'Lc_ds': fit_ds["ds_part/Lc"].value,
         'Lp_ds': fit_ds["ds_part/Lp"].value,
@@ -174,12 +191,12 @@ def fitting_ss(filename_i, input_settings, export_data, input_fitting, i_start,
 
     # calculate the integrals of the fitted functions
     distance_integral_fit_start = np.arange(min(Force_Distance[:, 1]), i_start)
-    ss_integral_start = model_ss(distance_integral_fit_start, fit_ss)
+    ss_integral_start = model_ss(distance_integral_fit_start, fit_ss.params)
     area_ss_fit_start = simps(ss_integral_start)
     print("area_ss_start = " + str(area_ss_fit_start))
 
     distance_integral_fit_end = np.arange(min(Force_Distance[:, 1]), i_end)
-    ss_integral_end = model_ss(distance_integral_fit_end, fit_ss)
+    ss_integral_end = model_ss(distance_integral_fit_end, fit_ss.params)
     area_ss_fit_end = simps(ss_integral_end)
     print("area_ss_end = " + str(area_ss_fit_end))
 
@@ -210,7 +227,7 @@ def fitting_ss(filename_i, input_settings, export_data, input_fitting, i_start,
 
 def plot_fit(fit, start_force_ss, start_distance_ss, Force_Distance, save_folder, filename_i, start_time):
     distance = np.arange(min(Force_Distance[:, 1]), max(Force_Distance[:, 1]) + 50, 2)
-    F_ds_model = model_ds(distance, fit_ds)
+    F_ds_model = model_ds(distance, fit_ds.params)
 
     legend_elements = [
         Line2D([0], [0], color='k', lw=1, alpha=0.85),
@@ -226,7 +243,7 @@ def plot_fit(fit, start_force_ss, start_distance_ss, Force_Distance, save_folder
     plt.legend(legend_elements, ['FD-Curve', 'Part used for fitting', 'Fitted WLC model'])
 
     for i in range(0, len(fit)):
-        F_ss_model = model_ss(distance, fit[i])
+        F_ss_model = model_ss(distance, fit[i].params)
         plt.scatter(start_distance_ss[i], start_force_ss[i], s=4)
         plt.plot(distance, F_ss_model, linestyle='dashed', color='gray')
 

POTATO_preprocessing.py

@@ -13,50 +13,37 @@ def preprocess_RAW(Force, Distance, input_settings):
     b, a = signal.butter(input_settings['filter_degree'], input_settings['filter_cut_off'])
     filteredForce = signal.filtfilt(b, a, Force_ds)
     filteredDistance = signal.filtfilt(b, a, Distance_ds)
-    filteredDistance_ready = filteredDistance * 1000
-
-    Force_Distance = np.column_stack((filteredForce, filteredDistance_ready))
-
-    if Force_Distance[0, 1] > Force_Distance[-1, 1]:  # reverse
-        Force_Distance = np.flipud(Force_Distance)
 
+    Force_Distance = np.column_stack((filteredForce, filteredDistance * 1000))
     Force_Distance_um = np.column_stack((filteredForce, filteredDistance))
 
     return Force_Distance, Force_Distance_um
 
 
 # creates a dataset from min force threshold to max force value
 def trim_data(FD_data, F_min):
-    global F_trimmed, PD_trimmed, F_low
-
-    F_trimmed = []
-    PD_trimmed = []
+    F_trimmed = np.array([])
+    PD_trimmed = np.array([])
+    F_low = np.array([])
 
     F_max = np.where(FD_data[:, 0] == max(FD_data[:, 0]))
     fi = F_max[0][0]
 
-    while FD_data[fi, 0] < FD_data[fi - 10, 0]:
-        fi = fi - 1
-
-    while FD_data[fi, 1] < FD_data[fi - 10, 1]:
+    while FD_data[fi, 0] > F_min and fi > 0:
         fi = fi - 1
-    fi0 = fi
-
-    fi = F_max[0][0]
-    # print(fi)
 
-    while FD_data[fi, 0] > F_min:
-        fi = fi - 1
-        # print(Force_Distance[fi, 0])
-    F_trimmed = FD_data[fi:fi0, 0]
-    PD_trimmed = FD_data[fi:fi0, 1]
-    F_low = FD_data[:fi, 0]
+    if not fi == F_max[0][0]:
+        F_trimmed = FD_data[fi:F_max[0][0], 0]
+        PD_trimmed = FD_data[fi:F_max[0][0], 1]
+        F_low = FD_data[:fi, 0]
+    elif fi == 0:
+        print('Could not trim this curve, data below minimum force threshold!')
 
     return F_trimmed, PD_trimmed, F_low
 
 
 # creates derivations for Force and Distance of the trimmed datasets
-def create_derivation(input_settings, Frequency_value):
+def create_derivation(input_settings, Frequency_value, F_trimmed, PD_trimmed, F_low):
     d_time = 1 / Frequency_value * input_settings['downsample_value'] * input_settings['step_d']
 
     x = input_settings['step_d']