fix(datasets): preserve zeros and propagate NaN in williams_2018

HISTORY.md

@@ -7,8 +7,29 @@ and this project adheres to
 
 ## [Unreleased]
 
+### Fixed
+
+- **Datasets** (`pr.datasets`) and **Download** (`pr.download`):
+  `williams_2018()` no longer conflates measured zeros with missing
+  values. Three defects, each masked by another:
+  - zeros in `.X` were coerced to `np.nan`, discarding 13,547 genuine
+    measurements
+  - charge-state summation used pandas' default `min_count=0`, so a
+    group with no measurements at all summed to `0.0`, inventing 3,324
+  - the same summation skipped `NaN` inside a *partially* measured
+    group, reporting a partial total as complete for 260 cells
+  - **this changes `.X`**, and therefore the output of
+    `pr.download.williams_2018()`. Verified against the PRIDE deposit
+    used by the original publication: the two are now bit-identical
+    with an identical missingness pattern.
+
 ### Changed
 
+- **Datasets** (`pr.datasets`) and **Download** (`pr.download`):
+  `williams_2018()` gained a `zero_to_na` parameter (default `False`),
+  for consistency with sibling functions. It is mutually exclusive
+  with `fill_na`. Note it governs zero semantics only and does not
+  restore the summation defects above.
 - **Preprocessing** (`pr.pp`): `normalize_median()`
   - now defaults to `log_space=True`
   - renamed the `batch_id` parameter to `group_by`

proteopy/datasets/williams_2018.py

@@ -16,6 +16,7 @@
 
 
 def williams_2018(
+    zero_to_na: bool = False,
     fill_na: float | int | None = None,
 ) -> ad.AnnData:
     """Load Williams 2018 mouse multi-tissue proteomics dataset.
@@ -30,6 +31,21 @@ def williams_2018(
     per peptide sequence. By default, missing values
     are represented as ``np.nan``.
 
+    Missing values and zeros are kept distinct, because
+    they mean different things:
+
+    - A **zero** is a measurement, and is preserved
+      as ``0.0``.
+    - A **missing** value is the absence of one.
+      Charge-state summation propagates it, so a peptide
+      is quantified in a sample only if *every* one of
+      its charge states was quantified there; a
+      partially measured group yields ``np.nan`` rather
+      than a partial total.
+
+    Pass ``zero_to_na=True`` to treat zeros as missing
+    instead.
+
     Sample annotation (``.obs``) includes:
         - ``sample_id``: Unique sample identifier
         - ``tissue``: Tissue type (Brain, BAT, Heart, Liver, Quad)
@@ -45,9 +61,14 @@ def williams_2018(
 
     Parameters
     ----------
+    zero_to_na : bool, optional
+        If True, zeros in ``.X`` are treated as missing
+        values (NaN). Mutually exclusive with
+        ``fill_na``.
     fill_na : float | int | None, optional
         If not ``None``, replace ``np.nan`` in ``.X``
-        with this value.
+        with this value. Mutually exclusive with
+        ``zero_to_na``.
 
     Returns
     -------
@@ -78,13 +99,20 @@ def williams_2018(
        Cellular Proteomics, 2018, 17(9):1766-1777.
        DOI: 10.1074/mcp.RA118.000554.
     """
+    if not isinstance(zero_to_na, bool):
+        raise TypeError(
+            f"zero_to_na must be bool, got {type(zero_to_na).__name__}"
+        )
     if fill_na is not None and not isinstance(
-        fill_na, (int, float),
+        fill_na,
+        (int, float),
     ):
         raise TypeError(
             f"fill_na must be float, int, or None, "
             f"got {type(fill_na).__name__}"
         )
