t-tests. Plot pos/neg values. Can use img_unique to get region list. …

…Guides updated

unravel/image_tools/uniq_intensities.py

@@ -1,19 +1,26 @@
 #!/usr/bin/env python3
 
 """
-Print list of IDs for clusters > X voxels
+Print list of unique intensities greater than 0
 
-Usage:
-    uniq_intensities.py -i path/input_img.nii.gz -m 100 -id -s
-    uniq_intensities.py -i path/input_img.nii.gz -m 100 -id
+Usage for printing all non-zero intensities: 
+    uniq_intensities.py -i path/input_img.nii.gz
+
+Usage for printing the number of voxels for each intensity that is present: 
+    uniq_intensities.py -i path/input_img.nii.gz
+
+Usage for checking which clusters are present if the min cluster size was 100 voxels: 
+    uniq_intensities.py -i path/input_img.nii.gz -m 100
 """
 
 import argparse
+import nibabel as nib
+import numpy as np
 from rich import print
 from rich.traceback import install
 
 from unravel.core.argparse_utils import SuppressMetavar, SM
-from unravel.core.config import Configuration
+# from unravel.core.config import Configuration
 from unravel.core.img_io import load_3D_img
 from unravel.core.img_tools import cluster_IDs
 from unravel.core.utils import print_cmd_and_times
@@ -23,21 +30,25 @@ def parse_args():
     parser = argparse.ArgumentParser(formatter_class=SuppressMetavar)
     parser.add_argument('-i', '--input', help='path/input_img.nii.gz', required=True, action=SM)
     parser.add_argument('-m', '--minextent', help='Min cluster size in voxels (Default: 1)', default=1, action=SM, type=int)
-    parser.add_argument('-id', '--print_IDs', help='Print cluster IDs. Default: True', default=True, action='store_true')
     parser.add_argument('-s', '--print_sizes', help='Print cluster IDs and sizes. Default: False', default=False, action='store_true')
-    parser.add_argument('-v', '--verbose', help='Increase verbosity. Default: False', action='store_true', default=False)
+    # parser.add_argument('-v', '--verbose', help='Increase verbosity. Default: False', action='store_true', default=False)
     parser.epilog = __doc__
     return parser.parse_args()
 
 def main():
     args = parse_args()
 
-    img = load_3D_img(args.input)
+    if str(args.input).endswith(".nii.gz"):
+        nii = nib.load(args.input)
+        img = np.asanyarray(nii.dataobj, dtype=np.uint16).squeeze()
+    else: 
+        img = load_3D_img(args.input)
 
-    cluster_IDs(img, min_extent=args.minextent, print_IDs=args.print_IDs, print_sizes=args.print_sizes)
+    cluster_IDs(img, min_extent=args.minextent, print_IDs=True, print_sizes=args.print_sizes)
 
 if __name__ == '__main__': 
     install()
     args = parse_args()
-    Configuration.verbose = args.verbose
-    print_cmd_and_times(main)()
+    # Configuration.verbose = args.verbose
+    # print_cmd_and_times(main)()
+    main()

unravel/region_stats/regional_IF_mean_intensities.py

@@ -1,7 +1,15 @@
 #!/usr/bin/env python3
 
 """
-Measure mean intensity of immunofluorescence staining in brain regions.
+Measure mean intensity of immunofluorescence staining in brain regions in atlas space.
+
+Usage: 
+    atlas=/SSD4/ET/atlas/gubra_ano_split_25um_w_A13_RH.nii.gz ; regions=$(img_unique -i $atlas) ; for i in <asterisk>.nii.gz ; do rstats_IF_mean -i $i -a $atlas -r $regions ; done
+    mkdir regional_mean_intensities
+    mv *_regional_mean_intensities.csv regional_mean_intensities/
+    cd regional_mean_intensities
+    rstats_IF_mean_summary 
+
 """
 
 import argparse 
@@ -13,12 +21,11 @@
 
 
 def parse_args():
-    DEFAULT_ATLAS = '/usr/local/unravel/atlases/gubra/gubra_ano_combined_25um.nii.gz'
     parser = argparse.ArgumentParser(formatter_class=SuppressMetavar)
