scverse / scverse/scanpy

`calculate_qc_metrics(use_raw=True)` labels `.raw`'s matrix with `adata.var`

Open
#4,347 0 comments 0 reactions 0 assignees View on GitHub

Nobody has claimed this yet.

Dominant language
Python
Stars
2.6k
Forks
779
Avg merge
1d 4h
Merged PRs (30d)
27

Description

Please make sure these conditions are met
  • I have checked that this issue has not already been reported.
  • I have confirmed this bug exists on the latest version of scanpy.
  • (optional) I have confirmed this bug exists on the main branch of scanpy.
What happened?

sc.pp.calculate_qc_metrics(adata, use_raw=True) takes the matrix from .raw and then labels and masks it with adata.var. Those two agree only while .raw and .X share a var axis, and they stop sharing one as soon as .X is subset to the highly variable genes.

The ordinary highly variable gene workflow leaves that shape. Set .raw on the log-normalised matrix, subset .X to the highly variable genes, then ask for QC over everything that was measured. This runs on sc.datasets.pbmc68k_reduced().raw.to_adata(), whose 765 genes are already log-normalised, so nothing is transformed again:

after adata.raw = adata       (700, 765) with .raw (700, 765)
after subsetting to 500 HVGs  (700, 500) with .raw (700, 765)
.raw still the untouched log-normalised matrix: True

>>> sc.pp.calculate_qc_metrics(adata, use_raw=True)
    ValueError: Length of values (765) does not match length of index (500)

The same shape turns up in scanpy's own docs/tutorials/trajectories/paga-paul15.ipynb. In [6] runs sc.pp.recipe_zheng17(adata), which subsets .X in place to n_top_genes=1000, and In [32] sets adata.raw = adata_raw holding all 3451 genes of sc.datasets.paul15(). That notebook never asks for QC metrics. Running those two cells and then asking:

sc.datasets.paul15()         (2730, 3451)
after sc.pp.recipe_zheng17   (2730, 1000)
adata                        (2730, 1000)
adata.raw                    (2730, 3451)

>>> sc.pp.calculate_qc_metrics(adata, use_raw=True)
    ValueError: Length of values (3451) does not match length of index (1000)

That run also prints a RuntimeWarning: invalid value encountered in log1p, because the notebook puts log1p'd and scaled values in .raw and 1488 of the 2730 per-cell totals fall below -1. Nothing to do with this.

inplace=True raises too, but only after writing four .raw-derived columns into adata.obs. On the sample object below it leaves total_counts [7., 4., 8.], which are the .raw row sums, on an object whose .X row sums are [3., 3., 5.], and no qc columns in .var at all.

There is a one line workaround for the raise: sc.pp.calculate_qc_metrics(adata.raw.to_adata(), percent_top=None) returns the right numbers, a (2730, 4) obs frame and a (3451, 6) var frame on the paul15 object. It has limits worth saying out loud. It copies the whole raw matrix. Its annotations land on the new object, so it cannot serve inplace=True. And qc_vars only survives the trip when the column was on .var before .raw was set, because that is what .raw.var keeps: annotate mito afterwards, which is the ordinary order, and the same call raises KeyError: 'mito'. The sample below runs both orders.

The half I think is worth more than the traceback is the quiet one. When .raw and .X hold the same genes in a different order, nothing raises, the numbers land on the wrong rows, and inplace=True leaves them in adata.var. In the sample below g2 is reported with a total_counts of 6, which is g0's, and total_counts_mito comes back as [0., 1., 2.], which is g2's column, when mito is g0 alone whose .raw column is [1., 0., 5.]. I know of no scanpy workflow that produces that order, so I am reporting it as the other end of one bug and not as a second one.

Where it comes from, at ec374022:

  • calculate_qc_metrics reads the matrix with use_raw at _qc.py:279, then calls describe_obs and describe_var passing x but not use_raw (_qc.py:287-304). Both helpers run with their default use_raw=False while x carries the .raw var axis.
  • describe_var builds its result as pd.DataFrame(index=adata.var_names) (_qc.py:177).
  • describe_obs masks with adata.var[qc_var].to_numpy() (_qc.py:114).

Three siblings switch the axis on use_raw and _qc.py takes adata.var either way: get.obs_df sets var = adata.raw.var when use_raw and adata.var otherwise (get.py:295-301), tl.score_genes uses adata.raw.var_names if use_raw else adata.var_names (_score_genes.py:236), and tl.dendrogram does the same (_dendrogram.py:145).

