Dimensional reduction

Linear first, then the embeddings you look at. Every one of these writes a DimReduc into obj.reductions under a key, and records the features it actually used — not the features it was asked for, which are not always the same set.

jack_straw is the permutation test for how many PCs to keep. It is worth reading its docstring before trusting the number: R's JackRandom seeds each replicate from its loop index and is therefore deterministic, while this one seeds from its seed argument and moves. Across 60 seeds on PBMC 3k it keeps 12–15 PCs, mode 13, which is R's answer. That spread is asserted as a band, not described in prose — see Fidelity.

Linear

run_pca

run_pca(seurat, n_pcs: int = 50, features: Optional[list[str]] = None, assay: Optional[str] = None, reduction_name: str = 'pca', reduction_key: str = 'PC_', seed: int = 42, layer: str = 'scale.data') -> None

Compute PCA on scaled data.

Mirrors R's RunPCA(pbmc, features = VariableFeatures(object = pbmc)). Stores a DimReduc in seurat.reductions[reduction_name].

Parameters:

• n_pcs (int, default: 50 ) –

  number of principal components

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

  genes to use (defaults to variable features)

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

  assay name (defaults to active assay)

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

  key for storage in seurat.reductions

• reduction_key (str, default: 'PC_' ) –

  prefix for dimension names (e.g. 'PC_')

• seed (int, default: 42 ) –

  random seed for reproducibility

• layer (str, default: 'scale.data' ) –

  which layer to take data from

Source code in truecell/reduction.py

def run_pca(
    seurat,
    n_pcs: int = 50,
    features: Optional[list[str]] = None,
    assay: Optional[str] = None,
    reduction_name: str = "pca",
    reduction_key: str = "PC_",
    seed: int = 42,
    layer: str = "scale.data",
) -> None:
    """Compute PCA on scaled data.

    Mirrors R's RunPCA(pbmc, features = VariableFeatures(object = pbmc)).
    Stores a DimReduc in seurat.reductions[reduction_name].

    Parameters
    ----------
    n_pcs          : number of principal components
    features       : genes to use (defaults to variable features)
    assay          : assay name (defaults to active assay)
    reduction_name : key for storage in seurat.reductions
    reduction_key  : prefix for dimension names (e.g. 'PC_')
    seed           : random seed for reproducibility
    layer          : which layer to take data from
    """
    assay_name = assay or seurat.active_assay
    assay_obj = seurat.assays[assay_name]

    features = _default_features(assay_obj, features)

    # Get scale.data for the selected features. `used` is what the layer
    # actually carried; it is what labels the loadings below, because a feature
    # the layer does not have has no row to be labelled.
    scaled, used = _prep_dr(assay_obj, features, layer)
    # scaled shape: (n_features_selected × n_cells)
    # irlba is handed t(scale.data), so transpose to (n_cells × n_features)
    data_t = scaled.T  # (cells × features)

    n_pcs = min(n_pcs, min(data_t.shape) - 1)

    np.random.seed(seed)
    u, d, vt = _truncated_svd(data_t, n_pcs, seed)
    embeddings = u * d                    # RunPCA.default: u %*% diag(d)
    loadings = vt.T                       # (features × n_pcs), pca.results$v
    # sdev <- pca.results$d / sqrt(max(1, ncol(object) - 1)), where `object` is
    # features × cells — so the denominator counts cells.
    stdev = d / np.sqrt(max(1, data_t.shape[0] - 1))

    # Cell names
    cells = seurat.cell_names()

    dr = DimReduc(
        cell_embeddings=embeddings,
        feature_loadings=loadings,
        assay_used=assay_name,
        stdev=stdev,
        key=reduction_key,
        cell_names=cells,
        feature_names=used,
    )

    seurat.reductions[reduction_name] = dr
    log_truecell_command(
        seurat, "RunPCA", assay=assay_name,
        params={"n_pcs": n_pcs, "reduction_name": reduction_name, "seed": seed},
    )

run_spca

run_spca(seurat, graph: str, npcs: int = 50, features: Optional[list[str]] = None, assay: Optional[str] = None, reduction_name: str = 'spca', reduction_key: str = 'SPC_', seed: int = 42, layer: str = 'scale.data') -> None

Supervised PCA — the gene axes that best explain a cell-cell graph.