-    parser.add_argument('-a', '--atlas', help='path/atlas.nii.gz', default=DEFAULT_ATLAS, action=SM)
     parser.add_argument('-i', '--input', help='path/atlas_space_immunofluorescence_image.nii.gz', required=True, action=SM)
-    parser.add_argument('-o', '--output', help='path/name.csv', default=None, action=SM)
+    parser.add_argument('-a', '--atlas', help='path/atlas.nii.gz', required=True, action=SM)
     parser.add_argument('-r', '--regions', nargs='*', type=int, help='Space-separated list of region intensities to process', default=None)
+    parser.add_argument('-o', '--output', help='path/name.csv', default=None, action=SM)
     parser.epilog = __doc__
     return parser.parse_args()
 

unravel/region_stats/regional_IF_mean_intensities_summary.py

@@ -1,10 +1,13 @@
 #!/usr/bin/env python3
 
 """
-Plot mean IF intensity for a given region intensity ID for 3+ groups (only works for positive data)
+Plot mean IF intensities and comparisons for each region intensity ID (only works for positive data).
 
-Usage:
-    regional_IF_mean_intensities_summary.py -r 1 --order group3 group2 group1 --labels Group_3 Group_2 Group_1 -t dunnnett
+Usage for Tukey's tests:
+    atlas=/SSD4/ET/atlas/gubra_ano_split_25um_w_A13_RH.nii.gz ; rstats_IF_mean_summary -r $(img_unique -i $atlas) --order group3 group2 group1 --labels Group_3 Group_2 Group_1
+
+Usage for t-test:
+    atlas=/SSD4/ET/atlas/gubra_ano_split_25um_w_A13_RH.nii.gz ; rstats_IF_mean_summary --region_ids $(img_unique -i $atlas) --order group2 group1 --labels Group_2 Group_1 -t ttest
 """
 
 import argparse
@@ -18,19 +21,20 @@
 from rich import print
 from rich.traceback import install
 from pathlib import Path
-from scipy.stats import dunnett
+from scipy.stats import ttest_ind, dunnett
 from statsmodels.stats.multicomp import pairwise_tukeyhsd
 
 from unravel.core.argparse_utils import SuppressMetavar, SM
 
 
+
 def parse_args():
     parser = argparse.ArgumentParser(formatter_class=SuppressMetavar)
     parser.add_argument('--region_ids', nargs='*', type=int, help='List of region intensity IDs (Default: process all regions from the lut CSV)', action=SM)
     parser.add_argument('-l', '--lut', help='LUT csv name (in unravel/core/csvs/). Default: gubra__region_ID_side_name_abbr.csv', default="gubra__region_ID_side_name_abbr.csv", action=SM)
     parser.add_argument('--order', nargs='*', help='Group Order for plotting (must match 1st word of CSVs)', action=SM)
     parser.add_argument('--labels', nargs='*', help='Group Labels in same order', action=SM)
-    parser.add_argument('-t', '--test', help='Choose between "tukey" and "dunnett" post-hoc tests. (Default: tukey)', default='tukey', choices=['tukey', 'dunnett'], action=SM)
+    parser.add_argument('-t', '--test', help='Choose between "tukey", "dunnett", and "ttest" post-hoc tests. (Default: tukey)', default='tukey', choices=['tukey', 'dunnett', 'ttest'], action=SM)
     parser.add_argument('-alt', "--alternate", help="Number of tails and direction for Dunnett's test {'two-sided', 'less' (means < ctrl), 'greater'}. Default: two-sided", default='two-sided', action=SM)
     parser.add_argument('-s', '--show_plot', help='Show plot if flag is provided', action='store_true')
     parser.epilog = __doc__
@@ -70,6 +74,22 @@ def get_all_region_ids(csv_path):
     region_df = pd.read_csv(csv_path)
     return region_df["Region_ID"].tolist()
 