+    if zero_to_na and fill_na is not None:
+        raise ValueError("`zero_to_na` and `fill_na` are mutually exclusive.")
 
     url = (
         "https://ars.els-cdn.com/content/image/"
@@ -124,7 +152,8 @@ def williams_2018(
     # Select intensity columns: named cols where row 0 == "Intensity",
     # excluding _mito fractions
     intensity_cols = [
-        c for c in df.columns
+        c
+        for c in df.columns
         if "Unnamed" not in str(c)
         and df[c].iloc[0] == "Intensity"
         and "_mito" not in str(c)
@@ -133,29 +162,22 @@ def williams_2018(
     df = df[list(meta_cols.keys()) + intensity_cols]
 
     # Remove _WholeCell suffix from sample column names
-    df = df.rename(columns={
-        c: c.replace("_WholeCell", "")
-        for c in intensity_cols
-    })
+    df = df.rename(
+        columns={c: c.replace("_WholeCell", "") for c in intensity_cols}
+    )
     df = df.rename(columns=meta_cols)
 
     # Drop the first row (secondary header)
     df = df.iloc[1:].reset_index(drop=True)
 
     # Extract peptide sequence (remove prefixes and suffixes)
-    df["peptide_id"] = (
-        df["peptide_id"].str.split("_").str[1]
-    )
+    df["peptide_id"] = df["peptide_id"].str.split("_").str[1]
 
     # Verify protein_id and gene_id are consistent
     # across charge states of the same peptide
-    meta_check = (
-        df.groupby("peptide_id")[["protein_id", "gene_id"]]
-        .nunique()
-    )
+    meta_check = df.groupby("peptide_id")[["protein_id", "gene_id"]].nunique()
     inconsistent = meta_check[
-        (meta_check["protein_id"] > 1)
-        | (meta_check["gene_id"] > 1)
+        (meta_check["protein_id"] > 1) | (meta_check["gene_id"] > 1)
     ]
     if not inconsistent.empty:
         raise ValueError(
@@ -164,45 +186,69 @@ def williams_2018(
             f"{inconsistent.index.tolist()}"
         )
 
-    # Sum intensities across charge states of the same peptide
+    # Sum intensities across charge states of the same peptide.
+    #
+    # The sum PROPAGATES missing values: a peptide is quantified in a
+    # sample only if every one of its charge states was quantified
+    # there. `DataFrameGroupBy.sum()` cannot express this on its own --
+    # it has no `skipna` argument, and `min_count` governs only the
+    # all-missing case -- so the sum is masked on the count.
+    #
+    # Both halves matter, and neither is sufficient alone:
+    #
+    #   min_count=1     without it, a group whose charge states are ALL
+    #                   missing sums to 0.0, inventing a measurement
+    #                   that was never made.
+    #   .where(complete) without it, a PARTIALLY measured group reports
+    #                   a partial total as though it were complete --
+    #                   e.g. [NaN, 5000] -> 5000.
     sample_cols = [
-        c for c in df.columns
+        c
+        for c in df.columns
         if c not in ("peptide_id", "protein_id", "gene_id")
     ]
     df[sample_cols] = df[sample_cols].astype(float)
-    var = (
-        df.groupby("peptide_id")[["protein_id", "gene_id"]]
-        .first()
-    )
+    grouped = df.groupby("peptide_id")
+    var = grouped[["protein_id", "gene_id"]].first()
     var["peptide_id"] = var.index
-    X = (
-        df.groupby("peptide_id")[sample_cols]
-        .sum()
-        .values.T
-    )
+
+    intensities = grouped[sample_cols]
+    complete = intensities.count().eq(grouped.size(), axis=0)
+    X = intensities.sum(min_count=1).where(complete).values.T
 
     # Build obs annotation with tissue and mouse_id
     obs = pd.DataFrame({"sample_id": sample_cols})
     parts = obs["sample_id"].str.split(
-        "_", n=1, expand=True,
+        "_",
+        n=1,
+        expand=True,
     )
     parts.columns = ["p1", "p2"]
-    tissue_first = parts["p1"].str.fullmatch(
-        r"Brain|BAT|Heart|Liver|Quad"
-    )
+    tissue_first = parts["p1"].str.fullmatch(r"Brain|BAT|Heart|Liver|Quad")
     obs["tissue"] = np.where(
-        tissue_first, parts["p1"], parts["p2"],
+        tissue_first,
+        parts["p1"],
+        parts["p2"],
     )
     obs["mouse_id"] = np.where(
-        tissue_first, parts["p2"], parts["p1"],
+        tissue_first,
+        parts["p2"],
+        parts["p1"],
     )
     obs = obs.set_index("sample_id")
     obs.index.name = None
     obs["sample_id"] = obs.index
 
-    # Construct anndata
+    # Construct anndata.
+    #
+    # NOTE: by default zeros are NOT coerced to NaN. A zero in this
+    # dataset is a measurement -- the peptide was looked for and its
+    # intensity was zero -- and is not interchangeable with "not
+    # measured". 13,547 of the 1,307,600 cells are genuine zeros.
     adata = ad.AnnData(X=X, obs=obs, var=var)
-    adata.X[adata.X == 0] = np.nan
+
+    if zero_to_na:
+        adata.X[adata.X == 0] = np.nan
 
     if fill_na is not None:
         adata.X[np.isnan(adata.X)] = fill_na

proteopy/download/williams_2018.py

@@ -7,16 +7,13 @@
 
 
 _DEFAULT_INTENSITIES = (
-    "williams-2018_ms-proteomics"
-    "_mouse-tissue_intensities.tsv"
+    "williams-2018_ms-proteomics" "_mouse-tissue_intensities.tsv"
 )
 _DEFAULT_VAR = (
-    "williams-2018_ms-proteomics"
-    "_mouse-tissue_peptide-annotation.tsv"
+    "williams-2018_ms-proteomics" "_mouse-tissue_peptide-annotation.tsv"
 )
 _DEFAULT_SAMPLE = (
-    "williams-2018_ms-proteomics"
-    "_mouse-tissue_sample-annotation.tsv"
+    "williams-2018_ms-proteomics" "_mouse-tissue_sample-annotation.tsv"
 )
 
 
@@ -25,6 +22,7 @@ def _check_williams_2018_types(
     var_annotation_path,
     sample_annotation_path,
     sep,
+    zero_to_na,
     fill_na,
     force,
 ):
@@ -36,27 +34,27 @@ def _check_williams_2018_types(
     ):
         if not isinstance(value, (str, Path)):
             raise TypeError(
-                f"{name} must be str or Path, "
-                f"got {type(value).__name__}"
+                f"{name} must be str or Path, " f"got {type(value).__name__}"
             )
     if sep is not None and not isinstance(sep, str):
         raise TypeError(
-            f"sep must be str or None, "
-            f"got {type(sep).__name__}"
+            f"sep must be str or None, " f"got {type(sep).__name__}"
+        )
+    if not isinstance(zero_to_na, bool):
+        raise TypeError(
+            f"zero_to_na must be bool, " f"got {type(zero_to_na).__name__}"
         )
     if fill_na is not None and (
-        isinstance(fill_na, bool)
-        or not isinstance(fill_na, (int, float))
+        isinstance(fill_na, bool) or not isinstance(fill_na, (int, float))
     ):
         raise TypeError(
             f"fill_na must be float, int, or None, "
             f"got {type(fill_na).__name__}"
         )
+    if zero_to_na and fill_na is not None:
+        raise ValueError("`zero_to_na` and `fill_na` are mutually exclusive.")
     if not isinstance(force, bool):
-        raise TypeError(
-            f"force must be bool, "
-            f"got {type(force).__name__}"
-        )
+        raise TypeError(f"force must be bool, " f"got {type(force).__name__}")
 
 
 def _check_williams_2018_paths(
@@ -101,6 +99,7 @@ def williams_2018(
     sample_annotation_path: str | Path = _DEFAULT_SAMPLE,
     *,
     sep: str | None = None,
+    zero_to_na: bool = False,
     fill_na: float | int | None = None,
     force: bool = False,
 ) -> None:
@@ -135,9 +134,13 @@ def williams_2018(
         separator is inferred from each file extension via
         ``detect_separator_from_extension()``
         (``.tsv`` → tab, ``.csv`` → comma).
+    zero_to_na : bool, optional
+        If True, zeros in the intensity matrix are treated as missing
+        values (NaN). Mutually exclusive with ``fill_na``.
     fill_na : float | int | None, optional
         If not ``None``, replace NaN values in the long-format
-        intensities DataFrame with this value before saving.
+        intensities DataFrame with this value before saving. Mutually
+        exclusive with ``zero_to_na``.
     force : bool, optional
         If ``True``, overwrite existing files at the output
         paths. Otherwise, raise ``FileExistsError`` when a
@@ -171,6 +174,7 @@ def williams_2018(
         var_annotation_path,
         sample_annotation_path,
         sep,
+        zero_to_na,
         fill_na,
         force,
     )
@@ -183,7 +187,7 @@ def williams_2018(
         )
     )
 
-    adata = _load_williams_2018()
+    adata = _load_williams_2018(zero_to_na=zero_to_na)
 
     # Auto-detect separator from file extension if not provided
     if sep is None:
@@ -222,33 +226,31 @@ def williams_2018(
         intensities_path,
         sep=sep_intensities,
         index=False,
-        lineterminator='\n',
+        lineterminator="\n",
     )
 
     # Save .var annotation
-    df_var = adata.var[
-        ["peptide_id", "protein_id", "gene_id"]
-    ].copy()
+    df_var = adata.var[["peptide_id", "protein_id", "gene_id"]].copy()
     var_annotation_path.parent.mkdir(
-        parents=True, exist_ok=True,
+        parents=True,
+        exist_ok=True,
     )
     df_var.to_csv(
         var_annotation_path,
         sep=sep_var,
         index=False,
-        lineterminator='\n',
+        lineterminator="\n",
     )
 
     # Save .obs annotation
-    df_obs = adata.obs[
-        ["sample_id", "tissue", "mouse_id"]
-    ].copy()
+    df_obs = adata.obs[["sample_id", "tissue", "mouse_id"]].copy()
     sample_annotation_path.parent.mkdir(
-        parents=True, exist_ok=True,
+        parents=True,
+        exist_ok=True,
     )
     df_obs.to_csv(
         sample_annotation_path,
         sep=sep_sample,
         index=False,
-        lineterminator='\n',
+        lineterminator="\n",
     )