Preprocessing

Counts to something you can do statistics on. Two routes, both of Seurat's: log-normalize → find variable features → scale, or sctransform in one call.

find_variable_features has two selectors and they honour different arguments — "vst" and "disp" take nfeatures, "mvp" takes the mean and dispersion cutoffs. That is Seurat's behaviour, not a quirk of the port; the docstring says which is which.

Verified against Seurat in PBMC 3k and, for the regularized-NB route, per fitted gene in SCTransform.

Log-normalize workflow

normalize_data

normalize_data(seurat, normalization_method: str = 'LogNormalize', scale_factor: float = 10000.0, assay: Optional[str] = None, margin: int = 1) -> None

Log-normalize counts.

Mirrors R's NormalizeData(pbmc, normalization.method = "LogNormalize", scale.factor = 10000). Modifies the assay in-place by adding / updating the 'data' layer.

Parameters:

• normalization_method (str, default: 'LogNormalize' ) –

  'LogNormalize', 'CLR', or 'RC'

• margin (int, default: 1 ) –

  for CLR only, and matching Seurat's flag exactly — normalize each feature across cells (1, Seurat's default) or each cell across its features (2). ADT/CITE-seq panels typically use margin=2.

Source code in truecell/preprocessing.py

def normalize_data(
    seurat,
    normalization_method: str = "LogNormalize",
    scale_factor: float = 10000.0,
    assay: Optional[str] = None,
    margin: int = 1,
) -> None:
    """Log-normalize counts.

    Mirrors R's NormalizeData(pbmc, normalization.method = "LogNormalize",
    scale.factor = 10000). Modifies the assay in-place by adding / updating
    the 'data' layer.

    Parameters
    ----------
    normalization_method : 'LogNormalize', 'CLR', or 'RC'
    margin               : for CLR only, and matching Seurat's flag exactly —
                           normalize each feature across cells (1, Seurat's
                           default) or each cell across its features (2).
                           ADT/CITE-seq panels typically use margin=2.
    """
    assay_name = assay or seurat.active_assay
    assay_obj = seurat.assays[assay_name]
    counts = _get_layer(assay_obj, "counts")

    if normalization_method == "LogNormalize":
        normed = _log_normalize(counts, scale_factor)
    elif normalization_method == "CLR":
        normed = _clr_normalize(counts, margin=margin)
    elif normalization_method == "RC":
        normed = _rc_normalize(counts, scale_factor)
    else:
        raise ValueError(f"Unknown normalization_method: {normalization_method!r}")

    _set_layer(assay_obj, "data", normed)
    log_truecell_command(
        seurat, "NormalizeData", assay=assay or seurat.active_assay,
        params={"normalization_method": normalization_method,
                "scale_factor": scale_factor, "margin": margin},
    )

find_variable_features

find_variable_features(seurat, selection_method: str = 'vst', nfeatures: int = 2000, assay: Optional[str] = None, layer: Optional[str] = None, mean_cutoff: tuple = (0.1, 8), dispersion_cutoff: tuple = (1, float('inf')), num_bin: int = 20, binning_method: str = 'equal_width') -> None

Select highly variable features.

Mirrors R's FindVariableFeatures(pbmc, selection.method = "vst", nfeatures = 2000). Modifies the assay in-place by setting var_features (Assay) or highly_variable in meta_data (Assay5).

On an Assay5 the per-feature statistics land in assay.meta_data under the names HVFInfo() uses, so a column reads the same in either language:

• selection_method="vst" — mean, variance, variance.expected, variance.standardized
• "mvp" / "mean.var.plot" / "dispersion" / "disp" — mvp.mean, mvp.dispersion, mvp.dispersion.scaled

plus highly_variable, a boolean flag that has no Seurat counterpart of its own (HVFInfo(status = TRUE) spells it variable); it is the fallback Assay5.variable_features reads when the ordered list is empty.

The two dispersion spellings are not synonyms, however much they share. Seurat routes them to different selectors — MVP for "mvp" / "mean.var.plot", DISP for "dispersion" / "disp" — and only the second honours nfeatures. mean_cutoff and dispersion_cutoff apply to the first, and to nothing else; they were accepted and discarded before, so mean.var.plot returned a top-nfeatures list under a name that promises a cutoff.

Source code in truecell/preprocessing.py

