Differential expression

find_markers implements all eight of Seurat's tests: wilcox (tie-corrected), t, bimod, LR, negbinom, roc, mast and deseq2. Seven of them are per-cell and reproduce Seurat's top 50 genes exactly on PBMC 3k; deseq2 is pseudobulk and deliberately does not, because it is answering a different question — see the DE vignette.

Two numbers to know before reading a result table:

• avg_log2FC carries Seurat's pseudocount on the group sum, not the group mean. Getting that backwards shifts every fold change and also changes which genes clear logfc_threshold, so it silently changes the returned gene set, not just a column.
• pct.1 and pct.2 are rounded to three decimals, by Seurat, inside FindMarkers. Anything comparing two runs gene-by-gene should not expect them closer than 5e-4.

Per-cluster and per-pair tests

find_markers

find_markers(seurat, ident_1: Union[str, list[str]], ident_2: Optional[Union[str, list[str]]] = None, assay: Optional[str] = None, layer: Optional[str] = None, test_use: str = 'wilcox', only_pos: bool = False, min_pct: float = 0.1, logfc_threshold: float = 0.25, features: Optional[list[str]] = None, latent_vars: Optional[list[str]] = None, sample_col: Optional[str] = None, max_cells_per_ident: Optional[int] = None, random_seed: int = 1) -> DataFrame

Find differentially expressed marker genes.

Mirrors R's FindMarkers(pbmc, ident.1 = 2).

Parameters:

• ident_1 (Union[str, list[str]]) –

  cluster label(s) for group 1

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

  cluster label(s) for group 2 (None = all others)

• test_use (str, default: 'wilcox' ) –

  statistical test — 'wilcox' (default), 't', 'bimod' (McDavid 2013 bimodal LRT), 'LR' (logistic-regression LRT), 'negbinom' (negative-binomial GLM LRT on counts), 'mast' (MAST two-part hurdle LRT on log-normalized data), 'deseq2' (pseudobulk DESeq2 — sums counts per sample then tests sample-level, requires sample_col; needs pip install truecell[deseq2]), or 'roc' (AUC classifier power).

• only_pos (bool, default: False ) –

  only return positive markers

• min_pct (float, default: 0.1 ) –

  minimum fraction cells expressing gene in either group

