diff --git a/unravel/image_tools/uniq_intensities.py b/unravel/image_tools/uniq_intensities.py index 8b5191b9..6ac6e761 100755 --- a/unravel/image_tools/uniq_intensities.py +++ b/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)() \ No newline at end of file + # Configuration.verbose = args.verbose + # print_cmd_and_times(main)() + main() \ No newline at end of file diff --git a/unravel/region_stats/regional_IF_mean_intensities.py b/unravel/region_stats/regional_IF_mean_intensities.py index f689c513..00d7ccc6 100755 --- a/unravel/region_stats/regional_IF_mean_intensities.py +++ b/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 .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() diff --git a/unravel/region_stats/regional_IF_mean_intensities_summary.py b/unravel/region_stats/regional_IF_mean_intensities_summary.py index a297e227..8122b3e5 100755 --- a/unravel/region_stats/regional_IF_mean_intensities_summary.py +++ b/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)