diff --git a/README.md b/README.md index 1e63437..6e375ad 100644 --- a/README.md +++ b/README.md @@ -8,6 +8,17 @@ The aim is to evaluate the impact of denoising strategy on functional connectivi **Preprint of the manuscript is now on [biorxiv](https://www.biorxiv.org/content/10.1101/2023.04.18.537240). The reporducible Jupyter Book preprint is on [NeuroLibre](https://neurolibre.org/papers/10.55458/neurolibre.00012).** +## Recommandations for those who thought this project is a software + +Bad news, this is not a software but a research project. +It's more similar to your regular data science project. +In other words, the code in this repository reflects the research done for the manuscript, and is not suitable for production level application. + +Some useful part of the code has been extracted and further reviewed within SIMEXP lab for deplyment on generic fmriprep derivatives as docker images. + + - *time series and connectome workflow*: [`giga_connectome`](https://github.com/SIMEXP/giga_connectome). + - *motion quality control metrics*: [`giga_auto_qc`](https://github.com/SIMEXP/giga_auto_qc). + ## Quick start ```bash @@ -37,12 +48,3 @@ make book - Custom code is located in `fmriprep_denoise/`. This project is installable. - Preprocessing SLURM scripts, and scripts for creating figure for manuscript are in `scripts/`. - - -## Poster and presentations - -The results will be presented at QBIN science day 2023 as a flash talk, and OHBM 2023 Montreal as a poster :tada:. - -The preliminary results were presented at OHBM 2022 as a poster. - -![spoiler](./content/images/ohbm2022_abstract_head.png) diff --git a/content/images/denoise-benchmark-workflow.drawio b/content/images/denoise-benchmark-workflow.drawio index ee33526..b727b76 100644 --- a/content/images/denoise-benchmark-workflow.drawio +++ b/content/images/denoise-benchmark-workflow.drawio @@ -1,6 +1,6 @@ - + - + @@ -34,7 +34,7 @@ - + @@ -220,38 +220,38 @@ - + - + - + - + - + - + - + - + - + - + @@ -306,7 +306,7 @@ - + @@ -321,10 +321,10 @@ - + - - + + @@ -335,7 +335,7 @@ - + @@ -366,7 +366,7 @@ - + diff --git a/content/images/fig-1-masker.png b/content/images/fig-1-masker.png index a80a2f6..cca8121 100644 Binary files a/content/images/fig-1-masker.png and b/content/images/fig-1-masker.png differ diff --git a/fmriprep_denoise/visualization/degrees_of_freedom_loss.py b/fmriprep_denoise/visualization/degrees_of_freedom_loss.py index fa18de6..b94dfb1 100644 --- a/fmriprep_denoise/visualization/degrees_of_freedom_loss.py +++ b/fmriprep_denoise/visualization/degrees_of_freedom_loss.py @@ -24,15 +24,80 @@ def load_data(path_root, datasets, criteria_name, fmriprep_version): fd_thresh=criteria["fd_thresh"], proportion_thresh=criteria["proportion_thresh"], ) - confounds_phenotypes[dataset] = confounds_phenotype[strategy_order] + df_plotting, full_length = _organise_data(confounds_phenotype[strategy_order]) + confounds_phenotypes[dataset] = { + "group_values": groups, + "participant_labels": participant_groups, + "confounds_phenotype": df_plotting, + "confounds_stats": confounds_phenotype[strategy_order], + "full_length": full_length, + } return confounds_phenotypes -def plot_stats(confounds_phenotypes): +def _plot_single_report(confounds_phenotype, full_length, palette, ax, title): + # change up the data a bit for plotting + sns.barplot( + x="strategy", + y="value", + data=confounds_phenotype[confounds_phenotype["type"] == "total"], + ci=95, + color=palette[4], + linewidth=1, + ax=ax, + ) + sns.barplot( + x="strategy", + y="value", + data=confounds_phenotype[confounds_phenotype["type"] == "compcor"], + ci=95, + color=palette[0], + linewidth=1, + ax=ax, + ) + sns.barplot( + x="strategy", + y="value", + data=confounds_phenotype[confounds_phenotype["type"] == "aroma"], + ci=95, + color=palette[2], + linewidth=1, + ax=ax, + ) + sns.barplot( + x="strategy", + y="value", + data=confounds_phenotype[confounds_phenotype["type"] == "fixed_regressors"], + ci=95, + color=palette[1], + linewidth=1, + ax=ax, + ) + sns.barplot( + x="strategy", + y="value", + data=confounds_phenotype[confounds_phenotype["type"] == "high_pass"], + ci=95, + color=palette[3], + linewidth=1, + ax=ax, + ) + ax.set_ylim(0, 100) + ax.set_ylabel(f"% Degrees of freedom loss\n(Full length: {full_length})") + ax.set_title(title) + ax.set_xticklabels( + strategy_order, rotation=45, ha="right", rotation_mode="anchor" + ) + + return confounds_phenotype + + +def plot_stats(confounds_phenotypes, plot_subgroup=False): sns.set_palette("colorblind") palette = sns.color_palette(n_colors=5) - fig = plt.figure(constrained_layout=True, figsize=(11, 5)) - axs = fig.subplots(1, 2, sharey=True) + figsize = (11, 13) if plot_subgroup else (11, 5) + fig = plt.figure(constrained_layout=True, figsize=figsize) + axs = fig.subplots(3, 3, sharey=True) if plot_subgroup else fig.subplots(1, 2, sharey=True) fig.suptitle( "Loss of temporal degrees of freedom", weight="heavy", @@ -40,81 +105,30 @@ def plot_stats(confounds_phenotypes): ) print("Generating new graph") - for ax, dataset in zip(axs, confounds_phenotypes): - confounds_phenotype = confounds_phenotypes[dataset] - _descriptive_stats(dataset, confounds_phenotype) - - # change up the data a bit for plotting - full_length = confounds_phenotype.iloc[0, -1] - confounds_phenotype.loc[:, ("aroma", "aroma")] += confounds_phenotype.loc[ - :, ("aroma", "fixed_regressors") - ] - confounds_phenotype.loc[:, ("compcor", "compcor")] += confounds_phenotype.loc[ - :, ("compcor", "fixed_regressors") - ] - confounds_phenotype.loc[:, ("compcor6", "compcor")] += confounds_phenotype.loc[ - :, ("compcor6", "fixed_regressors") - ] - - confounds_phenotype = confounds_phenotype.reset_index() - confounds_phenotype = confounds_phenotype.melt( - id_vars=["index"], - var_name=["strategy", "type"], - ) - confounds_phenotype["value"] /= full_length - confounds_phenotype["value"] *= 100 - sns.barplot( - x="strategy", - y="value", - data=confounds_phenotype[confounds_phenotype["type"] == "total"], - ci=95, - color=palette[4], - linewidth=1, - ax=ax, - ) - sns.barplot( - x="strategy", - y="value", - data=confounds_phenotype[confounds_phenotype["type"] == "compcor"], - ci=95, - color=palette[0], - linewidth=1, - ax=ax, - ) - sns.barplot( - x="strategy", - y="value", - data=confounds_phenotype[confounds_phenotype["type"] == "aroma"], - ci=95, - color=palette[2], - linewidth=1, - ax=ax, - ) - sns.barplot( - x="strategy", - y="value", - data=confounds_phenotype[confounds_phenotype["type"] == "fixed_regressors"], - ci=95, - color=palette[1], - linewidth=1, - ax=ax, - ) - sns.barplot( - x="strategy", - y="value", - data=confounds_phenotype[confounds_phenotype["type"] == "high_pass"], - ci=95, - color=palette[3], - linewidth=1, - ax=ax, - ) - ax.set_ylim(0, 100) - ax.set_ylabel(f"% Degrees of freedom loss\n(Full length: {full_length})") - ax.set_title(dataset) - ax.set_xticklabels( - strategy_order, rotation=45, ha="right", rotation_mode="anchor" - ) - + key_axs = [axs[0, 0], axs[1, 0]] if plot_subgroup else [axs[0], axs[1]] + for ax, dataset in zip(key_axs, confounds_phenotypes): + confounds_phenotype = confounds_phenotypes[dataset]["confounds_phenotype"] + _descriptive_stats(dataset, confounds_phenotypes[dataset]["confounds_stats"]) + n = confounds_phenotypes[dataset]["participant_labels"].shape[0] + _plot_single_report(confounds_phenotype, + confounds_phenotypes[dataset]["full_length"], + palette, ax, f"{dataset}(N={n})") + if plot_subgroup: + starting_x = 0 + for dataset in confounds_phenotypes: + subjects = confounds_phenotypes[dataset]["participant_labels"].index + for i, group in enumerate(confounds_phenotypes[dataset]["group_values"]): + # return index from participant_labels if group is in participant_labels + selected = confounds_phenotypes[dataset]["participant_labels"] == group + selected = subjects[selected] + confounds_phenotype = confounds_phenotypes[dataset]["confounds_phenotype"].loc[selected, :] + ax = axs[starting_x, (i % 2) + 1] + if i % 2 == 1: + starting_x += 1 + # _descriptive_stats(dataset, confounds_phenotype) + _plot_single_report(confounds_phenotype, + confounds_phenotypes[dataset]["full_length"], + palette, ax, f"{dataset}: {group} (N={selected.shape[0]})") colors = [palette[4], palette[0], palette[2], palette[1], palette[3]] labels = [ "Censored volumes", @@ -124,10 +138,94 @@ def plot_stats(confounds_phenotypes): "Discrete cosine-basis \nregressors", ] handles = [mpatches.Patch(color=c, label=l) for c, l in zip(colors, labels)] - axs[1].legend(handles=handles, loc=0) + if plot_subgroup: + axs[2,0].legend(handles=handles, loc=0) + axs[2,0].axis('off') + else: + axs[1].legend(handles=handles, loc=0) return fig +def _organise_data(confounds_phenotype): + full_length = confounds_phenotype.iloc[0, -1] + confounds_phenotype.loc[:, ("aroma", "aroma")] += confounds_phenotype.loc[ + :, ("aroma", "fixed_regressors") + ] + confounds_phenotype.loc[:, ("compcor", "compcor")] += confounds_phenotype.loc[ + :, ("compcor", "fixed_regressors") + ] + confounds_phenotype.loc[:, ("compcor6", "compcor")] += confounds_phenotype.loc[ + :, ("compcor6", "fixed_regressors") + ] + + confounds_phenotype = confounds_phenotype.reset_index() + confounds_phenotype = confounds_phenotype.melt( + id_vars=["index"], + var_name=["strategy", "type"], + ) + confounds_phenotype["value"] /= full_length + confounds_phenotype["value"] *= 100 + confounds_phenotype = confounds_phenotype.set_index("index") + return confounds_phenotype, full_length + + +def _plot_single_report(confounds_phenotype, full_length, palette, ax, title): + # change up the data a bit for plotting + sns.barplot( + x="strategy", + y="value", + data=confounds_phenotype[confounds_phenotype["type"] == "total"], + ci=95, + color=palette[4], + linewidth=1, + ax=ax, + ) + sns.barplot( + x="strategy", + y="value", + data=confounds_phenotype[confounds_phenotype["type"] == "compcor"], + ci=95, + color=palette[0], + linewidth=1, + ax=ax, + ) + sns.barplot( + x="strategy", + y="value", + data=confounds_phenotype[confounds_phenotype["type"] == "aroma"], + ci=95, + color=palette[2], + linewidth=1, + ax=ax, + ) + sns.barplot( + x="strategy", + y="value", + data=confounds_phenotype[confounds_phenotype["type"] == "fixed_regressors"], + ci=95, + color=palette[1], + linewidth=1, + ax=ax, + ) + sns.barplot( + x="strategy", + y="value", + data=confounds_phenotype[confounds_phenotype["type"] == "high_pass"], + ci=95, + color=palette[3], + linewidth=1, + ax=ax, + ) + ax.set_ylim(0, 100) + ax.set_ylabel(f"% Degrees of freedom loss\n(Full length: {full_length})") + ax.set_title(title) + ax.set_xticklabels( + strategy_order, rotation=45, ha="right", rotation_mode="anchor" + ) + + return confounds_phenotype + + def _descriptive_stats(dataset, confounds_phenotype): print(dataset) _report_descriptive_stats( diff --git a/fmriprep_denoise/visualization/strategy_ranking.py b/fmriprep_denoise/visualization/strategy_ranking.py index a06b63b..9ba5eb7 100644 --- a/fmriprep_denoise/visualization/strategy_ranking.py +++ b/fmriprep_denoise/visualization/strategy_ranking.py @@ -36,7 +36,7 @@ def _rank_degrees_of_freedom(path_root, criteria_name, v, d): """Rank the loss of temporal degrees of freedom.""" dof = degrees_of_freedom_loss.load_data(path_root, [d], criteria_name, v) current_ranking = { - s: [dof[d].loc[:, (s, "total")].mean()] for s in strategy_order + s: [dof[d]["confounds_stats"].loc[:, (s, "total")].mean()] for s in strategy_order } order = pd.DataFrame(current_ranking).T.sort_values(0) order.index = order.index.set_names(["strategy"]) diff --git a/scripts/make_manuscript_figures.py b/scripts/make_manuscript_figures.py index d71c230..1f15d26 100644 --- a/scripts/make_manuscript_figures.py +++ b/scripts/make_manuscript_figures.py @@ -56,7 +56,27 @@ Path(__file__).parents[1] / "outputs" / "loss_degrees_of_freedom.png", transparent=True, ) - + fig_degrees_of_freedom = degrees_of_freedom_loss.plot_stats(data, plot_subgroup=True) + fig_degrees_of_freedom.savefig( + Path(__file__).parents[1] / "outputs" / "loss_degrees_of_freedom_subgroups.png", + transparent=True, + ) + data = degrees_of_freedom_loss.load_data( + path_root, datasets, "minimal", fmriprep_version + ) + fig_degrees_of_freedom = degrees_of_freedom_loss.plot_stats(data, plot_subgroup=True) + fig_degrees_of_freedom.savefig( + Path(__file__).parents[1] / "outputs" / "loss_degrees_of_freedom_subgroups_qc-minimal.png", + transparent=True, + ) + data = degrees_of_freedom_loss.load_data( + path_root, datasets, None, fmriprep_version + ) + fig_degrees_of_freedom = degrees_of_freedom_loss.plot_stats(data, plot_subgroup=True) + fig_degrees_of_freedom.savefig( + Path(__file__).parents[1] / "outputs" / "loss_degrees_of_freedom_subgroups_qc-none.png", + transparent=True, + ) # Plotting metrics metrics = { "p_values": "sig_qcfc",