Mirrors R's RunSPCA(obj, assay = "SCT", graph = "wsnn"). Ordinary PCA picks the directions of greatest variance and knows nothing about which cells you consider neighbours. sPCA is handed a graph you already trust — typically the WNN graph from find_multi_modal_neighbors, which knows about protein as well as RNA — and finds the directions in gene space that best reproduce it. Where PCA maximises vᵀXᵀXv, sPCA maximises vᵀXᵀGXv: the same problem with the identity swapped for the graph. Set G = I and you get PCA back exactly.

The point is the loadings. Because sPCA is still a linear map from genes to components, a query dataset can be pushed into a reference's graph-defined space with a single matrix multiply, which is what makes it the reduction Azimuth maps onto.

Parameters:

• graph (str) –

  key in seurat.graphs — a cell × cell graph (e.g. "wsnn")

• npcs (int, default: 50 ) –

  number of components

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

  genes to use (defaults to variable features)

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

  assay name (defaults to active assay)

• reduction_name (str, default: 'spca' ) –

  key for storage in seurat.reductions

• reduction_key (str, default: 'SPC_' ) –

  prefix for dimension names (e.g. 'SPC_')

• seed (int, default: 42 ) –

  random seed (the eigensolver starts from a random vector)

• layer (str, default: 'scale.data' ) –

  which layer to take data from

Notes

Seurat runs irlba on XᵀGX, which is an SVD and so ranks components by |λ|; we take the largest eigenvalues themselves, since vᵀXᵀGXv is what is being maximised and a graph can push some eigenvalues negative. With non-negative edge weights the leading eigenvalues are positive and the two orderings agree, so this only ever differs in the tail.

Source code in truecell/reduction.py

def run_spca(
    seurat,
    graph: str,
    npcs: int = 50,
    features: Optional[list[str]] = None,
    assay: Optional[str] = None,
    reduction_name: str = "spca",
    reduction_key: str = "SPC_",
    seed: int = 42,
    layer: str = "scale.data",
) -> None:
    """Supervised PCA — the gene axes that best explain a cell-cell graph.

    Mirrors R's ``RunSPCA(obj, assay = "SCT", graph = "wsnn")``. Ordinary PCA
    picks the directions of greatest variance and knows nothing about which
    cells you consider neighbours. sPCA is handed a graph you already trust —
    typically the WNN graph from `find_multi_modal_neighbors`, which knows about
    protein as well as RNA — and finds the directions *in gene space* that best
    reproduce it. Where PCA maximises ``vᵀXᵀXv``, sPCA maximises ``vᵀXᵀGXv``:
    the same problem with the identity swapped for the graph. Set ``G = I`` and
    you get PCA back exactly.

    The point is the loadings. Because sPCA is still a linear map from genes to
    components, a query dataset can be pushed into a reference's graph-defined
    space with a single matrix multiply, which is what makes it the reduction
    Azimuth maps onto.

    Parameters
    ----------
    graph          : key in ``seurat.graphs`` — a cell × cell graph (e.g. "wsnn")
    npcs           : number of components
    features       : genes to use (defaults to variable features)
    assay          : assay name (defaults to active assay)
    reduction_name : key for storage in seurat.reductions
    reduction_key  : prefix for dimension names (e.g. 'SPC_')
    seed           : random seed (the eigensolver starts from a random vector)
    layer          : which layer to take data from

    Notes
    -----
    Seurat runs `irlba` on ``XᵀGX``, which is an SVD and so ranks components by
    ``|λ|``; we take the largest eigenvalues themselves, since ``vᵀXᵀGXv`` is
    what is being maximised and a graph can push some eigenvalues negative. With
    non-negative edge weights the leading eigenvalues are positive and the two
    orderings agree, so this only ever differs in the tail.
    """
    assay_name = assay or seurat.active_assay
    assay_obj = seurat.assays[assay_name]

    if graph not in seurat.graphs:
        raise KeyError(
            f"Graph '{graph}' not found. Run find_neighbors() or "
            f"find_multi_modal_neighbors() first. "
            f"Available graphs: {list(seurat.graphs)}"
        )

    features = _default_features(assay_obj, features)

    # `_prep_dr` densifies on the way out, so there is nothing left to convert
    # here — the `sp.issparse` branch that used to stand between these two lines
    # could not run.
    scaled, used = _prep_dr(assay_obj, features, layer)
    X = np.asarray(scaled.T, dtype=float)
    n_cells, n_features = X.shape

    G = seurat.graphs[graph].tocsr()
    if G.shape != (n_cells, n_cells):
        raise ValueError(
            f"Graph '{graph}' is {G.shape[0]}×{G.shape[1]} but the assay has "
            f"{n_cells} cells."
        )
    # A KNN graph need not be symmetric; the eigendecomposition needs it to be,
    # and (G + Gᵀ)/2 leaves the quadratic form vᵀXᵀGXv untouched anyway.
    G = (G + G.T) * 0.5

    npcs = min(npcs, n_features - 1)
    if npcs < 1:
        raise ValueError(
            f"Need at least 2 features to run sPCA; got {n_features}.")

    np.random.seed(seed)
    Z = X.T @ (G @ X)                                       # features × features
    Z = np.asarray(Z)
    Z = (Z + Z.T) * 0.5                                     # float asymmetry
    eigenvalues, loadings = _top_eigenvectors(Z, npcs, seed=seed)
    loadings = _flip_signs(loadings)                        # reproducible signs

    embeddings = X @ loadings                               # cells × npcs
    stdev = np.sqrt(np.var(embeddings, axis=0, ddof=1))

    seurat.reductions[reduction_name] = DimReduc(
        cell_embeddings=embeddings,
        feature_loadings=loadings,
        assay_used=assay_name,
        stdev=stdev,
        key=reduction_key,
        cell_names=seurat.cell_names(),
        feature_names=used,
        misc={"spca_graph": graph, "eigenvalues": eigenvalues},
    )

