Scanpy: sc.tl.umap numba error when used with init_pos="paga"

[jenzopr]

Hi Alex,

UMAP throws an error if I use scanpy.tl.ump with initial positions from sc.tl.paga. Based on the error (see below) I thought it was a problem of UMAP itself. However, the error is not thrown when called without initial positions from paga.

Here is the output / error:

sc.tl.umap(adata, init_pos='paga')

computing UMAP
    using 'X_pca' with n_pcs = 50

---------------------------------------------------------------------------
TypingError                               Traceback (most recent call last)
<ipython-input-35-924452b37e5b> in <module>
----> 1 sc.tl.umap(adata, init_pos='paga')

/opt/conda/lib/python3.7/site-packages/scanpy/tools/_umap.py in umap(adata, min_dist, spread, n_components, maxiter, alpha, gamma, negative_sample_rate, init_pos, random_state, a, b, copy)
    137         neigh_params.get('metric', 'euclidean'),
    138         neigh_params.get('metric_kwds', {}),
--> 139         verbose=max(0, verbosity-3))
    140     adata.obsm['X_umap'] = X_umap  # annotate samples with UMAP coordinates
    141     logg.info('    finished', time=True, end=' ' if _settings_verbosity_greater_or_equal_than(3) else '\n')

/opt/conda/lib/python3.7/site-packages/umap/umap_.py in simplicial_set_embedding(data, graph, n_components, initial_alpha, a, b, gamma, negative_sample_rate, n_epochs, init, random_state, metric, metric_kwds, verbose)
    984         initial_alpha,
    985         negative_sample_rate,
--> 986         verbose=verbose,
    987     )
    988 

/opt/conda/lib/python3.7/site-packages/numba/dispatcher.py in _compile_for_args(self, *args, **kws)
    348                 e.patch_message(msg)
    349 
--> 350             error_rewrite(e, 'typing')
    351         except errors.UnsupportedError as e:
    352             # Something unsupported is present in the user code, add help info

/opt/conda/lib/python3.7/site-packages/numba/dispatcher.py in error_rewrite(e, issue_type)
    315                 raise e
    316             else:
--> 317                 reraise(type(e), e, None)
    318 
    319         argtypes = []

/opt/conda/lib/python3.7/site-packages/numba/six.py in reraise(tp, value, tb)
    656             value = tp()
    657         if value.__traceback__ is not tb:
--> 658             raise value.with_traceback(tb)
    659         raise value
    660 

TypingError: Failed in nopython mode pipeline (step: nopython frontend)
Invalid use of type(CPUDispatcher(<function rdist at 0x7fb3c827f840>)) with parameters (array(float64, 1d, C), array(float64, 1d, C))
Known signatures:
 * (array(float32, 1d, A), array(float32, 1d, A)) -> float32
 * parameterized
[1] During: resolving callee type: type(CPUDispatcher(<function rdist at 0x7fb3c827f840>))
[2] During: typing of call at /opt/conda/lib/python3.7/site-packages/umap/umap_.py (776)


