Suggested edits to address reviews of Chris' transient growth rate branch

nuc_morph_analysis/analyses/evaluate_filters_and_outliers/evaluate_time_binning.py

@@ -2,6 +2,3 @@
 from nuc_morph_analysis.lib.preprocessing import add_times
 fig,ax = add_times.validate_dig_time_with_plot()
 fig.show()
-#%%
-fig,ax = add_times.validate_dig_time_with_plot(old_method=True)
-fig.show()

nuc_morph_analysis/analyses/neighbor_of_X/misc_neighbor_helper_functions.py

@@ -59,11 +59,6 @@ def get_a_cells_neighbors_as_track_id_list(df0,main_track_id,TIMEPOINT,return_se
     # find all neighbors that are immediate neighbors (and have passed the transition point)
     immediate_neighbor_track_ids = find_immediate_neighbors(df_after_transition_only,main_track_id,TIMEPOINT)
 
-    # # compute distances between nuclei and sort them by distance to main track
-    # dftime = dfcolony[dfcolony["index_sequence"] == TIMEPOINT]
-    # dist = compute_distances_between_nuclei(dftime,MAIN_TRACK_ID)
-    # sorted_dist,sorted_cell_ids,sorted_track_ids,sorted_index = get_ordered_list_of_nuclei_by_distance_to_track(dist,dftime)
-
     # combine the main track with the immediate neighbors into a list
     track_id_list = list(immediate_neighbor_track_ids)
     if return_self:

nuc_morph_analysis/analyses/volume/FigS10_workflow.py

@@ -4,14 +4,14 @@
 import matplotlib.pyplot as plt
 from pathlib import Path
 from nuc_morph_analysis.lib.preprocessing.global_dataset_filtering import load_dataset_with_features
-from nuc_morph_analysis.analyses.volume import plot_help
 
 from nuc_morph_analysis.analyses.volume.plot_help import (
     plot_neighbors_volume_over_time, plot_tracks_aligned_at_volume_drop_onset,update_plotting_params,
     adjust_axis_positions, plot_track_with_fit_line, plot_track_with_volume_dip,
     plot_dip_detection_validation, plot_dxdt_over_time_by_cell_cycle
 )
-from nuc_morph_analysis.lib.preprocessing import filter_data, compute_change_over_time, add_times
+from nuc_morph_analysis.lib.preprocessing import filter_data, add_times
+from nuc_morph_analysis.lib.visualization.matplotlib_to_axlist import type_axlist
 from nuc_morph_analysis.lib.visualization.notebook_tools import save_and_show_plot 
 from nuc_morph_analysis.analyses.volume import filter_out_dips
 from nuc_morph_analysis.analyses.neighbor_of_X.misc_neighbor_helper_functions import get_a_cells_neighbors_as_track_id_list
@@ -20,6 +20,8 @@
 from nuc_morph_analysis.lib.visualization import plotting_tools
 from nuc_morph_analysis.analyses.volume_variation import plot_features
 
+from nuc_morph_analysis.lib.visualization.example_tracks import EXAMPLE_TRACKS
+
 #%%
 # now load data with growth outliers
 df_outliers = load_dataset_with_features('all_baseline', remove_growth_outliers=False)
@@ -36,56 +38,45 @@
 # %%
 # S10 panel A, left
 # choose one track and its neighbors (at a given time) to plot over time
-MAIN_TRACK_ID = 75725
+MAIN_TRACK_ID = EXAMPLE_TRACKS['volume_dip_example']
 TIMEPOINT = 239
 track_id_list = get_a_cells_neighbors_as_track_id_list(df_outliers,MAIN_TRACK_ID,TIMEPOINT)
-fig,axlist = plot_neighbors_volume_over_time(df_outliers,track_id_list)
+fig,_ = plot_neighbors_volume_over_time(df_outliers,track_id_list)
 
-save_name = f"S10_A_left-immediate_neighbors_of_main_track_{MAIN_TRACK_ID}"
-save_path = str(save_dir / save_name)
+save_name: str = f"S10_A_left-immediate_neighbors_of_main_track_{MAIN_TRACK_ID}"
+save_path: Path = save_dir / save_name
 for ext in ['.png','.pdf']:
     save_and_show_plot(save_path,ext,fig,transparent=False,keep_open=True)
 plt.show()
 
 #%%
 # S10 panel A, middle
-fig,axlist = plot_tracks_aligned_at_volume_drop_onset(df_outliers,track_id_list,MAIN_TRACK_ID,TIMEPOINT)
+fig,_ = plot_tracks_aligned_at_volume_drop_onset(df_outliers,track_id_list,MAIN_TRACK_ID,TIMEPOINT)
 
 save_name = f"S10_A_middle-dip_shape_{MAIN_TRACK_ID}"
-save_path = str(save_dir / save_name)
+save_path = Path( save_dir / save_name )
 for ext in ['.png','.pdf']:
     save_and_show_plot(save_path,ext,fig,transparent=False,keep_open=True)
 plt.show()
 
-#%%
-# S10 panel A, right
-# view data in timelapse feature explorer (TFE) with the following link 
-# TODO: update data path to be the final TFE url
-datapath = "https%3A%2F%2Fdev-aics-dtp-001.int.allencell.org%2Fassay-dev%2Fusers%2FFrick%2FPythonProjects%2Frepos%2Flocal_storage%2Ftimelapse_feature_explorer_datasets%2FTFE_new%2Fexploratory_dataset%2Fsmall%2Fmanifest.json"
-url = f"https://timelapse.allencell.org/viewer?dataset={datapath}&feature=change_in_volume_in_25_minute_window&t=239&filters=growth_outlier_filter%3A%3A3%2Cbaseline_colonies_dataset_filter%3A%3A3%2Cfullinterphase_dataset_filter%3A%3A3%2Clineageannotated_dataset_filter%3A%3A3%2Cvolume_jumps_right_magnitude_mask%3A%3A43.916%3A254.426&range=-200%2C200&color=matplotlib-purple_orange&palette-key=adobe&bg-sat=100&bg-brightness=100&fg-alpha=100&outlier-color=c0c0c0&outlier-mode=1&filter-color=dddddd&filter-mode=1&tab=scatter_plot&scalebar=1&timestamp=1&path=1&keep-range=1&scatter-range=all&scatter-x=scatterplot_time&scatter-y=volume"
-print(url)
-# save url as text file in figure folder
-with open(save_dir / 'S10_A_right_url.txt','w') as f:
-    f.write(url)
 # %%
 # S10 panel B, illustrate the effect of volume dip on transient growth rate
 
 df_full = filter_data.all_timepoints_full_tracks(df) # filter to only full tracks
 
-fig,ax = plt.subplots(nrows=2,ncols=1,figsize=(6.5,8))
-axlist = np.asarray([ax]) if type(ax) != np.ndarray else ax
-assert type(ax) == np.ndarray # for mypy
+fig,axlist_untyped = plt.subplots(nrows=2,ncols=1,figsize=(6.5,8))
+axlist = type_axlist(axlist_untyped)
 
 # add_time_point_lines=False,timepoint=None
 volume_dip_example_track = 86570
 
-main_track_list = [(volume_dip_example_track, 263)] #[(86570, 263),(75725, 239), (71532,131)]
+main_track_list = [(volume_dip_example_track, 263)]
 for main_track_id, timepoint in main_track_list:
     ax = axlist[0]
     ax = plot_track_with_volume_dip(ax,df_full,main_track_id,add_time_point_lines=True,timepoint=timepoint)
     ax = axlist[1]
     ax = plot_track_with_volume_dip(ax,df_full,main_track_id,xcol='index_sequence',ycol='dxdt_48_volume')
-    fig,axlist = adjust_axis_positions(fig,axlist,curr_pos=None,width=0.6,height=0.6,space=0.2,horizontal=False)
+    fig,axlist_untyped = adjust_axis_positions(fig,axlist,curr_pos=None,width=0.6,height=0.6,space=0.2,horizontal=False)
     for ext in ['.png','.pdf']:
         savepath = save_dir / f"S10B_track_{main_track_id}_volume_dip{ext}"
         save_and_show_plot(str(savepath),ext,fig,transparent=False,keep_open=True)
@@ -109,8 +100,8 @@
 #%%
 # S10 panel C step3
 df_track = df_full[df_full.track_id == volume_dip_example_track]
-fig,axlist = plt.subplots(2,1,figsize=(fw,fh),sharey=False)
-axlist = np.asarray([axlist]) if type(axlist) != np.ndarray else axlist
+fig,axlist_untyped = plt.subplots(2,1,figsize=(fw,fh),sharey=False)
+axlist = type_axlist(axlist_untyped)
 
 _ = plot_track_with_fit_line(df_track,axlist[0],
                                 ycol1='volume',
@@ -128,7 +119,8 @@
     ax.set_xticks(np.arange(0,20,4))
     ax.set_xlim(-2,xlimmax)
 
-fig,axlist = adjust_axis_positions(fig,axlist,curr_pos=None,width=0.6,height=0.6,space=0.2,horizontal=False)
+fig,axlist_untyped = adjust_axis_positions(fig,axlist,curr_pos=None,width=0.6,height=0.6,space=0.2,horizontal=False)
+axlist = type_axlist(axlist_untyped)
 axlist[0].text(0.05,0.99,f"track {volume_dip_example_track}",transform=axlist[0].transAxes,
         ha = 'left',va='top',fontsize=fs)
 axlist[0].legend(loc='lower left',bbox_to_anchor=(1.05,0.0),
@@ -138,9 +130,9 @@
                     )
 
 # now save
-savename = f"S10C_right-volume_fit_volume_fit_track{volume_dip_example_track}"
+savename: str = f"S10C_right-volume_fit_volume_fit_track{volume_dip_example_track}"
 for ext in ['.png','.pdf']:
-    save_path = str(save_dir / savename)
+    save_path = save_dir / savename
     save_and_show_plot(str(save_path),ext,fig,transparent=False,keep_open=True)
 plt.show()
 
@@ -151,9 +143,8 @@
 colony_list = ['small','medium','large']
 
 for threshold in [-50,0]:
-    fig,axlist = plt.subplots(1,1,figsize=(fw,fh))
-    axlist = np.asarray([axlist]) if type(axlist) != np.ndarray else axlist # for mypy
-    assert type(axlist) == np.ndarray # for mypy
+    fig,axlist_untyped = plt.subplots(1,1,figsize=(fw,fh))
+    axlist = type_axlist(axlist_untyped)
     ax = axlist[0]
     for ci,colony in enumerate(colony_list):
         dfcolony = df[df['colony'] == colony]
@@ -173,8 +164,6 @@
         yn = df_all['number_of_nuclei']
         y = yd / yn *100
 
-        print(np.where(df_all['number_of_dips'] > 5))
-
         zorderval = 1 if threshold !=0 else -1 # to ensure large colony is in front when it has fewer peaks
         ax.plot(x,y,label=colony,color=plotting_tools.COLONY_COLORS[colony],zorder=ci*1000*zorderval)
     ax.set_xlabel(f"{xlabel} {xunit}")
@@ -187,15 +176,12 @@
     ax.text(0.05,0.99,text_str,transform=ax.transAxes,
             ha = 'left',va='top',fontsize=fs)
 
-    fig,axlist = adjust_axis_positions(fig,axlist,curr_pos=None,width=0.9,height=0.5,space=0.075)
+    fig,axlist_untyped = adjust_axis_positions(fig,axlist,curr_pos=None,width=0.9,height=0.5,space=0.075)
     ax.legend(loc='center left', bbox_to_anchor=(1.05, 0.5),
                 fontsize=fs,frameon=False,
                 markerscale=1,handlelength=1,
                 labelspacing=0,
                 )
-    # if threshold == -50:
-    #     ax.set_yticks(np.arange(0,110,10))
-    #     ax.set_ylim(0,30)
     if threshold != -100:
         curr_ylim = ax.get_ylim()
 
@@ -209,9 +195,9 @@
     ax.set_xlabel('Movie time (hr)')
 
 
-    savename = save_dir / f"S10_D-volume_dips_over_time_all_colonies-{ycol}-{threshold}"
+    savename = f"S10_D-volume_dips_over_time_all_colonies-{ycol}-{threshold}"
     for ext in ['.png','.pdf']:
-        savepath = str(save_dir / f"{savename}")
+        savepath = save_dir / savename
         save_and_show_plot(str(savepath),ext,fig,transparent=False,keep_open=True)
 
 
@@ -228,7 +214,6 @@
     fig,ax = adjust_axis_positions(fig,ax,curr_pos=None,width=0.9,height=0.6,space=0.075)
 
     plt.suptitle(f"{ycol}")
-    # savepath = figdir / f"cell_cycle_bins_{ycol}_{xcol1}_{plot_type}.png"
     for ext in ['.png','.pdf']:
         savepath = save_dir / f"S10_E-cell_cycle_bins_for_only_{colony}_{ycol}_{ext}"
         save_and_show_plot(str(savepath),ext,fig,transparent=False,keep_open=True)
@@ -237,7 +222,6 @@
 #%%
 # S10 panel F and G
 colony='all_baseline'
-# color = "colony" if colony != "all_baseline" else "#808080"
 
 dfc = df_full[df_full["colony"] == colony] if colony != "all_baseline" else df_full
 for local_radius_str in ["90um", "whole_colony"]:

nuc_morph_analysis/analyses/volume/extra_dxdt_analysis/correlating_transient_growth_rates.py

@@ -61,7 +61,7 @@
         dfcat = pd.merge(df1,df2,on='index_sequence',suffixes=('_1','_2'),how='inner')
         dfcat.dropna(subset=[f'nanmean_1',f'nanmean_2'],inplace=True)
 
-        x = dfcat[f'nanmean_1']
+        x = dfcat[f'nanmean_1'].astype('array')
         y = dfcat[f'nanmean_2']
 
         reg = LinearRegression().fit(x.values.reshape(-1,1),y)

nuc_morph_analysis/analyses/volume/figure_5_s9_workflow.py

@@ -248,7 +248,6 @@
     fig,ax = adjust_axis_positions(fig,ax,curr_pos=None,width=0.9,height=0.6,space=0.075)
 
     plt.suptitle(f"{ycol}")
-    # savepath = figdir / f"cell_cycle_bins_{ycol}_{xcol1}_{plot_type}.png"
     for ext in ['.png','.pdf']:
         savepath = Path(figdir) / f"S5_E-cell_cycle_bins_for_only_{colony}_{ycol}_{ext}"
         save_and_show_plot(str(savepath),ext,fig,transparent=False,keep_open=True)

nuc_morph_analysis/analyses/volume/filter_out_dips.py

@@ -159,7 +159,6 @@ def find_peaks_and_collect_features(vol_det_array, vol_array, index_sequence_vec
     # convert the output to a dataframe
     dfout_list = [pd.DataFrame(x.values(),columns = index_sequence_vec, index=x.keys()).T for x in out]
     dfout_list = [x.reset_index().rename(columns={'index':'index_sequence'}).set_index('index_sequence') for x in dfout_list]
-    # keys = ['volume_dips_peak_mask_at_region','volume_dips_centers','volume_dips_has_peak','volume_dips_prom','volume_dips_left_base','volume_dips_right_base','volume_dips_y2_magnitude']
     dfout = pd.concat(dfout_list,axis=0,keys=track_id_vec, names=['track_id']).reset_index()
 
     return dfout
@@ -335,7 +334,6 @@ def filter_out_volume_dips(dfd, volume_cols, find_dips=True, use_detrended=True,
     # they are added during the pivot operation
     # we will drop these rows
     not_in_dfdi = dfmi.index.difference(dfdi.index)
-    # print(f"dropping {len(not_in_dfdi)} rows")
     dfmi.drop(not_in_dfdi, inplace=True)
 
     dfmi.loc[dfmi.index.values, "CellId"] = dfdi.loc[dfmi.index.values, "CellId"]

nuc_morph_analysis/analyses/volume/plot_help.py

@@ -34,7 +34,6 @@ def adjust_axis_positions(fig,ax,curr_pos=None,width=1,height=0.7,space=0.075,ke
     fw,fh = fig.get_size_inches()
     for ci,cax in enumerate(ax):
         # make the axis = 1.0" wide x 0.7" tall
-        print(width,height,fw,fh)
         if curr_pos is None:
             curr_pos = [1,1 +  height,width,height]
         else:
@@ -119,10 +118,7 @@ def plot_dfg(dfcc,xcol,ycol,labelstr,curr_ax,plot_type='mean',colorby=None,requi
     # remove rows with less than 10 counts
     dfg = group_and_extract(dfcc,xcol,ycol)
     dfgindex = dfg['count']<required_N
-    print(f" timepoints with less than {required_N} counts: {dfg[dfgindex].index.values}")
-    dfg= dfg[dfg['count'] >= required_N]
-    print(labelstr,dfg['count'].min(),dfg['count'].max(),dfg['count'].mean(),dfg['count'].sum(), "t=",dfg.shape[0])
-
+    dfg= dfg[dfg['count'] >= required_N]    
 
     xscale,xlabel,xunit,_ = get_plot_labels_for_metric(xcol)
     yscale,ylabel,yunit,_ = get_plot_labels_for_metric(ycol)

nuc_morph_analysis/lib/preprocessing/add_features.py

@@ -1,5 +1,6 @@
 from nuc_morph_analysis.analyses.lineage.get_features import lineage_trees
 from nuc_morph_analysis.lib.visualization.plotting_tools import get_plot_labels_for_metric
+from nuc_morph_analysis.lib.preprocessing.compute_change_over_time import run_script
 import numpy as np
 
 FRAME_COL = {"Ff": "A", "frame_transition": "B", "Fb": "C"}
@@ -683,3 +684,41 @@ def add_mean_features(df,
     multiplier_list = [get_plot_labels_for_metric(x)[0] for x in feature_list]
     df = add_mean_feature_over_trajectory(df, feature_list, multiplier_list)
     return df
+
+def add_volume_change_over_25_minute_window(df, bin_interval=5):
+    """
+    Adds a new column to the dataframe that quantifies how much the volume has changed relative to 
+    25 minutes in the past (units are pixels^3)
+    this is useful for identifying volume dips in all tracks (see Fig S10)
+
+    Parameters
+    ----------
+    df : pandas.DataFrame
+        The input dataframe.
+    bin_interval : int
+        represents the number of frames to compute change in volume over
+        default is 5 frames, which is 25 minutes
+
+    Returns
+    -------
+    df : pandas.DataFrame
+        The dataframe with the new column 'volume_change_over_25_minutes' added.
+        (units are pixels^3)
+    """
+    dfm = df.copy()
+    # run the compute_change_over_time workflow for a given bin_interval
+    dfm = run_script(dfm,['volume'], [bin_interval], time_location='end')
+    dfm['volume_change_over_25_minutes'] = dfm['dxdt_5_volume_end']*5
+
+    # now check that all columns in df have the same dtype as columns in dfm
+    for col in df.columns:
+        if dfm[col].dtype != df[col].dtype:
+            print(f"column {col} has dtype {dfm[col].dtype} in dfm and {df[col].dtype} in df")
+
+    if dfm.shape[0] != df.shape[0]:
+        raise Exception(
+            f"The loaded manifest has {df.shape[0]} rows and your \
+            final manifest has {dfm.shape[0]} rows.\
+            Please revise code to leave manifest rows unchanged."
+        )
+    return dfm

nuc_morph_analysis/lib/preprocessing/add_times.py

@@ -512,25 +512,15 @@ def digitize_time_column(df, minval, maxval, number_of_bins=None, step_size=None
     df[new_col] = dig_time_array
     return df
 
-def validate_dig_time_with_plot(time_array = np.linspace(0,1,1000), number_of_bins=6, old_method=False):
+def validate_dig_time_with_plot(time_array = np.linspace(0,1,1000), number_of_bins=6):
     """
     this visualizes how the input array is binned by plotting
     the input array (x-axis) vs the digitized array (y-axis)
     """
 
-    if old_method:
-        time_array = np.linspace(0,1,1000)
-        TIME_BIN = 1/number_of_bins
-        df_agg = pd.DataFrame({'normalized_time':time_array})
-        timedig_bins = np.arange(0, 1 + TIME_BIN, TIME_BIN)
-        inds = np.digitize(df_agg["normalized_time"], timedig_bins)
-        df_agg["dig_time"] = timedig_bins[inds - 1]
-        dig_time = df_agg['dig_time'].values
-        extrastr = '\n(old method)'
-    else:
-        bin_centers = determine_bin_centers(0,1,number_of_bins=number_of_bins)
-        dig_time = digitize_time_array(time_array,bin_centers)
-        extrastr = ''
+    bin_centers = determine_bin_centers(0,1,number_of_bins=number_of_bins)
+    dig_time = digitize_time_array(time_array,bin_centers)
+    extrastr = ''
     fig, ax = plt.subplots(figsize=(3,3))
     plt.plot(time_array,dig_time)
     plt.xlabel('time')

nuc_morph_analysis/lib/preprocessing/compute_change_over_time.py

@@ -131,28 +131,3 @@ def run_script(df=None, dxdt_feature_list = None, bin_interval_list=None, exclud
     assert dforig.index.name == "CellId"
     return dforig
 
-
-def add_dvdt_over_V(df,columns=None,volume_col = 'volume'):
-    """
-    adds dvdt over V for all timepoints if dxdt_{time}_volume columns exist
-
-    Parameters
-    ----------
-    df : pd.DataFrame
-        dataframe with columns  columns + ['volume']
-    columns : list
-        list of columns to compute change normalized by volume
-    volume_col : str
-        name of the volume column, default is 'volume'
-
-    Returns
-    -------
-    df : pd.DataFrame
-        dataframe with columns ['{column}_per_V'] added
-    """
-    if columns is None:
-        columns = [f"{DXDT_PREFIX}{bin_interval}_{feature}" for bin_interval in BIN_INTERVAL_LIST for feature in DXDT_FEATURE_LIST]
-
-    for col in columns:
-        df[f"{col}_per_V"] = df[col] / df[volume_col]
-    return df