Add command to extract mean timeseries activity from ABA regions

src/mesoscopy/export/__init__.py

@@ -28,8 +28,8 @@
 import numpy as np
 
 import mesoscopy.export.nwb as exp_nwb
-import mesoscopy.process as proc
 
+from mesoscopy import io
 from tqdm import tqdm
 
 
@@ -97,7 +97,7 @@ def export_deltaf(path: str, out_path: str) -> str:
     nwb = bool(path.endswith(".nwb"))
 
     click.echo("Loading delta F data...")
-    _, deltaf, _ = proc.load_deltaf(path, nwb)
+    _, deltaf, _ = io.load_deltaf(path, nwb)
 
     # Rescale deltaf to positive integer values between 0 and 255
     click.echo("Rescaling delta F for export...")

src/mesoscopy/io.py

@@ -23,6 +23,7 @@
 from collections import OrderedDict
 
 import h5py
+import numpy as np
 import numpy.typing as npt
 import xmltodict
 import zarr
@@ -142,6 +143,30 @@ def store_interim(
     return zarr.load(interim_path)
 
 
+def load_deltaf(path: str, nwb: bool = False) -> tuple[str, np.ndarray, np.ndarray]:
+    """Load preprocessed dF/F from an HDF5 or NWB file.
+
+    Args:
+        path (str): Path to the preprocessed file.
+        nwb (bool, optional): Whether the file is an NWB file. Defaults to False.
+
+    Returns:
+        tuple[str, np.ndarray, np.ndarray]: Session identifier, dF/F series, and timestamps.
+    """
+    if nwb:
+        nwbfile = read_nwb(path)
+        session_id = nwbfile.identifier
+        deltaf_series = np.array(nwbfile.processing["ophys"]["DeltaFSeries"].data)
+        timestamps = np.array(nwbfile.processing["ophys"]["DeltaFSeries"].timestamps)
+    else:
+        session_id = path.split("/")[-1].replace(".h5", "")
+        with h5py.File(path, "r") as f:
+            deltaf_series = f["/F"][:]
+            timestamps = f["/timestamps"][:]
+
+    return session_id, deltaf_series, timestamps
+
+
 def read_points(path: str) -> dict[str, tuple[float, float]]:
     """Read a landmark points file.
 

src/mesoscopy/process/__init__.py

@@ -22,16 +22,17 @@
 """Processing submodule."""
 
 import os
+from pathlib import Path
+
 import click
-import numpy as np
-import dask.array as da
+import pandas as pd
 from pynwb.image import ImageSeries
 
-import mesoscopy.timer as timer
-
-import mesoscopy.io as io
+import mesoscopy.process.region as pr
 import mesoscopy.process.smooth as psm
 import mesoscopy.process.zscore as pzs
+from mesoscopy import io
+from mesoscopy import timer
 
 
 @click.group("process")
@@ -61,12 +62,12 @@ def smooth_cmd(path: str, out_dir: str, sigma: int = 2) -> None:
     """Generate a smoothed DeltaF/F recording using a Laplace of Gaussian filter."""
     if not os.path.exists(out_dir):
         click.echo(f"Creating output directory {out_dir}...")
-        os.makedirs(out_dir)
+        Path(out_dir).mkdir(parents=True)
 
     click.echo(f"Loading preprocessed recording from {path}...")
     # Determine whether we're working with an NWB file
     nwb = bool(path.endswith(".nwb"))
-    session_id, deltaf_series, timestamps = load_deltaf(path, nwb=nwb)
+    session_id, deltaf_series, timestamps = io.load_deltaf(path, nwb=nwb)
 
     outpath = out_dir + os.sep + session_id + "_smoothed.h5"
 
@@ -97,14 +98,14 @@ def smooth_cmd(path: str, out_dir: str, sigma: int = 2) -> None:
 )
 def zscore_cmd(path: str, out_dir: str) -> None:
     """Pixel-wise z-score ∆F/F signal."""
-    if not os.path.exists(out_dir):
+    if not Path(out_dir).exists():
         click.echo(f"Creating output directory {out_dir}...")
-        os.makedirs(out_dir)
+        Path(out_dir).mkdir(parents=True)
 
     click.echo(f"Loading preprocessed recording from {path}...")
     # Determine whether we're working with an NWB file
     nwb = bool(path.endswith(".nwb"))
-    session_id, deltaf_series, timestamps = load_deltaf(path, nwb=nwb)
+    session_id, deltaf_series, timestamps = io.load_deltaf(path, nwb=nwb)
 
     h5_outpath = out_dir + os.sep + session_id + "_zscored.h5"
     with timer.Timer(message="Z-scoring DeltaF/F"):
@@ -145,25 +146,37 @@ def zscore_cmd(path: str, out_dir: str) -> None:
         io.write_nwb(path, nwbfile, io=nwbio)
 
 
-def load_deltaf(path: str, nwb: bool = False) -> tuple[str, np.ndarray, np.ndarray]:
-    """Load preprocessed deltaf from an HDF5 or NWB file.
+@process_cmd.command("regions")
+@click.argument(
+    "path",
+    type=click.Path(exists=True),
+)
+@click.option(
+    "-o",
+    "--out_dir",
+    type=click.Path(dir_okay=True),
+    default="./",
+    help="Output directory for smoothed recording.",
+)
+def regions_cmd(path: str, out_dir: str) -> None:
+    """Extract ∆F signal averages from ABA-defined regions."""
+    if not Path(out_dir).exists():
+        click.echo(f"Creating output directory {out_dir}...")
+        Path(out_dir).mkdir(parents=True)
+
+    click.echo(f"Loading preprocessed recording from {path}...")
+    # Determine whether we're working with an NWB file
+    nwb = bool(path.endswith(".nwb"))
+    session_id, deltaf_series, timestamps = io.load_deltaf(path, nwb=nwb)
 
-    Args:
-        path (str): Path to the preprocessed file.
-        nwb (bool, optional): Whether the file is an NWB file. Defaults to False.
+    outpath = out_dir + os.sep + session_id + "_regions.csv"
 
-    Returns:
-        tuple[str, np.ndarray, np.ndarray]: Session identifier, dF/F series, and timestamps.
-    """
-    if nwb:
-        nwbfile = io.read_nwb(path)
-        session_id = nwbfile.identifier
-        deltaf_series = nwbfile.processing["ophys"]["DeltaFSeries"].data
-        timestamps = nwbfile.processing["ophys"]["DeltaFSeries"].timestamps
-    else:
-        session_id = path.split("/")[-1].replace(".h5", "")
-        f_preproc = io.read_h5(path)
-        deltaf_series = f_preproc["/F"]
-        timestamps = f_preproc["/timestamps"]
-
-    return session_id, np.array(deltaf_series), np.array(timestamps)
+    with timer.Timer(message="Extracting region activity"):
+        region_activity = pd.DataFrame(pr.extract_all_regions(deltaf_series, as_dataframe=True))
+        region_activity["time_idx"] = [
+            str(timestamp, encoding="utf-8") for timestamp in timestamps[region_activity["time_idx"]]
+        ]
+        region_activity.rename(columns={"time_idx": "timestamp"}, inplace=True)
+        region_activity.to_csv(outpath, index=False)
+
+    click.echo(f"Saved region activity at {outpath}")

src/mesoscopy/process/region.py

@@ -0,0 +1,136 @@
+#  Copyright (c) 2026 Constantinos Eleftheriou <Constantinos.Eleftheriou@ed.ac.uk>.
+#
+#   Permission is hereby granted, free of charge, to any person obtaining a copy of this
+#   software and associated documentation files (the "Software"), to deal in the
+#   Software without restriction, including without limitation the rights to use, copy,
+#   modify, merge, publish, distribute, sublicense, and/or sell copies of the Software,
+#   and to permit persons to whom the Software is furnished to do so, subject to the
+#  following conditions:
+#
+#  The above copyright notice and this permission notice shall be included in all copies
+#  or substantial portions of the Software
+#
+#  THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
+#  EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
+#  MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
+#  NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS
+#  BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER
+#  IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR
+#  IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
+#  SOFTWARE.
+
+from typing import Literal
+
+import numpy as np
+import numpy.typing as npt
+import pandas as pd
+
+from mesoscopy import resources
+
+DEFAULT_EXCLUDE = [
+    "FRP1",
+    "VISpl1",
+    "VISpor1",
+    "VISli1",
+    "TEa1",
+    "AUDd1",
+    "AUDp1",
+    "AUDpo1",
+    "AUDv1",
+    "ORBm1",
+]
+
+
+def extract_region_activity(
+    deltaf_series: npt.NDArray, region_acronym: str, hemisphere: Literal["left", "right", "both"]
+) -> npt.NDArray:
+    """Extracts the mean ∆F/F signal of a specific cortical region in the Allen Brain Atlas from a registered recording.
+
+    Args:
+        deltaf_series (npt.NDArray): A 3D array of shape (time, height, width) representing the DeltaF/F signal.
+            Should be registered to the Allen Brain Atlas.
+        region_acronym (str): The acronym of the cortical region to extract (e.g., 'VISp', 'MOp etc.).
+        hemisphere (Literal["left", "right", "both"]): The hemisphere to extract the activity from.
+            Options are 'left', 'right', or 'both'.
+
+    Returns:
+        npt.NDArray: A 1D array of shape (time,) representing the mean ∆F/F signal of the specified cortical
+            region over time.
+
+    Raises:
+        ValueError: If the specified region acronym is not recognized or if the hemisphere option is invalid.
+    """
+    annotations = resources.get_atlas_annotations()
+    region_id = annotations.loc[annotations["acronym"] == region_acronym, "id"].values[0]
+
+    if not region_id:
+        msg = f"Region acronym {region_acronym} not recognised."
+        raise ValueError(msg)
+
+    left_aba, right_aba = resources.get_atlas()
+
+    if hemisphere.lower() == "left":
+        region_mask = np.broadcast_to(left_aba == region_id, deltaf_series.shape)
+    elif hemisphere.lower() == "right":
+        region_mask = np.broadcast_to(right_aba == region_id, deltaf_series.shape)
+    elif hemisphere.lower() == "both":
+        aba = left_aba + right_aba
+        region_mask = np.broadcast_to(aba == region_id, deltaf_series.shape)
+    else:
+        msg = f"Could not recognise hemisphere option {hemisphere}, select left, right or both."
+        raise ValueError(msg)
+
+    return np.ma.array(deltaf_series, mask=~region_mask).mean(axis=(1, 2))
+
+
+def extract_all_regions(
+    deltaf_series: npt.NDArray,
+    exclude: list | None = None,
+    ignore_default_exclude: bool = False,
+    as_dataframe: bool = False,
+) -> dict | pd.DataFrame:
+    """Extracts the mean ∆F/F signal of all cortical regions in the Allen Brain Atlas from a registered recording.
+
+    Args:
+        deltaf_series (npt.NDArray): A 3D array of shape (time, height, width) representing the DeltaF/F signal.
+        exclude (list, optional): A list of region acronyms to exclude from the extraction. Defaults to None.
+        ignore_default_exclude (bool, optional): If True, ignores the default excluded regions. Defaults to False.
+        as_dataframe (bool, optional): If True, returns the result as a pandas DataFrame. Defaults to False.
+
+    Returns:
+        dict | pd.DataFrame: A dictionary with region acronyms as keys and their mean activity as values,
+            or a DataFrame with columns 'region', 'time_idx', and 'F'.
+    """
+    annotations = resources.get_atlas_annotations()
+    left_aba, right_aba = resources.get_atlas()
+
+    if exclude is None:
+        exclude = []
+
+    excluded_regions = list(DEFAULT_EXCLUDE) if not ignore_default_exclude else []
+    excluded_regions.extend(exclude)
+    regions = [region for region in annotations.acronym.unique() if region not in excluded_regions]
+
+    region_ids = {region: annotations.loc[annotations["acronym"] == region, "id"].values[0] for region in regions}
+
+    hw = left_aba.shape[0] * left_aba.shape[1]
+    left_masks = np.array([(left_aba == region_ids[r]).ravel() for r in regions])   # (n_regions, H*W)
+    right_masks = np.array([(right_aba == region_ids[r]).ravel() for r in regions])
+
+    left_counts = left_masks.sum(axis=1).astype(float)   # (n_regions,)
+    right_counts = right_masks.sum(axis=1).astype(float)
+
+    data_flat = deltaf_series.reshape(deltaf_series.shape[0], hw)  # (T, H*W)
+    left_activities = (data_flat @ left_masks.T) / left_counts     # (T, n_regions)
+    right_activities = (data_flat @ right_masks.T) / right_counts
+
+    region_activity = {}
+    for i, region in enumerate(regions):
+        region_activity[f"L_{region}"] = left_activities[:, i]
+        region_activity[f"R_{region}"] = right_activities[:, i]
+
+    if as_dataframe:
+        df = pd.DataFrame(region_activity)
+        return df.unstack().reset_index().rename(columns={"level_0": "region", "level_1": "time_idx", 0: "F"})
+
+    return region_activity

src/mesoscopy/register/__init__.py

@@ -176,7 +176,7 @@ def landmarks_cmd(
     # Determine whether we're working with an NWB file
     nwb = True if path.endswith(".nwb") else False
 
-    session_id, deltaf_series, timestamps = load_deltaf(path, nwb)
+    session_id, deltaf_series, timestamps = io.load_deltaf(path, nwb)
 
     click.echo("Loading landmarks...")
     template_landmarks = res.get_default_landmarks()
@@ -224,30 +224,6 @@ def landmarks_cmd(
     return outpath
 
 
-def load_deltaf(path: str, nwb: bool = False) -> tuple[str, np.ndarray, np.ndarray]:
-    """Load preprocessed deltaf from an HDF5 or NWB file.
-
-    Args:
-        path (str): Path to the preprocessed file.
-        nwb (bool, optional): Whether the file is an NWB file. Defaults to False.
-
-    Returns:
-        tuple[str, np.ndarray, np.ndarray]: Session identifier, dF/F series, and timestamps.
-    """
-    if nwb:
-        nwbfile = io.read_nwb(path)
-        session_id = nwbfile.identifier
-        deltaf_series = np.array(nwbfile.processing["ophys"]["DeltaFSeries"].data)
-        timestamps = np.array(nwbfile.processing["ophys"]["DeltaFSeries"].timestamps)
-    else:
-        session_id = path.split("/")[-1].replace(".h5", "")
-        with h5py.File(path, "r") as f_preproc:
-            deltaf_series = f_preproc["/F"][:]
-            timestamps = f_preproc["/timestamps"][:]
-
-    return session_id, deltaf_series, timestamps
-
-
 def load_maxips(path: str) -> tuple[np.ndarray, np.ndarray]:
     """Load maximum intensity projections from a preprocessed HDF5 file.
 

src/mesoscopy/resources/__init__.py

@@ -20,6 +20,10 @@
 #  SOFTWARE.
 from importlib import resources
 
+import imageio
+import numpy as np
+import pandas as pd
+
 import mesoscopy.resources
 from mesoscopy import io
 
@@ -30,10 +34,18 @@ def get_default_landmarks() -> dict:
     Returns:
         dict: Default landmark coordinates.
     """
-    return io.read_points(
-        str(
-            resources.files(mesoscopy.resources).joinpath(
-                "ccf_template_landmarks_140x142.csv"
-            )
-        )
+    return io.read_points(str(resources.files(mesoscopy.resources).joinpath("ccf_template_landmarks_140x142.csv")))
+
+
+def get_atlas():
+    left_hemisphere_aba = imageio.imread(
+        str(resources.files(mesoscopy.resources).joinpath("ccf_template_top_140x142.tiff"))
+    )
+    right_hemisphere_aba = np.flip(left_hemisphere_aba, axis=1)
+    return left_hemisphere_aba, right_hemisphere_aba
+
+
+def get_atlas_annotations():
+    return pd.read_csv(
+        str(resources.files(mesoscopy.resources).joinpath("ccf_annotations.csv")), delimiter=", ", engine="python"
     )