File "../../../opt/conda/lib/python3.7/site-packages/umap/umap_.py", line 776:
def optimize_layout(
    <source elided>

                dist_squared = rdist(current, other)
                ^

This is not usually a problem with Numba itself but instead often caused by
the use of unsupported features or an issue in resolving types.

To see Python/NumPy features supported by the latest release of Numba visit:
http://numba.pydata.org/numba-doc/dev/reference/pysupported.html
and
http://numba.pydata.org/numba-doc/dev/reference/numpysupported.html

For more information about typing errors and how to debug them visit:
http://numba.pydata.org/numba-doc/latest/user/troubleshoot.html#my-code-doesn-t-compile

If you think your code should work with Numba, please report the error message
and traceback, along with a minimal reproducer at:
https://github.com/numba/numba/issues/new

What I basically do from raw UMI counts:

1. total counts normalization / logarithmization
2. PCA, bbknn, louvain
3. combat, HVG, PCA, UMAP (works well)
4. Paga (with louvain from 2., works well)
5. UMAP (with positions from 4., does not work)

Any idea? Any further info needed?
Best,
Jens

[jenzopr]

Issue #666 👹 🙈

[falexwolf]

Hi! Are you beyond 1.4.2, I solved a bug around this here (https://github.com/theislab/scanpy/commit/ec2d12608dae77c68811d7742075079664eb0e4c)?

[jenzopr]

This is the output of sc.logging.print_versions():

scanpy==1.4.3 anndata==0.6.20 umap==0.3.8 numpy==1.16.3 scipy==1.2.1 pandas==0.24.2 scikit-learn==0.21.2 statsmodels==0.9.0 python-igraph==0.7.1 louvain==0.6.1

Probably downgrading of scikit-learn might help?

[jenzopr]

Downgrading scikit-learn to 0.20.3 does not resolve the error...

[LuckyMD]

What I basically do from raw UMI counts:
1. total counts normalization / logarithmization
2. PCA, bbknn, louvain
3. combat, HVG, PCA, UMAP (works well)

Just out of curiousity, you use both BBKNN and combat? Does Louvain after ComBat, HVG, and PCA not work as well for you? It's interesting that you go with two different knn graphs for clustering and visualization.

[flying-sheep]

The error seems to be that rdist is called from umap.umap_.optimize_layout with float64 arrays while it can only handle float32 arrays.

There seem to have been a few changes in umap between 0.3.8 and 0.3.9 maybe you should try 0.3.9.

[jenzopr]

Just out of curiousity, you use both BBKNN and combat? Does Louvain after ComBat, HVG, and PCA not work as well for you? It's interesting that you go with two different knn graphs for clustering and visualization.

@LuckyMD I found that the clustering using the bbknn kNN graph is much cleaner on UMAP, compared to e.g. sc.pp.neighbors. From combat, I just obtain the adjusted data, not a kNN graph.

There seem to have been a few changes in umap between 0.3.8 and 0.3.9 maybe you should try 0.3.9.

@flying-sheep Thanks! With scikit-learn pinned to 0.20.3 the umap version was updated to 0.3.9 (I use a container and rebuilt it).

I guess I will try to create a reproducible example and open an issue in umap.
Edit: https://github.com/lmcinnes/umap/issues/179

[flying-sheep]

Great, thank you. Seems like the issue is well placed over there.

[LuckyMD]

@LuckyMD I found that the clustering using the bbknn kNN graph is much cleaner on UMAP, compared to e.g. sc.pp.neighbors. From combat, I just obtain the adjusted data, not a kNN graph.

I see. So you don't recalculate sc.pp.neighbors after the combat run? Otherwise your umap is not using the combat corrected counts at all, but instead the KNN graph from bbknn as far as I understand. Then you are basically using bbknn for the embedding, and combat only for the visualization of expression values in your plots. Do your batches have similar cell type compositions?

[jenzopr]

Exactly, I do not recalculate sc.pp.neighbors after bbknn. In my understanding, UMAP uses the corrected counts via the preceding PCA step?

Yes, my batches have very similar (if not the same) cell type composition.

[LuckyMD]

I think umap works on the connectivities matrix generated from sc.pp.neighbors. I guess you can test this by omitting sc.pp.neighbors before running umap. Or just running umap after your step 2 and look at the difference.

Yes, my batches have very similar (if not the same) cell type composition.

In that case it would make sense for both bbknn and combat to work quite well.

[jenzopr]

@falexwolf @flying-sheep FYI it seems that sc.pp.combat returns a dtype('float64'). When this ends up in adata.X, problems occur during sc.tl.umap with init_pos. Specifically, because the type of get_init_pos_from_paga is set to the type from adata.X here.

I could work around the problem using

# Correct batch effect
adjusted_data = sc.pp.combat(adata, key = "Batch", inplace=False)

# Add adjusted data as data
adata.X = adjusted_data.astype(adata.X.dtype)

[jenzopr]

@LuckyMD Indeed, UMAP does not use information from X_pca when bbknn was run before, although it states so:

computing UMAP
    using 'X_pca' with n_pcs = 50
    finished (0:00:19.09) --> added
    'X_umap', UMAP coordinates (adata.obsm)

Any idea why?
Then my setup uses combat batch correction only for visualization of expression values.

[LuckyMD]

@jenzopr I think this is legacy from when sc.pp.neighbors was not a separate function, but sc.tl.umap instead had options to recompute the neighborhood graph. This should probably be changed. I don't know where the info about the 50 pcs comes from... is this stored somewhere related to the connectivities matrix? Do you know @flying-sheep?

Also, I wonder if the clustering looks better for you only because you do both the umap and louvain with bbknn. It should also look good if both is done via combat in my opinion. Would be curious about this, as I recommend ComBat in my best-practices tutorial (based on several papers, like kBet).

[jenzopr]

@LuckyMD I can put together some slides for you :-)

[falexwolf]

Although it states so:

UMAP only uses the representation of a data matrix for determining the number of connected components of the graph for the init conditions [if these aren't explicitly defined (they are if choosing init_pos='paga'): https://github.com/lmcinnes/umap/blob/948f60ff0caf7ccef0ab68626c7b99a11e66f1bb/umap/umap_.py#L958-L965

In fact, the only place where it enters is for the computation of the mean positions of the disconnected components: https://github.com/lmcinnes/umap/blob/948f60ff0caf7ccef0ab68626c7b99a11e66f1bb/umap/spectral.py#L50

Implementation-wise, it's a bit unfortunate that the data matrix is carried through all these functions just for that reason... But it's not a problem for the results.

The confusing logging is fixed via
https://github.com/theislab/scanpy/commit/a5bd1ecd8ab04ec79369f60d3656f578a4cde40c

[LuckyMD]

@LuckyMD I can put together some slides for you :-)

that would be awesome :)