optional QC filtering module

aegle-docs/docs/Main/Cell_filter_qc.md

@@ -0,0 +1,64 @@
+# Optional Cell Level QC Filtering
+
+## Overview
+During the cell profiling stage of the pipeline, various morphological metrics (such as eccentricity, convex area, area, and solidity) are calculated for both the cell and its nucleus. **In addition, the pipeline calculates a "retained fraction" (e.g., `nucleus_retained_fraction`, `cell_retained_fraction`), representing the proportion of original segmentation pixels kept after the mask repair and matching stage.**
+
+The **Optional Cell Level QC Filter** module (`aegle/qc_filter_cells.py`) utilizes these specific metrics to systematically identify and remove anomalous cells. This helps to clean up segmentation artifacts such as:
+* Cells that are unrealistically small or large.
+* Highly irregular or fragmented cells (low solidity).
+* Obvious artifacts (like perfectly circular false positives).
+* **Heavily trimmed nuclei or mismatched cells** that lost too much of their original area during the repair stage (i.e. nucleus leaked outside the cell membrane).
+
+By filtering out these low-quality segmentations, you can significantly improve the reliability of downstream spatial and single-cell expression analyses.
+
+## How it works
+The filter acts as a post-processing step. It takes the un-filtered output CSVs from the `cell_profiling` directory (specifically `cell_metadata.csv` and `cell_by_marker.csv`) and applies a set of configurable threshold boundaries to whichever metrics you specify.
+
+Because the filtering logic dynamically reads the columns in `cell_metadata.csv`, you can filter on literally *any* column present in that file (including the newly added `nucleus_retained_fraction`). Cells that meet all criteria are kept, and the module outputs new, leaner CSV files prefixed with `filtered_`.
+
+## Usage
+
+### 1. Standalone Script (Recommended)
+Because tuning QC thresholds often requires iterative trial and error, it is highly recommended to run this filter as a standalone script on the outputs of a completed pipeline run. This way, you do not need to re-run the heavy segmentation or feature extraction steps.
+
+```bash
+python aegle/qc_filter_cells.py \
+  --input_dir out/your_run/cell_profiling \
+  --config exps/configs/your_filter_config.yaml
+```
+
+* **`--input_dir`**: The directory containing `cell_metadata.csv` and `cell_by_marker.csv`.
+* **`--config`**: A YAML file containing your chosen filtering thresholds.
+
+### 2. Integration with `pipeline.py`
+If your QC thresholds are entirely finalized and you want the filtering to happen automatically at the end of a big batch-processing job, you can trigger the module directly in Python:
+
+```python
+from aegle.qc_filter_cells import apply_morphology_filters
+
+# In pipeline.py, assuming profiling_out_dir is defined
+if config.get("qc_filtering"):
+    apply_morphology_filters(profiling_out_dir, config["qc_filtering"])
+```
+
+## Configuration Reference
+Your YAML configuration file must nest its rules under the `qc_filtering` key. For each metric, you can specify a `min` bound, a `max` bound, or both.
+
+**Example `your_filter_config.yaml`:**
+```yaml
+qc_filtering:
+  rules:
+    nucleus_area:
+      min: 15       # Remove tiny artifacts
+      max: 1000     # Remove massive clumps of unresolved nuclei
+    nucleus_solidity:
+      min: 0.8      # Ensure the nucleus has a regular, non-fragmented shape
+    cell_eccentricity:
+      max: 0.95     # Remove perfectly linear streaks or artifacts
+    nucleus_retained_fraction:
+      min: 0.7      # Remove nuclei where more than 30% of pixels were trimmed during repair
+    cell_retained_fraction:
+      min: 0.8      # Remove cells heavily modified during mask alignment
+```
+
+*Note: The keys under `rules` must exactly match the morphology column names found in your `cell_metadata.csv`.*

aegle/cell_profiling.py

@@ -171,6 +171,40 @@ def run_cell_profiling(codex_patches, config, args):
             # Clean up image_dict after CPU extraction
             del image_dict
         logger.info(f"Extracted features for patch {patch_idx} with shape: {exp_df.shape}")
+
+        # Calculate per-cell retained fraction after segmentation repair
+        try:
+            orig_seg_result = codex_patches.original_seg_res_batch[seg_batch_idx]
+            orig_nucleus_mask = orig_seg_result.get("nucleus")
+            orig_cell_mask = orig_seg_result.get("cell")
+
+            cell_ids = metadata_df.index.values.astype(int)
+            orig_nuc_areas = np.bincount(orig_nucleus_mask.ravel()) if orig_nucleus_mask is not None else np.array([])
+            orig_cell_areas = np.bincount(orig_cell_mask.ravel()) if orig_cell_mask is not None else np.array([])
+
+            def map_areas(area_array, ids):
+                out = np.zeros(len(ids), dtype=float)
+                valid_mask = ids < len(area_array)
+                out[valid_mask] = area_array[ids[valid_mask]]
+                return out
+            # add to metadata_df
+            metadata_df["nucleus_area_before"] = map_areas(orig_nuc_areas, cell_ids)
+            metadata_df["nucleus_retained_fraction"] = (
+                metadata_df["nucleus_area"] / metadata_df["nucleus_area_before"].clip(lower=1e-6)
+            ).fillna(0)
+
+            metadata_df["cell_area_before"] = map_areas(orig_cell_areas, cell_ids)
+            metadata_df["cell_retained_fraction"] = (
+                metadata_df.get("cell_area", metadata_df.get("area", 0)) / 
+                metadata_df["cell_area_before"].clip(lower=1e-6)
+            ).fillna(0)
+
+            # Optional flags
+            metadata_df["nucleus_was_trimmed"] = metadata_df["nucleus_retained_fraction"] < 0.99
+
+        except Exception as e:
+            logger.warning(f"Failed to calculate cell-level repair QC metrics: {e}")
+
         logger.info(f"Metadata for patch {patch_idx} with shape: {metadata_df.shape}")
         logger.info(f"Exp DataFrame: {exp_df.head()}")
         logger.info(f"Metadata DataFrame: {metadata_df.head()}")