• logfc_threshold (float, default: 0.25 ) –

  minimum log2 fold-change filter

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

  restrict to these genes (default: all)

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

  metadata columns to regress out as covariates in the 'LR', 'negbinom', and 'mast' models (Seurat's latent.vars). Note that Seurat's MASTDETest fits ~ condition alone — it adds no cellular detection rate term unless you pass one — so leaving this empty is what matches Seurat's default. Passing CDR is the MAST paper's advice, and a deliberate departure from Seurat.

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

  metadata column identifying pseudobulk replicates (donor / sample); required for test_use='deseq2', ignored otherwise.

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

  downsample each group to this many cells

Returns:

• For 'wilcox' / 't' / 'LR' / 'negbinom': DataFrame with columns –
• p_val, avg_log2FC, pct.1, pct.2, p_val_adj (sorted by p_val). –
• For 'roc': columns myAUC, avg_diff, power, avg_log2FC, pct.1, pct.2 –
• (sorted by power), with no p-value — matching Seurat. –

Source code in truecell/markers.py

def find_markers(
    seurat,
    ident_1: Union[str, list[str]],
    ident_2: Optional[Union[str, list[str]]] = None,
    assay: Optional[str] = None,
    layer: Optional[str] = None,
    test_use: str = "wilcox",
    only_pos: bool = False,
    min_pct: float = 0.1,
    logfc_threshold: float = 0.25,
    features: Optional[list[str]] = None,
    latent_vars: Optional[list[str]] = None,
    sample_col: Optional[str] = None,
    max_cells_per_ident: Optional[int] = None,
    random_seed: int = 1,
) -> pd.DataFrame:
    """Find differentially expressed marker genes.

    Mirrors R's FindMarkers(pbmc, ident.1 = 2).

    Parameters
    ----------
    ident_1         : cluster label(s) for group 1
    ident_2         : cluster label(s) for group 2 (None = all others)
    test_use        : statistical test — 'wilcox' (default), 't', 'bimod'
                      (McDavid 2013 bimodal LRT), 'LR' (logistic-regression LRT),
                      'negbinom' (negative-binomial GLM LRT on counts), 'mast'
                      (MAST two-part hurdle LRT on log-normalized data), 'deseq2'
                      (pseudobulk DESeq2 — sums counts per sample then tests
                      sample-level, requires ``sample_col``; needs
                      ``pip install truecell[deseq2]``), or 'roc' (AUC classifier
                      power).
    only_pos        : only return positive markers
    min_pct         : minimum fraction cells expressing gene in either group
    logfc_threshold : minimum log2 fold-change filter
    features        : restrict to these genes (default: all)
    latent_vars     : metadata columns to regress out as covariates in the
                      'LR', 'negbinom', and 'mast' models (Seurat's latent.vars).
                      Note that Seurat's ``MASTDETest`` fits ``~ condition``
                      alone — it adds **no** cellular detection rate term unless
                      you pass one — so leaving this empty is what matches
                      Seurat's default. Passing CDR is the MAST paper's advice,
                      and a deliberate departure from Seurat.
    sample_col      : metadata column identifying pseudobulk replicates (donor /
                      sample); required for ``test_use='deseq2'``, ignored
                      otherwise.
    max_cells_per_ident : downsample each group to this many cells

    Returns
    -------
    For 'wilcox' / 't' / 'LR' / 'negbinom': DataFrame with columns
    p_val, avg_log2FC, pct.1, pct.2, p_val_adj (sorted by p_val).
    For 'roc': columns myAUC, avg_diff, power, avg_log2FC, pct.1, pct.2
    (sorted by power), with no p-value — matching Seurat.
    """
    assay_name = assay or seurat.active_assay
    assay_obj = seurat.assays[assay_name]
    cells = seurat.cell_names()
    idents = list(seurat.idents)

    # Resolve ident strings
    ident_1_set = {str(ident_1)} if isinstance(ident_1, str) else {str(i) for i in ident_1}
    if ident_2 is None:
        ident_2_set = {str(i) for i in set(idents) if str(i) not in ident_1_set}
    else:
        ident_2_set = {str(ident_2)} if isinstance(ident_2, str) else {str(i) for i in ident_2}

    # Cell indices for each group
    cells_1 = [c for c, i in zip(cells, idents) if str(i) in ident_1_set]
    cells_2 = [c for c, i in zip(cells, idents) if str(i) in ident_2_set]

    if not cells_1:
        raise ValueError(f"No cells found with ident {ident_1}.")
    if not cells_2:
        # `ident_2=None` means "every other ident", so name what was searched
        # rather than a parameter the caller never passed.
        raise ValueError(
            f"No cells found with ident {ident_2}." if ident_2 is not None
            else f"No cells left outside ident {ident_1} to compare against."
        )

    # Optional downsampling
    if max_cells_per_ident is not None:
        rng = np.random.default_rng(random_seed)
        if len(cells_1) > max_cells_per_ident:
            cells_1 = list(rng.choice(cells_1, max_cells_per_ident, replace=False))
        if len(cells_2) > max_cells_per_ident:
            cells_2 = list(rng.choice(cells_2, max_cells_per_ident, replace=False))

    # Get expression matrix for all genes (features × cells)
    data, feature_names = _get_expression_matrix(assay_obj, layer)

    cell_idx_map = {c: i for i, c in enumerate(cells)}
    idx_1 = [cell_idx_map[c] for c in cells_1]
    idx_2 = [cell_idx_map[c] for c in cells_2]

    # Column-slice each group and leave it sparse. A lazy layer indexes like a
    # sparse one and hands back a scipy CSC, so it takes the same branch -- the
    # dense fallback would materialise the whole store twice, once per group.
    #
    # Nothing here is densified. Both pre-filters below (`min_pct` and
    # `logfc_threshold`) are computable on the sparse matrices, and they are
    # what decides the handful of genes that a dense array is finally built for.
    if sp.issparse(data) or is_lazy(data):
        sub1 = data[:, idx_1]  # (features × n1)
        sub2 = data[:, idx_2]  # (features × n2)
    else:
        sub1 = np.asarray(data)[:, idx_1].astype(float)
        sub2 = np.asarray(data)[:, idx_2].astype(float)
    n1, n2 = sub1.shape[1], sub2.shape[1]

    # Restrict features
    if features is not None:
        feat_set = set(features)
        feat_mask = np.array([f in feat_set for f in feature_names])
    else:
        feat_mask = np.ones(len(feature_names), dtype=bool)

    # Percent cells expressing (> 0)
    pct1 = _row_pct_positive(sub1)
    pct2 = _row_pct_positive(sub2)

    # Pre-filter: gene must be expressed in at least min_pct of either group
    pct_mask = (pct1 >= min_pct) | (pct2 >= min_pct)
    combined_mask = feat_mask & pct_mask

    # Log2 fold change, matching Seurat 5's `log1pdata.mean.fxn` exactly:
    #   log2((sum(expm1(x)) + pseudocount) / n)
    # Data is log1p-normalized, so each cell is un-logged (expm1) before averaging.
    #
    # The pseudocount goes on the **sum**, not on the mean — it is one count added
    # to the whole group, worth 1/n on the mean scale, not 1. Adding it to the mean
    # instead (Seurat 4's `rowMeans(expm1(x)) + pseudocount`) floors every fold
    # change near zero: a gene seen in 0 % of one group and 24 % of the other
    # reads -1.26 that way against Seurat 5's -9.92. On pbmc3k that moved 98.9 %
    # of genes, and because `logfc_threshold` filters on this value it also
    # changed *which* genes came back — 2,298 against Seurat's 11,931 at 0.25.
    #
    # NOTE: the mean must be taken AFTER expm1, not before — expm1(mean(x)) is the
    # geometric-style mean and systematically compresses fold-changes (Jensen).
    group1_mean = (_row_expm1_sum(sub1) + PSEUDOCOUNT) / n1
    group2_mean = (_row_expm1_sum(sub2) + PSEUDOCOUNT) / n2
    avg_log2fc = np.log2(group1_mean) - np.log2(group2_mean)

    # Pre-filter by logfc_threshold
    if logfc_threshold > 0:
        fc_mask_arr = np.abs(avg_log2fc) >= logfc_threshold
        combined_mask = combined_mask & fc_mask_arr

    test_indices = np.where(combined_mask)[0]

    # The pre-filters are done, so this is where the dense arrays the per-gene
    # tests need get built — for the surviving genes only. `deseq2` is excluded
    # because it never looks at them; it aggregates the counts layer itself.
    if test_use != "deseq2" and len(test_indices) > 0:
        mat1 = _dense_rows(sub1, test_indices)  # (tested genes × n1)
        mat2 = _dense_rows(sub2, test_indices)  # (tested genes × n2)
    # Nothing below reads the sparse slices, and the per-gene loops they would
    # otherwise sit through run for seconds. Between them they hold about one
    # copy of the layer, so dropping them here is worth the line.
    del sub1, sub2

    # Per-cell covariates for the regression-based tests (LR / negbinom).
    latent = None
    if latent_vars and test_use in ("LR", "negbinom", "mast"):
        lat1 = seurat.meta_data.loc[cells_1, latent_vars].to_numpy(dtype=float)
        lat2 = seurat.meta_data.loc[cells_2, latent_vars].to_numpy(dtype=float)
        latent = np.vstack([lat1, lat2])

    # ---- ROC test: returns AUC / power, no p-value (matches Seurat) ----------
    if test_use == "roc":
        if len(test_indices) == 0:
            return pd.DataFrame(
                columns=["myAUC", "avg_diff", "power", "avg_log2FC", "pct.1", "pct.2"]
            )
        aucs = np.empty(len(test_indices))
        powers = np.empty(len(test_indices))
        avg_diff = np.empty(len(test_indices))
        for i in range(len(test_indices)):
            auc, power = _roc_auc(mat1[i, :], mat2[i, :])
            aucs[i] = auc
            powers[i] = power
            avg_diff[i] = mat1[i, :].mean() - mat2[i, :].mean()
        roc_res = pd.DataFrame(
            {
                "myAUC": aucs,
                "avg_diff": avg_diff,
                "power": powers,
                "avg_log2FC": avg_log2fc[test_indices],
                "pct.1": pct1[test_indices],
                "pct.2": pct2[test_indices],
            },
            index=[feature_names[i] for i in test_indices],
        )
        if only_pos:
            roc_res = roc_res[roc_res["avg_log2FC"] > 0]
        return roc_res.sort_values("power", ascending=False)

    if len(test_indices) == 0:
        return pd.DataFrame(
            columns=["p_val", "avg_log2FC", "pct.1", "pct.2", "p_val_adj"]
        )

    # ---- pseudobulk DESeq2: sample-level test, not per-cell -------------------
    if test_use == "deseq2":
        return _deseq2_pseudobulk(
            seurat, assay_obj, cells_1, cells_2, sample_col,
            feature_names, test_indices, pct1, pct2, only_pos,
        )

    # ---- p-value-based tests -------------------------------------------------
    p_vals = np.ones(len(test_indices))

    if test_use == "wilcox":
        for i in range(len(test_indices)):
            x1 = mat1[i, :]
            x2 = mat2[i, :]
            if x1.sum() == 0 and x2.sum() == 0:
                p_vals[i] = 1.0
            else:
                # mannwhitneyu (asymptotic) applies the tie correction and
                # continuity correction that base-R wilcox.test / presto use —
                # essential for scRNA data, which is dominated by zero ties.
                # scipy.stats.ranksums does NOT correct for ties.
                try:
                    _, p = mannwhitneyu(
                        x1, x2, alternative="two-sided",
                        use_continuity=True, method="asymptotic",
                    )
                except ValueError:
                    # Raised only when every value in both groups is identical.
                    p = 1.0
                p_vals[i] = p if not np.isnan(p) else 1.0
    elif test_use == "t":
        from scipy.stats import ttest_ind
        for i in range(len(test_indices)):
            x1 = mat1[i, :]
            x2 = mat2[i, :]
            _, p = ttest_ind(x1, x2, equal_var=False)
            p_vals[i] = p if not np.isnan(p) else 1.0
    elif test_use == "bimod":
        for i in range(len(test_indices)):
            p_vals[i] = _bimod_pvalue(mat1[i, :], mat2[i, :])
    elif test_use == "LR":
        group = np.concatenate([np.ones(n1), np.zeros(n2)])
        for i in range(len(test_indices)):
            expr = np.concatenate([mat1[i, :], mat2[i, :]])
            p_vals[i] = _lr_pvalue(expr, group, latent)
    elif test_use == "mast":
        group = np.concatenate([np.ones(n1), np.zeros(n2)])
        for i in range(len(test_indices)):
            expr = np.concatenate([mat1[i, :], mat2[i, :]])
            p_vals[i] = _mast_pvalue(expr, group, latent)
    elif test_use == "negbinom":
        # Counts, not the data layer — and restricted to the tested genes for
        # the same reason as above.
        counts_mat, _ = _get_expression_matrix(assay_obj, "counts")
        if sp.issparse(counts_mat) or is_lazy(counts_mat):
            c1 = _dense_rows(counts_mat[:, idx_1], test_indices)
            c2 = _dense_rows(counts_mat[:, idx_2], test_indices)
        else:
            dense_counts = np.asarray(counts_mat)
            c1 = dense_counts[np.ix_(test_indices, np.asarray(idx_1))]
            c2 = dense_counts[np.ix_(test_indices, np.asarray(idx_2))]
        group = np.concatenate([np.ones(n1), np.zeros(n2)])
        for i in range(len(test_indices)):
            cnts = np.concatenate([c1[i, :], c2[i, :]])
            p_vals[i] = _negbinom_pvalue(cnts, group, latent)
    else:
        raise ValueError(
            f"Unsupported test_use: {test_use!r}. "
            "Use 'wilcox', 't', 'bimod', 'LR', 'negbinom', 'mast', 'deseq2', or 'roc'."
        )

    # Bonferroni correction (Seurat default: multiply by total gene count)
    n_total = len(feature_names)
    p_val_adj = np.minimum(p_vals * n_total, 1.0)

    results = pd.DataFrame(
        {
            "p_val": p_vals,
            "avg_log2FC": avg_log2fc[test_indices],
            "pct.1": pct1[test_indices],
            "pct.2": pct2[test_indices],
            "p_val_adj": p_val_adj,
        },
        index=[feature_names[i] for i in test_indices],
    )

    if only_pos:
        results = results[results["avg_log2FC"] > 0]

    return results.sort_values("p_val")

find_all_markers

find_all_markers(seurat, assay: Optional[str] = None, layer: Optional[str] = None, test_use: str = 'wilcox', only_pos: bool = False, min_pct: float = 0.1, logfc_threshold: float = 0.25, sample_col: Optional[str] = None, max_cells_per_ident: Optional[int] = None, random_seed: int = 1, return_thresh: float = 0.01) -> DataFrame

Find marker genes for each cluster vs all others.

Mirrors R's FindAllMarkers(pbmc, only.pos = TRUE).

Returns a single DataFrame with an extra 'cluster' column.

return_thresh is Seurat's return.thresh: only genes with p_val < return_thresh are returned (for test_use="roc", only genes whose myAUC is further than return_thresh from 0.5 in either direction, since ROC reports no p-value). Pass None for the unfiltered table. Without it truecell returned every gene that survived the pct and logfc pre-filters, including plainly non-significant ones: on PBMC 3k that was 3,036 rows against Seurat's 3,446 spread over one fewer cluster, and on the two clusters whose membership matched Seurat exactly the filtered table reproduces Seurat's gene set exactly (151 and 242 genes).

Rows are ordered by p_val ascending and then avg_log2FC descending within each cluster, matching Seurat's order(gde$p_val, -gde[, 2]). The tie-break matters: Wilcoxon p-values tie at 0 for the strongest markers, so without it "the top 10 markers" depends on incoming row order.

Source code in truecell/markers.py

def find_all_markers(
    seurat,
    assay: Optional[str] = None,
    layer: Optional[str] = None,
    test_use: str = "wilcox",
    only_pos: bool = False,
    min_pct: float = 0.1,
    logfc_threshold: float = 0.25,
    sample_col: Optional[str] = None,
    max_cells_per_ident: Optional[int] = None,
    random_seed: int = 1,
    return_thresh: float = 1e-2,
) -> pd.DataFrame:
    """Find marker genes for each cluster vs all others.

    Mirrors R's FindAllMarkers(pbmc, only.pos = TRUE).

    Returns a single DataFrame with an extra 'cluster' column.

    ``return_thresh`` is Seurat's ``return.thresh``: only genes with
    ``p_val < return_thresh`` are returned (for ``test_use="roc"``, only genes
    whose ``myAUC`` is further than ``return_thresh`` from 0.5 in either
    direction, since ROC reports no p-value). Pass ``None`` for the unfiltered
    table. Without it truecell returned every gene that survived the pct and
    logfc pre-filters, including plainly non-significant ones: on PBMC 3k that
    was 3,036 rows against Seurat's 3,446 spread over one fewer cluster, and on
    the two clusters whose membership matched Seurat exactly the filtered table
    reproduces Seurat's gene set exactly (151 and 242 genes).

    Rows are ordered by ``p_val`` ascending and then ``avg_log2FC`` descending
    within each cluster, matching Seurat's ``order(gde$p_val, -gde[, 2])``.
    The tie-break matters: Wilcoxon p-values tie at 0 for the strongest
    markers, so without it "the top 10 markers" depends on incoming row order.
    """
    clusters = sorted(set(str(i) for i in seurat.idents), key=_ident_sort_key)
    all_results = []

    for cluster in clusters:
        try:
            df = find_markers(
                seurat,
                ident_1=cluster,
                ident_2=None,
                assay=assay,
                layer=layer,
                test_use=test_use,
                only_pos=only_pos,
                min_pct=min_pct,
                logfc_threshold=logfc_threshold,
                sample_col=sample_col,
                max_cells_per_ident=max_cells_per_ident,
                random_seed=random_seed,
            )
            if len(df) > 0:
                df = df.copy()
                df["cluster"] = cluster
                df["gene"] = df.index
                all_results.append(df)
        except Exception as e:
            print(f"Warning: cluster {cluster} marker finding failed: {e}")

    if not all_results:
        return pd.DataFrame(
            columns=["p_val", "avg_log2FC", "pct.1", "pct.2", "p_val_adj", "cluster", "gene"]
        )

    combined = pd.concat(all_results, axis=0)
    if return_thresh is not None:
        if test_use == "roc":
            # ROC has no p-value; Seurat thresholds on distance from chance.
            auc = combined["myAUC"]
            combined = combined[(auc > return_thresh) | (auc < 1 - return_thresh)]
        else:
            combined = combined[combined["p_val"] < return_thresh]
    if test_use == "roc":
        cols = ["cluster", "gene", "myAUC", "avg_diff", "power", "avg_log2FC",
                "pct.1", "pct.2"]
        cols = [c for c in cols if c in combined.columns]
        return combined[cols].sort_values(["cluster", "myAUC"], ascending=[True, False])
    combined = combined[["cluster", "gene", "p_val", "avg_log2FC", "pct.1", "pct.2", "p_val_adj"]]
    return combined.sort_values(["cluster", "p_val", "avg_log2FC"],
                                ascending=[True, True, False])

find_conserved_markers

find_conserved_markers(seurat, ident_1: Union[str, list[str]], grouping_var: str, ident_2: Optional[Union[str, list[str]]] = None, assay: Optional[str] = None, layer: Optional[str] = None, test_use: str = 'wilcox', only_pos: bool = False, min_pct: float = 0.1, logfc_threshold: float = 0.25, features: Optional[list[str]] = None) -> DataFrame

Find markers conserved across the levels of a grouping variable.

Mirrors R's FindConservedMarkers(obj, ident.1, grouping.var = "stim"): runs find_markers for ident_1 vs ident_2 independently within each level of grouping_var, keeps only genes detected as markers in every level, and combines their per-level p-values with Fisher's method (scipy.stats.combine_pvalues).

Every argument not listed below is forwarded verbatim to find_markers.

Parameters:

• ident_1 (Union[str, list[str]]) –

  cluster label(s) for group 1.

• grouping_var (str) –

  metadata column whose levels define the independent comparisons (e.g. condition, batch, donor).

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

  cluster label(s) for group 2 (None = all other cells).

Returns:

• DataFrame indexed by gene with, for each level ``g``, the columns –
• ``{g}_p_val, {g}_avg_log2FC, {g}_pct.1, {g}_pct.2, {g}_p_val_adj`` plus –
• ``max_pval`` (worst per-level p-value) and ``combined_p_val`` (Fisher-combined –
• across levels), sorted by ``combined_p_val``. Only genes that are markers in –
• all levels are returned. –

Source code in truecell/markers.py

def find_conserved_markers(
    seurat,
    ident_1: Union[str, list[str]],
    grouping_var: str,
    ident_2: Optional[Union[str, list[str]]] = None,
    assay: Optional[str] = None,
    layer: Optional[str] = None,
    test_use: str = "wilcox",
    only_pos: bool = False,
    min_pct: float = 0.1,
    logfc_threshold: float = 0.25,
    features: Optional[list[str]] = None,
) -> pd.DataFrame:
    """Find markers conserved across the levels of a grouping variable.

    Mirrors R's ``FindConservedMarkers(obj, ident.1, grouping.var = "stim")``:
    runs :func:`find_markers` for ``ident_1`` vs ``ident_2`` independently within
    each level of ``grouping_var``, keeps only genes detected as markers in
    *every* level, and combines their per-level p-values with Fisher's method
    (:func:`scipy.stats.combine_pvalues`).

    Every argument not listed below is forwarded verbatim to
    :func:`find_markers`.

    Parameters
    ----------
    ident_1      : cluster label(s) for group 1.
    grouping_var : metadata column whose levels define the independent
                   comparisons (e.g. condition, batch, donor).
    ident_2      : cluster label(s) for group 2 (None = all other cells).

    Returns
    -------
    DataFrame indexed by gene with, for each level ``g``, the columns
    ``{g}_p_val, {g}_avg_log2FC, {g}_pct.1, {g}_pct.2, {g}_p_val_adj`` plus
    ``max_pval`` (worst per-level p-value) and ``combined_p_val`` (Fisher-combined
    across levels), sorted by ``combined_p_val``. Only genes that are markers in
    all levels are returned.
    """
    from scipy.stats import combine_pvalues

    if grouping_var not in seurat.meta_data.columns:
        raise KeyError(
            f"grouping_var {grouping_var!r} not found in meta_data "
            f"(columns: {list(seurat.meta_data.columns)})."
        )

    cells = seurat.cell_names()
    group_of = seurat.meta_data.loc[cells, grouping_var].astype(str)
    levels = sorted(group_of.unique())

    per_level: dict[str, pd.DataFrame] = {}
    for level in levels:
        level_cells = [c for c, g in zip(cells, group_of) if g == level]
        sub = seurat.subset(cells=level_cells)
        try:
            df = find_markers(
                sub,
                ident_1=ident_1,
                ident_2=ident_2,
                assay=assay,
                layer=layer,
                test_use=test_use,
                only_pos=only_pos,
                min_pct=min_pct,
                logfc_threshold=logfc_threshold,
                features=features,
            )
        except ValueError as e:
            warnings.warn(
                f"Skipping {grouping_var}={level!r}: {e}", RuntimeWarning, stacklevel=2
            )
            continue
        if len(df) > 0:
            per_level[level] = df

    if not per_level:
        raise ValueError(
            f"No level of {grouping_var!r} yielded markers for the requested comparison."
        )

    # Genes must be markers in every retained level.
    common = set.intersection(*(set(df.index) for df in per_level.values()))

    used = list(per_level)
    cols: dict[str, pd.Series] = {}
    for level in used:
        df = per_level[level].loc[list(common)]
        for c in df.columns:
            cols[f"{level}_{c}"] = df[c]
    result = pd.DataFrame(cols, index=list(common))

    if test_use == "roc":
        # ROC has no p-value; conservation is summarised by the min power.
        power_cols = [f"{level}_power" for level in used]
        result["min_power"] = result[power_cols].min(axis=1)
        return result.sort_values("min_power", ascending=False)

    pval_cols = [f"{level}_p_val" for level in used]
    result["max_pval"] = result[pval_cols].max(axis=1)
    if len(used) == 1:
        result["combined_p_val"] = result[pval_cols[0]]
    else:
        result["combined_p_val"] = [
            combine_pvalues(result.loc[g, pval_cols].to_numpy(dtype=float),
                            method="fisher").pvalue
            for g in result.index
        ]
    return result.sort_values("combined_p_val")

Group summaries

AggregateExpression sums raw counts and is what pseudobulk differential expression wants. AverageExpression means the back-transformed values and is what a per-group expression summary wants. They are different functions, not two scalings of one — see each docstring.

aggregate_expression

aggregate_expression(seurat, group_by: Union[str, list[str]] = 'ident', assays: Optional[Union[str, list[str]]] = None, features: Optional[list[str]] = None, layer: str = 'counts', return_object: bool = False, normalization_method: str = 'LogNormalize', scale_factor: float = 10000.0)

Sum counts within cell groups to form a pseudobulk profile.

Mirrors R's AggregateExpression(obj, group.by = c("celltype", "donor")).

Parameters:

• group_by (Union[str, list[str]], default: 'ident' ) –

  metadata column(s) defining the groups. Multiple columns are combined into a single label joined by "_" (as Seurat does). "ident" uses the object's active identities.

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

  assay name(s) to aggregate (default: the active assay).

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

  restrict to these features (default: all).

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

  layer to aggregate (default "counts" — pseudobulk is defined on raw counts). Seurat's AggregateExpression has no such argument and always sums counts; this is a superset with a matching default.

• return_object (bool, default: False ) –

  if True, return a new Truecell object with one "cell" per group; if False (default), return a pd.DataFrame (features × groups), or a dict of them when several assays are requested.

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

  how to fill the returned object's data layer when return_object=True. "LogNormalize" (Seurat's default) or None to leave data as the raw sums.