run_ica

run_ica(seurat, nics: int = 50, features: Optional[list[str]] = None, assay: Optional[str] = None, reduction_name: str = 'ica', reduction_key: str = 'ICA_', seed: int = 42, layer: str = 'scale.data', max_iter: int = 200) -> None

Independent Component Analysis on scaled data.

Mirrors R's RunICA(obj, nics = 50). Stores a DimReduc (embeddings + loadings) under reduction_name; find_neighbors / run_umap already accept reduction="ica".

Source code in truecell/reduction.py

def run_ica(
    seurat,
    nics: int = 50,
    features: Optional[list[str]] = None,
    assay: Optional[str] = None,
    reduction_name: str = "ica",
    reduction_key: str = "ICA_",
    seed: int = 42,
    layer: str = "scale.data",
    max_iter: int = 200,
) -> None:
    """Independent Component Analysis on scaled data.

    Mirrors R's ``RunICA(obj, nics = 50)``. Stores a DimReduc (embeddings +
    loadings) under ``reduction_name``; ``find_neighbors`` / ``run_umap``
    already accept ``reduction="ica"``.
    """
    from sklearn.decomposition import FastICA

    assay_name = assay or seurat.active_assay
    assay_obj = seurat.assays[assay_name]

    features = _default_features(assay_obj, features)

    scaled, used = _prep_dr(assay_obj, features, layer)  # (features × cells)
    data_t = scaled.T  # (cells × features) — already dense, see run_pca

    nics = min(nics, min(data_t.shape))

    ica = FastICA(n_components=nics, random_state=seed, max_iter=max_iter)
    embeddings = ica.fit_transform(data_t)  # (cells × nics)
    loadings = ica.components_.T  # (features × nics)

    cells = seurat.cell_names()
    seurat.reductions[reduction_name] = DimReduc(
        cell_embeddings=embeddings,
        feature_loadings=loadings,
        assay_used=assay_name,
        key=reduction_key,
        cell_names=cells,
        feature_names=used,
    )

glm_pca

glm_pca(seurat, n_components: int = 10, features: Optional[list[str]] = None, assay: Optional[str] = None, reduction_name: str = 'glmpca', reduction_key: str = 'GLMPC_', family: str = 'poisson', layer: str = 'counts', max_iter: int = 100, tol: float = 0.0001, penalty: float = 1.0, learning_rate: float = 0.1, theta: float = 100.0, optimize_theta: bool = True, seed: int = 42) -> None

Fit a Poisson GLM-PCA and store it as a DimReduc.

Mirrors R's RunGLMPCA(obj, L = 10). Takes raw counts, not normalised or scaled data — the whole point is to model the counts as counts. Stores factors as cell_embeddings and loadings as feature_loadings, so find_neighbors(obj, reduction="glmpca") and run_umap work downstream exactly as they do off PCA.

Parameters:

• n_components (int, default: 10 ) –

  rank of the fit (L) — the number of factors

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

  genes to use (defaults to variable features)

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

  assay name (defaults to active assay)

• family (str, default: 'poisson' ) –

  noise model — "poisson" or "nb" (negative binomial)

• layer (str, default: 'counts' ) –

  layer to read counts from

