Merge pull request #398 from AllenInstitute/analysis

Figure script and statistics

Author: alexpiet

visual_behavior_glm/GLM_analysis_tools.py

@@ -898,7 +898,7 @@ def build_weights_df(run_params,results_pivoted, cache_results=False,load_cache=
 
     # Make a metric of omission excitation/inhibition
     #weights_df['omission_excited'] = [np.sum(x[0:24]) for x in weights_df['omissions_weights']
-    weights_df['omissions_excited'] = weights_df_24.apply(lambda x: omission_excitation(x['omissions_weights']),axis=1)
+    weights_df['omissions_excited'] = weights_df.apply(lambda x: omission_excitation(x['omissions_weights']),axis=1)
 
     # Return weights_df
     return weights_df 

visual_behavior_glm/GLM_cell_metrics.py

@@ -120,7 +120,7 @@ def evaluate_against_metrics(results_metrics, ymetric='variance_explained_full',
         plt.savefig(filepath)
         filepath = os.path.join(run_params['figure_dir'], 'variance_explained_vs_'+xmetric+'_'+label+'.svg') 
         plt.savefig(filepath)
-        print(filepath)
+        print('Figure saved to: '+filepath)
 
 
 

visual_behavior_glm/GLM_fit_dev.py

@@ -186,9 +186,7 @@ def make_glm(fit, run_params, design, session):
 
     # Make Nested Model plot (rainbow plot)
     # A couple versions with more or less detail
-    schematic_df = gsm.plot_all_dropouts(run_params['version'])
-    schematic_df = gsm.plot_high_level_dropouts(run_params['version'])
-    schematic_df = gsm.plot_nice_dropouts(run_params['version'])
+    schematic_df = gsm.plot_dropouts_final(run_params['version'])
 
     # Make the platform paper schematic examples
     oeid = 967008471
@@ -259,7 +257,7 @@ def make_glm(fit, run_params, design, session):
     gvt.plot_population_averages_by_area(results_pivoted, run_params, dropouts_to_show=['omissions','omissions_positive','omissions_negative'],extra='_omissions_excitation')
     gvt.plot_population_averages_by_depth(results_pivoted, run_params, dropouts_to_show=['omissions','omissions_positive','omissions_negative'],extra='_omissions_excitation',area='VISp')
     gvt.plot_population_averages_by_depth(results_pivoted, run_params, dropouts_to_show=['omissions','omissions_positive','omissions_negative'],extra='_omissions_excitation',area='VISl')
-    gvt.plot_kernel_comparison_by_omission_excitation(weights_df_24, run_params_24)
+    gvt.plot_kernel_comparison_by_omission_excitation(weights_df, run_params)
 
 def get_analysis_dfs(VERSION):
     run_params = glm_params.load_run_json(VERSION)

visual_behavior_glm/GLM_schematic_plots.py

@@ -25,7 +25,10 @@ def change_breakdown_schematic(run_params):
     ax.tick_params(axis='x',labelsize=style['fs2'])
     ax.set_yticks([])
     plt.tight_layout()
-    plt.savefig('/allen/programs/braintv/workgroups/nc-ophys/visual_behavior/ophys_glm/figures/change_breakdown_schematic.svg')
+    filename = '/allen/programs/braintv/workgroups/nc-ophys/visual_behavior/ophys_glm/figures/change_breakdown_schematic.svg'
+    print('Figure saved to: '+filename)
+    plt.savefig(filename)
+
 
 def omission_breakdown_schematic(run_params):
     plt.figure(figsize=(2.5,1.5))
@@ -41,7 +44,9 @@ def omission_breakdown_schematic(run_params):
     ax.tick_params(axis='x',labelsize=style['fs2'])
     ax.set_yticks([])
     plt.tight_layout()
-    plt.savefig('/allen/programs/braintv/workgroups/nc-ophys/visual_behavior/ophys_glm/figures/omission_breakdown_schematic.svg')
+    filename = '/allen/programs/braintv/workgroups/nc-ophys/visual_behavior/ophys_glm/figures/omission_breakdown_schematic.svg'
+    print('Figure saved to: '+filename)
+    plt.savefig(filename)
 
 def get_example_style():
     style={
@@ -230,9 +235,13 @@ def plot_glm_example_inputs(g,times,style,run_params, ax=None, add_stimulus=True
     ax.spines['top'].set_visible(False)
     ax.spines['right'].set_visible(False)
     #plt.tight_layout()
-    plt.savefig('/allen/programs/braintv/workgroups/nc-ophys/visual_behavior/ophys_glm/figures/example_inputs.svg')
+
     if add_stimulus:
-        plt.savefig('/allen/programs/braintv/workgroups/nc-ophys/visual_behavior/ophys_glm/figures/example_inputs_add_stimulus.svg')
+        filename ='/allen/programs/braintv/workgroups/nc-ophys/visual_behavior/ophys_glm/figures/example_inputs_add_stimulus.svg'
+    else:
+        filename = '/allen/programs/braintv/workgroups/nc-ophys/visual_behavior/ophys_glm/figures/example_inputs.svg'
+    print('Figure saved to: '+filename)
+    plt.savefig(filename)
 
 def plot_glm_example_components(g, cell_specimen_id, times, style):
     fig = plt.figure(figsize=(8,4))
@@ -339,11 +348,15 @@ def plot_glm_example_trace(g,cell_specimen_id,times,style,include_events=True,ax
     ax.set_ylim(-0.035,.9)
     ax.set_xlim(times)
     #plt.tight_layout()
+    
     if model is not None:
-        plt.savefig('/allen/programs/braintv/workgroups/nc-ophys/visual_behavior/ophys_glm/figures/example_trace_'+model+'.svg')
+        filename='/allen/programs/braintv/workgroups/nc-ophys/visual_behavior/ophys_glm/figures/example_trace_'+model+'.svg'
         plt.savefig('/allen/programs/braintv/workgroups/nc-ophys/visual_behavior/ophys_glm/figures/example_trace_'+model+'.png')
     else:
-        plt.savefig('/allen/programs/braintv/workgroups/nc-ophys/visual_behavior/ophys_glm/figures/example_trace.svg')
+        filename='/allen/programs/braintv/workgroups/nc-ophys/visual_behavior/ophys_glm/figures/example_trace.svg'
+
+    print('Figure saved to: '+filename)
+    plt.savefig(filename)
     return
 
 def plot_all_dropouts(VERSION):

visual_behavior_glm/GLM_visualization_tools.py

@@ -594,11 +594,16 @@ def var_explained_by_experience(results_pivoted, run_params,threshold = 0):
     plt.tick_params(axis='both',labelsize=14)
     plt.tight_layout() 
     if threshold !=0:
-        plt.savefig(run_params['figure_dir']+'/variance_explained_by_experience_filtered.svg')
+        filename = run_params['figure_dir']+'/variance_explained_by_experience_filtered.svg'
         plt.savefig(run_params['figure_dir']+'/variance_explained_by_experience_filtered.png')
+        print('Figure saved to: ' + filename)
+        plt.savefig(filename)
     else:
-        plt.savefig(run_params['figure_dir']+'/variance_explained_by_experience.svg')
+        filename = run_params['figure_dir']+'/variance_explained_by_experience.svg'
         plt.savefig(run_params['figure_dir']+'/variance_explained_by_experience.png')
+        print('Figure saved to: ' + filename)
+        plt.savefig(filename)
+    return results_pivoted.groupby(['cell_type','experience_level'])['variance_explained_percent'].describe()
 
 def compare_var_explained_by_version(results=None, fig=None, ax=None, test_data=True, figsize=(9,5), use_violin=True,cre=None,metric='Full',show_equipment=True,zoom_xlim=True,sort_by_signal=True):
     '''
@@ -1850,13 +1855,19 @@ def plot_kernel_comparison_by_experience(weights_df, run_params, kernel,threshol
     fig_f.savefig(run_params['fig_kernels_dir']+'/'+kernel+'_familiar_kernel'+extra+'.svg')
     fig_n.savefig(run_params['fig_kernels_dir']+'/'+kernel+'_novel1_kernel'+extra+'.svg')
     fig_np.savefig(run_params['fig_kernels_dir']+'/'+kernel+'_novelp1_kernel'+extra+'.svg')
+    print('Figure saved to: '+run_params['fig_kernels_dir']+'/'+kernel+'_familiar_kernel'+extra+'.svg')
+    print('Figure saved to: '+run_params['fig_kernels_dir']+'/'+kernel+'_novel1_kernel'+extra+'.svg')
+    print('Figure saved to: '+run_params['fig_kernels_dir']+'/'+kernel+'_novelp1_kernel'+extra+'.svg')
 
     k, fig_v , ax_v =plot_kernel_comparison(weights_df,run_params,kernel,save_results=False,session_filter=['Familiar','Novel 1','Novel >1'],cell_filter='Vip-IRES-Cre',compare=['experience_level'],threshold=threshold,drop_threshold=drop_threshold)   
     k, fig_s , ax_s =plot_kernel_comparison(weights_df,run_params,kernel,save_results=False,session_filter=['Familiar','Novel 1','Novel >1'],cell_filter='Sst-IRES-Cre',compare=['experience_level'],threshold=threshold,drop_threshold=drop_threshold)   
     k, fig_e , ax_e =plot_kernel_comparison(weights_df,run_params,kernel,save_results=False,session_filter=['Familiar','Novel 1','Novel >1'],cell_filter='Slc17a7-IRES2-Cre',compare=['experience_level'],threshold=threshold,drop_threshold=drop_threshold)   
     fig_v.savefig(run_params['fig_kernels_dir']+'/'+kernel+'_vip_kernel'+extra+'.svg')
     fig_s.savefig(run_params['fig_kernels_dir']+'/'+kernel+'_sst_kernel'+extra+'.svg')
     fig_e.savefig(run_params['fig_kernels_dir']+'/'+kernel+'_exc_kernel'+extra+'.svg')
+    print('Figure saved to: '+run_params['fig_kernels_dir']+'/'+kernel+'_vip_kernel'+extra+'.svg')
+    print('Figure saved to: '+run_params['fig_kernels_dir']+'/'+kernel+'_sst_kernel'+extra+'.svg')
+    print('Figure saved to: '+run_params['fig_kernels_dir']+'/'+kernel+'_exc_kernel'+extra+'.svg')
 
 def plot_kernel_comparison_by_omission_excitation(weights_df, run_params):
     plot_kernel_comparison(weights_df,run_params,'omissions',session_filter=['Familiar'],cell_filter='Slc17a7-IRES2-Cre',compare=['omissions_excited'])
@@ -2884,6 +2895,7 @@ def plot_kernel_heatmap_with_dropout(vip_table, sst_table, slc_table, time_vec,k
     #dax3.set_xticklabels(['Coding\n Score'],rotation=-70,fontsize=16)
     filename = os.path.join(run_params['fig_kernels_dir'],kernel+'_heatmap_with_dropout'+extra+'.svg')
     plt.savefig(filename)
+    print('Figure saved to: '+filename)
     filename = os.path.join(run_params['fig_kernels_dir'],kernel+'_heatmap_with_dropout'+extra+'.png')
     plt.savefig(filename)
     return zlims
@@ -3497,6 +3509,7 @@ def plot_population_averages_by_depth(results_pivoted, run_params, dropouts_to_s
     experience_levels = np.sort(results_pivoted.experience_level.unique())
     colors = project_colors()
 
+    summary = {}
     # Iterate cell types and make a plot for each
     for cell_type in cell_types:
         if len(dropouts_to_show) ==3:
@@ -3505,9 +3518,11 @@ def plot_population_averages_by_depth(results_pivoted, run_params, dropouts_to_s
             fig, ax = plt.subplots(1,len(dropouts_to_show),figsize=(10.8,4), sharey=sharey)
         all_data = results_pivoted.query('cell_type ==@cell_type')
         stats = {}
+        summary[cell_type + ' data'] = {}
         # Iterate dropouts and plot each by experience
         for index, feature in enumerate(dropouts_to_show):
             stats[feature] = test_significant_dropout_averages_by_depth(all_data,feature)
+            summary[cell_type+' data'][feature] = all_data.groupby(['experience_level','coarse_binned_depth'])[feature].describe()
             # Plot all cells in active sessions 
             ax[index] = sns.pointplot(
                 data = all_data,
@@ -3560,10 +3575,12 @@ def plot_population_averages_by_depth(results_pivoted, run_params, dropouts_to_s
         ax[0].set_ylabel('Coding Score',fontsize=20)
         plt.suptitle(cell_type+', '+area,fontsize=20)
         fig.tight_layout() 
-        plt.savefig(run_params['figure_dir']+'/dropout_average_by_depth_'+cell_type[0:3]+extra+'.svg')
+        filename = run_params['figure_dir']+'/dropout_average_by_depth_'+cell_type[0:3]+extra+'.svg'
         plt.savefig(run_params['figure_dir']+'/dropout_average_by_depth_'+cell_type[0:3]+extra+'.png')
-
-
+        print('Figure saved to: '+filename)
+        plt.savefig(filename)
+        summary[cell_type + ' stats'] = stats
+    return summary
 
 def plot_population_averages_by_area(results_pivoted, run_params, dropouts_to_show = ['all-images','omissions','behavioral','task'],sharey=False,include_zero_cells=True,add_stats=True,extra='',equipment="mesoscope"):
     '''
@@ -3609,7 +3626,7 @@ def plot_population_averages_by_area(results_pivoted, run_params, dropouts_to_sh
     cell_types = results_pivoted.cell_type.unique()
     experience_levels = np.sort(results_pivoted.experience_level.unique())
     colors = project_colors()
-
+    summary = {}
     # Iterate cell types and make a plot for each
     for cell_type in cell_types:
         if len(dropouts_to_show) == 3:
@@ -3618,10 +3635,12 @@ def plot_population_averages_by_area(results_pivoted, run_params, dropouts_to_sh
             fig, ax = plt.subplots(1,len(dropouts_to_show),figsize=(10.8,4), sharey=sharey)
         all_data = results_pivoted.query('cell_type ==@cell_type')
         stats = {}
+        summary[cell_type + ' data'] = {}
         # Iterate dropouts and plot each by experience
         for index, feature in enumerate(dropouts_to_show):
             stats[feature] = test_significant_dropout_averages_by_area(all_data,feature)
             # Plot all cells in active sessions 
+            summary[cell_type+' data'][feature] = all_data.groupby(['experience_level','targeted_structure'])[feature].describe()
             ax[index] = sns.pointplot(
                 data = all_data,
                 x = 'experience_level',
@@ -3672,9 +3691,12 @@ def plot_population_averages_by_area(results_pivoted, run_params, dropouts_to_sh
         ax[0].set_ylabel('Coding Score',fontsize=20)
         plt.suptitle(cell_type,fontsize=20)
         fig.tight_layout() 
-        plt.savefig(run_params['figure_dir']+'/dropout_average_by_area_'+cell_type[0:3]+extra+'.svg')
+        filename = run_params['figure_dir']+'/dropout_average_by_area_'+cell_type[0:3]+extra+'.svg'
         plt.savefig(run_params['figure_dir']+'/dropout_average_by_area_'+cell_type[0:3]+extra+'.png')
-
+        plt.savefig(filename)
+        print('Figure saved to: '+filename)
+        summary[cell_type + ' stats'] = stats
+    return summary
 
 
 def plot_population_averages(results_pivoted, run_params, dropouts_to_show = ['all-images','omissions','behavioral','task'],sharey=True,include_zero_cells=True,boxplot=False,add_stats=True,extra='',strict_experience_matching=False,plot_by_cell_type=False):
@@ -3755,7 +3777,9 @@ def plot_population_averages(results_pivoted, run_params, dropouts_to_show = ['a
             ax[index].tick_params(axis='y',labelsize=16)
         ax[0].set_ylabel('Coding Score',fontsize=20)
         plt.tight_layout()
-        plt.savefig(run_params['figure_dir']+'/dropout_average_combined'+extra+'.svg')  
+        filename = run_params['figure_dir']+'/dropout_average_combined'+extra+'.svg'
+        print('Figure saved to: '+filename)
+        plt.savefig(filename)  
         #plt.savefig(run_params['figure_dir']+'/dropout_average_combined'+extra+'.png')  
 
         # Iterate cell types and make a plot for each
@@ -3862,11 +3886,14 @@ def plot_population_averages(results_pivoted, run_params, dropouts_to_show = ['a
                     ax[index].set_ylim(0,ytop*1.2)
             ax[0].set_ylabel('Coding Score',fontsize=18)
             plt.suptitle(cell_type,fontsize=18)
-            fig.tight_layout() 
-            plt.savefig(run_params['figure_dir']+'/dropout_average_'+cell_type[0:3]+extra+'.svg')
+            fig.tight_layout()
+            filename = run_params['figure_dir']+'/dropout_average_'+cell_type[0:3]+extra+'.svg' 
+            plt.savefig(filename)
+            print('Figure saved to: '+filename)
 
     # Repeat the plots but transposed
     # Iterate cell types and make a plot for each
+    summary_data = {}
     for index, feature in enumerate(dropouts_to_show):   
         fig, ax = plt.subplots(1,4,figsize=(10.8,4), sharey=sharey) 
 
@@ -3893,6 +3920,7 @@ def plot_population_averages(results_pivoted, run_params, dropouts_to_show = ['a
         ax[3].spines['top'].set_visible(False)
         ax[3].spines['right'].set_visible(False)
 
+        summary_data[feature + ' data'] = {}
         stats = {}
         # Iterate dropouts and plot each by experience
         for cindex, cell_type in enumerate(cell_types):
@@ -3902,6 +3930,8 @@ def plot_population_averages(results_pivoted, run_params, dropouts_to_show = ['a
                 strict_matched_data = all_data.query('cell_specimen_id in @strict_matched_cells') 
             anova, tukey = test_significant_dropout_averages(all_data,feature)
             stats[cell_type]=(anova, tukey)
+            summary_data[feature+' data'][cell_type+' all data'] = all_data.groupby(['experience_level'])[feature].describe()
+            summary_data[feature+' data'][cell_type+' matched data'] = matched_data.groupby(['experience_level'])[feature].describe()
             # Plot all cells in active sessions
             if boxplot:
                 ax[cindex] = sns.boxplot(
@@ -3945,6 +3975,7 @@ def plot_population_averages(results_pivoted, run_params, dropouts_to_show = ['a
 
             if strict_experience_matching:
                 # Plot cells in matched active sessions
+                summary_data[feature+' data'][cell_type+' strict matched data'] = strict_matched_data.groupby(['experience_level'])[feature].describe()
                 ax[cindex] = sns.pointplot(
                     data = strict_matched_data,
                     x = 'experience_level',
@@ -4002,8 +4033,12 @@ def plot_population_averages(results_pivoted, run_params, dropouts_to_show = ['a
         ax[2].set_ylim(bottom=0)
         ax[3].set_ylim(bottom=0)
         fig.tight_layout() 
-        plt.savefig(run_params['figure_dir']+'/dropout_average_'+clean_feature.replace(' ','_')+extra+'.svg')
+        filename = run_params['figure_dir']+'/dropout_average_'+clean_feature.replace(' ','_')+extra+'.svg'
+        plt.savefig(filename)
+        print('Figure saved to: '+filename)
+        summary_data[feature+' stats'] = stats
 
+    return summary_data
 
 
 def test_significant_dropout_averages(data,feature):
@@ -4166,13 +4201,16 @@ def plot_dropout_individual_population(results, run_params,ax=None,palette=None,
     if add_title:
         plt.title(run_params['version'])
     if use_violin:
-        plt.savefig(run_params['figure_dir']+'/dropout_individual.svg')
+        filename = run_params['figure_dir']+'/dropout_individual.svg'
     elif use_single:
-        plt.savefig(run_params['figure_dir']+'/dropout_individual_boxplot_single.svg')
+        filename = run_params['figure_dir']+'/dropout_individual_boxplot_single.svg'
     else:
-        plt.savefig(run_params['figure_dir']+'/dropout_individual_boxplot.svg')
+        filename = run_params['figure_dir']+'/dropout_individual_boxplot.svg'
         plt.savefig(run_params['figure_dir']+'/dropout_individual_boxplot.png')
 
+    plt.savefig(filename)
+    print('Figure saved to: '+filename)
+    return data_to_plot.groupby(['cre_line','dropout'])['explained_variance'].describe()
 
 def plot_dropout_summary_population(results, run_params,dropouts_to_show =  ['all-images','omissions','behavioral','task'],ax=None,palette=None,use_violin=False,add_median=True,include_zero_cells=True,add_title=False): 
     '''
@@ -4281,8 +4319,11 @@ def plot_dropout_summary_population(results, run_params,dropouts_to_show =  ['al
     if use_violin:
         plt.savefig(run_params['figure_dir']+'/dropout_summary.svg')
     else:
-        plt.savefig(run_params['figure_dir']+'/dropout_summary_boxplot.svg')
+        filename = run_params['figure_dir']+'/dropout_summary_boxplot.svg'
+        print('Figure saved to: '+filename)
+        plt.savefig(filename)
         plt.savefig(run_params['figure_dir']+'/dropout_summary_boxplot.png')
+    return data_to_plot.groupby(['cre_line','dropout'])['explained_variance'].describe() 
 
 def plot_fraction_summary_population(results_pivoted, run_params,sharey=True,omissions_excitation=False):
     # compute coding fractions
@@ -4352,12 +4393,16 @@ def plot_fraction_summary_population(results_pivoted, run_params,sharey=True,omi
     ax[0].set_ylabel('Fraction of cells \n coding for ',fontsize=20)
     plt.tight_layout()
     if omissions_excitation:
-        plt.savefig(run_params['figure_dir']+'/coding_fraction_omissions_summary.svg')  
-        plt.savefig(run_params['figure_dir']+'/coding_fraction_omissions_summary.png')  
+        filename = run_params['figure_dir']+'/coding_fraction_omissions_summary.svg'  
+        plt.savefig(filename)  
+        plt.savefig(run_params['figure_dir']+'/coding_fraction_omissions_summary.png') 
+        print('Figure saved to: '+filename) 
     else:
-        plt.savefig(run_params['figure_dir']+'/coding_fraction_summary.svg')  
+        filename = run_params['figure_dir']+'/coding_fraction_summary.svg'
+        plt.savefig(filename)  
         plt.savefig(run_params['figure_dir']+'/coding_fraction_summary.png')  
- 
+        print('Figure saved to: '+filename) 
+    return summary_df 
 
 def make_cosyne_schematic(glm,cell=1028768972,t_range=5,time_to_plot=3291,alpha=.25):
     '''