add cosmx ingestion

CHANGELOG.MD

@@ -1,5 +1,9 @@
 # openpipeline_qc 0.2.0
 
+## NEW FEATURES
+
+* Add support for CosMx data in the QC Reporting workflow (PR #12).
+
 ## MINOR CHANGES
 
 * `src/workflows/generate_qc_report`: Add min and max constraints to `--max_samples_per_report`, to limit the number of samples per report between 1 and 20 (PR #11).

resources_test_scripts/qc_sample_data.sh

@@ -147,6 +147,14 @@ viash run src/ingestion_qc/h5mu_to_qc_json/config.vsh.yaml --engine docker -- \
   --output "$OUT_DIR_SPATIAL"/xenium_dataset.json \
   --output_reporting_json "$OUT_DIR_SPATIAL"/xenium_report_structure.json
 
+viash run src/ingestion_qc/h5mu_to_qc_json/config.vsh.yaml --engine docker -- \
+  --input "$OUT_DIR_SPATIAL"/Lung5_Rep2_tiny.qc.h5mu \
+  --input "$OUT_DIR_SPATIAL"/Lung5_Rep2_tiny.qc.h5mu \
+  --ingestion_method cosmx \
+  --min_num_nonzero_vars 1 \
+  --output "$OUT_DIR_SPATIAL"/cosmx_dataset.json \
+  --output_reporting_json "$OUT_DIR_SPATIAL"/cosmx_report_structure.json
+
 # remove all state yaml files
 rm "$OUT_DIR"/*.yaml
 rm "$OUT_DIR_SPATIAL"/*.yaml

src/ingestion_qc/h5mu_to_qc_json/config.vsh.yaml

@@ -34,6 +34,7 @@ argument_groups:
         choices:
           - cellranger_multi
           - xenium
+          - cosmx
         description: Method that was used to ingest the data - this will define the structure of the report that is generated.
       - name: --obs_sample_id
         type: string
@@ -132,7 +133,41 @@ argument_groups:
         description: The key in the h5mu .obs field that contains the number of control codewords.
         default: control_codeword_counts  
 
-  # - name: Options for CosMx reports
+  - name: Options for CosMx reports
+    arguments:
+      - name: --obs_area
+        type: string
+        description: The key in the h5mu .obs field that contains the cell area.
+        default: Area
+      - name: --obs_aspect_ratio
+        type: string
+        description: The key in the h5mu .obs field that contains the aspect ratio.
+        default: AspectRatio
+      - name: --obs_dapi_stain
+        type: string
+        multiple: true
+        description: The keys in the h5mu .obs field that contains the DAPI stain intensity.
+        default: [Mean.DAPI, Max.DAPI]
+      - name: --obs_membrane_stain
+        type: string
+        multiple: true
+        description: The keys in the h5mu .obs field that contains the membrane stain intensity.
+        default: [Mean.MembraneStain, Max.MembraneStain]
+      - name: --obs_panck_stain
+        type: string
+        multiple: true
+        description: The keys in the h5mu .obs field that contains the PanCK stain intensity.
+        default: [Mean.PanCK, Max.PanCK]
+      - name: --obs_cd45_stain
+        type: string
+        multiple: true
+        description: The keys in the h5mu .obs field that contains the CD45 stain intensity.
+        default: [Mean.CD45, Max.CD45]
+      - name: --obs_cd3_stain
+        type: string
+        multiple: true
+        description: The keys in the h5mu .obs field that contains the CD3 stain intensity.
+        default: [Mean.CD3, Max.CD3]
 
 resources:
   - type: python_script

src/ingestion_qc/h5mu_to_qc_json/report_structure/cosmx.json

@@ -0,0 +1,151 @@
+{
+  "categories": [
+    {
+      "name": "Sample QC",
+      "key": "sample_summary_stats",
+      "additionalAxes": false,
+      "defaultFilters": []
+    },
+    {
+      "name": "Cell RNA QC",
+      "key": "cell_rna_stats",
+      "additionalAxes": true,
+      "defaultFilters": [
+        {
+          "type": "histogram",
+          "visualizationType": "histogram",
+          "field": "total_counts",
+          "label": "Total UMI per cell",
+          "description": "Total number of RNA molecules detected per cell. Low values typically indicate empty droplets or low-quality cells that should be filtered out.",
+          "cutoffMin": null,
+          "cutoffMax": null,
+          "zoomMax": null,
+          "nBins": 50,
+          "groupBy": "sample_id",
+          "yAxisType": "linear"
+        },
+        {
+          "type": "histogram",
+          "visualizationType": "histogram",
+          "field": "num_nonzero_vars",
+          "label": "Number of non-zero genes per cell",
+          "description": "Count of unique genes detected in each cell. Low gene counts often indicate poor-quality cells.",
+          "cutoffMin": null,
+          "cutoffMax": null,
+          "zoomMax": null,
+          "nBins": 50,
+          "groupBy": "sample_id",
+          "yAxisType": "linear"
+        },
+        {
+          "type": "histogram",
+          "visualizationType": "histogram",
+          "field": "fraction_mitochondrial",
+          "label": "Fraction UMI of mitochondrial genes per cell",
+          "description": "Proportion of cell's RNA from mitochondrial genes.",
+          "cutoffMin": null,
+          "cutoffMax": null,
+          "nBins": 50,
+          "groupBy": "sample_id",
+          "yAxisType": "linear"
+        },
+        {
+          "type": "histogram",
+          "visualizationType": "histogram",
+          "field": "fraction_ribosomal",
+          "label": "Fraction UMI of ribosomal genes per cell",
+          "description": "Proportion of cell's RNA from ribosomal protein genes. Extreme values may indicate stress responses or cell cycle abnormalities.",
+          "cutoffMin": null,
+          "cutoffMax": null,
+          "nBins": 50,
+          "groupBy": "sample_id",
+          "yAxisType": "linear"
+        },
+        {
+          "type": "histogram",
+          "visualizationType": "histogram",
+          "field": "Area",
+          "label": "Cell Area",
+          "description": "Area of the segmented cells.",
+          "cutoffMin": null,
+          "cutoffMax": null,
+          "nBins": 50,
+          "groupBy": "sample_id",
+          "yAxisType": "linear"
+        },
+        {
+          "type": "histogram",
+          "visualizationType": "histogram",
+          "field": "AspectRatio",
+          "label": "Aspect Ratio",
+          "description": "Aspect ratio of the segmented cells.",
+          "cutoffMin": null,
+          "cutoffMax": null,
+          "nBins": 50,
+          "groupBy": "sample_id",
+          "yAxisType": "linear"
+        },
+        {
+          "type": "histogram",
+          "visualizationType": "histogram",
+          "field": "Mean.DAPI",
+          "label": "Mean DAPI Intensity",
+          "description": "Mean intensity of DAPI stain.",
+          "cutoffMin": null,
+          "cutoffMax": null,
+          "nBins": 50,
+          "groupBy": "sample_id",
+          "yAxisType": "linear"
+        },
+        {
+          "type": "histogram",
+          "visualizationType": "histogram",
+          "field": "Mean.MembraneStain",
+          "label": "Mean Membrane Stain Intensity",
+          "description": "Mean intensity of Membrane stain.",
+          "cutoffMin": null,
+          "cutoffMax": null,
+          "nBins": 50,
+          "groupBy": "sample_id",
+          "yAxisType": "linear"
+        },
+        {
+          "type": "histogram",
+          "visualizationType": "histogram",
+          "field": "Mean.PanCK",
+          "label": "Mean PanCK Intensity",
+          "description": "Mean intensity of PanCK stain.",
+          "cutoffMin": null,
+          "cutoffMax": null,
+          "nBins": 50,
+          "groupBy": "sample_id",
+          "yAxisType": "linear"
+        },
+        {
+          "type": "histogram",
+          "visualizationType": "histogram",
+          "field": "Mean.CD45",
+          "label": "Mean CD45 Intensity",
+          "description": "Mean intensity of CD45 stain.",
+          "cutoffMin": null,
+          "cutoffMax": null,
+          "nBins": 50,
+          "groupBy": "sample_id",
+          "yAxisType": "linear"
+        },
+        {
+          "type": "histogram",
+          "visualizationType": "histogram",
+          "field": "Mean.CD3",
+          "label": "Mean CD3 Intensity",
+          "description": "Mean intensity of CD3 stain.",
+          "cutoffMin": null,
+          "cutoffMax": null,
+          "nBins": 50,
+          "groupBy": "sample_id",
+          "yAxisType": "linear"
+        }
+      ]
+    }
+  ]
+}

src/ingestion_qc/h5mu_to_qc_json/script.py

@@ -39,7 +39,14 @@
     "obs_x_coord": "x_coord",
     "obs_y_coord": "y_coord",
     "obs_control_probe_counts": "control_probe_counts",
-    "obs_control_codeword_counts": "control_codeword_counts"
+    "obs_control_codeword_counts": "control_codeword_counts",
+    "obs_area": "Area",
+    "obs_aspect_ratio": "AspectRatio",
+    "obs_dapi_stain": ["Mean.DAPI", "Max.DAPI"],
+    "obs_membrane_stain": ["Mean.MembraneStain", "Max.MembraneStain"],
+    "obs_panck_stain": ["Mean.PanCK", "Max.PanCK"],
+    "obs_cd45_stain": ["Mean.CD45", "Max.CD45"],
+    "obs_cd3_stain": ["Mean.CD3", "Max.CD3"],
 }
 meta = {
     "resources_dir": os.path.abspath("src/ingestion_qc/h5mu_to_qc_json"),
@@ -56,6 +63,11 @@
 logger = setup_logger()
 
 par["obs_cellbender"] = {} if not par["obs_cellbender"] else par["obs_cellbender"]
+par["obs_dapi_stain"] = [] if not par["obs_dapi_stain"] else par["obs_dapi_stain"]
+par["obs_membrane_stain"] = [] if not par["obs_membrane_stain"] else par["obs_membrane_stain"]
+par["obs_panck_stain"] = [] if not par["obs_panck_stain"] else par["obs_panck_stain"]
+par["obs_cd45_stain"] = [] if not par["obs_cd45_stain"] else par["obs_cd45_stain"]
+par["obs_cd3_stain"] = [] if not par["obs_cd3_stain"] else par["obs_cd3_stain"]
 
 
 def transform_df(df):
@@ -247,6 +259,60 @@ def generate_xenium_stats(mod_obs, sample_id, required_keys):
     return cell_rna_stats
 
 
+def format_cosmx_columns(mod_obs):
+
+    cosmx_obs_keys = [par["obs_area"], par["obs_aspect_ratio"]]
+
+    # Add multiple keys
+    for key_list in [par["obs_dapi_stain"], par["obs_membrane_stain"], par["obs_panck_stain"],
+                     par["obs_cd45_stain"], par["obs_cd3_stain"]]:
+                     if key_list:
+                         cosmx_obs_keys.extend(key_list)
+
+    # Check if keys are present
+    check_optional_obs_keys(mod_obs, cosmx_obs_keys, "CosMx metrics.")
+
+    # Filter only existing keys
+    cosmx_formatted_columns = [key for key in cosmx_obs_keys if key in mod_obs]
+
+    # Add coordinates
+    cosmx_formatted_columns.extend(["x_coord", "y_coord"])
+
+    for key in cosmx_formatted_columns:
+        if key in mod_obs:
+            if not pd.api.types.is_float_dtype(mod_obs[key]):
+                try:
+                    mod_obs[key] = mod_obs[key].astype("float16")
+                except ValueError:
+                    logger.warning(f"Could not convert column {key} to a float dtype.")
+
+    return mod_obs, cosmx_formatted_columns
+
+
+def generate_cosmx_stats(mod_obs, sample_id, required_keys):
+
+    # Format required columns
+    mod_obs = format_required_columns(required_keys, mod_obs)
+
+    # Format cosmx-specific columns
+    mod_obs, cosmx_formatted_columns = format_cosmx_columns(mod_obs)
+
+    # Fetch and format  all categorical columns for grouping
+    metadata_obs_keys, mod_obs = format_categorical_columns(mod_obs)
+
+    # Create cell RNA stats dataframe
+    cell_rna_stats = pd.DataFrame(
+        {
+            "sample_id": pd.Categorical(sample_id),
+            **{key: mod_obs[key] for key in required_keys},
+            **{key: mod_obs[key] for key in cosmx_formatted_columns},
+            **{key: mod_obs[key] for key in metadata_obs_keys}
+        }
+    )
+
+    return cell_rna_stats
+
+
 def concatenate_dataframes(dfs):
     '''Concatenates a list of dataframes into a single dataframe, preserving categorical columns.'''
     df = pd.concat(dfs, ignore_index=True)
@@ -286,7 +352,7 @@ def main(par):
         barcodes_original_count = mod_obs.shape[0]
 
         # Add coordinates to obs before filtering
-        if par["ingestion_method"] == "xenium":
+        if par["ingestion_method"] == "xenium" or par["ingestion_method"] == "cosmx":
             mod_obs["x_coord"] = mod_obsm["spatial"][:, 0]
             mod_obs["y_coord"] = mod_obsm["spatial"][:, 1]
 
@@ -340,6 +406,9 @@ def main(par):
 
         if par["ingestion_method"] == "xenium":
             cell_rna_stats = generate_xenium_stats(mod_obs, sample_id, required_keys)
+
+        if par["ingestion_method"] == "cosmx":
+            cell_rna_stats = generate_cosmx_stats(mod_obs, sample_id, required_keys)
 
         cell_stats_dfs.append(cell_rna_stats)
         sample_stats_dfs.append(sample_summary_stats)
@@ -374,7 +443,8 @@ def main(par):
 
     report_structures = {
         "cellranger_multi": os.path.join(meta["resources_dir"], "report_structure/cellranger.json"),
-        "xenium": os.path.join(meta["resources_dir"], "report_structure/xenium.json")
+        "xenium": os.path.join(meta["resources_dir"], "report_structure/xenium.json"),
+        "cosmx": os.path.join(meta["resources_dir"], "report_structure/cosmx.json")
     }
 
     logger.info(f"Writing output report structure json to {par['output_reporting_json']}")

src/ingestion_qc/h5mu_to_qc_json/test.py

@@ -140,5 +140,43 @@ def test_xenium_execution(run_component, tmp_path):
             assert {"name", "dtype", "data"}.issubset(col.keys())
 
 
+def test_cosmx_execution(run_component, tmp_path):
+    output_json_path = tmp_path / "output.json"
+    output_reporting_json_path = tmp_path / "output_reporting.json"
+
+    run_component(
+        [
+            "--input", meta["resources_dir"] + "/resources_test/spatial_qc_sample_data/Lung5_Rep2_tiny.qc.h5mu",
+            "--ingestion_method", "cosmx",
+            "--min_num_nonzero_vars", "1",
+            "--output", output_json_path,
+            "--output_reporting_json", output_reporting_json_path
+        ]
+    )
+
+    assert os.path.exists(output_json_path), "Output file was not created"
+
+    with open(output_json_path, "r") as f:
+        output_json_dict = json.load(f)
+
+    assert output_json_dict.keys() == {"cell_rna_stats", "sample_summary_stats"}
+    assert "metrics_cellranger_stats" not in output_json_dict.keys()
+
+    column_names_cell = [col["name"] for col in output_json_dict["cell_rna_stats"]["columns"]]
+    expected_column_names = [
+        "sample_id", "total_counts", "num_nonzero_vars",
+        "fraction_mitochondrial",  "fraction_ribosomal",
+        "Area", "AspectRatio",
+        "Mean.DAPI", "Mean.MembraneStain",
+        "x_coord", "y_coord", "fov"
+        ]
+    assert np.all([column in column_names_cell for column in expected_column_names])
+
+    for key in output_json_dict.keys():
+        assert output_json_dict[key].keys() == {"num_rows", "num_cols", "min_total_counts", "min_num_nonzero_vars", "columns"}
+        for col in output_json_dict[key]["columns"]:
+            assert {"name", "dtype", "data"}.issubset(col.keys())
+
+
 if __name__ == "__main__":
     sys.exit(pytest.main([__file__]))