• max_iter (int, default: 100 ) –

  maximum Fisher scoring iterations

• tol (float, default: 0.0001 ) –

  stop when the relative change in deviance falls below this

• penalty (float, default: 1.0 ) –

  L2 ridge on U and V. U and V can trade scale freely (U·Vᵀ = (cU)·(V/c)ᵀ); the ridge is what pins that down.

• learning_rate (float, default: 0.1 ) –

  initial Fisher step size. Halved on any step that fails to lower the deviance, so this is an opening bid, not a commitment.

• theta (float, default: 100.0 ) –

  negative-binomial dispersion (Var = μ + μ²/θ). Ignored for Poisson. With optimize_theta it is only the starting value; otherwise it is held fixed at this number for the whole fit.

• optimize_theta (bool, default: True ) –

  re-estimate θ by maximum likelihood between factor updates (NB only). Turn off to fit at a dispersion you already trust — that also restores strict monotone deviance, since a moving θ re-scales the deviance under it.

• seed (int, default: 42 ) –

  only used if the counts have no structure at all — the fit is started deterministically from the data (see _init_factors)

Notes

Deviance falls monotonically by construction: a step that raises it (or overflows) is rejected outright and retried at half the step size. The full trace is kept in reduction.misc["deviance"] — if it is still dropping steeply at the end, raise max_iter. misc["converged"] says whether the fit stopped because it was done or because it ran out of iterations.

Fitting is dense in genes × cells. Pass a few thousand variable features rather than the whole transcriptome, as you would to run_pca.

Source code in truecell/glmpca.py

def glm_pca(
    seurat,
    n_components: int = 10,
    features: Optional[list[str]] = None,
    assay: Optional[str] = None,
    reduction_name: str = "glmpca",
    reduction_key: str = "GLMPC_",
    family: str = "poisson",
    layer: str = "counts",
    max_iter: int = 100,
    tol: float = 1e-4,
    penalty: float = 1.0,
    learning_rate: float = 0.1,
    theta: float = 100.0,
    optimize_theta: bool = True,
    seed: int = 42,
) -> None:
    """Fit a Poisson GLM-PCA and store it as a DimReduc.

    Mirrors R's ``RunGLMPCA(obj, L = 10)``. Takes **raw counts**, not normalised
    or scaled data — the whole point is to model the counts as counts. Stores
    factors as ``cell_embeddings`` and loadings as ``feature_loadings``, so
    ``find_neighbors(obj, reduction="glmpca")`` and ``run_umap`` work downstream
    exactly as they do off PCA.

    Parameters
    ----------
    n_components  : rank of the fit (L) — the number of factors
    features      : genes to use (defaults to variable features)
    assay         : assay name (defaults to active assay)
    family        : noise model — ``"poisson"`` or ``"nb"`` (negative binomial)
    layer         : layer to read counts from
    max_iter      : maximum Fisher scoring iterations
    tol           : stop when the relative change in deviance falls below this
    penalty       : L2 ridge on U and V. U and V can trade scale freely
                    (``U·Vᵀ = (cU)·(V/c)ᵀ``); the ridge is what pins that down.
    learning_rate : initial Fisher step size. Halved on any step that fails to
                    lower the deviance, so this is an opening bid, not a
                    commitment.
    theta         : negative-binomial dispersion (``Var = μ + μ²/θ``). Ignored for
                    Poisson. With ``optimize_theta`` it is only the starting value;
                    otherwise it is held fixed at this number for the whole fit.
    optimize_theta: re-estimate ``θ`` by maximum likelihood between factor updates
                    (NB only). Turn off to fit at a dispersion you already trust —
                    that also restores strict monotone deviance, since a moving
                    ``θ`` re-scales the deviance under it.
    seed          : only used if the counts have no structure at all — the fit is
                    started deterministically from the data (see `_init_factors`)

    Notes
    -----
    Deviance falls monotonically by construction: a step that raises it (or
    overflows) is rejected outright and retried at half the step size. The full
    trace is kept in ``reduction.misc["deviance"]`` — if it is still dropping
    steeply at the end, raise ``max_iter``. ``misc["converged"]`` says whether the
    fit stopped because it was done or because it ran out of iterations.

    Fitting is dense in genes × cells. Pass a few thousand variable features
    rather than the whole transcriptome, as you would to `run_pca`.
    """
    if family not in FAMILIES:
        raise NotImplementedError(
            f"family={family!r} is not implemented; use one of {list(FAMILIES)}.")

    from .markers import _get_expression_matrix
    from .reduction import _default_features

    assay_name = assay or seurat.active_assay
    assay_obj = seurat.assays[assay_name]

    features = _default_features(assay_obj, features)
    Y = _counts_for(assay_obj, features, layer, _get_expression_matrix)

    n_genes, n_cells = Y.shape
    n_components = min(n_components, min(n_genes, n_cells) - 1)
    if n_components < 1:
        raise ValueError(
            f"Need at least 2 genes and 2 cells to fit GLM-PCA; got "
            f"{n_genes} × {n_cells}.")

    # Library size as a fixed offset: depth is known, not something to rediscover.
    totals = Y.sum(axis=0)
    if np.any(totals <= 0):
        raise ValueError(
            "Some cells have zero total counts, so they have no library size to "
            "offset by. Filter them out first.")
    offsets = np.log(totals / totals.mean())

    if family == "poisson":
        intercept, U, V, deviance, converged = _fit_poisson(
            Y, n_components, offsets,
            max_iter=max_iter, tol=tol, penalty=penalty,
            learning_rate=learning_rate, seed=seed,
        )
        fitted_theta = np.inf
    else:                                                   # family == "nb"
        intercept, U, V, deviance, converged, fitted_theta = _fit_nb(
            Y, n_components, offsets,
            max_iter=max_iter, tol=tol, penalty=penalty,
            learning_rate=learning_rate, theta=theta,
            optimize_theta=optimize_theta, seed=seed,
        )
    loadings, factors = _orthogonalize(U, V)

    seurat.reductions[reduction_name] = DimReduc(
        cell_embeddings=factors,
        feature_loadings=loadings,
        assay_used=assay_name,
        stdev=np.sqrt(np.var(factors, axis=0, ddof=1)),
        key=reduction_key,
        cell_names=seurat.cell_names(),
        feature_names=list(features),
        misc={
            "glmpca_family": family,
            "deviance": deviance,
            "converged": converged,
            "intercept": intercept,
            "theta": fitted_theta,
        },
    )

