From 996e66c238f3be0a11cf2e8e85b6500715332903 Mon Sep 17 00:00:00 2001 From: Talley Lambert Date: Thu, 6 Jun 2024 20:48:11 -0400 Subject: [PATCH 1/7] feat: pull neurons from allen brain --- src/microsim/allen/__init__.py | 3 + src/microsim/allen/fetch.py | 230 ++++++++++++++++++ .../schema/sample/matslines/_bresenham.py | 34 ++- src/microsim/util.py | 26 +- x.py | 9 + 5 files changed, 297 insertions(+), 5 deletions(-) create mode 100644 src/microsim/allen/__init__.py create mode 100644 src/microsim/allen/fetch.py create mode 100644 x.py diff --git a/src/microsim/allen/__init__.py b/src/microsim/allen/__init__.py new file mode 100644 index 0000000..d1b14d1 --- /dev/null +++ b/src/microsim/allen/__init__.py @@ -0,0 +1,3 @@ +from .fetch import ApiCellTypesSpecimenDetail, Specimen, get_reconstructions + +__all__ = ["get_reconstructions", "Specimen", "ApiCellTypesSpecimenDetail"] diff --git a/src/microsim/allen/fetch.py b/src/microsim/allen/fetch.py new file mode 100644 index 0000000..87fa370 --- /dev/null +++ b/src/microsim/allen/fetch.py @@ -0,0 +1,230 @@ +from collections import defaultdict +from collections.abc import Sequence +from enum import IntEnum +from functools import cache, cached_property +from typing import TYPE_CHECKING, Literal, NamedTuple, cast + +import numpy as np +import requests +from pydantic import BaseModel, Field + +if TYPE_CHECKING: + from collections.abc import Iterable + +ALLEN_ROOT = "http://api.brain-map.org" +ALLEN_V2_API = f"{ALLEN_ROOT}/api/v2/data" +ALLEN_V2_QUERY = ALLEN_V2_API + "/query.json" +SWC_FILE_TYPE = "3DNeuronReconstruction" + + +class WellKnownFileType(BaseModel): + id: int + name: str + + +class WellKnownFile(BaseModel): + attachable_id: int | None + attachable_type: str | None + download_link: str | None + id: int | None + path: str | None + well_known_file_type_id: int | None + well_known_file_type: WellKnownFileType | None + + +class SWCType(IntEnum): + UNDEFINED = 0 + SOMA = 1 + AXON = 2 + BASAL_DENDRITE = 3 + APICAL_DENDRITE = 4 + CUSTOM = 5 + UNSPECIFIED_NEURITE = 6 + GLIA_PROCESSES = 7 + + +class Compartment(NamedTuple): + id: int + t: int # type + x: float + y: float + z: float + r: float # radius + c: int # parent id or "connectivity" + + +class SWC: + def __init__(self, compartments: Sequence[Compartment] = ()): + self.compartments = compartments + + self.children: defaultdict[int, list[Compartment]] = defaultdict(list) + self.node_types: defaultdict[int, list[Compartment]] = defaultdict(list) + self._map: dict[int, Compartment] = {} + for comp in compartments: + self._map[comp.id] = comp + self.children[comp.c].append(comp) + self.node_types[comp.t].append(comp) + + def iter_pairs(self, starting_type: int): + seen = set() + for comp_id, children in self.children.items(): + if comp_id == -1: + continue + comp = self._map[comp_id] + if comp.t == starting_type: + for child in children: + if (comp, child) not in seen: + yield comp, child + seen.add((comp, child)) + + @classmethod + def from_string(cls, content: str | bytes) -> "SWC": + if isinstance(content, bytes): + content = content.decode() + + compartments = [] + for num, line in enumerate(content.splitlines()): + if line.startswith("#"): + continue + try: + a, b, c, d, e, f, g = line.split() + comp = Compartment( + int(a), int(b), float(c), float(d), float(e), float(f), int(g) + ) + except ValueError as e: + raise ValueError(f"Invalid SWC line {num}: {line}") from e + compartments.append(comp) + return cls(compartments) + + @classmethod + def from_url(cls, url: str) -> "SWC": + response = requests.get(url) + response.raise_for_status() + return cls.from_string(response.text) + + @cached_property + def coords(self) -> np.ndarray: + """Return the coordinates of the compartments as (M, 3) array.""" + return np.array([(c.z, c.y, c.x) for c in self.compartments]).astype("float32") + + def empty_grid(self, resolution: float = 1.0) -> np.ndarray: + extent = np.ptp(self.coords, axis=0) + size = np.ceil(extent / resolution).astype(int) + return np.zeros(size, dtype=np.uint8) + + def build_mask( + self, + resolution: float = 1.0, + dend_scale: float = 2, + global_scale_factor: float = 10, + ): + from microsim.schema.sample.matslines._bresenham import bres_draw_segment_3d + + grid = self.empty_grid(resolution) + for par, child in self.iter_pairs(3): + r = int(max(1, 0.5 * global_scale_factor * dend_scale * (par.r + child.r))) + bres_draw_segment_3d( + int(par.x), + int(par.y), + int(par.z), + int(child.x), + int(child.y), + int(child.z), + grid, + 5000, + width=r, + ) + return grid + + +class NeuronReconstruction(BaseModel): + id: int + specimen_id: int + well_known_files: list[WellKnownFile] = Field(default_factory=list) + + @property + def swc_path(self) -> SWC: + """The SWC file for this reconstruction.""" + for f in self.well_known_files: + if f.well_known_file_type.name == SWC_FILE_TYPE: + return ALLEN_ROOT + f.download_link + raise ValueError("No SWC file found for this reconstruction.") + + def load_swc(self) -> SWC: + """Load the SWC file for this reconstruction.""" + return SWC.from_url(self.swc_path) + + +class Specimen(BaseModel): + id: int + is_cell_specimen: bool + specimen_id_path: str + structure_id: int + neuron_reconstructions: list[NeuronReconstruction] = Field(default_factory=list) + + @classmethod + @cache + def fetch(cls, id: int) -> "Specimen": + """Fetch this specimen from the Allen brain map API.""" + q = [ + "model::Specimen", + f"rma::criteria[id$eq{id}],neuron_reconstructions(well_known_files)", + "rma::include,neuron_reconstructions(well_known_files(" + f"well_known_file_type[name$eq'{SWC_FILE_TYPE}']))", + "rma::options[num_rows$eq'all']", + ] + response = requests.get(ALLEN_V2_QUERY, params={"q": ",".join(q)}) + response.raise_for_status() + qr = QueryResponse.model_validate_json(response.content) + if not qr.success: + raise ValueError(qr.msg) + return cast("Specimen", qr.msg[0]) + + +class ApiCellTypesSpecimenDetail(BaseModel): + specimen__id: int + structure__name: str | None + structure__acronym: str | None + donor__species: Literal["Homo Sapiens", "Mus musculus"] + nr__reconstruction_type: str | None # probably just 'full' or 'dendrite-only' + nr__average_contraction: float | None + nr__average_parent_daughter_ratio: float | None + nr__max_euclidean_distance: float | None + nr__number_bifurcations: int | None + nr__number_stems: int | None + + def specimen(self) -> Specimen: + return Specimen.fetch(self.specimen__id) + + @classmethod + @cache + def all_reconstructions(cls) -> tuple["ApiCellTypesSpecimenDetail", ...]: + """Fetch details for all Specimens with reconstruction info.""" + q = ( + "model::ApiCellTypesSpecimenDetail", + "rma::criteria[nr__reconstruction_type$ne'null']", + "rma::options[num_rows$eq'all']", + ) + response = requests.get(ALLEN_V2_QUERY, params={"q": ",".join(q)}) + response.raise_for_status() + qr = QueryResponse.model_validate_json(response.content) + if not qr.success: + raise ValueError(qr.msg) + return tuple(qr.msg) # type: ignore[arg-type] + + +class QueryResponse(BaseModel): + success: bool + msg: list[Specimen] | list[ApiCellTypesSpecimenDetail] | str + + +def get_reconstructions( + species: Literal["Homo Sapiens", "Mus musculus"] | None, + reconstruction_type: Literal["full", "dendrite-only"] | None = None, +) -> tuple[ApiCellTypesSpecimenDetail, ...]: + recons: Iterable = ApiCellTypesSpecimenDetail.all_reconstructions() + if species is not None: + recons = (x for x in recons if x.donor__species == species) + if reconstruction_type is not None: + recons = (x for x in recons if x.nr__reconstruction_type == reconstruction_type) + return tuple(recons) diff --git a/src/microsim/schema/sample/matslines/_bresenham.py b/src/microsim/schema/sample/matslines/_bresenham.py index e04800b..fad1715 100644 --- a/src/microsim/schema/sample/matslines/_bresenham.py +++ b/src/microsim/schema/sample/matslines/_bresenham.py @@ -42,7 +42,15 @@ def bres_draw_segment_2d( def bres_draw_segment_3d( - x0: int, y0: int, z0: int, x1: int, y1: int, z1: int, grid: np.ndarray, max_r: float + x0: int, + y0: int, + z0: int, + x1: int, + y1: int, + z1: int, + grid: np.ndarray, + max_r: float, + width: float = 1.0, ) -> None: """Bresenham's algorithm. @@ -66,9 +74,14 @@ def bres_draw_segment_3d( xr = grid.shape[2] / 2 while True: - r = ((x0 - xr) / xr) ** 2 + ((y0 - yr) / yr) ** 2 + ((z0 - zr) / zr) ** 2 - if sqrt(r) <= max_r: - grid[z0, y0, x0] += 1 + if width != 1: + # Draw a sphere around the current point with the given width + draw_sphere(grid, x0, y0, z0, width) + else: + r = ((x0 - xr) / xr) ** 2 + ((y0 - yr) / yr) ** 2 + ((z0 - zr) / zr) ** 2 + if sqrt(r) <= max_r: + grid[z0, y0, x0] += 1 + if i == 0: break @@ -87,6 +100,18 @@ def bres_draw_segment_3d( i -= 1 +def draw_sphere(grid: np.ndarray, x0: int, y0: int, z0: int, radius: float) -> None: + """Draw a sphere of a given radius around a point in a 3D grid.""" + z_range = range(int(max(0, z0 - radius)), int(min(grid.shape[0], z0 + radius + 1))) + y_range = range(int(max(0, y0 - radius)), int(min(grid.shape[1], y0 + radius + 1))) + x_range = range(int(max(0, x0 - radius)), int(min(grid.shape[2], x0 + radius + 1))) + for z in z_range: + for y in y_range: + for x in x_range: + if (x - x0) ** 2 + (y - y0) ** 2 + (z - z0) ** 2 <= radius**2: + grid[z, y, x] = True + + try: from numba import jit except Exception: @@ -94,3 +119,4 @@ def bres_draw_segment_3d( else: bres_draw_segment_2d = jit(nopython=True)(bres_draw_segment_2d) bres_draw_segment_3d = jit(nopython=True)(bres_draw_segment_3d) + draw_sphere = jit(nopython=True)(draw_sphere) diff --git a/src/microsim/util.py b/src/microsim/util.py index 08766a9..935b21e 100644 --- a/src/microsim/util.py +++ b/src/microsim/util.py @@ -3,7 +3,7 @@ import itertools import shutil import warnings -from typing import TYPE_CHECKING, Protocol, TypeVar, cast +from typing import TYPE_CHECKING, Any, Protocol, TypeVar, cast import numpy as np import numpy.typing as npt @@ -288,3 +288,27 @@ def downsample( for d in range(array.ndim): reshaped = method(reshaped, -1 * (d + 1), dtype) return reshaped + + +def view_nd( + ary: Any, figsize: tuple[int, int] = (1280, 1000), **view_kwargs: Any +) -> None: + try: + from pymmcore_widgets._stack_viewer_v2 import StackViewer + except ImportError as e: + raise ImportError( + "This feature uses a not-yet published widget from pymmcore-widgets." + "It will eventually be made available outside of pymmcore-widgets..." + ) from e + from qtpy.QtWidgets import QApplication + + app = QApplication.instance() + if not (hadapp := app is not None): + app = QApplication([]) + + s = StackViewer(ary, **view_kwargs) + s.resize(*figsize) + s.show() + + if not hadapp: + app.exec() diff --git a/x.py b/x.py new file mode 100644 index 0000000..78a84c3 --- /dev/null +++ b/x.py @@ -0,0 +1,9 @@ +from microsim.allen import Specimen +from microsim.util import view_nd + +spec = Specimen.fetch(555241040) +swc = spec.neuron_reconstructions[0].load_swc() +print(spec.id) + +mask = swc.build_mask(global_scale_factor=2) +view_nd(mask) From 03849fe02455f81a556355b67bffc95c4c031052 Mon Sep 17 00:00:00 2001 From: Talley Lambert Date: Fri, 7 Jun 2024 09:36:58 -0400 Subject: [PATCH 2/7] update binary mask --- .../matslines/_bresenham.py => _draw.py} | 20 ++- src/microsim/allen/fetch.py | 170 ++++++++++++------ .../schema/sample/matslines/_matslines.py | 6 +- src/microsim/util.py | 7 +- x.py | 14 +- 5 files changed, 135 insertions(+), 82 deletions(-) rename src/microsim/{schema/sample/matslines/_bresenham.py => _draw.py} (87%) diff --git a/src/microsim/schema/sample/matslines/_bresenham.py b/src/microsim/_draw.py similarity index 87% rename from src/microsim/schema/sample/matslines/_bresenham.py rename to src/microsim/_draw.py index fad1715..2c17aaf 100644 --- a/src/microsim/schema/sample/matslines/_bresenham.py +++ b/src/microsim/_draw.py @@ -9,7 +9,7 @@ import numpy as np -def bres_draw_segment_2d( +def draw_line_2d( y0: int, x0: int, y1: int, x1: int, grid: np.ndarray, max_r: float ) -> None: """Bresenham's algorithm. @@ -41,7 +41,7 @@ def bres_draw_segment_2d( y0 += sy -def bres_draw_segment_3d( +def draw_line_3d( x0: int, y0: int, z0: int, @@ -73,6 +73,9 @@ def bres_draw_segment_3d( yr = grid.shape[1] / 2 xr = grid.shape[2] / 2 + if max_r < 0: + max_r = sqrt(zr**2 + yr**2 + xr**2) + while True: if width != 1: # Draw a sphere around the current point with the given width @@ -108,15 +111,16 @@ def draw_sphere(grid: np.ndarray, x0: int, y0: int, z0: int, radius: float) -> N for z in z_range: for y in y_range: for x in x_range: - if (x - x0) ** 2 + (y - y0) ** 2 + (z - z0) ** 2 <= radius**2: - grid[z, y, x] = True + distance = (x - x0) ** 2 + (y - y0) ** 2 + (z - z0) ** 2 + if distance <= radius**2: + grid[z, y, x] += 1 try: - from numba import jit + from numba import njit except Exception: pass else: - bres_draw_segment_2d = jit(nopython=True)(bres_draw_segment_2d) - bres_draw_segment_3d = jit(nopython=True)(bres_draw_segment_3d) - draw_sphere = jit(nopython=True)(draw_sphere) + draw_line_2d = njit(draw_line_2d) + draw_line_3d = njit(draw_line_3d) + draw_sphere = njit(draw_sphere) diff --git a/src/microsim/allen/fetch.py b/src/microsim/allen/fetch.py index 87fa370..c171f9f 100644 --- a/src/microsim/allen/fetch.py +++ b/src/microsim/allen/fetch.py @@ -1,7 +1,8 @@ from collections import defaultdict -from collections.abc import Sequence +from collections.abc import Iterator, Sequence from enum import IntEnum from functools import cache, cached_property +from pathlib import Path from typing import TYPE_CHECKING, Literal, NamedTuple, cast import numpy as np @@ -52,30 +53,26 @@ class Compartment(NamedTuple): r: float # radius c: int # parent id or "connectivity" + def coord(self) -> np.ndarray: + return np.array([self.z, self.y, self.x]) -class SWC: - def __init__(self, compartments: Sequence[Compartment] = ()): - self.compartments = compartments + def shifted_coord( + self, origin: np.ndarray, resolution: float = 1 + ) -> tuple[int, int, int]: + shifted_coord = (self.coord() - origin) / resolution + return tuple(c for c in shifted_coord.astype(int)) - self.children: defaultdict[int, list[Compartment]] = defaultdict(list) - self.node_types: defaultdict[int, list[Compartment]] = defaultdict(list) - self._map: dict[int, Compartment] = {} - for comp in compartments: - self._map[comp.id] = comp - self.children[comp.c].append(comp) - self.node_types[comp.t].append(comp) - def iter_pairs(self, starting_type: int): - seen = set() - for comp_id, children in self.children.items(): - if comp_id == -1: - continue - comp = self._map[comp_id] - if comp.t == starting_type: - for child in children: - if (comp, child) not in seen: - yield comp, child - seen.add((comp, child)) +class SWC: + @classmethod + def from_path(cls, path: str | Path) -> "SWC": + if str(path).startswith(("http://", "https://")): + response = requests.get(str(path)) + response.raise_for_status() + content = response.text + else: + content = Path(path).expanduser().read_text() + return cls.from_string(content) @classmethod def from_string(cls, content: str | bytes) -> "SWC": @@ -96,70 +93,117 @@ def from_string(cls, content: str | bytes) -> "SWC": compartments.append(comp) return cls(compartments) - @classmethod - def from_url(cls, url: str) -> "SWC": - response = requests.get(url) - response.raise_for_status() - return cls.from_string(response.text) + def __init__(self, compartments: Sequence[Compartment] = ()): + self.compartments = compartments + + self._id_map: dict[int, Compartment] = {} + self._children_of: defaultdict[int, list[Compartment]] = defaultdict(list) + self._node_types: defaultdict[int, list[Compartment]] = defaultdict(list) + for comp in compartments: + self._id_map[comp.id] = comp + self._children_of[comp.c].append(comp) + self._node_types[comp.t].append(comp) + + def iter_pairs(self, *types: int) -> Iterator[tuple[Compartment, Compartment]]: + seen = set() + for comp_id, children in self._children_of.items(): + if comp_id == -1: + continue + comp = self._id_map[comp_id] + if not types or comp.t in types: + for child in children: + if (comp, child) not in seen: + yield comp, child + seen.add((comp, child)) @cached_property def coords(self) -> np.ndarray: """Return the coordinates of the compartments as (M, 3) array.""" return np.array([(c.z, c.y, c.x) for c in self.compartments]).astype("float32") - def empty_grid(self, resolution: float = 1.0) -> np.ndarray: + def root(self) -> Compartment: + """Return the root compartment of the SWC.""" + return self._id_map[1] + + def origin(self) -> tuple[float, float, float]: + """Return the (Z,Y,X) coordinate of the root node.""" + root = self.root() + return root.z, root.y, root.x + + def _empty_grid(self, voxel_size: float = 1.0) -> np.ndarray: + """Create an empty 3D grid for the binary mask.""" extent = np.ptp(self.coords, axis=0) - size = np.ceil(extent / resolution).astype(int) - return np.zeros(size, dtype=np.uint8) - - def build_mask( - self, - resolution: float = 1.0, - dend_scale: float = 2, - global_scale_factor: float = 10, - ): - from microsim.schema.sample.matslines._bresenham import bres_draw_segment_3d - - grid = self.empty_grid(resolution) - for par, child in self.iter_pairs(3): - r = int(max(1, 0.5 * global_scale_factor * dend_scale * (par.r + child.r))) - bres_draw_segment_3d( - int(par.x), - int(par.y), - int(par.z), - int(child.x), - int(child.y), - int(child.z), - grid, - 5000, - width=r, - ) - return grid + size = np.ceil(extent / voxel_size).astype(int) + return np.zeros(size, dtype=bool) + + def binary_mask(self, voxel_size: float = 1, scale_factor: float = 3) -> np.ndarray: + from microsim._draw import draw_line_3d, draw_sphere + + grid = self._empty_grid(voxel_size) + origin = np.min(self.coords, axis=0) + + dend_scale: float = 1 + max_r = float(np.sum(grid.shape)) + for par, child in self.iter_pairs( + SWCType.BASAL_DENDRITE, SWCType.APICAL_DENDRITE + ): + r = int(max(1, 0.5 * scale_factor * dend_scale * (par.r + child.r))) + pz, py, px = par.shifted_coord(origin, voxel_size) + cz, cy, cx = child.shifted_coord(origin, voxel_size) + draw_line_3d(px, py, pz, cx, cy, cz, grid, max_r=max_r, width=r) + + soma_scale: float = 1.2 + for comp in self._node_types[SWCType.SOMA]: + z, y, x = comp.shifted_coord(origin, voxel_size) + r = int(0.5 * scale_factor * soma_scale * comp.r) + draw_sphere(grid, x, y, z, r) + + return grid.astype(np.uint8) class NeuronReconstruction(BaseModel): id: int specimen_id: int + number_nodes: int + number_branches: int + neuron_reconstruction_type: str + overall_height: float + overall_width: float + overall_depth: float + scale_factor_x: float + scale_factor_y: float + scale_factor_z: float well_known_files: list[WellKnownFile] = Field(default_factory=list) @property - def swc_path(self) -> SWC: + def swc_path(self) -> str: """The SWC file for this reconstruction.""" for f in self.well_known_files: - if f.well_known_file_type.name == SWC_FILE_TYPE: + if ( + getattr(f.well_known_file_type, "name", None) == SWC_FILE_TYPE + and f.download_link + ): return ALLEN_ROOT + f.download_link raise ValueError("No SWC file found for this reconstruction.") def load_swc(self) -> SWC: """Load the SWC file for this reconstruction.""" - return SWC.from_url(self.swc_path) + return SWC.from_path(self.swc_path) + + +class Structure(BaseModel): + id: int + name: str + acronym: str + structure_id_path: str class Specimen(BaseModel): id: int + name: str is_cell_specimen: bool specimen_id_path: str - structure_id: int + structure: Structure neuron_reconstructions: list[NeuronReconstruction] = Field(default_factory=list) @classmethod @@ -169,7 +213,7 @@ def fetch(cls, id: int) -> "Specimen": q = [ "model::Specimen", f"rma::criteria[id$eq{id}],neuron_reconstructions(well_known_files)", - "rma::include,neuron_reconstructions(well_known_files(" + "rma::include,structure,neuron_reconstructions(well_known_files(" f"well_known_file_type[name$eq'{SWC_FILE_TYPE}']))", "rma::options[num_rows$eq'all']", ] @@ -180,6 +224,14 @@ def fetch(cls, id: int) -> "Specimen": raise ValueError(qr.msg) return cast("Specimen", qr.msg[0]) + def binary_mask(self, voxel_size: float = 1, scale_factor: float = 1) -> np.ndarray: + """Return 3D binary mask for this specimen's neuron reconstructions.""" + for recon in self.neuron_reconstructions: + return recon.load_swc().binary_mask( + voxel_size=voxel_size, scale_factor=scale_factor + ) + raise ValueError("No neuron reconstructions found for this specimen.") + class ApiCellTypesSpecimenDetail(BaseModel): specimen__id: int diff --git a/src/microsim/schema/sample/matslines/_matslines.py b/src/microsim/schema/sample/matslines/_matslines.py index 951df9f..0af5f81 100644 --- a/src/microsim/schema/sample/matslines/_matslines.py +++ b/src/microsim/schema/sample/matslines/_matslines.py @@ -73,15 +73,15 @@ def render(self, space: xrDataArray, xp: NumpyAPI | None = None) -> xrDataArray: def drawlines_bresenham( segments: npt.NDArray, grid: npt.NDArray, max_r: float = 2.0 ) -> None: - from ._bresenham import bres_draw_segment_2d, bres_draw_segment_3d + from microsim._draw import draw_line_2d, draw_line_3d if grid.ndim == 2: for segment in segments: y0, x0, y1, x1 = (int(x) for x in segment) - bres_draw_segment_2d(x0, y0, x1, y1, grid, max_r) + draw_line_2d(x0, y0, x1, y1, grid, max_r) elif grid.ndim == 3: for segment in segments: z0, y0, x0, z1, y1, x1 = (int(x) for x in segment) - bres_draw_segment_3d(x0, y0, z0, x1, y1, z1, grid, max_r) + draw_line_3d(x0, y0, z0, x1, y1, z1, grid, max_r) else: raise ValueError(f"grid must be either 2 or 3 dimensional. Got {grid.ndim}") diff --git a/src/microsim/util.py b/src/microsim/util.py index 935b21e..ef970b1 100644 --- a/src/microsim/util.py +++ b/src/microsim/util.py @@ -302,13 +302,10 @@ def view_nd( ) from e from qtpy.QtWidgets import QApplication - app = QApplication.instance() - if not (hadapp := app is not None): - app = QApplication([]) + app = QApplication.instance() or QApplication([]) s = StackViewer(ary, **view_kwargs) s.resize(*figsize) s.show() - if not hadapp: - app.exec() + app.exec() diff --git a/x.py b/x.py index 78a84c3..6cb620e 100644 --- a/x.py +++ b/x.py @@ -1,9 +1,9 @@ -from microsim.allen import Specimen -from microsim.util import view_nd +from rich import print -spec = Specimen.fetch(555241040) -swc = spec.neuron_reconstructions[0].load_swc() -print(spec.id) +from microsim.allen import get_reconstructions +from microsim.util import view_nd -mask = swc.build_mask(global_scale_factor=2) -view_nd(mask) +for recon in get_reconstructions("Mus musculus", "full"): + spec = recon.specimen() + print("viewing", recon, spec) + view_nd(spec.binary_mask()) From ef0426a1f8f15c43cba8ca762024e800245b49d0 Mon Sep 17 00:00:00 2001 From: Talley Lambert Date: Fri, 7 Jun 2024 13:25:37 -0400 Subject: [PATCH 3/7] save stuff --- src/microsim/allen/fetch.py | 6 ++---- x.py | 3 ++- 2 files changed, 4 insertions(+), 5 deletions(-) diff --git a/src/microsim/allen/fetch.py b/src/microsim/allen/fetch.py index c171f9f..1b9cc03 100644 --- a/src/microsim/allen/fetch.py +++ b/src/microsim/allen/fetch.py @@ -152,7 +152,7 @@ def binary_mask(self, voxel_size: float = 1, scale_factor: float = 3) -> np.ndar cz, cy, cx = child.shifted_coord(origin, voxel_size) draw_line_3d(px, py, pz, cx, cy, cz, grid, max_r=max_r, width=r) - soma_scale: float = 1.2 + soma_scale: float = 1 for comp in self._node_types[SWCType.SOMA]: z, y, x = comp.shifted_coord(origin, voxel_size) r = int(0.5 * scale_factor * soma_scale * comp.r) @@ -224,7 +224,7 @@ def fetch(cls, id: int) -> "Specimen": raise ValueError(qr.msg) return cast("Specimen", qr.msg[0]) - def binary_mask(self, voxel_size: float = 1, scale_factor: float = 1) -> np.ndarray: + def binary_mask(self, voxel_size: float = 1, scale_factor: float = 3) -> np.ndarray: """Return 3D binary mask for this specimen's neuron reconstructions.""" for recon in self.neuron_reconstructions: return recon.load_swc().binary_mask( @@ -239,8 +239,6 @@ class ApiCellTypesSpecimenDetail(BaseModel): structure__acronym: str | None donor__species: Literal["Homo Sapiens", "Mus musculus"] nr__reconstruction_type: str | None # probably just 'full' or 'dendrite-only' - nr__average_contraction: float | None - nr__average_parent_daughter_ratio: float | None nr__max_euclidean_distance: float | None nr__number_bifurcations: int | None nr__number_stems: int | None diff --git a/x.py b/x.py index 6cb620e..50335e8 100644 --- a/x.py +++ b/x.py @@ -5,5 +5,6 @@ for recon in get_reconstructions("Mus musculus", "full"): spec = recon.specimen() - print("viewing", recon, spec) + print("viewing", spec) view_nd(spec.binary_mask()) + break From ce9b44b9cb0648cd861af0845491146146241dcc Mon Sep 17 00:00:00 2001 From: Talley Lambert Date: Sun, 16 Jun 2024 17:26:11 -0400 Subject: [PATCH 4/7] more docs and refactor --- src/microsim/allen/__init__.py | 16 +- src/microsim/allen/_fetch.py | 230 +++++++++++++++++++++++++++ src/microsim/allen/_swc.py | 160 +++++++++++++++++++ src/microsim/allen/fetch.py | 280 --------------------------------- tests/test_allen.py | 0 5 files changed, 404 insertions(+), 282 deletions(-) create mode 100644 src/microsim/allen/_fetch.py create mode 100644 src/microsim/allen/_swc.py delete mode 100644 src/microsim/allen/fetch.py create mode 100644 tests/test_allen.py diff --git a/src/microsim/allen/__init__.py b/src/microsim/allen/__init__.py index d1b14d1..787bd05 100644 --- a/src/microsim/allen/__init__.py +++ b/src/microsim/allen/__init__.py @@ -1,3 +1,15 @@ -from .fetch import ApiCellTypesSpecimenDetail, Specimen, get_reconstructions +from ._fetch import ( + ApiCellTypesSpecimenDetail, + NeuronReconstruction, + Specimen, + get_reconstructions, +) +from ._swc import SWC -__all__ = ["get_reconstructions", "Specimen", "ApiCellTypesSpecimenDetail"] +__all__ = [ + "ApiCellTypesSpecimenDetail", + "get_reconstructions", + "NeuronReconstruction", + "Specimen", + "SWC", +] diff --git a/src/microsim/allen/_fetch.py b/src/microsim/allen/_fetch.py new file mode 100644 index 0000000..8f4f5b5 --- /dev/null +++ b/src/microsim/allen/_fetch.py @@ -0,0 +1,230 @@ +from __future__ import annotations + +from functools import cache, cached_property +from typing import TYPE_CHECKING, Literal, cast + +import requests +from pydantic import BaseModel, Field + +if TYPE_CHECKING: + from collections.abc import Iterable + + import numpy as np + + from ._swc import SWC + +ALLEN_ROOT = "http://api.brain-map.org" +ALLEN_V2_API = f"{ALLEN_ROOT}/api/v2/data" +ALLEN_V2_QUERY = ALLEN_V2_API + "/query.json" +SWC_FILE_TYPE = "3DNeuronReconstruction" + + +class WellKnownFileType(BaseModel): + """Model representing a well-known file type in the Allen Brain Map API.""" + + id: int + name: str # something like '3DNeuronReconstruction' + + +class WellKnownFile(BaseModel): + """Model representing a file in the Allen Brain Map API.""" + + attachable_id: int | None + attachable_type: str | None + download_link: str | None + id: int | None + path: str | None + well_known_file_type_id: int | None + well_known_file_type: WellKnownFileType | None + + +class NeuronReconstruction(BaseModel): + """Model representing a neuron reconstruction in the Allen Brain Map API.""" + + id: int + specimen_id: int + number_nodes: int + number_branches: int + number_stems: int + number_bifurcations: int + max_euclidean_distance: float + neuron_reconstruction_type: str + overall_height: float + overall_width: float + overall_depth: float + scale_factor_x: float + scale_factor_y: float + scale_factor_z: float + total_length: float + total_surface: float + total_volume: float + well_known_files: list[WellKnownFile] = Field(default_factory=list) + + @property + def swc_path(self) -> str: + """The SWC file for this reconstruction.""" + for f in self.well_known_files: + if ( + getattr(f.well_known_file_type, "name", None) == SWC_FILE_TYPE + and f.download_link + ): + return ALLEN_ROOT + f.download_link + raise ValueError("No SWC file found for this reconstruction.") + + @cached_property + def swc(self) -> SWC: + """Load the SWC file for this reconstruction.""" + from ._swc import SWC + + return SWC.from_path(self.swc_path) + + def binary_mask(self, voxel_size: float = 1, scale_factor: float = 3) -> np.ndarray: + """Return 3D binary mask for this neuron reconstructions.""" + return self.swc.binary_mask(voxel_size=voxel_size, scale_factor=scale_factor) + + @classmethod + @cache + def fetch(cls, id: int) -> NeuronReconstruction: + """Fetch NeuronReconstruction by ID from the Allen brain map API.""" + q = [ + "model::NeuronReconstruction", + f"rma::criteria[id$eq{id}],well_known_files", + f"rma::include,well_known_files(well_known_file_type[name$eq'{SWC_FILE_TYPE}'])", + # get all rows + "rma::options[num_rows$eq'all']", + ] + response = requests.get(ALLEN_V2_QUERY, params={"q": ",".join(q)}) + response.raise_for_status() + qr = _QueryResponse.model_validate_json(response.content) + if not qr.success: + raise ValueError(qr.msg) + return cast("NeuronReconstruction", qr.msg[0]) + + def specimen(self) -> Specimen: + """Fetch the specimen that owns this neuron reconstruction.""" + return Specimen.fetch(self.specimen_id) + + +class Structure(BaseModel): + """Speciment structure model from the Allen Brain Map API.""" + + id: int + name: str + acronym: str + structure_id_path: str + + +class Specimen(BaseModel): + """Model representing a specimen in the Allen Brain Map API.""" + + id: int + name: str + is_cell_specimen: bool + specimen_id_path: str + structure: Structure + neuron_reconstructions: list[NeuronReconstruction] = Field(default_factory=list) + + @classmethod + @cache + def fetch(cls, id: int) -> Specimen: + """Fetch this specimen from the Allen brain map API.""" + q = [ + # query the Specimen model + "model::Specimen", + # limit to the specimen with the given ID + # and join on NeuronReconstruction and WellKnownFile + f"rma::criteria[id$eq{id}],neuron_reconstructions(well_known_files)", + # include structure + # and neuron_reconstructions where the well_known_file_type is SWC + "rma::include,structure,neuron_reconstructions(well_known_files(" + f"well_known_file_type[name$eq'{SWC_FILE_TYPE}']))", + # get all rows + "rma::options[num_rows$eq'all']", + ] + response = requests.get(ALLEN_V2_QUERY, params={"q": ",".join(q)}) + response.raise_for_status() + qr = _QueryResponse.model_validate_json(response.content) + if not qr.success: + raise ValueError(qr.msg) + return cast("Specimen", qr.msg[0]) + + def binary_masks( + self, voxel_size: float = 1, scale_factor: float = 3 + ) -> list[np.ndarray]: + """Return all binary masks for this specimen's neuron reconstructions.""" + masks = [] + for recon in self.neuron_reconstructions: + masks.append( + recon.binary_mask(voxel_size=voxel_size, scale_factor=scale_factor) + ) + return masks + + @property + def url(self) -> str: + """Return the URL for this specimen on the Allen Brain Map.""" + return f"http://celltypes.brain-map.org/experiment/morphology/{self.id}" + + def open_webpage(self) -> None: + """Open the webpage for this specimen in the Allen Brain Map.""" + import webbrowser + + webbrowser.open(self.url) + + +class ApiCellTypesSpecimenDetail(BaseModel): + """Model representing Specimen details from the Allen Brain Map API.""" + + specimen__id: int + structure__name: str | None + structure__acronym: str | None + donor__species: Literal["Homo Sapiens", "Mus musculus"] + nr__reconstruction_type: str | None # probably just 'full' or 'dendrite-only' + nr__max_euclidean_distance: float | None + nr__number_bifurcations: int | None + nr__number_stems: int | None + + @classmethod + @cache + def all_reconstructions(cls) -> tuple[ApiCellTypesSpecimenDetail, ...]: + """Fetch details for all Specimens with reconstruction info.""" + q = ( + "model::ApiCellTypesSpecimenDetail", + # limit to specimens that have neuron reconstructions + "rma::criteria[nr__reconstruction_type$ne'null']", + # get all rows + "rma::options[num_rows$eq'all']", + ) + response = requests.get(ALLEN_V2_QUERY, params={"q": ",".join(q)}) + response.raise_for_status() + qr = _QueryResponse.model_validate_json(response.content) + if not qr.success: + raise ValueError(qr.msg) + return tuple(qr.msg) # type: ignore[arg-type] + + def specimen(self) -> Specimen: + """Return associated Specimen object.""" + return Specimen.fetch(self.specimen__id) + + +class _QueryResponse(BaseModel): + """Query response from the Allen Brain Map API.""" + + success: bool + msg: ( + list[NeuronReconstruction] + | list[Specimen] + | list[ApiCellTypesSpecimenDetail] + | str + ) + + +def get_reconstructions( + species: Literal["Homo Sapiens", "Mus musculus"] | None = None, + reconstruction_type: Literal["full", "dendrite-only"] | None = None, +) -> tuple[ApiCellTypesSpecimenDetail, ...]: + recons: Iterable = ApiCellTypesSpecimenDetail.all_reconstructions() + if species is not None: + recons = (x for x in recons if x.donor__species == species) + if reconstruction_type is not None: + recons = (x for x in recons if x.nr__reconstruction_type == reconstruction_type) + return tuple(recons) diff --git a/src/microsim/allen/_swc.py b/src/microsim/allen/_swc.py new file mode 100644 index 0000000..3f52f15 --- /dev/null +++ b/src/microsim/allen/_swc.py @@ -0,0 +1,160 @@ +from __future__ import annotations + +from collections import defaultdict +from enum import IntEnum +from functools import cached_property +from pathlib import Path +from typing import TYPE_CHECKING, NamedTuple + +import numpy as np +import requests + +if TYPE_CHECKING: + from collections.abc import Iterable, Iterator, Sequence + + +class SWCType(IntEnum): + """Constants for the SWC node types.""" + + UNDEFINED = 0 + SOMA = 1 + AXON = 2 + BASAL_DENDRITE = 3 + APICAL_DENDRITE = 4 + CUSTOM = 5 + UNSPECIFIED_NEURITE = 6 + GLIA_PROCESSES = 7 + + +class Compartment(NamedTuple): + """A compartment in an SWC file.""" + + id: int + t: int # type + x: float + y: float + z: float + r: float # radius + c: int # parent id or "connectivity" + + def coord(self) -> np.ndarray: + return np.array([self.z, self.y, self.x]) + + def shifted_coord( + self, origin: np.ndarray, resolution: float = 1 + ) -> tuple[int, int, int]: + shifted_coord = (self.coord() - origin) / resolution + return tuple(c for c in shifted_coord.astype(int)) + + +class SWC: + """Object representing an SWC file. + + Provides methods to parse SWC files and render them as binary masks. + """ + + @classmethod + def from_path(cls, path: str | Path) -> SWC: + if str(path).startswith(("http://", "https://")): + response = requests.get(str(path)) + response.raise_for_status() + content = response.text + else: + content = Path(path).expanduser().read_text() + return cls.from_string(content) + + @classmethod + def from_string(cls, content: str | bytes) -> SWC: + if isinstance(content, bytes): + content = content.decode() + + compartments = [] + for num, line in enumerate(content.splitlines()): + if line.startswith("#"): + continue + try: + a, b, c, d, e, f, g = line.split() + comp = Compartment( + int(a), int(b), float(c), float(d), float(e), float(f), int(g) + ) + except ValueError as e: + raise ValueError(f"Invalid SWC line {num}: {line}") from e + compartments.append(comp) + return cls(compartments) + + def __init__(self, compartments: Sequence[Compartment] = ()): + self.compartments = compartments + + self._id_map: dict[int, Compartment] = {} + self._children_of: defaultdict[int, list[Compartment]] = defaultdict(list) + self._node_types: defaultdict[int, list[Compartment]] = defaultdict(list) + for comp in compartments: + self._id_map[comp.id] = comp + self._children_of[comp.c].append(comp) + self._node_types[comp.t].append(comp) + + def iter_pairs(self, *types: int) -> Iterator[tuple[Compartment, Compartment]]: + seen = set() + for comp_id, children in self._children_of.items(): + if comp_id == -1: + continue + comp = self._id_map[comp_id] + if not types or comp.t in types: + for child in children: + if (comp, child) not in seen: + yield comp, child + seen.add((comp, child)) + + @cached_property + def coords(self) -> np.ndarray: + """Return the coordinates of the compartments as (M, 3) array.""" + return np.array([(c.z, c.y, c.x) for c in self.compartments]).astype("float32") + + def root(self) -> Compartment: + """Return the root compartment of the SWC.""" + return self._id_map[1] + + def origin(self) -> tuple[float, float, float]: + """Return the (Z,Y,X) coordinate of the root node.""" + root = self.root() + return root.z, root.y, root.x + + def _empty_grid(self, voxel_size: float = 1.0) -> np.ndarray: + """Create an empty 3D grid for the binary mask.""" + extent = np.ptp(self.coords, axis=0) + size = np.ceil(extent / voxel_size).astype(int) + return np.zeros(size, dtype=bool) + + def binary_mask( + self, + voxel_size: float = 1, + scale_factor: float = 3, + *, + include_types: Iterable[int] = ( + SWCType.BASAL_DENDRITE, + SWCType.APICAL_DENDRITE, + SWCType.AXON, + ), + ) -> np.ndarray: + """Render a binary mask of the neuron reconstruction.""" + from microsim._draw import draw_line_3d, draw_sphere + + grid = self._empty_grid(voxel_size) + origin = np.min(self.coords, axis=0) + + dend_scale: float = 1 + max_r = float(np.sum(grid.shape)) + + for par, child in self.iter_pairs(*include_types): + r = int(max(1, 0.5 * scale_factor * dend_scale * (par.r + child.r))) + pz, py, px = par.shifted_coord(origin, voxel_size) + cz, cy, cx = child.shifted_coord(origin, voxel_size) + draw_line_3d(px, py, pz, cx, cy, cz, grid, max_r=max_r, width=r) + + soma_scale: float = 1 + for comp in self._node_types[SWCType.SOMA]: + z, y, x = comp.shifted_coord(origin, voxel_size) + r = int(0.5 * scale_factor * soma_scale * comp.r) + draw_sphere(grid, x, y, z, r) + + return grid.astype(np.uint8) diff --git a/src/microsim/allen/fetch.py b/src/microsim/allen/fetch.py deleted file mode 100644 index 1b9cc03..0000000 --- a/src/microsim/allen/fetch.py +++ /dev/null @@ -1,280 +0,0 @@ -from collections import defaultdict -from collections.abc import Iterator, Sequence -from enum import IntEnum -from functools import cache, cached_property -from pathlib import Path -from typing import TYPE_CHECKING, Literal, NamedTuple, cast - -import numpy as np -import requests -from pydantic import BaseModel, Field - -if TYPE_CHECKING: - from collections.abc import Iterable - -ALLEN_ROOT = "http://api.brain-map.org" -ALLEN_V2_API = f"{ALLEN_ROOT}/api/v2/data" -ALLEN_V2_QUERY = ALLEN_V2_API + "/query.json" -SWC_FILE_TYPE = "3DNeuronReconstruction" - - -class WellKnownFileType(BaseModel): - id: int - name: str - - -class WellKnownFile(BaseModel): - attachable_id: int | None - attachable_type: str | None - download_link: str | None - id: int | None - path: str | None - well_known_file_type_id: int | None - well_known_file_type: WellKnownFileType | None - - -class SWCType(IntEnum): - UNDEFINED = 0 - SOMA = 1 - AXON = 2 - BASAL_DENDRITE = 3 - APICAL_DENDRITE = 4 - CUSTOM = 5 - UNSPECIFIED_NEURITE = 6 - GLIA_PROCESSES = 7 - - -class Compartment(NamedTuple): - id: int - t: int # type - x: float - y: float - z: float - r: float # radius - c: int # parent id or "connectivity" - - def coord(self) -> np.ndarray: - return np.array([self.z, self.y, self.x]) - - def shifted_coord( - self, origin: np.ndarray, resolution: float = 1 - ) -> tuple[int, int, int]: - shifted_coord = (self.coord() - origin) / resolution - return tuple(c for c in shifted_coord.astype(int)) - - -class SWC: - @classmethod - def from_path(cls, path: str | Path) -> "SWC": - if str(path).startswith(("http://", "https://")): - response = requests.get(str(path)) - response.raise_for_status() - content = response.text - else: - content = Path(path).expanduser().read_text() - return cls.from_string(content) - - @classmethod - def from_string(cls, content: str | bytes) -> "SWC": - if isinstance(content, bytes): - content = content.decode() - - compartments = [] - for num, line in enumerate(content.splitlines()): - if line.startswith("#"): - continue - try: - a, b, c, d, e, f, g = line.split() - comp = Compartment( - int(a), int(b), float(c), float(d), float(e), float(f), int(g) - ) - except ValueError as e: - raise ValueError(f"Invalid SWC line {num}: {line}") from e - compartments.append(comp) - return cls(compartments) - - def __init__(self, compartments: Sequence[Compartment] = ()): - self.compartments = compartments - - self._id_map: dict[int, Compartment] = {} - self._children_of: defaultdict[int, list[Compartment]] = defaultdict(list) - self._node_types: defaultdict[int, list[Compartment]] = defaultdict(list) - for comp in compartments: - self._id_map[comp.id] = comp - self._children_of[comp.c].append(comp) - self._node_types[comp.t].append(comp) - - def iter_pairs(self, *types: int) -> Iterator[tuple[Compartment, Compartment]]: - seen = set() - for comp_id, children in self._children_of.items(): - if comp_id == -1: - continue - comp = self._id_map[comp_id] - if not types or comp.t in types: - for child in children: - if (comp, child) not in seen: - yield comp, child - seen.add((comp, child)) - - @cached_property - def coords(self) -> np.ndarray: - """Return the coordinates of the compartments as (M, 3) array.""" - return np.array([(c.z, c.y, c.x) for c in self.compartments]).astype("float32") - - def root(self) -> Compartment: - """Return the root compartment of the SWC.""" - return self._id_map[1] - - def origin(self) -> tuple[float, float, float]: - """Return the (Z,Y,X) coordinate of the root node.""" - root = self.root() - return root.z, root.y, root.x - - def _empty_grid(self, voxel_size: float = 1.0) -> np.ndarray: - """Create an empty 3D grid for the binary mask.""" - extent = np.ptp(self.coords, axis=0) - size = np.ceil(extent / voxel_size).astype(int) - return np.zeros(size, dtype=bool) - - def binary_mask(self, voxel_size: float = 1, scale_factor: float = 3) -> np.ndarray: - from microsim._draw import draw_line_3d, draw_sphere - - grid = self._empty_grid(voxel_size) - origin = np.min(self.coords, axis=0) - - dend_scale: float = 1 - max_r = float(np.sum(grid.shape)) - for par, child in self.iter_pairs( - SWCType.BASAL_DENDRITE, SWCType.APICAL_DENDRITE - ): - r = int(max(1, 0.5 * scale_factor * dend_scale * (par.r + child.r))) - pz, py, px = par.shifted_coord(origin, voxel_size) - cz, cy, cx = child.shifted_coord(origin, voxel_size) - draw_line_3d(px, py, pz, cx, cy, cz, grid, max_r=max_r, width=r) - - soma_scale: float = 1 - for comp in self._node_types[SWCType.SOMA]: - z, y, x = comp.shifted_coord(origin, voxel_size) - r = int(0.5 * scale_factor * soma_scale * comp.r) - draw_sphere(grid, x, y, z, r) - - return grid.astype(np.uint8) - - -class NeuronReconstruction(BaseModel): - id: int - specimen_id: int - number_nodes: int - number_branches: int - neuron_reconstruction_type: str - overall_height: float - overall_width: float - overall_depth: float - scale_factor_x: float - scale_factor_y: float - scale_factor_z: float - well_known_files: list[WellKnownFile] = Field(default_factory=list) - - @property - def swc_path(self) -> str: - """The SWC file for this reconstruction.""" - for f in self.well_known_files: - if ( - getattr(f.well_known_file_type, "name", None) == SWC_FILE_TYPE - and f.download_link - ): - return ALLEN_ROOT + f.download_link - raise ValueError("No SWC file found for this reconstruction.") - - def load_swc(self) -> SWC: - """Load the SWC file for this reconstruction.""" - return SWC.from_path(self.swc_path) - - -class Structure(BaseModel): - id: int - name: str - acronym: str - structure_id_path: str - - -class Specimen(BaseModel): - id: int - name: str - is_cell_specimen: bool - specimen_id_path: str - structure: Structure - neuron_reconstructions: list[NeuronReconstruction] = Field(default_factory=list) - - @classmethod - @cache - def fetch(cls, id: int) -> "Specimen": - """Fetch this specimen from the Allen brain map API.""" - q = [ - "model::Specimen", - f"rma::criteria[id$eq{id}],neuron_reconstructions(well_known_files)", - "rma::include,structure,neuron_reconstructions(well_known_files(" - f"well_known_file_type[name$eq'{SWC_FILE_TYPE}']))", - "rma::options[num_rows$eq'all']", - ] - response = requests.get(ALLEN_V2_QUERY, params={"q": ",".join(q)}) - response.raise_for_status() - qr = QueryResponse.model_validate_json(response.content) - if not qr.success: - raise ValueError(qr.msg) - return cast("Specimen", qr.msg[0]) - - def binary_mask(self, voxel_size: float = 1, scale_factor: float = 3) -> np.ndarray: - """Return 3D binary mask for this specimen's neuron reconstructions.""" - for recon in self.neuron_reconstructions: - return recon.load_swc().binary_mask( - voxel_size=voxel_size, scale_factor=scale_factor - ) - raise ValueError("No neuron reconstructions found for this specimen.") - - -class ApiCellTypesSpecimenDetail(BaseModel): - specimen__id: int - structure__name: str | None - structure__acronym: str | None - donor__species: Literal["Homo Sapiens", "Mus musculus"] - nr__reconstruction_type: str | None # probably just 'full' or 'dendrite-only' - nr__max_euclidean_distance: float | None - nr__number_bifurcations: int | None - nr__number_stems: int | None - - def specimen(self) -> Specimen: - return Specimen.fetch(self.specimen__id) - - @classmethod - @cache - def all_reconstructions(cls) -> tuple["ApiCellTypesSpecimenDetail", ...]: - """Fetch details for all Specimens with reconstruction info.""" - q = ( - "model::ApiCellTypesSpecimenDetail", - "rma::criteria[nr__reconstruction_type$ne'null']", - "rma::options[num_rows$eq'all']", - ) - response = requests.get(ALLEN_V2_QUERY, params={"q": ",".join(q)}) - response.raise_for_status() - qr = QueryResponse.model_validate_json(response.content) - if not qr.success: - raise ValueError(qr.msg) - return tuple(qr.msg) # type: ignore[arg-type] - - -class QueryResponse(BaseModel): - success: bool - msg: list[Specimen] | list[ApiCellTypesSpecimenDetail] | str - - -def get_reconstructions( - species: Literal["Homo Sapiens", "Mus musculus"] | None, - reconstruction_type: Literal["full", "dendrite-only"] | None = None, -) -> tuple[ApiCellTypesSpecimenDetail, ...]: - recons: Iterable = ApiCellTypesSpecimenDetail.all_reconstructions() - if species is not None: - recons = (x for x in recons if x.donor__species == species) - if reconstruction_type is not None: - recons = (x for x in recons if x.nr__reconstruction_type == reconstruction_type) - return tuple(recons) diff --git a/tests/test_allen.py b/tests/test_allen.py new file mode 100644 index 0000000..e69de29 From 60bf868f61d95031f94291bef242be20e4dfbb17 Mon Sep 17 00:00:00 2001 From: Talley Lambert Date: Sun, 16 Jun 2024 18:03:17 -0400 Subject: [PATCH 5/7] add tests --- src/microsim/allen/_fetch.py | 14 +++++++------- tests/test_allen.py | 23 +++++++++++++++++++++++ x.py | 10 ---------- 3 files changed, 30 insertions(+), 17 deletions(-) delete mode 100644 x.py diff --git a/src/microsim/allen/_fetch.py b/src/microsim/allen/_fetch.py index 8f4f5b5..e12af24 100644 --- a/src/microsim/allen/_fetch.py +++ b/src/microsim/allen/_fetch.py @@ -69,7 +69,9 @@ def swc_path(self) -> str: and f.download_link ): return ALLEN_ROOT + f.download_link - raise ValueError("No SWC file found for this reconstruction.") + raise ValueError( + "No SWC file found for this reconstruction." + ) # pragma: no cover @cached_property def swc(self) -> SWC: @@ -96,7 +98,7 @@ def fetch(cls, id: int) -> NeuronReconstruction: response = requests.get(ALLEN_V2_QUERY, params={"q": ",".join(q)}) response.raise_for_status() qr = _QueryResponse.model_validate_json(response.content) - if not qr.success: + if not qr.success: # pragma: no cover raise ValueError(qr.msg) return cast("NeuronReconstruction", qr.msg[0]) @@ -144,7 +146,7 @@ def fetch(cls, id: int) -> Specimen: response = requests.get(ALLEN_V2_QUERY, params={"q": ",".join(q)}) response.raise_for_status() qr = _QueryResponse.model_validate_json(response.content) - if not qr.success: + if not qr.success: # pragma: no cover raise ValueError(qr.msg) return cast("Specimen", qr.msg[0]) @@ -164,7 +166,7 @@ def url(self) -> str: """Return the URL for this specimen on the Allen Brain Map.""" return f"http://celltypes.brain-map.org/experiment/morphology/{self.id}" - def open_webpage(self) -> None: + def open_webpage(self) -> None: # pragma: no cover """Open the webpage for this specimen in the Allen Brain Map.""" import webbrowser @@ -189,15 +191,13 @@ def all_reconstructions(cls) -> tuple[ApiCellTypesSpecimenDetail, ...]: """Fetch details for all Specimens with reconstruction info.""" q = ( "model::ApiCellTypesSpecimenDetail", - # limit to specimens that have neuron reconstructions "rma::criteria[nr__reconstruction_type$ne'null']", - # get all rows "rma::options[num_rows$eq'all']", ) response = requests.get(ALLEN_V2_QUERY, params={"q": ",".join(q)}) response.raise_for_status() qr = _QueryResponse.model_validate_json(response.content) - if not qr.success: + if not qr.success: # pragma: no cover raise ValueError(qr.msg) return tuple(qr.msg) # type: ignore[arg-type] diff --git a/tests/test_allen.py b/tests/test_allen.py index e69de29..b2cf159 100644 --- a/tests/test_allen.py +++ b/tests/test_allen.py @@ -0,0 +1,23 @@ +import numpy as np + +from microsim.allen import NeuronReconstruction, Specimen + + +def test_neuron_reconstruction(): + nr = NeuronReconstruction.fetch(638976782) + + mask = nr.binary_mask() + assert isinstance(mask, np.ndarray) + assert mask.ndim == 3 + + swc = nr.swc + root = swc.root() + assert root.id == 1 + assert swc.origin() == (root.z, root.y, root.x) + + +def test_specimen(): + spec = Specimen.fetch(586073850) + masks = spec.binary_masks() + assert len(masks) == len(spec.neuron_reconstructions) + assert isinstance(spec.url, str) diff --git a/x.py b/x.py deleted file mode 100644 index 50335e8..0000000 --- a/x.py +++ /dev/null @@ -1,10 +0,0 @@ -from rich import print - -from microsim.allen import get_reconstructions -from microsim.util import view_nd - -for recon in get_reconstructions("Mus musculus", "full"): - spec = recon.specimen() - print("viewing", spec) - view_nd(spec.binary_mask()) - break From 1d94993986c6ded93494c17dec9a6751f60f763b Mon Sep 17 00:00:00 2001 From: Talley Lambert Date: Sun, 16 Jun 2024 18:16:10 -0400 Subject: [PATCH 6/7] don't require requests --- src/microsim/allen/_fetch.py | 18 ++++++++---------- src/microsim/allen/_swc.py | 8 ++++---- src/microsim/util.py | 15 +++++++++++++++ 3 files changed, 27 insertions(+), 14 deletions(-) diff --git a/src/microsim/allen/_fetch.py b/src/microsim/allen/_fetch.py index e12af24..ed1beb8 100644 --- a/src/microsim/allen/_fetch.py +++ b/src/microsim/allen/_fetch.py @@ -3,9 +3,10 @@ from functools import cache, cached_property from typing import TYPE_CHECKING, Literal, cast -import requests from pydantic import BaseModel, Field +from microsim.util import http_get + if TYPE_CHECKING: from collections.abc import Iterable @@ -95,9 +96,8 @@ def fetch(cls, id: int) -> NeuronReconstruction: # get all rows "rma::options[num_rows$eq'all']", ] - response = requests.get(ALLEN_V2_QUERY, params={"q": ",".join(q)}) - response.raise_for_status() - qr = _QueryResponse.model_validate_json(response.content) + response = http_get(ALLEN_V2_QUERY, params={"q": ",".join(q)}) + qr = _QueryResponse.model_validate_json(response) if not qr.success: # pragma: no cover raise ValueError(qr.msg) return cast("NeuronReconstruction", qr.msg[0]) @@ -143,9 +143,8 @@ def fetch(cls, id: int) -> Specimen: # get all rows "rma::options[num_rows$eq'all']", ] - response = requests.get(ALLEN_V2_QUERY, params={"q": ",".join(q)}) - response.raise_for_status() - qr = _QueryResponse.model_validate_json(response.content) + response = http_get(ALLEN_V2_QUERY, params={"q": ",".join(q)}) + qr = _QueryResponse.model_validate_json(response) if not qr.success: # pragma: no cover raise ValueError(qr.msg) return cast("Specimen", qr.msg[0]) @@ -194,9 +193,8 @@ def all_reconstructions(cls) -> tuple[ApiCellTypesSpecimenDetail, ...]: "rma::criteria[nr__reconstruction_type$ne'null']", "rma::options[num_rows$eq'all']", ) - response = requests.get(ALLEN_V2_QUERY, params={"q": ",".join(q)}) - response.raise_for_status() - qr = _QueryResponse.model_validate_json(response.content) + response = http_get(ALLEN_V2_QUERY, params={"q": ",".join(q)}) + qr = _QueryResponse.model_validate_json(response) if not qr.success: # pragma: no cover raise ValueError(qr.msg) return tuple(qr.msg) # type: ignore[arg-type] diff --git a/src/microsim/allen/_swc.py b/src/microsim/allen/_swc.py index 3f52f15..d471083 100644 --- a/src/microsim/allen/_swc.py +++ b/src/microsim/allen/_swc.py @@ -7,7 +7,8 @@ from typing import TYPE_CHECKING, NamedTuple import numpy as np -import requests + +from microsim.util import http_get if TYPE_CHECKING: from collections.abc import Iterable, Iterator, Sequence @@ -56,9 +57,8 @@ class SWC: @classmethod def from_path(cls, path: str | Path) -> SWC: if str(path).startswith(("http://", "https://")): - response = requests.get(str(path)) - response.raise_for_status() - content = response.text + response = http_get(str(path)) + content = response.decode() else: content = Path(path).expanduser().read_text() return cls.from_string(content) diff --git a/src/microsim/util.py b/src/microsim/util.py index 8af29ed..2db8e8a 100644 --- a/src/microsim/util.py +++ b/src/microsim/util.py @@ -4,6 +4,8 @@ import shutil import warnings from typing import TYPE_CHECKING, Any, Protocol, TypeVar, cast +from urllib import parse, request +from urllib.error import HTTPError import numpy as np import numpy.typing as npt @@ -354,3 +356,16 @@ def norm_name(name: str) -> str: for char in " -/\\()[],;:!?@#$%^&*+=|<>'\"": name = name.replace(char, "_") return name + + +def http_get(url: str, params: dict | None = None) -> bytes: + """API like requests.get but with standard-library urllib.""" + if params: + url += "?" + parse.urlencode(params) + + with request.urlopen(url) as response: + if not 200 <= response.getcode() < 300: + raise HTTPError( + url, response.getcode(), "HTTP request failed", response.headers, None + ) + return response.read() From 552421e23ea1a8d7d91cf668bfd39a4b03d670f6 Mon Sep 17 00:00:00 2001 From: Talley Lambert Date: Sun, 16 Jun 2024 18:42:47 -0400 Subject: [PATCH 7/7] lint --- src/microsim/util.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/microsim/util.py b/src/microsim/util.py index 2db8e8a..586084c 100644 --- a/src/microsim/util.py +++ b/src/microsim/util.py @@ -368,4 +368,4 @@ def http_get(url: str, params: dict | None = None) -> bytes: raise HTTPError( url, response.getcode(), "HTTP request failed", response.headers, None ) - return response.read() + return cast(bytes, response.read())