from_cellranger_multi_to_h5mu: support converting cellranger filtered h5 to h5mu (#1170)

Author: dorien-er

CHANGELOG.md

@@ -6,6 +6,12 @@
 
 * `workflows/rna/rna_multisample`, `workflows/multiomics/process_samples`, `workflows/multiomics/process_batches`: the RNA scaling zero-center argument is now a regular `boolean` (default behaviour remains unaltered and is set explicitly to `true`) instead of `boolean_false` (PR #1216).
 
+## NEW FEATURES
+
+* `convert/from_cellranger_multi_to_h5mu`: add `--output_filtered_data` flag to convert the per-sample filtered count matrices instead of the aggregated raw count matrix (PR #1170).
+
+* `workflows/ingestion/cellranger_multi`: surface the `--output_filtered_data` flag to convert the per-sample filtered count matrices instead of the aggregated raw count matrix (PR #1170).
+
 # openpipelines 4.2.0
 
 ## NEW FEATURES

src/convert/from_cellranger_multi_to_h5mu/config.vsh.yaml

@@ -46,6 +46,13 @@ arguments:
     type: string
     description: Name of the .uns slot under which to QC metrics (if any).
     default: "metrics_cellranger"
+  - name: "--output_filtered_data"
+    type: boolean_true
+    description: |
+      If enabled, read the per-sample filtered count matrices
+      (per_sample_outs/{sample}/count/sample_filtered_feature_bc_matrix.h5)
+      instead of the aggregated raw count matrix
+      (multi/count/raw_feature_bc_matrix.h5).
 __merge__: [., /src/base/h5_compression_argument.yaml]
 
 resources:

src/convert/from_cellranger_multi_to_h5mu/script.py

@@ -20,6 +20,7 @@
     "uns_metrics": "metrics_cellranger",
     "output_compression": "gzip",
     "sample_csv": "samples.csv",
+    "output_filtered_data": False,
 }
 meta = {"resources_dir": "./src/utils/"}
 ## VIASH END
@@ -101,6 +102,7 @@ def gather_input_data(dir: Path):
     #           |   +-- per_barcode.csv
     #           |   +-- antigen_specificity_scores.csv
     #           +-- count
+    #           |   +-- sample_filtered_feature_bc_matrix.h5
     #           |   +-- antibody_analysis
     #           |   +-- crispr_analysis
     #           |       +-- perturbation_efficiencies_by_feature.csv
@@ -174,6 +176,26 @@ def gather_input_data(dir: Path):
                 file_name = found_file.name.removesuffix(".csv")
                 found_input.setdefault(file_name, {})[samples_dir.name] = found_file
 
+    if par["output_filtered_data"]:
+        for samples_dir in samples_dirs:
+            for file_part in (
+                "count/sample_filtered_feature_bc_matrix.h5",
+                "sample_filtered_feature_bc_matrix.h5",  # Cell Ranger v10
+            ):
+                found_file = samples_dir / file_part
+                if found_file.exists():
+                    found_input.setdefault("filtered_counts", {})[samples_dir.name] = (
+                        found_file
+                    )
+                    break
+            else:
+                raise ValueError(
+                    f"Expected a filtered count matrix under {samples_dir}, "
+                    "but none was found. Make sure the input directory is a "
+                    "valid cellranger multi output that contains per-sample "
+                    "filtered feature-barcode matrices."
+                )
+
     return found_input
 
 
@@ -219,22 +241,21 @@ def process_feature_reference(
     return mudatas
 
 
-def process_counts(counts_folder: Path, multiplexing_info, metrics_files):
-    counts_matrix_file = counts_folder / "raw_feature_bc_matrix.h5"
-    logger.info("Reading %s.", counts_matrix_file)
-    adata = scanpy.read_10x_h5(counts_matrix_file, gex_only=False)
+def _modality_name_factory(library_type):
+    return ("".join(library_type.replace("-", "_").split())).lower()
 
-    # set the gene ids as var_names
-    logger.info("Renaming var columns")
-    adata.var = adata.var.rename_axis("gene_symbol").reset_index().set_index("gene_ids")
 
-    # generate output
-    logger.info("Convert to mudata")
+def _rename_var_to_gene_ids(adata: anndata.AnnData):
+    # set the gene ids as var_names (unique Ensembl IDs); gene symbols, which
+    # scanpy.read_10x_h5 uses as var_names by default, are not unique
+    adata.var = adata.var.rename_axis("gene_symbol").reset_index().set_index("gene_ids")
 
-    def modality_name_factory(library_type):
-        return ("".join(library_type.replace("-", "_").split())).lower()
 
-    feature_types = defaultdict(modality_name_factory, FEATURE_TYPES_NAMES)
+def _aggregated_counts_to_per_sample_mudatas(
+    adata: anndata.AnnData, multiplexing_info, metrics_files
+):
+    logger.info("Convert to mudata")
+    feature_types = defaultdict(_modality_name_factory, FEATURE_TYPES_NAMES)
     mudata_all_samples = mudata.MuData(adata, feature_types_names=feature_types)
     if multiplexing_info:
         # Get the mapping between the barcode and the sample ID from one of the metrics files
@@ -267,6 +288,29 @@ def modality_name_factory(library_type):
     return {"run": mudata_all_samples}
 
 
+def process_counts_filtered(filtered_counts: dict[str, Path]):
+    # Unlike the raw matrix, per-sample filtered matrices are already
+    # demultiplexed by cellranger, so each h5 maps 1:1 to an output mudata.
+    feature_types = defaultdict(_modality_name_factory, FEATURE_TYPES_NAMES)
+    mudatas = {}
+    for sample_name, filtered_h5 in filtered_counts.items():
+        logger.info("Reading %s.", filtered_h5)
+        adata = scanpy.read_10x_h5(filtered_h5, gex_only=False)
+        _rename_var_to_gene_ids(adata)
+        mudatas[sample_name] = mudata.MuData(adata, feature_types_names=feature_types)
+    return mudatas
+
+
+def process_counts(counts_folder: Path, multiplexing_info, metrics_files):
+    counts_matrix_file = counts_folder / "raw_feature_bc_matrix.h5"
+    logger.info("Reading %s.", counts_matrix_file)
+    adata = scanpy.read_10x_h5(counts_matrix_file, gex_only=False)
+    _rename_var_to_gene_ids(adata)
+    return _aggregated_counts_to_per_sample_mudatas(
+        adata, multiplexing_info, metrics_files
+    )
+
+
 def split_samples(mudata_obj, multiplexing_analysis_folder, barcode_sample_mapping):
     result = {}
     cells_per_tag_file = multiplexing_analysis_folder / "cells_per_tag.json"
@@ -410,14 +454,19 @@ def get_modalities(input_data):
         ),
         "antigen_analysis": process_antigen_analysis,
     }
-    mudata_per_sample = process_counts(
-        input_data["count"],
-        input_data["multiplexing_analysis"],
-        input_data["metrics_summary"],
-    )
+    if input_data.get("filtered_counts"):
+        mudata_per_sample = process_counts_filtered(
+            input_data["filtered_counts"],
+        )
+    else:
+        mudata_per_sample = process_counts(
+            input_data["count"],
+            input_data["multiplexing_analysis"],
+            input_data["metrics_summary"],
+        )
     for modality_name, modality_data_path in input_data.items():
         if (
-            modality_name in ("count", "multiplexing_analysis")
+            modality_name in ("count", "multiplexing_analysis", "filtered_counts")
             or not modality_data_path
         ):
             continue

src/convert/from_cellranger_multi_to_h5mu/test.py

@@ -340,6 +340,85 @@ def test_no_feature_reference(run_component, tmp_path, input_4plex_dtc):
     assert list(converted_data.mod.keys()) == ["rna", "prot"]
 
 
+def test_cellranger_multi_output_filtered_data(run_component, tmp_path, input_anticmv):
+    output_dir = tmp_path / "converted"
+    output_path_template = output_dir / "*.h5mu"
+    samples_csv = tmp_path / "samples.csv"
+    run_component(
+        [
+            "--input",
+            input_anticmv,
+            "--output",
+            str(output_path_template),
+            "--output_compression",
+            "gzip",
+            "--sample_csv",
+            samples_csv,
+            "--output_filtered_data",
+        ]
+    )
+    assert output_dir.is_dir()
+    samples = [item for item in output_dir.iterdir() if item.is_file()]
+    assert len(samples) == 1
+    filtered_data = read_h5mu(samples[0])
+    assert list(filtered_data.mod.keys()) == ["rna", "prot", "vdj_t"]
+    assert filtered_data.mod["rna"].n_obs == 3798
+
+
+def test_cellranger_multi_output_filtered_data_multiplexed(
+    run_component, tmp_path, input_fixed_rna
+):
+    output_dir = tmp_path / "converted"
+    output_path_template = output_dir / "*.h5mu"
+    samples_csv = tmp_path / "samples.csv"
+    run_component(
+        [
+            "--input",
+            input_fixed_rna,
+            "--output",
+            str(output_path_template),
+            "--output_compression",
+            "gzip",
+            "--sample_csv",
+            samples_csv,
+            "--output_filtered_data",
+        ]
+    )
+    assert output_dir.is_dir()
+    samples = [item for item in output_dir.iterdir() if item.is_file()]
+    sample_names = {item.name.removesuffix(".h5mu") for item in samples}
+    assert sample_names == {
+        "Colorectal_BC3",
+        "Liver_BC1",
+        "Ovarian_BC2",
+        "Pancreas_BC4",
+    }
+    expected_n_obs_by_version = {
+        "v9": {
+            "Colorectal_BC3": 1910,
+            "Liver_BC1": 1779,
+            "Ovarian_BC2": 1941,
+            "Pancreas_BC4": 862,
+        },
+        "v10": {
+            "Colorectal_BC3": 1915,
+            "Liver_BC1": 1784,
+            "Ovarian_BC2": 1946,
+            "Pancreas_BC4": 866,
+        },
+    }
+    expected_n_obs = expected_n_obs_by_version[
+        "v10" if "_v10" in input_fixed_rna else "v9"
+    ]
+    actual_n_obs = {}
+    for output_path in samples:
+        converted_data = read_h5mu(output_path)
+        assert list(converted_data.mod.keys()) == ["rna", "prot"]
+        sample_name = output_path.name.removesuffix(".h5mu")
+        actual_n_obs[sample_name] = converted_data.mod["rna"].n_obs
+    assert actual_n_obs == expected_n_obs
+
+
 def test_vdj_no_cells(run_component, tmp_path, input_no_vdj_cells):
     """
     Test what happens when a VDJ analysis was performed by Cell Ranger,

src/workflows/ingestion/cellranger_multi/config.vsh.yaml

@@ -32,6 +32,13 @@ argument_groups:
         type: string
         description: Name of the .uns slot under which to QC metrics (if any).
         default: "metrics_cellranger"
+      - name: "--output_filtered_data"
+        type: boolean_true
+        description: |
+          If enabled, read the per-sample filtered count matrices
+          (per_sample_outs/{sample}/count/sample_filtered_feature_bc_matrix.h5)
+          instead of the aggregated raw count matrix
+          (multi/count/raw_feature_bc_matrix.h5).
 dependencies:
   - name: mapping/cellranger_multi
     alias: cellranger_multi_component

src/workflows/ingestion/cellranger_multi/main.nf

@@ -93,6 +93,7 @@ workflow run_wf {
         [
           "input": state.input,
           "uns_metrics": state.uns_metrics,
+          "output_filtered_data": state.output_filtered_data,
           "output_compression": "gzip"
         ]
       },

src/workflows/ingestion/cellranger_multi/test.nf

@@ -9,32 +9,44 @@ workflow test_wf {
 
   resources_test = file(params.resources_test)
 
-  output_ch = Channel.fromList([
-      [
-        id: "foo",
-        input:[
-          resources_test.resolve("10x_5k_anticmv/raw/5k_human_antiCMV_T_TBNK_connect_GEX_1_subset_S1_L001_R1_001.fastq.gz"),
-          resources_test.resolve("10x_5k_anticmv/raw/5k_human_antiCMV_T_TBNK_connect_GEX_1_subset_S1_L001_R2_001.fastq.gz"),
-          resources_test.resolve("10x_5k_anticmv/raw/5k_human_antiCMV_T_TBNK_connect_AB_subset_S2_L004_R1_001.fastq.gz"),
-          resources_test.resolve("10x_5k_anticmv/raw/5k_human_antiCMV_T_TBNK_connect_AB_subset_S2_L004_R2_001.fastq.gz"),
-          resources_test.resolve("10x_5k_anticmv/raw/5k_human_antiCMV_T_TBNK_connect_VDJ_subset_S1_L001_R1_001.fastq.gz"),
-          resources_test.resolve("10x_5k_anticmv/raw/5k_human_antiCMV_T_TBNK_connect_VDJ_subset_S1_L001_R2_001.fastq.gz")
-        ],
-        gex_reference: resources_test.resolve("reference_gencodev41_chr1/reference_cellranger.tar.gz"),
-        vdj_reference: resources_test.resolve("10x_5k_anticmv/raw/refdata-cellranger-vdj-GRCh38-alts-ensembl-7.0.0.tar.gz"),
-        feature_reference: resources_test.resolve("10x_5k_anticmv/raw/feature_reference.csv"),
-        library_id: [
-          "5k_human_antiCMV_T_TBNK_connect_GEX_1_subset",
-          "5k_human_antiCMV_T_TBNK_connect_AB_subset",
-          "5k_human_antiCMV_T_TBNK_connect_VDJ_subset"
-        ],
-        library_type: [
-          "Gene Expression",
-          "Antibody Capture",
-          "VDJ"
-        ]
-      ]
-    ])
+  base_state = [
+    id: "foo",
+    input:[
+      resources_test.resolve("10x_5k_anticmv/raw/5k_human_antiCMV_T_TBNK_connect_GEX_1_subset_S1_L001_R1_001.fastq.gz"),
+      resources_test.resolve("10x_5k_anticmv/raw/5k_human_antiCMV_T_TBNK_connect_GEX_1_subset_S1_L001_R2_001.fastq.gz"),
+      resources_test.resolve("10x_5k_anticmv/raw/5k_human_antiCMV_T_TBNK_connect_AB_subset_S2_L004_R1_001.fastq.gz"),
+      resources_test.resolve("10x_5k_anticmv/raw/5k_human_antiCMV_T_TBNK_connect_AB_subset_S2_L004_R2_001.fastq.gz"),
+      resources_test.resolve("10x_5k_anticmv/raw/5k_human_antiCMV_T_TBNK_connect_VDJ_subset_S1_L001_R1_001.fastq.gz"),
+      resources_test.resolve("10x_5k_anticmv/raw/5k_human_antiCMV_T_TBNK_connect_VDJ_subset_S1_L001_R2_001.fastq.gz")
+    ],
+    gex_reference: resources_test.resolve("reference_gencodev41_chr1/reference_cellranger.tar.gz"),
+    vdj_reference: resources_test.resolve("10x_5k_anticmv/raw/refdata-cellranger-vdj-GRCh38-alts-ensembl-7.0.0.tar.gz"),
+    feature_reference: resources_test.resolve("10x_5k_anticmv/raw/feature_reference.csv"),
+    library_id: [
+      "5k_human_antiCMV_T_TBNK_connect_GEX_1_subset",
+      "5k_human_antiCMV_T_TBNK_connect_AB_subset",
+      "5k_human_antiCMV_T_TBNK_connect_VDJ_subset"
+    ],
+    library_type: [
+      "Gene Expression",
+      "Antibody Capture",
+      "VDJ"
+    ]
+  ]
+
+  // Convert the aggregated raw count matrix (default behaviour).
+  raw_ch = Channel.fromList([base_state])
+    | map{ state -> [state.id, state] }
+    | cellranger_multi
+    | view { output ->
+      assert output.size() == 2 : "outputs should contain two elements; [id, out]"
+      assert output[1] instanceof Map : "Output should be a Map."
+      // todo: check whether output dir contains fastq files
+      "Output: $output"
+    }
+
+  // Convert the per-sample filtered count matrices instead.
+  filtered_ch = Channel.fromList([base_state + [output_filtered_data: true]])
     | map{ state -> [state.id, state] }
     | cellranger_multi
     | view { output ->
@@ -44,8 +56,13 @@ workflow test_wf {
       "Output: $output"
     }
 
+  raw_ch
+    | join(filtered_ch)
+    | map { id, raw_state, filtered_state ->
+      [id, ["input": raw_state.output_h5mu, "input_filtered": filtered_state.output_h5mu]]
+    }
     | cellranger_multi_test.run(
-      fromState: ["input": "output_h5mu"]
+      fromState: ["input": "input", "input_filtered": "input_filtered"]
     )
 
     | toSortedList()

src/workflows/test_workflows/ingestion/cellranger_multi/config.vsh.yaml

@@ -13,6 +13,15 @@ argument_groups:
         multiple: true
         description: Path to h5mu output.
         example: foo.final.h5mu
+      - name: "--input_filtered"
+        type: file
+        required: false
+        multiple: true
+        description: |
+          Path to h5mu output generated with --output_filtered_data enabled.
+          When provided, the per-sample filtered count matrices are checked to
+          contain fewer cells than the corresponding raw count matrices.
+        example: foo.filtered.h5mu
 resources:
   - type: python_script
     path: script.py

src/workflows/test_workflows/ingestion/cellranger_multi/script.py

@@ -3,7 +3,7 @@
 import pytest
 
 ##VIASH START
-par = {"input": "input.h5mu"}
+par = {"input": "input.h5mu", "input_filtered": None}
 
 meta = {"resources_dir": "resources_test"}
 ##VIASH END
@@ -28,5 +28,18 @@ def test_run():
         )
 
 
+def test_filtered_data_has_fewer_cells():
+    if not par.get("input_filtered"):
+        pytest.skip("No filtered input provided.")
+    for raw_path, filtered_path in zip(par["input"], par["input_filtered"]):
+        raw_mudata = read_h5mu(raw_path)
+        filtered_mudata = read_h5mu(filtered_path)
+
+        assert filtered_mudata.mod["rna"].n_obs < raw_mudata.mod["rna"].n_obs, (
+            "Expected the filtered count matrix to contain fewer cells than the "
+            "raw count matrix when --output_filtered_data is enabled."
+        )
+
+
 if __name__ == "__main__":
     sys.exit(pytest.main([__file__, "--import-mode=importlib"]))