Error: "invalid class dgRMatrix object"

[LuckyMD]

Hi Phil,

I'm just updating my case study notebook and updated anndata2ri for that. However, when converting my anndata object via:

%%R -i adata_ent
print('test')

I get the following error:

---------------------------------------------------------------------------
RRuntimeError                             Traceback (most recent call last)
<ipython-input-81-f8bea5081b84> in <module>
----> 1 get_ipython().run_cell_magic('R', '-i adata_ent', "print('test')\n")

~/anaconda3/envs/sc-tutorial/lib/python3.6/site-packages/IPython/core/interactiveshell.py in run_cell_magic(self, magic_name, line, cell)
   2321             magic_arg_s = self.var_expand(line, stack_depth)
   2322             with self.builtin_trap:
-> 2323                 result = fn(magic_arg_s, cell)
   2324             return result
   2325 

<decorator-gen-797> in R(self, line, cell, local_ns)

~/anaconda3/envs/sc-tutorial/lib/python3.6/site-packages/IPython/core/magic.py in <lambda>(f, *a, **k)
    185     # but it's overkill for just that one bit of state.
    186     def magic_deco(arg):
--> 187         call = lambda f, *a, **k: f(*a, **k)
    188 
    189         if callable(arg):

~/anaconda3/envs/sc-tutorial/lib/python3.6/site-packages/rpy2/ipython/rmagic.py in R(self, line, cell, local_ns)
    704                         raise NameError("name '%s' is not defined" % input)
    705                 with localconverter(converter) as cv:
--> 706                     ro.r.assign(input, val)
    707 
    708         tmpd = self.setup_graphics(args)

~/anaconda3/envs/sc-tutorial/lib/python3.6/site-packages/rpy2/robjects/functions.py in __call__(self, *args, **kwargs)
    190                 kwargs[r_k] = v
    191         return (super(SignatureTranslatedFunction, self)
--> 192                 .__call__(*args, **kwargs))
    193 
    194 

~/anaconda3/envs/sc-tutorial/lib/python3.6/site-packages/rpy2/robjects/functions.py in __call__(self, *args, **kwargs)
    111 
    112     def __call__(self, *args, **kwargs):
--> 113         new_args = [conversion.py2rpy(a) for a in args]
    114         new_kwargs = {}
    115         for k, v in kwargs.items():

~/anaconda3/envs/sc-tutorial/lib/python3.6/site-packages/rpy2/robjects/functions.py in <listcomp>(.0)
    111 
    112     def __call__(self, *args, **kwargs):
--> 113         new_args = [conversion.py2rpy(a) for a in args]
    114         new_kwargs = {}
    115         for k, v in kwargs.items():

~/anaconda3/envs/sc-tutorial/lib/python3.6/functools.py in wrapper(*args, **kw)
    805                             '1 positional argument')
    806 
--> 807         return dispatch(args[0].__class__)(*args, **kw)
    808 
    809     funcname = getattr(func, '__name__', 'singledispatch function')

~/anaconda3/envs/sc-tutorial/lib/python3.6/site-packages/anndata2ri/py2r.py in py2rpy_anndata(obj)
     69         # Convert everything we know
     70         with localconverter(full_converter() + dict_converter):
---> 71             metadata = ListVector(obj.uns.items())
     72 
     73         rd_args = {conv_name.scanpy2sce(k): mat_converter.py2rpy(obj.obsm[k]) for k in obj.obsm.keys()}

~/anaconda3/envs/sc-tutorial/lib/python3.6/site-packages/rpy2/robjects/vectors.py in __init__(self, tlist)
    628                 raise ValueError('tlist should have a /method/ __iter__ '
    629                                  '(not an attribute)')
--> 630             kv = [(str(k), conversion.py2rpy(v)) for k, v in tlist]
    631             kv = tuple(kv)
    632             df = baseenv_ri.find("list").rcall(kv, globalenv_ri)

~/anaconda3/envs/sc-tutorial/lib/python3.6/site-packages/rpy2/robjects/vectors.py in <listcomp>(.0)
    628                 raise ValueError('tlist should have a /method/ __iter__ '
    629                                  '(not an attribute)')