Non-linear embeddings

run_umap

run_umap(seurat, dims: Optional[Union[list[int], range]] = None, reduction: str = 'pca', graph: Optional[str] = None, n_components: int = 2, n_neighbors: int = 30, min_dist: float = 0.3, metric: str = 'euclidean', reduction_name: str = 'umap', reduction_key: str = 'UMAP_', seed: int = 42, assay: Optional[str] = None) -> None

Compute a UMAP embedding.

Mirrors R's RunUMAP(pbmc, dims = 1:10) and RunUMAP(pbmc, graph = "wsnn"). Stores a DimReduc in seurat.reductions[reduction_name].

Two input modes (mutually exclusive):

• reduction (default): embed the cells from a low-dimensional reduction (PCA/Harmony/ICA…). The fitted umap-learn model is stashed in dr.misc["umap_model"] for later transform-only projection.
• graph: embed a precomputed neighbour graph directly (e.g. the WNN "wsnn" graph from find_multi_modal_neighbors), via UMAP's simplicial_set_embedding. reduction/dims are ignored.

Parameters:

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

  which dimensions of 'reduction' to use (0-indexed)

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

  source reduction ('pca' by default)

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

  name of a precomputed graph in seurat.graphs to embed (takes precedence over reduction when given)

• n_components (int, default: 2 ) –

  output dimensions (2 for visualization)

• n_neighbors (int, default: 30 ) –

  UMAP n_neighbors (Seurat default 30)

• min_dist (float, default: 0.3 ) –

  UMAP min_dist (Seurat default 0.3)

• metric (str, default: 'euclidean' ) –

  distance metric (reduction mode only)

• reduction_name (str, default: 'umap' ) –

  storage key in seurat.reductions

• seed (int, default: 42 ) –

  random seed

Source code in truecell/umap.py

