Skip to content
Snippets Groups Projects

Compare revisions

Changes are shown as if the source revision was being merged into the target revision. Learn more about comparing revisions.

Source

Select target project
No results found

Target

Select target project
  • swain-lab/aliby/aliby-mirror
  • swain-lab/aliby/alibylite
2 results
Show changes
Commits on Source (55)
Hide whitespace changes
Inline Side-by-side
Showing
with 3280 additions and 2526 deletions
poetry.lock
View file @ b95e79fc
Source diff could not be displayed: it is too large. Options to address this: view the blob.
pyproject.toml
View file @ b95e79fc
...@@ -33,7 +33,7 @@ pathos = "^0.2.8" # Lambda-friendly multithreading ...@@ -33,7 +33,7 @@ pathos = "^0.2.8" # Lambda-friendly multithreading
p-tqdm = "^1.3.3" p-tqdm = "^1.3.3"
pandas = ">=1.3.3" pandas = ">=1.3.3"
py-find-1st = "^1.1.5" # Fast indexing py-find-1st = "^1.1.5" # Fast indexing
scikit-learn = ">=1.0.2" # Used for an extraction metric scikit-learn = ">=1.0.2, <1.3" # Used for an extraction metric
scipy = ">=1.7.3" scipy = ">=1.7.3"
# Pipeline + I/O # Pipeline + I/O
...@@ -46,14 +46,11 @@ xmltodict = "^0.13.0" # read ome-tiff metadata ...@@ -46,14 +46,11 @@ xmltodict = "^0.13.0" # read ome-tiff metadata
zarr = "^2.14.0" zarr = "^2.14.0"
GitPython = "^3.1.27" GitPython = "^3.1.27"
h5py = "2.10" # File I/O h5py = "2.10" # File I/O
aliby-baby = "^0.1.17"
# Networking # Networking
omero-py = { version = ">=5.6.2", optional = true } # contact omero server omero-py = { version = ">=5.6.2", optional = true } # contact omero server
# Baby segmentation
aliby-baby = {version = "^0.1.17", optional=true}
# Postprocessing # Postprocessing
[tool.poetry.group.pp.dependencies] [tool.poetry.group.pp.dependencies]
leidenalg = "^0.8.8" leidenalg = "^0.8.8"
...@@ -113,7 +110,6 @@ grid-strategy = {version = "^0.0.1", optional=true} ...@@ -113,7 +110,6 @@ grid-strategy = {version = "^0.0.1", optional=true}
[tool.poetry.extras] [tool.poetry.extras]
omero = ["omero-py"] omero = ["omero-py"]
baby = ["aliby-baby"]
[tool.black] [tool.black]
line-length = 79 line-length = 79
......
src/agora/abc.py
View file @ b95e79fc
...@@ -17,16 +17,14 @@ atomic = t.Union[int, float, str, bool] ...@@ -17,16 +17,14 @@ atomic = t.Union[int, float, str, bool]
class ParametersABC(ABC): class ParametersABC(ABC):
""" """
Defines parameters as attributes and allows parameters to Define parameters as attributes and allow parameters to
be converted to either a dictionary or to yaml. be converted to either a dictionary or to yaml.
No attribute should be called "parameters"! No attribute should be called "parameters"!
""" """
def __init__(self, **kwargs): def __init__(self, **kwargs):
""" """Define parameters as attributes."""
Defines parameters as attributes
"""
assert ( assert (
"parameters" not in kwargs "parameters" not in kwargs
), "No attribute should be named parameters" ), "No attribute should be named parameters"
...@@ -35,8 +33,9 @@ class ParametersABC(ABC): ...@@ -35,8 +33,9 @@ class ParametersABC(ABC):
def to_dict(self, iterable="null") -> t.Dict: def to_dict(self, iterable="null") -> t.Dict:
""" """
Recursive function to return a nested dictionary of the Return a nested dictionary of the attributes of the class instance.
attributes of the class instance.
Uses recursion.
""" """
if isinstance(iterable, dict): if isinstance(iterable, dict):
if any( if any(
...@@ -62,7 +61,8 @@ class ParametersABC(ABC): ...@@ -62,7 +61,8 @@ class ParametersABC(ABC):
def to_yaml(self, path: Union[Path, str] = None): def to_yaml(self, path: Union[Path, str] = None):
""" """
Returns a yaml stream of the attributes of the class instance. Return a yaml stream of the attributes of the class instance.
If path is provided, the yaml stream is saved there. If path is provided, the yaml stream is saved there.
Parameters Parameters
...@@ -81,9 +81,7 @@ class ParametersABC(ABC): ...@@ -81,9 +81,7 @@ class ParametersABC(ABC):
@classmethod @classmethod
def from_yaml(cls, source: Union[Path, str]): def from_yaml(cls, source: Union[Path, str]):
""" """Return instance from a yaml filename or stdin."""
Returns instance from a yaml filename or stdin
"""
is_buffer = True is_buffer = True
try: try:
if Path(source).exists(): if Path(source).exists():
...@@ -107,7 +105,8 @@ class ParametersABC(ABC): ...@@ -107,7 +105,8 @@ class ParametersABC(ABC):
def update(self, name: str, new_value): def update(self, name: str, new_value):
""" """
Update values recursively Update values recursively.
if name is a dictionary, replace data where existing found or add if not. if name is a dictionary, replace data where existing found or add if not.
It warns against type changes. It warns against type changes.
...@@ -179,7 +178,8 @@ def add_to_collection( ...@@ -179,7 +178,8 @@ def add_to_collection(
class ProcessABC(ABC): class ProcessABC(ABC):
""" """
Base class for processes. Base class for processes.
Defines parameters as attributes and requires run method to be defined.
Define parameters as attributes and requires a run method.
""" """
def __init__(self, parameters): def __init__(self, parameters):
...@@ -243,11 +243,9 @@ class StepABC(ProcessABC): ...@@ -243,11 +243,9 @@ class StepABC(ProcessABC):
@timer @timer
def run_tp(self, tp: int, **kwargs): def run_tp(self, tp: int, **kwargs):
""" """Time and log the timing of a step."""
Time and log the timing of a step.
"""
return self._run_tp(tp, **kwargs) return self._run_tp(tp, **kwargs)
def run(self): def run(self):
# Replace run with run_tp # Replace run with run_tp
raise Warning("Steps use run_tp instead of run") raise Warning("Steps use run_tp instead of run.")
src/agora/io/cells.py
View file @ b95e79fc
...@@ -14,183 +14,170 @@ from utils_find_1st import cmp_equal, find_1st ...@@ -14,183 +14,170 @@ from utils_find_1st import cmp_equal, find_1st
class Cells: class Cells:
""" """
Extracts information from an h5 file. This class accesses: Extract information from an h5 file.
Use output from BABY to find cells detected, get, and fill, edge masks
and retrieve mother-bud relationships.
This class accesses in the h5 file:
'cell_info', which contains 'angles', 'cell_label', 'centres', 'cell_info', which contains 'angles', 'cell_label', 'centres',
'edgemasks', 'ellipse_dims', 'mother_assign', 'mother_assign_dynamic', 'edgemasks', 'ellipse_dims', 'mother_assign', 'mother_assign_dynamic',
'radii', 'timepoint', 'trap'. 'radii', 'timepoint', and 'trap'. All of which except for 'edgemasks'
All of these except for 'edgemasks' are a 1D ndarray. are a 1D ndarray.
'trap_info', which contains 'drifts', 'trap_locations' 'trap_info', which contains 'drifts', and 'trap_locations'.
The "timepoint", "cell_label", and "trap" variables are mutually consistent
1D lists.
Examples are self["timepoint"][self.get_idx(1, 3)] to find the time points
where cell 1 was present in trap 3.
""" """
def __init__(self, filename, path="cell_info"): def __init__(self, filename, path="cell_info"):
"""Initialise from a filename."""
self.filename: t.Optional[t.Union[str, Path]] = filename self.filename: t.Optional[t.Union[str, Path]] = filename
self.cinfo_path: t.Optional[str] = path self.cinfo_path: t.Optional[str] = path
self._edgemasks: t.Optional[str] = None self._edgemasks: t.Optional[str] = None
self._tile_size: t.Optional[int] = None self._tile_size: t.Optional[int] = None
def __getitem__(self, item):
"""
Dynamically fetch data from the h5 file and save as an attribute.
These attributes are accessed like dict keys.
"""
assert item != "edgemasks", "Edgemasks must not be loaded as a whole"
_item = "_" + item
if not hasattr(self, _item):
setattr(self, _item, self.fetch(item))
return getattr(self, _item)
def fetch(self, path):
"""Get data from the h5 file."""
with h5py.File(self.filename, mode="r") as f:
return f[self.cinfo_path][path][()]
@classmethod @classmethod
def from_source(cls, source: t.Union[Path, str]): def from_source(cls, source: t.Union[Path, str]):
"""Ensure initiating file is a Path object."""
return cls(Path(source)) return cls(Path(source))
def _log(self, message: str, level: str = "warn"): def _log(self, message: str, level: str = "warn"):
# Log messages in the corresponding level """Log messages in the corresponding level."""
logger = logging.getLogger("aliby") logger = logging.getLogger("aliby")
getattr(logger, level)(f"{self.__class__.__name__}: {message}") getattr(logger, level)(f"{self.__class__.__name__}: {message}")
@staticmethod @staticmethod
def _asdense(array: np.ndarray): def asdense(array: np.ndarray):
"""Convert sparse array to dense array."""
if not isdense(array): if not isdense(array):
array = array.todense() array = array.todense()
return array return array
@staticmethod @staticmethod
def _astype(array: np.ndarray, kind: str): def astype(array: np.ndarray, kind: str):
# Convert sparse arrays if needed and if kind is 'mask' it fills the outline """Convert sparse arrays if needed; if kind is 'mask' fill the outline."""
array = Cells._asdense(array) array = Cells.asdense(array)
if kind == "mask": if kind == "mask":
array = ndimage.binary_fill_holes(array).astype(bool) array = ndimage.binary_fill_holes(array).astype(bool)
return array return array
def _get_idx(self, cell_id: int, trap_id: int): def get_idx(self, cell_id: int, trap_id: int):
# returns boolean array of time points where both the cell with cell_id and the trap with trap_id exist """Return boolean array giving indices for a cell_id and trap_id."""
return (self["cell_label"] == cell_id) & (self["trap"] == trap_id) return (self["cell_label"] == cell_id) & (self["trap"] == trap_id)
@property @property
def max_labels(self) -> t.List[int]: def max_labels(self) -> t.List[int]:
return [max((0, *self.labels_in_trap(i))) for i in range(self.ntraps)] """Return the maximum cell label per tile."""
return [
max((0, *self.cell_labels_in_trap(i))) for i in range(self.ntraps)
]
@property @property
def max_label(self) -> int: def max_label(self) -> int:
"""Return the maximum cell label over all tiles."""
return sum(self.max_labels) return sum(self.max_labels)
@property @property
def ntraps(self) -> int: def ntraps(self) -> int:
# find the number of traps from the h5 file """Find the number of tiles, or traps."""
with h5py.File(self.filename, mode="r") as f: with h5py.File(self.filename, mode="r") as f:
return len(f["trap_info/trap_locations"][()]) return len(f["trap_info/trap_locations"][()])
@property @property
def tinterval(self): def tinterval(self):
"""Return time interval in seconds."""
with h5py.File(self.filename, mode="r") as f: with h5py.File(self.filename, mode="r") as f:
return f.attrs["time_settings/timeinterval"] return f.attrs["time_settings/timeinterval"]
@property @property
def traps(self) -> t.List[int]: def traps(self) -> t.List[int]:
# returns a list of traps """List unique tile, or trap, IDs."""
return list(set(self["trap"])) return list(set(self["trap"]))
@property @property
def tile_size(self) -> t.Union[int, t.Tuple[int], None]: def tile_size(self) -> t.Union[int, t.Tuple[int], None]:
"""Give the x- and y- sizes of a tile."""
if self._tile_size is None: if self._tile_size is None:
with h5py.File(self.filename, mode="r") as f: with h5py.File(self.filename, mode="r") as f:
# self._tile_size = f["trap_info/tile_size"][0]
self._tile_size = f["cell_info/edgemasks"].shape[1:] self._tile_size = f["cell_info/edgemasks"].shape[1:]
return self._tile_size return self._tile_size
def nonempty_tp_in_trap(self, trap_id: int) -> set: def nonempty_tp_in_trap(self, trap_id: int) -> set:
# given a trap_id returns time points in which cells are available """Given a tile, return time points for which cells are available."""
return set(self["timepoint"][self["trap"] == trap_id]) return set(self["timepoint"][self["trap"] == trap_id])
@property @property
def edgemasks(self) -> t.List[np.ndarray]: def edgemasks(self) -> t.List[np.ndarray]:
# returns the masks per tile """Return a list of masks for every cell at every trap and time point."""
if self._edgemasks is None: if self._edgemasks is None:
edgem_path: str = "edgemasks" edgem_path: str = "edgemasks"
self._edgemasks = self._fetch(edgem_path) self._edgemasks = self.fetch(edgem_path)
return self._edgemasks return self._edgemasks
@property @property
def labels(self) -> t.List[t.List[int]]: def labels(self) -> t.List[t.List[int]]:
""" """Return all cell labels per tile as a set for all tiles."""
Return all cell labels in object return [self.cell_labels_in_trap(trap) for trap in range(self.ntraps)]
We use mother_assign to list traps because it is the only property that appears even
when no cells are found
"""
return [self.labels_in_trap(trap) for trap in range(self.ntraps)]
def max_labels_in_frame(self, frame: int) -> t.List[int]: def max_labels_in_frame(self, final_time_point: int) -> t.List[int]:
# Return the maximum label for each trap in the given frame """Get the maximal cell label for each tile within a frame of time."""
max_labels = [ max_labels = [
self["cell_label"][ self["cell_label"][
(self["timepoint"] <= frame) & (self["trap"] == trap_id) (self["timepoint"] <= final_time_point)
& (self["trap"] == trap_id)
] ]
for trap_id in range(self.ntraps) for trap_id in range(self.ntraps)
] ]
return [max([0, *labels]) for labels in max_labels] return [max([0, *labels]) for labels in max_labels]
def where(self, cell_id: int, trap_id: int): def where(self, cell_id: int, trap_id: int):
""" """Return time points, indices, and edge masks for a cell and trap."""
Parameters idx = self.get_idx(cell_id, trap_id)
----------
cell_id: int
Cell index
trap_id: int
Trap index
Returns
----------
indices int array
boolean mask array
edge_ix int array
"""
indices = self._get_idx(cell_id, trap_id)
edgem_ix = self._edgem_where(cell_id, trap_id)
return ( return (
self["timepoint"][indices], self["timepoint"][idx],
indices, idx,
edgem_ix, self.edgemasks_where(cell_id, trap_id),
) )
def mask(self, cell_id, trap_id): def mask(self, cell_id, trap_id):
""" """Return the times and the filled edge masks for a cell and trap."""
Returns the times and the binary masks of a given cell in a given tile.
Parameters
----------
cell_id : int
The unique ID of the cell.
tile_id : int
The unique ID of the tile.
Returns
-------
Tuple[np.ndarray, np.ndarray]
The times when the binary masks were taken and the binary masks of the given cell in the given tile.
"""
times, outlines = self.outline(cell_id, trap_id) times, outlines = self.outline(cell_id, trap_id)
return times, np.array( return times, np.array(
[ndimage.morphology.binary_fill_holes(o) for o in outlines] [ndimage.morphology.binary_fill_holes(o) for o in outlines]
) )
def at_time( def at_time(
self, timepoint: t.Iterable[int], kind="mask" self, timepoint: int, kind="mask"
) -> t.List[t.List[np.ndarray]]: ) -> t.List[t.List[np.ndarray]]:
""" """Return a dict with traps as keys and cell masks as values for a time point."""
Returns a list of lists of binary masks in a given list of time points. idx = self["timepoint"] == timepoint
traps = self["trap"][idx]
Parameters edgemasks = self.edgemasks_from_idx(idx)
----------
timepoints : Iterable[int]
The list of time points for which to return the binary masks.
kind : str, optional
The type of binary masks to return, by default "mask".
Returns
-------
List[List[np.ndarray]]
A list of lists with binary masks grouped by tile IDs.
"""
ix = self["timepoint"] == timepoint
traps = self["trap"][ix]
edgemasks = self._edgem_from_masking(ix)
masks = [ masks = [
self._astype(edgemask, kind) Cells.astype(edgemask, kind)
for edgemask in edgemasks for edgemask in edgemasks
if edgemask.any() if edgemask.any()
] ]
...@@ -199,22 +186,7 @@ class Cells: ...@@ -199,22 +186,7 @@ class Cells:
def at_times( def at_times(
self, timepoints: t.Iterable[int], kind="mask" self, timepoints: t.Iterable[int], kind="mask"
) -> t.List[t.List[np.ndarray]]: ) -> t.List[t.List[np.ndarray]]:
""" """Return a list of lists of cell masks one for specified time point."""
Returns a list of lists of binary masks for a given list of time points.
Parameters
----------
timepoints : Iterable[int]
The list of time points for which to return the binary masks.
kind : str, optional
The type of binary masks to return, by default "mask".
Returns
-------
List[List[np.ndarray]]
A list of lists with binary masks grouped by tile IDs.
"""
return [ return [
[ [
np.stack(tile_masks) if len(tile_masks) else [] np.stack(tile_masks) if len(tile_masks) else []
...@@ -226,91 +198,77 @@ class Cells: ...@@ -226,91 +198,77 @@ class Cells:
def group_by_traps( def group_by_traps(
self, traps: t.Collection, cell_labels: t.Collection self, traps: t.Collection, cell_labels: t.Collection
) -> t.Dict[int, t.List[int]]: ) -> t.Dict[int, t.List[int]]:
""" """Return a dict with traps as keys and a list of labels as values."""
Returns a dict with traps as keys and list of labels as value.
Note that the total number of traps are calculated from Cells.traps.
"""
iterator = groupby(zip(traps, cell_labels), lambda x: x[0]) iterator = groupby(zip(traps, cell_labels), lambda x: x[0])
d = {key: [x[1] for x in group] for key, group in iterator} d = {key: [x[1] for x in group] for key, group in iterator}
d = {i: d.get(i, []) for i in self.traps} d = {i: d.get(i, []) for i in self.traps}
return d return d
def labels_in_trap(self, trap_id: int) -> t.Set[int]: def cell_labels_in_trap(self, trap_id: int) -> t.Set[int]:
# return set of cell ids for a given trap """Return unique cell labels for a given trap."""
return set((self["cell_label"][self["trap"] == trap_id])) return set((self["cell_label"][self["trap"] == trap_id]))
def labels_at_time(self, timepoint: int) -> t.Dict[int, t.List[int]]: def labels_at_time(self, timepoint: int) -> t.Dict[int, t.List[int]]:
"""Return a dict with traps as keys and cell labels as values for a time point."""
labels = self["cell_label"][self["timepoint"] == timepoint] labels = self["cell_label"][self["timepoint"] == timepoint]
traps = self["trap"][self["timepoint"] == timepoint] traps = self["trap"][self["timepoint"] == timepoint]
return self.group_by_traps(traps, labels) return self.group_by_traps(traps, labels)
def __getitem__(self, item): def edgemasks_from_idx(self, idx):
assert item != "edgemasks", "Edgemasks must not be loaded as a whole" """Get edge masks from the h5 file."""
_item = "_" + item
if not hasattr(self, _item):
setattr(self, _item, self._fetch(item))
return getattr(self, _item)
def _fetch(self, path):
with h5py.File(self.filename, mode="r") as f:
return f[self.cinfo_path][path][()]
def _edgem_from_masking(self, mask):
with h5py.File(self.filename, mode="r") as f: with h5py.File(self.filename, mode="r") as f:
edgem = f[self.cinfo_path + "/edgemasks"][mask, ...] edgem = f[self.cinfo_path + "/edgemasks"][idx, ...]
return edgem return edgem
def _edgem_where(self, cell_id, trap_id): def edgemasks_where(self, cell_id, trap_id):
id_mask = self._get_idx(cell_id, trap_id) """Get the edge masks for a given cell and trap for all time points."""
edgem = self._edgem_from_masking(id_mask) idx = self.get_idx(cell_id, trap_id)
edgemasks = self.edgemasks_from_idx(idx)
return edgem return edgemasks
def outline(self, cell_id: int, trap_id: int): def outline(self, cell_id: int, trap_id: int):
id_mask = self._get_idx(cell_id, trap_id) """Get times and edge masks for a given cell and trap."""
times = self["timepoint"][id_mask] idx = self.get_idx(cell_id, trap_id)
times = self["timepoint"][idx]
return times, self._edgem_from_masking(id_mask) return times, self.edgemasks_from_idx(idx)
@property @property
def ntimepoints(self) -> int: def ntimepoints(self) -> int:
"""Return total number of time points in the experiment."""
return self["timepoint"].max() + 1 return self["timepoint"].max() + 1
@cached_property @cached_property
def _cells_vs_tps(self): def cells_vs_tps(self):
# Binary matrix showing the presence of all cells in all time points """Boolean matrix showing when cells are present for all time points."""
ncells_per_tile = [len(x) for x in self.labels] total_ncells = sum([len(x) for x in self.labels])
cells_vs_tps = np.zeros( cells_vs_tps = np.zeros((total_ncells, self.ntimepoints), dtype=bool)
(sum(ncells_per_tile), self.ntimepoints), dtype=bool
)
cells_vs_tps[ cells_vs_tps[
self._cell_cumsum[self["trap"]] + self["cell_label"] - 1, self.cell_cumlsum[self["trap"]] + self["cell_label"] - 1,
self["timepoint"], self["timepoint"],
] = True ] = True
return cells_vs_tps return cells_vs_tps
@cached_property @cached_property
def _cell_cumsum(self): def cell_cumlsum(self):
# Cumulative sum indicating the number of cells per tile """Find cumulative sum over tiles of the number of cells present."""
ncells_per_tile = [len(x) for x in self.labels] ncells_per_tile = [len(x) for x in self.labels]
cumsum = np.roll(np.cumsum(ncells_per_tile), shift=1) cumsum = np.roll(np.cumsum(ncells_per_tile), shift=1)
cumsum[0] = 0 cumsum[0] = 0
return cumsum return cumsum
def _flat_index_to_tuple_location(self, idx: int) -> t.Tuple[int, int]: def index_to_tile_and_cell(self, idx: int) -> t.Tuple[int, int]:
# Convert a cell index to a tuple """Convert an index to the equivalent pair of tile and cell IDs."""
# Note that it assumes tiles and cell labels are flattened, but tile_id = int(np.where(idx + 1 > self.cell_cumlsum)[0][-1])
# it is agnostic to tps cell_label = idx - self.cell_cumlsum[tile_id] + 1
tile_id = int(np.where(idx + 1 > self._cell_cumsum)[0][-1])
cell_label = idx - self._cell_cumsum[tile_id] + 1
return tile_id, cell_label return tile_id, cell_label
@property @property
def _tiles_vs_cells_vs_tps(self): def tiles_vs_cells_vs_tps(self):
"""
Boolean matrix showing if a cell is present.
The matrix is indexed by trap, cell label, and time point.
"""
ncells_mat = np.zeros( ncells_mat = np.zeros(
(self.ntraps, self["cell_label"].max(), self.ntimepoints), (self.ntraps, self["cell_label"].max(), self.ntimepoints),
dtype=bool, dtype=bool,
...@@ -325,32 +283,37 @@ class Cells: ...@@ -325,32 +283,37 @@ class Cells:
min_consecutive_tps: int = 15, min_consecutive_tps: int = 15,
interval: None or t.Tuple[int, int] = None, interval: None or t.Tuple[int, int] = None,
): ):
"""
Find cells present for all time points in a sliding window of time.
The result can be restricted to a particular interval of time.
"""
window = sliding_window_view( window = sliding_window_view(
self._cells_vs_tps, min_consecutive_tps, axis=1 self.cells_vs_tps, min_consecutive_tps, axis=1
) )
tp_min = window.sum(axis=-1) == min_consecutive_tps tp_min = window.sum(axis=-1) == min_consecutive_tps
# apply a filter to restrict to an interval of time
# Apply an interval filter to focucs on a slice
if interval is not None: if interval is not None:
interval = tuple(np.array(interval)) interval = tuple(np.array(interval))
else: else:
interval = (0, window.shape[1]) interval = (0, window.shape[1])
low_boundary, high_boundary = interval low_boundary, high_boundary = interval
tp_min[:, :low_boundary] = False tp_min[:, :low_boundary] = False
tp_min[:, high_boundary:] = False tp_min[:, high_boundary:] = False
return tp_min return tp_min
@lru_cache(20) @lru_cache(20)
def mothers_in_trap(self, trap_id: int): def mothers_in_trap(self, trap_id: int):
"""Return mothers at a trap."""
return self.mothers[trap_id] return self.mothers[trap_id]
@cached_property @cached_property
def mothers(self): def mothers(self):
""" """
Return nested list with final prediction of mother id for each cell Return a list of mother IDs for each cell in each tile.
Use Baby's "mother_assign_dynamic".
An ID of zero implies that no mother was assigned.
""" """
return self.mother_assign_from_dynamic( return self.mother_assign_from_dynamic(
self["mother_assign_dynamic"], self["mother_assign_dynamic"],
...@@ -362,73 +325,71 @@ class Cells: ...@@ -362,73 +325,71 @@ class Cells:
@cached_property @cached_property
def mothers_daughters(self) -> np.ndarray: def mothers_daughters(self) -> np.ndarray:
""" """
Return a single array with three columns, containing information about Return mother-daughter relationships for all tiles.
the mother-daughter relationships: tile, mothers and daughters.
Returns Returns
------- -------
np.ndarray mothers_daughters: np.ndarray
An array with shape (n, 3) where n is the number of mother-daughter pairs found. An array with shape (n, 3) where n is the number of mother-daughter
The columns contain: pairs found. The first column is the tile_id for the tile where the
- tile: the tile where the mother cell is located. mother cell is located. The second column is the cell index of a
- mothers: the index of the mother cell within the tile. mother cell in the tile. The third column is the index of the
- daughters: the index of the daughter cell within the tile. corresponding daughter cell.
""" """
nested_massign = self.mothers # list of arrays, one per tile, giving mothers of each cell in each tile
mothers = self.mothers
if sum([x for y in nested_massign for x in y]): if sum([x for y in mothers for x in y]):
mothers_daughters = np.array( mothers_daughters = np.array(
[ [
(tid, m, d) (trap_id, mother, bud)
for tid, trapcells in enumerate(nested_massign) for trap_id, trapcells in enumerate(mothers)
for d, m in enumerate(trapcells, 1) for bud, mother in enumerate(trapcells, start=1)
if m if mother
], ],
dtype=np.uint16, dtype=np.uint16,
) )
else: else:
mothers_daughters = np.array([]) mothers_daughters = np.array([])
self._log("No mother-daughters assigned") self._log("No mother-daughters assigned")
return mothers_daughters return mothers_daughters
@staticmethod @staticmethod
def mother_assign_to_mb_matrix(ma: t.List[np.array]): def mother_assign_to_mb_matrix(ma: t.List[np.array]):
""" """
Convert from a list of lists of mother-bud paired assignments to a Convert a list of mother-daughters into a boolean sparse matrix.
sparse matrix with a boolean dtype. The rows correspond to
to daughter buds. The values are boolean and indicate whether a Each row in the matrix correspond to daughter buds.
given cell is a mother cell and a given daughter bud is assigned If an entry is True, a given cell is a mother cell and a given
to the mother cell in the next timepoint. daughter bud is assigned to the mother cell in the next time point.
Parameters: Parameters:
----------- -----------
ma : list of lists of integers ma : list of lists of integers
A list of lists of mother-bud assignments. The i-th sublist contains the A list of lists of mother-bud assignments.
bud assignments for the i-th tile. The integers in each sublist The i-th sublist contains the bud assignments for the i-th tile.
represent the mother label, if it is zero no mother was found. The integers in each sublist represent the mother label, with zero
implying no mother found.
Returns: Returns:
-------- --------
mb_matrix : boolean numpy array of shape (n, m) mb_matrix : boolean numpy array of shape (n, m)
An n x m boolean numpy array where n is the total number of cells (sum An n x m array where n is the total number of cells (sum
of the lengths of all sublists in ma) and m is the maximum number of buds of the lengths of all sublists in ma) and m is the maximum
assigned to any mother cell in ma. The value at (i, j) is True if cell i number of buds assigned to any mother cell in ma.
is a daughter cell and cell j is its mother assigned to i. The value at (i, j) is True if cell i is a daughter cell and
cell j is its assigned mother.
Examples: Examples:
-------- --------
ma = [[0, 0, 1], [0, 1, 0]] >>> ma = [[0, 0, 1], [0, 1, 0]]
Cells(None).mother_assign_to_mb_matrix(ma) >>> Cells(None).mother_assign_to_mb_matrix(ma)
# array([[False, False, False, False, False, False], >>> array([[False, False, False, False, False, False],
# [False, False, False, False, False, False], [False, False, False, False, False, False],
# [ True, False, False, False, False, False], [ True, False, False, False, False, False],
# [False, False, False, False, False, False], [False, False, False, False, False, False],
# [False, False, False, True, False, False], [False, False, False, True, False, False],
# [False, False, False, False, False, False]]) [False, False, False, False, False, False]])
""" """
ncells = sum([len(t) for t in ma]) ncells = sum([len(t) for t in ma])
mb_matrix = np.zeros((ncells, ncells), dtype=bool) mb_matrix = np.zeros((ncells, ncells), dtype=bool)
c = 0 c = 0
...@@ -436,69 +397,78 @@ class Cells: ...@@ -436,69 +397,78 @@ class Cells:
for d, m in enumerate(cells): for d, m in enumerate(cells):
if m: if m:
mb_matrix[c + d, c + m - 1] = True mb_matrix[c + d, c + m - 1] = True
c += len(cells) c += len(cells)
return mb_matrix return mb_matrix
@staticmethod @staticmethod
def mother_assign_from_dynamic( def mother_assign_from_dynamic(
ma: np.ndarray, cell_label: t.List[int], trap: t.List[int], ntraps: int ma: np.ndarray,
cell_label: t.List[int],
trap: t.List[int],
ntraps: int,
) -> t.List[t.List[int]]: ) -> t.List[t.List[int]]:
""" """
Interpolate the associated mothers from the 'mother_assign_dynamic' feature. Find mothers from Baby's 'mother_assign_dynamic' variable.
Parameters Parameters
---------- ----------
ma: np.ndarray ma: np.ndarray
An array with shape (n_t, n_c) containing the 'mother_assign_dynamic' feature. An array with of length number of time points times number of cells
containing the 'mother_assign_dynamic' produced by Baby.
cell_label: List[int] cell_label: List[int]
A list containing the cell labels. A list of cell labels.
trap: List[int] trap: List[int]
A list containing the trap labels. A list of trap labels.
ntraps: int ntraps: int
The total number of traps. The total number of traps.
Returns Returns
------- -------
List[List[int]] List[List[int]]
A list of lists containing the interpolated mother assignment for each cell in each trap. A list giving the mothers for each cell at each trap.
""" """
idlist = list(zip(trap, cell_label)) ids = np.unique(list(zip(trap, cell_label)), axis=0)
cell_gid = np.unique(idlist, axis=0) # find when each cell last appeared at its trap
last_lin_preds = [ last_lin_preds = [
find_1st( find_1st(
((cell_label[::-1] == lbl) & (trap[::-1] == tr)), (
(cell_label[::-1] == cell_label_id)
& (trap[::-1] == trap_id)
),
True, True,
cmp_equal, cmp_equal,
) )
for tr, lbl in cell_gid for trap_id, cell_label_id in ids
] ]
# find the cell's mother using the latest prediction from Baby
mother_assign_sorted = ma[::-1][last_lin_preds] mother_assign_sorted = ma[::-1][last_lin_preds]
# rearrange as a list of mother IDs for each cell in each tile
traps = cell_gid[:, 0] traps = ids[:, 0]
iterator = groupby(zip(traps, mother_assign_sorted), lambda x: x[0]) iterator = groupby(zip(traps, mother_assign_sorted), lambda x: x[0])
d = {key: [x[1] for x in group] for key, group in iterator} d = {trap: [x[1] for x in mothers] for trap, mothers in iterator}
nested_massign = [d.get(i, []) for i in range(ntraps)] mothers = [d.get(i, []) for i in range(ntraps)]
return mothers
return nested_massign ###############################################################################
# Apparently unused below here
###############################################################################
@lru_cache(maxsize=200) @lru_cache(maxsize=200)
def labelled_in_frame( def labelled_in_frame(
self, frame: int, global_id: bool = False self, frame: int, global_id: bool = False
) -> np.ndarray: ) -> np.ndarray:
""" """
Returns labels in a 4D ndarray with the global ids with shape Return labels in a 4D ndarray with potentially global ids.
(ntraps, max_nlabels, ysize, xsize) at a given frame.
Use lru_cache to cache the results for speed.
Parameters Parameters
---------- ----------
frame : int frame : int
The frame number. The frame number (time point).
global_id : bool, optional global_id : bool, optional
If True, the returned array contains global ids, otherwise it If True, the returned array contains global ids, otherwise only
contains only the local ids of the labels. Default is False. the local ids of the labels.
Returns Returns
------- -------
...@@ -507,18 +477,12 @@ class Cells: ...@@ -507,18 +477,12 @@ class Cells:
The array has dimensions (ntraps, max_nlabels, ysize, xsize), The array has dimensions (ntraps, max_nlabels, ysize, xsize),
where max_nlabels is specific for this frame, not the entire where max_nlabels is specific for this frame, not the entire
experiment. experiment.
Notes
-----
This method uses lru_cache to cache the results for faster access.
""" """
labels_in_frame = self.labels_at_time(frame) labels_in_frame = self.labels_at_time(frame)
n_labels = [ n_labels = [
len(labels_in_frame.get(trap_id, [])) len(labels_in_frame.get(trap_id, []))
for trap_id in range(self.ntraps) for trap_id in range(self.ntraps)
] ]
# maxes = self.max_labels_in_frame(frame)
stacks_in_frame = self.get_stacks_in_frame(frame, self.tile_size) stacks_in_frame = self.get_stacks_in_frame(frame, self.tile_size)
first_id = np.cumsum([0, *n_labels]) first_id = np.cumsum([0, *n_labels])
labels_mat = np.zeros( labels_mat = np.zeros(
...@@ -552,7 +516,9 @@ class Cells: ...@@ -552,7 +516,9 @@ class Cells:
self, frame: int, tile_shape: t.Tuple[int] self, frame: int, tile_shape: t.Tuple[int]
) -> t.List[np.ndarray]: ) -> t.List[np.ndarray]:
""" """
Returns a list of stacked masks, each corresponding to a tile at a given timepoint. Return a list of stacked masks.
Each corresponds to a tile at a given time point.
Parameters Parameters
---------- ----------
...@@ -564,7 +530,7 @@ class Cells: ...@@ -564,7 +530,7 @@ class Cells:
Returns Returns
------- -------
List[np.ndarray] List[np.ndarray]
List of stacked masks for each tile at the given timepoint. List of stacked masks for each tile at the given time point.
""" """
masks = self.at_time(frame) masks = self.at_time(frame)
return [ return [
...@@ -574,7 +540,7 @@ class Cells: ...@@ -574,7 +540,7 @@ class Cells:
for trap_id in range(self.ntraps) for trap_id in range(self.ntraps)
] ]
def _sample_tiles_tps( def sample_tiles_tps(
self, self,
size=1, size=1,
min_consecutive_ntps: int = 15, min_consecutive_ntps: int = 15,
...@@ -582,7 +548,7 @@ class Cells: ...@@ -582,7 +548,7 @@ class Cells:
interval=None, interval=None,
) -> t.Tuple[np.ndarray, np.ndarray]: ) -> t.Tuple[np.ndarray, np.ndarray]:
""" """
Sample tiles that have a minimum number of cells and are occupied for at least a minimum number of consecutive timepoints. Sample tiles that have a minimum number of cells and are occupied for at least a minimum number of consecutive time points.
Parameters Parameters
---------- ----------
...@@ -591,7 +557,7 @@ class Cells: ...@@ -591,7 +557,7 @@ class Cells:
min_ncells: int, optional (default=2) min_ncells: int, optional (default=2)
The minimum number of cells per tile. The minimum number of cells per tile.
min_consecutive_ntps: int, optional (default=5) min_consecutive_ntps: int, optional (default=5)
The minimum number of consecutive timepoints a cell must be present in a trap. The minimum number of consecutive timep oints a cell must be present in a trap.
seed: int, optional (default=0) seed: int, optional (default=0)
Random seed value for reproducibility. Random seed value for reproducibility.
interval: None or Tuple(int,int), optional (default=None) interval: None or Tuple(int,int), optional (default=None)
...@@ -612,21 +578,15 @@ class Cells: ...@@ -612,21 +578,15 @@ class Cells:
min_consecutive_tps=min_consecutive_ntps, min_consecutive_tps=min_consecutive_ntps,
interval=interval, interval=interval,
) )
# Find all valid tiles with min_ncells for at least min_tps # Find all valid tiles with min_ncells for at least min_tps
index_id, tps = np.where(cell_availability_matrix) index_id, tps = np.where(cell_availability_matrix)
if interval is None: # Limit search if interval is None: # Limit search
interval = (0, cell_availability_matrix.shape[1]) interval = (0, cell_availability_matrix.shape[1])
np.random.seed(seed) np.random.seed(seed)
choices = np.random.randint(len(index_id), size=size) choices = np.random.randint(len(index_id), size=size)
linear_indices = np.zeros_like(self["cell_label"], dtype=bool) linear_indices = np.zeros_like(self["cell_label"], dtype=bool)
for cell_index_flat, tp in zip(index_id[choices], tps[choices]): for cell_index_flat, tp in zip(index_id[choices], tps[choices]):
tile_id, cell_label = self._flat_index_to_tuple_location( tile_id, cell_label = self.index_to_tile_and_cell(cell_index_flat)
cell_index_flat
)
linear_indices[ linear_indices[
( (
(self["cell_label"] == cell_label) (self["cell_label"] == cell_label)
...@@ -634,10 +594,9 @@ class Cells: ...@@ -634,10 +594,9 @@ class Cells:
& (self["timepoint"] == tp) & (self["timepoint"] == tp)
) )
] = True ] = True
return linear_indices return linear_indices
def _sample_masks( def sample_masks(
self, self,
size: int = 1, size: int = 1,
min_consecutive_ntps: int = 15, min_consecutive_ntps: int = 15,
...@@ -668,31 +627,28 @@ class Cells: ...@@ -668,31 +627,28 @@ class Cells:
The second tuple contains: The second tuple contains:
- `masks`: A list of 2D numpy arrays representing the binary masks of the sampled cells at each timepoint. - `masks`: A list of 2D numpy arrays representing the binary masks of the sampled cells at each timepoint.
""" """
sampled_bitmask = self._sample_tiles_tps( sampled_bitmask = self.sample_tiles_tps(
size=size, size=size,
min_consecutive_ntps=min_consecutive_ntps, min_consecutive_ntps=min_consecutive_ntps,
seed=seed, seed=seed,
interval=interval, interval=interval,
) )
# Sort sampled tiles to use automatic cache when possible # Sort sampled tiles to use automatic cache when possible
tile_ids = self["trap"][sampled_bitmask] tile_ids = self["trap"][sampled_bitmask]
cell_labels = self["cell_label"][sampled_bitmask] cell_labels = self["cell_label"][sampled_bitmask]
tps = self["timepoint"][sampled_bitmask] tps = self["timepoint"][sampled_bitmask]
masks = [] masks = []
for tile_id, cell_label, tp in zip(tile_ids, cell_labels, tps): for tile_id, cell_label, tp in zip(tile_ids, cell_labels, tps):
local_idx = self.labels_at_time(tp)[tile_id].index(cell_label) local_idx = self.labels_at_time(tp)[tile_id].index(cell_label)
tile_mask = self.at_time(tp)[tile_id][local_idx] tile_mask = self.at_time(tp)[tile_id][local_idx]
masks.append(tile_mask) masks.append(tile_mask)
return (tile_ids, cell_labels, tps), np.stack(masks) return (tile_ids, cell_labels, tps), np.stack(masks)
def matrix_trap_tp_where( def matrix_trap_tp_where(
self, min_ncells: int = 2, min_consecutive_tps: int = 5 self, min_ncells: int = 2, min_consecutive_tps: int = 5
): ):
""" """
NOTE CURRENLTY UNUSED WITHIN ALIBY THE MOMENT. MAY BE USEFUL IN THE FUTURE. NOTE CURRENTLY UNUSED BUT USEFUL.
Return a matrix of shape (ntraps x ntps - min_consecutive_tps) to Return a matrix of shape (ntraps x ntps - min_consecutive_tps) to
indicate traps and time-points where min_ncells are available for at least min_consecutive_tps indicate traps and time-points where min_ncells are available for at least min_consecutive_tps
...@@ -708,9 +664,8 @@ class Cells: ...@@ -708,9 +664,8 @@ class Cells:
(ntraps x ( ntps-min_consecutive_tps )) 2D boolean numpy array where rows are trap ids and columns are timepoint windows. (ntraps x ( ntps-min_consecutive_tps )) 2D boolean numpy array where rows are trap ids and columns are timepoint windows.
If the value in a cell is true its corresponding trap and timepoint contains more than min_ncells for at least min_consecutive time-points. If the value in a cell is true its corresponding trap and timepoint contains more than min_ncells for at least min_consecutive time-points.
""" """
window = sliding_window_view( window = sliding_window_view(
self._tiles_vs_cells_vs_tps, min_consecutive_tps, axis=2 self.tiles_vs_cells_vs_tps, min_consecutive_tps, axis=2
) )
tp_min = window.sum(axis=-1) == min_consecutive_tps tp_min = window.sum(axis=-1) == min_consecutive_tps
ncells_tp_min = tp_min.sum(axis=1) >= min_ncells ncells_tp_min = tp_min.sum(axis=1) >= min_ncells
...@@ -720,7 +675,7 @@ class Cells: ...@@ -720,7 +675,7 @@ class Cells:
def stack_masks_in_tile( def stack_masks_in_tile(
masks: t.List[np.ndarray], tile_shape: t.Tuple[int] masks: t.List[np.ndarray], tile_shape: t.Tuple[int]
) -> np.ndarray: ) -> np.ndarray:
# Stack all masks in a trap padding accordingly if no outlines found """Stack all masks in a trap, padding accordingly if no outlines found."""
result = np.zeros((0, *tile_shape), dtype=bool) result = np.zeros((0, *tile_shape), dtype=bool)
if len(masks): if len(masks):
result = np.stack(masks) result = np.stack(masks)
......
src/agora/io/decorators.py
View file @ b95e79fc
...@@ -6,17 +6,19 @@ import typing as t ...@@ -6,17 +6,19 @@ import typing as t
from functools import wraps from functools import wraps
def _first_arg_str_to_df( def _first_arg_str_to_raw_df(
fn: t.Callable, fn: t.Callable,
): ):
"""Enable Signal-like classes to convert strings to data sets.""" """Enable Signal-like classes to convert strings to data sets."""
@wraps(fn) @wraps(fn)
def format_input(*args, **kwargs): def format_input(*args, **kwargs):
cls = args[0] cls = args[0]
data = args[1] data = args[1]
if isinstance(data, str): if isinstance(data, str):
# get data from h5 file # get data from h5 file using Signal's get_raw
data = cls.get_raw(data) data = cls.get_raw(data)
# replace path in the undecorated function with data # replace path in the undecorated function with data
return fn(cls, data, *args[2:], **kwargs) return fn(cls, data, *args[2:], **kwargs)
return format_input return format_input
src/agora/io/metadata.py
View file @ b95e79fc
...@@ -66,7 +66,7 @@ class MetaData: ...@@ -66,7 +66,7 @@ class MetaData:
# Needed because HDF5 attributes do not support dictionaries # Needed because HDF5 attributes do not support dictionaries
def flatten_dict(nested_dict, separator="/"): def flatten_dict(nested_dict, separator="/"):
""" """
Flattens nested dictionary. If empty return as-is. Flatten nested dictionary. If empty return as-is.
""" """
flattened = {} flattened = {}
if nested_dict: if nested_dict:
...@@ -79,9 +79,7 @@ def flatten_dict(nested_dict, separator="/"): ...@@ -79,9 +79,7 @@ def flatten_dict(nested_dict, separator="/"):
# Needed because HDF5 attributes do not support datetime objects # Needed because HDF5 attributes do not support datetime objects
# Takes care of time zones & daylight saving # Takes care of time zones & daylight saving
def datetime_to_timestamp(time, locale="Europe/London"): def datetime_to_timestamp(time, locale="Europe/London"):
""" """Convert datetime object to UNIX timestamp."""
Convert datetime object to UNIX timestamp
"""
return timezone(locale).localize(time).timestamp() return timezone(locale).localize(time).timestamp()
...@@ -189,36 +187,37 @@ def parse_swainlab_metadata(filedir: t.Union[str, Path]): ...@@ -189,36 +187,37 @@ def parse_swainlab_metadata(filedir: t.Union[str, Path]):
Dictionary with minimal metadata Dictionary with minimal metadata
""" """
filedir = Path(filedir) filedir = Path(filedir)
filepath = find_file(filedir, "*.log") filepath = find_file(filedir, "*.log")
if filepath: if filepath:
# new log files
raw_parse = parse_from_swainlab_grammar(filepath) raw_parse = parse_from_swainlab_grammar(filepath)
minimal_meta = get_meta_swainlab(raw_parse) minimal_meta = get_meta_swainlab(raw_parse)
else: else:
# old log files
if filedir.is_file() or str(filedir).endswith(".zarr"): if filedir.is_file() or str(filedir).endswith(".zarr"):
# log file is in parent directory
filedir = filedir.parent filedir = filedir.parent
legacy_parse = parse_logfiles(filedir) legacy_parse = parse_logfiles(filedir)
minimal_meta = ( minimal_meta = (
get_meta_from_legacy(legacy_parse) if legacy_parse else {} get_meta_from_legacy(legacy_parse) if legacy_parse else {}
) )
return minimal_meta return minimal_meta
def dispatch_metadata_parser(filepath: t.Union[str, Path]): def dispatch_metadata_parser(filepath: t.Union[str, Path]):
""" """
Function to dispatch different metadata parsers that convert logfiles into a Dispatch different metadata parsers that convert logfiles into a dictionary.
basic metadata dictionary. Currently only contains the swainlab log parsers.
Currently only contains the swainlab log parsers.
Input: Input:
-------- --------
filepath: str existing file containing metadata, or folder containing naming conventions filepath: str existing file containing metadata, or folder containing naming
conventions
""" """
parsed_meta = parse_swainlab_metadata(filepath) parsed_meta = parse_swainlab_metadata(filepath)
if parsed_meta is None: if parsed_meta is None:
parsed_meta = dir_to_meta parsed_meta = dir_to_meta
return parsed_meta return parsed_meta
......
src/agora/io/signal.py
View file @ b95e79fc
...@@ -10,8 +10,8 @@ import numpy as np ...@@ -10,8 +10,8 @@ import numpy as np
import pandas as pd import pandas as pd
from agora.io.bridge import BridgeH5 from agora.io.bridge import BridgeH5
from agora.io.decorators import _first_arg_str_to_df from agora.io.decorators import _first_arg_str_to_raw_df
from agora.utils.indexing import validate_association from agora.utils.indexing import validate_lineage
from agora.utils.kymograph import add_index_levels from agora.utils.kymograph import add_index_levels
from agora.utils.merge import apply_merges from agora.utils.merge import apply_merges
...@@ -20,11 +20,14 @@ class Signal(BridgeH5): ...@@ -20,11 +20,14 @@ class Signal(BridgeH5):
""" """
Fetch data from h5 files for post-processing. Fetch data from h5 files for post-processing.
Signal assumes that the metadata and data are accessible to perform time-adjustments and apply previously recorded post-processes. Signal assumes that the metadata and data are accessible to
perform time-adjustments and apply previously recorded
post-processes.
""" """
def __init__(self, file: t.Union[str, Path]): def __init__(self, file: t.Union[str, Path]):
"""Define index_names for dataframes, candidate fluorescence channels, and composite statistics.""" """Define index_names for dataframes, candidate fluorescence channels,
and composite statistics."""
super().__init__(file, flag=None) super().__init__(file, flag=None)
self.index_names = ( self.index_names = (
"experiment", "experiment",
...@@ -46,9 +49,9 @@ class Signal(BridgeH5): ...@@ -46,9 +49,9 @@ class Signal(BridgeH5):
def __getitem__(self, dsets: t.Union[str, t.Collection]): def __getitem__(self, dsets: t.Union[str, t.Collection]):
"""Get and potentially pre-process data from h5 file and return as a dataframe.""" """Get and potentially pre-process data from h5 file and return as a dataframe."""
if isinstance(dsets, str): # no pre-processing if isinstance(dsets, str):
return self.get(dsets) return self.get(dsets)
elif isinstance(dsets, list): # pre-processing elif isinstance(dsets, list):
is_bgd = [dset.endswith("imBackground") for dset in dsets] is_bgd = [dset.endswith("imBackground") for dset in dsets]
# Check we are not comparing tile-indexed and cell-indexed data # Check we are not comparing tile-indexed and cell-indexed data
assert sum(is_bgd) == 0 or sum(is_bgd) == len( assert sum(is_bgd) == 0 or sum(is_bgd) == len(
...@@ -58,22 +61,23 @@ class Signal(BridgeH5): ...@@ -58,22 +61,23 @@ class Signal(BridgeH5):
else: else:
raise Exception(f"Invalid type {type(dsets)} to get datasets") raise Exception(f"Invalid type {type(dsets)} to get datasets")
def get(self, dsets: t.Union[str, t.Collection], **kwargs): def get(self, dset_name: t.Union[str, t.Collection], **kwargs):
"""Get and potentially pre-process data from h5 file and return as a dataframe.""" """Return pre-processed data as a dataframe."""
if isinstance(dsets, str): # no pre-processing if isinstance(dset_name, str):
df = self.get_raw(dsets, **kwargs) dsets = self.get_raw(dset_name, **kwargs)
prepost_applied = self.apply_prepost(dsets, **kwargs) prepost_applied = self.apply_prepost(dsets, **kwargs)
return self.add_name(prepost_applied, dset_name)
return self.add_name(prepost_applied, dsets) else:
raise Exception("Error in Signal.get")
@staticmethod @staticmethod
def add_name(df, name): def add_name(df, name):
"""Add column of identical strings to a dataframe.""" """Add name of the Signal as an attribute to its corresponding dataframe."""
df.name = name df.name = name
return df return df
def cols_in_mins(self, df: pd.DataFrame): def cols_in_mins(self, df: pd.DataFrame):
# Convert numerical columns in a dataframe to minutes """Convert numerical columns in a dataframe to minutes."""
try: try:
df.columns = (df.columns * self.tinterval // 60).astype(int) df.columns = (df.columns * self.tinterval // 60).astype(int)
except Exception as e: except Exception as e:
...@@ -94,14 +98,15 @@ class Signal(BridgeH5): ...@@ -94,14 +98,15 @@ class Signal(BridgeH5):
if tinterval_location in f.attrs: if tinterval_location in f.attrs:
return f.attrs[tinterval_location][0] return f.attrs[tinterval_location][0]
else: else:
logging.getlogger("aliby").warn( logging.getLogger("aliby").warn(
f"{str(self.filename).split('/')[-1]}: using default time interval of 5 minutes" f"{str(self.filename).split('/')[-1]}: using default time interval of 5 minutes"
) )
return 5 return 5
@staticmethod @staticmethod
def get_retained(df, cutoff): def get_retained(df, cutoff):
"""Return a fraction of the df, one without later time points.""" """Return rows of df with at least cutoff fraction of the total number
of time points."""
return df.loc[bn.nansum(df.notna(), axis=1) > df.shape[1] * cutoff] return df.loc[bn.nansum(df.notna(), axis=1) > df.shape[1] * cutoff]
@property @property
...@@ -110,15 +115,15 @@ class Signal(BridgeH5): ...@@ -110,15 +115,15 @@ class Signal(BridgeH5):
with h5py.File(self.filename, "r") as f: with h5py.File(self.filename, "r") as f:
return list(f.attrs["channels"]) return list(f.attrs["channels"])
@_first_arg_str_to_df
def retained(self, signal, cutoff=0.8): def retained(self, signal, cutoff=0.8):
""" """
Load data (via decorator) and reduce the resulting dataframe. Load data (via decorator) and reduce the resulting dataframe.
Load data for a signal or a list of signals and reduce the resulting Load data for a signal or a list of signals and reduce the resulting
dataframes to a fraction of their original size, losing late time dataframes to rows with sufficient numbers of time points.
points.
""" """
if isinstance(signal, str):
signal = self.get_raw(signal)
if isinstance(signal, pd.DataFrame): if isinstance(signal, pd.DataFrame):
return self.get_retained(signal, cutoff) return self.get_retained(signal, cutoff)
elif isinstance(signal, list): elif isinstance(signal, list):
...@@ -131,17 +136,15 @@ class Signal(BridgeH5): ...@@ -131,17 +136,15 @@ class Signal(BridgeH5):
""" """
Get lineage data from a given location in the h5 file. Get lineage data from a given location in the h5 file.
Returns an array with three columns: the tile id, the mother label, and the daughter label. Returns an array with three columns: the tile id, the mother label,
and the daughter label.
""" """
if lineage_location is None: if lineage_location is None:
lineage_location = "modifiers/lineage_merged" lineage_location = "modifiers/lineage_merged"
with h5py.File(self.filename, "r") as f: with h5py.File(self.filename, "r") as f:
# if lineage_location not in f:
# lineage_location = lineage_location.split("_")[0]
if lineage_location not in f: if lineage_location not in f:
lineage_location = "postprocessing/lineage" lineage_location = "postprocessing/lineage"
tile_mo_da = f[lineage_location] tile_mo_da = f[lineage_location]
if isinstance(tile_mo_da, h5py.Dataset): if isinstance(tile_mo_da, h5py.Dataset):
lineage = tile_mo_da[()] lineage = tile_mo_da[()]
else: else:
...@@ -154,7 +157,7 @@ class Signal(BridgeH5): ...@@ -154,7 +157,7 @@ class Signal(BridgeH5):
).T ).T
return lineage return lineage
@_first_arg_str_to_df @_first_arg_str_to_raw_df
def apply_prepost( def apply_prepost(
self, self,
data: t.Union[str, pd.DataFrame], data: t.Union[str, pd.DataFrame],
...@@ -162,57 +165,40 @@ class Signal(BridgeH5): ...@@ -162,57 +165,40 @@ class Signal(BridgeH5):
picks: t.Union[t.Collection, bool] = True, picks: t.Union[t.Collection, bool] = True,
): ):
""" """
Apply modifier operations (picker or merger) to a dataframe. Apply picking and merging to a Signal data frame.
Parameters Parameters
---------- ----------
data : t.Union[str, pd.DataFrame] data : t.Union[str, pd.DataFrame]
DataFrame or path to one. A data frame or a path to one.
merges : t.Union[np.ndarray, bool] merges : t.Union[np.ndarray, bool]
(optional) 2-D array with three columns: the tile id, the mother label, and the daughter id. (optional) An array of pairs of (trap, cell) indices to merge.
If True, fetch merges from file. If True, fetch merges from file.
picks : t.Union[np.ndarray, bool] picks : t.Union[np.ndarray, bool]
(optional) 2-D array with two columns: the tiles and (optional) An array of (trap, cell) indices.
the cell labels.
If True, fetch picks from file. If True, fetch picks from file.
Examples
--------
FIXME: Add docs.
""" """
if isinstance(merges, bool): if isinstance(merges, bool):
merges: np.ndarray = self.load_merges() if merges else np.array([]) merges = self.load_merges() if merges else np.array([])
if merges.any(): if merges.any():
merged = apply_merges(data, merges) merged = apply_merges(data, merges)
else: else:
merged = copy(data) merged = copy(data)
if isinstance(picks, bool): if isinstance(picks, bool):
picks = ( picks = (
self.get_picks(names=merged.index.names) self.get_picks(
names=merged.index.names, path="modifiers/picks/"
)
if picks if picks
else set(merged.index) else merged.index
) )
with h5py.File(self.filename, "r") as f: if picks:
if "modifiers/picks" in f and picks: picked_indices = set(picks).intersection(
if picks: [tuple(x) for x in merged.index]
return merged.loc[ )
set(picks).intersection( return merged.loc[picked_indices]
[tuple(x) for x in merged.index] else:
) return merged
]
else:
if isinstance(merged.index, pd.MultiIndex):
empty_lvls = [[] for i in merged.index.names]
index = pd.MultiIndex(
levels=empty_lvls,
codes=empty_lvls,
names=merged.index.names,
)
else:
index = pd.Index([], name=merged.index.name)
merged = pd.DataFrame([], index=index)
return merged
@cached_property @cached_property
def p_available(self): def p_available(self):
...@@ -272,10 +258,11 @@ class Signal(BridgeH5): ...@@ -272,10 +258,11 @@ class Signal(BridgeH5):
Parameters Parameters
---------- ----------
dataset: str or list of strs dataset: str or list of strs
The name of the h5 file or a list of h5 file names The name of the h5 file or a list of h5 file names.
in_minutes: boolean in_minutes: boolean
If True, If True, convert column headings to times in minutes.
lineage: boolean lineage: boolean
If True, add mother_label to index.
""" """
try: try:
if isinstance(dataset, str): if isinstance(dataset, str):
...@@ -288,15 +275,17 @@ class Signal(BridgeH5): ...@@ -288,15 +275,17 @@ class Signal(BridgeH5):
self.get_raw(dset, in_minutes=in_minutes, lineage=lineage) self.get_raw(dset, in_minutes=in_minutes, lineage=lineage)
for dset in dataset for dset in dataset
] ]
if lineage: # assume that df is sorted if lineage:
# assume that df is sorted
mother_label = np.zeros(len(df), dtype=int) mother_label = np.zeros(len(df), dtype=int)
lineage = self.lineage() lineage = self.lineage()
a, b = validate_association( # information on buds
valid_lineage, valid_indices = validate_lineage(
lineage, lineage,
np.array(df.index.to_list()), np.array(df.index.to_list()),
match_column=1, "daughters",
) )
mother_label[b] = lineage[a, 1] mother_label[valid_indices] = lineage[valid_lineage, 1]
df = add_index_levels(df, {"mother_label": mother_label}) df = add_index_levels(df, {"mother_label": mother_label})
return df return df
except Exception as e: except Exception as e:
...@@ -316,13 +305,14 @@ class Signal(BridgeH5): ...@@ -316,13 +305,14 @@ class Signal(BridgeH5):
names: t.Tuple[str, ...] = ("trap", "cell_label"), names: t.Tuple[str, ...] = ("trap", "cell_label"),
path: str = "modifiers/picks/", path: str = "modifiers/picks/",
) -> t.Set[t.Tuple[int, str]]: ) -> t.Set[t.Tuple[int, str]]:
"""Get the relevant picks based on names.""" """Get picks from the h5 file."""
with h5py.File(self.filename, "r") as f: with h5py.File(self.filename, "r") as f:
picks = set()
if path in f: if path in f:
picks = set( picks = set(
zip(*[f[path + name] for name in names if name in f[path]]) zip(*[f[path + name] for name in names if name in f[path]])
) )
else:
picks = set()
return picks return picks
def dataset_to_df(self, f: h5py.File, path: str) -> pd.DataFrame: def dataset_to_df(self, f: h5py.File, path: str) -> pd.DataFrame:
...@@ -353,10 +343,7 @@ class Signal(BridgeH5): ...@@ -353,10 +343,7 @@ class Signal(BridgeH5):
fullname: str, fullname: str,
node: t.Union[h5py.Dataset, h5py.Group], node: t.Union[h5py.Dataset, h5py.Group],
): ):
""" """Store the name of a signal if it is a leaf node and if it starts with extraction."""
Store the name of a signal if it is a leaf node
(a group with no more groups inside) and if it starts with extraction.
"""
if isinstance(node, h5py.Group) and np.all( if isinstance(node, h5py.Group) and np.all(
[isinstance(x, h5py.Dataset) for x in node.values()] [isinstance(x, h5py.Dataset) for x in node.values()]
): ):
......
src/agora/utils/indexing.py
View file @ b95e79fc
...@@ -9,6 +9,152 @@ This can be: ...@@ -9,6 +9,152 @@ This can be:
import numpy as np import numpy as np
import typing as t import typing as t
# data type to link together trap and cell ids
i_dtype = {"names": ["trap_id", "cell_id"], "formats": [np.int64, np.int64]}
def validate_lineage(
lineage: np.ndarray, indices: np.ndarray, how: str = "families"
):
"""
Identify mother-bud pairs that exist both in lineage and a Signal's
indices.
We expect the lineage information to be unique: a bud should not have
two mothers.
Parameters
----------
lineage : np.ndarray
2D array of lineage associations where columns are
(trap, mother, daughter)
or
a 3D array, which is an array of 2 X 2 arrays comprising
[[trap_id, mother_label], [trap_id, daughter_label]].
indices : np.ndarray
A 2D array of cell indices from a Signal, (trap_id, cell_label).
This array should not include mother_label.
how: str
If "mothers", matches indicate mothers from mother-bud pairs;
If "daughters", matches indicate daughters from mother-bud pairs;
If "families", matches indicate mothers and daughters in mother-bud pairs.
Returns
-------
valid_lineage: boolean np.ndarray
1D array indicating matched elements in lineage.
valid_indices: boolean np.ndarray
1D array indicating matched elements in indices.
Examples
--------
>>> import numpy as np
>>> from agora.utils.indexing import validate_lineage
>>> lineage = np.array([ [[0, 1], [0, 3]], [[0, 1], [0, 4]], [[0, 1], [0, 6]], [[0, 4], [0, 7]] ])
>>> indices = np.array([ [0, 1], [0, 2], [0, 3]])
>>> valid_lineage, valid_indices = validate_lineage(lineage, indices)
>>> print(valid_lineage)
array([ True, False, False, False])
>>> print(valid_indices)
array([ True, False, True])
and
>>> lineage = np.array([[[0,3], [0,1]], [[0,2], [0,4]]])
>>> indices = np.array([[0,1], [0,2], [0,3]])
>>> valid_lineage, valid_indices = validate_lineage(lineage, indices)
>>> print(valid_lineage)
array([ True, False])
>>> print(valid_indices)
array([ True, False, True])
"""
if lineage.ndim == 2:
# [trap, mother, daughter] becomes [[trap, mother], [trap, daughter]]
lineage = _assoc_indices_to_3d(lineage)
if how == "mothers":
c_index = 0
elif how == "daughters":
c_index = 1
# find valid lineage
valid_lineages = index_isin(lineage, indices)
if how == "families":
# both mother and bud must be in indices
valid_lineage = valid_lineages.all(axis=1)
else:
valid_lineage = valid_lineages[:, c_index, :]
flat_valid_lineage = valid_lineage.flatten()
# find valid indices
selected_lineages = lineage[flat_valid_lineage, ...]
if how == "families":
# select only pairs of mother and bud indices
valid_indices = index_isin(indices, selected_lineages)
else:
valid_indices = index_isin(indices, selected_lineages[:, c_index, :])
flat_valid_indices = valid_indices.flatten()
if (
indices[flat_valid_indices, :].size
!= np.unique(
lineage[flat_valid_lineage, :].reshape(-1, 2), axis=0
).size
):
# all unique indices in valid_lineages should be in valid_indices
raise Exception(
"Error in validate_lineage: "
"lineage information is likely not unique."
)
return flat_valid_lineage, flat_valid_indices
def index_isin(x: np.ndarray, y: np.ndarray) -> np.ndarray:
"""
Find those elements of x that are in y.
Both arrays must be arrays of integer indices,
such as (trap_id, cell_id).
"""
x = np.ascontiguousarray(x, dtype=np.int64)
y = np.ascontiguousarray(y, dtype=np.int64)
xv = x.view(i_dtype)
inboth = np.intersect1d(xv, y.view(i_dtype))
x_bool = np.isin(xv, inboth)
return x_bool
def _assoc_indices_to_3d(ndarray: np.ndarray):
"""
Convert the last column to a new row and repeat first column's values.
For example: [trap, mother, daughter] becomes
[[trap, mother], [trap, daughter]].
Assumes the input array has shape (N,3).
"""
result = ndarray
if len(ndarray) and ndarray.ndim > 1:
# faster indexing for single positions
if ndarray.shape[1] == 3:
result = np.transpose(
np.hstack((ndarray[:, [0]], ndarray)).reshape(-1, 2, 2),
axes=[0, 2, 1],
)
else:
# 20% slower but more general indexing
columns = np.arange(ndarray.shape[1])
result = np.stack(
(
ndarray[:, np.delete(columns, -1)],
ndarray[:, np.delete(columns, -2)],
),
axis=1,
)
return result
###################################################################
def validate_association( def validate_association(
association: np.ndarray, association: np.ndarray,
...@@ -104,38 +250,8 @@ def validate_association( ...@@ -104,38 +250,8 @@ def validate_association(
return valid_association, valid_indices return valid_association, valid_indices
def _assoc_indices_to_3d(ndarray: np.ndarray):
"""
Convert the last column to a new row while repeating all previous indices.
This is useful when converting a signal multiindex before comparing association.
Assumes the input array has shape (N,3)
"""
result = ndarray
if len(ndarray) and ndarray.ndim > 1:
if ndarray.shape[1] == 3: # Faster indexing for single positions
result = np.transpose(
np.hstack((ndarray[:, [0]], ndarray)).reshape(-1, 2, 2),
axes=[0, 2, 1],
)
else: # 20% slower but more general indexing
columns = np.arange(ndarray.shape[1])
result = np.stack(
(
ndarray[:, np.delete(columns, -1)],
ndarray[:, np.delete(columns, -2)],
),
axis=1,
)
return result
def _3d_index_to_2d(array: np.ndarray): def _3d_index_to_2d(array: np.ndarray):
""" """Revert _assoc_indices_to_3d."""
Opposite to _assoc_indices_to_3d.
"""
result = array result = array
if len(array): if len(array):
result = np.concatenate( result = np.concatenate(
......
src/agora/utils/indexing_new.py 0 → 100644
View file @ b95e79fc
#!/usr/bin/env jupyter
"""
Utilities based on association are used to efficiently acquire indices of
tracklets with some kind of relationship.
This can be:
- Cells that are to be merged.
- Cells that have a lineage relationship.
"""
import numpy as np
import typing as t
def validate_association(
association: np.ndarray,
indices: np.ndarray,
match_column: t.Optional[int] = None,
) -> t.Tuple[np.ndarray, np.ndarray]:
"""
Identify mother-bud pairs that exist both in lineage and a Signal's indices.
"""
if association.ndim == 2:
# reshape into 3D array for broadcasting
# for each trap, [trap, mother, daughter] becomes
# [[trap, mother], [trap, daughter]]
association = _assoc_indices_to_3d(association)
valid_association, valid_indices = validate_lineage(association, indices)
# Alan's working code
# Compare existing association with available indices
# Swap trap and label axes for the association array to correctly cast
valid_ndassociation_a = association[..., None] == indices.T[None, ...]
# Broadcasting is confusing (but efficient):
# First we check the dimension across trap and cell id, to ensure both match
valid_cell_ids_a = valid_ndassociation_a.all(axis=2)
if match_column is None:
# Then we check the merge tuples to check which cases have both target and source
valid_association_a = valid_cell_ids_a.any(axis=2).all(axis=1)
# Finally we check the dimension that crosses all indices, to ensure the pair
# is present in a valid merge event.
valid_indices_a = (
valid_ndassociation_a[valid_association_a]
.all(axis=2)
.any(axis=(0, 1))
)
else: # We fetch specific indices if we aim for the ones with one present
valid_indices_a = valid_cell_ids_a[:, match_column].any(axis=0)
# Valid association then becomes a boolean array, true means that there is a
# match (match_column) between that cell and the index
valid_association_a = (
valid_cell_ids_a[:, match_column] & valid_indices
).any(axis=1)
assert np.array_equal(
valid_association, valid_association_a
), "valid_association error"
assert np.array_equal(
valid_indices, valid_indices_a
), "valid_indices error"
return valid_association, valid_indices
def validate_association_old(
association: np.ndarray,
indices: np.ndarray,
match_column: t.Optional[int] = None,
) -> t.Tuple[np.ndarray, np.ndarray]:
"""
Identify mother-bud pairs that exist both in lineage and a Signal's indices.
Parameters
----------
association : np.ndarray
2D array of lineage associations where columns are (trap, mother, daughter)
or
a 3D array, which is an array of 2 X 2 arrays comprising [[trap_id, mother_label], [trap_id, daughter_label]].
indices : np.ndarray
A 2D array where each column is a different level, such as (trap_id, cell_label), which typically is an index of a Signal
dataframe. This array should not include mother_label.
match_column: int
If 0, matches indicate mothers from mother-bud pairs;
If 1, matches indicate daughters from mother-bud pairs;
If None, matches indicate either mothers or daughters in mother-bud pairs.
Returns
-------
valid_association: boolean np.ndarray
1D array indicating elements in association with matches.
valid_indices: boolean np.ndarray
1D array indicating elements in indices with matches.
Examples
--------
>>> import numpy as np
>>> from agora.utils.indexing import validate_association
>>> association = np.array([ [[0, 1], [0, 3]], [[0, 1], [0, 4]], [[0, 1], [0, 6]], [[0, 4], [0, 7]] ])
>>> indices = np.array([ [0, 1], [0, 2], [0, 3]])
>>> print(indices.T)
>>> valid_association, valid_indices = validate_association(association, indices)
>>> print(valid_association)
array([ True, False, False, False])
>>> print(valid_indices)
array([ True, False, True])
and
>>> association = np.array([[[0,3], [0,1]], [[0,2], [0,4]]])
>>> indices = np.array([[0,1], [0,2], [0,3]])
>>> valid_association, valid_indices = validate_association(association, indices)
>>> print(valid_association)
array([ True, False])
>>> print(valid_indices)
array([ True, False, True])
"""
if association.ndim == 2:
# reshape into 3D array for broadcasting
# for each trap, [trap, mother, daughter] becomes
# [[trap, mother], [trap, daughter]]
association = _assoc_indices_to_3d(association)
# use broadcasting to compare association with indices
# swap trap and cell_label axes for correct broadcasting
indicesT = indices.T
# compare each of [[trap, mother], [trap, daughter]] for all traps
# in association with [trap, cell_label] for all traps in indices
# association is no_traps x 2 x 2; indices is no_traps X 2
# valid_ndassociation is no_traps_association x 2 x 2 x no_traps_indices
valid_ndassociation = (
association[..., np.newaxis] == indicesT[np.newaxis, ...]
)
# find matches in association
###
# make True comparisons with both trap_ids and cell labels matching
# compare trap_ids and cell_ids for each pair of traps
valid_cell_ids = valid_ndassociation.all(axis=2)
if match_column is None:
# make True comparisons match at least one row in indices
# at least one cell_id matches
va_intermediate = valid_cell_ids.any(axis=2)
# make True comparisons have both mother and bud matching rows in indices
valid_association = va_intermediate.all(axis=1)
else:
# match_column selects mothers if 0 and daughters if 1
# make True match at least one row in indices
valid_association = valid_cell_ids[:, match_column].any(axis=1)
# find matches in indices
###
# make True comparisons have a validated association for both the mother and bud
# make True comparisons have both trap_ids and cell labels matching
valid_cell_ids_va = valid_ndassociation[valid_association].all(axis=2)
if match_column is None:
# make True comparisons match either a mother or a bud in association
valid_indices = valid_cell_ids_va.any(axis=(0, 1))
else:
valid_indices = valid_cell_ids_va[:, match_column][0]
# Alan's working code
# Compare existing association with available indices
# Swap trap and label axes for the association array to correctly cast
valid_ndassociation_a = association[..., None] == indices.T[None, ...]
# Broadcasting is confusing (but efficient):
# First we check the dimension across trap and cell id, to ensure both match
valid_cell_ids_a = valid_ndassociation_a.all(axis=2)
if match_column is None:
# Then we check the merge tuples to check which cases have both target and source
valid_association_a = valid_cell_ids_a.any(axis=2).all(axis=1)
# Finally we check the dimension that crosses all indices, to ensure the pair
# is present in a valid merge event.
valid_indices_a = (
valid_ndassociation_a[valid_association_a]
.all(axis=2)
.any(axis=(0, 1))
)
else: # We fetch specific indices if we aim for the ones with one present
valid_indices_a = valid_cell_ids_a[:, match_column].any(axis=0)
# Valid association then becomes a boolean array, true means that there is a
# match (match_column) between that cell and the index
valid_association_a = (
valid_cell_ids_a[:, match_column] & valid_indices
).any(axis=1)
assert np.array_equal(
valid_association, valid_association_a
), "valid_association error"
assert np.array_equal(
valid_indices, valid_indices_a
), "valid_indices error"
return valid_association, valid_indices
def _assoc_indices_to_3d(ndarray: np.ndarray):
"""
Reorganise an array of shape (N, 3) into one of shape (N, 2, 2).
Reorganise an array so that the last entry of each row is removed
and generates a new row. This new row retains all other entries of
the original row.
Example:
[ [0, 1, 3], [0, 1, 4] ]
becomes
[ [[0, 1], [0, 3]], [[0, 1], [0, 4]] ]
"""
result = ndarray
if len(ndarray) and ndarray.ndim > 1:
if ndarray.shape[1] == 3:
# faster indexing for single positions
result = np.transpose(
np.hstack((ndarray[:, [0]], ndarray)).reshape(-1, 2, 2),
axes=[0, 2, 1],
)
else:
# 20% slower, but more general indexing
columns = np.arange(ndarray.shape[1])
result = np.stack(
(
ndarray[:, np.delete(columns, -1)],
ndarray[:, np.delete(columns, -2)],
),
axis=1,
)
return result
def _3d_index_to_2d(array: np.ndarray):
"""Revert switch from _assoc_indices_to_3d."""
result = array
if len(array):
result = np.concatenate(
(array[:, 0, :], array[:, 1, 1, np.newaxis]), axis=1
)
return result
def compare_indices(x: np.ndarray, y: np.ndarray) -> np.ndarray:
"""
Compare two 2D arrays using broadcasting.
Return a binary array where a True value links two cells where
all cells are the same.
"""
return (x[..., np.newaxis] == y.T[np.newaxis, ...]).all(axis=1)
src/agora/utils/kymograph.py
View file @ b95e79fc
...@@ -86,16 +86,19 @@ def bidirectional_retainment_filter( ...@@ -86,16 +86,19 @@ def bidirectional_retainment_filter(
daughters_thresh: int = 7, daughters_thresh: int = 7,
) -> pd.DataFrame: ) -> pd.DataFrame:
""" """
Retrieve families where mothers are present for more than a fraction of the experiment, and daughters for longer than some number of time-points. Retrieve families where mothers are present for more than a fraction
of the experiment and daughters for longer than some number of
time-points.
Parameters Parameters
---------- ----------
df: pd.DataFrame df: pd.DataFrame
Data Data
mothers_thresh: float mothers_thresh: float
Minimum fraction of experiment's total duration for which mothers must be present. Minimum fraction of experiment's total duration for which mothers
must be present.
daughters_thresh: int daughters_thresh: int
Minimum number of time points for which daughters must be observed Minimum number of time points for which daughters must be observed.
""" """
# daughters # daughters
all_daughters = df.loc[df.index.get_level_values("mother_label") > 0] all_daughters = df.loc[df.index.get_level_values("mother_label") > 0]
...@@ -170,6 +173,7 @@ def slices_from_spans(spans: t.Tuple[int], df: pd.DataFrame) -> t.List[slice]: ...@@ -170,6 +173,7 @@ def slices_from_spans(spans: t.Tuple[int], df: pd.DataFrame) -> t.List[slice]:
def drop_mother_label(index: pd.MultiIndex) -> np.ndarray: def drop_mother_label(index: pd.MultiIndex) -> np.ndarray:
"""Remove mother_label level from a MultiIndex."""
no_mother_label = index no_mother_label = index
if "mother_label" in index.names: if "mother_label" in index.names:
no_mother_label = index.droplevel("mother_label") no_mother_label = index.droplevel("mother_label")
......
src/agora/utils/lineage.py
View file @ b95e79fc
#!/usr/bin/env python3 #!/usr/bin/env python3
import re
import typing as t
import numpy as np import numpy as np
import pandas as pd
from agora.io.bridge import groupsort from agora.io.bridge import groupsort
from itertools import groupby
def mb_array_to_dict(mb_array: np.ndarray): def mb_array_to_dict(mb_array: np.ndarray):
...@@ -19,4 +15,3 @@ def mb_array_to_dict(mb_array: np.ndarray): ...@@ -19,4 +15,3 @@ def mb_array_to_dict(mb_array: np.ndarray):
for trap, mo_da in groupsort(mb_array).items() for trap, mo_da in groupsort(mb_array).items()
for mo, daughters in groupsort(mo_da).items() for mo, daughters in groupsort(mo_da).items()
} }
src/agora/utils/merge.py
View file @ b95e79fc
...@@ -3,90 +3,161 @@ ...@@ -3,90 +3,161 @@
Functions to efficiently merge rows in DataFrames. Functions to efficiently merge rows in DataFrames.
""" """
import typing as t import typing as t
from copy import copy
import numpy as np import numpy as np
import pandas as pd import pandas as pd
from utils_find_1st import cmp_larger, find_1st from utils_find_1st import cmp_larger, find_1st
from agora.utils.indexing import compare_indices, validate_association from agora.utils.indexing import index_isin
def group_merges(merges: np.ndarray) -> t.List[t.Tuple]:
"""
Convert merges into a list of merges for traps requiring multiple
merges and then for traps requiring single merges.
"""
left_tracks = merges[:, 0]
right_tracks = merges[:, 1]
# find traps requiring multiple merges
linr = merges[index_isin(left_tracks, right_tracks).flatten(), :]
rinl = merges[index_isin(right_tracks, left_tracks).flatten(), :]
# make unique and order merges for each trap
multi_merge = np.unique(np.concatenate((linr, rinl)), axis=0)
# find traps requiring a singe merge
single_merge = merges[
~index_isin(merges, multi_merge).all(axis=1).flatten(), :
]
# convert to lists of arrays
single_merge_list = [[sm] for sm in single_merge]
multi_merge_list = [
multi_merge[multi_merge[:, 0, 0] == trap_id, ...]
for trap_id in np.unique(multi_merge[:, 0, 0])
]
res = [*multi_merge_list, *single_merge_list]
return res
def merge_lineage(
lineage: np.ndarray, merges: np.ndarray
) -> (np.ndarray, np.ndarray):
"""
Use merges to update lineage information.
Check if merging causes any buds to have multiple mothers and discard
those incorrect merges.
Return updated lineage and merge arrays.
"""
flat_lineage = lineage.reshape(-1, 2)
bud_mother_dict = {
tuple(bud): mother for bud, mother in zip(lineage[:, 1], lineage[:, 0])
}
left_tracks = merges[:, 0]
# find left tracks that are in lineages
valid_lineages = index_isin(flat_lineage, left_tracks).flatten()
# group into multi- and then single merges
grouped_merges = group_merges(merges)
# perform merges
if valid_lineages.any():
# indices of each left track -> indices of rightmost right track
replacement_dict = {
tuple(contig_pair[0]): merge[-1][1]
for merge in grouped_merges
for contig_pair in merge
}
# if both key and value are buds, they must have the same mother
buds = lineage[:, 1]
incorrect_merges = [
key
for key in replacement_dict
if np.any(index_isin(buds, replacement_dict[key]).flatten())
and np.any(index_isin(buds, key).flatten())
and not np.array_equal(
bud_mother_dict[key],
bud_mother_dict[tuple(replacement_dict[key])],
)
]
if incorrect_merges:
# reassign incorrect merges so that they have no affect
for key in incorrect_merges:
replacement_dict[key] = key
# find only correct merges
new_merges = merges[
~index_isin(
merges[:, 0], np.array(incorrect_merges)
).flatten(),
...,
]
else:
new_merges = merges
# correct lineage information
# replace mother or bud index with index of rightmost track
flat_lineage[valid_lineages] = [
replacement_dict[tuple(index)]
for index in flat_lineage[valid_lineages]
]
else:
new_merges = merges
# reverse flattening
new_lineage = flat_lineage.reshape(-1, 2, 2)
# remove any duplicates
new_lineage = np.unique(new_lineage, axis=0)
return new_lineage, new_merges
def apply_merges(data: pd.DataFrame, merges: np.ndarray): def apply_merges(data: pd.DataFrame, merges: np.ndarray):
"""Split data in two, one subset for rows relevant for merging and one """
without them. It uses an array of source tracklets and target tracklets Generate a new data frame containing merged tracks.
to efficiently merge them.
Parameters Parameters
---------- ----------
data : pd.DataFrame data : pd.DataFrame
Input DataFrame. A Signal data frame.
merges : np.ndarray merges : np.ndarray
3-D ndarray where dimensions are (X,2,2): nmerges, source-target An array of pairs of (trap, cell) indices to merge.
pair and single-cell identifiers, respectively.
Examples
--------
FIXME: Add docs.
""" """
indices = data.index indices = data.index
if "mother_label" in indices.names: if "mother_label" in indices.names:
indices = indices.droplevel("mother_label") indices = indices.droplevel("mother_label")
valid_merges, indices = validate_association( indices = np.array(list(indices))
merges, np.array(list(indices)) # merges in the data frame's indices
) valid_merges = index_isin(merges, indices).all(axis=1).flatten()
# corresponding indices for the data frame in merges
# Assign non-merged selected_merges = merges[valid_merges, ...]
merged = data.loc[~indices] valid_indices = index_isin(indices, selected_merges).flatten()
# data not requiring merging
# Implement the merges and drop source rows. merged = data.loc[~valid_indices]
# TODO Use matrices to perform merges in batch # merge tracks
# for ecficiency
if valid_merges.any(): if valid_merges.any():
to_merge = data.loc[indices] to_merge = data.loc[valid_indices].copy()
targets, sources = zip(*merges[valid_merges]) left_indices = merges[:, 0]
for source, target in zip(sources, targets): right_indices = merges[:, 1]
target = tuple(target) # join left track with right track
to_merge.loc[target] = join_tracks_pair( for left_index, right_index in zip(left_indices, right_indices):
to_merge.loc[target].values, to_merge.loc[tuple(left_index)] = join_two_tracks(
to_merge.loc[tuple(source)].values, to_merge.loc[tuple(left_index)].values,
to_merge.loc[tuple(right_index)].values,
) )
to_merge.drop(map(tuple, sources), inplace=True) # drop indices for right tracks
to_merge.drop(map(tuple, right_indices), inplace=True)
# add to data not requiring merges
merged = pd.concat((merged, to_merge), names=data.index.names) merged = pd.concat((merged, to_merge), names=data.index.names)
return merged return merged
def join_tracks_pair(target: np.ndarray, source: np.ndarray) -> np.ndarray: def join_two_tracks(
""" left_track: np.ndarray, right_track: np.ndarray
Join two tracks and return the new value of the target. ) -> np.ndarray:
""" """Join two tracks and return the new one."""
target_copy = target new_track = left_track.copy()
end = find_1st(target_copy[::-1], 0, cmp_larger) # find last positive element by inverting track
target_copy[-end:] = source[-end:] end = find_1st(left_track[::-1], 0, cmp_larger)
return target_copy # merge tracks into one
new_track[-end:] = right_track[-end:]
return new_track
def group_merges(merges: np.ndarray) -> t.List[t.Tuple]:
# Return a list where the cell is present as source and target
# (multimerges)
sources_targets = compare_indices(merges[:, 0, :], merges[:, 1, :])
is_multimerge = sources_targets.any(axis=0) | sources_targets.any(axis=1)
is_monomerge = ~is_multimerge
multimerge_subsets = union_find(zip(*np.where(sources_targets)))
merge_groups = [merges[np.array(tuple(x))] for x in multimerge_subsets]
sorted_merges = list(map(sort_association, merge_groups))
# Ensure that source and target are at the edges ##################################################################
return [
*sorted_merges,
*[[event] for event in merges[is_monomerge]],
]
def union_find(lsts): def union_find(lsts):
...@@ -120,27 +191,3 @@ def sort_association(array: np.ndarray): ...@@ -120,27 +191,3 @@ def sort_association(array: np.ndarray):
[res.append(x) for x in np.flip(order).flatten() if x not in res] [res.append(x) for x in np.flip(order).flatten() if x not in res]
sorted_array = array[np.array(res)] sorted_array = array[np.array(res)]
return sorted_array return sorted_array
def merge_association(
association: np.ndarray, merges: np.ndarray
) -> np.ndarray:
grouped_merges = group_merges(merges)
flat_indices = association.reshape(-1, 2)
comparison_mat = compare_indices(merges[:, 0], flat_indices)
valid_indices = comparison_mat.any(axis=0)
if valid_indices.any(): # Where valid, perform transformation
replacement_d = {}
for dataset in grouped_merges:
for k in dataset:
replacement_d[tuple(k[0])] = dataset[-1][1]
flat_indices[valid_indices] = [
replacement_d[tuple(i)] for i in flat_indices[valid_indices]
]
merged_indices = flat_indices.reshape(-1, 2, 2)
return merged_indices
src/aliby/io/dataset.py
View file @ b95e79fc
...@@ -54,7 +54,7 @@ class DatasetLocalABC(ABC): ...@@ -54,7 +54,7 @@ class DatasetLocalABC(ABC):
Abstract Base class to find local files, either OME-XML or raw images. Abstract Base class to find local files, either OME-XML or raw images.
""" """
_valid_suffixes = ("tiff", "png", "zarr") _valid_suffixes = ("tiff", "png", "zarr", "tif")
_valid_meta_suffixes = ("txt", "log") _valid_meta_suffixes = ("txt", "log")
def __init__(self, dpath: t.Union[str, Path], *args, **kwargs): def __init__(self, dpath: t.Union[str, Path], *args, **kwargs):
......
src/aliby/io/image.py
View file @ b95e79fc
...@@ -30,14 +30,14 @@ from agora.io.metadata import dir_to_meta, dispatch_metadata_parser ...@@ -30,14 +30,14 @@ from agora.io.metadata import dir_to_meta, dispatch_metadata_parser
def get_examples_dir(): def get_examples_dir():
"""Get examples directory which stores dummy image for tiler""" """Get examples directory that stores dummy image for tiler."""
return files("aliby").parent.parent / "examples" / "tiler" return files("aliby").parent.parent / "examples" / "tiler"
def instantiate_image( def instantiate_image(
source: t.Union[str, int, t.Dict[str, str], Path], **kwargs source: t.Union[str, int, t.Dict[str, str], Path], **kwargs
): ):
"""Wrapper to instatiate the appropiate image """Wrapper to instantiate the appropriate image
Parameters Parameters
---------- ----------
...@@ -55,26 +55,26 @@ def instantiate_image( ...@@ -55,26 +55,26 @@ def instantiate_image(
def dispatch_image(source: t.Union[str, int, t.Dict[str, str], Path]): def dispatch_image(source: t.Union[str, int, t.Dict[str, str], Path]):
""" """Pick the appropriate Image class depending on the source of data."""
Wrapper to pick the appropiate Image class depending on the source of data.
"""
if isinstance(source, (int, np.int64)): if isinstance(source, (int, np.int64)):
from aliby.io.omero import Image from aliby.io.omero import Image
instatiator = Image instantiator = Image
elif isinstance(source, dict) or ( elif isinstance(source, dict) or (
isinstance(source, (str, Path)) and Path(source).is_dir() isinstance(source, (str, Path)) and Path(source).is_dir()
): ):
if Path(source).suffix == ".zarr": if Path(source).suffix == ".zarr":
instatiator = ImageZarr instantiator = ImageZarr
else: else:
instatiator = ImageDir instantiator = ImageDir
elif isinstance(source, Path) and source.is_file():
# my addition
instantiator = ImageLocalOME
elif isinstance(source, str) and Path(source).is_file(): elif isinstance(source, str) and Path(source).is_file():
instatiator = ImageLocalOME instantiator = ImageLocalOME
else: else:
raise Exception(f"Invalid data source at {source}") raise Exception(f"Invalid data source at {source}")
return instantiator
return instatiator
class BaseLocalImage(ABC): class BaseLocalImage(ABC):
...@@ -82,6 +82,7 @@ class BaseLocalImage(ABC): ...@@ -82,6 +82,7 @@ class BaseLocalImage(ABC):
Base Image class to set path and provide context management method. Base Image class to set path and provide context management method.
""" """
# default image order
_default_dimorder = "tczyx" _default_dimorder = "tczyx"
def __init__(self, path: t.Union[str, Path]): def __init__(self, path: t.Union[str, Path]):
...@@ -98,8 +99,7 @@ class BaseLocalImage(ABC): ...@@ -98,8 +99,7 @@ class BaseLocalImage(ABC):
return False return False
def rechunk_data(self, img): def rechunk_data(self, img):
# Format image using x and y size from metadata. """Format image using x and y size from metadata."""
self._rechunked_img = da.rechunk( self._rechunked_img = da.rechunk(
img, img,
chunks=( chunks=(
...@@ -145,16 +145,16 @@ class ImageLocalOME(BaseLocalImage): ...@@ -145,16 +145,16 @@ class ImageLocalOME(BaseLocalImage):
in which a multidimensional tiff image contains the metadata. in which a multidimensional tiff image contains the metadata.
""" """
def __init__(self, path: str, dimorder=None): def __init__(self, path: str, dimorder=None, **kwargs):
super().__init__(path) super().__init__(path)
self._id = str(path) self._id = str(path)
self.set_meta(str(path))
def set_meta(self): def set_meta(self, path):
meta = dict() meta = dict()
try: try:
with TiffFile(path) as f: with TiffFile(path) as f:
self._meta = xmltodict.parse(f.ome_metadata)["OME"] self._meta = xmltodict.parse(f.ome_metadata)["OME"]
for dim in self.dimorder: for dim in self.dimorder:
meta["size_" + dim.lower()] = int( meta["size_" + dim.lower()] = int(
self._meta["Image"]["Pixels"]["@Size" + dim] self._meta["Image"]["Pixels"]["@Size" + dim]
...@@ -165,21 +165,19 @@ class ImageLocalOME(BaseLocalImage): ...@@ -165,21 +165,19 @@ class ImageLocalOME(BaseLocalImage):
] ]
meta["name"] = self._meta["Image"]["@Name"] meta["name"] = self._meta["Image"]["@Name"]
meta["type"] = self._meta["Image"]["Pixels"]["@Type"] meta["type"] = self._meta["Image"]["Pixels"]["@Type"]
except Exception as e:
except Exception as e: # Images not in OMEXML # images not in OMEXML
print("Warning:Metadata not found: {}".format(e)) print("Warning:Metadata not found: {}".format(e))
print( print(
f"Warning: No dimensional info provided. Assuming {self._default_dimorder}" "Warning: No dimensional info provided. "
f"Assuming {self._default_dimorder}"
) )
# mark non-existent dimensions for padding
# Mark non-existent dimensions for padding
self.base = self._default_dimorder self.base = self._default_dimorder
# self.ids = [self.index(i) for i in dimorder] # self.ids = [self.index(i) for i in dimorder]
self._dimorder = self.base
self._dimorder = base
self._meta = meta self._meta = meta
# self._meta["name"] = Path(path).name.split(".")[0]
@property @property
def name(self): def name(self):
...@@ -246,7 +244,7 @@ class ImageDir(BaseLocalImage): ...@@ -246,7 +244,7 @@ class ImageDir(BaseLocalImage):
It inherits from BaseLocalImage so we only override methods that are critical. It inherits from BaseLocalImage so we only override methods that are critical.
Assumptions: Assumptions:
- One folders per position. - One folder per position.
- Images are flat. - Images are flat.
- Channel, Time, z-stack and the others are determined by filenames. - Channel, Time, z-stack and the others are determined by filenames.
- Provides Dimorder as it is set in the filenames, or expects order during instatiation - Provides Dimorder as it is set in the filenames, or expects order during instatiation
...@@ -318,7 +316,7 @@ class ImageZarr(BaseLocalImage): ...@@ -318,7 +316,7 @@ class ImageZarr(BaseLocalImage):
print(f"Could not add size info to metadata: {e}") print(f"Could not add size info to metadata: {e}")
def get_data_lazy(self) -> da.Array: def get_data_lazy(self) -> da.Array:
"""Return 5D dask array. For lazy-loading local multidimensional zarr files""" """Return 5D dask array for lazy-loading local multidimensional zarr files."""
return self._img return self._img
def add_size_to_meta(self): def add_size_to_meta(self):
......
src/aliby/pipeline.py
View file @ b95e79fc
...@@ -7,10 +7,10 @@ import typing as t ...@@ -7,10 +7,10 @@ import typing as t
from copy import copy from copy import copy
from importlib.metadata import version from importlib.metadata import version
from pathlib import Path from pathlib import Path
from pprint import pprint
import h5py import h5py
import numpy as np import numpy as np
import pandas as pd
from pathos.multiprocessing import Pool from pathos.multiprocessing import Pool
from tqdm import tqdm from tqdm import tqdm
...@@ -18,11 +18,7 @@ from agora.abc import ParametersABC, ProcessABC ...@@ -18,11 +18,7 @@ from agora.abc import ParametersABC, ProcessABC
from agora.io.metadata import MetaData, parse_logfiles from agora.io.metadata import MetaData, parse_logfiles
from agora.io.reader import StateReader from agora.io.reader import StateReader
from agora.io.signal import Signal from agora.io.signal import Signal
from agora.io.writer import ( from agora.io.writer import LinearBabyWriter, StateWriter, TilerWriter
LinearBabyWriter,
StateWriter,
TilerWriter,
)
from aliby.baby_client import BabyParameters, BabyRunner from aliby.baby_client import BabyParameters, BabyRunner
from aliby.haystack import initialise_tf from aliby.haystack import initialise_tf
from aliby.io.dataset import dispatch_dataset from aliby.io.dataset import dispatch_dataset
...@@ -32,6 +28,10 @@ from extraction.core.extractor import Extractor, ExtractorParameters ...@@ -32,6 +28,10 @@ from extraction.core.extractor import Extractor, ExtractorParameters
from extraction.core.functions.defaults import exparams_from_meta from extraction.core.functions.defaults import exparams_from_meta
from postprocessor.core.processor import PostProcessor, PostProcessorParameters from postprocessor.core.processor import PostProcessor, PostProcessorParameters
# stop warnings from TensorFlow
os.environ["TF_CPP_MIN_LOG_LEVEL"] = "2"
logging.getLogger("tensorflow").setLevel(logging.ERROR)
class PipelineParameters(ParametersABC): class PipelineParameters(ParametersABC):
"""Define parameters for the steps of the pipeline.""" """Define parameters for the steps of the pipeline."""
...@@ -39,7 +39,12 @@ class PipelineParameters(ParametersABC): ...@@ -39,7 +39,12 @@ class PipelineParameters(ParametersABC):
_pool_index = None _pool_index = None
def __init__( def __init__(
self, general, tiler, baby, extraction, postprocessing, reporting self,
general,
tiler,
baby,
extraction,
postprocessing,
): ):
"""Initialise, but called by a class method - not directly.""" """Initialise, but called by a class method - not directly."""
self.general = general self.general = general
...@@ -47,7 +52,6 @@ class PipelineParameters(ParametersABC): ...@@ -47,7 +52,6 @@ class PipelineParameters(ParametersABC):
self.baby = baby self.baby = baby
self.extraction = extraction self.extraction = extraction
self.postprocessing = postprocessing self.postprocessing = postprocessing
self.reporting = reporting
@classmethod @classmethod
def default( def default(
...@@ -76,16 +80,15 @@ class PipelineParameters(ParametersABC): ...@@ -76,16 +80,15 @@ class PipelineParameters(ParametersABC):
postprocessing: dict (optional) postprocessing: dict (optional)
Parameters for post-processing. Parameters for post-processing.
""" """
expt_id = general.get("expt_id", 19993) if (
if isinstance(expt_id, Path): isinstance(general["expt_id"], Path)
assert expt_id.exists() and general["expt_id"].exists()
):
expt_id = str(expt_id) expt_id = str(general["expt_id"])
general["expt_id"] = expt_id else:
expt_id = general["expt_id"]
directory = Path(general["directory"]) directory = Path(general["directory"])
# get metadata from log files either locally or via OMERO
# get log files, either locally or via OMERO
with dispatch_dataset( with dispatch_dataset(
expt_id, expt_id,
**{k: general.get(k) for k in ("host", "username", "password")}, **{k: general.get(k) for k in ("host", "username", "password")},
...@@ -107,7 +110,6 @@ class PipelineParameters(ParametersABC): ...@@ -107,7 +110,6 @@ class PipelineParameters(ParametersABC):
} }
# set minimal metadata # set minimal metadata
meta_d = minimal_default_meta meta_d = minimal_default_meta
# define default values for general parameters # define default values for general parameters
tps = meta_d.get("ntps", 2000) tps = meta_d.get("ntps", 2000)
defaults = { defaults = {
...@@ -128,8 +130,7 @@ class PipelineParameters(ParametersABC): ...@@ -128,8 +130,7 @@ class PipelineParameters(ParametersABC):
use_explog=True, use_explog=True,
) )
} }
# update default values for general using inputs
# update default values using inputs
for k, v in general.items(): for k, v in general.items():
if k not in defaults["general"]: if k not in defaults["general"]:
defaults["general"][k] = v defaults["general"][k] = v
...@@ -138,11 +139,9 @@ class PipelineParameters(ParametersABC): ...@@ -138,11 +139,9 @@ class PipelineParameters(ParametersABC):
defaults["general"][k][k2] = v2 defaults["general"][k][k2] = v2
else: else:
defaults["general"][k] = v defaults["general"][k] = v
# default Tiler parameters
# define defaults and update with any inputs
defaults["tiler"] = TilerParameters.default(**tiler).to_dict() defaults["tiler"] = TilerParameters.default(**tiler).to_dict()
# generate a backup channel for when logfile meta is available
# Generate a backup channel, for when logfile meta is available
# but not image metadata. # but not image metadata.
backup_ref_channel = None backup_ref_channel = None
if "channels" in meta_d and isinstance( if "channels" in meta_d and isinstance(
...@@ -152,17 +151,14 @@ class PipelineParameters(ParametersABC): ...@@ -152,17 +151,14 @@ class PipelineParameters(ParametersABC):
defaults["tiler"]["ref_channel"] defaults["tiler"]["ref_channel"]
) )
defaults["tiler"]["backup_ref_channel"] = backup_ref_channel defaults["tiler"]["backup_ref_channel"] = backup_ref_channel
# default BABY parameters
defaults["baby"] = BabyParameters.default(**baby).to_dict() defaults["baby"] = BabyParameters.default(**baby).to_dict()
defaults["extraction"] = ( # default Extraction parmeters
exparams_from_meta(meta_d) defaults["extraction"] = exparams_from_meta(meta_d)
or BabyParameters.default(**extraction).to_dict() # default PostProcessing parameters
)
defaults["postprocessing"] = PostProcessorParameters.default( defaults["postprocessing"] = PostProcessorParameters.default(
**postprocessing **postprocessing
).to_dict() ).to_dict()
defaults["reporting"] = {}
return cls(**{k: v for k, v in defaults.items()}) return cls(**{k: v for k, v in defaults.items()})
def load_logs(self): def load_logs(self):
...@@ -187,7 +183,7 @@ class Pipeline(ProcessABC): ...@@ -187,7 +183,7 @@ class Pipeline(ProcessABC):
"postprocessing", "postprocessing",
] ]
# Specify the group in the h5 files written by each step # specify the group in the h5 files written by each step
writer_groups = { writer_groups = {
"tiler": ["trap_info"], "tiler": ["trap_info"],
"baby": ["cell_info"], "baby": ["cell_info"],
...@@ -230,8 +226,7 @@ class Pipeline(ProcessABC): ...@@ -230,8 +226,7 @@ class Pipeline(ProcessABC):
@classmethod @classmethod
def from_yaml(cls, fpath): def from_yaml(cls, fpath):
# This is just a convenience function, think before implementing # an unfinished convenience function
# for other processes
return cls(parameters=PipelineParameters.from_yaml(fpath)) return cls(parameters=PipelineParameters.from_yaml(fpath))
@classmethod @classmethod
...@@ -304,11 +299,16 @@ class Pipeline(ProcessABC): ...@@ -304,11 +299,16 @@ class Pipeline(ProcessABC):
def run(self): def run(self):
"""Run separate pipelines for all positions in an experiment.""" """Run separate pipelines for all positions in an experiment."""
# general information in config # display configuration
config = self.parameters.to_dict() config = self.parameters.to_dict()
for step in config:
print("\n---\n" + step + "\n---")
pprint(config[step])
print()
# extract from configuration
expt_id = config["general"]["id"] expt_id = config["general"]["id"]
distributed = config["general"]["distributed"] distributed = config["general"]["distributed"]
pos_filter = config["general"]["filter"] position_filter = config["general"]["filter"]
root_dir = Path(config["general"]["directory"]) root_dir = Path(config["general"]["directory"])
self.server_info = { self.server_info = {
k: config["general"].get(k) k: config["general"].get(k)
...@@ -320,56 +320,73 @@ class Pipeline(ProcessABC): ...@@ -320,56 +320,73 @@ class Pipeline(ProcessABC):
) )
# get log files, either locally or via OMERO # get log files, either locally or via OMERO
with dispatcher as conn: with dispatcher as conn:
image_ids = conn.get_images() position_ids = conn.get_images()
directory = self.store or root_dir / conn.unique_name directory = self.store or root_dir / conn.unique_name
if not directory.exists(): if not directory.exists():
directory.mkdir(parents=True) directory.mkdir(parents=True)
# download logs to use for metadata # download logs to use for metadata
conn.cache_logs(directory) conn.cache_logs(directory)
print("Positions available:")
for i, pos in enumerate(position_ids.keys()):
print("\t" + f"{i}: " + pos.split(".")[0])
# update configuration # update configuration
self.parameters.general["directory"] = str(directory) self.parameters.general["directory"] = str(directory)
config["general"]["directory"] = directory config["general"]["directory"] = directory
self.setLogger(directory) self.setLogger(directory)
# pick particular images if desired # pick particular positions if desired
if pos_filter is not None: if position_filter is not None:
if isinstance(pos_filter, list): if isinstance(position_filter, list):
image_ids = { position_ids = {
k: v k: v
for filt in pos_filter for filt in position_filter
for k, v in self.apply_filter(image_ids, filt).items() for k, v in self.apply_filter(position_ids, filt).items()
} }
else: else:
image_ids = self.apply_filter(image_ids, pos_filter) position_ids = self.apply_filter(position_ids, position_filter)
assert len(image_ids), "No images to segment" if not len(position_ids):
# create pipelines raise Exception("No images to segment.")
else:
print("\nPositions selected:")
for pos in position_ids:
print("\t" + pos.split(".")[0])
# create and run pipelines
if distributed != 0: if distributed != 0:
# multiple cores # multiple cores
with Pool(distributed) as p: with Pool(distributed) as p:
results = p.map( results = p.map(
lambda x: self.run_one_position(*x), lambda x: self.run_one_position(*x),
[(k, i) for i, k in enumerate(image_ids.items())], [(k, i) for i, k in enumerate(position_ids.items())],
) )
else: else:
# single core # single core
results = [] results = [
for k, v in tqdm(image_ids.items()): self.run_one_position((k, v), 1)
r = self.run_one_position((k, v), 1) for k, v in tqdm(position_ids.items())
results.append(r) ]
return results return results
def apply_filter(self, image_ids: dict, filt: int or str): def apply_filter(self, position_ids: dict, position_filter: int or str):
"""Select images by picking a particular one or by using a regular expression to parse their file names.""" """
if isinstance(filt, str): Select positions.
# pick images using a regular expression
image_ids = { Either pick a particular position or use a regular expression
k: v for k, v in image_ids.items() if re.search(filt, k) to parse their file names.
"""
if isinstance(position_filter, str):
# pick positions using a regular expression
position_ids = {
k: v
for k, v in position_ids.items()
if re.search(position_filter, k)
} }
elif isinstance(filt, int): elif isinstance(position_filter, int):
# pick the filt'th image # pick a particular position
image_ids = { position_ids = {
k: v for i, (k, v) in enumerate(image_ids.items()) if i == filt k: v
for i, (k, v) in enumerate(position_ids.items())
if i == position_filter
} }
return image_ids return position_ids
def run_one_position( def run_one_position(
self, self,
...@@ -379,120 +396,98 @@ class Pipeline(ProcessABC): ...@@ -379,120 +396,98 @@ class Pipeline(ProcessABC):
"""Set up and run a pipeline for one position.""" """Set up and run a pipeline for one position."""
self._pool_index = index self._pool_index = index
name, image_id = name_image_id name, image_id = name_image_id
# session and filename are defined by calling setup_pipeline. # session is defined by calling pipe_pipeline.
# can they be deleted here? # can it be deleted here?
session = None session = None
filename = None run_kwargs = {"extraction": {"cell_labels": None, "masks": None}}
#
run_kwargs = {"extraction": {"labels": None, "masks": None}}
try: try:
( pipe, session = self.setup_pipeline(image_id)
filename,
meta,
config,
process_from,
tps,
steps,
earlystop,
session,
trackers_state,
) = self._setup_pipeline(image_id)
loaded_writers = { loaded_writers = {
name: writer(filename) name: writer(pipe["filename"])
for k in self.step_sequence for k in self.step_sequence
if k in self.writers if k in self.writers
for name, writer in self.writers[k] for name, writer in self.writers[k]
} }
writer_ow_kwargs = { writer_overwrite_kwargs = {
"state": loaded_writers["state"].datatypes.keys(), "state": loaded_writers["state"].datatypes.keys(),
"baby": ["mother_assign"], "baby": ["mother_assign"],
} }
# START PIPELINE # START PIPELINE
frac_clogged_traps = 0.0 frac_clogged_traps = 0.0
min_process_from = min(process_from.values()) min_process_from = min(pipe["process_from"].values())
with dispatch_image(image_id)( with dispatch_image(image_id)(
image_id, **self.server_info image_id, **self.server_info
) as image: ) as image:
# initialise steps # initialise steps
if "tiler" not in steps: if "tiler" not in pipe["steps"]:
steps["tiler"] = Tiler.from_image( pipe["steps"]["tiler"] = Tiler.from_image(
image, TilerParameters.from_dict(config["tiler"]) image,
TilerParameters.from_dict(pipe["config"]["tiler"]),
) )
if process_from["baby"] < tps: if pipe["process_from"]["baby"] < pipe["tps"]:
session = initialise_tf(2) session = initialise_tf(2)
steps["baby"] = BabyRunner.from_tiler( pipe["steps"]["baby"] = BabyRunner.from_tiler(
BabyParameters.from_dict(config["baby"]), BabyParameters.from_dict(pipe["config"]["baby"]),
steps["tiler"], pipe["steps"]["tiler"],
) )
if trackers_state: if pipe["trackers_state"]:
steps["baby"].crawler.tracker_states = trackers_state pipe["steps"]["baby"].crawler.tracker_states = pipe[
# limit extraction parameters using the available channels in tiler "trackers_state"
if process_from["extraction"] < tps: ]
# TODO Move this parameter validation into Extractor if pipe["process_from"]["extraction"] < pipe["tps"]:
av_channels = set((*steps["tiler"].channels, "general"))
config["extraction"]["tree"] = {
k: v
for k, v in config["extraction"]["tree"].items()
if k in av_channels
}
config["extraction"]["sub_bg"] = av_channels.intersection(
config["extraction"]["sub_bg"]
)
av_channels_wsub = av_channels.union(
[c + "_bgsub" for c in config["extraction"]["sub_bg"]]
)
tmp = copy(config["extraction"]["multichannel_ops"])
for op, (input_ch, _, _) in tmp.items():
if not set(input_ch).issubset(av_channels_wsub):
del config["extraction"]["multichannel_ops"][op]
exparams = ExtractorParameters.from_dict( exparams = ExtractorParameters.from_dict(
config["extraction"] pipe["config"]["extraction"]
) )
steps["extraction"] = Extractor.from_tiler( pipe["steps"]["extraction"] = Extractor.from_tiler(
exparams, store=filename, tiler=steps["tiler"] exparams,
store=pipe["filename"],
tiler=pipe["steps"]["tiler"],
) )
# set up progress meter # initiate progress bar
pbar = tqdm( pbar = tqdm(
range(min_process_from, tps), range(min_process_from, pipe["tps"]),
desc=image.name, desc=image.name,
initial=min_process_from, initial=min_process_from,
total=tps, total=pipe["tps"],
) )
# run through time points
for i in pbar: for i in pbar:
if ( if (
frac_clogged_traps frac_clogged_traps
< earlystop["thresh_pos_clogged"] < pipe["earlystop"]["thresh_pos_clogged"]
or i < earlystop["min_tp"] or i < pipe["earlystop"]["min_tp"]
): ):
# run through steps # run through steps
for step in self.pipeline_steps: for step in self.pipeline_steps:
if i >= process_from[step]: if i >= pipe["process_from"][step]:
result = steps[step].run_tp( # perform step
result = pipe["steps"][step].run_tp(
i, **run_kwargs.get(step, {}) i, **run_kwargs.get(step, {})
) )
# write to h5 file using writers
# extractor writes to h5 itself
if step in loaded_writers: if step in loaded_writers:
loaded_writers[step].write( loaded_writers[step].write(
data=result, data=result,
overwrite=writer_ow_kwargs.get( overwrite=writer_overwrite_kwargs.get(
step, [] step, []
), ),
tp=i, tp=i,
meta={"last_processed": i}, meta={"last_processed": i},
) )
# perform step # clean up
if ( if (
step == "tiler" step == "tiler"
and i == min_process_from and i == min_process_from
): ):
logging.getLogger("aliby").info( logging.getLogger("aliby").info(
f"Found {steps['tiler'].n_tiles} traps in {image.name}" f"Found {pipe['steps']['tiler'].n_tiles} traps in {image.name}"
) )
elif step == "baby": elif step == "baby":
# write state and pass info to Extractor # write state
loaded_writers["state"].write( loaded_writers["state"].write(
data=steps[ data=pipe["steps"][
step step
].crawler.tracker_states, ].crawler.tracker_states,
overwrite=loaded_writers[ overwrite=loaded_writers[
...@@ -501,12 +496,14 @@ class Pipeline(ProcessABC): ...@@ -501,12 +496,14 @@ class Pipeline(ProcessABC):
tp=i, tp=i,
) )
elif step == "extraction": elif step == "extraction":
# remove mask/label after extraction # remove masks and labels after extraction
for k in ["masks", "labels"]: for k in ["masks", "cell_labels"]:
run_kwargs[step][k] = None run_kwargs[step][k] = None
# check and report clogging # check and report clogging
frac_clogged_traps = self.check_earlystop( frac_clogged_traps = self.check_earlystop(
filename, earlystop, steps["tiler"].tile_size pipe["filename"],
pipe["earlystop"],
pipe["steps"]["tiler"].tile_size,
) )
if frac_clogged_traps > 0.3: if frac_clogged_traps > 0.3:
self._log( self._log(
...@@ -519,18 +516,17 @@ class Pipeline(ProcessABC): ...@@ -519,18 +516,17 @@ class Pipeline(ProcessABC):
self._log( self._log(
f"{name}:Stopped early at time {i} with {frac_clogged_traps} clogged traps" f"{name}:Stopped early at time {i} with {frac_clogged_traps} clogged traps"
) )
meta.add_fields({"end_status": "Clogged"}) pipe["meta"].add_fields({"end_status": "Clogged"})
break break
meta.add_fields({"last_processed": i}) pipe["meta"].add_fields({"last_processed": i})
pipe["meta"].add_fields({"end_status": "Success"})
# run post-processing # run post-processing
meta.add_fields({"end_status": "Success"})
post_proc_params = PostProcessorParameters.from_dict( post_proc_params = PostProcessorParameters.from_dict(
config["postprocessing"] pipe["config"]["postprocessing"]
) )
PostProcessor(filename, post_proc_params).run() PostProcessor(pipe["filename"], post_proc_params).run()
self._log("Analysis finished successfully.", "info") self._log("Analysis finished successfully.", "info")
return 1 return 1
except Exception as e: except Exception as e:
# catch bugs during setup or run time # catch bugs during setup or run time
logging.exception( logging.exception(
...@@ -541,87 +537,9 @@ class Pipeline(ProcessABC): ...@@ -541,87 +537,9 @@ class Pipeline(ProcessABC):
traceback.print_exc() traceback.print_exc()
raise e raise e
finally: finally:
_close_session(session) close_session(session)
@staticmethod
def check_earlystop(filename: str, es_parameters: dict, tile_size: int):
"""
Check recent time points for tiles with too many cells.
Returns the fraction of clogged tiles, where clogged tiles have
too many cells or too much of their area covered by cells.
Parameters def setup_pipeline(
----------
filename: str
Name of h5 file.
es_parameters: dict
Parameters defining when early stopping should happen.
For example:
{'min_tp': 100,
'thresh_pos_clogged': 0.4,
'thresh_trap_ncells': 8,
'thresh_trap_area': 0.9,
'ntps_to_eval': 5}
tile_size: int
Size of tile.
"""
# get the area of the cells organised by trap and cell number
s = Signal(filename)
df = s.get_raw("/extraction/general/None/area")
# check the latest time points only
cells_used = df[
df.columns[-1 - es_parameters["ntps_to_eval"] : -1]
].dropna(how="all")
# find tiles with too many cells
traps_above_nthresh = (
cells_used.groupby("trap").count().apply(np.mean, axis=1)
> es_parameters["thresh_trap_ncells"]
)
# find tiles with cells covering too great a fraction of the tiles' area
traps_above_athresh = (
cells_used.groupby("trap").sum().apply(np.mean, axis=1)
/ tile_size**2
> es_parameters["thresh_trap_area"]
)
return (traps_above_nthresh & traps_above_athresh).mean()
# FIXME: Remove this functionality. It used to be for
# older hdf5 file formats.
def _load_config_from_file(
self,
filename: Path,
process_from: t.Dict[str, int],
trackers_state: t.List,
overwrite: t.Dict[str, bool],
):
with h5py.File(filename, "r") as f:
for k in process_from.keys():
if not overwrite[k]:
process_from[k] = self.legacy_get_last_tp[k](f)
process_from[k] += 1
return process_from, trackers_state, overwrite
# FIXME: Remove this functionality. It used to be for
# older hdf5 file formats.
@staticmethod
def legacy_get_last_tp(step: str) -> t.Callable:
"""Get last time-point in different ways depending
on which step we are using
To support segmentation in aliby < v0.24
TODO Deprecate and replace with State method
"""
switch_case = {
"tiler": lambda f: f["trap_info/drifts"].shape[0] - 1,
"baby": lambda f: f["cell_info/timepoint"][-1],
"extraction": lambda f: f[
"extraction/general/None/area/timepoint"
][-1],
}
return switch_case[step]
def _setup_pipeline(
self, image_id: int self, image_id: int
) -> t.Tuple[ ) -> t.Tuple[
Path, Path,
...@@ -645,83 +563,97 @@ class Pipeline(ProcessABC): ...@@ -645,83 +563,97 @@ class Pipeline(ProcessABC):
Returns Returns
------- -------
filename: str pipe: dict
Path to a h5 file to write to. With keys
meta: object filename: str
agora.io.metadata.MetaData object Path to a h5 file to write to.
config: dict meta: object
Configuration parameters. agora.io.metadata.MetaData object
process_from: dict config: dict
Gives from which time point each step of the pipeline should start. Configuration parameters.
tps: int process_from: dict
Number of time points. Gives time points from which each step of the
steps: dict pipeline should start.
earlystop: dict tps: int
Parameters to check whether the pipeline should be stopped. Number of time points.
steps: dict
earlystop: dict
Parameters to check whether the pipeline should
be stopped.
trackers_state: list
States of any trackers from earlier runs.
session: None session: None
trackers_state: list
States of any trackers from earlier runs.
""" """
pipe = {}
config = self.parameters.to_dict() config = self.parameters.to_dict()
# TODO Alan: Verify if session must be passed # TODO Alan: Verify if session must be passed
session = None session = None
earlystop = config["general"].get("earlystop", None) pipe["earlystop"] = config["general"].get("earlystop", None)
process_from = {k: 0 for k in self.pipeline_steps} pipe["process_from"] = {k: 0 for k in self.pipeline_steps}
steps = {} pipe["steps"] = {}
# check overwriting # check overwriting
ow_id = config["general"].get("overwrite", 0) overwrite_id = config["general"].get("overwrite", 0)
ow = {step: True for step in self.step_sequence} overwrite = {step: True for step in self.step_sequence}
if ow_id and ow_id is not True: if overwrite_id and overwrite_id is not True:
ow = { overwrite = {
step: self.step_sequence.index(ow_id) < i step: self.step_sequence.index(overwrite_id) < i
for i, step in enumerate(self.step_sequence, 1) for i, step in enumerate(self.step_sequence, 1)
} }
# set up
# Set up
directory = config["general"]["directory"] directory = config["general"]["directory"]
pipe["trackers_state"] = []
trackers_state: t.List[np.ndarray] = []
with dispatch_image(image_id)(image_id, **self.server_info) as image: with dispatch_image(image_id)(image_id, **self.server_info) as image:
filename = Path(f"{directory}/{image.name}.h5") pipe["filename"] = Path(f"{directory}/{image.name}.h5")
meta = MetaData(directory, filename) # load metadata from h5 file
from_start = True if np.any(ow.values()) else False pipe["meta"] = MetaData(directory, pipe["filename"])
# remove existing file if overwriting from_start = True if np.any(overwrite.values()) else False
# remove existing h5 file if overwriting
if ( if (
from_start from_start
and ( and (
config["general"].get("overwrite", False) config["general"].get("overwrite", False)
or np.all(list(ow.values())) or np.all(list(overwrite.values()))
) )
and filename.exists() and pipe["filename"].exists()
): ):
os.remove(filename) os.remove(pipe["filename"])
# if the file exists with no previous segmentation use its tiler # if the file exists with no previous segmentation use its tiler
if filename.exists(): if pipe["filename"].exists():
self._log("Result file exists.", "info") self._log("Result file exists.", "info")
if not ow["tiler"]: if not overwrite["tiler"]:
steps["tiler"] = Tiler.from_hdf5(image, filename) pipe["steps"]["tiler"] = Tiler.from_h5(
image, pipe["filename"]
)
try: try:
( (
process_from, process_from,
trackers_state, trackers_state,
ow, overwrite,
) = self._load_config_from_file( ) = self._load_config_from_file(
filename, process_from, trackers_state, ow pipe["filename"],
pipe["process_from"],
pipe["trackers_state"],
overwrite,
) )
# get state array # get state array
trackers_state = ( pipe["trackers_state"] = (
[] []
if ow["baby"] if overwrite["baby"]
else StateReader(filename).get_formatted_states() else StateReader(
pipe["filename"]
).get_formatted_states()
) )
config["tiler"] = steps["tiler"].parameters.to_dict() config["tiler"] = pipe["steps"][
"tiler"
].parameters.to_dict()
except Exception: except Exception:
self._log(f"Overwriting tiling data") self._log("Overwriting tiling data")
if config["general"]["use_explog"]: if config["general"]["use_explog"]:
meta.run() pipe["meta"].run()
pipe["config"] = config
# add metadata not in the log file # add metadata not in the log file
meta.add_fields( pipe["meta"].add_fields(
{ {
"aliby_version": version("aliby"), "aliby_version": version("aliby"),
"baby_version": version("aliby-baby"), "baby_version": version("aliby-baby"),
...@@ -734,20 +666,53 @@ class Pipeline(ProcessABC): ...@@ -734,20 +666,53 @@ class Pipeline(ProcessABC):
).to_yaml(), ).to_yaml(),
} }
) )
tps = min(config["general"]["tps"], image.data.shape[0]) pipe["tps"] = min(config["general"]["tps"], image.data.shape[0])
return ( return pipe, session
filename,
meta, @staticmethod
config, def check_earlystop(filename: str, es_parameters: dict, tile_size: int):
process_from, """
tps, Check recent time points for tiles with too many cells.
steps,
earlystop, Returns the fraction of clogged tiles, where clogged tiles have
session, too many cells or too much of their area covered by cells.
trackers_state,
) Parameters
----------
filename: str
Name of h5 file.
es_parameters: dict
Parameters defining when early stopping should happen.
For example:
{'min_tp': 100,
'thresh_pos_clogged': 0.4,
'thresh_trap_ncells': 8,
'thresh_trap_area': 0.9,
'ntps_to_eval': 5}
tile_size: int
Size of tile.
"""
# get the area of the cells organised by trap and cell number
s = Signal(filename)
df = s.get_raw("/extraction/general/None/area")
# check the latest time points only
cells_used = df[
df.columns[-1 - es_parameters["ntps_to_eval"] : -1]
].dropna(how="all")
# find tiles with too many cells
traps_above_nthresh = (
cells_used.groupby("trap").count().apply(np.mean, axis=1)
> es_parameters["thresh_trap_ncells"]
)
# find tiles with cells covering too great a fraction of the tiles' area
traps_above_athresh = (
cells_used.groupby("trap").sum().apply(np.mean, axis=1)
/ tile_size**2
> es_parameters["thresh_trap_area"]
)
return (traps_above_nthresh & traps_above_athresh).mean()
def _close_session(session): def close_session(session):
if session: if session:
session.close() session.close()
src/aliby/tile/tiler.py
View file @ b95e79fc
""" """
Tiler: Divides images into smaller tiles. Tiler: Divides images into smaller tiles.
The tasks of the Tiler are selecting regions of interest, or tiles, of images - with one trap per tile, correcting for the drift of the microscope stage over time, and handling errors and bridging between the image data and Aliby’s image-processing steps. The tasks of the Tiler are selecting regions of interest, or tiles, of
images - with one trap per tile, correcting for the drift of the microscope
stage over time, and handling errors and bridging between the image data
and Aliby’s image-processing steps.
Tiler subclasses deal with either network connections or local files. Tiler subclasses deal with either network connections or local files.
To find tiles, we use a two-step process: we analyse the bright-field image to produce the template of a trap, and we fit this template to the image to find the tiles' centres. To find tiles, we use a two-step process: we analyse the bright-field image
to produce the template of a trap, and we fit this template to the image to
find the tiles' centres.
We use texture-based segmentation (entropy) to split the image into foreground -- cells and traps -- and background, which we then identify with an Otsu filter. Two methods are used to produce a template trap from these regions: pick the trap with the smallest minor axis length and average over all validated traps. We use texture-based segmentation (entropy) to split the image into
foreground -- cells and traps -- and background, which we then identify with
an Otsu filter. Two methods are used to produce a template trap from these
regions: pick the trap with the smallest minor axis length and average over
all validated traps.
A peak-identifying algorithm recovers the x and y-axis location of traps in the original image, and we choose the approach to template that identifies the most tiles. A peak-identifying algorithm recovers the x and y-axis location of traps in
the original image, and we choose the approach to template that identifies
the most tiles.
The experiment is stored as an array with a standard indexing order of (Time, Channels, Z-stack, X, Y). The experiment is stored as an array with a standard indexing order of
(Time, Channels, Z-stack, X, Y).
""" """
import logging import logging
import re import re
...@@ -169,7 +181,7 @@ class TileLocations: ...@@ -169,7 +181,7 @@ class TileLocations:
return cls(initial_location, tile_size, max_size, drifts=[]) return cls(initial_location, tile_size, max_size, drifts=[])
@classmethod @classmethod
def read_hdf5(cls, file): def read_h5(cls, file):
"""Instantiate from a h5 file.""" """Instantiate from a h5 file."""
with h5py.File(file, "r") as hfile: with h5py.File(file, "r") as hfile:
tile_info = hfile["trap_info"] tile_info = hfile["trap_info"]
...@@ -234,12 +246,9 @@ class Tiler(StepABC): ...@@ -234,12 +246,9 @@ class Tiler(StepABC):
"channels", "channels",
list(range(metadata.get("size_c", 0))), list(range(metadata.get("size_c", 0))),
) )
self.ref_channel = self.get_channel_index(parameters.ref_channel) self.ref_channel = self.get_channel_index(parameters.ref_channel)
if self.ref_channel is None: if self.ref_channel is None:
self.ref_channel = self.backup_ref_channel self.ref_channel = self.backup_ref_channel
self.ref_channel = self.get_channel_index(parameters.ref_channel)
self.tile_locs = tile_locs self.tile_locs = tile_locs
try: try:
self.z_perchannel = { self.z_perchannel = {
...@@ -316,7 +325,7 @@ class Tiler(StepABC): ...@@ -316,7 +325,7 @@ class Tiler(StepABC):
Path to a directory of h5 files Path to a directory of h5 files
parameters: an instance of TileParameters (optional) parameters: an instance of TileParameters (optional)
""" """
tile_locs = TileLocations.read_hdf5(filepath) tile_locs = TileLocations.read_h5(filepath)
metadata = BridgeH5(filepath).meta_h5 metadata = BridgeH5(filepath).meta_h5
metadata["channels"] = image.metadata["channels"] metadata["channels"] = image.metadata["channels"]
if parameters is None: if parameters is None:
...@@ -332,7 +341,7 @@ class Tiler(StepABC): ...@@ -332,7 +341,7 @@ class Tiler(StepABC):
return tiler return tiler
@lru_cache(maxsize=2) @lru_cache(maxsize=2)
def get_tc(self, t: int, c: int) -> np.ndarray: def get_tc(self, tp: int, c: int) -> np.ndarray:
""" """
Load image using dask. Load image using dask.
...@@ -345,7 +354,7 @@ class Tiler(StepABC): ...@@ -345,7 +354,7 @@ class Tiler(StepABC):
Parameters Parameters
---------- ----------
t: integer tp: integer
An index for a time point An index for a time point
c: integer c: integer
An index for a channel An index for a channel
...@@ -354,10 +363,10 @@ class Tiler(StepABC): ...@@ -354,10 +363,10 @@ class Tiler(StepABC):
------- -------
full: an array of images full: an array of images
""" """
full = self.image[t, c] full = self.image[tp, c]
if hasattr(full, "compute"): # If using dask fetch images here if hasattr(full, "compute"):
# if using dask fetch images
full = full.compute(scheduler="synchronous") full = full.compute(scheduler="synchronous")
return full return full
@property @property
...@@ -558,9 +567,8 @@ class Tiler(StepABC): ...@@ -558,9 +567,8 @@ class Tiler(StepABC):
Returns Returns
------- -------
res: array res: array
Data arranged as (tiles, channels, time points, X, Y, Z) Data arranged as (tiles, channels, Z, X, Y)
""" """
# FIXME add support for sub-tiling a tile
# FIXME can we ignore z # FIXME can we ignore z
if channels is None: if channels is None:
channels = [0] channels = [0]
...@@ -571,8 +579,7 @@ class Tiler(StepABC): ...@@ -571,8 +579,7 @@ class Tiler(StepABC):
for c in channels: for c in channels:
# only return requested z # only return requested z
val = self.get_tp_data(tp, c)[:, z] val = self.get_tp_data(tp, c)[:, z]
# starts with the order: tiles, z, y, x # starts with the order: tiles, Z, Y, X
# returns the order: tiles, C, T, Z, X, Y
val = np.expand_dims(val, axis=1) val = np.expand_dims(val, axis=1)
res.append(val) res.append(val)
if tile_shape is not None: if tile_shape is not None:
...@@ -584,7 +591,10 @@ class Tiler(StepABC): ...@@ -584,7 +591,10 @@ class Tiler(StepABC):
for tile_size, ax in zip(tile_shape, res[0].shape[-3:-2]) for tile_size, ax in zip(tile_shape, res[0].shape[-3:-2])
] ]
) )
return np.stack(res, axis=1) # convert to array with channels as first column
# final has dimensions (tiles, channels, 1, Z, X, Y)
final = np.stack(res, axis=1)
return final
@property @property
def ref_channel_index(self): def ref_channel_index(self):
...@@ -593,7 +603,10 @@ class Tiler(StepABC): ...@@ -593,7 +603,10 @@ class Tiler(StepABC):
def get_channel_index(self, channel: str or int) -> int or None: def get_channel_index(self, channel: str or int) -> int or None:
""" """
Find index for channel using regex. Returns the first matched string. Find index for channel using regex.
Return the first matched string.
If self.channels is integers (no image metadata) it returns None. If self.channels is integers (no image metadata) it returns None.
If channel is integer If channel is integer
...@@ -602,10 +615,8 @@ class Tiler(StepABC): ...@@ -602,10 +615,8 @@ class Tiler(StepABC):
channel: string or int channel: string or int
The channel or index to be used. The channel or index to be used.
""" """
if all(map(lambda x: isinstance(x, int), self.channels)): if all(map(lambda x: isinstance(x, int), self.channels)):
channel = channel if isinstance(channel, int) else None channel = channel if isinstance(channel, int) else None
if isinstance(channel, str): if isinstance(channel, str):
channel = find_channel_index(self.channels, channel) channel = find_channel_index(self.channels, channel)
return channel return channel
......
src/extraction/core/extractor.py
View file @ b95e79fc
...@@ -26,14 +26,14 @@ extraction_result = t.Dict[ ...@@ -26,14 +26,14 @@ extraction_result = t.Dict[
str, t.Dict[reduction_method, t.Dict[str, t.Dict[str, pd.Series]]] str, t.Dict[reduction_method, t.Dict[str, t.Dict[str, pd.Series]]]
] ]
# Global variables used to load functions that either analyse cells or their background. These global variables both allow the functions to be stored in a dictionary for access only on demand and to be defined simply in extraction/core/functions. # Global variables used to load functions that either analyse cells
CELL_FUNS, TRAPFUNS, FUNS = load_funs() # or their background. These global variables both allow the functions
# to be stored in a dictionary for access only on demand and to be
# defined simply in extraction/core/functions.
CELL_FUNS, TRAP_FUNS, ALL_FUNS = load_funs()
CUSTOM_FUNS, CUSTOM_ARGS = load_custom_args() CUSTOM_FUNS, CUSTOM_ARGS = load_custom_args()
RED_FUNS = load_redfuns() RED_FUNS = load_redfuns()
# Assign datatype depending on the metric used
# m2type = {"mean": np.float32, "median": np.ubyte, "imBackground": np.ubyte}
class ExtractorParameters(ParametersABC): class ExtractorParameters(ParametersABC):
"""Base class to define parameters for extraction.""" """Base class to define parameters for extraction."""
...@@ -74,13 +74,19 @@ class Extractor(StepABC): ...@@ -74,13 +74,19 @@ class Extractor(StepABC):
""" """
Apply a metric to cells identified in the tiles. Apply a metric to cells identified in the tiles.
Using the cell masks, the Extractor applies a metric, such as area or median, to cells identified in the image tiles. Using the cell masks, the Extractor applies a metric, such as
area or median, to cells identified in the image tiles.
Its methods require both tile images and masks. Its methods require both tile images and masks.
Usually the metric is applied to only a tile's masked area, but some metrics depend on the whole tile. Usually the metric is applied to only a tile's masked area, but
some metrics depend on the whole tile.
Extraction follows a three-level tree structure. Channels, such as GFP, are the root level; the reduction algorithm, such as maximum projection, is the second level; the specific metric, or operation, to apply to the masks, such as mean, is the third level. Extraction follows a three-level tree structure. Channels, such
as GFP, are the root level; the reduction algorithm, such as
maximum projection, is the second level; the specific metric,
or operation, to apply to the masks, such as mean, is the third
or leaf level.
""" """
# TODO Alan: Move this to a location with the SwainLab defaults # TODO Alan: Move this to a location with the SwainLab defaults
...@@ -107,7 +113,8 @@ class Extractor(StepABC): ...@@ -107,7 +113,8 @@ class Extractor(StepABC):
store: str store: str
Path to the h5 file containing the cell masks. Path to the h5 file containing the cell masks.
tiler: pipeline-core.core.segmentation tiler tiler: pipeline-core.core.segmentation tiler
Class that contains or fetches the images used for segmentation. Class that contains or fetches the images used for
segmentation.
""" """
self.params = parameters self.params = parameters
if store: if store:
...@@ -118,6 +125,24 @@ class Extractor(StepABC): ...@@ -118,6 +125,24 @@ class Extractor(StepABC):
self.meta = {"channel": parameters.to_dict()["tree"].keys()} self.meta = {"channel": parameters.to_dict()["tree"].keys()}
if tiler: if tiler:
self.tiler = tiler self.tiler = tiler
available_channels = set((*tiler.channels, "general"))
# only extract for channels available
self.params.tree = {
k: v
for k, v in self.params.tree.items()
if k in available_channels
}
self.params.sub_bg = available_channels.intersection(
self.params.sub_bg
)
# add background subtracted channels to those available
available_channels_bgsub = available_channels.union(
[c + "_bgsub" for c in self.params.sub_bg]
)
# remove any multichannel operations requiring a missing channel
for op, (input_ch, _, _) in self.params.multichannel_ops.items():
if not set(input_ch).issubset(available_channels_bgsub):
self.params.multichannel_ops.pop(op)
self.load_funs() self.load_funs()
@classmethod @classmethod
...@@ -161,13 +186,16 @@ class Extractor(StepABC): ...@@ -161,13 +186,16 @@ class Extractor(StepABC):
def load_custom_funs(self): def load_custom_funs(self):
""" """
Incorporate the extra arguments of custom functions into their definitions. Incorporate the extra arguments of custom functions into their
definitions.
Normal functions only have cell_masks and trap_image as their Normal functions only have cell_masks and trap_image as their
arguments, and here custom functions are made the same by arguments, and here custom functions are made the same by
setting the values of their extra arguments. setting the values of their extra arguments.
Any other parameters are taken from the experiment's metadata and automatically applied. These parameters therefore must be loaded within an Extractor instance. Any other parameters are taken from the experiment's metadata
and automatically applied. These parameters therefore must be
loaded within an Extractor instance.
""" """
# find functions specified in params.tree # find functions specified in params.tree
funs = set( funs = set(
...@@ -202,11 +230,11 @@ class Extractor(StepABC): ...@@ -202,11 +230,11 @@ class Extractor(StepABC):
self._custom_funs[k] = tmp(f) self._custom_funs[k] = tmp(f)
def load_funs(self): def load_funs(self):
"""Define all functions, including custum ones.""" """Define all functions, including custom ones."""
self.load_custom_funs() self.load_custom_funs()
self._all_cell_funs = set(self._custom_funs.keys()).union(CELL_FUNS) self._all_cell_funs = set(self._custom_funs.keys()).union(CELL_FUNS)
# merge the two dicts # merge the two dicts
self._all_funs = {**self._custom_funs, **FUNS} self._all_funs = {**self._custom_funs, **ALL_FUNS}
def load_meta(self): def load_meta(self):
"""Load metadata from h5 file.""" """Load metadata from h5 file."""
...@@ -222,7 +250,8 @@ class Extractor(StepABC): ...@@ -222,7 +250,8 @@ class Extractor(StepABC):
""" """
Find tiles for a given time point, channels, and z-stacks. Find tiles for a given time point, channels, and z-stacks.
Any additional keyword arguments are passed to tiler.get_tiles_timepoint Any additional keyword arguments are passed to
tiler.get_tiles_timepoint
Parameters Parameters
---------- ----------
...@@ -243,8 +272,9 @@ class Extractor(StepABC): ...@@ -243,8 +272,9 @@ class Extractor(StepABC):
# a list of the indices of the z stacks # a list of the indices of the z stacks
channel_ids = None channel_ids = None
if z is None: if z is None:
# gets the tiles data via tiler # include all Z channels
z = list(range(self.tiler.shape[-3])) z = list(range(self.tiler.shape[-3]))
# get the image data via tiler
res = ( res = (
self.tiler.get_tiles_timepoint( self.tiler.get_tiles_timepoint(
tp, channels=channel_ids, z=z, **kwargs tp, channels=channel_ids, z=z, **kwargs
...@@ -252,15 +282,15 @@ class Extractor(StepABC): ...@@ -252,15 +282,15 @@ class Extractor(StepABC):
if channel_ids if channel_ids
else None else None
) )
# data arranged as (tiles, channels, time points, X, Y, Z) # res has dimensions (tiles, channels, 1, Z, X, Y)
return res return res
def extract_traps( def extract_traps(
self, self,
traps: t.List[np.ndarray], traps: t.List[np.ndarray],
masks: t.List[np.ndarray], masks: t.List[np.ndarray],
metric: str, cell_property: str,
labels: t.Dict[int, t.List[int]], cell_labels: t.Dict[int, t.List[int]],
) -> t.Tuple[t.Union[t.Tuple[float], t.Tuple[t.Tuple[int]]]]: ) -> t.Tuple[t.Union[t.Tuple[float], t.Tuple[t.Tuple[int]]]]:
""" """
Apply a function to a whole position. Apply a function to a whole position.
...@@ -271,35 +301,35 @@ class Extractor(StepABC): ...@@ -271,35 +301,35 @@ class Extractor(StepABC):
t.List of images. t.List of images.
masks: list of arrays masks: list of arrays
t.List of masks. t.List of masks.
metric: str cell_property: str
Metric to extract. Property to extract, including imBackground.
labels: dict cell_labels: dict
A dict of cell labels with trap_ids as keys and a list of cell labels as values. A dict of cell labels with trap_ids as keys and a list
pos_info: bool of cell labels as values.
Whether to add the position as an index or not.
Returns Returns
------- -------
res_idx: a tuple of tuples res_idx: a tuple of tuples
A two-tuple comprising a tuple of results and a tuple of the tile_id and cell labels A two-tuple comprising a tuple of results and a tuple of
the tile_id and cell labels
""" """
if labels is None: if cell_labels is None:
self._log("No labels given. Sorting cells using index.") self._log("No cell labels given. Sorting cells using index.")
cell_fun = True if metric in self._all_cell_funs else False cell_fun = True if cell_property in self._all_cell_funs else False
idx = [] idx = []
results = [] results = []
for trap_id, (mask_set, trap, lbl_set) in enumerate( for trap_id, (mask_set, trap, local_cell_labels) in enumerate(
zip(masks, traps, labels.values()) zip(masks, traps, cell_labels.values())
): ):
# ignore empty traps # ignore empty traps
if len(mask_set): if len(mask_set):
# apply metric either a cell function or otherwise # find property from the tile
result = self._all_funs[metric](mask_set, trap) result = self._all_funs[cell_property](mask_set, trap)
if cell_fun: if cell_fun:
# store results for each cell separately # store results for each cell separately
for lbl, val in zip(lbl_set, result): for cell_label, val in zip(local_cell_labels, result):
results.append(val) results.append(val)
idx.append((trap_id, lbl)) idx.append((trap_id, cell_label))
else: else:
# background (trap) function # background (trap) function
results.append(result) results.append(result)
...@@ -311,19 +341,19 @@ class Extractor(StepABC): ...@@ -311,19 +341,19 @@ class Extractor(StepABC):
self, self,
traps: t.List[np.array], traps: t.List[np.array],
masks: t.List[np.array], masks: t.List[np.array],
metrics: t.List[str], cell_properties: t.List[str],
**kwargs, **kwargs,
) -> t.Dict[str, pd.Series]: ) -> t.Dict[str, pd.Series]:
""" """
Return dict with metrics as key and metrics applied to data as values. Return dict with metrics as key and cell_properties as values.
Data from one time point is used. Data from one time point is used.
""" """
d = { d = {
metric: self.extract_traps( cell_property: self.extract_traps(
traps=traps, masks=masks, metric=metric, **kwargs traps=traps, masks=masks, cell_property=cell_property, **kwargs
) )
for metric in metrics for cell_property in cell_properties
} }
return d return d
...@@ -331,11 +361,11 @@ class Extractor(StepABC): ...@@ -331,11 +361,11 @@ class Extractor(StepABC):
self, self,
traps: np.ndarray, traps: np.ndarray,
masks: t.List[np.ndarray], masks: t.List[np.ndarray],
red_metrics: t.Dict[reduction_method, t.Collection[str]], tree_branch: t.Dict[reduction_method, t.Collection[str]],
**kwargs, **kwargs,
) -> t.Dict[str, t.Dict[reduction_method, t.Dict[str, pd.Series]]]: ) -> t.Dict[str, t.Dict[reduction_method, t.Dict[str, pd.Series]]]:
""" """
Wrapper to apply reduction and then extraction. Wrapper to reduce to a 2D image and then extract.
Parameters Parameters
---------- ----------
...@@ -343,8 +373,10 @@ class Extractor(StepABC): ...@@ -343,8 +373,10 @@ class Extractor(StepABC):
An array of image data arranged as (tiles, X, Y, Z) An array of image data arranged as (tiles, X, Y, Z)
masks: list of arrays masks: list of arrays
An array of masks for each trap: one per cell at the trap An array of masks for each trap: one per cell at the trap
red_metrics: dict tree_branch: dict
dict for which keys are reduction functions and values are either a list or a set of strings giving the metric functions. An upper branch of the extraction tree: a dict for which
keys are reduction functions and values are either a list
or a set of strings giving the cell properties to be found.
For example: {'np_max': {'max5px', 'mean', 'median'}} For example: {'np_max': {'max5px', 'mean', 'median'}}
**kwargs: dict **kwargs: dict
All other arguments passed to Extractor.extract_funs. All other arguments passed to Extractor.extract_funs.
...@@ -353,22 +385,27 @@ class Extractor(StepABC): ...@@ -353,22 +385,27 @@ class Extractor(StepABC):
------ ------
Dict of dataframes with the corresponding reductions and metrics nested. Dict of dataframes with the corresponding reductions and metrics nested.
""" """
# create dict with keys naming the reduction in the z-direction and the reduced data as values # FIXME hack to pass tests
reduced_tiles_data = {} if "labels" in kwargs:
kwargs["cell_labels"] = kwargs.pop("labels")
# create dict with keys naming the reduction in the z-direction
# and the reduced data as values
reduced_tiles = {}
if traps is not None: if traps is not None:
for red_fun in red_metrics.keys(): for red_fun in tree_branch.keys():
reduced_tiles_data[red_fun] = [ reduced_tiles[red_fun] = [
self.reduce_dims(tile_data, method=RED_FUNS[red_fun]) self.reduce_dims(tile_data, method=RED_FUNS[red_fun])
for tile_data in traps for tile_data in traps
] ]
# calculate cell and tile properties
d = { d = {
red_fun: self.extract_funs( red_fun: self.extract_funs(
metrics=metrics, cell_properties=cell_properties,
traps=reduced_tiles_data.get(red_fun, [None for _ in masks]), traps=reduced_tiles.get(red_fun, [None for _ in masks]),
masks=masks, masks=masks,
**kwargs, **kwargs,
) )
for red_fun, metrics in red_metrics.items() for red_fun, cell_properties in tree_branch.items()
} }
return d return d
...@@ -392,64 +429,24 @@ class Extractor(StepABC): ...@@ -392,64 +429,24 @@ class Extractor(StepABC):
reduced = reduce_z(img, method) reduced = reduce_z(img, method)
return reduced return reduced
def extract_tp( def make_tree_bits(self, tree):
self, """Put extraction tree and information for the channels into a dict."""
tp: int,
tree: t.Optional[extraction_tree] = None,
tile_size: int = 117,
masks: t.Optional[t.List[np.ndarray]] = None,
labels: t.Optional[t.List[int]] = None,
**kwargs,
) -> t.Dict[str, t.Dict[str, t.Dict[str, tuple]]]:
"""
Extract for an individual time point.
Parameters
----------
tp : int
Time point being analysed.
tree : dict
Nested dictionary indicating channels, reduction functions and
metrics to be used.
For example: {'general': {'None': ['area', 'volume', 'eccentricity']}}
tile_size : int
Size of the tile to be extracted.
masks : list of arrays
A list of masks per trap with each mask having dimensions (ncells, tile_size,
tile_size).
labels : dict
A dictionary with trap_ids as keys and cell_labels as values.
**kwargs : keyword arguments
Passed to extractor.reduce_extract.
Returns
-------
d: dict
Dictionary of the results with three levels of dictionaries.
The first level has channels as keys.
The second level has reduction metrics as keys.
The third level has cell or background metrics as keys and a two-tuple as values.
The first tuple is the result of applying the metrics to a particular cell or trap; the second tuple is either (trap_id, cell_label) for a metric applied to a cell or a trap_id for a metric applied to a trap.
An example is d["GFP"]["np_max"]["mean"][0], which gives a tuple of the calculated mean GFP fluorescence for all cells.
"""
# TODO Can we split the different extraction types into sub-methods to make this easier to read?
if tree is None: if tree is None:
# use default # use default
tree: extraction_tree = self.params.tree tree: extraction_tree = self.params.tree
# dictionary with channel: {reduction algorithm : metric} tree_bits = {
ch_tree = {ch: v for ch, v in tree.items() if ch != "general"} "tree": tree,
# tuple of the channels # dictionary with channel: {reduction algorithm : metric}
tree_chs = (*ch_tree,) "channel_tree": {
# create a Cells object to extract information from the h5 file ch: v for ch, v in tree.items() if ch != "general"
cells = Cells(self.local) },
# find the cell labels and store as dict with trap_ids as keys }
if labels is None: # tuple of the fluorescence channels
raw_labels = cells.labels_at_time(tp) tree_bits["tree_channels"] = (*tree_bits["channel_tree"],)
labels = { return tree_bits
trap_id: raw_labels.get(trap_id, [])
for trap_id in range(cells.ntraps) def get_masks(self, tp, masks, cells):
} """Get the masks as a list with an array of masks for each trap."""
# find the cell masks for a given trap as a dict with trap_ids as keys # find the cell masks for a given trap as a dict with trap_ids as keys
if masks is None: if masks is None:
raw_masks = cells.at_time(tp, kind="mask") raw_masks = cells.at_time(tp, kind="mask")
...@@ -458,16 +455,31 @@ class Extractor(StepABC): ...@@ -458,16 +455,31 @@ class Extractor(StepABC):
if len(cells): if len(cells):
masks[trap_id] = np.stack(np.array(cells)).astype(bool) masks[trap_id] = np.stack(np.array(cells)).astype(bool)
# convert to a list of masks # convert to a list of masks
# one array of size (no cells, tile_size, tile_size) per trap
masks = [np.array(v) for v in masks.values()] masks = [np.array(v) for v in masks.values()]
# find image data at the time point return masks
# stored as an array arranged as (traps, channels, time points, X, Y, Z)
tiles = self.get_tiles(tp, tile_shape=tile_size, channels=tree_chs) def get_cell_labels(self, tp, cell_labels, cells):
# generate boolean masks for background as a list with one mask per trap """Get the cell labels per trap as a dict with trap_ids as keys."""
bgs = np.array([]) if cell_labels is None:
raw_cell_labels = cells.labels_at_time(tp)
cell_labels = {
trap_id: raw_cell_labels.get(trap_id, [])
for trap_id in range(cells.ntraps)
}
return cell_labels
def get_background_masks(self, masks, tile_size):
"""
Generate boolean background masks.
Combine masks per trap and then take the logical inverse.
"""
if self.params.sub_bg: if self.params.sub_bg:
bgs = ~np.array( bgs = ~np.array(
list( list(
map( map(
# sum over masks for each cell
lambda x: np.sum(x, axis=0) lambda x: np.sum(x, axis=0)
if np.any(x) if np.any(x)
else np.zeros((tile_size, tile_size)), else np.zeros((tile_size, tile_size)),
...@@ -475,77 +487,179 @@ class Extractor(StepABC): ...@@ -475,77 +487,179 @@ class Extractor(StepABC):
) )
) )
).astype(bool) ).astype(bool)
# perform extraction by applying metrics else:
bgs = np.array([])
return bgs
def extract_one_channel(
self, tree_bits, cell_labels, tiles, masks, bgs, **kwargs
):
"""
Extract using all metrics requiring a single channel.
Apply first without and then with background subtraction.
Return the extraction results and a dict of background
corrected images.
"""
d = {} d = {}
self.img_bgsub = {} img_bgsub = {}
for ch, red_metrics in tree.items(): for ch, tree_branch in tree_bits["tree"].items():
# NB ch != is necessary for threading # NB ch != is necessary for threading
if ch != "general" and tiles is not None and len(tiles): if ch != "general" and tiles is not None and len(tiles):
# image data for all traps and z sections for a particular channel # image data for all traps for a particular channel and time point
# as an array arranged as (tiles, Z, X, Y, ) # arranged as (traps, Z, X, Y)
img = tiles[:, tree_chs.index(ch), 0] # we use 0 here to access the single time point available
img = tiles[:, tree_bits["tree_channels"].index(ch), 0]
else: else:
# no reduction applied to bright-field images
img = None img = None
# apply metrics to image data # apply metrics to image data
d[ch] = self.reduce_extract( d[ch] = self.reduce_extract(
traps=img, traps=img,
masks=masks, masks=masks,
red_metrics=red_metrics, tree_branch=tree_branch,
labels=labels, cell_labels=cell_labels,
**kwargs, **kwargs,
) )
# apply metrics to image data with the background subtracted # apply metrics to image data with the background subtracted
if bgs.any() and ch in self.params.sub_bg and img is not None: if bgs.any() and ch in self.params.sub_bg and img is not None:
# calculate metrics with subtracted bg # calculate metrics with background subtracted
ch_bs = ch + "_bgsub" ch_bs = ch + "_bgsub"
# subtract median background # subtract median background
bgsub_mapping = map(
self.img_bgsub[ch_bs] = np.moveaxis( # move Z to last column to allow subtraction
np.stack( lambda img, bgs: np.moveaxis(img, 0, -1)
list( # median of background over all pixels for each Z section
map( - bn.median(img[:, bgs], axis=1),
lambda tile, mask: np.moveaxis(tile, 0, -1) img,
- bn.median(tile[:, mask], axis=1), bgs,
img, )
bgs, # apply map and convert to array
) mapping_result = np.stack(list(bgsub_mapping))
) # move Z axis back to the second column
), img_bgsub[ch_bs] = np.moveaxis(mapping_result, -1, 1)
-1,
1,
) # End with tiles, z, y, x
# apply metrics to background-corrected data # apply metrics to background-corrected data
d[ch_bs] = self.reduce_extract( d[ch_bs] = self.reduce_extract(
red_metrics=ch_tree[ch], tree_branch=tree_bits["channel_tree"][ch],
traps=self.img_bgsub[ch_bs], traps=img_bgsub[ch_bs],
masks=masks, masks=masks,
labels=labels, cell_labels=cell_labels,
**kwargs, **kwargs,
) )
# apply any metrics using multiple channels, such as pH calculations return d, img_bgsub
def extract_multiple_channels(
self, tree_bits, cell_labels, tiles, masks, **kwargs
):
"""
Extract using all metrics requiring multiple channels.
"""
available_chs = set(self.img_bgsub.keys()).union(
tree_bits["tree_channels"]
)
d = {}
for name, ( for name, (
chs, chs,
merge_fun, reduction_fun,
red_metrics, op,
) in self.params.multichannel_ops.items(): ) in self.params.multichannel_ops.items():
if len( common_chs = set(chs).intersection(available_chs)
set(chs).intersection( # all required channels should be available
set(self.img_bgsub.keys()).union(tree_chs) if len(common_chs) == len(chs):
)
) == len(chs):
channels_stack = np.stack( channels_stack = np.stack(
[self.get_imgs(ch, tiles, tree_chs) for ch in chs], axis=-1 [
self.get_imgs(ch, tiles, tree_bits["tree_channels"])
for ch in chs
],
axis=-1,
) )
merged = RED_FUNS[merge_fun](channels_stack, axis=-1) # reduce in Z
d[name] = self.reduce_extract( traps = RED_FUNS[reduction_fun](channels_stack, axis=1)
red_metrics=red_metrics, # evaluate multichannel op
traps=merged, if name not in d:
masks=masks, d[name] = {}
labels=labels, if reduction_fun not in d[name]:
**kwargs, d[name][reduction_fun] = {}
d[name][reduction_fun][op] = self.extract_traps(
traps,
masks,
op,
cell_labels,
) )
return d return d
def extract_tp(
self,
tp: int,
tree: t.Optional[extraction_tree] = None,
tile_size: int = 117,
masks: t.Optional[t.List[np.ndarray]] = None,
cell_labels: t.Optional[t.List[int]] = None,
**kwargs,
) -> t.Dict[str, t.Dict[str, t.Dict[str, tuple]]]:
"""
Extract for an individual time point.
Parameters
----------
tp : int
Time point being analysed.
tree : dict
Nested dictionary indicating channels, reduction functions
and metrics to be used.
For example: {'general': {'None': ['area', 'volume', 'eccentricity']}}
tile_size : int
Size of the tile to be extracted.
masks : list of arrays
A list of masks per trap with each mask having dimensions
(ncells, tile_size, tile_size) and with one mask per cell.
cell_labels : dict
A dictionary with trap_ids as keys and cell_labels as values.
**kwargs : keyword arguments
Passed to extractor.reduce_extract.
Returns
-------
d: dict
Dictionary of the results with three levels of dictionaries.
The first level has channels as keys.
The second level has reduction metrics as keys.
The third level has cell or background metrics as keys and a
two-tuple as values.
The first tuple is the result of applying the metrics to a
particular cell or trap; the second tuple is either
(trap_id, cell_label) for a metric applied to a cell or a
trap_id for a metric applied to a trap.
An example is d["GFP"]["np_max"]["mean"][0], which gives a tuple
of the calculated mean GFP fluorescence for all cells.
"""
# dict of information from extraction tree
tree_bits = self.make_tree_bits(tree)
# create a Cells object to extract information from the h5 file
cells = Cells(self.local)
# find the cell labels as dict with trap_ids as keys
cell_labels = self.get_cell_labels(tp, cell_labels, cells)
# get masks one per cell per trap
masks = self.get_masks(tp, masks, cells)
# find image data at the time point
# stored as an array arranged as (traps, channels, 1, Z, X, Y)
tiles = self.get_tiles(
tp, tile_shape=tile_size, channels=tree_bits["tree_channels"]
)
# generate boolean masks for background for each trap
bgs = self.get_background_masks(masks, tile_size)
# perform extraction
res_one, self.img_bgsub = self.extract_one_channel(
tree_bits, cell_labels, tiles, masks, bgs, **kwargs
)
res_multiple = self.extract_multiple_channels(
tree_bits, cell_labels, tiles, masks, **kwargs
)
res = {**res_one, **res_multiple}
return res
def get_imgs(self, channel: t.Optional[str], tiles, channels=None): def get_imgs(self, channel: t.Optional[str], tiles, channels=None):
""" """
Return image from a correct source, either raw or bgsub. Return image from a correct source, either raw or bgsub.
...@@ -555,14 +669,16 @@ class Extractor(StepABC): ...@@ -555,14 +669,16 @@ class Extractor(StepABC):
channel: str channel: str
Name of channel to get. Name of channel to get.
tiles: ndarray tiles: ndarray
An array of the image data having dimensions of (tile_id, channel, tp, tile_size, tile_size, n_zstacks). An array of the image data having dimensions of
(tile_id, channel, tp, tile_size, tile_size, n_zstacks).
channels: list of str (optional) channels: list of str (optional)
t.List of available channels. t.List of available channels.
Returns Returns
------- -------
img: ndarray img: ndarray
An array of image data with dimensions (no tiles, X, Y, no Z channels) An array of image data with dimensions
(no tiles, X, Y, no Z channels)
""" """
if channels is None: if channels is None:
channels = (*self.params.tree,) channels = (*self.params.tree,)
...@@ -579,7 +695,9 @@ class Extractor(StepABC): ...@@ -579,7 +695,9 @@ class Extractor(StepABC):
**kwargs, **kwargs,
) -> dict: ) -> dict:
""" """
Wrapper to add compatibility with other steps of the pipeline. Run extraction for one position and for the specified time points.
Save the results to a h5 file.
Parameters Parameters
---------- ----------
...@@ -597,7 +715,9 @@ class Extractor(StepABC): ...@@ -597,7 +715,9 @@ class Extractor(StepABC):
Returns Returns
------- -------
d: dict d: dict
A dict of the extracted data with a concatenated string of channel, reduction metric, and cell metric as keys and pd.Series of the extracted data as values. A dict of the extracted data for one position with a concatenated
string of channel, reduction metric, and cell metric as keys and
pd.DataFrame of the extracted data for all time points as values.
""" """
if tree is None: if tree is None:
tree = self.params.tree tree = self.params.tree
...@@ -628,12 +748,12 @@ class Extractor(StepABC): ...@@ -628,12 +748,12 @@ class Extractor(StepABC):
d[k].index.names = idx d[k].index.names = idx
# save # save
if save: if save:
self.save_to_hdf(d) self.save_to_h5(d)
return d return d
def save_to_hdf(self, dict_series, path=None): def save_to_h5(self, dict_series, path=None):
""" """
Save the extracted data to the h5 file. Save the extracted data for one position to the h5 file.
Parameters Parameters
---------- ----------
...@@ -672,14 +792,17 @@ def flatten_nesteddict( ...@@ -672,14 +792,17 @@ def flatten_nesteddict(
nest: dict of dicts nest: dict of dicts
Contains the nested results of extraction. Contains the nested results of extraction.
to: str (optional) to: str (optional)
Specifies the format of the output, either pd.Series (default) or a list Specifies the format of the output, either pd.Series (default)
or a list
tp: int tp: int
Time point used to name the pd.Series Time point used to name the pd.Series
Returns Returns
------- -------
d: dict d: dict
A dict with a concatenated string of channel, reduction metric, and cell metric as keys and either a pd.Series or a list of the corresponding extracted data as values. A dict with a concatenated string of channel, reduction metric,
and cell metric as keys and either a pd.Series or a list of the
corresponding extracted data as values.
""" """
d = {} d = {}
for k0, v0 in nest.items(): for k0, v0 in nest.items():
...@@ -689,14 +812,3 @@ def flatten_nesteddict( ...@@ -689,14 +812,3 @@ def flatten_nesteddict(
pd.Series(*v2, name=tp) if to == "series" else v2 pd.Series(*v2, name=tp) if to == "series" else v2
) )
return d return d
class hollowExtractor(Extractor):
"""
Extractor that only cares about receiving images and masks.
Used for testing.
"""
def __init__(self, parameters):
self.params = parameters
src/extraction/core/functions/cell.py
View file @ b95e79fc
...@@ -90,9 +90,10 @@ def max2p5pc(cell_mask, trap_image) -> float: ...@@ -90,9 +90,10 @@ def max2p5pc(cell_mask, trap_image) -> float:
return np.mean(top_values) return np.mean(top_values)
def max5px(cell_mask, trap_image) -> float: def max5px_median(cell_mask, trap_image) -> float:
""" """
Find the mean of the five brightest pixels in the cell. Find the mean of the five brightest pixels in the cell divided by the
median of all pixels.
Parameters Parameters
---------- ----------
...@@ -105,7 +106,7 @@ def max5px(cell_mask, trap_image) -> float: ...@@ -105,7 +106,7 @@ def max5px(cell_mask, trap_image) -> float:
top_values = bn.partition(pixels, len(pixels) - 5)[-5:] top_values = bn.partition(pixels, len(pixels) - 5)[-5:]
# find mean of five brightest pixels # find mean of five brightest pixels
max5px = np.mean(top_values) max5px = np.mean(top_values)
return max5px return max5px / np.median(pixels)
def std(cell_mask, trap_image): def std(cell_mask, trap_image):
...@@ -193,3 +194,50 @@ def min_maj_approximation(cell_mask) -> t.Tuple[int]: ...@@ -193,3 +194,50 @@ def min_maj_approximation(cell_mask) -> t.Tuple[int]:
# + distance from the center of cone top to edge of cone top # + distance from the center of cone top to edge of cone top
maj_ax = np.round(np.max(dn) + np.sum(cone_top) / 2) maj_ax = np.round(np.max(dn) + np.sum(cone_top) / 2)
return min_ax, maj_ax return min_ax, maj_ax
def moment_of_inertia(cell_mask, trap_image):
"""
Find moment of inertia - a measure of homogeneity.
From iopscience.iop.org/article/10.1088/1742-6596/1962/1/012028
which cites ieeexplore.ieee.org/document/1057692.
"""
# set pixels not in cell to zero
trap_image[~cell_mask] = 0
x = trap_image
if np.any(x):
# x-axis : column=x-axis
columnvec = np.arange(1, x.shape[1] + 1, 1)[:, None].T
# y-axis : row=y-axis
rowvec = np.arange(1, x.shape[0] + 1, 1)[:, None]
# find raw moments
M00 = np.sum(x)
M10 = np.sum(np.multiply(x, columnvec))
M01 = np.sum(np.multiply(x, rowvec))
# find centroid
Xm = M10 / M00
Ym = M01 / M00
# find central moments
Mu00 = M00
Mu20 = np.sum(np.multiply(x, (columnvec - Xm) ** 2))
Mu02 = np.sum(np.multiply(x, (rowvec - Ym) ** 2))
# find invariants
Eta20 = Mu20 / Mu00 ** (1 + (2 + 0) / 2)
Eta02 = Mu02 / Mu00 ** (1 + (0 + 2) / 2)
# find moments of inertia
moi = Eta20 + Eta02
return moi
else:
return np.nan
def ratio(cell_mask, trap_image):
"""Find the median ratio between two fluorescence channels."""
if trap_image.ndim == 3 and trap_image.shape[-1] == 2:
fl_1 = trap_image[..., 0][cell_mask]
fl_2 = trap_image[..., 1][cell_mask]
div = np.median(fl_1 / fl_2)
else:
div = np.nan
return div
src/extraction/core/functions/custom/localisation.py
View file @ b95e79fc
""" How to do the nuc Est Conv from MATLAB """
How to do the nuc Est Conv from MATLAB
Based on the code in MattSegCode/Matt Seg Based on the code in MattSegCode/Matt Seg
GUI/@timelapseTraps/extractCellDataStacksParfor.m GUI/@timelapseTraps/extractCellDataStacksParfor.m
Especially lines 342 to 399. Especially lines 342 to 399.
This part only replicates the method to get the nuc_est_conv values This part only replicates the method to get the nuc_est_conv values
""" """
import typing as t import typing as t
......
src/extraction/core/functions/defaults.py
View file @ b95e79fc
# File with defaults for ease of use # File with defaults for ease of use
import re
import typing as t import typing as t
from pathlib import Path from pathlib import Path
import h5py import h5py
# should we move these functions here?
from aliby.tile.tiler import find_channel_name from aliby.tile.tiler import find_channel_name
...@@ -59,6 +58,7 @@ def exparams_from_meta( ...@@ -59,6 +58,7 @@ def exparams_from_meta(
for ch in extant_fluorescence_ch: for ch in extant_fluorescence_ch:
base["tree"][ch] = default_reduction_metrics base["tree"][ch] = default_reduction_metrics
base["sub_bg"] = extant_fluorescence_ch base["sub_bg"] = extant_fluorescence_ch
# additional extraction defaults if the channels are available # additional extraction defaults if the channels are available
if "ph" in extras: if "ph" in extras:
# SWAINLAB specific names # SWAINLAB specific names
......