tests/test_qc_metrics.py::test_layer_raw does cover use_raw=True, but line 286 is adata.raw = adata.copy(), so the two axes are identical there and no length can disagree.

Two ways to close it, and I don't know which you would want:

  1. Follow get.py:295-301 and take both the index and the qc_vars mask from adata.raw. I tried it. Every use_raw=True call above returns and the reordering case comes out right. It leaves inplace=True writing only the rows whose names are also in adata.var, two of the four on the small object, dropped silently because pandas aligns on the index. And it inherits the workaround's limitation from the same root cause, .raw.var never having carried the column: a mito column added to adata.var after .raw was set becomes KeyError: 'mito' on this path too. The late.raw.to_adata() call is untouched by this option and still raises.
  2. Raise when use_raw=True and adata.raw.var_names.equals(adata.var_names) is False. That invents no contract and turns both halves above into one sentence. anndata.acc.A exposes no raw, so if .raw is not coming along to the new accessor API, this is probably the one that fits. Nothing shipped would start raising: no caller in src/, no tutorial and no docstring example passes use_raw=True to this function.

The output below is from uv run on the sample file, which built scanpy from main at ec374022 and ran it on python 3.14. The same file with the dependency changed to "scanpy" installs 1.12.4 and prints the same values and the same ValueError. The pbmc68k and paul15 runs above were both measured in an editable checkout of ec374022 on python 3.12, which is what the Versions block reports.

Minimal code sample
# /// script
# requires-python = ">=3.12"
# dependencies = [
#   "scanpy@git+https://github.com/scverse/scanpy.git@main",
# ]
# ///

import anndata as ad
import numpy as np
import pandas as pd

import scanpy as sc

counts = np.array([[1, 2, 0, 4], [0, 3, 1, 0], [5, 0, 2, 1]], dtype=np.float32)
full = ad.AnnData(
    counts,
    obs=pd.DataFrame(index=["c0", "c1", "c2"]),
    var=pd.DataFrame(index=["g0", "g1", "g2", "g3"]),
)
full.var["mito"] = [True, False, False, False]  # g0 alone, .raw column [1., 0., 5.]
full.raw = full

# percent_top=None throughout only because the default asks for the top 500
# of 4 genes, which is an unrelated IndexError.
subset = full[:, ["g0", "g1"]].copy()  # what an HVG subset leaves behind
reordered = full[:, ["g2", "g3", "g0", "g1"]].copy()  # same genes, other order


def show(label, fn):
    print(f">>> {label}")
    try:
        return fn()
    except Exception as e:  # noqa: BLE001
        print(f"    {type(e).__name__}: {e}")
        return None


print("subset  ", subset.shape, "with .raw", subset.raw.shape)
show(
    "calculate_qc_metrics(subset, use_raw=True, percent_top=None)",
    lambda: sc.pp.calculate_qc_metrics(subset, use_raw=True, percent_top=None),
)
show(
    'calculate_qc_metrics(subset, use_raw=True, percent_top=None, qc_vars=["mito"])',
    lambda: sc.pp.calculate_qc_metrics(
        subset, use_raw=True, percent_top=None, qc_vars=["mito"]
    ),
)
show(
    "calculate_qc_metrics(subset, use_raw=True, percent_top=None, inplace=True)",
    lambda: sc.pp.calculate_qc_metrics(
        subset, use_raw=True, percent_top=None, inplace=True
    ),
)
print("    adata.obs now has ", list(subset.obs.columns))
print("    obs total_counts  ", subset.obs["total_counts"].to_numpy())
print("    .X row sums       ", np.asarray(subset.X).sum(axis=1))
print("    .raw row sums     ", np.asarray(subset.raw.X).sum(axis=1))
print("    adata.var now has ", list(subset.var.columns))

print()
obs, var = show(
    'calculate_qc_metrics(subset.raw.to_adata(), percent_top=None, qc_vars=["mito"])',
    lambda: sc.pp.calculate_qc_metrics(
        subset.raw.to_adata(), percent_top=None, qc_vars=["mito"]
    ),
)
print("    var index         ", list(var.index))
print("    total_counts      ", var["total_counts"].to_numpy())
print("    total_counts_mito ", obs["total_counts_mito"].to_numpy())