def find_variable_features(
    seurat,
    selection_method: str = "vst",
    nfeatures: int = 2000,
    assay: Optional[str] = None,
    layer: Optional[str] = None,
    mean_cutoff: tuple = (0.1, 8),
    dispersion_cutoff: tuple = (1, float("inf")),
    num_bin: int = 20,
    binning_method: str = "equal_width",
) -> None:
    """Select highly variable features.

    Mirrors R's FindVariableFeatures(pbmc, selection.method = "vst",
    nfeatures = 2000). Modifies the assay in-place by setting
    var_features (Assay) or highly_variable in meta_data (Assay5).

    On an Assay5 the per-feature statistics land in ``assay.meta_data`` under
    the names `HVFInfo()` uses, so a column reads the same in either language:

    * ``selection_method="vst"`` — ``mean``, ``variance``,
      ``variance.expected``, ``variance.standardized``
    * ``"mvp"`` / ``"mean.var.plot"`` / ``"dispersion"`` / ``"disp"`` —
      ``mvp.mean``, ``mvp.dispersion``, ``mvp.dispersion.scaled``

    plus ``highly_variable``, a boolean flag that has no Seurat counterpart of
    its own (`HVFInfo(status = TRUE)` spells it ``variable``); it is the
    fallback `Assay5.variable_features` reads when the ordered list is empty.

    The two dispersion spellings are **not** synonyms, however much they share.
    Seurat routes them to different selectors — `MVP` for ``"mvp"`` /
    ``"mean.var.plot"``, `DISP` for ``"dispersion"`` / ``"disp"`` — and only
    the second honours ``nfeatures``. ``mean_cutoff`` and ``dispersion_cutoff``
    apply to the first, and to nothing else; they were accepted and discarded
    before, so ``mean.var.plot`` returned a top-``nfeatures`` list under a name
    that promises a cutoff.
    """
    assay_name = assay or seurat.active_assay
    assay_obj = seurat.assays[assay_name]

    if layer is not None:
        data = _get_layer(assay_obj, layer)
    else:
        from .assay5 import Assay5
        from ._sparse import is_matrix_empty
        if selection_method == "vst":
            # Seurat fits the vst mean-variance LOESS on RAW COUNTS, regardless of
            # whether NormalizeData() has run. Using normalized data here would
            # change which features are selected.
            if isinstance(assay_obj, Assay5):
                data = assay_obj.layers.get("counts")
                if data is None:
                    data = assay_obj.layers.get("data")
            else:
                data = assay_obj.counts
        else:
            # dispersion / mean.var.plot methods operate on log-normalized data
            if isinstance(assay_obj, Assay5):
                data = assay_obj.layers.get("data")
                if data is None:
                    data = assay_obj.layers.get("counts")
            else:
                data = assay_obj.data if not is_matrix_empty(assay_obj.data) else assay_obj.counts

    # `stats` becomes the per-feature columns written to meta_data below. The
    # keys are Seurat's, taken from `HVFInfo()` rather than from the raw slot:
    # SeuratObject stores these prefixed by method and layer
    # (`vf_vst_counts_mean`), and strips the prefix on the way out. `HVFInfo()`
    # is the name a Seurat user actually types, so it is the one to match.
    if selection_method == "vst":
        hvg_indices, means, variances, expected_var, var_std = _vst_hvg(data, nfeatures)
        stats = {
            "mean": means,
            "variance": variances,
            "variance.expected": expected_var,
            "variance.standardized": var_std,
        }
    elif selection_method in ("dispersion", "disp", "mvp", "mean.var.plot"):
        hvg_indices, means, dispersion, dispersion_scaled = _dispersion_hvg(
            data, nfeatures, mean_cutoff, dispersion_cutoff,
            num_bin=num_bin, binning_method=binning_method,
            select="disp" if selection_method in ("dispersion", "disp") else "mvp",
        )
        # A different method produces different quantities, and Seurat names
        # them so: `HVFInfo(method = "mvp")` returns mvp.mean / mvp.dispersion /
        # mvp.dispersion.scaled, never `variance.standardized`. Writing scaled
        # dispersions into a column called `variance.standardized` — which is
        # what sharing one set of column names across both methods amounts to —
        # is a mislabelling that no downstream reader can detect.
        stats = {
            "mvp.mean": means,
            "mvp.dispersion": dispersion,
            "mvp.dispersion.scaled": dispersion_scaled,
        }
    else:
        raise ValueError(f"Unknown selection_method: {selection_method!r}")

    from .assay5 import Assay5
    if isinstance(assay_obj, Assay5):
        feature_names = assay_obj._all_feature_names
    else:
        feature_names = assay_obj._feature_names

    hvg_names = [feature_names[i] for i in hvg_indices]

    # Store results
    if isinstance(assay_obj, Assay5):
        # Store HVF info in assay meta_data
        hvf_df = pd.DataFrame(
            {**stats, "highly_variable": np.zeros(len(feature_names), dtype=bool)},
            index=feature_names,
        )
        hvf_df.loc[hvg_names, "highly_variable"] = True
        # Retire the other method's columns. SeuratObject can keep both — it
        # namespaces them by method and layer (`vf_vst_counts_mean`) and
        # `HVFInfo(method = )` picks — but truecell's meta_data *is* the
        # user-facing table, with one flat name per statistic. Leaving the
        # previous method's columns in place next to a `variable_features` list
        # they no longer describe is how `variable_feature_plot` came to draw
        # standardized variances over an mvp selection.
        for col in _HVF_COLUMNS - set(hvf_df.columns):
            assay_obj.meta_data.drop(columns=col, inplace=True, errors="ignore")
        for col in hvf_df.columns:
            assay_obj.meta_data[col] = hvf_df[col]
        assay_obj.variable_features = hvg_names
    else:
        assay_obj.var_features = hvg_names
    log_truecell_command(
        seurat, "FindVariableFeatures", assay=assay or seurat.active_assay,
        params={"selection_method": selection_method, "nfeatures": nfeatures},
    )