• scale_factor (float, default: 10000.0 ) –

  the scale factor for that normalization (Seurat: 10000).

Notes

return_object=True normalizes, which is easy to miss and was wrong here until it was checked against R: Seurat's return.seurat = TRUE runs NormalizeData over the pseudobulk, so data holds log1p(sums / colSums × 10000) and not the sums. Leaving the sums in data — which is what this did — hands every downstream function that reads that layer un-normalized library-size-confounded values.

This is the one place AggregateExpression and AverageExpression diverge on their object output: average_expression writes plain log1p of the averages, with no library-size step. Verified against Seurat 5.5.1 for both.

Returns:

• ``pd.DataFrame`` | ``dict[str, pd.DataFrame]`` | ``Truecell`` –

  A single DataFrame when one assay is aggregated, a dict keyed by assay name when several are, or a Truecell object when return_object=True.

Source code in truecell/aggregate.py

def aggregate_expression(
    seurat,
    group_by: Union[str, list[str]] = "ident",
    assays: Optional[Union[str, list[str]]] = None,
    features: Optional[list[str]] = None,
    layer: str = "counts",
    return_object: bool = False,
    normalization_method: str = "LogNormalize",
    scale_factor: float = 10000.0,
):
    """Sum counts within cell groups to form a pseudobulk profile.

    Mirrors R's ``AggregateExpression(obj, group.by = c("celltype", "donor"))``.

    Parameters
    ----------
    group_by      : metadata column(s) defining the groups. Multiple columns are
                    combined into a single label joined by ``"_"`` (as Seurat
                    does). ``"ident"`` uses the object's active identities.
    assays        : assay name(s) to aggregate (default: the active assay).
    features      : restrict to these features (default: all).
    layer         : layer to aggregate (default ``"counts"`` — pseudobulk is
                    defined on raw counts). Seurat's ``AggregateExpression`` has
                    no such argument and always sums ``counts``; this is a
                    superset with a matching default.
    return_object : if True, return a new :class:`Truecell` object with one "cell"
                    per group; if False (default), return a ``pd.DataFrame``
                    (features × groups), or a ``dict`` of them when several
                    assays are requested.
    normalization_method : how to fill the returned object's ``data`` layer when
                    ``return_object=True``. ``"LogNormalize"`` (Seurat's default)
                    or ``None`` to leave ``data`` as the raw sums.
    scale_factor  : the scale factor for that normalization (Seurat: 10000).

    Notes
    -----
    ``return_object=True`` **normalizes**, which is easy to miss and was wrong
    here until it was checked against R: Seurat's ``return.seurat = TRUE`` runs
    ``NormalizeData`` over the pseudobulk, so ``data`` holds
    ``log1p(sums / colSums × 10000)`` and not the sums. Leaving the sums in
    ``data`` — which is what this did — hands every downstream function that
    reads that layer un-normalized library-size-confounded values.

    This is the one place ``AggregateExpression`` and ``AverageExpression``
    diverge on their object output: :func:`average_expression` writes plain
    ``log1p`` of the averages, with no library-size step. Verified against
    Seurat 5.5.1 for both.

    Returns
    -------
    ``pd.DataFrame`` | ``dict[str, pd.DataFrame]`` | ``Truecell``
        A single DataFrame when one assay is aggregated, a dict keyed by assay
        name when several are, or a Truecell object when ``return_object=True``.
    """
    group_cols = [group_by] if isinstance(group_by, str) else list(group_by)
    if assays is None:
        assay_names = [seurat.active_assay]
    else:
        assay_names = [assays] if isinstance(assays, str) else list(assays)

    labels = _group_labels(seurat, group_cols)
    groups = sorted(labels.unique())
    n_cells = len(labels)
    n_groups = len(groups)

    # One-hot cells × groups indicator; counts(features×cells) @ indicator
    # sums each group's columns in a single sparse matmul.
    group_index = {g: j for j, g in enumerate(groups)}
    col_idx = np.fromiter((group_index[g] for g in labels), dtype=int, count=n_cells)
    indicator = sp.csr_matrix(
        (np.ones(n_cells), (np.arange(n_cells), col_idx)),
        shape=(n_cells, n_groups),
    )

    agg_frames: dict[str, pd.DataFrame] = {}
    for name in assay_names:
        assay_obj = seurat.assays[name]
        data, feature_names = _get_expression_matrix(assay_obj, layer)
        if features is not None:
            feat_set = set(features)
            keep = np.array([f in feat_set for f in feature_names])
            feature_names = [f for f, k in zip(feature_names, keep) if k]
            data = data[keep, :]
        summed = data @ indicator  # features × groups
        if sp.issparse(summed):
            summed = summed.toarray()
        agg_frames[name] = pd.DataFrame(
            np.asarray(summed), index=feature_names, columns=groups
        )

    if return_object:
        obj = _to_truecell(seurat, agg_frames, labels, groups, group_cols)
        if normalization_method is not None:
            from .preprocessing import normalize_data
            for name in assay_names:
                normalize_data(obj, assay=name,
                               normalization_method=normalization_method,
                               scale_factor=scale_factor)
        return obj

    if len(assay_names) == 1:
        return agg_frames[assay_names[0]]
    return agg_frames

