Merge pull request #3 from czbiohub-sf/marika

Adding filter_min_consecutive_fractions function and updating documentation.

Author: mffrank

docs/source/api/preprocessing.md

@@ -28,6 +28,7 @@ Any transformation of the data matrix that is not a *tool*. Other than *tools*,
    pp.normalize_total
    pp.drop_excess_MQ_metadata
    pp.remove_contaminants
+   pp.filter_min_consecutive_fractions
 ```
 
 ### Imputation

docs/source/conf.py

@@ -1,4 +1,3 @@
-import os
 import sys
 from pathlib import Path
 

docs/source/generated/grassp.pp.filter_min_consecutive_fractions.rst

@@ -0,0 +1,6 @@
+grassp.pp.filter\_min\_consecutive\_fractions
+=============================================
+
+.. currentmodule:: grassp.pp
+
+.. autofunction:: filter_min_consecutive_fractions

docs/source/installation.md

@@ -6,7 +6,7 @@ We recommend using a conda environment to install the package.
 ```{code-block} bash
 git clone https://github.com/czbiohub-sf/grassp.git
 cd grassp
-conda create -n grassp python=3.10
+conda create -n grassp python=3.12
 conda activate grassp
 pip install -e .
 ```

grassp/preprocessing/__init__.py

@@ -11,4 +11,5 @@
     highly_variable_proteins,
     normalize_total,
     remove_contaminants,
+    filter_min_consecutive_fractions,
 )

grassp/preprocessing/enrichment.py

@@ -2,7 +2,7 @@
 from typing import TYPE_CHECKING
 
 if TYPE_CHECKING:
-    from typing import Optional
+    from typing import Optional, Literal
 
 import warnings
 
@@ -144,9 +144,11 @@ def calculate_enrichment_vs_all(
     adata: AnnData,
     covariates: Optional[list[str]] = None,
     subcellular_enrichment_column: str = "subcellular_enrichment",
+    enrichment_method: Literal["lfc", "proportion"] = "lfc",
     correlation_threshold: float = 1.0,
     original_intensities_key: str | None = "original_intensities",
     keep_raw: bool = True,
+    min_comparison_warning: int | None = None,
 ) -> AnnData:
     """Calculate enrichment of each subcellular enrichment vs all other samples as the background.
 
@@ -159,10 +161,16 @@ def calculate_enrichment_vs_all(
         If None, uses columns starting with "covariate_"
     subcellular_enrichment_column
         Column in adata.var containing subcellular enrichment labels
+    enrichment_method
+        Calculating enrichment based on Log Fold Change (lfc) or Proportion-based analysis.
+        Must be either "proportion" or "lfc"
     original_intensities_key
         If provided, store the original intensities in this layer
     keep_raw
         Whether to keep the unaggregated data in the .raw attribute of the returned AnnData object
+    min_comparison_warning
+        The minimum number of control samples required before issuing a warning about low statistical power.
+
 
     Returns
     -------
@@ -171,27 +179,31 @@ def calculate_enrichment_vs_all(
     Raw values are stored in .layers[original_intensities_key] if provided.
     """
 
+    if enrichment_method not in ["lfc", "proportion"]:
+        raise ValueError("enrichment_method must be either 'lfc' or 'proportion'")
+
     data = adata.copy()
 
     if covariates is None:
         covariates = data.var.columns[data.var.columns.str.startswith("covariate_")].tolist()
-    # Check that all covariates are in the data
+    if not isinstance(covariates, list):
+        covariates = [covariates]
+
     for c in covariates:
         if c not in data.var.columns:
             raise ValueError(f"Covariate {c} not found in data.var.columns")
 
-    if not isinstance(covariates, list):
-        covariates = [covariates]
-    # Create aggregated data with the desired output shape
     grouping_columns = [subcellular_enrichment_column] + covariates
-    # Create a temporary column that contains the experimental conditions
+
     data.var["_experimental_condition"] = data.var[grouping_columns].apply(
         lambda x: "_".join(x.dropna().astype(str)),
         axis=1,
     )
+    data.var["_covariates"] = data.var[covariates].apply(
+        lambda x: "_".join(x.dropna().astype(str)), axis=1
+    )
 
     data_aggr = aggregate_samples(data, grouping_columns=grouping_columns, keep_raw=False)
-    data_aggr.var_names = data_aggr.var_names.str.replace(r"_\d+", "", regex=True)
 
     if original_intensities_key is not None:
         data_aggr.layers[original_intensities_key] = data_aggr.X
@@ -203,22 +215,40 @@ def calculate_enrichment_vs_all(
 
     for experimental_condition in data_aggr.var["_experimental_condition"].unique():
         mask = data_aggr.var["_experimental_condition"] == experimental_condition
-        corr_mat_sub = corr_matrix[mask, :].mean(axis=0)
-        control_mask = ~mask & (corr_mat_sub < correlation_threshold)
-        if control_mask.sum() < 10:
-            warnings.warn(
-                f"Less than 10 ({control_mask.sum()}) control samples found for condition: {experimental_condition}"
-            )
-        intensities_control = intensities[:, control_mask]
+
         intensities_ip = intensities[:, mask]
+        covariate = data.var.loc[
+            data.var._experimental_condition == experimental_condition, "_covariates"
+        ].values[0]
+        covariate_mask = data_aggr.var["_covariates"] == covariate
+        control_mask = ~mask & covariate_mask
+        corr_mat_sub = corr_matrix[mask, control_mask].mean(axis=0)
+        control_mask = control_mask & (corr_mat_sub < correlation_threshold)
+        intensities_control = intensities[:, control_mask]
+        if min_comparison_warning is not None:
+            if control_mask.sum() < min_comparison_warning:
+                warnings.warn(
+                    f"Less than {min_comparison_warning} ({control_mask.sum()}) control samples found for condition: {experimental_condition}"
+                )  # Check for statistical power (if fewer than 10 samples selected )
+
         scores, pv = stats.ttest_ind(intensities_ip.T, intensities_control.T)
-        lfc = np.median(intensities_ip, axis=1) - np.median(intensities_control, axis=1)
+
+        if enrichment_method == "lfc":
+            enrichment_values = np.median(intensities_ip, axis=1) - np.median(
+                intensities_control, axis=1
+            )
+        else:
+            enrichment_values = np.nansum(intensities_ip, axis=1) / (
+                np.nansum(intensities_ip, axis=1) + np.nansum(intensities_control, axis=1)
+            )
+
         aggr_mask = data_aggr.var["_experimental_condition"] == experimental_condition
         data_aggr.layers["pvals"][:, aggr_mask] = pv[:, None]
-        data_aggr[:, aggr_mask].X = lfc[:, None]
+        data_aggr.X[:, aggr_mask] = enrichment_values[:, None]
         data_aggr.var.loc[aggr_mask, "enriched_vs"] = ",".join(
             data_aggr.var_names[control_mask]
         )
+
     data_aggr.var.drop(columns=["_experimental_condition"], inplace=True)
     if keep_raw:
         data_aggr.raw = data.copy()