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}")

    umi_assay_name = assay or seurat.active_assay
    src = seurat.assays[umi_assay_name]
    # 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)

    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 this object's median depth.
    # `median(log10(umi))` rather than `log10(median(umi))` because that is what
    # sctransform's `correct` does when no `scale_factor` is passed — it takes
    # the median of the latent variable itself. `prep_sct_find_markers` is the
    # one caller that *does* pass a scale factor, and there R uses
    # `log10(median(umi))`; the two differ whenever the cell count is even.
    corrected = _correct_counts(counts_csr, b0r, b1r, theta_r, log10_umi,
                                median_log10_umi, gene_chunk)

    # ---- 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)

    # The fitted model, recorded whole. `meta_data` above is per *assay*, so a
    # merge of two SCTransformed objects collapses it to one table and the fact
    # that the two halves were corrected at different sequencing depths is
    # lost. This is Seurat's `SCTModel.list`: keyed by model, carrying the cells
    # it was fitted on and the median UMI it corrected to, which is exactly what
    # `prep_sct_find_markers` needs to put them back on a common scale.
    sct.misc["SCTModel.list"] = {
        "model1": {
            "feature_attributes": pd.DataFrame(
                {"theta": theta_r, "(Intercept)": b0r, "log_umi": b1r},
                index=list(genes),
            ),
            "cell_attributes": pd.DataFrame(
                {"umi": total_umi, "log_umi": log10_umi}, index=list(cell_names),
            ),
            # R's `median(cell.attributes[, "umi"])` — the median on the *count*
            # scale, which is what `PrepSCTFindMarkers` minimises over. Not the
            # same as `median_log10_umi` above for an even number of cells.
            "median_umi": float(np.median(total_umi)),
            "umi_assay": umi_assay_name,
        }
    }

    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

Run prep_sct_find_markers before differential expression on an object that carries more than one SCT model — the state you get by SCTransforming several objects separately and merging them. Without it, a fold change across the merge partly measures how deeply each batch was sequenced.

prep_sct_find_markers

prep_sct_find_markers(seurat, assay: str = 'SCT', umi_assay: Optional[str] = None, verbose: bool = True)

Put a merged object's SCT counts on one scale, before differential expression.

Mirrors R's PrepSCTFindMarkers(object). Run it once on an object that carries more than one SCT model — the state you get by SCTransforming several objects separately and merging them — and before any find_markers call on the SCT assay.

Why it is needed

sctransform corrects each object's counts to that object's median sequencing depth. Merge two such objects and the two halves of the SCT counts layer are expressed at two different depths, so a fold change across them partly measures how deeply each batch happened to be sequenced. This re-corrects every cell to the minimum median UMI across the models, which is the deepest common scale all of them can reach without extrapolating.

What it does

For each model, the Pearson residual is taken at each cell's own depth from the raw UMI counts (never from the already-corrected ones, so the operation is idempotent) and read back out at log10(min_median_umi). counts becomes the recorrected matrix and data its log1p.

It returns early, unchanged, in the two cases R does: when only one model is stored, and when every model's recorded median_umi already sits above the minimum observed one.

Parameters:

• assay (str, default: 'SCT' ) –

  the SCT assay to re-correct (default "SCT").

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

  assay holding the raw counts. Defaults to the one each model recorded at fit time, and raises if the models disagree — as R does, since a single corrected matrix cannot come from two different count matrices.

Returns:

• ``seurat``, with the SCT assay's ``counts`` and ``data`` layers replaced. –

Notes

Scale factor. The target depth is log10(median(umi)) — the median on the count scale — where an ordinary sctransform call uses median(log10(umi)). That asymmetry is R's, not a slip: correct_counts takes the median of the latent variable when no scale_factor is passed and log10 of the supplied one when there is. The two differ whenever a model has an even number of cells.

Source code in truecell/sctransform.py

