Scanpy: Plot of ranked genes groups with non-raw data or layers

[SamueleSoraggi]

Hej,

I have been trying to plot the ranked gene groups with a dotplot in this way

sc.plotting.tools.rank_genes_groups_dotplot(all_data_flt_clst, n_genes=5, save='.pdf', groupby='batch', dendrogram=False, layer='imputed') 

or from non-raw data

sc.plotting.tools.rank_genes_groups_dotplot(all_data_flt_clst, n_genes=5, save='.pdf', groupby='batch', dendrogram=False, use_raw=False) 

I get an error regarding the index names of the genes. Other commands for plotting gene rank groups (such as heatmap) show the same problem. The error goes as follows:

KeyError: 'Indices "[\'PRM2\', \'ACAP1\', \'SPEM1\', \'SPATA3\', \'C10orf62\', \'TNP1\', \'MIR193BHG\', \'PRM1\', \'CCDC179\', \'AC007557.1\', \'SPACA1\', \'ERICH2\', \'RP11-360D2.1\', \'TIPARP-AS1\', \'GS1-124K5.4\', \'DCN\', \'C1S\', \'SERPING1\', \'C1R\', \'SERPINF1\', \'GAGE2A\', \'PTMA\', \'HMGB1\', \'VCX2\', \'ERP29\', \'ZCWPW1\', \'SMC1B\', \'DPH7\', \'SCML1\', \'CLSPN\', \'CSAD\', \'C1QBP\', \'DNAJB6\', \'TCF3\', \'RRBP1\', \'HSP90AA1\', \'TMED10\', \'ART3\', \'BUB1\', \'KRBOX1\', \'B2M\', \'IFITM3\', \'GNG11\', \'IFITM2\', \'IFI27\', \'TYROBP\', \'S100A4\', \'FCER1G\', \'CD163\', \'CYBA\', \'CALD1\', \'IGFBP7\', \'TIMP3\', \'PTGDS\', \'TSHZ2\', \'MT-ND3\', \'MT-ND1\', \'MT-ND2\', \'MT-ATP6\', \'MT-ND4\', \'MT-ND1\', \'HMGN5\', \'MT-ND3\', \'ALDH1A1\', \'MT-ATP6\']" contain invalid observation/variables names/indices.'

I checked the var_names and obs_names of my object, but they seem totally fine. Do you have any ideas about the origin of the problem? I can post the backtracking of the error if needed :)

Cheers,
Samuele

[falexwolf]

Sorry, hard to figure out what's going on. Yes, please paste the traceback.

[SamueleSoraggi]

Yeah it is pretty cryptical, but I hoped it happened already to some of you :) Here is the traceback.

---------------------------------------------------------------------------
KeyError                                  Traceback (most recent call last)
<ipython-input-25-c1200217276b> in <module>
----> 1 sc.plotting.tools.rank_genes_groups_stacked_violin(all_data_bygroup, n_genes=5, save='.pdf', use_raw=True, groupby='batch', stripplot=True)

~/miniconda3/envs/scRNA/lib/python3.6/site-packages/scanpy/plotting/tools/__init__.py in rank_genes_groups_stacked_violin(adata, groups, n_genes, groupby, key, show, save, **kwds)
    439 
    440     _rank_genes_groups_plot(adata, plot_type='stacked_violin', groups=groups, n_genes=n_genes,
--> 441                             groupby=groupby, key=key, show=show, save=save, **kwds)
    442 
    443 

~/miniconda3/envs/scRNA/lib/python3.6/site-packages/scanpy/plotting/tools/__init__.py in _rank_genes_groups_plot(adata, plot_type, groups, n_genes, groupby, key, show, save, **kwds)
    305         from ..anndata import stacked_violin
    306         return stacked_violin(adata, gene_names, groupby, var_group_labels=group_names,
--> 307                        var_group_positions=group_positions, show=show, save=save, **kwds)
    308 
    309     elif plot_type == 'tracksplot':

~/miniconda3/envs/scRNA/lib/python3.6/site-packages/scanpy/plotting/anndata.py in stacked_violin(adata, var_names, groupby, log, use_raw, num_categories, figsize, dendrogram, var_group_positions, var_group_labels, var_group_rotation, layer, stripplot, jitter, size, scale, order, swap_axes, show, save, row_palette, **kwds)
    819     if isinstance(var_names, str):
    820         var_names = [var_names]