+def perform_t_tests(df, order):
+    """Perform t-tests between groups in the DataFrame."""
+    comparisons = []
+    for i in range(len(order)):
+        for j in range(i + 1, len(order)):
+            group1, group2 = order[i], order[j]
+            data1 = df[df['group'] == group1]['mean_intensity']
+            data2 = df[df['group'] == group2]['mean_intensity']
+            t_stat, p_value = ttest_ind(data1, data2)
+            comparisons.append({
+                'group1': group1,
+                'group2': group2,
+                'p-adj': p_value
+            })
+    return pd.DataFrame(comparisons)
+
 def plot_data(region_id, order=None, labels=None, csv_path=None, test_type='tukey', show_plot=False, alt='two-sided'):
     df = load_data(region_id)
     region_name, region_abbr = get_region_details(region_id, csv_path)
@@ -109,6 +129,7 @@ def plot_data(region_id, order=None, labels=None, csv_path=None, test_type='tuke
 
     # Formatting
     ax.set_ylabel('Mean IF Intensity', weight='bold')
+    ax.set_xticks(np.arange(len(df['group_label'].unique())))
     ax.set_xticklabels(ax.get_xticklabels(), weight='bold')
     ax.tick_params(axis='both', which='major', width=2)
     ax.spines['top'].set_visible(False)
@@ -120,7 +141,8 @@ def plot_data(region_id, order=None, labels=None, csv_path=None, test_type='tuke
     sns.swarmplot(x='group_label', y='mean_intensity', hue='group', data=df, palette=group_colors, size=8, linewidth=1, edgecolor='black')
 
     # Remove the legend created by hue
-    ax.legend_.remove()
+    if ax.legend_:
+        ax.legend_.remove()
 
     # Perform the chosen post-hoc test
     if test_type == 'tukey':
@@ -138,8 +160,11 @@ def plot_data(region_id, order=None, labels=None, csv_path=None, test_type='tuke
             'p-adj': test_stats.pvalue
         })
         test_df['reject'] = test_df['p-adj'] < 0.05
-    significant_comparisons = test_df[test_df['reject'] == True]
+    elif test_type == 'ttest':
+        test_df = perform_t_tests(df, order)
+        test_df['reject'] = test_df['p-adj'] < 0.05
 
+    significant_comparisons = test_df[test_df['reject'] == True]
     y_max = df['mean_intensity'].max()
     y_min = df['mean_intensity'].min()
     height_diff = (y_max - y_min) * 0.1
@@ -166,16 +191,37 @@ def plot_data(region_id, order=None, labels=None, csv_path=None, test_type='tuke
         plt.text((x1+x2)*.5, y_pos + 0.8*height_diff, sig, horizontalalignment='center', size='xx-large', color='black', weight='bold')
         y_pos += 3 * height_diff
 
-    plt.ylim(0, y_pos + 2*height_diff)
+
+
+    # Calculate y-axis limits
+    y_max = df['mean_intensity'].max()
+    y_min = df['mean_intensity'].min()
+    height_diff = (y_max - y_min) * 0.1
+    y_pos = y_max + 0.5 * height_diff
+
+    # Ensure the y-axis starts from the minimum value, allowing for negative values
+    plt.ylim(y_min - 2 * height_diff, y_pos + 2 * height_diff)
+
+
+
+    # plt.ylim(0, y_pos + 2*height_diff)
     ax.set_xlabel(None)
 
     # Save the plot
+    output_folder = Path('regional_IF_mean_summary')  # Replace 'output_folder_name' with your desired folder name
+    output_folder.mkdir(parents=True, exist_ok=True)
+
     title = f"{region_name} ({region_abbr})"
     wrapped_title = textwrap.fill(title, 42)  # wraps at x characters. Adjust as needed.
     plt.title(wrapped_title)
     plt.tight_layout()
     region_abbr = region_abbr.replace("/", "-") # Replace problematic characters for file paths
-    plt.savefig(f'region_{region_id}_{region_abbr}.pdf')
+
+    is_significant = not significant_comparisons.empty
+    file_prefix = '_' if is_significant else ''
+    file_name = f"{file_prefix}region_{region_id}_{region_abbr}.pdf"
+    plt.savefig(output_folder / file_name)
+
     plt.close()
 
     if show_plot:
@@ -197,7 +243,6 @@ def main():
 
     # Process each region ID
     for region_id in region_ids_to_process:
-        csv_path = Path(__file__).parent.parent / 'csvs' / args.lut
 
         plot_data(region_id, args.order, args.labels, csv_path=lut, test_type=args.test, show_plot=args.show_plot, alt=args.alternate)