diff --git a/scanpy/neighbors/__init__.py b/scanpy/neighbors/__init__.py
index 55d0845aaa..133ab7786e 100644
--- a/scanpy/neighbors/__init__.py
+++ b/scanpy/neighbors/__init__.py
@@ -1,14 +1,13 @@
-from typing import Union, Optional, Any, Mapping, Callable
+from typing import Union, Optional, Any, Mapping, Callable, NamedTuple, Generator, Tuple
 
 import numpy as np
 import scipy
 from anndata import AnnData
 from numpy.random import RandomState
-from scipy.sparse import issparse, coo_matrix
+from scipy.sparse import issparse, coo_matrix, csr_matrix
 from sklearn.metrics import pairwise_distances
 from sklearn.utils import check_random_state
 
-from .._settings import settings
 from .. import logging as logg
 from .. import utils
 from ..utils import doc_params
@@ -19,6 +18,9 @@
 N_PCS = 50  # default number of PCs
 
 
+Metric = Callable[[np.ndarray, np.ndarray], float]
+
+
 @doc_params(n_pcs=doc_n_pcs, use_rep=doc_use_rep)
 def neighbors(
     adata: AnnData,
@@ -28,7 +30,7 @@ def neighbors(
     knn: bool = True,
     random_state: Optional[Union[int, RandomState]] = 0,
     method: str = 'umap',
-    metric: Union[str, Callable[[np.ndarray, np.ndarray], float]] = 'euclidean',
+    metric: Union[str, Metric] = 'euclidean',
     metric_kwds: Mapping[str, Any] = {},
     copy: bool = False
 ) -> Optional[AnnData]:
@@ -91,7 +93,8 @@ def neighbors(
     neighbors.compute_neighbors(
         n_neighbors=n_neighbors, knn=knn, n_pcs=n_pcs, use_rep=use_rep,
         method=method, metric=metric, metric_kwds=metric_kwds,
-        random_state=random_state)
+        random_state=random_state,
+    )
     adata.uns['neighbors'] = {}
     adata.uns['neighbors']['params'] = {'n_neighbors': n_neighbors, 'method': method}
     adata.uns['neighbors']['params']['metric'] = metric
@@ -113,12 +116,19 @@ def neighbors(
     return adata if copy else None
 
 
-def _rp_forest_generate(rp_forest_dict):
-    from collections import namedtuple
+class FlatTree(NamedTuple):
+    hyperplanes: None
+    offsets: None
+    children: None
+    indices: None
+
+
+RPForestDict = Mapping[str, Mapping[str, np.ndarray]]
 
-    props = ['hyperplanes', 'offsets', 'children', 'indices']
-    FlatTree = namedtuple('FlatTree', props)
-    num_trees = len(rp_forest_dict[props[0]]['start'])-1
+
+def _rp_forest_generate(rp_forest_dict: RPForestDict) -> Generator[FlatTree, None, None]:
+    props = FlatTree._fields[0]
+    num_trees = len(rp_forest_dict[props]['start'])-1
 
     for i in range(num_trees):
         tree = []
@@ -174,9 +184,14 @@ def neighbors_update(adata, adata_new, k=10, queue_size=5, random_state=0):
 
 
 def compute_neighbors_umap(
-        X, n_neighbors, random_state=None,
-        metric='euclidean', metric_kwds={}, angular=False,
-        verbose=False):
+    X: Union[np.ndarray, csr_matrix],
+    n_neighbors: int,
+    random_state: Optional[Union[int, RandomState]] = None,
+    metric: Union[str, Metric] = 'euclidean',
+    metric_kwds: Mapping[str, Any] = {},
+    angular: bool = False,
+    verbose: bool = False,
+):
     """This is from umap.fuzzy_simplicial_set [McInnes18]_.
 
     Given a set of data X, a neighborhood size, and a measure of distance
@@ -190,14 +205,14 @@ def compute_neighbors_umap(
     ----------
     X: array of shape (n_samples, n_features)
         The data to be modelled as a fuzzy simplicial set.
-    n_neighbors: int
+    n_neighbors
         The number of neighbors to use to approximate geodesic distance.
         Larger numbers induce more global estimates of the manifold that can
         miss finer detail, while smaller values will focus on fine manifold
         structure to the detriment of the larger picture.
-    random_state: numpy RandomState or equivalent
+    random_state
         A state capable being used as a numpy random state.
-    metric: string or function (optional, default 'euclidean')
+    metric
         The metric to use to compute distances in high dimensional space.
         If a string is passed it must match a valid predefined metric. If
         a general metric is required a function that takes two 1d arrays and
@@ -229,14 +244,14 @@ def compute_neighbors_umap(
         can have arguments passed via the metric_kwds dictionary. At this
         time care must be taken and dictionary elements must be ordered
         appropriately; this will hopefully be fixed in the future.
-    metric_kwds: dict (optional, default {})
+    metric_kwds
         Arguments to pass on to the metric, such as the ``p`` value for
         Minkowski distance.
-    angular: bool (optional, default False)
+    angular
         Whether to use angular/cosine distance for the random projection
         forest for seeding NN-descent to determine approximate nearest
         neighbors.
-    verbose: bool (optional, default False)
+    verbose
         Whether to report information on the current progress of the algorithm.
 
     Returns
@@ -247,9 +262,11 @@ def compute_neighbors_umap(
 
     random_state = check_random_state(random_state)
 
-    knn_indices, knn_dists, forest = nearest_neighbors(X, n_neighbors, random_state=random_state,
-                                                       metric=metric, metric_kwds=metric_kwds,
-                                                       angular=angular, verbose=verbose)
+    knn_indices, knn_dists, forest = nearest_neighbors(
+        X, n_neighbors, random_state=random_state,
+        metric=metric, metric_kwds=metric_kwds,
+        angular=angular, verbose=verbose,
+    )
 
     return knn_indices, knn_dists, forest
 
@@ -278,9 +295,11 @@ def get_sparse_matrix_from_indices_distances_umap(knn_indices, knn_dists, n_obs,
     return result.tocsr()
 
 
-def compute_connectivities_umap(knn_indices, knn_dists,
-        n_obs, n_neighbors, set_op_mix_ratio=1.0,
-        local_connectivity=1.0):
+def compute_connectivities_umap(
+    knn_indices, knn_dists,
+    n_obs, n_neighbors, set_op_mix_ratio=1.0,
+    local_connectivity=1.0,
+):
     """This is from umap.fuzzy_simplicial_set [McInnes18]_.
 
     Given a set of data X, a neighborhood size, and a measure of distance
@@ -305,15 +324,16 @@ def compute_connectivities_umap(knn_indices, knn_dists,
 def get_sparse_matrix_from_indices_distances_numpy(indices, distances, n_obs, n_neighbors):
     n_nonzero = n_obs * n_neighbors
     indptr = np.arange(0, n_nonzero + 1, n_neighbors)
-    D = scipy.sparse.csr_matrix((distances.copy().ravel(),  # copy the data, otherwise strange behavior here
-                                indices.copy().ravel(),
-                                indptr),
-                                shape=(n_obs, n_obs))
+    D = csr_matrix((
+        distances.copy().ravel(),  # copy the data, otherwise strange behavior here
+        indices.copy().ravel(),
+        indptr,
+    ), shape=(n_obs, n_obs))
     D.eliminate_zeros()
     return D
 
 
-def get_indices_distances_from_sparse_matrix(D, n_neighbors):
+def get_indices_distances_from_sparse_matrix(D, n_neighbors: int):
     indices = np.zeros((D.shape[0], n_neighbors), dtype=int)
     distances = np.zeros((D.shape[0], n_neighbors), dtype=D.dtype)
     n_neighbors_m1 = n_neighbors - 1
@@ -335,7 +355,7 @@ def get_indices_distances_from_sparse_matrix(D, n_neighbors):
     return indices, distances
 
 
-def get_indices_distances_from_dense_matrix(D, n_neighbors):
+def get_indices_distances_from_dense_matrix(D, n_neighbors: int):
     sample_range = np.arange(D.shape[0])[:, None]
     indices = np.argpartition(D, n_neighbors-1, axis=1)[:, :n_neighbors]
     indices = indices[sample_range, np.argsort(D[sample_range, indices])]
@@ -343,7 +363,7 @@ def get_indices_distances_from_dense_matrix(D, n_neighbors):
     return indices, distances
 
 
-def _backwards_compat_get_full_X_diffmap(adata):
+def _backwards_compat_get_full_X_diffmap(adata: AnnData) -> np.ndarray:
     if 'X_diffmap0' in adata.obs:
         return np.c_[adata.obs['X_diffmap0'].values[:, None],
                      adata.obsm['X_diffmap']]
@@ -351,7 +371,7 @@ def _backwards_compat_get_full_X_diffmap(adata):
         return adata.obsm['X_diffmap']
 
 
-def _backwards_compat_get_full_eval(adata):
+def _backwards_compat_get_full_eval(adata: AnnData):
     if 'X_diffmap0' in adata.obs:
         return np.r_[1, adata.uns['diffmap_evals']]
     else:
@@ -384,7 +404,15 @@ def _make_forest_dict(forest):
 class OnFlySymMatrix:
     """Emulate a matrix where elements are calculated on the fly.
     """
-    def __init__(self, get_row, shape, DC_start=0, DC_end=-1, rows=None, restrict_array=None):
+    def __init__(
+        self,
+        get_row: Callable[[Any], np.ndarray],
+        shape: Tuple[int, int],
+        DC_start: int = 0,
+        DC_end: int = -1,
+        rows: Optional[Mapping[Any, np.ndarray]] = None,
+        restrict_array: Optional[np.ndarray] = None,
+    ):
         self.get_row = get_row
         self.shape = shape
         self.DC_start = DC_start
@@ -393,7 +421,7 @@ def __init__(self, get_row, shape, DC_start=0, DC_end=-1, rows=None, restrict_ar
         self.restrict_array = restrict_array  # restrict the array to a subset
 
     def __getitem__(self, index):
-        if isinstance(index, int) or isinstance(index, np.integer):
+        if isinstance(index, (int, np.integer)):
             if self.restrict_array is None:
                 glob_index = index
             else:
@@ -420,9 +448,11 @@ def restrict(self, index_array):
         """Generate a view restricted to a subset of indices.
         """
         new_shape = index_array.shape[0], index_array.shape[0]
-        return OnFlySymMatrix(self.get_row, new_shape, DC_start=self.DC_start,
-                              DC_end=self.DC_end,
-                              rows=self.rows, restrict_array=index_array)
+        return OnFlySymMatrix(
+            self.get_row, new_shape, DC_start=self.DC_start,
+            DC_end=self.DC_end,
+            rows=self.rows, restrict_array=index_array,
+        )
 
 
 class Neighbors:
@@ -444,11 +474,12 @@ def __init__(self, adata: AnnData, n_dcs: Optional[int] = None):
         self._init_iroot()
         # use the graph in adata
         info_str = ''
-        self.knn = None
-        self._distances = None
-        self._connectivities = None
-        self._number_connected_components = None
-        self._rp_forest = None
+        self.knn: Optional[bool] = None
+        self._distances: Union[np.ndarray, csr_matrix, None] = None
+        self._connectivities: Union[np.ndarray, csr_matrix, None] = None
+        self._transitions_sym: Union[np.ndarray, csr_matrix, None] = None
+        self._number_connected_components: Optional[int] = None
+        self._rp_forest: Optional[RPForestDict] = None
         if 'neighbors' in adata.uns:
             if 'distances' in adata.uns['neighbors']:
                 self.knn = issparse(adata.uns['neighbors']['distances'])
@@ -461,11 +492,13 @@ def __init__(self, adata: AnnData, n_dcs: Optional[int] = None):
             if 'params' in adata.uns['neighbors']:
                 self.n_neighbors = adata.uns['neighbors']['params']['n_neighbors']
             else:
+                def count_nonzero(a: Union[np.ndarray, csr_matrix]) -> int:
+                    return a.count_nonzero() if issparse(a) else np.count_nonzero(a)
                 # estimating n_neighbors
                 if self._connectivities is None:
-                    self.n_neighbors = int(self._distances.count_nonzero() / self._distances.shape[0])
+                    self.n_neighbors = int(count_nonzero(self._distances) / self._distances.shape[0])
                 else:
-                    self.n_neighbors = int(self._connectivities.count_nonzero() / self._connectivities.shape[0] / 2)
+                    self.n_neighbors = int(count_nonzero(self._connectivities) / self._connectivities.shape[0] / 2)
             info_str += '`.distances` `.connectivities` '
             self._number_connected_components = 1
             if issparse(self._connectivities):
@@ -493,23 +526,23 @@ def __init__(self, adata: AnnData, n_dcs: Optional[int] = None):
             logg.msg('    initialized {}'.format(info_str), v=4)
 
     @property
-    def rp_forest(self):
+    def rp_forest(self) -> Optional[RPForestDict]:
         return self._rp_forest
 
     @property
-    def distances(self):
+    def distances(self) -> Union[np.ndarray, csr_matrix, None]:
         """Distances between data points (sparse matrix).
         """
         return self._distances
 
     @property
-    def connectivities(self):
+    def connectivities(self) -> Union[np.ndarray, csr_matrix, None]:
         """Connectivities between data points (sparse matrix).
         """
         return self._connectivities
 
     @property
-    def transitions(self):
+    def transitions(self) -> Union[np.ndarray, csr_matrix]:
         """Transition matrix (sparse matrix).
 
         Is conjugate to the symmetrized transition matrix via::
@@ -527,10 +560,10 @@ def transitions(self):
             Zinv = self.Z.power(-1)
         else:
             Zinv = np.diag(1./np.diag(self.Z))
-        return self.Z.dot(self.transitions_sym).dot(Zinv)
+        return self.Z @ self.transitions_sym @ Zinv
 
     @property
-    def transitions_sym(self):
+    def transitions_sym(self) -> Union[np.ndarray, csr_matrix, None]:
         """Symmetrized transition matrix (sparse matrix).
 
         Is conjugate to the transition matrix via::
@@ -554,12 +587,6 @@ def eigen_basis(self):
         """
         return self._eigen_basis
 
-    @property
-    def laplacian(self):
-        """Graph laplacian (sparse matrix).
-        """
-        return self._laplacian
-
     @property
     def distances_dpt(self):
         """DPT distances (on-fly matrix).
@@ -648,7 +675,11 @@ def compute_neighbors(
             # we need self._distances also for method == 'gauss' if we didn't
             # use dense distances
             self._distances, self._connectivities = compute_connectivities_umap(
-                knn_indices, knn_distances, self._adata.shape[0], self.n_neighbors)
+                knn_indices,
+                knn_distances,
+                self._adata.shape[0],
+                self.n_neighbors,
+            )
         # overwrite the umap connectivities if method is 'gauss'
         # self._distances is unaffected by this
         if method == 'gauss':
@@ -749,7 +780,7 @@ def compute_transitions(self, density_normalize=True):
                 Q = np.diag(1.0/q)
             else:
                 Q = scipy.sparse.spdiags(1.0/q, 0, W.shape[0], W.shape[0])
-            K = Q.dot(W).dot(Q)
+            K = Q @ W @ Q
         else:
             K = W
 
@@ -759,7 +790,7 @@ def compute_transitions(self, density_normalize=True):
             self.Z = np.diag(1.0/z)
         else:
             self.Z = scipy.sparse.spdiags(1.0/z, 0, K.shape[0], K.shape[0])
-        self._transitions_sym = self.Z.dot(K).dot(self.Z)
+        self._transitions_sym = self.Z @ K @ self.Z
         logg.msg('computed transitions', v=4, time=True)
 
     def compute_eigen(self, n_comps=15, sym=None, sort='decrease'):
@@ -837,77 +868,32 @@ def _init_iroot(self):
             self._set_iroot_via_xroot(xroot)
 
     def _get_dpt_row(self, i):
-        use_mask = False
+        mask = None
         if self._number_connected_components > 1:
-            use_mask = True
             label = self._connected_components[1][i]
             mask = self._connected_components[1] == label
-        row = sum([(self.eigen_values[l]/(1-self.eigen_values[l])
-                     * (self.eigen_basis[i, l] - self.eigen_basis[:, l]))**2
-                   # account for float32 precision
-                    for l in range(0, self.eigen_values.size) if self.eigen_values[l] < 0.9994])
+        row = sum(
+            (
+                self.eigen_values[l] / (1-self.eigen_values[l])
+                * (self.eigen_basis[i, l] - self.eigen_basis[:, l])
+            )**2
+            # account for float32 precision
+            for l in range(0, self.eigen_values.size)
+            if self.eigen_values[l] < 0.9994
+        )
         # thanks to Marius Lange for pointing Alex to this:
         # we will likely remove the contributions from the stationary state below when making
         # backwards compat breaking changes, they originate from an early implementation in 2015
         # they never seem to have deteriorated results, but also other distance measures (see e.g.
         # PAGA paper) don't have it, which makes sense
-        row += sum([(self.eigen_basis[i, l] - self.eigen_basis[:, l])**2
-                    for l in range(0, self.eigen_values.size) if self.eigen_values[l] >= 0.9994])
-        if not use_mask:
-            return np.sqrt(row)
-        else:
+        row += sum(
+            (self.eigen_basis[i, l] - self.eigen_basis[:, l])**2
+            for l in range(0, self.eigen_values.size)
+            if self.eigen_values[l] >= 0.9994
+        )
+        if mask is not None:
             row[~mask] = np.inf
-            return np.sqrt(row)
-
-    def _compute_Lp_matrix(self):
-        """See Fouss et al. (2006) and von Luxburg et al. (2007).
-
-        See Proposition 6 in von Luxburg (2007) and the inline equations
-        right in the text above.
-        """
-        self.Lp = sum([1/self.eigen_values[i]
-                      * np.outer(self.eigen_basis[:, i], self.eigen_basis[:, i])
-                      for i in range(1, self.eigen_values.size)])
-
-    def _compute_C_matrix(self):
-        """See Fouss et al. (2006) and von Luxburg et al. (2007).
-
-        This is the commute-time matrix. It's a squared-euclidian distance
-        matrix in :math:`\\mathbb{R}^n`.
-        """
-        self.C = np.repeat(np.diag(self.Lp)[:, np.newaxis],
-                           self.Lp.shape[0], axis=1)
-        self.C += np.repeat(np.diag(self.Lp)[np.newaxis, :],
-                            self.Lp.shape[0], axis=0)
-        self.C -= 2*self.Lp
-        # the following is much slower
-        # self.C = np.zeros(self.Lp.shape)
-        # for i in range(self.Lp.shape[0]):
-        #     for j in range(self.Lp.shape[1]):
-        #         self.C[i, j] = self.Lp[i, i] + self.Lp[j, j] - 2*self.Lp[i, j]
-        volG = np.sum(self.z)
-        self.C *= volG
-        settings.mt(0, 'computed commute distance matrix')
-        self.distances_dpt = self.C
-
-    def _compute_MFP_matrix(self):
-        """See Fouss et al. (2006).
-
-        This is the mean-first passage time matrix. It's not a distance.
-
-        Mfp[i, k] := m(k|i) in the notation of Fouss et al. (2006). This
-        corresponds to the standard notation for transition matrices (left index
-        initial state, right index final state, i.e. a right-stochastic
-        matrix, with each row summing to one).
-        """
-        self.MFP = np.zeros(self.Lp.shape)
-        for i in range(self.Lp.shape[0]):
-            for k in range(self.Lp.shape[1]):
-                for j in range(self.Lp.shape[1]):
-                    self.MFP[i, k] += (self.Lp[i, j] - self.Lp[i, k]
-                                       - self.Lp[k, j] + self.Lp[k, k]) * self.z[j]
-        settings.mt(0, 'computed mean first passage time matrix')
-        self.distances_dpt = self.MFP
+        return np.sqrt(row)
 
     def _set_pseudotime(self):
         """Return pseudotime with respect to root point.
@@ -936,7 +922,7 @@ def _set_iroot_via_xroot(self, xroot):
         iroot = 0
         for i in range(self._adata.shape[0]):
             diff = self._adata.X[i, :] - xroot
-            dsq = diff.dot(diff)
+            dsq = diff @ diff
             if dsq < dsqroot:
                 dsqroot = dsq
                 iroot = i
@@ -945,18 +931,3 @@ def _set_iroot_via_xroot(self, xroot):
         if self.iroot is not None and iroot != self.iroot:
             logg.warn('Changing index of iroot from {} to {}.'.format(self.iroot, iroot))
         self.iroot = iroot
-
-    def _test_embed(self):
-        """
-        Checks and tests for embed.
-        """
-        # pl.semilogy(w,'x',label=r'$ \widetilde K$')
-        # pl.show()
-        if _settings_verbosity_greater_or_equal_than(3):
-            # output of spectrum of K for comparison
-            w, v = np.linalg.eigh(self.K)
-            logg.msg('spectrum of K (kernel)')
-        if _settings_verbosity_greater_or_equal_than(4):
-            # direct computation of spectrum of T
-            w, vl, vr = scipy.linalg.eig(self.T, left=True)
-            logg.msg('spectrum of transition matrix (should be same as of Ktilde)')
diff --git a/scanpy/tests/test_binary.py b/scanpy/tests/test_binary.py
index 0b6af3e74a..f773be585d 100644
--- a/scanpy/tests/test_binary.py
+++ b/scanpy/tests/test_binary.py
@@ -40,7 +40,7 @@ def test_help_output(set_path: type(None), capsys: CaptureFixture):
     with pytest.raises(SystemExit, match='^0$'):
         main(['-h'])
     captured = capsys.readouterr()
-    assert re.search(r'^positional arguments:\n\s+\{settings,testbin[^}]*\}$', captured.out, re.MULTILINE)
+    assert re.search(r'^positional arguments:\n\s+\{settings,[\w,-]*testbin[\w,-]*\}$', captured.out, re.MULTILINE)
 
 
 def test_external(set_path: type(None)):
diff --git a/scanpy/tests/test_neighbors.py b/scanpy/tests/test_neighbors.py
index b88044f8ca..e37c371335 100644
--- a/scanpy/tests/test_neighbors.py
+++ b/scanpy/tests/test_neighbors.py
@@ -1,6 +1,7 @@
 import numpy as np
 import pytest
 from anndata import AnnData
+from scipy.sparse import csr_matrix
 
 from scanpy import Neighbors
 
@@ -13,13 +14,15 @@
     [0.0, 2.0, 0.0, 4.123105525970459],
     [2.0, 0.0, 0.0, 5.0],
     [0.0, 6.324555397033691, 0.0, 5.385164737701416],
-    [4.123105525970459, 5.0, 0.0, 0.0]]
+    [4.123105525970459, 5.0, 0.0, 0.0],
+]
 
 distances_euclidean_all = [
     [0.0, 2.0, 7.211102485656738, 4.123105525970459],
     [2.0, 0.0, 6.324555397033691, 5.0],
     [7.211102485656738, 6.324555397033691, 0.0, 5.385164737701416],
-    [4.123105525970459, 5.0, 5.385164737701416, 0.0]]
+    [4.123105525970459, 5.0, 5.385164737701416, 0.0],
+]
 
 
 # umap "kernel" - only knn results
@@ -27,19 +30,22 @@
     [0.0, 1.0, 0.0, 1.0],
     [1.0, 0.0, 0.5849691143165735, 0.8277419907567016],
     [0.0, 0.5849691143165735, 0.0, 1.0],
-    [1.0, 0.8277419907567016, 1.0, 0.0]]
+    [1.0, 0.8277419907567016, 1.0, 0.0],
+]
 
 transitions_sym_umap = [
     [0.0, 0.4818987107873648, 0.0, 0.3951883393150153],
     [0.48189871078736474, 0.0, 0.3594582764005241, 0.24216345431293487],
     [0.0, 0.3594582764005241, 0.0, 0.5039226836320637],
-    [0.39518833931501524, 0.24216345431293487, 0.5039226836320637, 0.0]]
+    [0.39518833931501524, 0.24216345431293487, 0.5039226836320637, 0.0],
+]
 
 transitions_umap = [
     [0.0, 0.5395987596963403, 0.0, 0.4604012403036599],
     [0.430368608684738, 0.0, 0.3176747629691457, 0.2519566283461165],
     [0.0, 0.40673754271561435, 0.0, 0.5932624572843856],
-    [0.33921243006981133, 0.23275092618009624, 0.42803664375009237, 0.0]]
+    [0.33921243006981133, 0.23275092618009624, 0.42803664375009237, 0.0],
+]
 
 
 # gauss kernel [diffmap, dpt] - knn and dense results
@@ -53,31 +59,36 @@
     [1.0, 0.676927387714386, 0.024883469566702843, 0.1962655782699585],
     [0.676927387714386, 1.0, 0.08414449542760849, 0.1353352814912796],
     [0.024883469566702843, 0.08414449542760849, 1.0, 0.16558068990707397],
-    [0.1962655782699585, 0.1353352814912796, 0.16558068990707397, 1.0]]
+    [0.1962655782699585, 0.1353352814912796, 0.16558068990707397, 1.0],
+]
 
 transitions_sym_gauss_knn = [
     [0.0, 0.5146393179893494, 0.0, 0.36445462703704834],
     [0.5146393179893494, 0.0, 0.3581143319606781, 0.2239987552165985],
     [0.0, 0.3581143319606781, 0.0, 0.5245543718338013],
-    [0.36445462703704834, 0.2239987552165985, 0.5245543718338013, 0.0]]
+    [0.36445462703704834, 0.2239987552165985, 0.5245543718338013, 0.0],
+]
 
 transitions_sym_gauss_noknn = [
     [0.5093212127685547, 0.34393802285194397, 0.016115963459014893, 0.11607448011636734],
     [0.34393805265426636, 0.506855845451355, 0.054364752024412155, 0.07984541356563568],
     [0.016115965321660042, 0.054364752024412155, 0.8235670328140259, 0.12452481687068939],
-    [0.11607448011636734, 0.07984541356563568, 0.1245248094201088, 0.6867417693138123]]
+    [0.11607448011636734, 0.07984541356563568, 0.1245248094201088, 0.6867417693138123],
+]
 
 transitions_gauss_knn = [
     [0.0, 0.5824036598205566, 0.0, 0.4175964295864105],
     [0.4547595679759979, 0.0, 0.3184431493282318, 0.22679725289344788],
     [0.0, 0.4027276933193207, 0.0, 0.5972723364830017],
-    [0.3180755078792572, 0.22123482823371887, 0.46068981289863586, 0.0]]
+    [0.3180755078792572, 0.22123482823371887, 0.46068981289863586, 0.0],
+]
 
 transitions_gauss_noknn = [
     [0.5093212127685547, 0.3450769782066345, 0.01887294091284275, 0.12672874331474304],
     [0.34280285239219666, 0.506855845451355, 0.06345486640930176, 0.08688655495643616],
     [0.01376173086464405, 0.04657683148980141, 0.8235670328140259, 0.11609435081481934],
-    [0.10631592571735382, 0.07337487488985062, 0.13356748223304749, 0.6867417693138123]]
+    [0.10631592571735382, 0.07337487488985062, 0.13356748223304749, 0.6867417693138123],
+]
 
 
 def get_neighbors() -> Neighbors:
@@ -130,3 +141,13 @@ def test_metrics_argument():
     no_knn_manhattan.compute_neighbors(method="gauss", knn=False,
         n_neighbors=n_neighbors, metric="manhattan")
     assert not np.allclose(no_knn_euclidean.distances, no_knn_manhattan.distances)
+
+
+@pytest.mark.parametrize('conv', [csr_matrix.toarray, csr_matrix])
+def test_restore_n_neighbors(neigh, conv):
+    neigh.compute_neighbors(method='gauss', n_neighbors=n_neighbors)
+
+    ad = AnnData(np.array(X))
+    ad.uns['neighbors'] = dict(connectivities=conv(neigh.connectivities))
+    neigh_restored = Neighbors(ad)
+    assert neigh_restored.n_neighbors == 1