diff --git a/src/pixelator/pna/anndata.py b/src/pixelator/pna/anndata.py index 94815b9f7..5914c56b5 100644 --- a/src/pixelator/pna/anndata.py +++ b/src/pixelator/pna/anndata.py @@ -41,7 +41,9 @@ def add_panel_information(adata: AnnData, panel: PNAAntibodyPanel) -> AnnData: def pna_edgelist_to_anndata( - pixel_connection: duckdb.DuckDBPyConnection, panel: PNAAntibodyPanel + pixel_connection: duckdb.DuckDBPyConnection, + panel: PNAAntibodyPanel, + markers: list[str] | None = None, ) -> AnnData: """Build an AnnData object from a DuckDB connection to a pixel file and a panel object. @@ -49,6 +51,9 @@ def pna_edgelist_to_anndata( pixel_connection: A DuckDB connection to a pixel file. The connection must contain an 'edgelist' table with the required columns (e.g., component, marker_1, marker_2, umi1, umi2, read_count). panel: The antibody panel object containing marker metadata. + markers: Marker ids to use as the count-matrix columns. Defaults to + every marker on ``panel``. Sample calling passes the markers left + after hashing clones collapse to their base name. Returns: An AnnData object with counts and panel information. @@ -73,10 +78,11 @@ def pna_edgelist_to_anndata( .tolist() ) + marker_list = list(panel.markers if markers is None else markers) n_components = len(components) - n_markers = len(panel.markers) + n_markers = len(marker_list) component_to_idx = {c: i for i, c in enumerate(components)} - marker_to_idx = {m: i for i, m in enumerate(panel.markers)} + marker_to_idx = {m: i for i, m in enumerate(marker_list)} X = np.zeros((n_components, n_markers), dtype=np.uint32) n_umi1_arr = np.zeros(n_components, dtype=np.uint64) @@ -152,7 +158,7 @@ def pna_edgelist_to_anndata( node_counts_df = pd.DataFrame( X, index=component_index, - columns=pd.Index(panel.markers, name="marker_id"), + columns=pd.Index(marker_list, name="marker_id"), ) logger.debug("Computing component metrics.") @@ -182,7 +188,7 @@ def pna_edgelist_to_anndata( logger.debug("Computing antibody metrics.") antibody_metrics_df = calculate_antibody_metrics(counts_df=node_counts_df) - antibody_metrics_df = antibody_metrics_df.reindex(index=panel.markers, fill_value=0) + antibody_metrics_df = antibody_metrics_df.reindex(index=marker_list, fill_value=0) antibody_metrics_df.index.name = "marker_id" # Do a dtype conversion of the columns here since AnnData cannot handle # a pyarrow arrays. @@ -201,7 +207,7 @@ def pna_edgelist_to_anndata( adata = add_panel_information(adata, panel) total_marker_counts = node_counts_df.sum(axis=1) - isotype_markers = adata.var[adata.var["control"]].index + isotype_markers = panel.df.index[panel.df["control"]] isotype_counts = node_counts_df[isotype_markers].sum(axis=1) adata.obs["isotype_fraction"] = isotype_counts / total_marker_counts diff --git a/src/pixelator/pna/cli/sample_calling.py b/src/pixelator/pna/cli/sample_calling.py index df24c2c7f..d16039956 100644 --- a/src/pixelator/pna/cli/sample_calling.py +++ b/src/pixelator/pna/cli/sample_calling.py @@ -118,9 +118,13 @@ def sample_calling_cli( pool_metadata=PxlFile(Path(input_pxl_file)).metadata(), ) return - hashing_antibodies_in_panel = set( - panel_info.df[panel_info.df["sample_hashing"] == "yes"].index.to_list() - ) + if "sample_hashing" not in panel_info.df.columns: + raise ValueError( + "Sample calling requires a sample_hashing column on the panel " + "so hashing markers can be identified. This panel has no " + "sample_hashing column." + ) + hashing_antibodies_in_panel = panel_info.hashing_marker_ids samplesheet_df = pl.read_csv(samplesheet) _reject_reserved_samplesheet_names( samplesheet_df["sample"].to_list(), diff --git a/src/pixelator/pna/config/panel.py b/src/pixelator/pna/config/panel.py deleted file mode 100644 index cf89c7a1e..000000000 --- a/src/pixelator/pna/config/panel.py +++ /dev/null @@ -1,623 +0,0 @@ -"""Marker panel management for different PNA assays. - -Copyright © 2022 Pixelgen Technologies AB. -""" - -from __future__ import annotations - -import warnings -from functools import cached_property -from pathlib import Path -from typing import TYPE_CHECKING, List, Optional, Set - -from anndata import AnnData - -try: - from typing import Self -except ImportError: - from typing_extensions import Self - -import re - -import pandas as pd -import polars as pl - -from pixelator.common.config.panel import ( - AntibodyPanelMetadata, - parse_panel_header_metadata, -) -from pixelator.common.types import PathType -from pixelator.common.utils import logger - -if TYPE_CHECKING: - from pixelator.pna.config.config_class import PNAConfig - from pixelator.pna.pixeldataset.dataset import PNAPixelDataset - - -class PNAAntibodyPanel: - """Class representing a PNA antibody panel.""" - - # required columns - _INDEX_COLUMN = "marker_id" - _INDEX_COLUMN_TYPE = str - _REQUIRED_COLUMNS = { - "control": bool, - "sequence_1": str, - "sequence_2": str, - } - - # and these should have unique values - _UNIQUE_COLUMNS = ["sequence_1", "sequence_2"] - - def __init__( - self, - df: pd.DataFrame, - metadata: AntibodyPanelMetadata, - file_name: Optional[str] = None, - filepath: Optional[PathType] = None, - ) -> None: - """Load a panel from a dataframe and metadata. - - Args: - df: The dataframe containing the panel information. - metadata: The metadata for the panel. - file_name: The optional basename of the file from which the panel is loaded. - filepath: The optional full path of the file from which the panel is loaded. - - Returns: - None - Raises: - AssertionError: exception if panel file is missing, invalid or with incorrect format - """ - self._filename = file_name - self._filepath: Optional[Path] = Path(filepath).resolve() if filepath else None - self.metadata = metadata - self._df = df - - # validate the panel - errors = self.validate_antibody_panel(df) - if len(errors) > 0: - msg_str = "\n".join(errors) - raise AssertionError( - f"The following errors were found validating the panel: {msg_str}" - ) - - @classmethod - def from_csv(cls, filename: PathType) -> Self: - """Create an AntibodyPanel from a csv panel file. - - Args: - filename: The path to the panel file. - - Returns: - The AntibodyPanel object. (AntibodyPanel) - - Raises: - AssertionError: exception if panel file is missing, - """ - panel_file = Path(filename) - - if not panel_file.is_file() or panel_file.suffix != ".csv": - raise AssertionError( - f"Panel file {filename} not found or has an incorrect format" - ) - - logger.debug("Creating Antibody panel from file %s", filename) - - df = cls._parse_panel(panel_file) - metadata = cls._parse_header(panel_file) - - logger.debug("Antibody panel from file %s created", filename) - - return cls(df, metadata, file_name=panel_file.name, filepath=panel_file) - - @classmethod - def from_pxl_dataset( - cls, pxl_data: PNAPixelDataset, file_name: Optional[str] = None - ) -> Self: - """Create an AntibodyPanel from a pxl dataset. - - Args: - pxl_data: A PNAPixelDataset object. - file_name: The optional name of the file from which the pxl dataset was loaded. - - Returns: - The AntibodyPanel object. (AntibodyPanel) - - Raises: - KeyError: exception if panel information is missing in the pxl dataset, - """ - logger.debug("Creating Antibody panel from PNAPixelDataset object") - adata = pxl_data.adata() - panel = cls.from_adata(adata, file_name=file_name) - logger.debug("Antibody panel from PNAPixelDataset created") - return panel - - @classmethod - def from_adata(cls, adata: AnnData, file_name: Optional[str] = None) -> Self: - """Create an AntibodyPanel from an AnnData object. - - Args: - adata: An AnnData object containing panel information. - file_name: The optional name of the file from which the AnnData object was loaded. - - Returns: - The AntibodyPanel object. (AntibodyPanel) - - Raises: - KeyError: exception if panel information is missing in the AnnData object. - """ - logger.debug("Creating Antibody panel from AnnData object") - try: - panel_metadata = adata.uns["panel_metadata"] - except KeyError as err: - logger.error( # pylint: disable=logging-not-lazy - f"The provided AnnData object does not contain {err}. " - + "Please, regenerate your data with the most recent version of pixelator." - ) - raise - panel_columns = panel_metadata.get("panel_columns") - if not panel_columns: - raise KeyError( - "The provided AnnData object does not contain panel columns information in the metadata. " - + "Please, regenerate your data with the most recent version of pixelator." - ) - df = adata.var[panel_columns] - metadata = AntibodyPanelMetadata.model_validate(panel_metadata) - - logger.debug("Antibody panel from AnnData object created") - return cls(df, metadata, file_name=file_name) - - @property - def name(self) -> str: - """Panel name from metadata. - - Returns: - The panel name. - """ - return self.metadata.name - - @property - def product(self) -> Optional[str]: - """Product identifier from metadata, if present. - - Returns: - Product name, or None when not provided in panel metadata. - """ - return self.metadata.product - - @property - def version(self) -> str: - """Panel version from metadata. - - Returns: - Semantic version string for this panel. - """ - return self.metadata.version - - @property - def description(self) -> Optional[str]: - """Return the panel file description.""" - return self.metadata.description - - @property - def aliases(self) -> list[str]: - """Return the (optional) list of panel file aliases.""" - return self.metadata.aliases - - @property - def archived(self) -> Optional[bool]: - """Return whether the panel is marked as archived.""" - return self.metadata.archived - - @classmethod - def _parse_header(cls, file: Path) -> AntibodyPanelMetadata: - """Parse front-matter YAML metadata from a panel file. - - Args: - file: Panel CSV file whose leading comment block contains YAML metadata. - - Returns: - Parsed panel metadata. - - Raises: - ValueError: If no metadata header is present in the file. - """ - return parse_panel_header_metadata(file) - - @classmethod - def _parse_panel(cls, panel_file: Path) -> pd.DataFrame: - panel = pd.read_csv(str(panel_file), comment="#", index_col="marker_id").fillna( - "" - ) - - panel["control"] = panel["control"].map(lambda s: s.lower() == "yes") - - return panel.copy() - - @cached_property - def markers_control(self) -> List[str]: - """Return a list of marker control (names).""" - return list(self._df[self._df["control"]].index) - - @cached_property - def markers(self) -> List[str]: - """Return the list of unique markers in the panel.""" - return list(self._df.index.unique()) - - @property - def df(self) -> pd.DataFrame: - """Return the panel dataframe.""" - return self._df - - @property - def filename(self) -> Optional[str]: - """Return the filename of the marker panel.""" - return self._filename - - @property - def filepath(self) -> Optional[Path]: - """Return the full path of the marker panel file, if any.""" - return self._filepath - - @cached_property - def size(self) -> int: - """Return the size of the marker panel.""" - return self._df.shape[0] - - @staticmethod - def _validate_sequences(panel_df, sequence_col): - errors = [] - sequences = panel_df[sequence_col] - ref_length = len(sequences.iloc[0]) - if not sequences.apply(lambda x: len(x) == ref_length).all(): - errors.append(f"All {sequence_col} values must have the same length.") - - if not sequences.str.match("^[ATCG]*$").all(): - errors.append( - f"All {sequence_col} values must only contain ATCG characters. Offending values: " - f"{sequences[~sequences.str.match('^[ATCG]*$')].tolist()}" - ) - - return errors - - @staticmethod - def _validate_marker_names(panel_df): - errors = [] - if any(panel_df.index.str.contains("_")): - # Markers should not contain underscores since this messes - # things up with Seurat on the R side - errors.append( - "The marker_id column should not contain underscores. " - "Please use dashes instead. Offending values: " - f"{panel_df.index[panel_df.index.str.contains('_')]}" - ) - if any(panel_df.index.str.contains(r"\s")): - # Markers should not contain white-spaces since this causes - # issues in the demultiplexing step (and other places that - # might assume that marker names are single tokens) - problematic_lines = panel_df.index[panel_df.index.str.contains(r"\s")] - errors.append( - "The marker_id column should not contain white-spaces. " - "Please use dashes instead or remove the white-spaces. Offending values: " - f"{problematic_lines}" - ) - return errors - - @classmethod - def validate_antibody_panel( - cls, panel_df: pd.DataFrame, validate_types: bool = True - ) -> list[str]: - """Validate antibody panel schema and content. - - Args: - panel_df: Dataframe containing panel markers and sequences. - validate_types: If True, validate dataframe column types. - - Returns: - A list of validation error messages. Empty means valid input. - """ - errors = [] - - # some basic sanity check on the panel size and columns - if not set(cls._REQUIRED_COLUMNS).issubset(set(panel_df.columns)): - missing_columns = set(cls._REQUIRED_COLUMNS) - set(panel_df.columns) - errors.append(f"Panel has missing required columns: {missing_columns}") - return errors - - if validate_types: - panel_pl_df = pl.from_pandas(panel_df, include_index=True) - for col, expected_type in ( - cls._REQUIRED_COLUMNS | {cls._INDEX_COLUMN: cls._INDEX_COLUMN_TYPE} - ).items(): - found_type = panel_pl_df[col].dtype.to_python() - if not found_type == expected_type: - errors.append( - f"Column {col} has incorrect type. Expected {expected_type}, got {found_type}" - ) - - if panel_df.shape[0] == 0: - errors.append("Panel file is empty") - return errors - - # sanity check on the unique columns - for col in cls._UNIQUE_COLUMNS: - if not len(panel_df[col].unique()) == len(panel_df[col]): - errors.append(f"All values in column: {col} were not unique") - - if panel_df.index.name != cls._INDEX_COLUMN: - errors.append(f"`{cls._INDEX_COLUMN}` is missing or is not set as index") - return errors - - errors += cls._validate_marker_names(panel_df) - - if panel_df["control"].dtype != bool: - errors.append("`control` column is not boolean") - - # Check UniProt IDs format conforming to the UniProt naming convention. Empty IDs are allowed. - if "uniprot_id" in panel_df.columns: - # Pattern for valid UniProt IDs - pattern = r"^[OPQ][0-9][A-Z0-9]{3}[0-9]|[A-NR-Z][0-9]([A-Z][A-Z0-9]{2}[0-9]){1,2}|$" - - def check_id(id_str): - """Check id. - - Args: - id_str: id str. - """ - return all( - bool(re.match(pattern, id_)) for id_ in str(id_str).split(";") - ) - - bad_ids = panel_df[~panel_df["uniprot_id"].apply(check_id)]["uniprot_id"] - - if len(bad_ids) > 0: - errors.append( - "Invalid UniProt IDs found." - "Please conform to the naming convention or remove the following IDs:" - f"{bad_ids.tolist()}" - ) - - errors += cls._validate_sequences(panel_df, "sequence_1") - errors += cls._validate_sequences(panel_df, "sequence_2") - - return errors - - def to_polars(self) -> pl.DataFrame: - """Convert the panel to a Polars DataFrame.""" - return pl.from_pandas(self.df, include_index=True) - - def __eq__(self, other: object) -> bool: - """Check if two panels are equal based on their dataframes and metadata. - - Args: - other: Panel to compare for equality. - """ - if not isinstance(other, PNAAntibodyPanel): - raise ValueError("Can only compare with another PNAAntibodyPanel") - return self.df.equals(other.df) and self.metadata == other.metadata - - -def load_antibody_panel(config: PNAConfig, panel: PathType) -> PNAAntibodyPanel: - """Load an antibody panel from a file or from the config file. - - Args: - config: the config object - panel: the path to the panel file or the name of the panel in the config file - Returns: - the antibody panel (PNAAntibodyPanel) - """ - panel_str = str(panel) - panel_from_config = config.get_panel(panel_str) - - if panel_from_config is not None: - logger.info("Found panel in config file: %s", panel_from_config.name) - return panel_from_config - - panel_obj = PNAAntibodyPanel.from_csv(panel) - return panel_obj - - -class PNAAntibodyPanelDiff: - """Class representing the differences between two PNAAntibodyPanel objects.""" - - join_on_columns: list[str] = ["sequence_1", "sequence_2"] - - def __init__(self, panel_1: PNAAntibodyPanel, panel_2: PNAAntibodyPanel) -> None: - """Initialize the PNAAntibodyPanelDiff object. - - Args: - panel_1: The first panel to compare. - panel_2: The second panel to compare. - """ - self.panel_1 = panel_1 - self.panel_2 = panel_2 - - logger.debug( - "Comparing panels %s v%s and %s v%s", - panel_1.name, - panel_1.version, - panel_2.name, - panel_2.version, - ) - - self.joined = self.panel_1.to_polars().join( - self.panel_2.to_polars(), - on=self.join_on_columns, - how="full", - suffix="_panel_2", - ) - - self._identical_columns: List[str] | None = None - self._changed_columns: List[str] | None = None - self._removed_columns: Set[str] | None = None - self._added_columns: Set[str] | None = None - - @property - def col_names_in_both_panels(self) -> List[str]: - """Return a list of column names that are present in both panels.""" - return list( - set(self.panel_1.to_polars().columns).intersection( - set(self.panel_2.to_polars().columns) - ) - ) - - @property - def identical_columns(self) -> List[str]: - """Return a list of columns that are identical between the two panels.""" - return [ - col_name - for col_name in self.col_names_in_both_panels - if self.joined[col_name] - .eq_missing(self.joined[col_name + "_panel_2"]) - .all() - ] - - @cached_property - def changed_columns(self) -> List[str]: - """Return a list of columns that are different between the two panels.""" - changed_columns = [ - col_name - for col_name in set(self.col_names_in_both_panels).difference( - set(self.join_on_columns) - ) - if not self.joined[col_name] - .eq_missing(self.joined[col_name + "_panel_2"]) - .all() - ] - for col_name in changed_columns: - diff_count = self.joined.filter( - pl.col(col_name).ne_missing(pl.col(col_name + "_panel_2")) - ).shape[0] - logger.debug( - "Column %s is different between the two panels %s and %s (%d differing entries).", - col_name, - self.panel_1.name, - self.panel_2.name, - diff_count, - ) - return changed_columns - - @cached_property - def removed_columns(self) -> List[str]: - """Return a list of columns that are present in panel 1 but not in panel 2.""" - removed_columns = set(self.panel_1.to_polars().columns).difference( - set(self.panel_2.to_polars().columns) - ) - for col_name in removed_columns: - logger.debug( - "Column %s is present in panel %s but not in panel %s.", - col_name, - self.panel_1.name, - self.panel_2.name, - ) - return sorted(removed_columns) - - @cached_property - def added_columns(self) -> List[str]: - """Return a list of columns that are present in panel 2 but not in panel 1.""" - added_columns = set(self.panel_2.to_polars().columns).difference( - set(self.panel_1.to_polars().columns) - ) - for col_name in added_columns: - logger.debug( - "Column %s is present in panel %s but not in panel %s.", - col_name, - self.panel_2.name, - self.panel_1.name, - ) - return sorted(added_columns) - - @property - def added_clones(self) -> pl.DataFrame: - """Return a dataframe with the clones that are present in panel 2 but not in panel 1.""" - return ( - self.joined.filter( - pl.any_horizontal( - pl.col(col_name).is_null() - & pl.col(col_name + "_panel_2").is_not_null() - for col_name in self.join_on_columns - ) - ) - .drop([col_name for col_name in self.panel_1.to_polars().columns]) - .rename( - { - col_name + "_panel_2": col_name - for col_name in self.panel_2.to_polars().columns - if col_name + "_panel_2" in self.joined.columns - } - ) - ) - - @property - def removed_clones(self) -> pl.DataFrame: - """Return a dataframe with the clones that are present in panel 1 but not in panel 2.""" - return self.joined.filter( - pl.any_horizontal( - pl.col(col_name).is_not_null() & pl.col(col_name + "_panel_2").is_null() - for col_name in self.join_on_columns - ) - ).drop( - [ - col_name + "_panel_2" - if col_name in self.joined.columns - and col_name not in self.added_columns - else col_name - for col_name in self.panel_2.to_polars().columns - ] - ) - - def upgrade_adata(self, adata: AnnData) -> AnnData: - """Upgrade an AnnData object with the changes between the two panels. - - Args: - adata: An AnnData object containing panel information. - """ - adata_panel = PNAAntibodyPanel.from_adata(adata) - if self.panel_1 != adata_panel: - raise ValueError( - "The provided AnnData object does not match the panel. Cannot upgrade." - f"Expected panel {self.panel_2.name} v{self.panel_2.version}, but got panel {adata_panel.name} v{adata_panel.version}." - ) - - non_panel_columns = adata.var.copy()[ - [ - col - for col in adata.var.columns - if col not in adata.uns["panel_metadata"]["panel_columns"] - ] - + self.join_on_columns - ] - adata.var = ( - self.joined.select( - list( - set( - self.join_on_columns - + self.identical_columns - + [f"{col}_panel_2" for col in self.changed_columns] - + self.added_columns - ) - ) - ) - .rename({f"{col}_panel_2": col for col in self.changed_columns}) - # keep order and append new to the end - .select( - ["marker_id"] # index not in panel_metadata panel_columns below - + adata.uns["panel_metadata"]["panel_columns"] - + self.added_columns - ) - .to_pandas() - .set_index("marker_id") - ) - if adata.var.shape[0] != non_panel_columns.shape[0]: - raise ValueError( - "Row count mismatch in automatic patch panel patch version bump." - ) - adata.var = adata.var.join( - non_panel_columns.set_index(self.join_on_columns), - how="outer", - on=self.join_on_columns, - ) - adata.uns["panel_metadata"] = self.panel_2.metadata.model_dump() - adata.uns["panel_metadata"]["panel_columns"] = self.panel_2.df.columns.tolist() - return adata diff --git a/src/pixelator/pna/config/panel/__init__.py b/src/pixelator/pna/config/panel/__init__.py new file mode 100644 index 000000000..feb9b2277 --- /dev/null +++ b/src/pixelator/pna/config/panel/__init__.py @@ -0,0 +1,24 @@ +"""Marker panel management for different PNA assays. + +Copyright © 2022 Pixelgen Technologies AB. +""" + +from pixelator.pna.config.panel.antibody_panel import ( + PNAAntibodyPanel, + load_antibody_panel, +) +from pixelator.pna.config.panel.diff import PNAAntibodyPanelDiff +from pixelator.pna.config.panel.hashing import ( + collapsed_hashing_marker_id, + sample_hashing_mask, + split_hashing_marker_id, +) + +__all__ = [ + "PNAAntibodyPanel", + "PNAAntibodyPanelDiff", + "collapsed_hashing_marker_id", + "load_antibody_panel", + "sample_hashing_mask", + "split_hashing_marker_id", +] diff --git a/src/pixelator/pna/config/panel/antibody_panel.py b/src/pixelator/pna/config/panel/antibody_panel.py new file mode 100644 index 000000000..bf96e5b78 --- /dev/null +++ b/src/pixelator/pna/config/panel/antibody_panel.py @@ -0,0 +1,318 @@ +"""PNA antibody panel loading. + +Copyright © 2022 Pixelgen Technologies AB. +""" + +from __future__ import annotations + +from functools import cached_property +from pathlib import Path +from typing import TYPE_CHECKING, List, Optional + +from anndata import AnnData + +try: + from typing import Self +except ImportError: + from typing_extensions import Self + +import pandas as pd +import polars as pl + +from pixelator.common.config.panel import ( + AntibodyPanelMetadata, + parse_panel_header_metadata, +) +from pixelator.common.types import PathType +from pixelator.common.utils import logger +from pixelator.pna.config.panel.hashing import _hashing_marker_ids +from pixelator.pna.config.panel.validation import ( + _INDEX_COLUMN, + _INDEX_COLUMN_TYPE, + _REQUIRED_COLUMNS, + _UNIQUE_COLUMNS, + validate_antibody_panel, +) + +if TYPE_CHECKING: + from pixelator.pna.config.config_class import PNAConfig + from pixelator.pna.pixeldataset.dataset import PNAPixelDataset + + +class PNAAntibodyPanel: + """Class representing a PNA antibody panel.""" + + _INDEX_COLUMN = _INDEX_COLUMN + _INDEX_COLUMN_TYPE = _INDEX_COLUMN_TYPE + _REQUIRED_COLUMNS = _REQUIRED_COLUMNS + _UNIQUE_COLUMNS = _UNIQUE_COLUMNS + + def __init__( + self, + df: pd.DataFrame, + metadata: AntibodyPanelMetadata, + file_name: Optional[str] = None, + filepath: Optional[PathType] = None, + ) -> None: + """Load a panel from a dataframe and metadata. + + Args: + df: The dataframe containing the panel information. + metadata: The metadata for the panel. + file_name: The optional basename of the file from which the panel is loaded. + filepath: The optional full path of the file from which the panel is loaded. + + Returns: + None + Raises: + AssertionError: exception if panel file is missing, invalid or with incorrect format + """ + self._filename = file_name + self._filepath: Optional[Path] = Path(filepath).resolve() if filepath else None + self.metadata = metadata + self._df = df + + # validate the panel + errors = self.validate_antibody_panel(df) + if len(errors) > 0: + msg_str = "\n".join(errors) + raise AssertionError( + f"The following errors were found validating the panel: {msg_str}" + ) + + @classmethod + def from_csv(cls, filename: PathType) -> Self: + """Create an AntibodyPanel from a csv panel file. + + Args: + filename: The path to the panel file. + + Returns: + The AntibodyPanel object. (AntibodyPanel) + + Raises: + AssertionError: exception if panel file is missing, + """ + panel_file = Path(filename) + + if not panel_file.is_file() or panel_file.suffix != ".csv": + raise AssertionError( + f"Panel file {filename} not found or has an incorrect format" + ) + + logger.debug("Creating Antibody panel from file %s", filename) + + df = cls._parse_panel(panel_file) + metadata = cls._parse_header(panel_file) + + logger.debug("Antibody panel from file %s created", filename) + + return cls(df, metadata, file_name=panel_file.name, filepath=panel_file) + + @classmethod + def from_pxl_dataset( + cls, pxl_data: PNAPixelDataset, file_name: Optional[str] = None + ) -> Self: + """Create an AntibodyPanel from a pxl dataset. + + Args: + pxl_data: A PNAPixelDataset object. + file_name: The optional name of the file from which the pxl dataset was loaded. + + Returns: + The AntibodyPanel object. (AntibodyPanel) + + Raises: + KeyError: exception if panel information is missing in the pxl dataset, + """ + logger.debug("Creating Antibody panel from PNAPixelDataset object") + adata = pxl_data.adata() + panel = cls.from_adata(adata, file_name=file_name) + logger.debug("Antibody panel from PNAPixelDataset created") + return panel + + @classmethod + def from_adata(cls, adata: AnnData, file_name: Optional[str] = None) -> Self: + """Create an AntibodyPanel from an AnnData object. + + Args: + adata: An AnnData object containing panel information. + file_name: The optional name of the file from which the AnnData object was loaded. + + Returns: + The AntibodyPanel object. (AntibodyPanel) + + Raises: + KeyError: exception if panel information is missing in the AnnData object. + """ + logger.debug("Creating Antibody panel from AnnData object") + try: + panel_metadata = adata.uns["panel_metadata"] + except KeyError as err: + logger.error( # pylint: disable=logging-not-lazy + f"The provided AnnData object does not contain {err}. " + + "Please, regenerate your data with the most recent version of pixelator." + ) + raise + panel_columns = panel_metadata.get("panel_columns") + if not panel_columns: + raise KeyError( + "The provided AnnData object does not contain panel columns information in the metadata. " + + "Please, regenerate your data with the most recent version of pixelator." + ) + df = adata.var[panel_columns] + metadata = AntibodyPanelMetadata.model_validate(panel_metadata) + + logger.debug("Antibody panel from AnnData object created") + return cls(df, metadata, file_name=file_name) + + @property + def name(self) -> str: + """Panel name from metadata. + + Returns: + The panel name. + """ + return self.metadata.name + + @property + def product(self) -> Optional[str]: + """Product identifier from metadata, if present. + + Returns: + Product name, or None when not provided in panel metadata. + """ + return self.metadata.product + + @property + def version(self) -> str: + """Panel version from metadata. + + Returns: + Semantic version string for this panel. + """ + return self.metadata.version + + @property + def description(self) -> Optional[str]: + """Return the panel file description.""" + return self.metadata.description + + @property + def aliases(self) -> list[str]: + """Return the (optional) list of panel file aliases.""" + return self.metadata.aliases + + @property + def archived(self) -> Optional[bool]: + """Return whether the panel is marked as archived.""" + return self.metadata.archived + + @property + def hashing_marker_ids(self) -> set[str]: + """Return marker ids flagged by the ``sample_hashing`` column.""" + return _hashing_marker_ids(self.df) + + @classmethod + def _parse_header(cls, file: Path) -> AntibodyPanelMetadata: + """Parse front-matter YAML metadata from a panel file. + + Args: + file: Panel CSV file whose leading comment block contains YAML metadata. + + Returns: + Parsed panel metadata. + + Raises: + ValueError: If no metadata header is present in the file. + """ + return parse_panel_header_metadata(file) + + @classmethod + def _parse_panel(cls, panel_file: Path) -> pd.DataFrame: + panel = pd.read_csv(str(panel_file), comment="#", index_col="marker_id").fillna( + "" + ) + + panel["control"] = panel["control"].map(lambda s: s.lower() == "yes") + + return panel.copy() + + @cached_property + def markers_control(self) -> List[str]: + """Return a list of marker control (names).""" + return list(self._df[self._df["control"]].index) + + @cached_property + def markers(self) -> List[str]: + """Return the list of unique markers in the panel.""" + return list(self._df.index.unique()) + + @property + def df(self) -> pd.DataFrame: + """Return the panel dataframe.""" + return self._df + + @property + def filename(self) -> Optional[str]: + """Return the filename of the marker panel.""" + return self._filename + + @property + def filepath(self) -> Optional[Path]: + """Return the full path of the marker panel file, if any.""" + return self._filepath + + @cached_property + def size(self) -> int: + """Return the size of the marker panel.""" + return self._df.shape[0] + + @classmethod + def validate_antibody_panel( + cls, panel_df: pd.DataFrame, validate_types: bool = True + ) -> list[str]: + """Validate antibody panel schema and content. + + Args: + panel_df: Dataframe containing panel markers and sequences. + validate_types: If True, validate dataframe column types. + + Returns: + A list of validation error messages. Empty means valid input. + """ + return validate_antibody_panel(panel_df, validate_types=validate_types) + + def to_polars(self) -> pl.DataFrame: + """Convert the panel to a Polars DataFrame.""" + return pl.from_pandas(self.df, include_index=True) + + def __eq__(self, other: object) -> bool: + """Check if two panels are equal based on their dataframes and metadata. + + Args: + other: Panel to compare for equality. + """ + if not isinstance(other, PNAAntibodyPanel): + raise ValueError("Can only compare with another PNAAntibodyPanel") + return self.df.equals(other.df) and self.metadata == other.metadata + + +def load_antibody_panel(config: PNAConfig, panel: PathType) -> PNAAntibodyPanel: + """Load an antibody panel from a file or from the config file. + + Args: + config: the config object + panel: the path to the panel file or the name of the panel in the config file + Returns: + the antibody panel (PNAAntibodyPanel) + """ + panel_str = str(panel) + panel_from_config = config.get_panel(panel_str) + + if panel_from_config is not None: + logger.info("Found panel in config file: %s", panel_from_config.name) + return panel_from_config + + panel_obj = PNAAntibodyPanel.from_csv(panel) + return panel_obj diff --git a/src/pixelator/pna/config/panel/diff.py b/src/pixelator/pna/config/panel/diff.py new file mode 100644 index 000000000..055e8b848 --- /dev/null +++ b/src/pixelator/pna/config/panel/diff.py @@ -0,0 +1,220 @@ +"""Differences between two PNA antibody panels. + +Copyright © 2022 Pixelgen Technologies AB. +""" + +from __future__ import annotations + +from functools import cached_property +from typing import List, Set + +import polars as pl +from anndata import AnnData + +from pixelator.common.utils import logger +from pixelator.pna.config.panel.antibody_panel import PNAAntibodyPanel + + +class PNAAntibodyPanelDiff: + """Class representing the differences between two PNAAntibodyPanel objects.""" + + join_on_columns: list[str] = ["sequence_1", "sequence_2"] + + def __init__(self, panel_1: PNAAntibodyPanel, panel_2: PNAAntibodyPanel) -> None: + """Initialize the PNAAntibodyPanelDiff object. + + Args: + panel_1: The first panel to compare. + panel_2: The second panel to compare. + """ + self.panel_1 = panel_1 + self.panel_2 = panel_2 + + logger.debug( + "Comparing panels %s v%s and %s v%s", + panel_1.name, + panel_1.version, + panel_2.name, + panel_2.version, + ) + + self.joined = self.panel_1.to_polars().join( + self.panel_2.to_polars(), + on=self.join_on_columns, + how="full", + suffix="_panel_2", + ) + + self._identical_columns: List[str] | None = None + self._changed_columns: List[str] | None = None + self._removed_columns: Set[str] | None = None + self._added_columns: Set[str] | None = None + + @property + def col_names_in_both_panels(self) -> List[str]: + """Return a list of column names that are present in both panels.""" + return list( + set(self.panel_1.to_polars().columns).intersection( + set(self.panel_2.to_polars().columns) + ) + ) + + @property + def identical_columns(self) -> List[str]: + """Return a list of columns that are identical between the two panels.""" + return [ + col_name + for col_name in self.col_names_in_both_panels + if self.joined[col_name] + .eq_missing(self.joined[col_name + "_panel_2"]) + .all() + ] + + @cached_property + def changed_columns(self) -> List[str]: + """Return a list of columns that are different between the two panels.""" + changed_columns = [ + col_name + for col_name in set(self.col_names_in_both_panels).difference( + set(self.join_on_columns) + ) + if not self.joined[col_name] + .eq_missing(self.joined[col_name + "_panel_2"]) + .all() + ] + for col_name in changed_columns: + diff_count = self.joined.filter( + pl.col(col_name).ne_missing(pl.col(col_name + "_panel_2")) + ).shape[0] + logger.debug( + "Column %s is different between the two panels %s and %s (%d differing entries).", + col_name, + self.panel_1.name, + self.panel_2.name, + diff_count, + ) + return changed_columns + + @cached_property + def removed_columns(self) -> List[str]: + """Return a list of columns that are present in panel 1 but not in panel 2.""" + removed_columns = set(self.panel_1.to_polars().columns).difference( + set(self.panel_2.to_polars().columns) + ) + for col_name in removed_columns: + logger.debug( + "Column %s is present in panel %s but not in panel %s.", + col_name, + self.panel_1.name, + self.panel_2.name, + ) + return sorted(removed_columns) + + @cached_property + def added_columns(self) -> List[str]: + """Return a list of columns that are present in panel 2 but not in panel 1.""" + added_columns = set(self.panel_2.to_polars().columns).difference( + set(self.panel_1.to_polars().columns) + ) + for col_name in added_columns: + logger.debug( + "Column %s is present in panel %s but not in panel %s.", + col_name, + self.panel_2.name, + self.panel_1.name, + ) + return sorted(added_columns) + + @property + def added_clones(self) -> pl.DataFrame: + """Return a dataframe with the clones that are present in panel 2 but not in panel 1.""" + return ( + self.joined.filter( + pl.any_horizontal( + pl.col(col_name).is_null() + & pl.col(col_name + "_panel_2").is_not_null() + for col_name in self.join_on_columns + ) + ) + .drop([col_name for col_name in self.panel_1.to_polars().columns]) + .rename( + { + col_name + "_panel_2": col_name + for col_name in self.panel_2.to_polars().columns + if col_name + "_panel_2" in self.joined.columns + } + ) + ) + + @property + def removed_clones(self) -> pl.DataFrame: + """Return a dataframe with the clones that are present in panel 1 but not in panel 2.""" + return self.joined.filter( + pl.any_horizontal( + pl.col(col_name).is_not_null() & pl.col(col_name + "_panel_2").is_null() + for col_name in self.join_on_columns + ) + ).drop( + [ + col_name + "_panel_2" + if col_name in self.joined.columns + and col_name not in self.added_columns + else col_name + for col_name in self.panel_2.to_polars().columns + ] + ) + + def upgrade_adata(self, adata: AnnData) -> AnnData: + """Upgrade an AnnData object with the changes between the two panels. + + Args: + adata: An AnnData object containing panel information. + """ + adata_panel = PNAAntibodyPanel.from_adata(adata) + if self.panel_1 != adata_panel: + raise ValueError( + "The provided AnnData object does not match the panel. Cannot upgrade." + f"Expected panel {self.panel_2.name} v{self.panel_2.version}, but got panel {adata_panel.name} v{adata_panel.version}." + ) + + non_panel_columns = adata.var.copy()[ + [ + col + for col in adata.var.columns + if col not in adata.uns["panel_metadata"]["panel_columns"] + ] + + self.join_on_columns + ] + adata.var = ( + self.joined.select( + list( + set( + self.join_on_columns + + self.identical_columns + + [f"{col}_panel_2" for col in self.changed_columns] + + self.added_columns + ) + ) + ) + .rename({f"{col}_panel_2": col for col in self.changed_columns}) + # keep order and append new to the end + .select( + ["marker_id"] # index not in panel_metadata panel_columns below + + adata.uns["panel_metadata"]["panel_columns"] + + self.added_columns + ) + .to_pandas() + .set_index("marker_id") + ) + if adata.var.shape[0] != non_panel_columns.shape[0]: + raise ValueError( + "Row count mismatch in automatic patch panel patch version bump." + ) + adata.var = adata.var.join( + non_panel_columns.set_index(self.join_on_columns), + how="outer", + on=self.join_on_columns, + ) + adata.uns["panel_metadata"] = self.panel_2.metadata.model_dump() + adata.uns["panel_metadata"]["panel_columns"] = self.panel_2.df.columns.tolist() + return adata diff --git a/src/pixelator/pna/config/panel/hashing.py b/src/pixelator/pna/config/panel/hashing.py new file mode 100644 index 000000000..e9a1f1bdd --- /dev/null +++ b/src/pixelator/pna/config/panel/hashing.py @@ -0,0 +1,47 @@ +"""Helpers for sample-hashing marker ids. + +Copyright © 2022 Pixelgen Technologies AB. +""" + +from __future__ import annotations + +import re + +import pandas as pd + +# Trailing ``-`` is the hash group (``B2M-1`` → ``B2M``). The same +# pattern matches ordinary names such as ``PD-1``, so it is only applied to +# rows already flagged by ``sample_hashing``. +_HASHING_MARKER_ID_RE = re.compile(r"^(?P.+)-(?P\d+)$") + + +def sample_hashing_mask(sample_hashing: pd.Series) -> pd.Series: + """Return a boolean mask for values that flag a hashing marker.""" + if pd.api.types.is_bool_dtype(sample_hashing): + return sample_hashing.fillna(False).astype(bool) + if pd.api.types.is_numeric_dtype(sample_hashing): + return sample_hashing.fillna(0).astype(bool) + normalized = sample_hashing.astype(str).str.strip().str.lower() + return normalized.isin(["yes", "true"]) + + +def split_hashing_marker_id(marker_id: str) -> tuple[str, str] | None: + """Return ``(base, index)`` for a hashing id such as ``B2M-1``.""" + match = _HASHING_MARKER_ID_RE.fullmatch(str(marker_id)) + if match is None: + return None + return match.group("base"), match.group("index") + + +def collapsed_hashing_marker_id(marker_id: str) -> str: + """Return the marker id sample calling stores for a hashing antibody.""" + parts = split_hashing_marker_id(marker_id) + return parts[0] if parts is not None else str(marker_id) + + +def _hashing_marker_ids(panel_df: pd.DataFrame) -> set[str]: + """Return hashing marker ids, or an empty set when the column is absent.""" + if "sample_hashing" not in panel_df.columns: + return set() + mask = sample_hashing_mask(panel_df["sample_hashing"]) + return {str(marker_id) for marker_id in panel_df.index[mask]} diff --git a/src/pixelator/pna/config/panel/validation.py b/src/pixelator/pna/config/panel/validation.py new file mode 100644 index 000000000..21cccf571 --- /dev/null +++ b/src/pixelator/pna/config/panel/validation.py @@ -0,0 +1,227 @@ +"""Schema checks for a PNA antibody panel dataframe. + +Copyright © 2022 Pixelgen Technologies AB. +""" + +from __future__ import annotations + +import re + +import pandas as pd +import polars as pl + +from pixelator.pna.config.panel.hashing import ( + _hashing_marker_ids, + collapsed_hashing_marker_id, + sample_hashing_mask, + split_hashing_marker_id, +) + +# required columns +_INDEX_COLUMN = "marker_id" +_INDEX_COLUMN_TYPE = str +_REQUIRED_COLUMNS = { + "control": bool, + "sequence_1": str, + "sequence_2": str, +} + +# and these should have unique values +_UNIQUE_COLUMNS = ["sequence_1", "sequence_2"] + +# UniProt accession naming convention. The trailing alternative allows an empty id. +_UNIPROT_ID_RE = re.compile( + r"^[OPQ][0-9][A-Z0-9]{3}[0-9]|[A-NR-Z][0-9]([A-Z][A-Z0-9]{2}[0-9]){1,2}|$" +) + + +def _uniprot_ids_are_valid(id_str: object) -> bool: + """Return whether every semicolon-separated UniProt id matches the convention.""" + return all(bool(_UNIPROT_ID_RE.match(part)) for part in str(id_str).split(";")) + + +def _validate_sequences(panel_df, sequence_col): + errors = [] + sequences = panel_df[sequence_col] + ref_length = len(sequences.iloc[0]) + if not sequences.apply(lambda x: len(x) == ref_length).all(): + errors.append(f"All {sequence_col} values must have the same length.") + + if not sequences.str.match("^[ATCG]*$").all(): + errors.append( + f"All {sequence_col} values must only contain ATCG characters. Offending values: " + f"{sequences[~sequences.str.match('^[ATCG]*$')].tolist()}" + ) + + return errors + + +def _validate_marker_names(panel_df): + errors = [] + if any(panel_df.index.str.contains("_")): + # Markers should not contain underscores since this messes + # things up with Seurat on the R side + errors.append( + "The marker_id column should not contain underscores. " + "Please use dashes instead. Offending values: " + f"{panel_df.index[panel_df.index.str.contains('_')]}" + ) + if any(panel_df.index.str.contains(r"\s")): + # Markers should not contain white-spaces since this causes + # issues in the demultiplexing step (and other places that + # might assume that marker names are single tokens) + problematic_lines = panel_df.index[panel_df.index.str.contains(r"\s")] + errors.append( + "The marker_id column should not contain white-spaces. " + "Please use dashes instead or remove the white-spaces. Offending values: " + f"{problematic_lines}" + ) + return errors + + +def _validate_structure(panel_df: pd.DataFrame) -> list[str]: + """Return the first error that makes the remaining checks meaningless.""" + if not set(_REQUIRED_COLUMNS).issubset(set(panel_df.columns)): + missing_columns = set(_REQUIRED_COLUMNS) - set(panel_df.columns) + return [f"Panel has missing required columns: {missing_columns}"] + + if panel_df.shape[0] == 0: + return ["Panel file is empty"] + + if panel_df.index.name != _INDEX_COLUMN: + return [f"`{_INDEX_COLUMN}` is missing or is not set as index"] + + return [] + + +def _validate_column_types(panel_df: pd.DataFrame) -> list[str]: + """Return errors for required columns whose dtype does not match the schema.""" + errors = [] + panel_pl_df = pl.from_pandas(panel_df, include_index=True) + for col, expected_type in ( + _REQUIRED_COLUMNS | {_INDEX_COLUMN: _INDEX_COLUMN_TYPE} + ).items(): + found_type = panel_pl_df[col].dtype.to_python() + if not found_type == expected_type: + errors.append( + f"Column {col} has incorrect type. Expected {expected_type}, got {found_type}" + ) + return errors + + +def _validate_unique_values(panel_df: pd.DataFrame) -> list[str]: + """Return errors when a unique column or the marker id repeats.""" + errors = [] + for col in _UNIQUE_COLUMNS: + if not len(panel_df[col].unique()) == len(panel_df[col]): + errors.append(f"All values in column: {col} were not unique") + + if panel_df.index.duplicated().any(): + duplicated = panel_df.index[panel_df.index.duplicated()].unique().tolist() + errors.append( + "All values in column: marker_id were not unique. " + f"Offending values: {duplicated}" + ) + return errors + + +def _validate_control_column(panel_df: pd.DataFrame) -> list[str]: + """Return an error when the control column is not boolean.""" + if panel_df["control"].dtype != bool: + return ["`control` column is not boolean"] + return [] + + +def _validate_uniprot_ids(panel_df: pd.DataFrame) -> list[str]: + """Return an error when a present UniProt id breaks the naming convention.""" + if "uniprot_id" not in panel_df.columns: + return [] + + bad_ids = panel_df[~panel_df["uniprot_id"].apply(_uniprot_ids_are_valid)][ + "uniprot_id" + ] + if len(bad_ids) > 0: + return [ + "Invalid UniProt IDs found." + "Please conform to the naming convention or remove the following IDs:" + f"{bad_ids.tolist()}" + ] + return [] + + +def _validate_hashing_marker_ids(panel_df: pd.DataFrame) -> list[str]: + """Return errors when hashing ids lack a ``-`` suffix or nest.""" + hashing_ids = _hashing_marker_ids(panel_df) + if not hashing_ids: + return [] + errors: list[str] = [] + missing_suffix = [ + marker_id + for marker_id in sorted(hashing_ids) + if split_hashing_marker_id(marker_id) is None + ] + if missing_suffix: + errors.append( + "Hashing marker ids must end with - (e.g. B2M-1). " + f"Offending values: {missing_suffix}" + ) + nested = sorted( + marker_id + for marker_id in hashing_ids + if (base := collapsed_hashing_marker_id(marker_id)) != marker_id + and base in hashing_ids + ) + if nested: + errors.append( + "Hashing marker ids must not collapse to another hashing id " + f"(e.g. B2M-1-1 next to B2M-1). Offending values: {nested}" + ) + return errors + + +def _validate_control_is_not_hashing(panel_df: pd.DataFrame) -> list[str]: + """Return an error when a marker is both a control and a hashing marker.""" + if ( + "sample_hashing" in panel_df.columns + and panel_df["control"].dtype == bool + and len( + both := panel_df.index[ + panel_df["control"] & sample_hashing_mask(panel_df["sample_hashing"]) + ] + ) + > 0 + ): + return [ + "Control markers cannot be hashing markers. " + f"Offending values: {sorted(map(str, both))}" + ] + else: + return [] + + +def validate_antibody_panel( + panel_df: pd.DataFrame, validate_types: bool = True +) -> list[str]: + """Validate antibody panel schema and content. + + Args: + panel_df: Dataframe containing panel markers and sequences. + validate_types: If True, validate dataframe column types. + + Returns: + A list of validation error messages. Empty means valid input. + """ + errors = _validate_structure(panel_df) + if errors: + return errors + if validate_types: + errors += _validate_column_types(panel_df) + errors += _validate_unique_values(panel_df) + errors += _validate_marker_names(panel_df) + errors += _validate_control_column(panel_df) + errors += _validate_uniprot_ids(panel_df) + errors += _validate_sequences(panel_df, "sequence_1") + errors += _validate_sequences(panel_df, "sequence_2") + errors += _validate_hashing_marker_ids(panel_df) + errors += _validate_control_is_not_hashing(panel_df) + return errors diff --git a/src/pixelator/pna/sample_calling/__init__.py b/src/pixelator/pna/sample_calling/__init__.py index fdcffc229..62ece7a4b 100644 --- a/src/pixelator/pna/sample_calling/__init__.py +++ b/src/pixelator/pna/sample_calling/__init__.py @@ -4,7 +4,6 @@ """ import logging -import re import tempfile from itertools import chain from pathlib import Path @@ -20,7 +19,7 @@ from pixelator import __version__ from pixelator.pna.analysis_engine import AnalysisManager, PerComponentTask from pixelator.pna.anndata import add_missing_adata_info, pna_edgelist_to_anndata -from pixelator.pna.config.panel import PNAAntibodyPanel +from pixelator.pna.config.panel import PNAAntibodyPanel, collapsed_hashing_marker_id from pixelator.pna.pixeldataset import PNAPixelDataset from pixelator.pna.pixeldataset.io import PixelFileWriter from pixelator.pna.sample_calling.hash_antibodies import HashedAntibodyMapping @@ -209,6 +208,19 @@ def _add_original_hash_counts_to_obs( old_adata.obs[f"original_hash_counts_{ab}"] = 0 +def _markers_after_hash_collapse( + panel: PNAAntibodyPanel, hashing_antibodies: set[str] +) -> list[str]: + """Return count-matrix markers after hashing ids collapse to their base name. + + Hashing clones are left out. Each collapsed base that is not already on the + panel is appended in sorted order, so edgelist counts under that name are kept. + """ + kept = [marker for marker in panel.markers if marker not in hashing_antibodies] + bases = {collapsed_hashing_marker_id(marker) for marker in hashing_antibodies} + return kept + sorted(bases - set(kept)) + + def _build_post_sample_calling_anndata( con: duckdb.DuckDBPyConnection, old_adata: anndata.AnnData, @@ -235,14 +247,15 @@ def _build_post_sample_calling_anndata( old_adata, hashing_antibody_mapping.hashing_antibodies ) - # Create the anndata object and remove all panel hashing markers from var - new_adata = pna_edgelist_to_anndata(con, panel) - non_hashing_markers = [ - marker - for marker in new_adata.var.index - if marker not in hashing_antibody_mapping.hashing_antibodies - ] - new_adata = new_adata[:, non_hashing_markers].copy() + # The edgelist already uses collapsed hashing ids (B2M-1 -> B2M). Include + # a collapsed base that the panel does not define, and leave the clones out. + new_adata = pna_edgelist_to_anndata( + con, + panel, + markers=_markers_after_hash_collapse( + panel, hashing_antibody_mapping.hashing_antibodies + ), + ) new_adata = add_missing_adata_info(new_adata, old_adata) # `sample` is a reserved, transient column added when reading a @@ -390,7 +403,7 @@ def sample_calling( { "hashed_marker": hashed_markers, "base_marker": [ - re.sub(r"-\d+$", "", marker) for marker in hashed_markers + collapsed_hashing_marker_id(marker) for marker in hashed_markers ], } ) diff --git a/tests/common/data_generator/molecules.py b/tests/common/data_generator/molecules.py index e58d175c0..7b56c0567 100644 --- a/tests/common/data_generator/molecules.py +++ b/tests/common/data_generator/molecules.py @@ -11,12 +11,20 @@ import numpy as np import polars as pl +from pixelator.pna.config.panel import split_hashing_marker_id from tests.common.data_generator.topology import generate_cell_graph if TYPE_CHECKING: from pixelator.pna.config.panel import PNAAntibodyPanel +def _hashing_index(marker_id: str) -> int: + parts = split_hashing_marker_id(marker_id) + if parts is None: + raise ValueError(f"Hashing marker {marker_id!r} must end with -.") + return int(parts[1]) + + def generate_edgelist( n_cells: int, n_nodes: int, @@ -197,7 +205,7 @@ def _hashing_indices_per_cell( if hashing.size == 0: return [None] * n_cells if hashing_indices is None: - indices = np.unique([int(m.rsplit("-", 1)[-1]) for m in hashing]) + indices = np.unique([_hashing_index(str(m)) for m in hashing]) else: indices = np.unique(hashing_indices) cell_indices = np.resize(indices, n_cells) @@ -324,7 +332,7 @@ def _assign_markers( return node_umi_map hashing_markers = markers[is_hashing] - index = np.array([int(m.rsplit("-", 1)[-1]) for m in hashing_markers]) + index = np.array([_hashing_index(str(m)) for m in hashing_markers]) chosen = hashing_markers[index == hashing_index] if chosen.size == 0: return node_umi_map diff --git a/tests/common/data_generator/test_molecules.py b/tests/common/data_generator/test_molecules.py index 003ffa261..fe0814a79 100644 --- a/tests/common/data_generator/test_molecules.py +++ b/tests/common/data_generator/test_molecules.py @@ -11,6 +11,7 @@ from polars.testing import assert_frame_equal from pixelator.pna.config import pna_config +from pixelator.pna.config.panel import split_hashing_marker_id from tests.common.data_generator.molecules import ( _assign_markers, _assign_umis, @@ -154,7 +155,7 @@ def test_assign_markers_hashing_index(marker_panel): by_index = defaultdict(list) for marker, hashing in zip(markers, is_hashing): if hashing: - by_index[int(marker.rsplit("-", 1)[-1])].append(marker) + by_index[int(split_hashing_marker_id(marker)[1])].append(marker) # the chosen index's markers share the extra ~f and sit above the base chosen = by_index[hashing_index] diff --git a/tests/pna/config/test_panel.py b/tests/pna/config/test_panel.py index 6a2b7416b..9bfeabfdc 100644 --- a/tests/pna/config/test_panel.py +++ b/tests/pna/config/test_panel.py @@ -234,3 +234,108 @@ def test_panel_header_non_recoverable_yaml_still_fails(): with pytest.raises(yaml.YAMLError): PNAAntibodyPanel.from_csv(tmp_file.name) + + +def _marker_row(marker_id: str, sequence: str, **extra) -> dict: + row = { + "marker_id": marker_id, + "control": False, + "sequence_1": sequence, + "sequence_2": sequence, + } + row.update(extra) + return row + + +def test_duplicate_marker_ids_are_rejected(): + """Marker ids must be unique even when sequences differ.""" + frame = pd.DataFrame( + [_marker_row("CD3", "AAAA"), _marker_row("CD3", "CCCC")] + ).set_index("marker_id") + with pytest.raises(AssertionError, match="marker_id were not unique"): + PNAAntibodyPanel( + frame, + AntibodyPanelMetadata(name="base", version="1.0.0"), + ) + + +def test_hashing_marker_ids_need_a_numeric_suffix_and_must_not_nest(): + """Hashing ids end with - and must not collapse onto each other.""" + missing_suffix = pd.DataFrame( + [_marker_row("B2M", "AAAA", sample_hashing=True)] + ).set_index("marker_id") + with pytest.raises(AssertionError, match="must end with -"): + PNAAntibodyPanel( + missing_suffix, + AntibodyPanelMetadata(name="base", version="1.0.0"), + ) + + nested = pd.DataFrame( + [ + _marker_row("B2M-1", "AAAA", sample_hashing=True), + _marker_row("B2M-1-1", "CCCC", sample_hashing=True), + ] + ).set_index("marker_id") + with pytest.raises(AssertionError, match="must not collapse"): + PNAAntibodyPanel( + nested, + AntibodyPanelMetadata(name="base", version="1.0.0"), + ) + + +def test_control_markers_cannot_be_hashing_markers(): + """A marker cannot be both a control and a hashing marker.""" + frame = pd.DataFrame( + [_marker_row("B2M-1", "AAAA", control=True, sample_hashing=True)] + ).set_index("marker_id") + with pytest.raises(AssertionError, match="cannot be hashing markers"): + PNAAntibodyPanel( + frame, + AntibodyPanelMetadata(name="base", version="1.0.0"), + ) + + +def test_panel_validation_reports_missing_required_columns(panel_df): + """A panel missing a required column is rejected before later checks.""" + errors = PNAAntibodyPanel.validate_antibody_panel( + panel_df.drop(columns=["sequence_1"]) + ) + assert len(errors) == 1 + assert "missing required columns" in errors[0] + assert "sequence_1" in errors[0] + + +def test_panel_validation_reports_an_empty_panel(panel_df): + """An empty panel is rejected.""" + errors = PNAAntibodyPanel.validate_antibody_panel(panel_df.iloc[0:0]) + assert "Panel file is empty" in errors + + +def test_panel_validation_requires_marker_id_index(panel_df): + """The marker id must be the dataframe index.""" + panel_df.index.name = "other" + errors = PNAAntibodyPanel.validate_antibody_panel(panel_df, validate_types=False) + assert "`marker_id` is missing or is not set as index" in errors + + +def test_panel_validation_requires_unique_sequences(panel_df): + """Sequence columns must contain unique values.""" + panel_df.loc["marker2", "sequence_1"] = panel_df.loc["marker1", "sequence_1"] + errors = PNAAntibodyPanel.validate_antibody_panel(panel_df) + assert "All values in column: sequence_1 were not unique" in errors + + +def test_panel_validation_requires_boolean_control(panel_df): + """The control column must be boolean once types are not re-checked.""" + panel_df["control"] = ["no", "yes", "no"] + errors = PNAAntibodyPanel.validate_antibody_panel(panel_df, validate_types=False) + assert "`control` column is not boolean" in errors + + +def test_panel_validation_reports_column_type_mismatch(panel_df): + """Column types are checked when validate_types is left on.""" + panel_df["control"] = ["no", "yes", "no"] + errors = PNAAntibodyPanel.validate_antibody_panel(panel_df) + assert any( + error.startswith("Column control has incorrect type.") for error in errors + ) diff --git a/tests/pna/sample_calling/test_sample_calling.py b/tests/pna/sample_calling/test_sample_calling.py index f8752eba5..a9d0330bf 100644 --- a/tests/pna/sample_calling/test_sample_calling.py +++ b/tests/pna/sample_calling/test_sample_calling.py @@ -15,11 +15,12 @@ from pixelator.common.config import AntibodyPanelMetadata from pixelator.pna import read from pixelator.pna.anndata import pna_edgelist_to_anndata -from pixelator.pna.config.panel import PNAAntibodyPanel +from pixelator.pna.config.panel import PNAAntibodyPanel, collapsed_hashing_marker_id from pixelator.pna.pixeldataset.io import PixelFileWriter from pixelator.pna.sample_calling import ( _add_original_hash_counts_to_obs, _collect_nodes_to_remove, + _markers_after_hash_collapse, collect_hash_info, create_final_report, sample_calling, @@ -28,6 +29,35 @@ from pixelator.pna.sample_calling.hash_antibodies import HashedAntibodyMapping +def _panel(marker_ids: list[str], hashing: set[str]) -> PNAAntibodyPanel: + frame = pd.DataFrame( + { + "marker_id": marker_ids, + "control": [False] * len(marker_ids), + "sequence_1": ["AAAA", "AAAC", "AAAG", "AAAT", "AACA"][: len(marker_ids)], + "sequence_2": ["CCCC", "CCCA", "CCCG", "CCCT", "CCAC"][: len(marker_ids)], + "sample_hashing": [marker_id in hashing for marker_id in marker_ids], + } + ).set_index("marker_id") + return PNAAntibodyPanel( + frame, AntibodyPanelMetadata(name="hashing", version="1.0.0") + ) + + +def test_markers_after_hash_collapse_appends_missing_bases(): + """Hashing clones drop out, and each new collapsed base is appended once, sorted.""" + panel = _panel( + ["PD-1", "B2M", "CD29-1", "CD29-2", "ACTB-1"], + {"CD29-1", "CD29-2", "ACTB-1"}, + ) + assert _markers_after_hash_collapse(panel, {"CD29-1", "CD29-2", "ACTB-1"}) == [ + "PD-1", + "B2M", + "ACTB", + "CD29", + ] + + def test_add_original_hash_counts_includes_all_panel_antibodies(): """All panel hashing antibodies get original_hash_counts_* in obs, including those not in adata. @@ -386,7 +416,7 @@ def test_sample_calling(sample_hashed_pixel_files, tmp_path): assert dehashed_counts.shape[0] == 10 sample_name = list(dehashed_pxl.metadata().keys())[0] for ab in hashed_antibodies[sample_name]: - base_name = ab.split("-")[0] + base_name = collapsed_hashing_marker_id(ab) assert all( original_counts.loc[ dehashed_counts.index, @@ -520,6 +550,70 @@ def test_sample_calling_does_not_strip_suffix_from_non_hash_markers( assert dehashed_counts.loc["c1", "PD-1"] == original_counts.loc["c1", "PD-1"] +def test_sample_calling_keeps_collapsed_hash_marker_missing_from_panel( + tmp_path: Path, +): + """A collapsed hashing id stays in var when the panel has no such marker.""" + panel_df = pd.DataFrame( + [ + { + "marker_id": "PD-1", + "control": False, + "uniprot_id": "P00001", + "sequence_1": "ATCGATCGAA", + "sequence_2": "ATCGATCGAC", + }, + { + "marker_id": "HashA-1", + "control": False, + "uniprot_id": "P00003", + "sequence_1": "ATCGATCGTT", + "sequence_2": "ATCGATCGTG", + }, + ] + ).set_index("marker_id") + panel = PNAAntibodyPanel( + df=panel_df, + metadata=AntibodyPanelMetadata(name="test-panel", version="0.1.0"), + ) + edgelist = pl.DataFrame( + { + "umi1": pl.Series([1, 2, 3], dtype=pl.UInt64), + "umi2": pl.Series([2, 3, 4], dtype=pl.UInt64), + "read_count": pl.Series([10, 10, 10], dtype=pl.UInt32), + "uei_count": pl.Series([5, 5, 5], dtype=pl.UInt32), + "marker_1": ["PD-1", "PD-1", "HashA-1"], + "marker_2": ["PD-1", "HashA-1", "HashA-1"], + "component": ["c1", "c1", "c1"], + } + ) + target = tmp_path / "input.pxl" + with PixelFileWriter(target) as writer: + writer.write_edgelist(edgelist) + con = writer.get_connection() + adata = pna_edgelist_to_anndata(con, panel=panel) + writer.write_adata(adata) + writer.write_metadata({"sample_name": "input", "version": "0.1.0"}) + + out_dir = tmp_path / "out" + out_dir.mkdir() + sample_calling( + input_pxl=read(target), + hashing_antibody_mapping=HashedAntibodyMapping( + mapping={"S1": ["HashA-1"]}, + all_hashing_antibodies=["HashA-1"], + ), + output_folder=out_dir, + remove_incompatible=True, + enrichment_threshold=1.5, + ) + + dehashed_counts = read(next(out_dir.glob("*.dehashed.pxl"))).adata().to_df() + assert "HashA" not in panel.markers + assert "HashA-1" not in dehashed_counts.columns + assert dehashed_counts.loc["c1", "HashA"] == 3 + + @pytest.mark.slow def test_sample_calling_without_uei_count(tmp_path: Path): """sample_calling succeeds when the input edgelist has no uei_count column."""