def run_umap(
    seurat,
    dims: Optional[Union[list[int], range]] = None,
    reduction: str = "pca",
    graph: Optional[str] = None,
    n_components: int = 2,
    n_neighbors: int = 30,
    min_dist: float = 0.3,
    metric: str = "euclidean",
    reduction_name: str = "umap",
    reduction_key: str = "UMAP_",
    seed: int = 42,
    assay: Optional[str] = None,
) -> None:
    """Compute a UMAP embedding.

    Mirrors R's RunUMAP(pbmc, dims = 1:10) and RunUMAP(pbmc, graph = "wsnn").
    Stores a DimReduc in seurat.reductions[reduction_name].

    Two input modes (mutually exclusive):

    * ``reduction`` (default): embed the cells from a low-dimensional
      reduction (PCA/Harmony/ICA…). The fitted ``umap-learn`` model is stashed
      in ``dr.misc["umap_model"]`` for later transform-only projection.
    * ``graph``: embed a precomputed neighbour graph directly (e.g. the WNN
      ``"wsnn"`` graph from ``find_multi_modal_neighbors``), via UMAP's
      ``simplicial_set_embedding``. ``reduction``/``dims`` are ignored.

    Parameters
    ----------
    dims           : which dimensions of 'reduction' to use (0-indexed)
    reduction      : source reduction ('pca' by default)
    graph          : name of a precomputed graph in seurat.graphs to embed
                     (takes precedence over ``reduction`` when given)
    n_components   : output dimensions (2 for visualization)
    n_neighbors    : UMAP n_neighbors (Seurat default 30)
    min_dist       : UMAP min_dist (Seurat default 0.3)
    metric         : distance metric (reduction mode only)
    reduction_name : storage key in seurat.reductions
    seed           : random seed
    """
    assay_name = assay or seurat.active_assay
    cells = seurat.cell_names()
    dim_names = [f"{reduction_key}{i + 1}" for i in range(n_components)]

    if graph is not None:
        coords = _umap_from_graph(seurat, graph, n_components, min_dist, seed)
        seurat.reductions[reduction_name] = DimReduc(
            cell_embeddings=coords,
            assay_used=assay_name,
            key=reduction_key,
            cell_names=cells,
            feature_names=dim_names,
            misc={"umap_graph": graph},
        )
        return

    from umap import UMAP

    if reduction not in seurat.reductions:
        raise KeyError(f"Reduction '{reduction}' not found. Run run_pca() first.")

    embeddings = seurat.reductions[reduction].cell_embeddings  # (cells × n_dims)
    if dims is None:
        emb = embeddings
    else:
        emb = embeddings[:, list(dims)]

    reducer = UMAP(
        n_components=n_components,
        n_neighbors=n_neighbors,
        min_dist=min_dist,
        metric=metric,
        random_state=seed,
    )
    umap_coords = reducer.fit_transform(emb)  # (cells × n_components)

    seurat.reductions[reduction_name] = DimReduc(
        cell_embeddings=umap_coords,
        assay_used=assay_name,
        key=reduction_key,
        cell_names=cells,
        feature_names=dim_names,
        misc={"umap_model": reducer},
    )

run_tsne

run_tsne(seurat, dims: Optional[list[int]] = None, reduction: str = 'pca', n_components: int = 2, perplexity: float = 30.0, reduction_name: str = 'tsne', reduction_key: str = 'tSNE_', seed: int = 42, assay: Optional[str] = None) -> None

t-SNE embedding from an existing reduction.

Mirrors R's RunTSNE(obj, dims = 1:10). Stores a DimReduc under reduction_name.

Source code in truecell/reduction.py

def run_tsne(
    seurat,
    dims: Optional[list[int]] = None,
    reduction: str = "pca",
    n_components: int = 2,
    perplexity: float = 30.0,
    reduction_name: str = "tsne",
    reduction_key: str = "tSNE_",
    seed: int = 42,
    assay: Optional[str] = None,
) -> None:
    """t-SNE embedding from an existing reduction.

    Mirrors R's ``RunTSNE(obj, dims = 1:10)``. Stores a DimReduc under
    ``reduction_name``.
    """
    from sklearn.manifold import TSNE

    assay_name = assay or seurat.active_assay

    if reduction not in seurat.reductions:
        raise KeyError(f"Reduction '{reduction}' not found. Run run_pca() first.")
    emb = seurat.reductions[reduction].cell_embeddings
    if dims is not None:
        emb = emb[:, list(dims)]

    tsne = TSNE(
        n_components=n_components,
        perplexity=perplexity,
        random_state=seed,
        init="pca",
    )
    coords = tsne.fit_transform(emb)

    cells = seurat.cell_names()
    dim_names = [f"{reduction_key}{i + 1}" for i in range(n_components)]
    seurat.reductions[reduction_name] = DimReduc(
        cell_embeddings=coords,
        assay_used=assay_name,
        key=reduction_key,
        cell_names=cells,
        feature_names=dim_names,
    )

How many components to keep

jack_straw