scale_data

scale_data(seurat, features: Optional[list[str]] = None, vars_to_regress: Optional[list[str]] = None, assay: Optional[str] = None, do_scale: bool = True, do_center: bool = True, scale_max: float = 10.0, layer: str = 'data') -> None

Scale and optionally center expression data.

Mirrors R's ScaleData(). Stores result in the 'scale.data' layer (Assay5) or scale_data slot (Assay v3).

Parameters:

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

  genes to scale (defaults to variable features)

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

  metadata columns to regress out before scaling

• do_scale (bool, default: True ) –

  standardize variance to 1

• do_center (bool, default: True ) –

  subtract mean

• scale_max (float, default: 10.0 ) –

  clip scaled values at this magnitude

Source code in truecell/preprocessing.py

def scale_data(
    seurat,
    features: Optional[list[str]] = None,
    vars_to_regress: Optional[list[str]] = None,
    assay: Optional[str] = None,
    do_scale: bool = True,
    do_center: bool = True,
    scale_max: float = 10.0,
    layer: str = "data",
) -> None:
    """Scale and optionally center expression data.

    Mirrors R's ScaleData(). Stores result in the 'scale.data' layer
    (Assay5) or scale_data slot (Assay v3).

    Parameters
    ----------
    features      : genes to scale (defaults to variable features)
    vars_to_regress : metadata columns to regress out before scaling
    do_scale      : standardize variance to 1
    do_center     : subtract mean
    scale_max     : clip scaled values at this magnitude
    """
    assay_name = assay or seurat.active_assay
    assay_obj = seurat.assays[assay_name]

    from .assay5 import Assay5

    if isinstance(assay_obj, Assay5):
        all_features = assay_obj._all_feature_names
        feat_idx_map = {f: i for i, f in enumerate(all_features)}
    else:
        all_features = assay_obj._feature_names
        feat_idx_map = {f: i for i, f in enumerate(all_features)}

    # Select features to scale
    if features is None:
        if isinstance(assay_obj, Assay5):
            features = assay_obj.variable_features or all_features
        else:
            features = assay_obj.var_features or all_features

    feat_idx = [feat_idx_map[f] for f in features if f in feat_idx_map]
    features_present = [all_features[i] for i in feat_idx]

    # Get log-normalized data for the selected features (features × cells)
    data = _get_layer(assay_obj, layer)
    if sp.issparse(data) or is_lazy(data):
        # Subset first, densify second -- on a lazy layer the reverse would
        # read the whole store off disk to keep a few thousand rows of it.
        sub = data[feat_idx, :].toarray().astype(float)
    else:
        sub = np.asarray(data)[feat_idx, :].astype(float)

    # Regress out covariates
    if vars_to_regress:
        sub = _regress_out(sub, seurat.meta_data, vars_to_regress, features_present)

    # Center
    if do_center:
        gene_means = sub.mean(axis=1, keepdims=True)
        sub = sub - gene_means

    # Scale (sample SD, ddof=1, matching Seurat's ScaleData)
    if do_scale:
        gene_stds = sub.std(axis=1, ddof=1, keepdims=True)
        gene_stds[gene_stds == 0] = 1.0
        sub = sub / gene_stds

    # Clip
    sub = np.clip(sub, -scale_max, scale_max)

    # Store scaled data
    scaled_sparse = sp.csc_matrix(sub)
    if isinstance(assay_obj, Assay5):
        # Add/replace scale.data layer; update internal cell/feature maps
        assay_obj.set_layer_data("scale.data", scaled_sparse, feature_names=features_present)
        # Store which features are scaled (needed for PCA)
        assay_obj._scaled_features = features_present
    else:
        # Through `_set_layer` so the row labels are written with the matrix.
        # Assigning `assay_obj.scale_data` directly leaves `_scaled_features`
        # describing whatever was there before.
        assay_obj._set_layer("scale_data", sub, features_present)
    log_truecell_command(
        seurat, "ScaleData", assay=assay or seurat.active_assay,
        params={"do_center": do_center, "do_scale": do_scale,
                "n_features": len(features_present)},
    )

