Skip to content

Graphs and clustering

find_neighbors builds both graphs Seurat builds — the directed KNN graph and the shared-nearest-neighbour graph on top of it — and find_clusters runs community detection over the SNN.

Four details here were wrong before the integration tutorial went looking, and each mattered downstream: the KNN graph is stored directed rather than symmetrized, the SNN keeps its diagonal and is computed in float64, run_umap zeroes the diagonal when handed a graph, and singletons are folded into their nearest community the way GroupSingletons does. They are noted in the docstrings because each one changes cluster assignments, not just internals.

Neighbour graphs

find_neighbors 

find_neighbors(seurat, dims: Optional[Union[list[int], range]] = None, k_param: int = 20, assay: Optional[str] = None, reduction: str = 'pca', graph_name: Optional[str] = None, nn_name: Optional[str] = None, prune_snn: float = 1 / 15, seed: int = 42) -> None

Build KNN and SNN graphs from a low-dimensional embedding.

Mirrors R's FindNeighbors(pbmc, dims = 1:10). Stores Graph objects in seurat.graphs[graph_name + '_nn'] and seurat.graphs[graph_name + '_snn'].

Parameters:

  • dims (Optional[Union[list[int], range]], default: None ) –

    which PCs to use (0-indexed; default all available)

  • k_param (int, default: 20 ) –

    number of nearest neighbors

  • reduction (str, default: 'pca' ) –

    which reduction to use ('pca' by default)

  • graph_name (Optional[str], default: None ) –

    prefix for graph names (defaults to active assay name)

  • prune_snn (float, default: 1 / 15 ) –

    edges with Jaccard index below this are pruned (Seurat default 1/15)

Source code in truecell/neighbors.py 
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
def find_neighbors(
    seurat,
    dims: Optional[Union[list[int], range]] = None,
    k_param: int = 20,
    assay: Optional[str] = None,
    reduction: str = "pca",
    graph_name: Optional[str] = None,
    nn_name: Optional[str] = None,
    prune_snn: float = 1 / 15,
    seed: int = 42,
) -> None:
    """Build KNN and SNN graphs from a low-dimensional embedding.

    Mirrors R's FindNeighbors(pbmc, dims = 1:10).
    Stores Graph objects in seurat.graphs[graph_name + '_nn'] and
    seurat.graphs[graph_name + '_snn'].

    Parameters
    ----------
    dims        : which PCs to use (0-indexed; default all available)
    k_param     : number of nearest neighbors
    reduction   : which reduction to use ('pca' by default)
    graph_name  : prefix for graph names (defaults to active assay name)
    prune_snn   : edges with Jaccard index below this are pruned (Seurat default 1/15)
    """
    assay_name = assay or seurat.active_assay

    # Get embeddings
    if reduction not in seurat.reductions:
        raise KeyError(f"Reduction '{reduction}' not found. Run run_pca() first.")
    dr = seurat.reductions[reduction]
    embeddings = dr.cell_embeddings  # (cells × dims)

    if dims is None:
        emb = embeddings
    else:
        dims_list = list(dims)
        emb = embeddings[:, dims_list]

    cells = seurat.cell_names()
    n_cells = len(cells)

    # Build KNN
    nn_idx, nn_dist = _build_knn(emb, k_param, seed)

    # Build KNN sparse graph (symmetric)
    knn_mat = _knn_to_sparse(nn_idx, n_cells)

    # Build SNN (shared nearest neighbor) sparse graph with Jaccard weights
    snn_mat = _build_snn(nn_idx, n_cells, k_param, prune_snn)

    prefix = graph_name or assay_name
    knn_name = f"{prefix}_nn"
    snn_name = f"{prefix}_snn"

    seurat.graphs[knn_name] = Graph(
        matrix=knn_mat, cell_names=cells, assay_used=assay_name
    )
    seurat.graphs[snn_name] = Graph(
        matrix=snn_mat, cell_names=cells, assay_used=assay_name
    )
    log_truecell_command(
        seurat, "FindNeighbors", assay=assay_name, reduction=reduction,
        params={"k_param": k_param, "prune_snn": prune_snn,
                "dims": list(dims) if dims is not None else None},
    )

find_multi_modal_neighbors 

find_multi_modal_neighbors(seurat, reduction_list: Sequence[str] = ('pca', 'apca'), dims_list: Optional[Sequence[Sequence[int]]] = None, k_nn: int = 20, l2_norm: bool = True, knn_graph_name: str = 'wknn', snn_graph_name: str = 'wsnn', knn_range: int = 200, prune_snn: float = 1 / 15, sd_scale: float = 1.0, cross_constant: Optional[float] = None, smooth: bool = False, seed: int = 42) -> None