jack_straw(seurat, reduction: str = 'pca', dims: int = 20, num_replicate: int = 100, prop_freq: float = 0.01, layer: str = 'scale.data', seed: int = 42) -> 'JackStrawData'

Permutation test for the significance of PCA dimensions.

Mirrors R's JackStraw(): a small fraction (prop_freq) of features is permuted across cells, the PCA is re-run on the permuted matrix, and the permuted features' loadings in that refit basis form the null distribution per PC. Each observed loading is then assigned an empirical p-value. Results are stored on seurat.reductions[reduction].jackstraw and returned.

The refit is the expensive part and it is not optional: an earlier version of this function built the null by projecting the permuted rows onto the fixed original embedding, which is far cheaper but produces a much tighter null — a fixed basis cannot rotate to absorb the scrambled signal, so the permuted loadings come out too small and ordinary noise features look extreme against them. On pbmc3k that inflated the count of "significant" features on the pure-noise PCs 14-20 from R's 0-5 to 109-203, and left score_jackstraw unable to reject any PC at all.

Cost scales as num_replicate full PCAs; ~1-2 minutes for the Seurat defaults on a 2000-feature, 2700-cell object. Lower num_replicate when iterating, but note it also sets the p-value resolution: the smallest non-zero empirical p is 1 / (num_replicate * n_permuted).

Parameters:

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

  reduction to test (default 'pca')

• dims (int, default: 20 ) –

  number of PCs to score

• num_replicate (int, default: 100 ) –

  permutation replicates (Seurat default 100)

• prop_freq (float, default: 0.01 ) –

  fraction of features permuted per replicate (default 0.01)

• layer (str, default: 'scale.data' ) –

  scaled layer feeding the reduction

• seed (int, default: 42 ) –

  RNG seed

Source code in truecell/jackstraw.py

def jack_straw(
    seurat,
    reduction: str = "pca",
    dims: int = 20,
    num_replicate: int = 100,
    prop_freq: float = 0.01,
    layer: str = "scale.data",
    seed: int = 42,
) -> "JackStrawData":
    """Permutation test for the significance of PCA dimensions.

    Mirrors R's ``JackStraw()``: a small fraction (``prop_freq``) of features is
    permuted across cells, the **PCA is re-run on the permuted matrix**, and the
    permuted features' loadings in that refit basis form the null distribution
    per PC. Each observed loading is then assigned an empirical p-value. Results
    are stored on ``seurat.reductions[reduction].jackstraw`` and returned.

    The refit is the expensive part and it is not optional: an earlier version of
    this function built the null by projecting the permuted rows onto the
    *fixed* original embedding, which is far cheaper but produces a much tighter
    null — a fixed basis cannot rotate to absorb the scrambled signal, so the
    permuted loadings come out too small and ordinary noise features look
    extreme against them. On pbmc3k that inflated the count of "significant"
    features on the pure-noise PCs 14-20 from R's 0-5 to 109-203, and left
    :func:`score_jackstraw` unable to reject any PC at all.

    Cost scales as ``num_replicate`` full PCAs; ~1-2 minutes for the Seurat
    defaults on a 2000-feature, 2700-cell object. Lower ``num_replicate`` when
    iterating, but note it also sets the p-value resolution: the smallest
    non-zero empirical p is ``1 / (num_replicate * n_permuted)``.

    Parameters
    ----------
    reduction      : reduction to test (default 'pca')
    dims           : number of PCs to score
    num_replicate  : permutation replicates (Seurat default 100)
    prop_freq      : fraction of features permuted per replicate (default 0.01)
    layer          : scaled layer feeding the reduction
    seed           : RNG seed
    """
    if reduction not in seurat.reductions:
        raise KeyError(f"Reduction '{reduction}' not found. Run run_pca() first.")
    dr = seurat.reductions[reduction]

    X, features = _scaled_matrix_for_reduction(seurat, dr, layer)
    n_features, n_cells = X.shape
    ndims = int(min(dims, dr.cell_embeddings.shape[1]))

    # The observed statistic is the reduction's own feature loadings, exactly as
    # R takes Loadings(object[[reduction]], projected = FALSE).
    loadings = np.asarray(dr.feature_loadings, dtype=float)
    if loadings.size == 0:
        raise ValueError(
            f"Reduction '{reduction}' has no feature loadings; JackStraw needs "
            "them as the observed statistic. Re-run run_pca().")
    if loadings.shape[0] != n_features:
        raise ValueError(
            f"Reduction '{reduction}' has {loadings.shape[0]} loadings for "
            f"{n_features} scaled features — the reduction and the layer "
            "disagree about the feature set.")
    obs_stat = np.abs(loadings[:, :ndims])                         # features × ndims

    # R: sample(rownames, size = nrow * prop.use), floored, with a hard floor of 3.
    n_perm = max(3, int(n_features * prop_freq))
    if n_perm > n_features:
        raise ValueError(
            f"prop_freq={prop_freq} selects {n_perm} of only {n_features} features")

    rng = np.random.default_rng(seed)
    null_chunks = []
    for _ in range(num_replicate):
        idx = rng.choice(n_features, size=n_perm, replace=False)
        # Scramble in place and restore afterwards: copying the whole matrix per
        # replicate would dominate the runtime on a real object.
        saved = X[idx, :].copy()
        for r in idx:
            X[r, :] = rng.permutation(X[r, :])
        refit = _refit_loadings(X, ndims, seed)
        null_chunks.append(np.abs(refit[idx, :]))                  # n_perm × ndims
        X[idx, :] = saved
    null_all = np.vstack(null_chunks)                              # (R·n_perm) × ndims

    # R's EmpiricalP: the fraction of null loadings STRICTLY greater than the
    # observed one. searchsorted(side='right') counts null <= obs.
    empirical = np.empty((n_features, ndims))
    n_null = null_all.shape[0]
    for j in range(ndims):
        col = np.sort(null_all[:, j])
        ranks = np.searchsorted(col, obs_stat[:, j], side="right")
        empirical[:, j] = (n_null - ranks) / n_null

    js = JackStrawData(
        empirical_p_values=empirical,
        fake_reduction_scores=null_all,
        overall_p_values=None,
        score=obs_stat,
        method="jackstraw",
    )
    dr.jackstraw = js
    return js