average_expression

average_expression(seurat, group_by: Union[str, list[str]] = 'ident', assays: Optional[Union[str, list[str]]] = None, features: Optional[list[str]] = None, layer: str = 'data', return_object: bool = False)

Mean expression within cell groups.

Mirrors R's AverageExpression(obj, group.by = "celltype"), which is a different function from aggregate_expression and not a rescaling of it. Two things differ, and both matter.

It averages rather than sums, and on the data layer it averages the back-transformed values: mean(expm1(x)), not mean(x) and not expm1(mean(x)). The data layer holds log1p-normalized expression, so a mean taken on it is a geometric-ish mean of nothing in particular; Seurat undoes the log first, averages on the linear scale, and returns that. The difference is not cosmetic — on a small Poisson object the first gene reads 332.84 under AverageExpression against a count mean of 3.17.

The back-transform applies to the data layer only. counts and scale.data are not log-normalized, so those are averaged as they stand. Verified against Seurat 5.5.1 for all three layers.

Parameters:

• group_by (Union[str, list[str]], default: 'ident' ) –

  metadata column(s) defining the groups; "ident" uses the object's active identities. Several are joined by "_".

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

  assay name(s) to average (default: the active assay).

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

  restrict to these features (default: all).

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

  layer to average (default "data", as Seurat's is).

• return_object (bool, default: False ) –

  if True, return a Truecell with one "cell" per group. Seurat puts the averaged matrix in that object's counts layer and log1p of it in data; this does the same.

Returns:

• ``pd.DataFrame`` | ``dict[str, pd.DataFrame]`` | ``Truecell`` –

See Also

aggregate_expression : sums raw counts, which is what pseudobulk DE wants.

Source code in truecell/aggregate.py

def average_expression(
    seurat,
    group_by: Union[str, list[str]] = "ident",
    assays: Optional[Union[str, list[str]]] = None,
    features: Optional[list[str]] = None,
    layer: str = "data",
    return_object: bool = False,
):
    """Mean expression within cell groups.

    Mirrors R's ``AverageExpression(obj, group.by = "celltype")``, which is a
    different function from :func:`aggregate_expression` and not a rescaling of
    it. Two things differ, and both matter.

    **It averages rather than sums**, and **on the ``data`` layer it averages the
    back-transformed values**: ``mean(expm1(x))``, not ``mean(x)`` and not
    ``expm1(mean(x))``. The ``data`` layer holds log1p-normalized expression, so
    a mean taken on it is a geometric-ish mean of nothing in particular; Seurat
    undoes the log first, averages on the linear scale, and returns that. The
    difference is not cosmetic — on a small Poisson object the first gene reads
    **332.84** under ``AverageExpression`` against a count mean of **3.17**.

    The back-transform applies to the ``data`` layer **only**. ``counts`` and
    ``scale.data`` are not log-normalized, so those are averaged as they stand.
    Verified against Seurat 5.5.1 for all three layers.

    Parameters
    ----------
    group_by      : metadata column(s) defining the groups; ``"ident"`` uses the
                    object's active identities. Several are joined by ``"_"``.
    assays        : assay name(s) to average (default: the active assay).
    features      : restrict to these features (default: all).
    layer         : layer to average (default ``"data"``, as Seurat's is).
    return_object : if True, return a :class:`Truecell` with one "cell" per group.
                    Seurat puts the averaged matrix in that object's ``counts``
                    layer and ``log1p`` of it in ``data``; this does the same.

    Returns
    -------
    ``pd.DataFrame`` | ``dict[str, pd.DataFrame]`` | ``Truecell``

    See Also
    --------
    aggregate_expression : sums raw counts, which is what pseudobulk DE wants.
    """
    from .markers import expm1_keeping_sparsity

    group_cols = [group_by] if isinstance(group_by, str) else list(group_by)
    if assays is None:
        assay_names = [seurat.active_assay]
    else:
        assay_names = [assays] if isinstance(assays, str) else list(assays)

    labels = _group_labels(seurat, group_cols)
    groups = sorted(labels.unique())
    n_cells = len(labels)

    group_index = {g: j for j, g in enumerate(groups)}
    col_idx = np.fromiter((group_index[g] for g in labels), dtype=int, count=n_cells)
    indicator = sp.csr_matrix(
        (np.ones(n_cells), (np.arange(n_cells), col_idx)),
        shape=(n_cells, len(groups)),
    )
    # Divide by the group's own size, not by the mean size: the groups are
    # rarely balanced and Seurat averages within each.
    sizes = np.asarray(indicator.sum(axis=0)).ravel()

    avg_frames: dict[str, pd.DataFrame] = {}
    for name in assay_names:
        assay_obj = seurat.assays[name]
        data, feature_names = _get_expression_matrix(assay_obj, layer)
        if features is not None:
            feat_set = set(features)
            keep = np.array([f in feat_set for f in feature_names])
            feature_names = [f for f, k in zip(feature_names, keep) if k]
            data = data[keep, :]
        if layer == "data":
            data = expm1_keeping_sparsity(data)
        summed = data @ indicator
        if sp.issparse(summed):
            summed = summed.toarray()
        avg_frames[name] = pd.DataFrame(
            np.asarray(summed) / sizes, index=feature_names, columns=groups
        )

    if return_object:
        obj = _to_truecell(seurat, avg_frames, labels, groups, group_cols)
        # Seurat leaves the averaged values in `counts` and writes log1p of them
        # to `data`, rather than running NormalizeData over a matrix that is
        # already a per-group average.
        from .preprocessing import _set_layer
        for name in assay_names:
            _set_layer(obj.assays[name], "data",
                       np.log1p(avg_frames[name].to_numpy()))
        return obj

    if len(assay_names) == 1:
        return avg_frames[assay_names[0]]
    return avg_frames