Problem at reproducibility of UMAP / leiden

Hi,

I am working on a project with a labmate and we are using the same dataset. We have found that, when running the same pipeline on the same adata the neighbors / bbknn + UMAP + leiden results, even with the same seed, the clustering solution and UMAP are considerably different. This renders the analysis unreproducible and makes the downstream analysis far more difficult to do, since I have to map my clustering solutions and UMAP plots with hers using markers, and it is quite impractical.

We have the same versions of scanpy, leiden, umap, and bbknn on the two computers:

• scanpy==1.4.5.post2
• umap-learn==0.3.10
• leidenalg==0.7.0
• bbknn==1.3.6

To try to reproduce the issue, we have created a random matrix with the same seed (10), and create one annData with sc.pp.neighbours, and another one with bbknn. We have made the adatas to have two batches, so that we can use bbknn.

seed = 10
np.random.seed(seed)
a = np.random.rand(100, 100)
b = np.random.rand(100, 100)
print(np.sum(a), np.sum(b))

adata = sc.AnnData.concatenate(sc.AnnData(X=a), sc.AnnData(X=b), batch_categories=['a', 'b'])
sc.tl.pca(adata)
sce.pp.bbknn(adata, metric='angular')
sc.tl.umap(adata, random_state=seed)
sc.tl.leiden(adata, resolution=0.5, random_state=seed)
sc.pl.umap(adata, color=['batch', 'leiden'], alpha=0.3)
print(adata.uns['neighbors']['connectivities'].sum())

adata_neigh = adata.copy()
sc.pp.neighbors(adata_neigh, metric='cosine', random_state=seed)
sc.tl.umap(adata_neigh, random_state=seed)
sc.tl.leiden(adata_neigh, resolution=0.6, random_state=seed)
sc.pl.umap(adata_neigh, color=['batch', 'leiden'], alpha=0.3)
print(adata_neigh.uns['neighbors']['connectivities'].sum())

Our matrices are the same (the sums are 4918.370372081173 and 5005.088472351332), so the random generation works, but then the UMAPs and clustering solutions are different.

For the adata run with sc.pp.neighbors (left are batches and right are leiden cluster labels):
Mine

Hers

For the adata run with sce.pp.bbknn:
Mine

Hers

Our PCA decomposition has the same coordinates, so we discard the PCA as the source of variability. Also, since both UMAP and leiden look different, we think the source might come from the neighbor calculation.

When running adata.uns['neighbors']['connectivities'].sum() I get 801.5580058219996 and 1204.5274490986717 for adata and adata_neigh. I don't have her values now, but the values using the real dataset were in the order of 5000, and they were of by less than 0.001; so we are confused that with such a small difference on the sum, the results can be so different.

I attach the adatas for you to inspect them if you need more info.

adatas.zip

Thanks for the help!

Start with the reproduction script in the issue and the attached adatas.zip, using the listed scanpy, umap-learn, leidenalg, and bbknn versions. Compare the neighbor connectivities and subsequent UMAP and leiden outputs between environments; done would require identifying why identical inputs and seeds diverge and documenting a reproducible resolution.