percentage_feature_set

percentage_feature_set(seurat, pattern: str, col_name: Optional[str] = None, assay: Optional[str] = None, layer: str = 'counts') -> None

Add a metadata column with % of counts matching a gene name pattern.

Mirrors R's PercentageFeatureSet(pbmc, pattern = "^MT-"). Modifies seurat.meta_data in-place.

Source code in truecell/preprocessing.py

def percentage_feature_set(
    seurat,
    pattern: str,
    col_name: Optional[str] = None,
    assay: Optional[str] = None,
    layer: str = "counts",
) -> None:
    """Add a metadata column with % of counts matching a gene name pattern.

    Mirrors R's PercentageFeatureSet(pbmc, pattern = "^MT-").
    Modifies seurat.meta_data in-place.
    """
    assay_obj = _get_assay(seurat, assay)
    mat = _get_layer(assay_obj, layer)

    from .assay5 import Assay5
    if isinstance(assay_obj, Assay5):
        feature_names = assay_obj._all_feature_names
    else:
        feature_names = assay_obj._feature_names

    # Find matching features
    rx = re.compile(pattern)
    match_mask = np.array([bool(rx.search(f)) for f in feature_names])

    if not match_mask.any():
        pct = np.zeros(mat.shape[1])
    else:
        # A lazy layer belongs on the sparse branch: indexing it returns a
        # sparse block, so it satisfies this branch's contract exactly, while
        # the dense one would `np.asarray` the whole store just to sum it.
        if sp.issparse(mat) or is_lazy(mat):
            total = np.array(mat.sum(axis=0)).flatten()
            matching = np.array(mat[match_mask, :].sum(axis=0)).flatten()
        else:
            total = mat.sum(axis=0)
            matching = mat[match_mask, :].sum(axis=0)
        total[total == 0] = 1
        pct = (matching / total) * 100.0

    if col_name is None:
        col_name = "percent.mt" if re.search(r"mt|mito", pattern, re.I) else "percent_feature"

    cells = seurat.cell_names()
    if isinstance(assay_obj, Assay5):
        cell_list = assay_obj._all_cell_names
    else:
        cell_list = assay_obj._cell_names

    seurat.meta_data[col_name] = pd.Series(pct, index=cell_list).reindex(cells).values

Regularized negative binomial

sctransform

sctransform(seurat, assay: Optional[str] = None, new_assay_name: str = 'SCT', n_cells: int = 5000, n_genes: int = 2000, n_features: int = 3000, min_cells: int = 5, vars_to_regress: Optional[list[str]] = None, clip_range: Optional[tuple] = None, gene_chunk: int = 500, seed: int = 42, set_default: bool = True, vst_flavor: str = 'v2', bw_adjust: float = 3.0, verbose: bool = False) -> 'object'

Run SCTransform and attach a normalized assay.

Mirrors R's SCTransform(object). Fits the regularized NB model on the active assay's counts and stores the result as new_assay_name ("SCT").

Parameters:

• n_cells (int, default: 5000 ) –

  cells subsampled for parameter estimation (Seurat default 5000).

