adding boxplot analysis + bugfix in unique track selection

README.md

@@ -73,7 +73,7 @@ To run track-analyzer
 
 
 ## Installation-free version
-The installation-free online version is available [here](https://galaxy.pasteur.fr/root?tool_id=toolshed.pasteur.fr/repos/rplanel/track_analyzer/track-analyzer/0.1.0). It is run on the web-base platform [Galaxy](https://galaxyproject.org/), which is easy to use (some documentation regarding Galaxy is available [here](https://training.galaxyproject.org/training-material/)). A quickstart tutorial to Galaxy's interface is presented in Track Analyzer [documentation](https://track-analyzer.pages.pasteur.fr/track-analyzer/). 
+The installation-free online version is available [here](https://galaxy.pasteur.fr/root?tool_id=toolshed.pasteur.fr/repos/rplanel/track_analyzer/track-analyzer/0.1.0). It is run on the web-base platform [Galaxy](https://galaxyproject.org/), which is easy to use (some documentation regarding Galaxy is available [here](https://training.galaxyproject.org/training-material/)). A quickstart tutorial to Galaxy's interface is presented in Track Analyzer's [documentation](https://track-analyzer.pages.pasteur.fr/track-analyzer/). 
 
 ## Documentation
 You can find a complete documentation [here](https://track-analyzer.pages.pasteur.fr/track-analyzer/).

track_analyzer/calculate.py

@@ -26,6 +26,8 @@
 import sys
 import os.path as osp
 
+from itertools import product
+
 import numpy as np
 import pandas as pd
 import scipy.interpolate as sci
@@ -169,36 +171,48 @@ def recompute_pos(df, dimensions=['x', 'y', 'z'], lengthscale=1.):
 def compute_track_prop(df, dimensions=['x', 'y', 'z'],custom_var={}):
     """Compute track properties by averaging all data along a traj"""
 
-    groups = df.groupby('track')
-
     # colmuns to average
     scaled_dimensions = [dim + '_scaled' for dim in dimensions]
     vel = ['v' + dim for dim in dimensions]
     acc = ['a' + dim for dim in dimensions]
-    mean_col = ['t'] + dimensions + scaled_dimensions + vel + acc + ['v', 'a'] + list(custom_var.keys())
+    col_to_avg = ['t'] + dimensions + scaled_dimensions + vel + acc + ['v', 'a'] + list(custom_var.keys())
     
     # all columns
-    columns = ['track', 'track_length'] + mean_col
+    columns = ['track', 'track_length'] + col_to_avg
     if 'subset' in df.columns:
         columns.append('subset')
-    df_out = pd.DataFrame(columns=columns)
+    df_mean = pd.DataFrame(columns=columns)
+    df_std = pd.DataFrame(columns=columns)
 
     #check custom_var are numeric
     for var in custom_var.keys():
         df[var] = pd.to_numeric(df[var],errors='coerce')
 
     # average
-    for i, track in enumerate(df['track'].unique()):
-        traj = groups.get_group(track)
+    groups, group_values = tpr.get_unique_track_groups(df)
+
+    for i, group_value in enumerate(group_values):
+        traj = groups.get_group(group_value)
         track_length = traj['t'].max() - traj['t'].min()
-        val = [track, track_length]
-        for col in mean_col:
-            val.append(traj[col].mean())
+        if type(group_value) is tuple:
+            subset = group_value[0] 
+            track = group_value[1] 
+        else: 
+            track = group_value
+
+        mean_val = {'track':track, 'track_length':track_length}
+        std_val = {'track':track, 'track_length':track_length}
+        for col in col_to_avg:
+            mean_val[col] = traj[col].mean()
+        for col in col_to_avg:
+            std_val[col] = traj[col].std()
         if 'subset' in df.columns:
-            val.append(traj['subset'].values[0])
-        df_out.loc[i, columns] = val
+            mean_val['subset'] = subset
+            std_val['subset'] = subset
+        df_mean.loc[i,:] = mean_val
+        df_std.loc[i,:] = std_val
 
-    return df_out
+    return df_mean, df_std
 
 
 def interpolate_field(data_dir, df, groups, grids, frame, field_values=['vx', 'vy', 'vz', 'ax', 'ay', 'az', 'v', 'a'],
@@ -594,10 +608,10 @@ def compute_all_Voronoi(data_dir, df, outdir=None, compute_local_area=True, area
     # calculate track average
     if df_mean is not None:
         df_mean['area'] = np.nan
-        groups = df.groupby('track')
-        for track in df['track'].unique():
-            traj = groups.get_group(track)
-            ind = df_mean['track'] == track
+        track_groups, track_groups_values = tpr.get_unique_track_groups(df)
+        for tgv in track_groups_values:
+            traj = tpr.get_traj(track_groups, tgv)
+            ind = ((df_mean['subset'] == tgv[0]) & (df_mean['track'] == tgv[1])) if type(tgv) is tuple else df_mean['track'] == tgv
             df_mean.loc[ind, 'area'] = traj['area'].mean()
 
     return data_dict

track_analyzer/plotting.py

@@ -212,29 +212,29 @@ def plot_traj(df, frame, data_dir, groups=None, image={'image_fn': None, 't_dim'
 
     #### plotting trajectories as lines
     if show_tail:
-        track_list = group['track'].values
-        track_groups = df.groupby(['track'])
-        for track in track_list:
-            traj = tpr.get_traj(track_groups, track, max_frame=frame)
+        track_groups, track_groups_values = tpr.get_unique_track_groups(df)
+        for tgv in track_groups_values:
+            track = tgv[1] if type(tgv) is tuple else tgv
+            traj = tpr.get_traj(track_groups, tgv, max_frame=frame)
             traj_length = traj.shape[0]
             X = traj['x'].values
             Y = traj['y'].values
 
-            ## color code
-            if color_code == "group":
-                if 'subset' in traj.columns:
-                    colors = color_list[traj['subset_order'].values[0] % len(color_list)]
+            if traj_length > 1:
+                ## color code
+                if color_code == "group":
+                    if 'subset' in traj.columns:
+                        colors = color_list[traj['subset_order'].values[0] % len(color_list)]
+                    else:
+                        colors = color_list[0]  # if color_coded by group but there's none, use only one color
+                elif color_code == "random":
+                    colors = color_list[track % len(color_list)]
+                elif color_code == "none":
+                    colors = color_list[0]
                 else:
-                    colors = color_list[0]  # if color_coded by group but there's none, use only one color
-            elif color_code == "random":
-                colors = color_list[track % len(color_list)]
-            elif color_code == "none":
-                colors = color_list[0]
-            else:
-                Val = traj[color_code].values
-                colors = tpr.get_cmap_color(Val, cmap, vmin=cmap_lim[0], vmax=cmap_lim[1])
+                    Val = traj[color_code].values
+                    colors = tpr.get_cmap_color(Val, cmap, vmin=cmap_lim[0], vmax=cmap_lim[1])
 
-            if traj_length > 1:
                 if color_code not in ['group','random','none']:
                     for j in range(1, traj_length):
                         ax.plot([X[j - 1], X[j]], [Y[j - 1], Y[j]], lw=lw, ls='-', color=colors[j])
@@ -563,7 +563,7 @@ def plot_MSD(data_dir, track, track_groups=None, df=None, df_out=None, fit_model
         df = data['df']
 
     if track_groups is None:
-        track_groups = df.groupby(['track'])
+        track_groups, track_groups_values = tpr.get_unique_track_groups(df)
     traj = tpr.get_traj(track_groups, track)
 
     if frame_subset is not None:
@@ -634,16 +634,17 @@ def plot_param_vs_param(data_dir, x_param, y_param, df=None, hue=None, hue_order
     if df is None:
         data = tpr.get_data(data_dir)
         df = data['df']
+    df_ = df.copy()  # copy to make edit df locally
 
     # define plotting directory
     if plot_dir is None:
         plot_dir = data_dir
 
     # make sure params are in df
-    if x_param not in df.columns: 
+    if x_param not in df_.columns: 
         print("Warning: parameter {} does not exist".format(x_param))
         return -1
-    if y_param not in df.columns: 
+    if y_param not in df_.columns: 
         print("Warning: parameter {} does not exist".format(y_param))
         return -1
 
@@ -660,23 +661,23 @@ def plot_param_vs_param(data_dir, x_param, y_param, df=None, hue=None, hue_order
 
     # make sure data is float and finite
     for p in [x_param,y_param]:
-        df[p] = pd.to_numeric(df[p],errors='coerce')
-        df = df[np.isfinite(df[p])]
+        df_[p] = pd.to_numeric(df_[p],errors='coerce')
+        df_ = df_[np.isfinite(df_[p])]
 
     # make sure that sns.lmplot does not use the continuous colormap
     if hue is not None:
-        df[hue] = df[hue].astype('category')  
+        df_[hue] = df_[hue].astype('category')  
 
     # make x_bins
     if x_bin_num is None:
         x_bins = None 
     else: 
         # make evenly spaced bins along x_param and exclude min and max
-        x_bins = np.linspace(df[x_param].min(),df[x_param].max(),x_bin_num+2)  # +2 so len(x_bins)=x_bin_num after removing min and max
+        x_bins = np.linspace(df_[x_param].min(),df_[x_param].max(),x_bin_num+2)  # +2 so len(x_bins)=x_bin_num after removing min and max
         x_bins = x_bins[1:-1]
 
     # plot
-    g = sns.lmplot(x=x_param, y=y_param, hue=hue, hue_order=hue_order, data=df, ci=ci, fit_reg=fit_reg, x_bins=x_bins, scatter=scatter,
+    g = sns.lmplot(x=x_param, y=y_param, hue=hue, hue_order=hue_order, data=df_, ci=ci, fit_reg=fit_reg, x_bins=x_bins, scatter=scatter,
         facet_kws={'despine':plot_config['despine']})
     if hue is not None:
         sns.move_legend(g, "right", frameon=False, title=None, bbox_to_anchor=(1.05, 0.5))
@@ -690,8 +691,8 @@ def plot_param_vs_param(data_dir, x_param, y_param, df=None, hue=None, hue_order
         ax.set_xlim(set_axis_lim[0], set_axis_lim[1])
         ax.set_ylim(set_axis_lim[2], set_axis_lim[3])
     else:  # recalculate because sometimes matplotlib auto limit fails
-        xlim_ = [df[x_param].min(), df[x_param].max()]
-        ylim_ = [df[y_param].min(), df[y_param].max()]
+        xlim_ = [df_[x_param].min(), df_[x_param].max()]
+        ylim_ = [df_[y_param].min(), df_[y_param].max()]
         xlim_ = [xlim_[0] - 0.05 * (xlim_[1] - xlim_[0]), xlim_[1] + 0.05 * (xlim_[1] - xlim_[0])]
         ylim_ = [ylim_[0] - 0.05 * (ylim_[1] - ylim_[0]), ylim_[1] + 0.05 * (ylim_[1] - ylim_[0])]
         ax.set_xlim(xlim_[0], xlim_[1])
@@ -707,7 +708,7 @@ def plot_param_vs_param(data_dir, x_param, y_param, df=None, hue=None, hue_order
     if export_data_pts:
         cols = [x_param, y_param, hue] if hue is not None else [x_param, y_param]
         fn = osp.join(plot_dir, prefix + '{}_vs_{}{}{}'.format(y_param, x_param, suffix, '.csv'))
-        df[cols].to_csv(fn)
+        df_[cols].to_csv(fn)
 
 
 def plot_param_hist(data_dir, param, df=None, hue=None, hue_order=None, hist=True, kde=True,
@@ -723,13 +724,14 @@ def plot_param_hist(data_dir, param, df=None, hue=None, hue_order=None, hist=Tru
     if df is None:
         data = tpr.get_data(data_dir)
         df = data['df']
+    df_ = df.copy()  # copy to make edit df locally
 
     # define plotting directory
     if plot_dir is None:
         plot_dir = data_dir
 
     # make sure param is in df
-    if param not in df.columns: 
+    if param not in df_.columns: 
         print("Warning: parameter {} does not exist".format(param))
         return -1
 
@@ -741,11 +743,11 @@ def plot_param_hist(data_dir, param, df=None, hue=None, hue_order=None, hist=Tru
         param_label = tpr.make_param_label(param, l_unit=info['length_unit'], t_unit=info['time_unit'])
 
     # make sure data is float and finite
-    df[param] = pd.to_numeric(df[param],errors='coerce')
-    df = df[np.isfinite(df[param])]
+    df_[param] = pd.to_numeric(df_[param],errors='coerce')
+    df_ = df_[np.isfinite(df_[param])]
 
     kind = "hist" if hist else "kde"
-    g = sns.displot(data=df, x=param, hue=hue, kind=kind, kde=kde,facet_kws={'despine':plot_config['despine']})
+    g = sns.displot(data=df_, x=param, hue=hue, kind=kind, kde=kde,facet_kws={'despine':plot_config['despine']})
     if hue is not None:
         sns.move_legend(g, "right", frameon=False, title=None, bbox_to_anchor=(1.05, 0.5))
     fig = g.figure
@@ -760,12 +762,12 @@ def plot_param_hist(data_dir, param, df=None, hue=None, hue_order=None, hist=Tru
     if export_data_pts:
         cols = [param, hue] if hue is not None else param
         fn = osp.join(plot_dir, prefix + '{}_hist{}{}'.format(param, suffix, '.csv'))
-        df[cols].to_csv(fn)
+        df_[cols].to_csv(fn)
 
 
 def plot_param_boxplot(data_dir, df, x_param, param, order=None, hue=None, save_stat=False, hue_order=None,
                        boxplot=True, swarmplot=True, plot_config=None, plot_dir=None, prefix='', suffix='',
-                       leg_lab=None):
+                       leg_lab=None, custom_var={}):
     """Plot boxplot between categories (given by x_param). Sub-categories can be plotted too (given by hue). A ttest can be performed has well"""
 
     plot_config = make_plot_config() if plot_config is None else plot_config
@@ -783,7 +785,9 @@ def plot_param_boxplot(data_dir, df, x_param, param, order=None, hue=None, save_
         plot_dir = data_dir
 
     # make sure params are in df
-    if x_param not in df.columns: 
+    if x_param is None:
+        pass
+    elif x_param not in df.columns: 
         print("Warning: parameter {} does not exist".format(x_param))
         return -1
     if param not in df.columns: 
@@ -791,7 +795,10 @@ def plot_param_boxplot(data_dir, df, x_param, param, order=None, hue=None, save_
         return -1
 
     info = tpr.get_info(data_dir)
-    param_label = tpr.make_param_label(param, l_unit=info['length_unit'], t_unit=info['time_unit'])
+    if param in custom_var.keys():
+        param_label = tpr.make_param_label(None, manual_symbol=custom_var[param]['name'], manual_unit=custom_var[param]['unit'])
+    else: 
+        param_label = tpr.make_param_label(param, l_unit=info['length_unit'], t_unit=info['time_unit'])
 
     # plotting
     fig, ax = plt.subplots(1, 1, figsize=figsize)
@@ -826,7 +833,7 @@ def plot_param_boxplot(data_dir, df, x_param, param, order=None, hue=None, save_
     if plot_config['despine']:
         sns.despine(fig)
 
-    filename_basis = osp.join(plot_dir, prefix + '{}_box{}'.format(param, suffix))
+    filename_basis = osp.join(plot_dir, prefix + '{}_boxplot{}'.format(param, suffix))
     fig.savefig(filename_basis + plot_config['format'], dpi=plot_config['dpi'], bbox_inches='tight')
     plt.close(fig)
 
@@ -868,7 +875,10 @@ def plot_param_boxplot(data_dir, df, x_param, param, order=None, hue=None, save_
         df_pval.to_csv(filename_basis + '_pvalue.csv')
 
     if export_data_pts:
-        cols = [x_param, param, hue] if hue is not None else [x_param, param]
+        cols = [x_param, param, hue]
+        for p in [x_param, hue]: 
+            if p is None: 
+                cols.remove(p)
         df[cols].to_csv(filename_basis + '.csv')
 
 
@@ -1246,7 +1256,7 @@ def plot_all_MSD(data_dir, df=None, df_out=None, fit_model="biased_diff", msd_al
 
     if df_out is None and fit_model is not None:  # compute track properties if MSD are to be fitted
         dimensions = ['x', 'y', 'z'] if 'z' in df.columns else ['x', 'y']
-        df_out = tca.compute_track_prop(df, dimensions)
+        df_out,df_std = tca.compute_track_prop(df, dimensions)
 
     if refresh:  # if refresh erase msd_all
         msd_all = None
@@ -1257,20 +1267,19 @@ def plot_all_MSD(data_dir, df=None, df_out=None, fit_model="biased_diff", msd_al
 
     # compute MSD
     if refresh:
-        track_groups = df.groupby(['track'])
-        track_list = df['track'].unique()
+        track_groups, track_groups_values = tpr.get_unique_track_groups(df)
         msd_list = []
-        for track in track_list:
-            msd = plot_MSD(data_dir, track=track, track_groups=track_groups, df=df, df_out=df_out, fit_model=fit_model,
+        for tgv in track_groups_values:
+            msd = plot_MSD(data_dir, track=tgv, track_groups=track_groups, df=df, df_out=df_out, fit_model=fit_model,
                            dim=dim, save_plot=plot_single_MSD, fitrange=fitrange, plot_dir=plot_dir,
                            plot_config=plot_config, logx=logx, logy=logy)
-            msd['track'] = track
+            
+            if type(tgv) is tuple: # if both subset and track are needed to identify tracks
+                msd['subset'] = tgv[0]
+                msd['track'] = tgv[1]
+            else: 
+                msd['track'] = tgv
 
-            # get subset if exists
-            if 'subset' in df.columns:
-                t = track_groups.get_group(track)
-                subset = t['subset'].values[0]
-                msd['subset'] = subset
             msd_list.append(msd)
         #concatenate
         msd_all = pd.concat(msd_list)
@@ -1400,7 +1409,7 @@ def plot_total_traj(data_dir, df, dim=3, plot_dir=None, plot_fn=None, plot_confi
     tpr.write_dict(specific_config, fn)
 
     # group by tracks
-    track_groups = df.groupby(['track'])
+    track_groups, track_groups_values = tpr.get_unique_track_groups(df)
 
     # color code  
     cmap_lim,color_code = tpr.get_cmap_lim(df,color_code,cmap_lim,dim)
@@ -1428,8 +1437,9 @@ def plot_total_traj(data_dir, df, dim=3, plot_dir=None, plot_fn=None, plot_confi
 
     label_list = []  # to store label that have already been used 
     fig, ax = plt.subplots(1, 1, figsize=plot_config['figsize'])
-    for i, track in enumerate(df['track'].unique()):
-        traj = tpr.get_traj(track_groups, track)
+    for tgv in track_groups_values:
+        track = tgv[1] if type(tgv) is tuple else tgv
+        traj = tpr.get_traj(track_groups, tgv)
         traj_length = traj.shape[0]
         first_frame = traj['frame'].min()
         x0, y0 = traj[traj['frame'] == first_frame][['x_scaled', 'y_scaled']].values[0]

track_analyzer/prepare.py

@@ -785,10 +785,10 @@ def filter_by_traj_len(df, min_traj_len=1, max_traj_len=None):
         max_traj_len = df['frame'].max() - df['frame'].min() + 1
     min_traj_len = 1 if min_traj_len is None else min_traj_len  # assign 1, if not given
 
-    tracks = df.groupby('track')
+    groups, group_values = get_unique_track_groups(df)
     df_list = []
-    for t in df['track'].unique():
-        track = tracks.get_group(t)
+    for group_value in group_values:
+        track = groups.get_group(group_value)
         if (track.shape[0] >= min_traj_len) & (track.shape[0] <= max_traj_len):
             df_list.append(track)
 
@@ -1406,10 +1406,20 @@ def make_traj_config(data_dir=None, export_config=True):
 
     hist_config = {'run': True,  # run plot_param_hist
                    'var_list': ['v'],  # list of variables to be plotted in histogram
-                   'mean_var_list': ['v']
-                   # list of variables averaged along the whole track to be plotted in histogram
+                   'mean_var_list': ['v']  # list of variables averaged along the whole track to be plotted in histogram
                    }
 
+    boxplot_config = {'run': True,  # run plot_param_hist
+                       'var_list': ['v'],  # list of variables to be plotted 
+                       'mean_var_list': ['v'],  # list of variables averaged along the whole track to be plotted
+                       'save_stat': True ,  # run statistical test (scipy.stats.ttest_ind(equal_var=False))
+                       'boxplot': True,  # show boxplot for variables in var_list
+                       'swarmplot': True,  # show swarmplot for variables in var_list
+                       'mean_boxplot': True,  # show boxplot for variables in mean_var_list
+                       'mean_swarmplot': True,  # show swarmplot for variables in mean_var_list
+                       'mean_save_stat': True ,  # run statistical test (scipy.stats.ttest_ind(equal_var=False)) for wholetrack parameters
+                       }
+
     total_traj_config = {'run': True,  # run plot_centered_traj
                          'hide_labels': False,  # hide trajectory ID
                          'label_size': 5,  # label size in points if hide_labels is False
@@ -1445,6 +1455,7 @@ def make_traj_config(data_dir=None, export_config=True):
               'hist_config': hist_config,
               'total_traj_config': total_traj_config,
               'voronoi_config': voronoi,
+              'boxplot_config': boxplot_config,
               }
 
     if export_config:
@@ -1547,3 +1558,19 @@ def make_all_config(data_dir=None, export_config=True):
                     "plot_config": plot_config,
                     }
     return all_config
+
+
+def get_unique_track_groups(df):
+    """
+    Split a df in track groups making sure they are unique tracks (splitting by subsets too if exist)
+    Return pandas groups and the values to get them
+    """
+    if 'subset' in df.columns: 
+        groups = df.groupby(by=['subset','track'])
+        group_values = zip(df['subset'].tolist(), df['track'].tolist())
+    else: 
+        groups = df.groupby(by='track')
+        group_values = df['track'].unique()
+
+    return groups, group_values
+

track_analyzer/scripts/analyze_maps.py

@@ -75,7 +75,7 @@ def map_analysis(data_dir, data=None, image=None, refresh=False, parallelize=Fal
     map_config = map_config_default if map_config is None else map_config
 
     # check that all configs are in map_config, if not load default
-    for key in ["grid_param", "map_param", "scalar_fields", "vector_fields", "vector_mean"]:
+    for key in map_config_default.keys():
         if key not in map_config.keys():
             map_config[key] = map_config_default[key]
         # add default param if missing

track_analyzer/scripts/analyze_tracks.py

@@ -75,15 +75,15 @@ def traj_analysis(data_dir, data=None, image=None, refresh=False, parallelize=Fa
     traj_config = traj_config_default if traj_config is None else traj_config
 
     # check that all configs are in traj_confign, if not load default
-    for key in ["traj_config_", "MSD_config", "scatter_config", "hist_config", "total_traj_config", "voronoi_config"]:
-        traj_config_default_key = traj_config_default[key]
-        # if don't exist, add
+    for key in traj_config_default.keys():
+        traj_config_default_ = traj_config_default[key]
+        # if doesn't exist, add
         if key not in traj_config.keys():
-            traj_config[key] = traj_config_default_key
+            traj_config[key] = traj_config_default_
         # check all defaut param are present
-        for k in traj_config_default_key.keys():
+        for k in traj_config_default_.keys():
             if k not in traj_config[key].keys():
-                traj_config[key][k] = traj_config_default_key[k] 
+                traj_config[key][k] = traj_config_default_[k] 
 
     traj_config_ = traj_config["traj_config_"]
     MSD_config = traj_config["MSD_config"]
@@ -91,6 +91,7 @@ def traj_analysis(data_dir, data=None, image=None, refresh=False, parallelize=Fa
     hist_config = traj_config["hist_config"]
     total_traj_config = traj_config["total_traj_config"]
     voronoi_config = traj_config["voronoi_config"]
+    boxplot_config = traj_config["boxplot_config"]
 
     # Filter data
     filters = tpr.init_filters(data_dir=data_dir, export_config=False) if filters is None else filters
@@ -147,9 +148,11 @@ def traj_analysis(data_dir, data=None, image=None, refresh=False, parallelize=Fa
         tpr.write_dict(filters_[i], fn)
 
         # compute mean track properties
-        mean_fn = osp.join(sub_dir, 'track_prop.csv')
-        df_prop = tca.compute_track_prop(df, dimensions,custom_var)
-        df_prop.to_csv(mean_fn)
+        mean_fn = osp.join(sub_dir, 'mean_track_data.csv')
+        std_fn = osp.join(sub_dir, 'std_track_data.csv')
+        df_mean, df_std = tca.compute_track_prop(df, dimensions,custom_var)
+        df_mean.to_csv(mean_fn)
+        df_std.to_csv(std_fn)
 
         # analysis and plotting
         if subset_analysis == 'together':
@@ -178,19 +181,19 @@ def traj_analysis(data_dir, data=None, image=None, refresh=False, parallelize=Fa
         if MSD_config['run']:
             if quiet < 2:
                 print("MSD analysis...")
-            df_prop = tpl.plot_all_MSD(data_dir, df, df_out=df_prop, fit_model=MSD_config['MSD_model'],
+            df_mean = tpl.plot_all_MSD(data_dir, df, df_out=df_mean, fit_model=MSD_config['MSD_model'],
                                        MSD_parameters=MSD_config, plot_config=plot_config, plot_dir=sub_dir, hue=hue,
                                        hue_order=hue_order)
-            df_prop.to_csv(mean_fn)
+            df_mean.to_csv(mean_fn)
 
         # Voronoi analysis
         if voronoi_config['run']:
             if quiet < 2:
                 print("Voronoi analysis...")
             vor_data = tca.compute_all_Voronoi(data_dir, df, outdir=sub_dir,compute_local_area=voronoi_config['compute_local_area'],
-                                            area_threshold=voronoi_config['area_threshold'], df_mean=df_prop)
+                                            area_threshold=voronoi_config['area_threshold'], df_mean=df_mean)
             # update csv files with area data
-            df_prop.to_csv(mean_fn)
+            df_mean.to_csv(mean_fn)
             df.to_csv(csv_fn) 
             
             if voronoi_config['plot']:
@@ -211,7 +214,7 @@ def traj_analysis(data_dir, data=None, image=None, refresh=False, parallelize=Fa
             if 'var_list' in hist_config.keys():
                 if len(hist_config['var_list']) > 0:
                     if quiet < 2:
-                        print("Plotting parameters histograms...")
+                        print("Plotting instantaneous parameters histograms...")
                 for p in hist_config['var_list']:
                     tpl.plot_param_hist(data_dir, p, df, plot_config=plot_config, plot_dir=sub_dir, hue=hue,
                                         hue_order=hue_order,custom_var=custom_var)
@@ -219,9 +222,9 @@ def traj_analysis(data_dir, data=None, image=None, refresh=False, parallelize=Fa
             if 'mean_var_list' in hist_config.keys():
                 if len(hist_config['mean_var_list']) > 0:
                     if quiet < 2:
-                        print("Plotting whole-track histograms...")
+                        print("Plotting whole-track parameters histograms...")
                 for p in hist_config['mean_var_list']:
-                    tpl.plot_param_hist(data_dir, p, df_prop, plot_config=plot_config, plot_dir=sub_dir, prefix='track_',
+                    tpl.plot_param_hist(data_dir, p, df_mean, plot_config=plot_config, plot_dir=sub_dir, prefix='track_',
                                         hue=hue, hue_order=hue_order,custom_var=custom_var)
 
         if scatter_config['run']:
@@ -231,7 +234,7 @@ def traj_analysis(data_dir, data=None, image=None, refresh=False, parallelize=Fa
             if 'couple_list' in scatter_config.keys():
                 if len(scatter_config['couple_list']) > 0:
                     if quiet < 2:
-                        print("Plotting couples of parameters...")
+                        print("Plotting scatter plots of instantaneous parameters...")
                 for param_vs_param in scatter_config['couple_list']:
                     x_param, y_param = param_vs_param
                     tpl.plot_param_vs_param(data_dir, x_param, y_param, df, plot_dir=sub_dir, plot_config=plot_config,
@@ -241,13 +244,34 @@ def traj_analysis(data_dir, data=None, image=None, refresh=False, parallelize=Fa
             if 'mean_couple_list' in scatter_config.keys():
                 if len(scatter_config['mean_couple_list']) > 0:
                     if quiet < 2:
-                        print("Plotting couples of whole-track parameters...")
+                        print("Plotting scatter plots of whole-track parameters...")
                 for param_vs_param in scatter_config['mean_couple_list']:
                     x_param, y_param = param_vs_param
-                    tpl.plot_param_vs_param(data_dir, x_param, y_param, df_prop, plot_dir=sub_dir, plot_config=plot_config,
+                    tpl.plot_param_vs_param(data_dir, x_param, y_param, df_mean, plot_dir=sub_dir, plot_config=plot_config,
                                             prefix='track_', hue=hue, hue_order=hue_order, x_bin_num=scatter_config["x_bin_num"], ci=scatter_config["ci"], 
                                             fit_reg=scatter_config["fit_reg"], scatter=scatter_config["scatter"],custom_var=custom_var)
 
+        if boxplot_config['run']:
+            fn = osp.join(config_dir, 'boxplot_config.csv')
+            tpr.write_dict(boxplot_config, fn)
+
+            if 'var_list' in boxplot_config.keys():
+                if len(boxplot_config['var_list']) > 0:
+                    if quiet < 2:
+                        print("Plotting instantaneous parameters boxplots...")
+                for p in boxplot_config['var_list']:
+                    tpl.plot_param_boxplot(data_dir, df=df, x_param=hue, param=p, order=hue_order,save_stat=boxplot_config['save_stat'], boxplot=boxplot_config['boxplot'],
+                                   swarmplot=boxplot_config['swarmplot'], plot_config=plot_config, plot_dir=sub_dir,custom_var=custom_var)
+
+            if 'mean_var_list' in boxplot_config.keys():
+                if len(boxplot_config['mean_var_list']) > 0:
+                    if quiet < 2:
+                        print("Plotting whole-track boxplots...")
+                for p in boxplot_config['mean_var_list']:
+                    tpl.plot_param_boxplot(data_dir, df=df_mean, x_param=hue, param=p, order=hue_order, save_stat=boxplot_config['mean_save_stat'], boxplot=boxplot_config['mean_boxplot'],
+                                   swarmplot=boxplot_config['mean_swarmplot'], plot_config=plot_config, plot_dir=sub_dir, prefix='track_',custom_var=custom_var)
+
+
     return df_list