ITADN

CELLBENDER_MERGE fails for 10x Flex data due to gene-set mismatch between filtered and raw matrices

#289Opennictru 创建于 2026-05-28
N
nictrucommented
## Summary `CELLBENDER_MERGE` crashes with a `ValueError` when the pipeline is run on data generated by the **10x Chromium Fixed RNA Profiling (Flex)** protocol. The root cause is that Cell Ranger outputs the filtered and raw matrices with **different gene sets** for Flex data, which the merge template does not account for. ## Error ``` ValueError: Value passed for key 'ambient_corrected_cellbender' is of incorrect shape. Values of layers must match dimensions (0, 1) of parent. Value had shape (9802, 33774) while it should have had (9802, 19070). ``` ## Root cause In standard 10x scRNA-seq (Cell Ranger count, STARsolo), the filtered and raw/unfiltered matrices always share the same gene set — they differ only in which barcodes are retained. With Cell Ranger Flex, this assumption breaks down: | Matrix | Barcodes | Genes | |--------|----------|-------| | `sample_filtered_feature_bc_matrix.h5` | called cells only | **probe-targeted genes only** (e.g. 19,070) | | `sample_raw_feature_bc_matrix.h5` | all droplets | **all reference genes** (e.g. 33,774) | Cell Ranger includes every gene from the reference genome in the raw matrix, but restricts the filtered matrix to genes actually covered by probes in the probe set (here: Chromium Mouse Transcriptome Probe Set v1.1.1). The genes present in the filtered matrix are always a strict subset of the genes in the raw matrix. The current `merge.py` template subsets only barcodes (rows) after loading the CellBender output: ```python adata_cellbender = adata_cellbender[adata.obs_names] ``` When filtered and raw have different gene counts (Flex case), assigning `adata_cellbender.layers["cellbender"]` into `adata.layers[output_layer]` then fails because the gene dimensions disagree. ## Proposed fix Also subset genes (columns) to match the filtered AnnData: ```python adata_cellbender = adata_cellbender[adata.obs_names, adata.var_names] ``` This is the entire change needed in `modules/nf-core/cellbender/merge/templates/merge.py`. ## Why this is safe for standard (non-Flex) data In the standard case, filtered and raw matrices share an identical gene set. Subsetting `adata_cellbender` to `adata.var_names` is then a no-op — the gene order and content are the same, so no data is lost and behaviour is unchanged. For Flex, all probe-targeted genes (the filtered set) are guaranteed to be present in the raw matrix — there is no gene in filtered that is missing from raw, so the index lookup always succeeds. The ~14,700 extra genes dropped from the raw matrix carry negligible signal (empirically: 0.23% of total counts across 354 of 14,704 genes), as they cannot be reliably captured without a targeting probe. ## Reproduction Run `nf-core/scdownstream` with `ambient_correction = 'cellbender'` on a samplesheet pointing to per-sample Cell Ranger multi Flex outputs (`per_sample_outs/<sample>/count/sample_filtered_feature_bc_matrix.h5` as filtered, `sample_raw_feature_bc_matrix.h5` as unfiltered). The pipeline will fail at `CELLBENDER_MERGE` for every sample with the shape mismatch above. The presence of `probe_set.csv` and `sample_raw_probe_bc_matrix.h5` alongside the count matrices is a reliable indicator of Flex data.
0 条评论