diff --git a/aegle-docs/docs/Main/Cell_filter_qc.md b/aegle-docs/docs/Main/Cell_filter_qc.md new file mode 100644 index 0000000..4e9f3ce --- /dev/null +++ b/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`.* diff --git a/aegle/cell_profiling.py b/aegle/cell_profiling.py index a5070c9..e9f6be6 100644 --- a/aegle/cell_profiling.py +++ b/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()}") diff --git a/aegle/qc_filter_cells.py b/aegle/qc_filter_cells.py new file mode 100644 index 0000000..be9869c --- /dev/null +++ b/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)