"""JUMP Cell Painting profiles and their perturbation annotation.
A JUMP plate parquet carries three metadata columns (source, plate, well) and 4762 features.
The compound or gene each well received is recorded in a separate repository, keyed by ``Metadata_JCP2022``.
Sources, all public over HTTPS:
* profiles: Cell Painting Gallery, accession ``cpg0016-jump``
* annotation: ``jump-cellpainting/datasets`` on GitHub
References:
:cite:t:`Chandrasekaran_2023`, the JUMP Cell Painting datasets.
"""
from __future__ import annotations
import warnings
from collections.abc import Sequence
from pathlib import Path
from typing import Any
import numpy as np
import pandas as pd
from anndata import AnnData
from mantispy._core.logging import get_logger
from mantispy.io._profiles import _read_frame, read_profiles
#: Pinned by sha256 in the dataset registry because the upstream repository is mutable.
TABLES = (
"plate",
"well",
"compound",
"crispr",
"orf",
"perturbation_control",
"gene_chromosome_map",
)
#: JUMP's negative control: DMSO, under its JCP identifier.
NEGATIVE_CONTROL = "JCP2022_033924"
KINDS = ("compound", "crispr", "orf")
_JOIN_ON = ["Metadata_Source", "Metadata_Plate", "Metadata_Well"]
_ON_UNMAPPED = ("warn", "drop", "error")
def jump_metadata(name: str) -> pd.DataFrame:
"""Read one of JUMP's annotation tables, such as ``"well"``, ``"compound"`` or ``"crispr"``.
Args:
name: Which table to read, one of :data:`TABLES`.
Returns:
The table as JUMP publishes it, downloaded once into :attr:`mantispy.settings.cache_dir` and checked against the sha256 the dataset registry pins.
Raises:
ValueError: `name` is not one of :data:`TABLES`.
"""
if name not in TABLES:
raise ValueError(f"name must be one of {TABLES}, got {name!r}")
from mantispy.ds._datasets import _files
(path,) = _files("_jump_annotation", select=lambda file_name: file_name.split(".")[0] == f"jump_{name}")
return _read_frame(path)
[docs]
def read_jump(paths: str | Path | Sequence[str | Path], annotate: bool = True, **kwargs: Any) -> AnnData:
"""Read JUMP plate profiles, optionally joining the annotation.
Args:
paths: One or more ``{plate}.parquet`` files in the Cell Painting Gallery layout.
annotate: Join the well and compound tables, which map the three metadata columns to a perturbation.
Downloads once and caches it.
kwargs: Passed to :func:`~mantispy.io.read_profiles`.
``on_column_mismatch="intersect"`` is useful when plates come from different sources.
Returns:
An :class:`~anndata.AnnData` at well resolution.
With ``annotate`` it carries ``Metadata_JCP2022``, ``Metadata_Perturbation``, ``Metadata_Perturbation_Type`` (``"compound"``), ``Metadata_InChIKey`` and ``Metadata_Control``.
"""
adata = read_profiles(paths, resolution="well", **kwargs)
if annotate:
from mantispy.pp._annotate import annotate_jump
annotate_jump(adata)
return adata
def join_jump_annotation(obs: pd.DataFrame, kind: str = "compound", on_unmapped: str = "warn") -> pd.DataFrame:
"""Join the JUMP annotation onto an ``obs`` frame keyed by source, plate and well.
Args:
obs: A frame carrying ``Metadata_Source``, ``Metadata_Plate`` and ``Metadata_Well``, whose three key columns are cast to strings in place before the join, or one that already carries ``Metadata_JCP2022``, as JUMP's assembled profiles do.
kind: Which perturbation the annotation is read for, one of :data:`KINDS`.
on_unmapped: What to do with wells whose id the annotation cannot map: ``"warn"`` keeps them, ``"drop"`` removes them and warns, ``"error"`` raises.
Returns:
A new frame with ``Metadata_JCP2022``, ``Metadata_Perturbation``, ``Metadata_Perturbation_Type`` and ``Metadata_Control`` joined onto `obs`, missing on the wells the annotation does not cover.
For ``"compound"`` the perturbation is the ``Metadata_JCP2022`` compound id, ``Metadata_Perturbation_Type`` is ``"compound"``, it adds ``Metadata_InChIKey``, and ``Metadata_Control_Type`` is ``"negcon"`` for :data:`NEGATIVE_CONTROL` else ``"treatment"``, with ``Metadata_Control`` the boolean it implies.
For ``"crispr"`` and ``"orf"`` ``Metadata_Gene`` and ``Metadata_Perturbation`` are the gene symbol (its reagents are the replicates), ``Metadata_Perturbation_Type`` is the kind, ``Metadata_Control_Type`` is ``"negcon"``, ``"poscon"`` or ``"treatment"`` (JUMP's ``"trt"``), ``Metadata_Control`` marks any control row, and ``Metadata_ChromosomeArm`` is the arm the gene sits on, such as ``"1p"``, missing for a gene without a mapped locus.
Raises:
ValueError: `kind` is not one of :data:`KINDS`, `on_unmapped` is not ``"warn"``, ``"drop"`` or ``"error"``, or it is ``"error"`` and an id is unmapped.
KeyError: `obs` has no ``Metadata_JCP2022`` and is missing one of the three columns the annotation is keyed by.
"""
if kind not in KINDS:
raise ValueError(f"kind must be one of {KINDS}, got {kind!r}")
if on_unmapped not in _ON_UNMAPPED:
raise ValueError(f"on_unmapped must be one of {_ON_UNMAPPED}, got {on_unmapped!r}")
# The JUMP annotation is authoritative for control status; drop any incoming values so the controls
# merge below does not collide with them into Metadata_Control_Type_x / _y.
obs = obs.drop(columns=["Metadata_Control", "Metadata_Control_Type"], errors="ignore")
joined = obs if "Metadata_JCP2022" in obs else _join_wells(obs)
if kind == "compound":
compounds = jump_metadata("compound")[["Metadata_JCP2022", "Metadata_InChIKey"]]
joined = joined.merge(compounds, on="Metadata_JCP2022", how="left", validate="m:1")
joined["Metadata_Perturbation"] = joined["Metadata_JCP2022"].astype(str)
joined["Metadata_Perturbation_Type"] = "compound"
is_negcon = (joined["Metadata_JCP2022"] == NEGATIVE_CONTROL).to_numpy()
joined["Metadata_Control_Type"] = np.where(is_negcon, "negcon", "treatment")
joined["Metadata_Control"] = is_negcon
return _handle_unmapped(
joined, joined["Metadata_InChIKey"].isna().to_numpy() & ~is_negcon, "compound", on_unmapped
)
genes = jump_metadata(kind)[["Metadata_JCP2022", "Metadata_Symbol"]].rename(
columns={"Metadata_Symbol": "Metadata_Gene"}
)
controls = jump_metadata("perturbation_control")
controls = controls.loc[controls["Metadata_modality"] == kind, ["Metadata_JCP2022", "Metadata_pert_type"]].rename(
columns={"Metadata_pert_type": "Metadata_Control_Type"}
)
for table in (genes, controls):
joined = joined.merge(table, on="Metadata_JCP2022", how="left", validate="m:1")
# Spec 10.2 fixes the control-type vocabulary; JUMP's non-control pert_type "trt" is the schema's "treatment".
joined["Metadata_Control_Type"] = joined["Metadata_Control_Type"].fillna("treatment").replace({"trt": "treatment"})
# The gene is the replication unit here: JUMP's mAP benchmark scores the genetic arms at the gene level,
# treating the several reagents per gene as its replicates. The reagent stays in Metadata_JCP2022.
joined["Metadata_Perturbation"] = joined["Metadata_Gene"].fillna(joined["Metadata_JCP2022"]).astype(str)
joined["Metadata_Perturbation_Type"] = kind
joined["Metadata_Control"] = (joined["Metadata_Control_Type"] != "treatment").to_numpy()
joined = _handle_unmapped(joined, joined["Metadata_Gene"].isna().to_numpy(), "gene", on_unmapped)
joined["Metadata_ChromosomeArm"] = joined["Metadata_Gene"].map(_chromosome_arms())
return joined
def _handle_unmapped(joined: pd.DataFrame, unmapped: np.ndarray, what: str, on_unmapped: str) -> pd.DataFrame:
"""Warn, drop or raise for annotated ids that have no ``what`` in the JUMP tables; wells with no id are logged by the well join."""
unmapped = unmapped & joined["Metadata_JCP2022"].notna().to_numpy()
n_ids = len(joined.loc[unmapped, "Metadata_JCP2022"].unique())
if not n_ids:
return joined
message = f"{n_ids} JUMP id(s) could not map to a {what}"
if on_unmapped == "error":
raise ValueError(f"{message}; pass on_unmapped='drop' to remove their wells or 'warn' to keep them.")
if on_unmapped == "drop":
warnings.warn(f"{message}; their {int(unmapped.sum())} well(s) were dropped.", UserWarning, stacklevel=3)
return joined.loc[~unmapped]
warnings.warn(f"{message}; their wells are kept as treatments without one.", UserWarning, stacklevel=3)
return joined
def _join_wells(obs: pd.DataFrame) -> pd.DataFrame:
"""Map source, plate and well to ``Metadata_JCP2022`` through JUMP's well table."""
missing = [column for column in _JOIN_ON if column not in obs]
if missing:
raise KeyError(
f"obs is missing {missing}, the columns the JUMP annotation is keyed by. Read the plate "
"parquet with mt.io.read_jump to get them."
)
wells = jump_metadata("well")
for frame in (obs, wells):
for column in _JOIN_ON:
frame[column] = frame[column].astype(str)
joined = obs.merge(wells, on=_JOIN_ON, how="left", validate="m:1")
unannotated = int(joined["Metadata_JCP2022"].isna().sum())
if unannotated:
get_logger().warning(
"%d of %d wells have no JUMP annotation; their Metadata_JCP2022 is missing", unannotated, len(joined)
)
return joined
def _chromosome_arms() -> pd.Series:
"""The chromosome arm of every gene symbol, read off its cytogenetic locus as jump-profiling-recipe does."""
loci = jump_metadata("gene_chromosome_map").drop_duplicates("Approved_symbol").set_index("Approved_symbol")["Locus"]
return loci.astype(str).str.extract(r"^(\w+?[pq])", expand=False)
[docs]
def unexpressed_genes(
expression: str | Path | pd.DataFrame,
cell_line: str,
*,
tpm_cutoff: float = 0.0,
) -> set[str]:
"""Gene symbols a cell line does not express, read from DepMap as an empirical null for hit calling.
The reference is always a DepMap expression matrix you have downloaded, so the null is tied to the screened cell line rather than to any hard-coded reference.
Nothing is downloaded or re-hosted; point it at the file from the DepMap data page, such as ``OmicsExpressionProteinCodingGenesTPMLogp1.csv`` (one row per model, gene columns named ``"SYMBOL (ENTREZ)"``, holding log2(TPM+1)).
The cell line is named by its DepMap model id, DepMap's own key, so the lookup reuses their identifier directly.
A gene is unexpressed when its value is at or below `tpm_cutoff`, which is zero for the zero-TPM genes PERISCOPE uses.
Args:
expression: A DepMap expression matrix, as a path or a loaded frame indexed by model id.
cell_line: The DepMap model id to read, such as ``"ACH-000364"`` for U2OS.
tpm_cutoff: The log2(TPM+1) value at or below which a gene is called unexpressed.
Returns:
The set of unexpressed gene symbols.
Raises:
ValueError: `cell_line` is not a row of the expression matrix.
"""
frame = expression if isinstance(expression, pd.DataFrame) else pd.read_csv(expression, index_col=0)
if cell_line not in frame.index:
raise ValueError(f"cell_line={cell_line!r} is not a model id in the expression matrix")
symbols = frame.columns.to_series().str.replace(r"\s*\(\d+\)\s*$", "", regex=True)
row = pd.Series(np.asarray(frame.loc[cell_line]).ravel(), index=frame.columns)
values = pd.to_numeric(row, errors="coerce").to_numpy()
return set(symbols[values <= tpm_cutoff].astype(str))