--> 630             kv = [(str(k), conversion.py2rpy(v)) for k, v in tlist]
    631             kv = tuple(kv)
    632             df = baseenv_ri.find("list").rcall(kv, globalenv_ri)

~/anaconda3/envs/sc-tutorial/lib/python3.6/functools.py in wrapper(*args, **kw)
    805                             '1 positional argument')
    806 
--> 807         return dispatch(args[0].__class__)(*args, **kw)
    808 
    809     funcname = getattr(func, '__name__', 'singledispatch function')

~/anaconda3/envs/sc-tutorial/lib/python3.6/site-packages/anndata2ri/py2r.py in py2rpy_dict(obj)
     33     for k, v in obj.items():
     34         try:
---> 35             converted[str(k)] = conversion.py2rpy(v)
     36         except NotImplementedError as e:
     37             warn(str(e), NotConvertedWarning)

~/anaconda3/envs/sc-tutorial/lib/python3.6/functools.py in wrapper(*args, **kw)
    805                             '1 positional argument')
    806 
--> 807         return dispatch(args[0].__class__)(*args, **kw)
    808 
    809     funcname = getattr(func, '__name__', 'singledispatch function')

~/anaconda3/envs/sc-tutorial/lib/python3.6/site-packages/anndata2ri/scipy2ri/py2r.py in wrapper(obj)
     40 
     41         with localconverter(default_converter + numpy2ri.converter):
---> 42             return f(obj)
     43 
     44     return wrapper

~/anaconda3/envs/sc-tutorial/lib/python3.6/site-packages/anndata2ri/scipy2ri/py2r.py in csr_to_rmat(csr)
     67         p=as_integer(csr.indptr),
     68         x=conv_data(csr.data),
---> 69         Dim=as_integer(list(csr.shape)),
     70     )
     71 

~/anaconda3/envs/sc-tutorial/lib/python3.6/site-packages/rpy2/robjects/functions.py in __call__(self, *args, **kwargs)
    190                 kwargs[r_k] = v
    191         return (super(SignatureTranslatedFunction, self)
--> 192                 .__call__(*args, **kwargs))
    193 
    194 

~/anaconda3/envs/sc-tutorial/lib/python3.6/site-packages/rpy2/robjects/functions.py in __call__(self, *args, **kwargs)
    119             else:
    120                 new_kwargs[k] = conversion.py2rpy(v)
--> 121         res = super(Function, self).__call__(*new_args, **new_kwargs)
    122         res = conversion.rpy2py(res)
    123         return res

~/anaconda3/envs/sc-tutorial/lib/python3.6/site-packages/rpy2/rinterface_lib/conversion.py in _(*args, **kwargs)
     26 def _cdata_res_to_rinterface(function):
     27     def _(*args, **kwargs):
---> 28         cdata = function(*args, **kwargs)
     29         # TODO: test cdata is of the expected CType
     30         return _cdata_to_rinterface(cdata)

~/anaconda3/envs/sc-tutorial/lib/python3.6/site-packages/rpy2/rinterface.py in __call__(self, *args, **kwargs)
    771                     error_occured))
    772             if error_occured[0]:
--> 773                 raise embedded.RRuntimeError(_rinterface._geterrmessage())
    774         return res
    775 

RRuntimeError: Error in validObject(.Object) : 
  invalid class “dgRMatrix” object: slot j is not increasing inside a column
Calls: <Anonymous> ... initialize -> callNextMethod -> .nextMethod -> validObject

I'm pretty sure this has to do with versions in my conda environment as anndata2ri conversion worked fine when I worked locally on my laptop. Any idea if there are any R version dependencies that might not be met?

[LuckyMD]

Is this maybe falsely assuming that I am passing a sparse matrix object? That would explain why it works for a different case on my local machine.

[flying-sheep]

Anndata2ri only calls csr_to_rmat through singledispatch:

anndata2ri/anndata2ri/scipy2ri/py2r.py

Lines 60 to 62 in 3fe0c7d

	@converter.py2rpy.register(sparse.csr_matrix)
	@py2r_context
	def csr_to_rmat(csr: sparse.csr_matrix):

So I’m pretty sure it’s impossible that it “falsely assumes” that you pass a sparse matrix object. I’m, pretty sure you actually pass a scipy.sparse.csr_matrix. The traceback contains these lines:

anndata2ri/anndata2ri/py2r.py

Lines 70 to 71 in 3fe0c7d

	with localconverter(full_converter() + dict_converter):
	metadata = ListVector(obj.uns.items())

So I assume there’s a sparse matrix in adata_ent.uns. Can you please provide a minimal reproducible example?

[LuckyMD]

It's in the case study of the best practices notebook that worked before the recent update (not exactly a minimal example). I will see if there's a sparse matrix in there, and otherwise check when this breaks.

[LuckyMD]

You are correct. There are the following sparse matrices in adata_ent.uns:

  'distances': <9498x9498 sparse matrix of type '<class 'numpy.float64'>'
  	with 132569 stored elements in Compressed Sparse Row format>,
  'connectivities': <9498x9498 sparse matrix of type '<class 'numpy.float64'>'
  	with 215190 stored elements in Compressed Sparse Row format>},

Wasn't it the case before that you ignored sparse matrices in adata.uns? Or was that a general ignoring of sparse matrices that is no longer there as you now do convert adata.X if it's sparse?

Would it be possible to only ignore sparse matrices in adata.uns and put a warning out for that? Or where do SingleCellExperiment objects store the knn adjacency matrices?

[flying-sheep]

I think it’s a case of “couldn’t handle sparse matrices before, can do that now”

I don’t think SCEs have a canonical location for those, so I just put it into the generic “metadata” list.

I don’t think ignoring this is the right way, there seems to be a bug either in scipy or here. My code looks pretty OK I think:

anndata2ri/anndata2ri/scipy2ri/py2r.py

Lines 62 to 70 in 3fe0c7d

	def csr_to_rmat(csr: sparse.csr_matrix):
	t, conv_data, _ = get_type_conv(csr.dtype)
	return methods.new(
	f"{t}gRMatrix",
	j=as_integer(csr.indices),
	p=as_integer(csr.indptr),
	x=conv_data(csr.data),
	Dim=as_integer(list(csr.shape)),
	)

So why does R say “invalid class “dgRMatrix” object: slot j is not increasing inside a column”. Did whatever code created the csr_matrix make a mistake and csr.indices is actually faulty? Or does the data model differ and scipy accepts that structure while the Matrix package doesn’t? Or did my as_integer(csr.indices) return the wrong thing?

[LuckyMD]

The sparse matrices in adata.uns don't seem to be the issue. Just tested this with sc.datasets.blobs() and it works fine after sc.pp.neighbors().

[flying-sheep]

(one of) those specific sparse matrices’ memory layout seems to be the problem, not any sparse matrix stored in .uns

[LuckyMD]

I have an idea what it could be. In this case I subset an anndata object. I assume this doesn't subset the sparse matrices in adata.uns. Maybe that's the issue here. They probably don't have the same dimensions as the subsetted anndata object anymore. Do you assume that is the case for a neighbourhood graph?

[LuckyMD]

So it has to do with subsetting. Although the adata.uns['neighbors'] dictionary is also subsetted. So it's not a shape mismatch. I have a minimal reproducible example now:

import scanpy as sc
from rpy2.robjects import pandas2ri
import anndata2ri
import numpy as np

# Automatically convert rpy2 outputs to pandas dataframes
pandas2ri.activate()
anndata2ri.activate()
%load_ext rpy2.ipython

%%R
library(SingleCellExperiment)

adata = sc.datasets.blobs()
sc.pp.pca(adata)
sc.pp.neighbors(adata)
sc.tl.louvain(adata)
adata_sub = adata[np.isin(adata.obs['louvain'], ['0', '1']),:].copy()

This works:

%%R -i adata
print(adata)

This doesn't work:

%%R -i adata_sub
print(adata_sub)

[LuckyMD]

I added the .copy() in there to ensure it's not a COW/View, but it doesn't work in either case.

[flying-sheep]

Maybe the subsetting results in csr.indices being some slice structure that can’t be properly converted into an R integer vector.

Anyway, I can’t reproduce this:

%%R -i adata_sub
adata_sub

class: SingleCellExperiment 
dim: 11 256 
metadata(3): pca neighbors louvain
assays(1): X
rownames(11): 0 1 ... 9 10
rowData names(0):
colnames(256): 2 3 ... 635 638
colData names(2): blobs louvain
reducedDimNames(1): PCA
spikeNames(0):

