Repository: neurostuff/NiMARE

Issue:
Other annotation issues
## Summary
<!--What is the nature of the bug?-->

Moving this here so as not to clutter up neurostars with a bunch of bugs. I currently run into two issues:
1. If I add the sample sizes as described on neurostars, slicing the studyset does not retain that information and actually appears to lose some of the information that is present in the annotations and metadata.
2. If I go ahead an just run the ALE on the whole dataset (as proof of concept) the following happens
 
## Additional details
<!--Please fill in the following details-->
- NiMARE version: 0.4.2

### What were you trying to do?

### What did you expect to happen?

### What actually happened?

## Reproducing the bug
<!--Please share any steps you performed that revealed the bug-->
<!--Please include any code snippets.
Enclose them in triple back-ticks (```)
Like this:

```
<code>
```
-->
```
from importlib import import_module
from pathlib import Path
import gzip
import json
import io
from requests import request
from nimare.io import convert_nimads_to_dataset
from nimare.nimads import Studyset 

def download_file(url):
    """Download a file from NeuroStore."""
    response = request("GET", url)
    return response.json()

nimads_studyset = download_file("https://neurostore.org/api/studysets/8mMthr2GDTaC?nested=true")
nimads_annotation = download_file("https://neurostore.org/api/annotations/MBdEAkukAz7V")

from nimare.io import convert_nimads_to_dataset
from nimare.nimads import Studyset, Annotation
studyset = Studyset(nimads_studyset, annotations = nimads_annotation)
annotation = Annotation(nimads_annotation, studyset)

for study in studyset.studies:
    for analysis in study.analyses:
        analysis_annotation = list(analysis.annotations.values())[0]
        if analysis_annotation['included']:
            sample_sizes = [
                analysis_annotation['patient_samplesize'],
                analysis_annotation['control_samplesize'],
            ]
        else:
            sample_sizes = [10,10]
        print(sample_sizes)
        analysis.metadata['sample_sizes'] = sample_sizes
nimare_dset = studyset.to_dataset()
nimare_dset.coordinates.head()

from nimare.meta.cbma.ale import ALE
from nimare.correct import FWECorrector

out_dir = '../results/a2nd/2025-03-11_ALE/hypermetabolism'

ale = ALE()
results = ale.fit(nimare_dset)
corr = FWECorrector(
    method="montecarlo", n_iters=100, voxel_thresh=0.001, n_cores=12
)
cres = corr.transform(results)
cres.save_maps(output_dir = out_dir, prefix = "ALE_FWE")
```

Error message:
```
/opt/conda/lib/python3.12/site-packages/nimare/utils.py:210: RuntimeWarning: invalid value encountered in cast
  ijk = nib.affines.apply_affine(np.linalg.inv(affine), xyz).astype(int)

---------------------------------------------------------------------------
UFuncTypeError                            Traceback (most recent call last)
Cell In[138], line 4
      1 out_dir = '../results/a2nd/2025-03-11_ALE/hypermetabolism'
      3 ale = ALE()
----> 4 results = ale.fit(nimare_dset)
      5 corr = FWECorrector(
      6     method="montecarlo", n_iters=100, voxel_thresh=0.001, n_cores=12
      7 )
      8 cres = corr.transform(results)

File /opt/conda/lib/python3.12/site-packages/nimare/estimator.py:131, in Estimator.fit(self, dataset, drop_invalid)
    103 """Fit Estimator to Dataset.
    104 
    105 Parameters
   (...)
    128 call `fit`.
    129 """
    130 self._collect_inputs(dataset, drop_invalid=drop_invalid)
--> 131 self._preprocess_input(dataset)
    132 maps, tables, description = self._cache(self._fit, func_memory_level=1)(dataset)
    134 if hasattr(self, "masker") and self.masker is not None:

File /opt/conda/lib/python3.12/site-packages/nimare/meta/cbma/base.py:144, in CBMAEstimator._preprocess_input(self, dataset)
    141 # Integrate "sample_size" from metadata into DataFrame so that
    142 # kernel_transformer can access it.
    143 if "sample_size" in kt_args:
--> 144     self.inputs_["coordinates"] = _add_metadata_to_dataframe(
    145         dataset,
    146         self.inputs_["coordinates"],
    147         metadata_field="sample_sizes",
    148         target_column="sample_size",
    149         filter_func=np.mean,
    150     )

File /opt/conda/lib/python3.12/site-packages/nimare/utils.py:860, in _add_metadata_to_dataframe(dataset, dataframe, metadata_field, target_column, filter_func)
    854 metadata = pd.DataFrame(
    855     index=dataset.ids,
    856     data=metadata,
    857     columns=[metadata_field],
    858 )
    859 # Reduce the metadata (if in list/array format) to single values
--> 860 metadata[target_column] = metadata[metadata_field].apply(
    861     lambda x: None if x is None else filter_func(x)
    862 )
    863 # Merge metadata df into coordinates df
    864 dataframe = dataframe.merge(
    865     right=metadata,
    866     left_on="id",
   (...)
    871     how="left",
    872 )

File /opt/conda/lib/python3.12/site-packages/pandas/core/series.py:4924, in Series.apply(self, func, convert_dtype, args, by_row, **kwargs)
   4789 def apply(
   4790     self,
   4791     func: AggFuncType,
   (...)
   4796     **kwargs,
   4797 ) -> DataFrame | Series:
   4798     """
   4799     Invoke function on values of Series.
   4800 
   (...)
   4915     dtype: float64
   4916     """
   4917     return SeriesApply(
   4918         self,
   4919         func,
   4920         convert_dtype=convert_dtype,
   4921         by_row=by_row,
   4922         args=args,
   4923         kwargs=kwargs,
-> 4924     ).apply()

File /opt/conda/lib/python3.12/site-packages/pandas/core/apply.py:1427, in SeriesApply.apply(self)
   1424     return self.apply_compat()
   1426 # self.func is Callable
-> 1427 return self.apply_standard()

File /opt/conda/lib/python3.12/site-packages/pandas/core/apply.py:1507, in SeriesApply.apply_standard(self)
   1501 # row-wise access
   1502 # apply doesn't have a `na_action` keyword and for backward compat reasons
   1503 # we need to give `na_action="ignore"` for categorical data.
   1504 # TODO: remove the `na_action="ignore"` when that default has been changed in
   1505 #  Categorical (GH51645).
   1506 action = "ignore" if isinstance(obj.dtype, CategoricalDtype) else None
-> 1507 mapped = obj._map_values(
   1508     mapper=curried, na_action=action, convert=self.convert_dtype
   1509 )
   1511 if len(mapped) and isinstance(mapped[0], ABCSeries):
   1512     # GH#43986 Need to do list(mapped) in order to get treated as nested
   1513     #  See also GH#25959 regarding EA support
   1514     return obj._constructor_expanddim(list(mapped), index=obj.index)

File /opt/conda/lib/python3.12/site-packages/pandas/core/base.py:921, in IndexOpsMixin._map_values(self, mapper, na_action, convert)
    918 if isinstance(arr, ExtensionArray):
    919     return arr.map(mapper, na_action=na_action)
--> 921 return algorithms.map_array(arr, mapper, na_action=na_action, convert=convert)

File /opt/conda/lib/python3.12/site-packages/pandas/core/algorithms.py:1743, in map_array(arr, mapper, na_action, convert)
   1741 values = arr.astype(object, copy=False)
   1742 if na_action is None:
-> 1743     return lib.map_infer(values, mapper, convert=convert)
   1744 else:
   1745     return lib.map_infer_mask(
   1746         values, mapper, mask=isna(values).view(np.uint8), convert=convert
   1747     )

File lib.pyx:2972, in pandas._libs.lib.map_infer()

File /opt/conda/lib/python3.12/site-packages/nimare/utils.py:861, in _add_metadata_to_dataframe.<locals>.<lambda>(x)
    854 metadata = pd.DataFrame(
    855     index=dataset.ids,
    856     data=metadata,
    857     columns=[metadata_field],
    858 )
    859 # Reduce the metadata (if in list/array format) to single values
    860 metadata[target_column] = metadata[metadata_field].apply(
--> 861     lambda x: None if x is None else filter_func(x)
    862 )
    863 # Merge metadata df into coordinates df
    864 dataframe = dataframe.merge(
    865     right=metadata,
    866     left_on="id",
   (...)
    871     how="left",
    872 )

File /opt/conda/lib/python3.12/site-packages/numpy/core/fromnumeric.py:3504, in mean(a, axis, dtype, out, keepdims, where)
   3501     else:
   3502         return mean(axis=axis, dtype=dtype, out=out, **kwargs)
-> 3504 return _methods._mean(a, axis=axis, dtype=dtype,
   3505                       out=out, **kwargs)

File /opt/conda/lib/python3.12/site-packages/numpy/core/_methods.py:118, in _mean(a, axis, dtype, out, keepdims, where)
    115         dtype = mu.dtype('f4')
    116         is_float16_result = True
--> 118 ret = umr_sum(arr, axis, dtype, out, keepdims, where=where)
    119 if isinstance(ret, mu.ndarray):
    120     with _no_nep50_warning():

UFuncTypeError: ufunc 'add' did not contain a loop with signature matching types (dtype('<U2'), dtype('<U2')) -> None
```

Hints from issue comments:
@araikes: Mmm... solved it potentially...

```
for study in studyset.studies:
    for analysis in study.analyses:
        analysis_annotation = list(analysis.annotations.values())[0]
        if analysis_annotation['included']:
            analysis_ids.append(analysis.id)
            sample_sizes = [
                int(analysis_annotation['patient_samplesize']),
                int(analysis_annotation['control_samplesize']),
            ]
        else:
            sample_sizes = [10,10]
        print(sample_sizes)
        analysis.metadata['sample_sizes'] = sample_sizes
```

Annotations grab the sample sizes as strings, rather than `ints`.

Doesn't solve the slicing issue@jdkent: Thanks for reporting this!

yes, you are correct that annotations on the platform are either treated as strings or boolean values currently, we can expand the possible value types on the frontend.

Solution patch:
diff --git a/examples/01_datasets/05_plot_nimads.py b/examples/01_datasets/05_plot_nimads.py
index ddc19b5..264d518 100644
--- a/examples/01_datasets/05_plot_nimads.py
+++ b/examples/01_datasets/05_plot_nimads.py
@@ -6,14 +6,16 @@
 Using NIMADS with NiMARE
 ========================
 
-How to use the NeuroImaging Meta-Analysis Data Structure
-`(NIMADS) <https://neurostuff.github.io/NIMADS/>`_ with NiMARE.
+This example demonstrates the key functionality of the NeuroImaging Meta-Analysis Data Structure
+(NIMADS) with NiMARE, including working with StudySets, annotations, coordinates, and metadata.
 """
+
+from pprint import pprint
 from requests import request
 
-from nimare.io import convert_nimads_to_dataset
 from nimare.nimads import Studyset
 
+
 ###############################################################################
 # Download Data from NeuroStore
 # -----------------------------------------------------------------------------
@@ -25,32 +27,103 @@ def download_file(url):
     return response.json()
 
 
+# Download a studyset and its annotation
 nimads_studyset = download_file("https://neurostore.org/api/studysets/Cv2LLUqG76W9?nested=true")
 nimads_annotation = download_file("https://neurostore.org/api/annotations/76PyNqoTNEsE")
 
 
 ###############################################################################
-# Load Data
+# Create and Explore StudySet
 # -----------------------------------------------------------------------------
-# Load the json files into a NiMADS Studyset object.
+# Load the data into a NiMADS Studyset object and explore its contents
 
 studyset = Studyset(nimads_studyset, annotations=nimads_annotation)
 
+# Display basic information about the studyset
+print("\nStudySet Information:")
+print("-" * 50)
+print(f"ID: {studyset.id}")
+print(f"Name: {studyset.name}")
+print(f"Number of studies: {len(studyset.studies)}")
+print(f"Number of annotations: {len(studyset.annotations)}")
+
 
 ###############################################################################
-# Convert to NiMARE Dataset
+# Explore Studies and Analyses
 # -----------------------------------------------------------------------------
-# Convert the NiMADS Studyset object to a NiMARE Dataset object.
-# Then you can run NiMARE analyses on the Dataset object.
+# Look at the first study and its analyses in detail
+
+first_study = studyset.studies[0]
+print("\nFirst Study Details:")
+print("-" * 50)
+print(f"Study ID: {first_study.id}")
+print(f"Title: {first_study.name}")
+print(f"Authors: {first_study.authors}")
+print(f"Publication: {first_study.publication}")
+print(f"Number of analyses: {len(first_study.analyses)}")
+
+# Show details of the first analysis
+first_analysis = first_study.analyses[0]
+print("\nFirst Analysis Details:")
+print("-" * 50)
+print(f"Analysis ID: {first_analysis.id}")
+print(f"Analysis Name: {first_analysis.name}")
+print(f"Number of coordinates: {len(first_analysis.points)}")
+print(f"Number of conditions: {len(first_analysis.conditions)}")
+
+
+###############################################################################
+# Working with Coordinates
+# -----------------------------------------------------------------------------
+# Demonstrate coordinate-based queries
+
+# Example coordinate in MNI space
+example_coord = [-42, -58, -15]  # MNI coordinates
+print("\nCoordinate Search Results:")
+print("-" * 50)
+print(f"Searching near coordinate: {example_coord}")
+
+# Find analyses with coordinates within 10mm
+nearby_analyses = studyset.get_analyses_by_coordinates(example_coord, r=10)
+print(f"\nFound {len(nearby_analyses)} analyses within 10mm")
+
+# Find 5 closest analyses
+closest_analyses = studyset.get_analyses_by_coordinates(example_coord, n=5)
+print(f"\nClosest 5 analyses: {closest_analyses}")
 
-nimare_dset = studyset.to_dataset()
-nimare_dset.coordinates.head()
 
 ###############################################################################
-# Directly to NiMARE Dataset
+# Working with Annotations
 # -----------------------------------------------------------------------------
-# Alternatively, you can convert the NiMADS json files directly to a NiMARE Dataset object
-# if you wish to skip using the nimads studyset object directly.
+# Demonstrate how to work with study annotations
 
-nimare_dset_2 = convert_nimads_to_dataset(nimads_studyset, nimads_annotation)
-nimare_dset_2.coordinates.head()
+print("\nAnnotation Information:")
+print("-" * 50)
+for annotation in studyset.annotations:
+    print(f"\nAnnotation ID: {annotation.id}")
+    print(f"Annotation Name: {annotation.name}")
+    print(f"Number of notes: {len(annotation.notes)}")
+
+
+###############################################################################
+# Query Metadata
+# -----------------------------------------------------------------------------
+# Show how to query analyses based on metadata
+
+# Get all analyses that have a specific metadata field
+metadata_results = studyset.get_analyses_by_metadata("contrast_type")
+print("\nAnalyses with contrast_type metadata:")
+print("-" * 50)
+pprint(metadata_results)
+
+
+###############################################################################
+# Convert to NiMARE Dataset
+# -----------------------------------------------------------------------------
+# Convert the NiMADS Studyset to a NiMARE Dataset for further analysis
+
+nimare_dset = studyset.to_dataset()
+print("\nNiMARE Dataset Information:")
+print("-" * 50)
+print("Coordinates DataFrame Preview:")
+print(nimare_dset.coordinates.head())
diff --git a/nimare/nimads.py b/nimare/nimads.py
index 654c27a..13fd2c6 100644
--- a/nimare/nimads.py
+++ b/nimare/nimads.py
@@ -4,7 +4,11 @@ import json
 import weakref
 from copy import deepcopy
 
+import numpy as np
+from nilearn._utils import load_niimg
+
 from nimare.io import convert_nimads_to_dataset
+from nimare.utils import mm2vox
 
 
 class Studyset:
@@ -84,7 +88,7 @@ class Studyset:
         for study in studyset.studies:
             if len(study.analyses) > 1:
                 source_lst = [analysis.to_dict() for analysis in study.analyses]
-                ids, names, conditions, images, points, weights = [
+                ids, names, conditions, images, points, weights, metadata = [
                     [source[key] for source in source_lst] for key in source_lst[0]
                 ]
 
@@ -95,6 +99,7 @@ class Studyset:
                     "images": [image for i_list in images for image in i_list],
                     "points": [point for p_list in points for point in p_list],
                     "weights": [weight for w_list in weights for weight in w_list],
+                    "metadata": {k: v for m_dict in metadata for k, v in m_dict.items()},
                 }
                 study.analyses = [Analysis(new_source)]
 
@@ -118,12 +123,42 @@ class Studyset:
         return convert_nimads_to_dataset(self)
 
     def load(self, filename):
-        """Load a Studyset from a pickled file."""
-        raise NotImplementedError("Loading from pickled files is not yet supported.")
+        """Load a Studyset from a pickled file.
+
+        Parameters
+        ----------
+        filename : str
+            Path to the pickled file to load from.
+
+        Returns
+        -------
+        Studyset
+            The loaded Studyset object.
+        """
+        import pickle
+
+        with open(filename, "rb") as f:
+            loaded_data = pickle.load(f)
+
+        # Update current instance with loaded data
+        self.id = loaded_data.id
+        self.name = loaded_data.name
+        self.studies = loaded_data.studies
+        self._annotations = loaded_data._annotations
+        return self
 
     def save(self, filename):
-        """Write the Studyset to a pickled file."""
-        raise NotImplementedError("Saving to pickled files is not yet supported.")
+        """Write the Studyset to a pickled file.
+
+        Parameters
+        ----------
+        filename : str
+            Path where the pickled file should be saved.
+        """
+        import pickle
+
+        with open(filename, "wb") as f:
+            pickle.dump(self, f)
 
     def copy(self):
         """Create a copy of the Studyset."""
@@ -141,49 +176,221 @@ class Studyset:
 
         for annot in annotations:
             annot["notes"] = [n for n in annot["notes"] if n["analysis"] in analyses]
-            studyset.annotation = annot
+            studyset.annotations = annot
 
         return studyset
 
     def merge(self, right):
-        """Merge a separate Studyset into the current one."""
-        raise NotImplementedError("Merging Studysets is not yet supported.")
-
-    def update_image_path(self, new_path):
-        """Point to a new location for image files on the local filesystem."""
-        raise NotImplementedError("Updating image paths is not yet supported.")
+        """Merge a separate Studyset into the current one.
+
+        Parameters
+        ----------
+        right : Studyset
+            The other Studyset to merge with this one.
+
+        Returns
+        -------
+        Studyset
+            A new Studyset containing merged studies from both input Studysets.
+            For studies with the same ID, their analyses and metadata are combined,
+            with data from self (left) taking precedence in case of conflicts.
+        """
+        if not isinstance(right, Studyset):
+            raise ValueError("Can only merge with another Studyset")
+
+        # Create new source dictionary starting with left (self) studyset
+        merged_source = self.to_dict()
+        merged_source["id"] = f"{self.id}_{right.id}"
+        merged_source["name"] = f"Merged: {self.name} + {right.name}"
+
+        # Create lookup of existing studies by ID
+        left_studies = {study["id"]: study for study in merged_source["studies"]}
+
+        # Process studies from right studyset
+        right_dict = right.to_dict()
+        for right_study in right_dict["studies"]:
+            study_id = right_study["id"]
+
+            if study_id in left_studies:
+                # Merge study data
+                left_study = left_studies[study_id]
+
+                # Keep metadata from left unless missing
+                left_study["metadata"].update(
+                    {
+                        k: v
+                        for k, v in right_study["metadata"].items()
+                        if k not in left_study["metadata"]
+                    }
+                )
+
+                # Keep basic info from left unless empty
+                for field in ["name", "authors", "publication"]:
+                    if not left_study[field]:
+                        left_study[field] = right_study[field]
+
+                # Combine analyses, avoiding duplicates by ID
+                left_analyses = {a["id"]: a for a in left_study["analyses"]}
+                for right_analysis in right_study["analyses"]:
+                    if right_analysis["id"] not in left_analyses:
+                        left_study["analyses"].append(right_analysis)
+            else:
+                # Add new study
+                merged_source["studies"].append(right_study)
+
+        # Create new merged studyset
+        merged = self.__class__(source=merged_source)
+
+        # Merge annotations, preferring left's annotations for conflicts
+        existing_annot_ids = {a.id for a in self.annotations}
+        for right_annot in right.annotations:
+            if right_annot.id not in existing_annot_ids:
+                merged.annotations = right_annot.to_dict()
+
+        return merged
 
     def get_analyses_by_coordinates(self, xyz, r=None, n=None):
-        """Extract a list of Analyses with at least one Point near the requested coordinates."""
-        raise NotImplementedError("Getting analyses by coordinates is not yet supported.")
+        """Extract a list of Analyses with at least one Point near the requested coordinates.
+
+        Parameters
+        ----------
+        xyz : array_like
+            1 x 3 array of coordinates in mm space to search from
+        r : float, optional
+            Search radius in millimeters.
+            Mutually exclusive with n.
+        n : int, optional
+            Number of closest analyses to return.
+            Mutually exclusive with r.
+
+        Returns
+        -------
+        list[str]
+            A list of Analysis IDs with at least one point within the search criteria.
+
+        Notes
+        -----
+        Either r or n must be provided, but not both.
+        """
+        if (r is None and n is None) or (r is not None and n is not None):
+            raise ValueError("Exactly one of r or n must be provided.")
+
+        xyz = np.asarray(xyz).ravel()
+        if xyz.shape != (3,):
+            raise ValueError("xyz must be a 1 x 3 array-like object.")
+
+        # Extract all points from all analyses
+        all_points = []
+        analysis_ids = []
+        for study in self.studies:
+            for analysis in study.analyses:
+                for point in analysis.points:
+                    if hasattr(point, "x") and hasattr(point, "y") and hasattr(point, "z"):
+                        all_points.append([point.x, point.y, point.z])
+                        analysis_ids.append(analysis.id)
+
+        if not all_points:  # Return empty list if no coordinates found
+            return []
+
+        all_points = np.array(all_points)
+
+        # Calculate Euclidean distances to all points
+        distances = np.sqrt(np.sum((all_points - xyz) ** 2, axis=1))
+
+        if r is not None:
+            # Find analyses with points within radius r
+            within_radius = distances <= r
+            found_analyses = set(np.array(analysis_ids)[within_radius])
+        else:
+            # Find n closest analyses
+            closest_n_idx = np.argsort(distances)[:n]
+            found_analyses = set(np.array(analysis_ids)[closest_n_idx])
+
+        return list(found_analyses)
 
     def get_analyses_by_mask(self, img):
-        """Extract a list of Analyses with at least one Point in the specified mask."""
-        raise NotImplementedError("Getting analyses by mask is not yet supported.")
+        """Extract a list of Analyses with at least one Point in the specified mask.
 
-    def get_analyses_by_annotations(self):
-        """Extract a list of Analyses with a given label/annotation."""
-        raise NotImplementedError("Getting analyses by annotations is not yet supported.")
+        Parameters
+        ----------
+        img : img_like
+            Mask across which to search for coordinates.
 
-    def get_analyses_by_texts(self):
-        """Extract a list of Analyses with a given text."""
-        raise NotImplementedError("Getting analyses by texts is not yet supported.")
+        Returns
+        -------
+        list[str]
+            A list of Analysis IDs with at least one point in the mask.
+        """
+        # Load mask
+        mask = load_niimg(img)
+
+        # Extract all points from all analyses
+        all_points = []
+        analysis_ids = []
+        for study in self.studies:
+            for analysis in study.analyses:
+                for point in analysis.points:
+                    if hasattr(point, "x") and hasattr(point, "y") and hasattr(point, "z"):
+                        all_points.append([point.x, point.y, point.z])
+                        analysis_ids.append(analysis.id)
+
+        if not all_points:  # Return empty list if no coordinates found
+            return []
+
+        # Convert to voxel coordinates
+        all_points = np.array(all_points)
+        ijk = mm2vox(all_points, mask.affine)
+
+        # Get mask coordinates
+        mask_data = mask.get_fdata()
+        mask_coords = np.vstack(np.where(mask_data)).T
 
-    def get_analyses_by_images(self):
-        """Extract a list of Analyses with a given image."""
-        raise NotImplementedError("Getting analyses by images is not yet supported.")
+        # Check for presence of coordinates in mask
+        in_mask = np.any(np.all(ijk[:, None] == mask_coords[None, :], axis=-1), axis=-1)
 
-    def get_analyses_by_metadata(self):
+        # Get unique analysis IDs where points are in mask
+        found_analyses = set(np.array(analysis_ids)[in_mask])
+
+        return list(found_analyses)
+
+    def get_analyses_by_annotations(self, key, value=None):
+        """Extract a list of Analyses with a given label/annotation."""
+        annotations = {}
+        for study in self.studies:
+            for analysis in study.analyses:
+                a_annot = analysis.annotations
+                if key in a_annot and (value is None or a_annot[key] == value):
+                    annotations[analysis.id] = {key: a_annot[key]}
+        return annotations
+
+    def get_analyses_by_metadata(self, key, value=None):
         """Extract a list of Analyses with a metadata field/value."""
-        raise NotImplementedError("Getting analyses by metadata is not yet supported.")
+        metadata = {}
+        for study in self.studies:
+            for analysis in study.analyses:
+                a_metadata = analysis.metadata
+                if key in a_metadata and (value is None or a_metadata[key] == value):
+                    metadata[analysis.id] = {key: a_metadata[key]}
+        return metadata
 
     def get_points(self, analyses):
         """Collect Points associated with specified Analyses."""
-        raise NotImplementedError("Getting points is not yet supported.")
+        points = {}
+        for study in self.studies:
+            for analysis in study.analyses:
+                if analysis.id in analyses:
+                    points[analysis.id] = analysis.points
+        return points
 
     def get_annotations(self, analyses):
         """Collect Annotations associated with specified Analyses."""
-        raise NotImplementedError("Getting annotations is not yet supported.")
+        annotations = {}
+        for study in self.studies:
+            for analysis in study.analyses:
+                if analysis.id in analyses:
+                    annotations[analysis.id] = analysis.annotations
+
+        return annotations
 
     def get_texts(self, analyses):
         """Collect texts associated with specified Analyses."""
@@ -191,11 +398,32 @@ class Studyset:
 
     def get_images(self, analyses):
         """Collect image files associated with specified Analyses."""
-        raise NotImplementedError("Getting images is not yet supported.")
+        images = {}
+        for study in self.studies:
+            for analysis in study.analyses:
+                if analysis.id in analyses:
+                    images[analysis.id] = analysis.images
+        return images
 
     def get_metadata(self, analyses):
-        """Collect metadata associated with specified Analyses."""
-        raise NotImplementedError("Getting metadata is not yet supported.")
+        """Collect metadata associated with specified Analyses.
+
+        Parameters
+        ----------
+        analyses : list of str
+            List of Analysis IDs to get metadata for.
+
+        Returns
+        -------
+        dict[str, dict]
+            Dictionary mapping Analysis IDs to their combined metadata (including study metadata).
+        """
+        metadata = {}
+        for study in self.studies:
+            for analysis in study.analyses:
+                if analysis.id in analyses:
+                    metadata[analysis.id] = analysis.get_metadata()
+        return metadata
 
 
 class Study:
@@ -226,7 +454,7 @@ class Study:
         self.authors = source["authors"] or ""
         self.publication = source["publication"] or ""
         self.metadata = source.get("metadata", {}) or {}
-        self.analyses = [Analysis(a) for a in source["analyses"]]
+        self.analyses = [Analysis(a, study=self) for a in source["analyses"]]
 
     def __repr__(self):
         """My Simple representation."""
@@ -285,7 +513,7 @@ class Analysis:
     Should the images attribute be a list instead, if the Images contain type information?
     """
 
-    def __init__(self, source):
+    def __init__(self, source, study=None):
         self.id = source["id"]
         self.name = source["name"]
         self.conditions = [
@@ -295,6 +523,7 @@ class Analysis:
         self.points = [Point(p) for p in source["points"]]
         self.metadata = source.get("metadata", {}) or {}
         self.annotations = {}
+        self._study = weakref.proxy(study) if study else None
 
     def __repr__(self):
         """My Simple representation."""
@@ -306,6 +535,22 @@ class Analysis:
             " ".join([self.name, f"images: {len(self.images)}", f"points: {len(self.points)}"])
         )
 
+    def get_metadata(self) -> "dict[str, any]":
+        """Get combined metadata from both analysis and parent study.
+
+        Returns
+        -------
+        dict[str, any]
+            Combined metadata dictionary with analysis metadata taking precedence
+            over study metadata for any overlapping keys.
+        """
+        if self._study is None:
+            return self.metadata.copy()
+
+        combined_metadata = self._study.metadata.copy()
+        combined_metadata.update(self.metadata)
+        return combined_metadata
+
     def to_dict(self):
         """Convert the Analysis to a dictionary."""
         return {
@@ -318,6 +563,7 @@ class Analysis:
             "images": [i.to_dict() for i in self.images],
             "points": [p.to_dict() for p in self.points],
             "weights": [c.to_dict()["weight"] for c in self.conditions],
+            "metadata": self.metadata,
         }
 
 


Test patch:
diff --git a/nimare/tests/test_io.py b/nimare/tests/test_io.py
index d1207d0..833162a 100644
--- a/nimare/tests/test_io.py
+++ b/nimare/tests/test_io.py
@@ -64,18 +64,6 @@ def test_analysis_to_dict_invalid_sample_sizes_type(example_nimads_studyset):
         io.convert_nimads_to_dataset(studyset)
 
 
-def test_analysis_to_dict_invalid_annotations_format(example_nimads_studyset):
-    """Test _analysis_to_dict raises ValueError when annotations are in an invalid format."""
-    studyset = Studyset(example_nimads_studyset)
-    # Here we assume that the annotation is expected to be a dict
-    # Set annotation to an invalid format (e.g., a string)
-    for study in studyset.studies:
-        for analysis in study.analyses:
-            analysis.metadata["annotations"] = "invalid_format"
-    with pytest.raises(TypeError):
-        io.convert_nimads_to_dataset(studyset)
-
-
 def test_convert_sleuth_to_dataset_smoke():
     """Smoke test for Sleuth text file conversion."""
     sleuth_file = os.path.join(get_test_data_path(), "test_sleuth_file.txt")
diff --git a/nimare/tests/test_nimads.py b/nimare/tests/test_nimads.py
index 2261ffa..8fcafff 100644
--- a/nimare/tests/test_nimads.py
+++ b/nimare/tests/test_nimads.py
@@ -1,5 +1,11 @@
 """Test NiMADS functionality."""
 
+import json
+import os
+import tempfile
+
+import pytest
+
 from nimare import nimads
 from nimare.dataset import Dataset
 
@@ -19,3 +25,260 @@ def test_load_nimads(example_nimads_studyset, example_nimads_annotation):
     assert isinstance(filtered_studyset, nimads.Studyset)
     dataset = filtered_studyset.to_dataset()
     assert isinstance(dataset, Dataset)
+
+
+def test_slice_preserves_metadata_and_annotations(
+    example_nimads_studyset, example_nimads_annotation
+):
+    """Test that slicing preserves both metadata and annotations.
+
+    This test verifies that both metadata attached to analyses and annotation
+    notes are correctly preserved when slicing a studyset.
+    """
+    studyset = nimads.Studyset(example_nimads_studyset)
+    studyset.annotations = example_nimads_annotation
+
+    # Get analysis IDs from the first annotation
+    annotation = studyset.annotations[0]
+    analysis_ids = [n.analysis.id for n in annotation.notes if n.note["include"]]
+    selected_ids = analysis_ids[:2]  # Take first two analyses
+
+    # Add metadata to the analyses we'll keep
+    metadata_map = {}
+    for study in studyset.studies:
+        for analysis in study.analyses:
+            if analysis.id in selected_ids:
+                analysis.metadata = {
+                    "sample_size": 30,
+                    "contrast_type": "activation",
+                    "significance_threshold": 0.001,
+                }
+                metadata_map[analysis.id] = analysis.metadata
+
+    # Slice studyset
+    sliced_studyset = studyset.slice(analyses=selected_ids)
+
+    # Verify analyses and their metadata are preserved
+    for study in sliced_studyset.studies:
+        for analysis in study.analyses:
+            assert analysis.id in selected_ids
+            assert analysis.metadata == metadata_map[analysis.id]
+
+    # Verify annotations are preserved for remaining analyses
+    sliced_annotation = sliced_studyset.annotations[0]
+    sliced_analysis_ids = [n.analysis.id for n in sliced_annotation.notes]
+    sliced_annotation_notes = {n.analysis.id: n.note for n in sliced_annotation.notes}
+
+    # Check that notes exist only for remaining analyses
+    assert set(sliced_analysis_ids) == set(selected_ids)
+
+    # Check that annotation contents are preserved
+    for analysis_id in selected_ids:
+        original_note = next(n.note for n in annotation.notes if n.analysis.id == analysis_id)
+        assert sliced_annotation_notes[analysis_id] == original_note
+
+
+def test_studyset_init(example_nimads_studyset):
+    """Test Studyset initialization."""
+    # Test initialization with dict
+    studyset1 = nimads.Studyset(example_nimads_studyset)
+    assert studyset1.id == example_nimads_studyset["id"]
+    assert studyset1.name == example_nimads_studyset["name"]
+    assert len(studyset1.studies) == len(example_nimads_studyset["studies"])
+
+    # Test initialization with JSON file
+    with tempfile.NamedTemporaryFile(mode="w", suffix=".json", delete=False) as tmp:
+        json.dump(example_nimads_studyset, tmp)
+        tmp_path = tmp.name
+
+    try:
+        studyset2 = nimads.Studyset(tmp_path)
+        assert studyset2.id == example_nimads_studyset["id"]
+        assert studyset2.name == example_nimads_studyset["name"]
+        assert len(studyset2.studies) == len(example_nimads_studyset["studies"])
+    finally:
+        os.unlink(tmp_path)
+
+
+def test_studyset_string_methods(example_nimads_studyset):
+    """Test string representation methods."""
+    studyset = nimads.Studyset(example_nimads_studyset)
+
+    # Test __repr__
+    assert repr(studyset) == f"'<Studyset: {studyset.id}>'"
+
+    # Test __str__
+    expected_str = f"Studyset: {studyset.name} :: studies: {len(studyset.studies)}"
+    assert str(studyset) == expected_str
+
+
+def test_studyset_save_load(example_nimads_studyset):
+    """Test saving and loading Studyset."""
+    studyset = nimads.Studyset(example_nimads_studyset)
+
+    with tempfile.NamedTemporaryFile(suffix=".pkl", delete=False) as tmp:
+        tmp_path = tmp.name
+
+    try:
+        # Test save
+        studyset.save(tmp_path)
+        assert os.path.exists(tmp_path)
+
+        # Test load
+        new_studyset = nimads.Studyset({"id": "temp", "name": "", "studies": []})
+        new_studyset.load(tmp_path)
+
+        assert new_studyset.id == studyset.id
+        assert new_studyset.name == studyset.name
+        assert len(new_studyset.studies) == len(studyset.studies)
+    finally:
+        os.unlink(tmp_path)
+
+
+def test_studyset_to_dict(example_nimads_studyset):
+    """Test conversion to dictionary."""
+    studyset = nimads.Studyset(example_nimads_studyset)
+    result = studyset.to_dict()
+
+    assert isinstance(result, dict)
+    assert "id" in result
+    assert "name" in result
+    assert "studies" in result
+    assert len(result["studies"]) == len(studyset.studies)
+
+
+def test_studyset_to_nimads(example_nimads_studyset):
+    """Test saving to NIMADS format."""
+    studyset = nimads.Studyset(example_nimads_studyset)
+
+    with tempfile.NamedTemporaryFile(suffix=".json", delete=False) as tmp:
+        tmp_path = tmp.name
+
+    try:
+        studyset.to_nimads(tmp_path)
+        assert os.path.exists(tmp_path)
+
+        # Verify the saved file can be loaded
+        with open(tmp_path, "r") as f:
+            saved_data = json.load(f)
+
+        assert saved_data["id"] == studyset.id
+        assert saved_data["name"] == studyset.name
+        assert len(saved_data["studies"]) == len(studyset.studies)
+    finally:
+        os.unlink(tmp_path)
+
+
+def test_studyset_copy(example_nimads_studyset):
+    """Test copying of Studyset."""
+    studyset = nimads.Studyset(example_nimads_studyset)
+    copied = studyset.copy()
+
+    assert copied is not studyset
+    assert copied.id == studyset.id
+    assert copied.name == studyset.name
+    assert len(copied.studies) == len(studyset.studies)
+
+
+def test_studyset_merge(example_nimads_studyset):
+    """Test merging of Studysets."""
+    studyset1 = nimads.Studyset(example_nimads_studyset)
+
+    # Create a modified copy for merging
+    modified_data = example_nimads_studyset.copy()
+    modified_data["id"] = "other_id"
+    modified_data["name"] = "Other name"
+    studyset2 = nimads.Studyset(modified_data)
+
+    merged = studyset1.merge(studyset2)
+
+    assert isinstance(merged, nimads.Studyset)
+    assert merged.id == f"{studyset1.id}_{studyset2.id}"
+    assert merged.name == f"Merged: {studyset1.name} + {studyset2.name}"
+
+    # Test invalid merge
+    with pytest.raises(ValueError):
+        studyset1.merge("not a studyset")
+
+
+def test_get_analyses_by_coordinates(example_nimads_studyset):
+    """Test retrieving analyses by coordinates."""
+    studyset = nimads.Studyset(example_nimads_studyset)
+
+    # Test with radius
+    xyz = [0, 0, 0]
+    results_r = studyset.get_analyses_by_coordinates(xyz, r=10)
+    assert isinstance(results_r, list)
+
+    # Test with n nearest
+    results_n = studyset.get_analyses_by_coordinates(xyz, n=5)
+    assert isinstance(results_n, list)
+    assert len(results_n) <= 5
+
+    # Test invalid parameters
+    with pytest.raises(ValueError):
+        studyset.get_analyses_by_coordinates(xyz)  # Neither r nor n
+    with pytest.raises(ValueError):
+        studyset.get_analyses_by_coordinates(xyz, r=10, n=5)  # Both r and n
+    with pytest.raises(ValueError):
+        studyset.get_analyses_by_coordinates([0, 0])  # Invalid coordinates
+
+
+def test_get_analyses_by_mask(example_nimads_studyset, mni_mask):
+    """Test retrieving analyses by mask."""
+    studyset = nimads.Studyset(example_nimads_studyset)
+
+    results = studyset.get_analyses_by_mask(mni_mask)
+    assert isinstance(results, list)
+
+
+def test_get_analyses_by_metadata(example_nimads_studyset):
+    """Test retrieving analyses by metadata."""
+    studyset = nimads.Studyset(example_nimads_studyset)
+
+    # Add some metadata for testing
+    key = "test_key"
+    value = "test_value"
+    for study in studyset.studies:
+        for analysis in study.analyses:
+            analysis.metadata[key] = value
+
+    # Test with key only
+    results1 = studyset.get_analyses_by_metadata(key)
+    assert isinstance(results1, dict)
+
+    # Test with key and value
+    results2 = studyset.get_analyses_by_metadata(key, value)
+    assert isinstance(results2, dict)
+    assert all(list(d.values())[0] == value for d in results2.values())
+
+
+def test_data_retrieval_methods(example_nimads_studyset):
+    """Test methods that retrieve data for specified analyses."""
+    studyset = nimads.Studyset(example_nimads_studyset)
+
+    # Get some analysis IDs to test with
+    analysis_ids = []
+    for study in studyset.studies:
+        for analysis in study.analyses:
+            analysis_ids.append(analysis.id)
+            if len(analysis_ids) >= 2:  # Just test with first two analyses
+                break
+        if len(analysis_ids) >= 2:
+            break
+
+    # Test get_points
+    points = studyset.get_points(analysis_ids)
+    assert isinstance(points, dict)
+
+    # Test get_images
+    images = studyset.get_images(analysis_ids)
+    assert isinstance(images, dict)
+
+    # Test get_metadata
+    metadata = studyset.get_metadata(analysis_ids)
+    assert isinstance(metadata, dict)
+
+    # Test get_annotations
+    annotations = studyset.get_annotations(analysis_ids)
+    assert isinstance(annotations, dict)