add visium support

CHANGELOG.MD

@@ -2,6 +2,12 @@
 
 ## MINOR CHANGES
 
+* Added support for Visium data ingestion and QC reporting (PR #16).
+
+# openpipeline_qc 0.2.0
+
+## 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).
 
 * Update `openpipeline` dependency to v4.0.0 and `siqc` to version v0.2.0 (PR #15).

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"/visium_tiny.qc.h5mu \
+  --input "$OUT_DIR_SPATIAL"/visium_tiny.qc.h5mu \
+  --ingestion_method visium \
+  --min_num_nonzero_vars 1 \
+  --output "$OUT_DIR_SPATIAL"/visium_dataset.json \
+  --output_reporting_json "$OUT_DIR_SPATIAL"/visium_report_structure.json
+
 # remove all state yaml files
 rm "$OUT_DIR"/*.yaml
 rm "$OUT_DIR_SPATIAL"/*.yaml

resources_test_scripts/spatial_qc_sample_data.sh

@@ -11,6 +11,8 @@ param_list:
     input: s3://openpipelines-bio/openpipeline_spatial/resources_test/xenium/xenium_tiny.h5mu
   - id: Lung5_Rep2_tiny
     input: s3://openpipelines-bio/openpipeline_spatial/resources_test/cosmx/Lung5_Rep2_tiny.h5mu
+  - id: visium_tiny
+    input: s3://openpipelines-bio/openpipeline_spatial/resources_test/visium/Visium_FFPE_Human_Ovarian_Cancer_tiny.h5mu
 var_name_mitochondrial_genes: mitochondrial
 var_name_ribosomal_genes: ribosomal
 output: '\$id.qc.h5mu'

src/ingestion_qc/generate_html/config.vsh.yaml

@@ -52,7 +52,7 @@ engines:
           - git
       - type: docker
         run: |
-          cd /opt && git clone -b v0.2.0 https://github.com/openpipelines-bio/siqc.git \
+          cd /opt && git clone -b v0.2.1 https://github.com/openpipelines-bio/siqc.git \
           && cd siqc && npm install \
           && true
 

src/ingestion_qc/generate_html/script.sh

@@ -5,4 +5,4 @@ ABSOLUTE_OUTPUT=$(realpath $par_output_qc_report)
 cd /opt/siqc
 mkdir src/data
 
-npm run cli render -- --data "$ABSOLUTE_INPUT_DATA" --structure "$ABSOLUTE_INPUT_STRUCTURE" --output "$ABSOLUTE_OUTPUT"
+npm run cli render -- --data "$ABSOLUTE_INPUT_DATA" --structure "$ABSOLUTE_INPUT_STRUCTURE" --output "$ABSOLUTE_OUTPUT"

src/ingestion_qc/h5mu_to_qc_json/config.vsh.yaml

@@ -34,6 +34,7 @@ argument_groups:
         choices:
           - cellranger_multi
           - xenium
+          - visium
         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

src/ingestion_qc/h5mu_to_qc_json/report_structure/visium.json

@@ -0,0 +1,67 @@
+{
+  "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"
+        }
+      ]
+    }
+  ]
+}

src/ingestion_qc/h5mu_to_qc_json/script.py

@@ -243,6 +243,31 @@ def generate_xenium_stats(mod_obs, sample_id, required_keys):
             **{key: mod_obs[key] for key in metadata_obs_keys}
         }
     )
+
+    return cell_rna_stats
+
+def generate_visium_stats(mod_obs, sample_id, required_keys):
+
+    # Format required columns
+    mod_obs = format_required_columns(required_keys, mod_obs)
+
+    # Format visium-specific columns
+    visium_formatted_columns = ["x_coord", "y_coord"]
+    for key in visium_formatted_columns:
+        mod_obs[key] = mod_obs[key].astype("float16")
+
+    # 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 visium_formatted_columns},
+            **{key: mod_obs[key] for key in metadata_obs_keys}
+        }
+    )
 
     return cell_rna_stats
 
@@ -286,7 +311,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"] == "visium":
             mod_obs["x_coord"] = mod_obsm["spatial"][:, 0]
             mod_obs["y_coord"] = mod_obsm["spatial"][:, 1]
 
@@ -340,6 +365,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"] == "visium":
+            cell_rna_stats = generate_visium_stats(mod_obs, sample_id, required_keys)
 
         cell_stats_dfs.append(cell_rna_stats)
         sample_stats_dfs.append(sample_summary_stats)
@@ -374,7 +402,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"),
+        "visium": os.path.join(meta["resources_dir"], "report_structure/visium.json")
     }
 
     logger.info(f"Writing output report structure json to {par['output_reporting_json']}")

src/workflows/generate_qc_report/config.vsh.yaml

@@ -35,6 +35,7 @@ argument_groups:
         choices:
           - cellranger_multi
           - xenium
+          - visium
       - name: --sample_metadata
         type: file
         required: false
@@ -131,6 +132,12 @@ test_resources:
   - type: nextflow_script
     path: test.nf
     entrypoint: test_with_cellbender
+  - type: nextflow_script
+    path: test.nf
+    entrypoint: test_xenium
+  - type: nextflow_script
+    path: test.nf
+    entrypoint: test_visium
 
 dependencies:
   - name: metadata/add_id

src/workflows/generate_qc_report/integration_test.sh

@@ -35,3 +35,10 @@ nextflow run . \
   -entry test_multiple_reports \
   -c src/configs/labels_ci.config \
   -resume
+
+nextflow run . \
+  -main-script src/workflows/generate_qc_report/test.nf \
+  -profile docker,no_publish,local \
+  -entry test_visium \
+  -c src/configs/labels_ci.config \
+  -resume

src/workflows/generate_qc_report/test.nf

@@ -221,4 +221,48 @@ workflow test_multiple_reports {
         assert files.size() == 5 : "Output directory should contain exactly 5 files, but found ${files.size()} files"
         "Output: $output"
     }
-}
+}
+workflow test_visium {
+
+  resources_test_file = file(params.resources_test)
+
+  output_ch = Channel.fromList([
+      [
+        id: "sample_one",
+        input: resources_test_file.resolve("spatial_qc_sample_data/visium_tiny.qc.h5mu"),
+        run_cellbender: false,
+        ingestion_method: "visium",
+        min_num_nonzero_vars: "1",
+        output_html: "report.html",
+        publish_dir: "test_out"
+      ],
+      [
+        id: "sample_two",
+        input: resources_test_file.resolve("spatial_qc_sample_data/visium_tiny.qc.h5mu"),
+        ingestion_method: "visium",
+        min_num_nonzero_vars: "1",
+        run_cellbender: false,
+        output_html: "report.html",
+        publish_dir: "test_out"
+      ]
+    ])
+
+    | map{ state -> [state.id, state] }
+    | generate_qc_report
+
+    | view { output ->
+        assert output.size() == 2 : "Outputs should contain two elements; [id, state]"
+        def id = output[0]
+        def state = output [1]
+        assert id == "combined": "Output ID should be `combined`"
+        assert state instanceof Map : "State should be a map. Found: ${state}"
+        assert state.containsKey("output_qc_report"): "Output should contain key `output_qc_report`"
+        assert state.containsKey("output_processed_h5mu"): "Output should contain key `output_processed_h5mu`"
+        assert state.output_qc_report.size() == 1 : "Expected exactly one output HTML file to be generated"
+        assert state.output_qc_report.every { it.isFile()} : "All output HTML report files should exist"
+        assert state.output_processed_h5mu.isDirectory() : "Output directory should exist"
+        def files = state.output_processed_h5mu.listFiles().findAll { it.isFile() }
+        assert files.size() == 2 : "Output directory should contain exactly 2 files, but found ${files.size()} files"
+        "Output: $output"
+    }
+}