# The workaround's one limit: qc_vars only survives if the column was on .var
# before .raw was set, because that is what .raw.var keeps.
late = ad.AnnData(
    counts.copy(),
    obs=pd.DataFrame(index=["c0", "c1", "c2"]),
    var=pd.DataFrame(index=["g0", "g1", "g2", "g3"]),
)
late.raw = late
late.var["mito"] = [True, False, False, False]  # annotated after .raw was set
late = late[:, ["g0", "g1"]].copy()
print()
print("mito in .var", "mito" in late.var, "| mito in .raw.var", "mito" in late.raw.var)
show(
    'calculate_qc_metrics(late.raw.to_adata(), percent_top=None, qc_vars=["mito"])',
    lambda: sc.pp.calculate_qc_metrics(
        late.raw.to_adata(), percent_top=None, qc_vars=["mito"]
    ),
)

print()
print("reordered", reordered.shape, "with .raw", reordered.raw.shape)
_, var = show(
    "calculate_qc_metrics(reordered, use_raw=True, percent_top=None)",
    lambda: sc.pp.calculate_qc_metrics(reordered, use_raw=True, percent_top=None),
)
print("    var index         ", list(var.index))
print("    total_counts      ", var["total_counts"].to_numpy())
print("    .raw var order    ", list(reordered.raw.var_names))
print("    .raw column sums  ", np.asarray(reordered.raw.X).sum(axis=0))
obs, _ = show(
    'calculate_qc_metrics(reordered, use_raw=True, percent_top=None, qc_vars=["mito"])',
    lambda: sc.pp.calculate_qc_metrics(
        reordered, use_raw=True, percent_top=None, qc_vars=["mito"]
    ),
)
print("    total_counts_mito ", obs["total_counts_mito"].to_numpy())
show(
    "calculate_qc_metrics(reordered, use_raw=True, percent_top=None, inplace=True)",
    lambda: sc.pp.calculate_qc_metrics(
        reordered, use_raw=True, percent_top=None, inplace=True
    ),
)
print("    adata.var index   ", list(reordered.var.index))
print("    adata.var totals  ", reordered.var["total_counts"].to_numpy())

print()
sc.pp.calculate_qc_metrics(subset, use_raw=True, percent_top=None)
Error output
subset   (3, 2) with .raw (3, 4)
>>> calculate_qc_metrics(subset, use_raw=True, percent_top=None)
    ValueError: Length of values (4) does not match length of index (2)
>>> calculate_qc_metrics(subset, use_raw=True, percent_top=None, qc_vars=["mito"])
    IndexError: boolean index did not match indexed array along axis 1; size of axis is 4 but size of corresponding boolean axis is 2
>>> calculate_qc_metrics(subset, use_raw=True, percent_top=None, inplace=True)
    ValueError: Length of values (4) does not match length of index (2)
    adata.obs now has  ['n_genes_by_counts', 'log1p_n_genes_by_counts', 'total_counts', 'log1p_total_counts']
    obs total_counts   [7. 4. 8.]
    .X row sums        [3. 3. 5.]
    .raw row sums      [7. 4. 8.]
    adata.var now has  ['mito']

>>> calculate_qc_metrics(subset.raw.to_adata(), percent_top=None, qc_vars=["mito"])
    var index          ['g0', 'g1', 'g2', 'g3']
    total_counts       [6. 5. 3. 5.]
    total_counts_mito  [1. 0. 5.]

mito in .var True | mito in .raw.var False
>>> calculate_qc_metrics(late.raw.to_adata(), percent_top=None, qc_vars=["mito"])
    KeyError: 'mito'

reordered (3, 4) with .raw (3, 4)
>>> calculate_qc_metrics(reordered, use_raw=True, percent_top=None)
    var index          ['g2', 'g3', 'g0', 'g1']
    total_counts       [6. 5. 3. 5.]
    .raw var order     ['g0', 'g1', 'g2', 'g3']
    .raw column sums   [6. 5. 3. 5.]
>>> calculate_qc_metrics(reordered, use_raw=True, percent_top=None, qc_vars=["mito"])
    total_counts_mito  [0. 1. 2.]
>>> calculate_qc_metrics(reordered, use_raw=True, percent_top=None, inplace=True)
    adata.var index    ['g2', 'g3', 'g0', 'g1']
    adata.var totals   [6. 5. 3. 5.]

