diff --git a/docs/api.rst b/docs/api.rst index 9d5420fbb..eeaae3430 100644 --- a/docs/api.rst +++ b/docs/api.rst @@ -95,6 +95,7 @@ For more information about the components of coordinate-based meta-analysis in N :template: class.rst diagnostics.Jackknife + diagnostics.FocusCounter .. _api_annotate_ref: diff --git a/examples/02_meta-analyses/05_plot_cbma_subtraction_conjunction.py b/examples/02_meta-analyses/05_plot_cbma_subtraction_conjunction.py index a603b99ce..395413998 100644 --- a/examples/02_meta-analyses/05_plot_cbma_subtraction_conjunction.py +++ b/examples/02_meta-analyses/05_plot_cbma_subtraction_conjunction.py @@ -20,7 +20,7 @@ from nimare import io, utils from nimare.correct import FWECorrector -from nimare.diagnostics import Jackknife +from nimare.diagnostics import FocusCounter, Jackknife from nimare.meta.cbma import ALE, ALESubtraction ############################################################################### @@ -83,21 +83,41 @@ ############################################################################### # Characterize the relative contributions of experiments in the ALE results # ----------------------------------------------------------------------------- +# NiMARE contains two methods for this: :class:`~nimare.diagnostics.Jackknife` +# and :class:`~nimare.diagnostics.FocusCounter`. +# We will show both below. -jknife = Jackknife( +counter = FocusCounter( target_image="z_desc-size_level-cluster_corr-FWE_method-montecarlo", voxel_thresh=None, ) -knowledge_cluster_table, knowledge_cluster_img = jknife.transform(knowledge_corrected_results) -related_cluster_table, related_cluster_img = jknife.transform(related_corrected_results) +knowledge_count_table, knowledge_cluster_img = counter.transform(knowledge_corrected_results) +related_count_table, related_cluster_img = counter.transform(related_corrected_results) # %% # ############################################################################# -knowledge_cluster_table.head(10) +knowledge_count_table.head(10) # %% # ############################################################################# -related_cluster_table.head(10) +related_count_table.head(10) + +# %% +# ############################################################################# +jackknife = Jackknife( + target_image="z_desc-size_level-cluster_corr-FWE_method-montecarlo", + voxel_thresh=None, +) +knowledge_jackknife_table, _ = jackknife.transform(knowledge_corrected_results) +related_jackknife_table, _ = jackknife.transform(related_corrected_results) + +# %% +# ############################################################################# +knowledge_jackknife_table.head(10) + +# %% +# ############################################################################# +related_jackknife_table.head(10) ############################################################################### # Subtraction analysis diff --git a/examples/02_meta-analyses/08_plot_simulated_data.py b/examples/02_meta-analyses/08_plot_simulated_data.py index 0b638f3b6..d0873778d 100644 --- a/examples/02_meta-analyses/08_plot_simulated_data.py +++ b/examples/02_meta-analyses/08_plot_simulated_data.py @@ -110,7 +110,6 @@ def analyze_and_plot(dset, ground_truth_foci=None, correct=True, return_cres=Fal # ------------------------------- display = analyze_and_plot(manual_dset, ground_truth_foci) -display.frame_axes.figure.show() ############################################################################### # Control percentage of studies with the foci of interest @@ -129,7 +128,6 @@ def analyze_and_plot(dset, ground_truth_foci=None, correct=True, return_cres=Fal # ------------------------------------- display = analyze_and_plot(perc_foci_dset, ground_truth_foci[0:2]) -display.frame_axes.figure.show() ############################################################################### # Create a null dataset diff --git a/nimare/diagnostics.py b/nimare/diagnostics.py index 7347e1982..945726379 100644 --- a/nimare/diagnostics.py +++ b/nimare/diagnostics.py @@ -9,10 +9,11 @@ from nibabel.funcs import squeeze_image from nilearn import image, input_data from scipy import ndimage +from scipy.spatial.distance import cdist from tqdm.auto import tqdm -from .base import NiMAREBase -from .utils import tqdm_joblib, vox2mm +from nimare.base import NiMAREBase +from nimare.utils import mm2vox, tqdm_joblib, vox2mm LGR = logging.getLogger(__name__) @@ -20,6 +21,8 @@ class Jackknife(NiMAREBase): """Run a jackknife analysis on a meta-analysis result. + .. versionadded:: 0.0.11 + Parameters ---------- target_image : :obj:`str`, optional @@ -205,6 +208,8 @@ def _transform( original_masker, cluster_masker, ): + estimator = copy.deepcopy(estimator) + # Fit Estimator to all studies except the target study other_ids = [id_ for id_ in all_ids if id_ != expid] temp_dset = dset.slice(other_ids) @@ -232,3 +237,179 @@ def _transform( stat_prop_values = cluster_masker.transform(stat_prop_img) return expid, stat_prop_values + + +class FocusCounter(NiMAREBase): + """Run a focus-count analysis on a coordinate-based meta-analysis result. + + .. versionadded:: 0.0.12 + + Parameters + ---------- + target_image : :obj:`str`, optional + The meta-analytic map for which clusters will be characterized. + The default is z because log-p will not always have value of zero for non-cluster voxels. + voxel_thresh : :obj:`float` or None, optional + An optional voxel-level threshold that may be applied to the ``target_image`` to define + clusters. This can be None if the ``target_image`` is already thresholded + (e.g., a cluster-level corrected map). + Default is None. + n_cores : :obj:`int`, optional + Number of cores to use for parallelization. + If <=0, defaults to using all available cores. + Default is 1. + + Notes + ----- + This analysis characterizes the relative contribution of each experiment in a meta-analysis + to the resulting clusters by counting the number of peaks from each experiment that fall within + each significant cluster. + + Warnings + -------- + This method only works for coordinate-based meta-analyses. + + Pairwise meta-analyses, like ALESubtraction and MKDAChi2, are not yet supported in this + method. + """ + + def __init__( + self, + target_image="z_desc-size_level-cluster_corr-FWE_method-montecarlo", + voxel_thresh=None, + n_cores=1, + ): + self.target_image = target_image + self.voxel_thresh = voxel_thresh + self.n_cores = self._check_ncores(n_cores) + + def transform(self, result): + """Apply the analysis to a MetaResult. + + Parameters + ---------- + result : :obj:`~nimare.results.MetaResult` + A MetaResult produced by a coordinate- or image-based meta-analysis. + + Returns + ------- + contribution_table : :obj:`pandas.DataFrame` + A DataFrame with information about relative contributions of each experiment to each + cluster in the thresholded map. + There is one row for each experiment, as well as one more row at the top of the table + (below the header), which has the center of mass of each cluster. + The centers of mass are not guaranteed to fall within the actual clusters, but can + serve as a useful heuristic for identifying them. + There is one column for each cluster, with column names being integers indicating the + cluster's associated value in the ``labeled_cluster_img`` output. + labeled_cluster_img : :obj:`nibabel.nifti1.Nifti1Image` + The labeled, thresholded map that is used to identify clusters characterized by this + analysis. + Each cluster in the map has a single value, which corresponds to the cluster's column + name in ``contribution_table``. + """ + if not hasattr(result.estimator, "dataset"): + raise AttributeError( + "MetaResult was not generated by an Estimator with a `dataset` attribute. " + "This may be because the Estimator was a pairwise Estimator. " + "The Jackknife method does not currently work with pairwise Estimators." + ) + + # We need to copy the estimator because it will otherwise overwrite the original version + # with one missing a study in its inputs. + estimator = copy.deepcopy(result.estimator) + + # Collect the thresholded cluster map + if self.target_image in result.maps: + target_img = result.get_map(self.target_image, return_type="image") + else: + available_maps = [f"'{m}'" for m in result.maps.keys()] + raise ValueError( + f"Target image ('{self.target_image}') not present in result. " + f"Available maps in result are: {', '.join(available_maps)}." + ) + + if self.voxel_thresh: + thresh_img = image.threshold_img(target_img, self.voxel_thresh) + else: + thresh_img = target_img + + thresh_arr = thresh_img.get_fdata() + + # Use study IDs in inputs_ instead of dataset, because we don't want to try fitting the + # estimator to a study that might have been filtered out by the estimator's criteria. + meta_ids = estimator.inputs_["id"] + rows = ["Center of Mass"] + list(meta_ids) + + # Let's label the clusters in the thresholded map so we can use it as a NiftiLabelsMasker + # This won't work when the Estimator's masker isn't a NiftiMasker... :( + conn = np.zeros((3, 3, 3), int) + conn[:, :, 1] = 1 + conn[:, 1, :] = 1 + conn[1, :, :] = 1 + labeled_cluster_arr, n_clusters = ndimage.measurements.label(thresh_arr, conn) + labeled_cluster_img = nib.Nifti1Image( + labeled_cluster_arr, + affine=target_img.affine, + header=target_img.header, + ) + + if n_clusters == 0: + LGR.warning("No clusters found") + contribution_table = pd.DataFrame(index=rows) + return contribution_table, labeled_cluster_img + + # Identify center of mass for each cluster + # This COM may fall outside the cluster, but it is a useful heuristic for identifying them + cluster_ids = list(range(1, n_clusters + 1)) + cluster_coms = ndimage.center_of_mass( + labeled_cluster_arr, labeled_cluster_arr, cluster_ids + ) + cluster_coms = np.array(cluster_coms) + cluster_coms = vox2mm(cluster_coms, target_img.affine) + + cluster_com_strs = [] + for i_peak in range(len(cluster_ids)): + x, y, z = cluster_coms[i_peak, :].astype(int) + xyz_str = f"({x}, {y}, {z})" + cluster_com_strs.append(xyz_str) + + # Create empty contribution table + contribution_table = pd.DataFrame(index=rows, columns=cluster_ids) + contribution_table.index.name = "Cluster ID" + contribution_table.loc["Center of Mass"] = cluster_com_strs + + with tqdm_joblib(tqdm(total=len(meta_ids))): + jackknife_results = Parallel(n_jobs=self.n_cores)( + delayed(self._transform)( + study_id, + coordinates_df=estimator.inputs_["coordinates"], + labeled_cluster_map=labeled_cluster_arr, + affine=target_img.affine, + ) + for study_id in meta_ids + ) + + # Add the results to the table + for expid, focus_counts in jackknife_results: + contribution_table.loc[expid] = focus_counts + + return contribution_table, labeled_cluster_img + + def _transform(self, expid, coordinates_df, labeled_cluster_map, affine): + coords = coordinates_df.loc[coordinates_df["id"] == expid] + ijk = mm2vox(coords[["x", "y", "z"]], affine) + + clust_ids = sorted(list(np.unique(labeled_cluster_map)[1:])) + focus_counts = [] + + for i_cluster, c_val in enumerate(clust_ids): + cluster_mask = labeled_cluster_map == c_val + cluster_idx = np.vstack(np.where(cluster_mask)) + distances = cdist(cluster_idx.T, ijk) + distances = distances < 1 + distances = np.any(distances, axis=0) + n_included_voxels = np.sum(distances) + focus_counts.append(n_included_voxels) + + return expid, focus_counts diff --git a/nimare/tests/test_diagnostics.py b/nimare/tests/test_diagnostics.py index c0108b9e6..3502cf0d9 100644 --- a/nimare/tests/test_diagnostics.py +++ b/nimare/tests/test_diagnostics.py @@ -4,7 +4,7 @@ import pytest from nilearn.input_data import NiftiLabelsMasker -from nimare.diagnostics import Jackknife +from nimare.diagnostics import FocusCounter, Jackknife from nimare.meta import cbma, ibma from nimare.tests.utils import get_test_data_path @@ -72,3 +72,38 @@ def test_jackknife_with_custom_masker_smoke(testdata_ibma): with pytest.raises(ValueError): jackknife = Jackknife(target_image="doggy", voxel_thresh=0.5) jackknife.transform(res) + + +@pytest.mark.parametrize( + "estimator,meta_type,n_samples,target_image", + [ + (cbma.ALE, "cbma", "onesample", "z"), + (cbma.MKDADensity, "cbma", "onesample", "z"), + (cbma.KDA, "cbma", "onesample", "z"), + (cbma.MKDAChi2, "cbma", "twosample", "z_desc-consistency"), + ], +) +def test_focuscounter_smoke( + testdata_ibma, + testdata_cbma_full, + estimator, + meta_type, + n_samples, + target_image, +): + """Smoke test the FocusCounter method.""" + meta = estimator() + testdata = testdata_ibma if meta_type == "ibma" else testdata_cbma_full + if n_samples == "twosample": + res = meta.fit(testdata, testdata) + else: + res = meta.fit(testdata) + + counter = FocusCounter(target_image=target_image, voxel_thresh=1.65) + + if n_samples == "twosample": + with pytest.raises(AttributeError): + counter.transform(res) + else: + cluster_table, labeled_img = counter.transform(res) + assert cluster_table.shape[0] == len(meta.inputs_["id"]) + 1 diff --git a/nimare/workflows/ale.py b/nimare/workflows/ale.py index cb405b97c..acc572d0a 100644 --- a/nimare/workflows/ale.py +++ b/nimare/workflows/ale.py @@ -6,9 +6,10 @@ import numpy as np -from ..correct import FWECorrector -from ..io import convert_sleuth_to_dataset -from ..meta import ALE, ALESubtraction +from nimare.correct import FWECorrector +from nimare.diagnostics import FocusCounter +from nimare.io import convert_sleuth_to_dataset +from nimare.meta import ALE, ALESubtraction LGR = logging.getLogger(__name__) @@ -72,7 +73,7 @@ def ale_sleuth_workflow( ---------- - Eickhoff, S. B., Bzdok, D., Laird, A. R., Kurth, F., & Fox, P. T. (2012). Activation likelihood estimation meta-analysis revisited. NeuroImage, -59(3), 2349–2361. +59(3), 2349-2361. - Fonov, V., Evans, A. C., Botteron, K., Almli, C. R., McKinstry, R. C., Collins, D. L., & Brain Development Cooperative Group. (2011). Unbiased average age-appropriate atlases for pediatric studies. @@ -82,11 +83,11 @@ def ale_sleuth_workflow( to adulthood. NeuroImage, (47), S102. - Turkeltaub, P. E., Eden, G. F., Jones, K. M., & Zeffiro, T. A. (2002). Meta-analysis of the functional neuroanatomy of single-word reading: method -and validation. NeuroImage, 16(3 Pt 1), 765–780. +and validation. NeuroImage, 16(3 Pt 1), 765-780. - Turkeltaub, P. E., Eickhoff, S. B., Laird, A. R., Fox, M., Wiener, M., & Fox, P. (2012). Minimizing within-experiment and within-group effects in Activation Likelihood Estimation meta-analyses. Human Brain Mapping, -33(1), 1–13. +33(1), 1-13. """ ale = ALE(kernel__fwhm=fwhm) @@ -97,6 +98,11 @@ def ale_sleuth_workflow( method="montecarlo", n_iters=n_iters, voxel_thresh=v_thr, n_cores=n_cores ) cres = corr.transform(results) + fcounter = FocusCounter( + target_image="z_desc-size_level-cluster_corr-FWE_method-montecarlo", + voxel_thresh=None, + ) + count_df, _ = fcounter.transform(cres) boilerplate = boilerplate.format( n_exps=len(dset.ids), @@ -149,14 +155,14 @@ def ale_sleuth_workflow( ---------- - Turkeltaub, P. E., Eden, G. F., Jones, K. M., & Zeffiro, T. A. (2002). Meta-analysis of the functional neuroanatomy of single-word reading: method -and validation. NeuroImage, 16(3 Pt 1), 765–780. +and validation. NeuroImage, 16(3 Pt 1), 765-780. - Eickhoff, S. B., Bzdok, D., Laird, A. R., Kurth, F., & Fox, P. T. (2012). Activation likelihood estimation meta-analysis revisited. NeuroImage, -59(3), 2349–2361. +59(3), 2349-2361. - Turkeltaub, P. E., Eickhoff, S. B., Laird, A. R., Fox, M., Wiener, M., & Fox, P. (2012). Minimizing within-experiment and within-group effects in Activation Likelihood Estimation meta-analyses. Human Brain Mapping, -33(1), 1–13. +33(1), 1-13. - Fonov, V., Evans, A. C., Botteron, K., Almli, C. R., McKinstry, R. C., Collins, D. L., & Brain Development Cooperative Group. (2011). Unbiased average age-appropriate atlases for pediatric studies. @@ -180,7 +186,15 @@ def ale_sleuth_workflow( method="montecarlo", n_iters=n_iters, voxel_thresh=v_thr, n_cores=n_cores ) cres1 = corr.transform(res1) + fcounter = FocusCounter( + target_image="z_desc-size_level-cluster_corr-FWE_method-montecarlo", + voxel_thresh=None, + ) + count_df1, _ = fcounter.transform(cres1) + cres2 = corr.transform(res2) + count_df2, _ = fcounter.transform(cres2) + sub = ALESubtraction(n_iters=n_iters, kernel__fwhm=fwhm) sres = sub.fit(dset1, dset2) @@ -211,13 +225,16 @@ def ale_sleuth_workflow( LGR.info("Saving output maps...") if not sleuth_file2: cres.save_maps(output_dir=output_dir, prefix=prefix) + count_df.to_csv(os.path.join(output_dir, prefix + "_clust.tsv"), index=False, sep="\t") copyfile(sleuth_file, os.path.join(output_dir, prefix + "input_coordinates.txt")) else: prefix1 = os.path.splitext(os.path.basename(sleuth_file))[0] + "_" prefix2 = os.path.splitext(os.path.basename(sleuth_file2))[0] + "_" prefix3 = prefix + "subtraction_" cres1.save_maps(output_dir=output_dir, prefix=prefix1) + count_df1.to_csv(os.path.join(output_dir, prefix1 + "_clust.tsv"), index=False, sep="\t") cres2.save_maps(output_dir=output_dir, prefix=prefix2) + count_df2.to_csv(os.path.join(output_dir, prefix2 + "_clust.tsv"), index=False, sep="\t") sres.save_maps(output_dir=output_dir, prefix=prefix3) copyfile(sleuth_file, os.path.join(output_dir, prefix + "group1_input_coordinates.txt")) copyfile(sleuth_file2, os.path.join(output_dir, prefix + "group2_input_coordinates.txt"))