score_jackstraw

score_jackstraw(seurat, reduction: str = 'pca', dims: Optional[int] = None, score_thresh: float = 1e-05) -> ndarray

Aggregate per-feature JackStraw p-values into one p-value per PC.

Mirrors R's ScoreJackStraw(): count the features whose empirical p-value falls at or below score_thresh, and test that count against the number expected under a uniform null (floor(n_features * score_thresh)) with a two-proportion test. A small returned value marks a significant PC. A PC with no feature below the threshold scores exactly 1, as in R.

Not a distributional goodness-of-fit test. An earlier version used a one-sided KS test against Uniform(0, 1), which is enormously more sensitive: with thousands of features it returned p-values around 1e-112 or smaller for every PC on pbmc3k, including pure noise, so no PC ever failed and the function could not do the one job it exists for — telling you where to cut.

Source code in truecell/jackstraw.py

def score_jackstraw(
    seurat,
    reduction: str = "pca",
    dims: Optional[int] = None,
    score_thresh: float = 1e-5,
) -> np.ndarray:
    """Aggregate per-feature JackStraw p-values into one p-value per PC.

    Mirrors R's ``ScoreJackStraw()``: count the features whose empirical p-value
    falls at or below ``score_thresh``, and test that count against the number
    expected under a uniform null (``floor(n_features * score_thresh)``) with a
    two-proportion test. A *small* returned value marks a significant PC. A PC
    with no feature below the threshold scores exactly 1, as in R.

    Not a distributional goodness-of-fit test. An earlier version used a
    one-sided KS test against Uniform(0, 1), which is enormously more sensitive:
    with thousands of features it returned p-values around 1e-112 or smaller for
    *every* PC on pbmc3k, including pure noise, so no PC ever failed and the
    function could not do the one job it exists for — telling you where to cut.
    """
    if reduction not in seurat.reductions:
        raise KeyError(f"Reduction '{reduction}' not found.")
    dr = seurat.reductions[reduction]
    js = dr.jackstraw
    if js is None or js.is_empty():
        raise ValueError("Run jack_straw() before score_jackstraw().")

    emp = js.empirical_p_values
    ndims = emp.shape[1] if dims is None else int(min(dims, emp.shape[1]))
    n_features = emp.shape[0]
    expected = float(np.floor(n_features * score_thresh))

    overall = np.ones(ndims)
    for j in range(ndims):
        observed = int((np.clip(emp[:, j], 0.0, 1.0) <= score_thresh).sum())
        if observed == 0:
            overall[j] = 1.0                 # R's explicit guard
        else:
            overall[j] = _prop_test(observed, expected, n_features, n_features)

    js.overall_p_values = overall
    return overall