--> 821     categories, obs_tidy = _prepare_dataframe(adata, var_names, groupby, use_raw, log, num_categories, layer=layer)
    822 
    823     if 'color' in kwds:

~/miniconda3/envs/scRNA/lib/python3.6/site-packages/scanpy/plotting/anndata.py in _prepare_dataframe(adata, var_names, groupby, use_raw, log, num_categories, layer)
   1983         matrix = adata[:, var_names].layers[layer]
   1984     elif use_raw:
-> 1985         matrix = adata.raw[:, var_names].X
   1986     else:
   1987         matrix = adata[:, var_names].X

~/miniconda3/envs/scRNA/lib/python3.6/site-packages/anndata/base.py in __getitem__(self, index)
    510 
    511     def __getitem__(self, index):
--> 512         oidx, vidx = self._normalize_indices(index)
    513         if self._adata is not None or not self._adata.isbacked: X = self._X[oidx, vidx]
    514         else: X = self._adata.file['raw.X'][oidx, vidx]

~/miniconda3/envs/scRNA/lib/python3.6/site-packages/anndata/base.py in _normalize_indices(self, packed_index)
    538         obs, var = super(Raw, self)._unpack_index(packed_index)
    539         obs = _normalize_index(obs, self._adata.obs_names)
--> 540         var = _normalize_index(var, self.var_names)
    541         return obs, var
    542 

~/miniconda3/envs/scRNA/lib/python3.6/site-packages/anndata/base.py in _normalize_index(index, names)
    270             raise KeyError(
    271                 'Indices "{}" contain invalid observation/variables names/indices.'
--> 272                 .format(index))
    273         return positions.values
    274     else:

Cheers,
Samuele

[fidelram]

Looks like if for some reason the gene names that you want to plot are not
in your anndata object.

can you run sc.pl.umap(all_data_flt_clst, color=['PRM2']) ?

On Wed, Jan 23, 2019 at 9:07 AM Samuele Soraggi [e-mail]
wrote:

Yeah it is pretty cryptical, but I hoped it happened already to some of
you :) Here is the traceback.

---

KeyError Traceback (most recent call last)
in
----> 1 sc.plotting.tools.rank_genes_groups_stacked_violin(all_data_bygroup, n_genes=5, save='.pdf', use_raw=True, groupby='batch', stripplot=True)

~/miniconda3/envs/scRNA/lib/python3.6/site-packages/scanpy/plotting/tools/__init__.py in rank_genes_groups_stacked_violin(adata, groups, n_genes, groupby, key, show, save, *kwds)
439
440 _rank_genes_groups_plot(adata, plot_type='stacked_violin', groups=groups, n_genes=n_genes,
--> 441 groupby=groupby, key=key, show=show, save=save, *kwds)
442
443

~/miniconda3/envs/scRNA/lib/python3.6/site-packages/scanpy/plotting/tools/__init__.py in _rank_genes_groups_plot(adata, plot_type, groups, n_genes, groupby, key, show, save, *kwds)
305 from ..anndata import stacked_violin
306 return stacked_violin(adata, gene_names, groupby, var_group_labels=group_names,
--> 307 var_group_positions=group_positions, show=show, save=save, *kwds)
308
309 elif plot_type == 'tracksplot':

~/miniconda3/envs/scRNA/lib/python3.6/site-packages/scanpy/plotting/anndata.py in stacked_violin(adata, var_names, groupby, log, use_raw, num_categories, figsize, dendrogram, var_group_positions, var_group_labels, var_group_rotation, layer, stripplot, jitter, size, scale, order, swap_axes, show, save, row_palette, **kwds)
819 if isinstance(var_names, str):
820 var_names = [var_names]
--> 821 categories, obs_tidy = _prepare_dataframe(adata, var_names, groupby, use_raw, log, num_categories, layer=layer)
822
823 if 'color' in kwds:

~/miniconda3/envs/scRNA/lib/python3.6/site-packages/scanpy/plotting/anndata.py in _prepare_dataframe(adata, var_names, groupby, use_raw, log, num_categories, layer)
1983 matrix = adata[:, var_names].layers[layer]
1984 elif use_raw:
-> 1985 matrix = adata.raw[:, var_names].X
1986 else:
1987 matrix = adata[:, var_names].X