• n_genes (int, default: 2000 ) –

  genes used for step-1 estimation, sampled to spread evenly over expression (Seurat default 2000). None uses every gene.

• n_features (int, default: 3000 ) –

  number of variable features by residual variance (default 3000).

• min_cells (int, default: 5 ) –

  drop genes detected in fewer than this many cells (default 5).

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

  metadata columns (e.g. 'percent.mt') regressed out of the Pearson residuals, mirroring SCTransform's vars.to.regress.

• clip_range (Optional[tuple], default: None ) –

  residual clip for scale.data; default (-√(N/30), √(N/30)). Note this is not the clip used when ranking variable features — see below.

• vst_flavor (str, default: 'v2' ) –

  "v2" (default, as Seurat 5) or "v1". See the module docstring.

• bw_adjust (float, default: 3.0 ) –

  multiplier on the Sheather-Jones smoothing bandwidth (R's 3).

• set_default (bool, default: True ) –

  make the new assay the active assay.

Returns:

• ``seurat``, with the new assay added. –

Source code in truecell/sctransform.py

def sctransform(
    seurat,
    assay: Optional[str] = None,
    new_assay_name: str = "SCT",
    n_cells: int = 5000,
    n_genes: int = 2000,
    n_features: int = 3000,
    min_cells: int = 5,
    vars_to_regress: Optional[list[str]] = None,
    clip_range: Optional[tuple] = None,
    gene_chunk: int = 500,
    seed: int = 42,
    set_default: bool = True,
    vst_flavor: str = "v2",
    bw_adjust: float = 3.0,
    verbose: bool = False,
) -> "object":
    """Run SCTransform and attach a normalized assay.

    Mirrors R's ``SCTransform(object)``. Fits the regularized NB model on the
    active assay's counts and stores the result as ``new_assay_name`` ("SCT").

    Parameters
    ----------
    n_cells     : cells subsampled for parameter estimation (Seurat default 5000).
    n_genes     : genes used for step-1 estimation, sampled to spread evenly over
                  expression (Seurat default 2000). ``None`` uses every gene.
    n_features  : number of variable features by residual variance (default 3000).
    min_cells   : drop genes detected in fewer than this many cells (default 5).
    vars_to_regress : metadata columns (e.g. 'percent.mt') regressed out of the
                  Pearson residuals, mirroring SCTransform's vars.to.regress.
    clip_range  : residual clip for ``scale.data``; default (-√(N/30), √(N/30)).
                  Note this is *not* the clip used when ranking variable features
                  — see below.
    vst_flavor  : "v2" (default, as Seurat 5) or "v1". See the module docstring.
    bw_adjust   : multiplier on the Sheather-Jones smoothing bandwidth (R's 3).
    set_default : make the new assay the active assay.

    Returns
    -------
    ``seurat``, with the new assay added.
    """
    if vst_flavor not in ("v1", "v2"):
        raise ValueError(f"vst_flavor must be 'v1' or 'v2', got {vst_flavor!r}")

    src = seurat.assays[assay or seurat.active_assay]
    # A layer may be dense, sparse, or an on-disk LazyMatrix — `sp.csc_matrix`
    # is what accepts all three, so keep `issparse` as the discriminator below.
    counts: Any
    if isinstance(src, Assay5):
        counts = src.layers.get("counts")
        all_genes = src._all_feature_names
    else:
        counts = src.counts
        all_genes = src._feature_names
    if counts is None:
        raise ValueError("SCTransform requires a counts layer.")

    counts = counts.tocsc() if sp.issparse(counts) else sp.csc_matrix(counts)
    cell_names = seurat.cell_names()
    G_all, N = counts.shape

    # Drop genes detected in too few cells.
    nnz_per_gene = np.diff(counts.tocsr().indptr)
    keep = np.where(nnz_per_gene >= min_cells)[0]
    counts = counts[keep, :]
    genes = [all_genes[i] for i in keep]
    G = len(genes)
    counts_csr = counts.tocsr()
    if verbose:
        print(f"  SCTransform ({vst_flavor}): {G}/{G_all} genes kept "
              f"(>= {min_cells} cells), {N} cells")

    total_umi = np.asarray(counts.sum(axis=0)).ravel().astype(float)
    total_umi[total_umi == 0] = 1.0
    log10_umi = np.log10(total_umi)
    log10_gmean = np.log10(_row_gmean(counts_csr, eps=1.0))

    rng = np.random.default_rng(seed)

    # ---- step 1: which cells and genes estimate the model ----
    if N > n_cells:
        cells_step1 = np.sort(rng.choice(N, n_cells, replace=False))
        det = np.diff(counts_csr[:, cells_step1].tocsr().indptr)
        genes_step1 = np.where(det >= min_cells)[0]
    else:
        cells_step1 = np.arange(N)
        genes_step1 = np.arange(G)

    gene_amean = np.asarray(counts_csr.mean(axis=1)).ravel()
    gene_var = _row_var(counts_csr)
    overdispersion = gene_var - gene_amean
    poisson_genes = np.zeros(G, dtype=bool)

    if vst_flavor == "v2":
        # Genes whose variance does not exceed their mean carry no NB signal;
        # regularizing them drags the trend, so R models them as pure Poisson.
        poisson_genes = (overdispersion <= 0) | (gene_amean < 0.001)
        genes_step1 = genes_step1[overdispersion[genes_step1] > 0]
        if verbose:
            print(f"  {poisson_genes.sum()} poisson genes excluded from regularization")

    if n_genes is not None and n_genes < len(genes_step1):
        # Sample step-1 genes inversely to their density in log-gmean, so the
        # smoother sees the sparse tails, not just the crowded middle.
        x = log10_gmean[genes_step1]
        bw = _bw_nrd(x)
        dens = np.exp(-0.5 * ((x[:, None] - x[None, :]) / bw) ** 2).sum(axis=1)
        dens /= (len(x) * bw * _SQRT2PI)
        prob = 1.0 / (dens + np.finfo(float).eps)
        genes_step1 = rng.choice(genes_step1, size=n_genes, replace=False,
                                 p=prob / prob.sum())
        genes_step1 = np.sort(genes_step1)

    log10_gmean_step1 = log10_gmean[genes_step1]
    Y1 = counts_csr[genes_step1][:, cells_step1].toarray().astype(float)
    log10_umi_step1 = log10_umi[cells_step1]

    if verbose:
        print(f"  fitting NB model on {len(cells_step1)} cells x "
              f"{len(genes_step1)} genes ...")

    # ---- step 1: fit ----
    if vst_flavor == "v2":
        b0_s1, theta_s1 = _fit_nb_offset(Y1, log10_umi_step1, gene_chunk)
        b1_s1 = np.full(len(genes_step1), np.log(10.0))
    else:
        b0_s1, b1_s1, theta_s1 = _fit_poisson(Y1, log10_umi_step1, gene_chunk)
    del Y1

    # ---- step 2: regularize ----
    model_pars = np.column_stack([_dispersion_par(log10_gmean_step1, theta_s1),
                                  b0_s1, b1_s1])
    fit = _reg_model_pars(model_pars, log10_gmean_step1, log10_gmean, bw_adjust, verbose)

    theta_r = _theta_from_dispersion_par(log10_gmean, fit[:, 0])
    b0r, b1r = fit[:, 1], fit[:, 2]

    if vst_flavor == "v2":
        mean_cell_sum = total_umi.mean()
        with np.errstate(divide="ignore"):
            b0r = np.where(poisson_genes,
                           np.log(np.maximum(gene_amean, 1e-300)) - np.log(mean_cell_sum),
                           b0r)
        theta_r = np.where(poisson_genes, np.inf, theta_r)
        b1r = np.full(G, np.log(10.0))  # fix_slope

    bad = ~np.isfinite(theta_r) & ~poisson_genes
    theta_r = np.where(bad, 1e6, theta_r)
    theta_r = np.where(poisson_genes, np.inf, np.clip(theta_r, 1e-2, 1e6))

    # ---- step 3: residuals ----
    # Two different clips, as in R. Residual *variance* — which ranks the
    # variable features — comes from residuals clipped at sqrt(N); only the
    # stored scale.data is clipped to the much tighter sqrt(N/30). Using the
    # tight clip for both crushes exactly the marker genes that define rare
    # subsets, so the ranking loses them.
    res_clip = np.sqrt(N)
    if clip_range is None:
        c = np.sqrt(N / 30.0)
        clip_lo, clip_hi = -c, c
    else:
        clip_lo, clip_hi = clip_range

    min_var = -np.inf
    if vst_flavor == "v2":
        # R: (median non-zero UMI / 5)^2 — a floor on model variance that stops
        # near-zero fitted means from manufacturing enormous residuals.
        min_var = (np.median(counts_csr.data) / 5.0) ** 2
        if verbose:
            print(f"  min_variance = {min_var:.6g}")

    median_log10_umi = float(np.median(log10_umi))
    res_var = np.zeros(G)
    res_mean = np.zeros(G)
    corrected_blocks = []

    for start in range(0, G, gene_chunk):
        end = min(start + gene_chunk, G)
        y = counts_csr[start:end].toarray().astype(float)
        mu = np.exp(np.clip(b0r[start:end, None] + b1r[start:end, None] * log10_umi[None, :],
                            -30, 30))
        var = mu + mu * mu / theta_r[start:end, None]
        var = np.maximum(var, min_var)
        z = (y - mu) / np.sqrt(var)
        z_clipped = np.clip(z, -res_clip, res_clip)
        res_var[start:end] = z_clipped.var(axis=1, ddof=1)
        res_mean[start:end] = z_clipped.mean(axis=1)

        # Corrected counts: the residual re-expressed at the median depth. R
        # uses the unclipped residual and no variance floor here.
        mu_med = np.exp(b0r[start:end, None] + b1r[start:end, None] * median_log10_umi)
        var_med = mu_med + mu_med * mu_med / theta_r[start:end, None]
        z_raw = (y - mu) / np.sqrt(mu + mu * mu / theta_r[start:end, None])
        corr = np.clip(np.round(mu_med + z_raw * np.sqrt(var_med)), 0.0, None)
        corrected_blocks.append(sp.csr_matrix(corr))

    corrected = sp.vstack(corrected_blocks, format="csc")
    del corrected_blocks

    # ---- variable features by residual variance ----
    n_feat = min(n_features, G)
    top = np.argsort(res_var)[::-1][:n_feat]
    var_features = [genes[i] for i in top]

    # ---- scale.data: residuals for the variable features only ----
    top_sorted = np.sort(top)
    y = counts_csr[top_sorted].toarray().astype(float)
    mu = np.exp(np.clip(b0r[top_sorted, None] + b1r[top_sorted, None] * log10_umi[None, :],
                        -30, 30))
    var = np.maximum(mu + mu * mu / theta_r[top_sorted, None], min_var)
    resid = np.clip((y - mu) / np.sqrt(var), clip_lo, clip_hi)
    scale_feats = [genes[i] for i in top_sorted]

    if vars_to_regress:
        resid = _regress_out(resid, seurat.meta_data, list(vars_to_regress), scale_feats)

    sct = Assay5(
        layers={"counts": corrected, "data": _log1p_sparse(corrected)},
        feature_names=list(genes),
        cell_names=list(cell_names),
        key=f"{new_assay_name.lower()}_",
    )
    sct.set_layer_data("scale.data", sp.csc_matrix(resid), feature_names=scale_feats)
    sct._scaled_features = scale_feats
    sct.variable_features = var_features
    # The fitted model, per gene. Column names are Seurat's, from
    # `obj[["SCT"]]@SCTModel.list[[1]]@feature.attributes` — including the
    # awkward `(Intercept)`, because this table exists to be read alongside
    # Seurat's and renaming the columns would defeat that. `b0r`/`b1r` are the
    # *regularized* coefficients, which is what R stores here too (the
    # unregularized step-1 fits are separate `step1_*` columns there).
    detection_rate = np.diff(counts_csr.indptr) / float(N)
    sct.meta_data["residual_variance"] = pd.Series(res_var, index=genes)
    sct.meta_data["residual_mean"] = pd.Series(res_mean, index=genes)
    sct.meta_data["theta"] = pd.Series(theta_r, index=genes)
    sct.meta_data["gmean"] = pd.Series(10.0 ** log10_gmean, index=genes)
    sct.meta_data["detection_rate"] = pd.Series(detection_rate, index=genes)
    sct.meta_data["(Intercept)"] = pd.Series(b0r, index=genes)
    sct.meta_data["log_umi"] = pd.Series(b1r, index=genes)

    seurat.assays[new_assay_name] = sct
    if set_default:
        seurat.active_assay = new_assay_name
    if verbose:
        print(f"  SCT assay '{new_assay_name}' added: {G} genes, "
              f"{n_feat} variable features")
    return seurat