Traceback (most recent call last):
  File "/home/levi/scanpy-issue-repro/issue.py", line 118, in <module>
    sc.pp.calculate_qc_metrics(subset, use_raw=True, percent_top=None)
    ~~~~~~~~~~~~~~~~~~~~~~~~~~^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/levi/.cache/uv/environments-v2/issue-2c5ce2c7d0fd00b0/lib/python3.14/site-packages/scanpy/preprocessing/_qc.py", line 297, in calculate_qc_metrics
    var_metrics = describe_var(
        adata,
    ...<4 lines>...
        log1p=log1p,
    )
  File "/home/levi/.cache/uv/environments-v2/issue-2c5ce2c7d0fd00b0/lib/python3.14/site-packages/scanpy/preprocessing/_qc.py", line 178, in describe_var
    var_metrics[f"n_cells_by_{expr_type}"], var_metrics[f"mean_{expr_type}"] = (
    ~~~~~~~~~~~^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/levi/.cache/uv/environments-v2/issue-2c5ce2c7d0fd00b0/lib/python3.14/site-packages/pandas/core/frame.py", line 4672, in __setitem__
    self._set_item(key, value)
    ~~~~~~~~~~~~~~^^^^^^^^^^^^
  File "/home/levi/.cache/uv/environments-v2/issue-2c5ce2c7d0fd00b0/lib/python3.14/site-packages/pandas/core/frame.py", line 4874, in _set_item
    value, refs = self._sanitize_column(value)
                  ~~~~~~~~~~~~~~~~~~~~~^^^^^^^
  File "/home/levi/.cache/uv/environments-v2/issue-2c5ce2c7d0fd00b0/lib/python3.14/site-packages/pandas/core/frame.py", line 5756, in _sanitize_column
    com.require_length_match(value, self.index)
    ~~~~~~~~~~~~~~~~~~~~~~~~^^^^^^^^^^^^^^^^^^^
  File "/home/levi/.cache/uv/environments-v2/issue-2c5ce2c7d0fd00b0/lib/python3.14/site-packages/pandas/core/common.py", line 601, in require_length_match
    raise ValueError(
    ...<4 lines>...
    )
ValueError: Length of values (4) does not match length of index (2)
Versions
scanpy	1.14.0.dev4+gec3740223
----	----
charset-normalizer	3.5.1
legacy-api-wrap	1.5
cycler	0.12.1
psutil	7.2.2
numba	0.67.0
session-info2	0.4.2
scipy	1.18.1
narwhals	2.25.0
fonttools	4.64.0
pyparsing	3.3.2
PyYAML	6.0.3
pydantic	2.13.5
matplotlib	3.11.1
typing_extensions	4.16.0
annotated-types	0.8.0
pillow	12.3.0
pydantic-settings	2.15.0
pandas	3.0.5
donfig	0.8.1.post1
python-dateutil	2.9.0.post0
threadpoolctl	3.6.0
h5py	3.16.0
zarr	3.3.0
packaging	26.3
cloudpickle	3.1.2
scikit-learn	1.9.0
anndata	0.13.3.post0
pydantic_core	2.46.5
python-dotenv	1.2.3
fast-array-utils	1.5
natsort	8.4.0
llvmlite	0.49.0
six	1.17.0
joblib	1.6.0
typing-inspection	0.4.4
google-crc32c	1.8.0
coverage	7.16.0
numcodecs	0.16.5
kiwisolver	1.5.1
scverse-misc	0.1.4
numpy	2.5.2
----	----
Python	3.12.3 (main, Jun 19 2026, 12:46:00) [GCC 13.3.0]
OS	Linux-7.0.0-28-generic-x86_64-with-glibc2.39
CPU	12/12 logical CPU cores, x86_64
GPU	ID: 0, NVIDIA GeForce RTX 5090, Driver: 610.43.02, Memory: 32607 MiB
Updated	2026-09-06 14:13

Contributor guide

Open the contributing guide

First steps

  1. Read the whole issue, then the project's contributing guide.
  2. Comment on the issue to say you are picking it up — it saves two people doing the same work.
  3. Fork the repository and make your change on a branch.
  4. Open a pull request that references the issue number.

Research direction

Read calculate_qc_metrics, describe_obs, and describe_var in _qc.py, then compare their use_raw handling with get.py:295-301. Extend tests/test_qc_metrics.py beyond test_layer_raw to cover subset and reordered .raw axes, including qc_vars; done means the behavior is decided and these cases no longer mislabel or fail unexpectedly.

Written by the indexing model from the issue text.

Assessment

Tech stack
numpy, pandas, python
Domain
bioinformatics, data
Issue type
Bug
Difficulty
4/5
Estimated time
3-5 days
Activity status
Active
Clarity
Mostly clear
Newbie friendliness
58/100

Get new issues in your inbox

A short digest of beginner-friendly GitHub issues.