scanpy==1.4.3 anndata==0.6.20 umap==0.3.0 numpy==1.16.4 scipy==1.3.0 pandas==0.24.2 scikit-learn==0.21.2 statsmodels==0.9.0 python-igraph==0.7.1 louvain==0.6.1

[LuckyMD]

Just linking to the anndata issue report for this: scverse/anndata#165

Also, in case anyone else sees this issue: It has to do with scipy version 1.2.1

[gokceneraslan]

I also encountered the same error, it is about subsetting but it has nothing to do with .uns, just to report. It's fixed when I upgrade to scipy 1.3.0.

[evafast]

I just encountered the same error. It only occurs on the subsetted matrices. However my scipy and statsmodels packages are current.

These are all the versions that I am using:
scanpy==1.4.4.post1 anndata==0.6.22.post1 umap==0.3.10 numpy==1.17.4 scipy==1.3.3 pandas==0.25.3 scikit-learn==0.21.3 statsmodels==0.10.2 python-igraph==0.7.1 louvain==0.6.1

I did try and downgrade to scipy==1.3.0 and statsmodels=0.10.0 but I am still getting the same error.

Any idea what else to try?

Thanks!!!!

Eva

[flying-sheep]

And your anndata2ri version?

[evafast]

Version for anndata2ri is 1.0.1

[LuckyMD]

I recently had an error with ordering of sparse matrix indices after subsetting. R doesn't like it if the indices aren't contiguous anymore, so it gives an error when converting. You could check if all of your sparse matrices have .has_sorted_index == True.

[evafast]

Just ran adata.raw.X.has_sorted_indices and it says True. However I'm still getting the error. Any other ideas what to try?

[flying-sheep]

At this point, we need a minimal reproducible example to see what’s happening. Just a code block that can be executed in a fresh session of python and triggers the issue in as little code as possible.

[evafast]

I think I have resolved the issue with a workaround that is probably more ideal anyways - doing everything in a Docker container. Seems to be an OS specific issue that occurs while saving a subsetted .h5ad file. Thanks for all your help.

These are my original commands (did not work):

this was done on Windows 8.1 OS

adata = sc.read('./write/previouslysaved_adata.h5ad')
adata[
    (adata.obs['sample']=='ct') |  
    (adata.obs['sample']=='treatment_1')
].write('/data/subset_adata.h5ad')

this was done in Docker container

adata = sc.read('./data/subset_adata.h5ad')

%%R -i adata
adata

This works:

all done in Docker container:

adata = sc.read('./data/previouslysaved_adata.h5ad')
adata[
    (adata.obs['sample']=='ct') |
    (adata.obs['sample']=='treatment_1')
].write('/data/subset_adata.h5ad')

adata = sc.read('./data/subset_adata.h5ad')

%%R -i adata
adata

[LuckyMD]

So this sounds like it is a windows issue, no?

[dburkhardt]

I'm finding this same error in a project I'm working on with @LuckyMD (https://github.com/singlecellopenproblems/SingleCellOpenProblems), also involving subsetting a sparse matrix

from openproblems.tasks.label_projection.datasets import  zebrafish_labels
import scIB.preprocessing

adata = zebrafish_labels(test=True) # This uses scanpy.pp.subsample()
scIB.preprocessing.normalize(adata)

raises

---------------------------------------------------------------------------
RRuntimeError                             Traceback (most recent call last)
<ipython-input-1-6c0cfb2c1423> in <module>
      3
      4 adata = zebrafish_labels(test=True) # This uses scanpy.pp.subsample()
----> 5 scIB.preprocessing.normalize(adata)

~/.local/lib/python3.8/site-packages/scIB/preprocessing.py in normalize(adata, min_mean)
    166     sc.tl.louvain(adata_pp, key_added='groups', resolution=0.5)
    167
--> 168     ro.globalenv['data_mat'] = adata.X.T
    169     ro.globalenv['input_groups'] = adata_pp.obs['groups']
    170     size_factors = ro.r(f'computeSumFactors(data_mat, clusters = input_groups, min.mean = {min_mean})')

/usr/lib/python3.8/site-packages/rpy2/robjects/environments.py in __setitem__(self, item, value)
     30
     31     def __setitem__(self, item: str, value: typing.Any) -> None:
---> 32         robj = conversion.converter.py2rpy(value)
     33         super(Environment, self).__setitem__(item, robj)
     34

/usr/lib/python3.8/functools.py in wrapper(*args, **kw)
    873                             '1 positional argument')
    874