~/miniconda3/envs/scRNA/lib/python3.6/site-packages/anndata/base.py in __getitem__(self, index)
510
511 def __getitem__(self, index):
--> 512 oidx, vidx = self._normalize_indices(index)
513 if self._adata is not None or not self._adata.isbacked: X = self._X[oidx, vidx]
514 else: X = self._adata.file['raw.X'][oidx, vidx]

~/miniconda3/envs/scRNA/lib/python3.6/site-packages/anndata/base.py in _normalize_indices(self, packed_index)
538 obs, var = super(Raw, self)._unpack_index(packed_index)
539 obs = _normalize_index(obs, self._adata.obs_names)
--> 540 var = _normalize_index(var, self.var_names)
541 return obs, var
542

~/miniconda3/envs/scRNA/lib/python3.6/site-packages/anndata/base.py in _normalize_index(index, names)
270 raise KeyError(
271 'Indices "{}" contain invalid observation/variables names/indices.'
--> 272 .format(index))
273 return positions.values
274 else:

Cheers,
Samuele

—
You are receiving this because you are subscribed to this thread.
Reply to this email directly, view it on GitHub
https://github.com/theislab/scanpy/issues/438#issuecomment-456707222,
or mute the thread
https://github.com/notifications/unsubscribe-auth/AEu_1XEPWhKAyeK0sWLrAzmqJvm45H-vks5vGBhDgaJpZM4aMT_6
.

[SamueleSoraggi]

Yeah, the UMAP plots work alright, and I can recognize many of the genes I get in the message are indeed genes that are expressed in the data, for I can visualize them on my UMAP plots. It is kind of weird, I am wondering if somewhere I messed up the format of those names, but then why would they work on UMAP?

[LuckyMD]

Have you subsetted your AnnData object to highly variable genes, while keeping the full dataset in .raw? In that case it could be that genes that are found as markers via rank_genes_groups, are not in adata.var_names, but only in adata.raw.var_names and therefore cannot be found by the plotting function. I've previously encountered issues with this, but I thought it had been solved now.

[SamueleSoraggi]

Yeah, the raw data has been indeed assigned before I subsetted through highly variable genes. It might be this sort of mismatching that is problematic.
However I encounter the same problem when trying to plot a layer. The layers should contain the same set of genes as the data.

[LuckyMD]

I'm not entirely sure what you mean by plotting a layer here. You mean that you run rank_genes_groups on a layer? Or is plotting a layer completely unrelated to rank_genes_groups?

[fidelram]

@SamueleSoraggi you are right, the layers contain the same genes as the adata.X matrix. I assume that in your case, you did a highly variable gene selection which affects both adata.X and adata.layers but not adata.raw.

The solution is to mark highly variable genes without removing the other less variable genes. This functionality was added some few months ago and may not be properly reflected on the documentation.

[MichaelPeibo]

Hi @SamueleSoraggi @fidelram
I came to the same problem---unable to find genes after filtering by HVGs.
Did you figure it out and go through?

Or I can only use unfiltered matrix to plot? However, if I do not filter matrix with HVG, there is a slight difference between clustering and take a lot of time running sc.tl.rank_genes_groups.

Any suggestion?

[fidelram]

@MichaelPeibo how did you run the HVG? There is as parameter called subset_data or something like that, which will allow you to keep all genes and only flag the HVG.

[MichaelPeibo]

Hi @fidelram
Thanks for pointing it out!

I used default parameters of subset=False in scanpy.pp.highly_variable_genes.

Actually, I think I found where the trick is. I used .concatenate to combine multiple datasets, and default parameter of join is inner, which will only filter common colomns/genes? Following this, I got some interesting genes missing when I did violin plot;

In the second time, I tried to set join=outer , which will keep more genes I think? Following this, I got some interesting genes found when I did violin plot;

Am I doing the meaningful or rationale way?

[falexwolf]

Yes, what you describe makes sense! :)

[SamueleSoraggi]

I stilll encounter the problem, even when I set adata.raw=adata.
However, it works with use_raw=false or by choosing a layer when I run the scanpy.plotting.rank_genes_groups function and does not when I use other plotting functions like sc.plotting.rank_genes_groups_dotplot.
It is weird, I guess they are based on the same way of looking at the rank of genes.