Compute weighted-nearest-neighbour graphs across modalities.

Mirrors R's FindMultiModalNeighbors(obj, reduction.list = list("pca","apca"), dims.list = list(1:30, 1:18)). Stores:

  • seurat.graphs[knn_graph_name] — joint KNN graph
  • seurat.graphs[snn_graph_name] — joint SNN graph (feed to find_clusters(graph_name=...) or run_umap(graph=...))
  • seurat.meta_data["<modality>.weight"] — per-cell modality weights

Parameters:

  • reduction_list (Sequence[str], default: ('pca', 'apca') ) –

    reductions to combine, one per modality

  • dims_list (Optional[Sequence[Sequence[int]]], default: None ) –

    dims (0-indexed) to use from each reduction; default all

  • k_nn (int, default: 20 ) –

    number of joint neighbours to keep per cell

  • l2_norm (bool, default: True ) –

    L2-normalise each embedding first (R's l2.norm)

  • knn_range (int, default: 200 ) –

    candidate neighbours each modality nominates before the joint re-ranking (R's knn.range)

  • prune_snn (float, default: 1 / 15 ) –

    Jaccard prune threshold for the joint SNN graph

  • sd_scale (float, default: 1.0 ) –

    scaling on the per-cell kernel bandwidth (R's sd.scale)

  • cross_constant (Optional[float], default: None ) –

    denominator guard in the modality score; default 1e-4

  • smooth (bool, default: False ) –

    average each cell's modality score over its neighbours

Notes

The weights are stored under each reduction's assay_used name, so an RNA "pca" and an ADT "apca" produce RNA.weight and ADT.weight — the same columns Seurat writes, and readable by the plotting functions as if they were features.

Source code in truecell/multimodal.py 
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
def find_multi_modal_neighbors(
    seurat,
    reduction_list: Sequence[str] = ("pca", "apca"),
    dims_list: Optional[Sequence[Sequence[int]]] = None,
    k_nn: int = 20,
    l2_norm: bool = True,
    knn_graph_name: str = "wknn",
    snn_graph_name: str = "wsnn",
    knn_range: int = 200,
    prune_snn: float = 1 / 15,
    sd_scale: float = 1.0,
    cross_constant: Optional[float] = None,
    smooth: bool = False,
    seed: int = 42,
) -> None:
    """Compute weighted-nearest-neighbour graphs across modalities.

    Mirrors R's ``FindMultiModalNeighbors(obj, reduction.list =
    list("pca","apca"), dims.list = list(1:30, 1:18))``. Stores:

    * ``seurat.graphs[knn_graph_name]`` — joint KNN graph
    * ``seurat.graphs[snn_graph_name]`` — joint SNN graph (feed to
      ``find_clusters(graph_name=...)`` or ``run_umap(graph=...)``)
    * ``seurat.meta_data["<modality>.weight"]`` — per-cell modality weights

    Parameters
    ----------
    reduction_list : reductions to combine, one per modality
    dims_list      : dims (0-indexed) to use from each reduction; default all
    k_nn           : number of joint neighbours to keep per cell
    l2_norm        : L2-normalise each embedding first (R's ``l2.norm``)
    knn_range      : candidate neighbours each modality nominates before the
                     joint re-ranking (R's ``knn.range``)
    prune_snn      : Jaccard prune threshold for the joint SNN graph
    sd_scale       : scaling on the per-cell kernel bandwidth (R's ``sd.scale``)
    cross_constant : denominator guard in the modality score; default 1e-4
    smooth         : average each cell's modality score over its neighbours

    Notes
    -----
    The weights are stored under each reduction's ``assay_used`` name, so an RNA
    ``"pca"`` and an ADT ``"apca"`` produce ``RNA.weight`` and ``ADT.weight`` —
    the same columns Seurat writes, and readable by the plotting functions as if
    they were features.
    """
    if len(reduction_list) < 2:
        raise ValueError("find_multi_modal_neighbors needs at least 2 reductions.")

    for r in reduction_list:
        if r not in seurat.reductions:
            raise KeyError(
                f"Reduction '{r}' not found. Compute it before WNN "
                "(e.g. run_pca for RNA, run_pca(reduction_name='apca') for ADT)."
            )

    cells = seurat.cell_names()
    n_cells = len(cells)
    cross_constant = _CROSS_CONSTANT if cross_constant is None else cross_constant

    if k_nn >= n_cells:
        raise ValueError(
            f"k_nn ({k_nn}) must be smaller than the number of cells ({n_cells})."
        )

    # Per-modality embeddings (optionally dim-subset, then L2-normalised).
    embs: list[np.ndarray] = []
    for m, r in enumerate(reduction_list):
        e = seurat.reductions[r].cell_embeddings
        if dims_list is not None and dims_list[m] is not None:
            e = e[:, list(dims_list[m])]
        e = np.asarray(e, dtype=float)
        embs.append(_l2_norm(e) if l2_norm else e)

    # Stage 1 — per-cell modality weights, plus the bandwidths and nearest
    # distances stage 2 reuses (R threads these through ModalityWeights@params).
    weights, sigmas, nearest_dist = _modality_weights(
        embs, k_nn, sd_scale, cross_constant, smooth, seed,
    )

    # Stage 2 — joint neighbour search, then graphs off that single ranking.
    select_nn = _multi_modal_nn(
        embs, weights, sigmas, nearest_dist, k_nn, knn_range, seed,
    )

    wknn = _knn_union_graph(select_nn, n_cells)
    wsnn = _build_snn(select_nn, n_cells, select_nn.shape[1], prune_snn).tocsc()

    assay_name = seurat.active_assay
    seurat.graphs[knn_graph_name] = Graph(matrix=wknn, cell_names=cells, assay_used=assay_name)
    seurat.graphs[snn_graph_name] = Graph(matrix=wsnn, cell_names=cells, assay_used=assay_name)

    # Store per-cell modality weights in meta_data, named by each modality's assay.
    for m, r in enumerate(reduction_list):
        assay_used = seurat.reductions[r].assay_used or r
        seurat.meta_data[f"{assay_used}.weight"] = weights[:, m]

Community detection

find_clusters 

find_clusters(seurat, resolution: Union[float, Sequence[float]] = 0.5, algorithm: int = 1, graph_name: Optional[str] = None, random_seed: int = 0, n_iterations: int = -1, group_singletons: bool = True, cluster_name: Optional[Union[str, Sequence[str]]] = None) -> None

Apply Louvain or Leiden clustering on the SNN graph.

Mirrors R's FindClusters(pbmc, resolution = 0.5), including its vector form FindClusters(pbmc, resolution = c(0.4, 0.8, 1.2)).

Each resolution is written to its own metadata column, named {graph_name}_res.{resolution} as Seurat names it. seurat_clusters and the active identities are set from the last resolution in the sequence — last as given, not largest, which is what Seurat does.

Parameters:

  • resolution (Union[float, Sequence[float]], default: 0.5 ) –

    higher values give more / finer clusters. A sequence runs each in turn, as R's resolution = c(...) does.

  • algorithm (int, default: 1 ) –

    1 = Louvain, 2 = Louvain (multilevel, igraph's default), 4 = Leiden. (3 = SLM is not implemented.)

  • graph_name (Optional[str], default: None ) –

    SNN graph to use (defaults to '{assay}_snn')

  • random_seed (int, default: 0 ) –

    for reproducibility

  • n_iterations (int, default: -1 ) –

    Leiden iterations (-1 = until stable)

  • group_singletons (bool, default: True ) –

    absorb size-1 clusters into their best-connected neighbour, as Seurat's GroupSingletons does. With False they are all pooled into one "singleton" cluster instead — again matching Seurat.

  • cluster_name (Optional[Union[str, Sequence[str]]], default: None ) –

    override the generated column name(s); one name, or one per resolution. Seurat's cluster.name. seurat_clusters is still written either way.

Notes

Seurat runs its own modularity optimiser with n.start = 10 restarts and keeps the highest-modularity partition; this runs a single pass of igraph's multilevel Louvain. On the same graph that makes truecell's partition land in a slightly shallower optimum — measurably so, but not necessarily a worse one. See the clustering section of tutorials/integration_vignette.md.

Each resolution is clustered from the same seed rather than from a running RNG stream, so a partition does not depend on which resolutions preceded it or on the order they were given in. Verified against Seurat 5.5.1, where resolution 0.8 gives the same partition alone, in c(0.4, 0.8, 1.2), and in c(1.2, 0.8, 0.4).

Source code in truecell/clustering.py 
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
def find_clusters(
    seurat,
    resolution: Union[float, Sequence[float]] = 0.5,
    algorithm: int = 1,
    graph_name: Optional[str] = None,
    random_seed: int = 0,
    n_iterations: int = -1,
    group_singletons: bool = True,
    cluster_name: Optional[Union[str, Sequence[str]]] = None,
) -> None:
    """Apply Louvain or Leiden clustering on the SNN graph.

    Mirrors R's ``FindClusters(pbmc, resolution = 0.5)``, including its
    vector form ``FindClusters(pbmc, resolution = c(0.4, 0.8, 1.2))``.

    Each resolution is written to its own metadata column, named
    ``{graph_name}_res.{resolution}`` as Seurat names it. ``seurat_clusters`` and
    the active identities are set from the **last** resolution in the sequence —
    last as given, not largest, which is what Seurat does.

    Parameters
    ----------
    resolution   : higher values give more / finer clusters. A sequence runs each
                   in turn, as R's ``resolution = c(...)`` does.
    algorithm    : 1 = Louvain, 2 = Louvain (multilevel, igraph's default),
                   4 = Leiden. (3 = SLM is not implemented.)
    graph_name   : SNN graph to use (defaults to '{assay}_snn')
    random_seed  : for reproducibility
    n_iterations : Leiden iterations (-1 = until stable)
    group_singletons : absorb size-1 clusters into their best-connected
                   neighbour, as Seurat's ``GroupSingletons`` does. With
                   ``False`` they are all pooled into one ``"singleton"``
                   cluster instead — again matching Seurat.
    cluster_name : override the generated column name(s); one name, or one per
                   resolution. Seurat's ``cluster.name``. ``seurat_clusters`` is
                   still written either way.

    Notes
    -----
    Seurat runs its own modularity optimiser with ``n.start = 10`` restarts and
    keeps the highest-modularity partition; this runs a single pass of igraph's
    multilevel Louvain. On the same graph that makes truecell's partition land in
    a slightly shallower optimum — measurably so, but not necessarily a worse
    one. See the clustering section of ``tutorials/integration_vignette.md``.

    Each resolution is clustered from the same seed rather than from a running
    RNG stream, so a partition does not depend on which resolutions preceded it
    or on the order they were given in. Verified against Seurat 5.5.1, where
    resolution 0.8 gives the same partition alone, in ``c(0.4, 0.8, 1.2)``, and
    in ``c(1.2, 0.8, 0.4)``.
    """
    assay_name = seurat.active_assay
    if graph_name is None:
        graph_name = f"{assay_name}_snn"
        if graph_name not in seurat.graphs:
            # Try knn graph
            graph_name = f"{assay_name}_nn"

    if graph_name not in seurat.graphs:
        raise KeyError(
            f"Graph '{graph_name}' not found. Run find_neighbors() first."
        )

    graph = seurat.graphs[graph_name]
    mat = graph._matrix  # scipy sparse (cells × cells)

    # Validated once rather than per iteration. This does *not* prevent a partial
    # write — the dispatch sits at the top of the loop body, so a bad `algorithm`
    # would raise on the first resolution either way, before anything is stored.
    # It is here so the check does not depend on the loop at all.
    #
    # A partial write is still reachable: if a *later* resolution fails inside
    # igraph, the earlier columns are already on the object. That is Seurat's
    # behaviour too, and is left alone.
    if algorithm == 3:
        raise NotImplementedError(
            "algorithm=3 (SLM) is not implemented. Use 1 or 2 (Louvain) or "
            "4 (Leiden)."
        )
    if algorithm not in (1, 2, 4):
        raise ValueError(
            f"Unknown algorithm {algorithm!r}. Use 1 or 2 (Louvain) or 4 (Leiden)."
        )

    # `np.number` is here for the numpy scalars a caller gets out of an array;
    # np.float64 subclasses float but np.float32 does not, and iterating one
    # raises rather than falling through to the sequence branch.
    if isinstance(resolution, (int, float, np.number)):
        resolutions = [float(resolution)]
    else:
        resolutions = [float(r) for r in resolution]
    if not resolutions:
        raise ValueError("`resolution` is empty; give at least one value.")

    if cluster_name is None:
        names = [f"{graph_name}_res.{_res_label(r)}" for r in resolutions]
    else:
        names = [cluster_name] if isinstance(cluster_name, str) else list(cluster_name)
        if len(names) != len(resolutions):
            raise ValueError(
                f"`cluster_name` has {len(names)} name(s) for "
                f"{len(resolutions)} resolution(s); give one per resolution."
            )

    cluster_series = None
    for res, name in zip(resolutions, names):
        if algorithm == 4:
            labels = _leiden_clustering(mat, res, random_seed, n_iterations)
        else:
            # python-igraph's community_multilevel is the multilevel Louvain
            # algorithm (closest to Seurat's algorithm 1/2).
            labels = _louvain_clustering(mat, res, random_seed)

        str_labels = _group_singletons(
            np.asarray([str(c) for c in labels]), mat, group_singletons
        )
        present = sorted(set(str_labels),
                         key=lambda s: (not s.isdigit(), s.isdigit() and int(s), s))
        cluster_series = pd.Categorical(str_labels, categories=present)
        seurat.meta_data[name] = cluster_series

    # The last resolution given, not the largest — Seurat takes the last column
    # of its results frame, so `resolution = c(1.2, 0.8, 0.4)` leaves the object
    # sitting on 0.4.
    seurat.meta_data["seurat_clusters"] = cluster_series
    seurat.idents = cluster_series