Leiden at scale: the same clusters on any number of threads

piaso.tl.leiden clusters the neighbour graph with PIASO’s own Leiden. It is the algorithm of Traag, Waltman and van Eck (2019), with the same objective and the same randomised refinement as igraph’s, and it reaches the same modularity. Two things are different: it uses every core it is given, and the labels do not depend on how many that is.

This page clusters a 200,061-cell mouse developing visual cortex atlas, shows the clusters beside igraph’s, and reports what it costs. Timings are from a 10-core workstation (20 threads).

Cluster

import piaso

path = "allen_devvis_rna.cytome"

piaso.tl.infog(path, n_top_genes=3000, save_layer=True)
piaso.tl.runSVD(path, layer="infog", n_components=50, key_added="X_svd")
piaso.tl.neighbors(path, use_rep="X_svd", n_neighbors=15, key_added="svd")

piaso.tl.leiden(path, neighbors_key="svd", key_added="leiden")

piaso.tl.umap(path, use_rep="X_svd", key_added="X_umap", neighbors_key="svd")

The call is the one you know. On an AnnData it is the same: piaso.tl.leiden(adata, resolution=1.0) reads adata.obsp["connectivities"] and writes adata.obs["leiden"].

The graph here has 200,061 cells and 4.1 million entries. Clustering it takes 2.7 s. On the store the call also reads the graph and writes the labels, 3.7 s in all.

The number of threads changes the time, not the labels

import numpy as np
import cytome

ds = cytome.open(path)
graph = {"connectivities": ds.graphs["svd_connectivities"].to_sparse()}
ds.close()

labels = {t: piaso.tl.leiden(None, knn_result=graph, n_threads=t) for t in (1, 2, 4, 10, 20)}
all(np.array_equal(labels[1], labels[t]) for t in labels)

True

threads	1	2	4	10	20
time	5.4 s	4.2 s	3.4 s	2.8 s	2.8 s

One graph and one random_state give one partition: on a laptop, on a workstation, on a cluster node with whatever allocation the job was given. A parallel clustering usually gives that up, because its threads see each other’s moves in an order that changes from run to run, and two runs of the same program on the same data then disagree. Here they cannot.

n_threads=0, the default, uses all the cores the process may use; on a cluster that is the job’s allocation. Since the labels do not depend on it, there is nothing to record about it in your methods.

What it looks like, beside igraph’s

piaso.tl.leiden(path, neighbors_key="svd", key_added="leiden_igraph", backend="igraph")

Two partitions number their clusters independently, so cluster 7 of one is not cluster 7 of the other. To compare them by eye, igraph’s clusters are given the number of the PIASO cluster they share most cells with, where the two are the same cluster in the main, and a new number otherwise:

import pandas as pd

ds = cytome.open(path)
ours = np.asarray(ds.cells["leiden"]).astype(str)
theirs = np.asarray(ds.cells["leiden_igraph"]).astype(str)

ct = pd.crosstab(theirs, ours)
jaccard = ct / (ct.sum(1).to_numpy()[:, None] + ct.sum(0).to_numpy()[None, :] - ct)
best = jaccard.idxmax(1)
new = iter(str(k) for k in range(ct.shape[1], 10_000))
name = {c: (best[c] if jaccard.loc[c, best[c]] > 0.5 and jaccard[best[c]].idxmax() == c else next(new))
        for c in ct.index}

ds.cells["leiden_igraph_named"] = np.array([name[c] for c in theirs])
ds.flush()
ds.set_categories("leiden_igraph_named", order=sorted(set(name.values()), key=int))
ds.close()

# one colour for one number, in both panels
numbers = sorted(set(ours) | set(name.values()), key=int)
colours = dict(zip(numbers, piaso.pl.color.d_color4 * 3))

import matplotlib.pyplot as plt

piaso.settings.set_figure_params(style="cell")
fig, axes = plt.subplots(1, 3, figsize=(18, 6))
for ax, key, title in zip(axes, ("subclass_label", "leiden", "leiden_igraph_named"),
                          ("Annotated subclass", "piaso.tl.leiden", "igraph")):
    piaso.pl.plotEmbedding(path, basis="X_umap", color=key, ax=ax, show=False,
                           palette=None if key == "subclass_label" else colours,
                           legend_loc="on_data", legend_fontsize=6, frameon=False)
    ax.set_title(title)

Both give 42 clusters, and 35 of igraph’s are a PIASO cluster in the main: same number, same colour, same place. The two differ where a continuum has to be cut somewhere, in the immature neurons at the centre of the map (clusters 1, 7, 18 and 21 on the left; 42, 45, 46 and 47 on the right). Nothing in the data says where those cuts belong, and another seed of either program moves them as well.

The Sankey diagram shows every cell, not a picture of them:

piaso.pl.sankey(path, left="leiden", right="leiden_igraph_named", palette=colours)

Most clusters pass straight across. The thin bands that cross are the continuum again.

And against the annotation of the atlas:

piaso.pl.sankey(path, left="subclass_label", right="leiden", flow_scale="sqrt")

The same quality

	piaso.tl.leiden	igraph
clusters	42	42
agreement with the annotated subclasses (ARI)	0.354	0.352
two seeds of the same program agree (ARI, median of 5 seeds)	0.80	0.81
PIASO against igraph (ARI, median)	0.81

The two programs differ from each other exactly as much as each differs from itself with another seed. On this graph, at this resolution, that is what Leiden’s partitions are: any two of them agree at about 0.8.

What it costs

10 iterations	cells	piaso.tl.leiden	igraph
this dataset	200,061	2.7 s	17.3 s
whole mouse brain	2,341,350	26 s, 1.3 GB	554 s, 3.7 GB

Memory is what the call adds to the process. The whole-brain graph has 54 million entries.

n_iterations

An iteration is one run of the whole algorithm, starting from the partition the run before it ended with. The default is 10.

	2 iterations	10 iterations
this dataset	1.4 s	2.7 s
whole mouse brain	12 s	26 s

Two iterations give 98 % of the partition of ten here (ARI 0.978). Ten raise the modularity a little and make two seeds agree a little more. Keep ten unless you are clustering tens of millions of cells.

Reproducing a result of PIASO 1.2.5 or earlier

Before 1.2.6 piaso.tl.leiden ran igraph. Those partitions are still there:

piaso.tl.leiden(adata, backend="igraph")
piaso.tl.leiden_local(adata, groupby="leiden", backend="igraph")

leiden_local and runGDR cluster with piaso.tl.leiden, so their results follow the same change.

Good to know

• The graph must be symmetric. The connectivities of piaso.tl.neighbors are. A directed graph is refused, and the message names a pair of nodes whose two directions differ.
• random_state=None refines without any random choice.
• Ctrl-C stops a run.
• On a shared machine, n_threads=8 limits one call, and the environment variable RAYON_NUM_THREADS=8 limits every call.
• A graph of fewer than 32,768 cells is clustered on one thread whatever is asked, which is faster at that size.