def prep_sct_find_markers(
    seurat,
    assay: str = "SCT",
    umi_assay: Optional[str] = None,
    verbose: bool = True,
):
    """Put a merged object's SCT counts on one scale, before differential expression.

    Mirrors R's ``PrepSCTFindMarkers(object)``. Run it once on an object that
    carries **more than one** SCT model — the state you get by SCTransforming
    several objects separately and merging them — and before any
    ``find_markers`` call on the SCT assay.

    Why it is needed
    ----------------
    `sctransform` corrects each object's counts to *that object's* median
    sequencing depth. Merge two such objects and the two halves of the SCT
    ``counts`` layer are expressed at two different depths, so a fold change
    across them partly measures how deeply each batch happened to be
    sequenced. This re-corrects every cell to the **minimum** median UMI across
    the models, which is the deepest common scale all of them can reach without
    extrapolating.

    What it does
    ------------
    For each model, the Pearson residual is taken at each cell's own depth from
    the **raw** UMI counts (never from the already-corrected ones, so the
    operation is idempotent) and read back out at ``log10(min_median_umi)``.
    ``counts`` becomes the recorrected matrix and ``data`` its ``log1p``.

    It returns early, unchanged, in the two cases R does: when only one model is
    stored, and when every model's recorded ``median_umi`` already sits above
    the minimum observed one.

    Parameters
    ----------
    assay     : the SCT assay to re-correct (default ``"SCT"``).
    umi_assay : assay holding the raw counts. Defaults to the one each model
                recorded at fit time, and raises if the models disagree — as
                R does, since a single corrected matrix cannot come from two
                different count matrices.

    Returns
    -------
    ``seurat``, with the SCT assay's ``counts`` and ``data`` layers replaced.

    Notes
    -----
    **Scale factor.** The target depth is ``log10(median(umi))`` — the median on
    the count scale — where an ordinary `sctransform` call uses
    ``median(log10(umi))``. That asymmetry is R's, not a slip: ``correct_counts``
    takes the median of the latent variable when no ``scale_factor`` is passed
    and ``log10`` of the supplied one when there is. The two differ whenever a
    model has an even number of cells.
    """
    sct = seurat.assays[assay]
    models = dict(getattr(sct, "misc", {}).get("SCTModel.list") or {})

    if len(models) <= 1:
        if verbose:
            print("Only one SCT model is stored - skipping recalculating "
                  "corrected counts")
        return seurat

    # R keeps these apart on purpose. `observed` is recomputed from each model's
    # cell attributes every call; `stored` is the depth the counts were last
    # corrected *to*, which this function overwrites. They coincide until the
    # first run.
    observed = {n: float(np.median(m["cell_attributes"]["umi"]))
                for n, m in models.items()}
    min_median_umi = float(min(observed.values()))
    if all(float(m["median_umi"]) > min_median_umi for m in models.values()):
        if verbose:
            print("Minimum UMI unchanged. Skipping re-correction.")
        return seurat

    if umi_assay is None:
        named = {m.get("umi_assay") for m in models.values()}
        named.discard(None)
        if len(named) > 1:
            raise ValueError(
                "Multiple UMI assays are used for SCTransform: "
                + ", ".join(sorted(str(n) for n in named))
            )
        umi_assay = named.pop() if named else "RNA"
    if umi_assay not in seurat.assays:
        raise KeyError(
            f"assay {umi_assay!r} holds the raw counts these SCT models were "
            f"fitted from, but the object has no such assay"
        )

    if verbose:
        print(f"Found {len(models)} SCT models. Recorrecting SCT counts using "
              f"minimum median counts: {min_median_umi}")

    src = seurat.assays[umi_assay]
    # `: Any` for the same reason the fit path above declares it — the layer may
    # be dense, sparse or an on-disk LazyMatrix, and `issparse` is what picks.
    raw: Any = src.layers.get("counts") if isinstance(src, Assay5) else src.counts
    if raw is None:
        raise ValueError(f"assay {umi_assay!r} has no counts layer.")
    raw = raw.tocsc() if sp.issparse(raw) else sp.csc_matrix(raw)
    raw_genes = (src._all_feature_names if isinstance(src, Assay5)
                 else src._feature_names)
    raw_gene_pos = {g: i for i, g in enumerate(raw_genes)}
    raw_cell_pos = {c: i for i, c in enumerate(
        src.cells() if isinstance(src, Assay5) else src._cell_names)}

    all_genes = sct.features()
    all_cells = sct.cells()
    gene_pos = {g: i for i, g in enumerate(all_genes)}
    cell_pos = {c: i for i, c in enumerate(all_cells)}

    target = float(np.log10(min_median_umi))
    rows: list[np.ndarray] = []
    cols: list[np.ndarray] = []
    vals: list[np.ndarray] = []

    for model in models.values():
        fa = model["feature_attributes"]
        ca = model["cell_attributes"]
        # R: theta not NA and both coefficients finite. A gene the regularization
        # could not fit has no model to re-correct through.
        ok = (fa["theta"].notna()
              & np.isfinite(fa["(Intercept)"])
              & np.isfinite(fa["log_umi"]))
        genes = [g for g in fa.index[ok]
                 if g in raw_gene_pos and g in gene_pos]
        cells = [c for c in ca.index if c in raw_cell_pos and c in cell_pos]
        if not genes or not cells:
            continue

        sub = raw[[raw_gene_pos[g] for g in genes], :][:, [raw_cell_pos[c] for c in cells]]
        corrected = _correct_counts(
            sub.tocsr(),
            fa.loc[genes, "(Intercept)"].to_numpy(dtype=float),
            fa.loc[genes, "log_umi"].to_numpy(dtype=float),
            fa.loc[genes, "theta"].to_numpy(dtype=float),
            ca.loc[cells, "log_umi"].to_numpy(dtype=float),
            target,
        )
        coo = corrected.tocoo()
        finite = np.isfinite(coo.data)
        gidx = np.array([gene_pos[g] for g in genes])
        cidx = np.array([cell_pos[c] for c in cells])
        rows.append(gidx[coo.row[finite]])
        cols.append(cidx[coo.col[finite]])
        vals.append(coo.data[finite])
        model["median_umi"] = min_median_umi

    corrected_full = sp.csc_matrix(
        (np.concatenate(vals), (np.concatenate(rows), np.concatenate(cols))),
        shape=(len(all_genes), len(all_cells)),
    )
    sct.set_layer_data("counts", corrected_full)
    sct.set_layer_data("data", _log1p_sparse(corrected_full))
    sct.misc["SCTModel.list"] = models
    return seurat