Scanpy: normalize_per_cell does not remove zero-expression cells from adata.raw

[gokceneraslan]

When I run this snippet:

import scanpy.api as sc

adata = sc.datasets.paul15()
adata.X[0, :] = 0
adata.raw = adata.copy()

sc.pp.normalize_per_cell(adata)

print(adata.X.shape)
print(adata.raw.X.shape)

sc.pp.pca(adata)
sc.pl.pca_scatter(adata, color='Gata1')

I get:

(2729, 3451)
(2730, 3451)

and an IndexError: index 2729 is out of bounds for axis 0 with size 2729 exception. Actually X and raw.X size mismatch happens every time we filter observations. So what would be the ideal solution?

• Don't remove zero-expression cells in normalize_per_cell
• Remove zero cells also from raw in normalize_per_cell
• Print a warning if some observations are removed in normalize_per_cell
• Whenever some observations are subsetted in adata, also subset raw.X (solution from anndata side)

[wflynny]

My two cents: normalize_per_cell shouldn't remove zero-expression cells. That functionality should explicitly be taken care of by either filter_cells or specific filtering expressions like
{python} adata = adata[adata.obs['n_genes'] > 200, :]
Modifying raw, especially without a warning, shouldn't happen unless it's done directly by the user.

[LuckyMD]

This is very related to a discussion I've been having with Alex about e.g. rank_genes_groups() which is run on the full .raw dataset, but then the visualization happens on a possibly HVG filtered subset which does not have certain indices. So the problem is bigger than just normalize_per_cell().

I think subsetting raw is not really a viable option as you want especially rank_genes_groups() to be run on all genes rather than on a subset which may be chosen for some reason that is unrelated to the genes' importance as marker genes. Rather you will need a solution that revolves around selecting which processing/subselection level of data you want to use to visualize...

[falexwolf]

Thank you for your thoughts!

1. normalize_per_cell needs to remove zero-expression cells as these can't be normalized, the alternative would be to require it as a preprocessing step; but you're right, wflynny, I'll frame it as a fall-back for normalize_per_cell in the next version and output a warning... which will make things backwards compatible...
2. Any filtering operation on the cells/observations should also affect .raw. I'll look into this today.
3. Any filtering operation on the variables should not affect .raw. I didn't know that this gives problems in rank_genes_groups? Of course, you don't find everything in .X that you find in .raw.X and you'll get a key error if you try to; but is there a fundamental problem, @LuckyMD?

[LuckyMD]

I would say there is a more general problem. You should be able to expect that sc.tl.rank_genes_groups(adata) and then sc.pl.rank_genes_groups(adata) always works. However, if adata is subsetted to HVGs and then sc.tl.rank_genes_groups() is run with the default of use_raw=True, then you can get genes as top ranked markers which are not in the adata.X subset. Thus, sc.pl.rank_genes_groups() or sc.pl.rank_genes_groups_violin() will fail as the gene in adata.uns['rank_genes_groups']['name'] is not found in adata.var_names. This occurs as the plotting is done on adata.X and not adata.raw.X.

I've tested this when using the gene_symbols parameter, but it feels like it should also happen when just using regular var_names.

[LuckyMD]

The same would also true for any sc.tl function that has the option of use_raw=True and has a plotting equivalent.

[falexwolf]

@LuckyMD: of course, what you remarked should work and I thought it would. Thanks for remarking it again, I'll look into it.

[falexwolf]

1. for now, only partially fixed through more verbose logging output: https://github.com/theislab/scanpy/commit/8dcacda93c91a5466c2ff23bceb6120ce1d5e0cf
2. was a bug that only affected "inplace subsetting", fixed via https://github.com/theislab/anndata/commit/0becb7b068dda31a60fd0ecb24360d0b5e3d3d7f and in anndata 0.6.6
3. hm, of course you should be able to call tl.rank_genes_groups and pl.rank_genes_groups, but https://github.com/theislab/scanpy/blob/8dcacda93c91a5466c2ff23bceb6120ce1d5e0cf/scanpy/plotting/tools/__init__.py#L271-L282 should do the job for pl.rank_genes_groups and https://github.com/theislab/scanpy/blob/8dcacda93c91a5466c2ff23bceb6120ce1d5e0cf/scanpy/plotting/tools/__init__.py#L355-L358 for pl.rank_genes_groups_violin, right? I don't see how you can observe what you describe. let's look an example together.

[LuckyMD]

Ah, yes... this was a commit I made a month ago, which I had forgotten about that fixed the issue at least in part (caddf9b5934301f9cf2048e6bb947161fd84a210).

I recall that I found it harder to fix for pl.rank_genes_groups_violin and so I left that for the time being as I felt it was a longer discussion. Let's do this offline next week.

[umangv]

The fact that normalize_per_cell filters cells caused me some trouble recently, especially because there is nothing in documentation that says this will happen and there the subset indices are not returned.

My use case is a little different: I am running a K-fold cross validation and normalize_per_cell is a step in the process. Since this function filters out some cells, the array of predicted labels has fewer entries than I would expect for the fold.

Ideally, this behavior could be switched off. For example, you could pass a min_counts argument to normalize_per_cell.

I've made this change if this is an approach you would like to take: https://github.com/umangv/scanpy/tree/umangv-normalize

[falexwolf]

That's a very straight-forward solution viable for the current way in which this is written. A pull request would be very welcome! Thank you!

[umangv]

I made the PR. Thanks!