--> 875         return dispatch(args[0].__class__)(*args, **kw)
    876
    877     funcname = getattr(func, '__name__', 'singledispatch function')

~/.local/lib/python3.8/site-packages/anndata2ri/scipy2ri/py2r.py in wrapper(obj)
     41
     42         with localconverter(default_converter + numpy2ri.converter):
---> 43             return f(obj)
     44
     45     return wrapper

~/.local/lib/python3.8/site-packages/anndata2ri/scipy2ri/py2r.py in csr_to_rmat(csr)
     63 def csr_to_rmat(csr: sparse.csr_matrix):
     64     t, conv_data, _ = get_type_conv(csr.dtype)
---> 65     return methods.new(
     66         f"{t}gRMatrix",
     67         j=as_integer(csr.indices),

/usr/lib/python3.8/site-packages/rpy2/robjects/functions.py in __call__(self, *args, **kwargs)
    195                 v = kwargs.pop(k)
    196                 kwargs[r_k] = v
--> 197         return (super(SignatureTranslatedFunction, self)
    198                 .__call__(*args, **kwargs))
    199

/usr/lib/python3.8/site-packages/rpy2/robjects/functions.py in __call__(self, *args, **kwargs)
    123             else:
    124                 new_kwargs[k] = conversion.py2rpy(v)
--> 125         res = super(Function, self).__call__(*new_args, **new_kwargs)
    126         res = conversion.rpy2py(res)
    127         return res

/usr/lib/python3.8/site-packages/rpy2/rinterface_lib/conversion.py in _(*args, **kwargs)
     42 def _cdata_res_to_rinterface(function):
     43     def _(*args, **kwargs):
---> 44         cdata = function(*args, **kwargs)
     45         # TODO: test cdata is of the expected CType
     46         return _cdata_to_rinterface(cdata)

/usr/lib/python3.8/site-packages/rpy2/rinterface.py in __call__(self, *args, **kwargs)
    622                     error_occured))
    623             if error_occured[0]:
--> 624                 raise embedded.RRuntimeError(_rinterface._geterrmessage())
    625         return res
    626

RRuntimeError: Error in validObject(.Object) :
  invalid class “dgRMatrix” object: slot j is not increasing inside a column

I have the following package versions installed

anndata 0.7.4
scanpy 1.4.5.post3
rpy2 3.3.5
anndata2ri 1.0.3
scanorama 1.6
memory_profiler 0.57.0
networkx 2.4
igraph 0.8.2
louvain 0.6.1
scipy 1.4.1
statsmodels 0.11.0rc1

[flying-sheep]

So adata.X.T.indices isn’t sorted. How does that come to be, is that legal, and can you fix it at its source or do we have to handle it here somehow?

[dburkhardt]

I think this is happening during the call to sc.pp.subsample. I added a call to np.sort in scverse/scanpy#1382 that fixes this issue

[flying-sheep]

Seems like that’s the correct thing to do. I think correctly constructed matrices do have the property R expects.

[LuckyMD]

So I think that this has to do with sparse matrices in R expecting contiguous indices. After subsetting, that is no longer the case. Thus, whenever I move a sparse index into R, I do the following:

if sparse.issparse(adata.X):
    if not adata.X.has_sorted_indices:
        adata.X.sort_indices()

That should solve the problem I think.

[flying-sheep]

Yeah, see #51 and #52

I think the newest anndata2ri should fix this

[LuckyMD]

Ah, nice! I didn't see the update. Does this fix affect passing back and forth between R and python? For example, if you have a subsetted python matrix, pass it to R and then back to python, would it still be seen as similar/equal?

[flying-sheep]

Comparison operations should test if the represented data is equal, not implementation details. Therefore I’d say that if they aren’t, that could be considered a bug in scipy’s spmatrix.__eq__ method.

