CELLBENDER_MERGE fails for 10x Flex data due to gene-set mismatch between filtered and raw matrices
## 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 条评论