aegle/qc_filter_cells.py

@@ -0,0 +1,100 @@
+import os
+import argparse
+import logging
+import pandas as pd
+import yaml
+from typing import Dict, Any
+
+logger = logging.getLogger(__name__)
+
+def apply_morphology_filters(profiling_out_dir: str, qc_config: Dict[str, Any]):
+    """
+    Applies quality control filters to cell metadata and intensities.
+    Expects profiling_out_dir to contain 'cell_metadata.csv' and 'cell_by_marker.csv'.
+    """
+    metadata_path = os.path.join(profiling_out_dir, "cell_metadata.csv")
+    expression_path = os.path.join(profiling_out_dir, "cell_by_marker.csv")
+
+    if not os.path.exists(metadata_path):
+        logger.warning(f"Could not find {metadata_path}. Skipping QC filtering.")
+        return
+
+    logger.info(f"Loading metadata from {metadata_path}")
+    metadata_df = pd.read_csv(metadata_path)
+
+    original_count = len(metadata_df)
+    valid_mask = pd.Series(True, index=metadata_df.index)
+
+    rules = qc_config.get("rules", {})
+    if not rules:
+        logger.info("No QC rules defined in config. Returning original data.")
+        return
+
+    logger.info("Applying morphology QC rules...")
+    for column, bounds in rules.items():
+        if column not in metadata_df.columns:
+            logger.warning(f"Column '{column}' not found in metadata. Skipping rule.")
+            continue
+
+        if "min" in bounds:
+            min_val = bounds["min"]
+            kept = metadata_df[column] >= min_val
+            valid_mask &= kept
+            logger.info(f"Rule: {column} >= {min_val} (Kept {kept.sum()})")
+
+        if "max" in bounds:
+            max_val = bounds["max"]
+            kept = metadata_df[column] <= max_val
+            valid_mask &= kept
+            logger.info(f"Rule: {column} <= {max_val} (Kept {kept.sum()})")
+
+    passed_count = valid_mask.sum()
+    logger.info(f"QC Filtering Complete: {passed_count}/{original_count} cells passed ({(passed_count/original_count)*100:.2f}%).")
+
+    filtered_metadata = metadata_df[valid_mask]
+    filtered_meta_path = os.path.join(profiling_out_dir, "filtered_cell_metadata.csv")
+    filtered_metadata.to_csv(filtered_meta_path, index=False)
+    logger.info(f"Saved filtered metadata to {filtered_meta_path}")
+
+    # Filter and Save Expression Data if it exists
+    if os.path.exists(expression_path):
+        exp_df = pd.read_csv(expression_path)
+        if len(exp_df) == original_count:
+            filtered_exp = exp_df[valid_mask]
+            filtered_exp_path = os.path.join(profiling_out_dir, "filtered_cell_by_marker.csv")
+            filtered_exp.to_csv(filtered_exp_path, index=False)
+            logger.info(f"Saved filtered cell-by-marker to {filtered_exp_path}")
+        else:
+            logger.warning("Row count mismatch between metadata and expression data. Cannot safely filter expression data.")
+
+    # Also attempt to filter overview csvs
+    for overview_name in ["cell_overview.csv", "nucleus_overview.csv"]:
+        overview_path = os.path.join(profiling_out_dir, overview_name)
+        if os.path.exists(overview_path):
+            overview_df = pd.read_csv(overview_path)
+            if len(overview_df) == original_count:
+                filtered_overview = overview_df[valid_mask]
+                filtered_overview_path = os.path.join(profiling_out_dir, f"filtered_{overview_name}")
+                filtered_overview.to_csv(filtered_overview_path, index=False)
+                logger.info(f"Saved filtered overview to {filtered_overview_path}")
+
+
+if __name__ == "__main__":
+    # Ensure standard logging when run as standalone script
+    logging.basicConfig(level=logging.INFO, format="%(asctime)s - %(levelname)s - %(message)s")
+
+    parser = argparse.ArgumentParser(description="Filter cell metadata and marker profiles based on QC metrics (e.g. morphology, nucleus_retained_fraction).")
+    parser.add_argument("--input_dir", type=str, required=True, 
+                        help="Directory containing cell_metadata.csv and cell_by_marker.csv (Usually out_dir/cell_profiling)")
+    parser.add_argument("--config", type=str, required=True, 
+                        help="Path to YAML configuration file with filtering rules")
+
+    args = parser.parse_args()
+
+    with open(args.config, 'r') as f:
+        config = yaml.safe_load(f)
+
+    # The config will have the rules nested under 'qc_filtering'
+    qc_config = config.get("qc_filtering", config)
+
+    apply_morphology_filters(args.input_dir, qc_config)