[LuckyMD]

fair point.

[nroak]

I'm seeing this issue again with latest versions of scipy and anndata2ri

if issparse(adata.X):
    if not adata.X.has_sorted_indices:
        adata.X.sort_indices()
ro.globalenv["adata"] = adata

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
Cell In[49], line 4
      2     if not adata.X.has_sorted_indices:
      3         adata.X.sort_indices()
----> 4 ro.globalenv["adata"] = adata

File ~/opt/anaconda3/envs/velocyto/lib/python3.10/site-packages/rpy2/robjects/environments.py:35, in Environment.__setitem__(self, item, value)
     34 def __setitem__(self, item: str, value: typing.Any) -> None:
---> 35     robj = conversion.get_conversion().py2rpy(value)
     36     super(Environment, self).__setitem__(item, robj)

File ~/opt/anaconda3/envs/velocyto/lib/python3.10/functools.py:889, in singledispatch.<locals>.wrapper(*args, **kw)
    885 if not args:
    886     raise TypeError(f'{funcname} requires at least '
    887                     '1 positional argument')
--> 889 return dispatch(args[0].__class__)(*args, **kw)

File ~/opt/anaconda3/envs/velocyto/lib/python3.10/site-packages/anndata2ri/py2r.py:56, in py2rpy_anndata(obj)
     54 # TODO: sparse
     55 x = {} if obj.X is None else dict(X=mat_converter.py2rpy(obj.X.T))
---> 56 layers = {k: mat_converter.py2rpy(v.T) for k, v in obj.layers.items()}
     57 assays = ListVector({**x, **layers})
     59 row_args = {k: pandas2ri.py2rpy(v) for k, v in obj.var.items()}

File ~/opt/anaconda3/envs/velocyto/lib/python3.10/site-packages/anndata2ri/py2r.py:56, in <dictcomp>(.0)
     54 # TODO: sparse
     55 x = {} if obj.X is None else dict(X=mat_converter.py2rpy(obj.X.T))
---> 56 layers = {k: mat_converter.py2rpy(v.T) for k, v in obj.layers.items()}
     57 assays = ListVector({**x, **layers})
     59 row_args = {k: pandas2ri.py2rpy(v) for k, v in obj.var.items()}

File ~/opt/anaconda3/envs/velocyto/lib/python3.10/functools.py:889, in singledispatch.<locals>.wrapper(*args, **kw)
    885 if not args:
    886     raise TypeError(f'{funcname} requires at least '
    887                     '1 positional argument')
--> 889 return dispatch(args[0].__class__)(*args, **kw)

File ~/opt/anaconda3/envs/velocyto/lib/python3.10/site-packages/anndata2ri/scipy2ri/py2r.py:88, in py2r_context.<locals>.wrapper(obj)
     36     importr('Matrix')  # make class available
     37     matrix = SignatureTranslatedAnonymousPackage(
     38         """
     39         sparse_matrix <- function(x, conv_data, dims, ...) {
   (...)
     85         'matrix',
     86     )
---> 88 return f(obj)

File ~/opt/anaconda3/envs/velocyto/lib/python3.10/site-packages/anndata2ri/scipy2ri/py2r.py:97, in csc_to_rmat(csc)
     93 @converter.py2rpy.register(sparse.csc_matrix)
     94 @py2r_context
     95 def csc_to_rmat(csc: sparse.csc_matrix):
     96     csc.sort_indices()
---> 97     conv_data = get_type_conv(csc.dtype)
     98     with localconverter(default_converter + numpy2ri.converter):
     99         return matrix.from_csc(i=csc.indices, p=csc.indptr, x=csc.data, dims=list(csc.shape), conv_data=conv_data)

File ~/opt/anaconda3/envs/velocyto/lib/python3.10/site-packages/anndata2ri/scipy2ri/py2r.py:28, in get_type_conv(dtype)
     26     return base.as_logical
     27 else:
---> 28     raise ValueError(f'Unknown dtype {dtype!r} cannot be converted to ?gRMatrix.')

ValueError: Unknown dtype dtype('uint16') cannot be converted to ?gRMatrix.

anndata2ri==1.1, scipy==1.10.0, rpy2==3.5.8, scanpy==1.9.1