diff --git a/CHANGELOG.md b/CHANGELOG.md index 2ee6caa..6ec52bc 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -4,6 +4,20 @@ All notable changes to this project will be documented in this file. The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.1.0/), and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.html). +## [0.3.0] - 2024-04-15 + +### Changed +__Locus spots visualisation__ +* QC-failed points are now... + * only text, no actual point + * sky blue rather than standard blue +* QC-passed points are now... + * goldenrod/dandelion yellow rather than red +* Color changes have been made to be more [colorblind-friendly](https://davidmathlogic.com/colorblind/). +* Points (or text) are now visible throughout each $z$-stack, rather than just in the slice nearest the center of the Gaussian fit. +* Each QC-pass point's position in $xy$ is shown as a stars in the $z$-slice corresponding to the truncation of the $z$-coordinate of the center of its Gaussian fit; it's shown in other slices as a circle. +* Each QC-pass point is 50% larger in the $z$ slice corresponding to the center of its fit. + ## [0.2.0] - 2024-04-14 ### Added diff --git a/docs/user_docs/locus-spots.md b/docs/user_docs/locus-spots.md index 7fc01de..2650936 100644 --- a/docs/user_docs/locus-spots.md +++ b/docs/user_docs/locus-spots.md @@ -2,7 +2,7 @@ ### Quickstart To visualise the locus-specific spots masks, you need to __drag-and-drop files__ into an active Napari window. -1. First, drag one of the `PXXXX.zarr` files--from the `locus_spots_visualisation` subfolder in your analysis folder from `looptrace`--into Napari. +1. First, drag one of the `PXXXX.zarr` files--from the `locus_spots_visualisation` (previously `locus_spot_images`) subfolder in your analysis folder from `looptrace`--into Napari. 1. Choose to open the ZARR with the Napari native plugin(s) / readers, not this plugin's. 1. Select "continuous" for the value of the "auto-contrast" option in the upper-left panel of the Napari window. 1. Second, drag a the corresponding (by field of view) pair of `*qcpass.csv` and `*qcfail.csv` files into Napari. @@ -22,4 +22,12 @@ This can makes scrolling time-consuming and will be simplified in a future relea __Finding points (centers of Gaussian fits to pixel volume)__ +In version 0.3 of this plugin, points should be visible throughout each z-stack, and the color of point which passed QC has been changed from red to yellow. +The shape is an star in the $z$ slice closest to the center of the Gaussian fit to a pixel volume, and a circle in other $z$ slices. +In this version, then, shape relates to $z$ position rather than to QC status; we allow point color to disambiguate QC status. + In version 0.2 of this plugin, a point (blue or red indicator) can only be visible only in the z-slice closest to the z-coordinate of the the center of the Gaussian fit to a particular pixel volume. Generally, these will be more toward the middle of the z-stack rather than at either end, so focus scrolling mainly through the midrange of the z-axis. + +### Frequently asked questions (FAQ) +1. Why do point centers appear _slightly_ off (thinking about subpixel resolution)?\ + We think this comes from the fact that Napari is regarding the _center_ of the upper-left pixel as the origin, rather than the upper-left corner of the upper-left pixel. diff --git a/looptrace_napari/__init__.py b/looptrace_napari/__init__.py index 0b00e64..a60cda7 100644 --- a/looptrace_napari/__init__.py +++ b/looptrace_napari/__init__.py @@ -1,3 +1,3 @@ """Package-level attributes and functions""" -__version__ = "0.2.0" +__version__ = "0.3.0" diff --git a/looptrace_napari/_colors.py b/looptrace_napari/_colors.py new file mode 100644 index 0000000..96fe653 --- /dev/null +++ b/looptrace_napari/_colors.py @@ -0,0 +1,5 @@ +"""Color tools""" + +# See: https://davidmathlogic.com/colorblind/ +DEEP_SKY_BLUE = "#0C7BDC" +GOLDENROD = "#FFC20A" diff --git a/looptrace_napari/geometry.py b/looptrace_napari/geometry.py index 1f9e8dc..bee2d8a 100644 --- a/looptrace_napari/geometry.py +++ b/looptrace_napari/geometry.py @@ -2,13 +2,12 @@ from abc import abstractmethod from dataclasses import dataclass -from typing import Protocol +from typing import Protocol, Union +import numpy as np from numpydoc_decorator import doc # type: ignore[import-untyped] -# TODO: need Python >= 3.12 -# See: https://github.com/gerlichlab/looptrace-napari/issues/6 -# from typing import override +ZCoordinate = Union[int, float, np.float64] # int to accommodate notion of "z-slice" class LocatableXY(Protocol): @@ -23,6 +22,14 @@ def get_y_coordinate(self) -> float: raise NotImplementedError +class LocatableZ(Protocol): + """Something that admits z-coordinate.""" + + @abstractmethod + def get_z_coordinate(self) -> ZCoordinate: + raise NotImplementedError + + @doc( summary="Bundle x and y position to create point in 2D space.", parameters=dict( @@ -41,13 +48,34 @@ def __post_init__(self) -> None: if any(c < 0 for c in [self.x, self.y]): raise ValueError(f"At least one coordinate is negative! {self}") - # TODO: adopt @override once on Python >= 3.12 - # See: https://github.com/gerlichlab/looptrace-napari/issues/6 def get_x_coordinate(self) -> float: return self.x - # TODO: adopt @override once on Python >= 3.12 - # See: https://github.com/gerlichlab/looptrace-napari/issues/6 - # @override def get_y_coordinate(self) -> float: return self.y + + +@doc( + summary="Bundle x and y position to create point in 2D space.", + parameters=dict( + x="Position in x", + y="Position in y", + z="Position in z", + ), + see_also=dict( + ImagePoint2D="Simpler, non-z implementation of an image point", + ), +) +@dataclass(kw_only=True, frozen=True) +class ImagePoint3D(ImagePoint2D, LocatableZ): + z: ZCoordinate + + def __post_init__(self) -> None: + super().__post_init__() + if not isinstance(self.z, (int, float, np.float64)): + raise TypeError(f"Bad z ({type(self.z).__name__}: {self.z}") + if any(c < 0 for c in [self.x, self.y]): + raise ValueError(f"z-coordinate is negative! {self}") + + def get_z_coordinate(self) -> ZCoordinate: + return self.z diff --git a/looptrace_napari/locus_specific_points_reader.py b/looptrace_napari/locus_specific_points_reader.py index f8c843c..8283ed5 100644 --- a/looptrace_napari/locus_specific_points_reader.py +++ b/looptrace_napari/locus_specific_points_reader.py @@ -1,32 +1,31 @@ """Reading locus-specific spots and points data from looptrace for visualisation in napari""" import csv +import dataclasses import logging import os +from enum import Enum +from math import floor from pathlib import Path -from typing import Callable, Literal, Optional, Tuple, Union +from typing import Callable, Dict, Literal, Optional, Tuple, Union +import numpy as np from numpydoc_decorator import doc # type: ignore[import-untyped] +from ._colors import DEEP_SKY_BLUE, GOLDENROD from ._docs import CommonReaderDoc -from .types import ( # pylint: disable=unused-import - CsvRow, - LayerParams, - PathLike, - PointRecord, - Timepoint, - TraceId, -) +from .geometry import ImagePoint3D, LocatableXY, LocatableZ, ZCoordinate +from .types import CsvRow, LayerParams, PathLike +from .types import TimepointFrom0 as Timepoint # pylint: disable=unused-import +from .types import TraceIdFrom0 as TraceId __author__ = "Vince Reuter" __credits__ = ["Vince Reuter"] -FullLayerData = Tuple["RawLocusPointsLike", "LayerParams", "LayerTypeName"] +FlatPointRecord = list[Union[float, ZCoordinate]] +FullLayerData = Tuple[list[FlatPointRecord], "LayerParams", "LayerTypeName"] LayerTypeName = Literal["points"] QCFailReasons = str -RawLocusPointsLike = list[Union["TraceId", "Timepoint", float]] - -QC_FAIL_CODES_KEY = "failCodes" @doc( @@ -41,8 +40,7 @@ def get_reader( path: CommonReaderDoc.path_type, ) -> Optional[Callable[[PathLike], list[FullLayerData]]]: static_params = { - "size": 0.5, - "edge_width": 0.05, + "edge_width": 0.1, "edge_width_is_relative": True, "n_dimensional": False, } @@ -67,67 +65,98 @@ def do_not_parse(why: str, *, level: int = logging.DEBUG) -> None: # Determine how to read and display the points layer to be parsed. status_name = base.lower().split(".")[-1].lstrip("qc_").lstrip("qc") if status_name in {"pass", "passed"}: - logging.debug("Parsing as QC-pass: %s", path) - color = "red" - symbol = "*" + logging.debug("Will parse sas QC-pass: %s", path) + color = GOLDENROD read_rows = parse_passed elif status_name in {"fail", "failed"}: - logging.debug("Parsing as QC-fail: %s", path) - color = "blue" - symbol = "o" + logging.debug("Will parse as QC-fail: %s", path) + color = DEEP_SKY_BLUE read_rows = parse_failed else: return do_not_parse( # type: ignore[func-returns-value] f"Could not infer QC status from '{status_name}'", level=logging.WARNING ) + base_meta = {"edge_color": color, "face_color": color} # Build the actual parser function. - base_meta = {"edge_color": color, "face_color": color, "symbol": symbol} - def parse(p): with open(p, mode="r", newline="") as fh: rows = list(csv.reader(fh)) - data, extra_meta = read_rows(rows) - if not data: + point_records, center_flags, extra_meta = read_rows(rows) + if not point_records: logging.warning("No data rows parsed!") - params = {**static_params, **base_meta, **extra_meta} - return data, params, "points" + shape_meta = { + "symbol": ["*" if is_center else "o" for is_center in center_flags], + } + params = {**static_params, **base_meta, **extra_meta, **shape_meta} + return [pt_rec.flatten() for pt_rec in point_records], params, "points" return lambda p: [parse(p)] @doc( summary="Parse records from points which passed QC.", - parameters=dict(records="Records to parse"), + parameters=dict(rows="Records to parse"), returns=""" A pair in which the first element is the array-like of points coordinates, and the second element is the mapping from attribute name to list of values (1 per point). """, notes="https://napari.org/stable/plugins/guides.html#layer-data-tuples", ) -def parse_passed(records: list[CsvRow]) -> Tuple[RawLocusPointsLike, LayerParams]: - return [parse_simple_record(r, exp_num_fields=5) for r in records], {} +def parse_passed( + rows: list[CsvRow], +) -> Tuple[list["PointRecord"], list[bool], LayerParams]: + records = [parse_simple_record(r, exp_num_fields=5) for r in rows] + max_z = max(r.get_z_coordinate() for r in records) + points: list["PointRecord"] = [] + center_flags: list[bool] = [] + for rec in records: + new_points, new_flags = expand_along_z(rec, z_max=max_z) + points.extend(new_points) + center_flags.extend(new_flags) + sizes = [1.5 if is_center else 1.0 for is_center in center_flags] + return points, center_flags, {"size": sizes} @doc( summary="Parse records from points which failed QC.", - parameters=dict(records="Records to parse"), + parameters=dict(rows="Records to parse"), returns=""" A pair in which the first element is the array-like of points coordinates, and the second element is the mapping from attribute name to list of values (1 per point). """, notes="https://napari.org/stable/plugins/guides.html#layer-data-tuples", ) -def parse_failed(records: list[CsvRow]) -> Tuple[RawLocusPointsLike, LayerParams]: - if not records: - data = [] - codes = [] - else: - data_codes_pairs: list[Tuple[PointRecord, QCFailReasons]] = [ - (parse_simple_record(r, exp_num_fields=6), r[5]) for r in records - ] - data, codes = map(list, zip(*data_codes_pairs)) - return data, {"text": QC_FAIL_CODES_KEY, "properties": {QC_FAIL_CODES_KEY: codes}} +def parse_failed( + rows: list[CsvRow], +) -> Tuple[list["PointRecord"], list[bool], LayerParams]: + record_qc_pairs: list[Tuple[PointRecord, QCFailReasons]] = [] + for row in rows: + try: + qc = row[InputFileColumn.QC.get] + rec = parse_simple_record(row, exp_num_fields=6) + except IndexError: + logging.error("Bad row: %s", row) + raise + record_qc_pairs.append((rec, qc)) + max_z = max(r.get_z_coordinate() for r, _ in record_qc_pairs) + points: list["PointRecord"] = [] + center_flags: list[bool] = [] + codes: list[QCFailReasons] = [] + for rec, qc in record_qc_pairs: + new_points, new_flags = expand_along_z(rec, z_max=max_z) + points.extend(new_points) + center_flags.extend(new_flags) + codes.extend([qc] * len(new_points)) + params = { + "size": 0, # Make the point invisible and just use text. + "text": { + "string": "{failCodes}", + "color": DEEP_SKY_BLUE, + }, + "properties": {"failCodes": codes}, + } + return points, center_flags, params @doc( @@ -143,7 +172,7 @@ def parse_failed(records: list[CsvRow]) -> Tuple[RawLocusPointsLike, LayerParams and the second element represents the (z, y, x) coordinates of the centroid of the spot fit. """, ) -def parse_simple_record(r: CsvRow, *, exp_num_fields: int) -> RawLocusPointsLike: +def parse_simple_record(r: CsvRow, *, exp_num_fields: int) -> "PointRecord": """Parse a single line from an input CSV file.""" if not isinstance(r, list): raise TypeError(f"Record to parse must be list, not {type(r).__name__}") @@ -151,9 +180,114 @@ def parse_simple_record(r: CsvRow, *, exp_num_fields: int) -> RawLocusPointsLike raise ValueError( f"Expected record of length {exp_num_fields} but got {len(r)}: {r}" ) - trace = int(r[0]) - timepoint = int(r[1]) - z = float(r[2]) - y = float(r[3]) - x = float(r[4]) - return [trace, timepoint, z, y, x] + trace = TraceId(int(r[InputFileColumn.TRACE.get])) + timepoint = Timepoint(int(r[InputFileColumn.TIMEPOINT.get])) + z = float(r[InputFileColumn.Z.get]) + y = float(r[InputFileColumn.Y.get]) + x = float(r[InputFileColumn.X.get]) + point = ImagePoint3D(z=z, y=y, x=x) + return PointRecord(trace_id=trace, timepoint=timepoint, point=point) + + +@dataclasses.dataclass(frozen=True, kw_only=True) +class PointRecord(LocatableXY, LocatableZ): + trace_id: TraceId + timepoint: Timepoint + point: ImagePoint3D + + def __post_init__(self) -> None: + bads: Dict[str, object] = {} + if not isinstance(self.trace_id, TraceId): + bads["trace ID"] = self.trace_id + if not isinstance(self.timepoint, Timepoint): + bads["timepoint"] = self.timepoint + if not isinstance(self.point, ImagePoint3D): + bads["point"] = self.point + if bads: + messages = "; ".join( + f"Bad type ({type(v).__name__}) for {k}" for k, v in bads.items() + ) + raise TypeError(f"Cannot create point record: {messages}") + + @doc(summary="Flatten") + def flatten(self) -> FlatPointRecord: + return [ + self.trace_id.get, + self.timepoint.get, + self.get_z_coordinate(), + self.get_y_coordinate(), + self.get_x_coordinate(), + ] + + def get_x_coordinate(self) -> float: + return self.point.x + + def get_y_coordinate(self) -> float: + return self.point.y + + def get_z_coordinate(self) -> ZCoordinate: + return self.point.z + + @doc(summary="Round point position to nearest z-slice") + def with_truncated_z(self) -> "PointRecord": + return self.with_new_z(floor(self.get_z_coordinate())) + + @doc( + summary="Replace this instance's point with a copy with updated z.", + parameters=dict(z="New z-coordinate value"), + ) + def with_new_z(self, z: int) -> "PointRecord": + pt = ImagePoint3D(x=self.point.x, y=self.point.y, z=z) + return dataclasses.replace(self, point=pt) + + +@doc( + summary="Create ancillary points from main point", + parameters=dict( + r="The record to expand along z-axis", + z_max="The maximum z-coordinate", + ), +) +def expand_along_z( + r: PointRecord, *, z_max: Union[float, np.float64] +) -> Tuple[list[PointRecord], list[bool]]: + if not isinstance(z_max, (int, float, np.float64)): + raise TypeError(f"Bad type for z_max: {type(z_max).__name__}") + + r = r.with_truncated_z() # type: ignore[no-redef] + z_center = int(r.get_z_coordinate()) + z_max = int(floor(z_max)) # type: ignore[no-redef] + assert isinstance(z_center, int) and isinstance(z_max, int), ( + f"Z center and Z max must be int; got {type(z_center).__name__} and" + f" {type(z_max).__name__}" + ) + + # Check that max z and center z make sense together. + if z_max < z_center: + raise ValueError(f"Max z must be at least as great as central z ({z_center})") + + # Build the records and flags of where the center in z really is. + predecessors = [(r.with_new_z(i), False) for i in range(z_center)] + successors = [(r.with_new_z(i), False) for i in range(z_center + 1, z_max + 1)] + points, params = zip(*(predecessors + [(r, True)] + successors)) + + # Each record should give rise to a total of 1 + z_max records, since numbering from 0. + assert ( + len(points) == 1 + z_max + ), f"Point={r}, z_max={z_max}, len(points)={len(points)}" + return points, params # type: ignore[return-value] + + +class InputFileColumn(Enum): + """Indices of the different columns to parse as particular fields""" + + TRACE = 0 + TIMEPOINT = 1 + Z = 2 + Y = 3 + X = 4 + QC = 5 + + @property + def get(self) -> int: + return self.value diff --git a/looptrace_napari/types.py b/looptrace_napari/types.py index e61f501..9ce53a0 100644 --- a/looptrace_napari/types.py +++ b/looptrace_napari/types.py @@ -2,22 +2,14 @@ from dataclasses import dataclass from pathlib import Path -from typing import Dict, Tuple, Union +from typing import Dict, Union from numpydoc_decorator import doc # type: ignore[import-untyped] CsvRow = list[str] -FailCodesText = str LayerParams = Dict PathLike = Union[str, Path] PathOrPaths = Union[PathLike, list[PathLike]] -PointId = Tuple["TraceId", "Timepoint"] -PointRecord = Tuple[PointId, "Point3D"] -Point3D = Tuple[float, float, float] -QCFailRecord = Tuple[PointId, Point3D, FailCodesText] -QCPassRecord = PointRecord -TraceId = int -Timepoint = int @doc( @@ -33,7 +25,13 @@ class FieldOfViewFrom1: get: int def __post_init__(self) -> None: - _refine_as_pos_int(self.get, context="FOV") + if not isinstance(self.get, int): + raise TypeError( + f"Non-integer as 1-based FOV! {self.get} (type" + f" {type(self.get).__name__})" + ) + if self.get < 1: + raise ValueError(f"1-based FOV view must be positive int; got {self.get}") @doc( @@ -49,11 +47,52 @@ class NucleusNumber: get: int def __post_init__(self) -> None: - _refine_as_pos_int(self.get, context="nucleus number") + if not isinstance(self.get, int): + raise TypeError( + f"Non-integer as nucleus number! {self.get} (type" + f" {type(self.get).__name__})" + ) + if self.get < 1: + raise ValueError(f"Nucleus number must be positive int; got {self.get}") -def _refine_as_pos_int(x: object, *, context: str) -> None: - if not isinstance(x, int): - raise TypeError(f"Non-integer as {context}! {x} (type {type(x).__name__})") - if x < 1: - raise ValueError(f"1-based {context} must be positive int, not {x}") +@doc( + summary="Wrap an int as a 0-based timepoint.", + parameters=dict(get="The value to wrap"), + raises=dict( + TypeError="If the given value to wrap isn't an integer", + ValueError="If the given value to wrap is negative", + ), +) +@dataclass(frozen=True, order=True) +class TimepointFrom0: + get: int + + def __post_init__(self) -> None: + if not isinstance(self.get, int): + raise TypeError( + f"Non-integer as timepoint! {self.get} (type {type(self.get).__name__})" + ) + if self.get < 0: + raise ValueError(f"Timepoint must be nonnegative int; got {self.get}") + + +@doc( + summary="Wrap an int as a 0-based trace ID.", + parameters=dict(get="The value to wrap"), + raises=dict( + TypeError="If the given value to wrap isn't an integer", + ValueError="If the given value to wrap is negative", + ), +) +@dataclass(frozen=True, order=True) +class TraceIdFrom0: + get: int + + def __post_init__(self) -> None: + if not isinstance(self.get, int): + raise TypeError( + f"Non-integer as trace ID! {self.get} (type {type(self.get).__name__})" + ) + if self.get < 0: + raise ValueError(f"Trace ID must be nonnegative int; got {self.get}") diff --git a/pyproject.toml b/pyproject.toml index 12786af..87bdc7a 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -89,5 +89,5 @@ profile = "black" [tool.pylint] # Disable missing docstrings since in general they're handled by numpydoc_decorator's @doc. -"messages control".disable = "fixme,missing-class-docstring,missing-function-docstring,too-many-branches,unspecified-encoding,use-dict-literal" +"messages control".disable = "fixme,missing-class-docstring,missing-function-docstring,too-few-public-methods,too-many-branches,unspecified-encoding,use-dict-literal" diff --git a/tests/test_locus_points_smoke.py b/tests/test_locus_points_smoke.py new file mode 100644 index 0000000..8b60a7d --- /dev/null +++ b/tests/test_locus_points_smoke.py @@ -0,0 +1,24 @@ +"""Smoketests for locus-specific points""" + +import pytest + +from looptrace_napari.locus_specific_points_reader import parse_failed + +FAIL_LINES_SAMPLE = """0,13,5.880338307654485,12.20211975317036,10.728294496728491,S +0,17,10.594366532607864,10.95875680073854,20.711938561802768,R;S;xy;z +0,47,10.198132167665957,14.398450914314138,15.378219719077295,R +0,48,0.7173811741523715,13.999908344240598,2.011625563698183,R;z +0,49,6.274792365451074,14.440034853392085,8.81613597404698,S +0,69,10.03907938905064,18.453449673327853,7.594187495036839,R;S;xy;z +0,71,8.157075382512406,12.780500232574296,3.1916456736466685,R +0,74,7.03360935199292,15.324188332927145,3.6995859572616823,xy +0,76,2.426576702500344,20.546530442060508,2.151493771689803,R;S;xy;O +0,77,6.0415531254567885,13.910733825016758,10.238202728231837,S +""" + + +@pytest.mark.parametrize("keep_line_ends", [False, True]) +def test_failed_sample_line_count(keep_line_ends): + lines = FAIL_LINES_SAMPLE.splitlines(keepends=keep_line_ends) + records, _, _ = parse_failed(list(map(lambda l: l.split(","), lines))) + assert 110 == len(records)