Differential expression

find_markers implements all nine of Seurat's tests: wilcox (tie-corrected), t, bimod, LR, negbinom, poisson, roc, mast and deseq2. On PBMC 3k the eight that return a p-value reproduce Seurat's top 50 genes exactly, and roc its AUCs to Seurat's three-decimal rounding — see the DE vignette. deseq2 tests every cell as a replicate, as Seurat's DESeq2DETest does; sample_col sums each sample's cells first and makes it a pseudobulk test.

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 before min_pct filters on them, as Seurat's FoldChange rounds them, and by R's rule rather than NumPy's. A gene detected in 19 of 2,000 cells (0.0095) therefore passes min_pct=0.01 in both tools, and the two columns match Seurat's exactly.

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.01, logfc_threshold: float = 0.1, 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 Wald test on counts), 'poisson' (Poisson GLM Wald test on counts — fast, but anti-conservative on overdispersed UMI data; prefer 'negbinom'), 'mast' (MAST two-part hurdle LRT on log-normalized data), 'deseq2' (DESeq2 as Seurat runs it: size factors, local dispersion trend and Wald test on counts, every cell a replicate unless sample_col aggregates them; needs pip install truecell[deseq2]), or 'roc' (AUC classifier power).

• only_pos (bool, default: False ) –

  only return positive markers

• min_pct (float, default: 0.01 ) –

  minimum fraction of cells expressing the gene in either group, compared after rounding to three decimals as Seurat does. Seurat 5's default, 0.01.

• logfc_threshold (float, default: 0.1 ) –

  minimum absolute log2 fold change. Seurat 5's default, 0.1.

• 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', 'poisson' and 'mast' models — the same four Seurat's DEmethods_latent() names for its 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 ) –

  for test_use='deseq2', a metadata column naming each cell's replicate (donor / sample); counts are summed per replicate before the test. None (default) tests every cell as a replicate, as Seurat's DESeq2DETest does. Ignored by the other tests.

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

  downsample each group to this many cells

Returns:

• For 'wilcox' / 't' / 'bimod' / 'LR' / 'negbinom' / 'poisson' / 'mast' / 'deseq2': –
• DataFrame with columns p_val, avg_log2FC, pct.1, pct.2, p_val_adj, in –
• Seurat's order: by p_val, ties to the larger |pct.1 - pct.2|. –
• For 'roc': columns myAUC, avg_diff, power, avg_log2FC, pct.1, pct.2, by –
• power and then myAUC, both descending, with no p-value — matching Seurat. –

Notes

Divergence, 'negbinom' and 'poisson'. Seurat's GLMDETest drops a gene from its output when the gene is detected in fewer than min.cells (3) cells in both groups, or has zero variance across them; it flags those with a sentinel p-value of 2 and deletes the rows. Here they are returned with p_val = 1.0 instead — no evidence rather than no row — which keeps the frame's gene set identical across every test_use and keeps p_val a p-value. The min_pct pre-filter already removes most such genes before either rule could fire.

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.01,
    logfc_threshold: float = 0.1,
    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 Wald test on counts),
                      'poisson' (Poisson GLM Wald test on counts — fast, but
                      anti-conservative on overdispersed UMI data; prefer
                      'negbinom'), 'mast' (MAST two-part hurdle LRT on
                      log-normalized data), 'deseq2' (DESeq2 as Seurat runs it:
                      size factors, local dispersion trend and Wald test on
                      counts, every cell a replicate unless ``sample_col``
                      aggregates them; needs ``pip install truecell[deseq2]``),
                      or 'roc' (AUC classifier power).
    only_pos        : only return positive markers
    min_pct         : minimum fraction of cells expressing the gene in either
                      group, compared after rounding to three decimals as
                      Seurat does. Seurat 5's default, 0.01.
    logfc_threshold : minimum absolute log2 fold change. Seurat 5's default, 0.1.
    features        : restrict to these genes (default: all)
    latent_vars     : metadata columns to regress out as covariates in the
                      'LR', 'negbinom', 'poisson' and 'mast' models — the same
                      four Seurat's ``DEmethods_latent()`` names for its
                      ``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      : for ``test_use='deseq2'``, a metadata column naming each
                      cell's replicate (donor / sample); counts are summed per
                      replicate before the test. ``None`` (default) tests every
                      cell as a replicate, as Seurat's ``DESeq2DETest`` does.
                      Ignored by the other tests.
    max_cells_per_ident : downsample each group to this many cells

    Returns
    -------
    For 'wilcox' / 't' / 'bimod' / 'LR' / 'negbinom' / 'poisson' / 'mast' / 'deseq2':
    DataFrame with columns p_val, avg_log2FC, pct.1, pct.2, p_val_adj, in
    Seurat's order: by p_val, ties to the larger |pct.1 - pct.2|.
    For 'roc': columns myAUC, avg_diff, power, avg_log2FC, pct.1, pct.2, by
    power and then myAUC, both descending, with no p-value — matching Seurat.

    Notes
    -----
    **Divergence, 'negbinom' and 'poisson'.** Seurat's ``GLMDETest`` *drops* a
    gene from its output when the gene is detected in fewer than ``min.cells``
    (3) cells in **both** groups, or has zero variance across them; it flags
    those with a sentinel p-value of 2 and deletes the rows. Here they are
    returned with ``p_val = 1.0`` instead — no evidence rather than no row —
    which keeps the frame's gene set identical across every ``test_use`` and
    keeps ``p_val`` a p-value. The ``min_pct`` pre-filter already removes most
    such genes before either rule could fire.
    """
    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)

    # Fraction of cells expressing (> 0), rounded as Seurat's `FoldChange` rounds
    # it. `min_pct` below filters on the rounded values, so a gene at 19 of 2,000
    # cells (0.0095, which R rounds to 0.010) passes a 0.01 cutoff in both tools.
    pct1 = _round_like_r(_row_pct_positive(sub1))
    pct2 = _round_like_r(_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. Seurat skips it for DESeq2
    # (`DEmethods_noprefilter()`) and keeps min.pct.
    if logfc_threshold > 0 and test_use != "deseq2":
        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. Mirrors Seurat's
    # `DEmethods_latent()` — negbinom, poisson, MAST, LR.
    latent = None
    if latent_vars and test_use in ("LR", "negbinom", "poisson", "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]
        # Seurat: `order(-power, -myAUC)`.
        order = np.lexsort((-roc_res["myAUC"].to_numpy(float), -roc_res["power"].to_numpy(float)))
        return roc_res.iloc[order]

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

    # ---- 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 in ("negbinom", "poisson"):
        # Seurat's two GLMDETest families. Both read 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)])
        glm_test = _negbinom_pvalue if test_use == "negbinom" else _poisson_pvalue
        for i in range(len(test_indices)):
            cnts = np.concatenate([c1[i, :], c2[i, :]])
            p_vals[i] = glm_test(cnts, group, latent)
    elif test_use == "deseq2":
        p_vals = _deseq2_pvalues(
            seurat, assay_obj, cells_1, cells_2, idx_1, idx_2, sample_col,
            feature_names, test_indices,
        )
    else:
        raise ValueError(
            f"Unsupported test_use: {test_use!r}. Use 'wilcox', 't', 'bimod', "
            "'LR', 'negbinom', 'poisson', '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 _order_like_seurat(results)

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.01, logfc_threshold: float = 0.1, 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", which has no p-value, only genes with myAUC > return_thresh or myAUC < 1 - return_thresh). 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).

Within each cluster, rows follow Seurat 5.5.1's order(p_val, -abs(pct.1 - pct.2)), which is how find_markers returns them. 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. It used to be descending avg_log2FC here, older Seurat's -gde[, 2]; on a full tie Seurat 5 keeps feature order instead.

With test_use="roc", Seurat swaps the default return.thresh of 0.01, which means nothing for an AUC, for 0.7 before filtering, so only genes with myAUC above 0.7 or below 0.3 come back. It tests the value, not whether it was passed, so an explicit 0.01 is swapped too; truecell does the same.

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.01,
    logfc_threshold: float = 0.1,
    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"``, which has no
    p-value, only genes with ``myAUC > return_thresh`` or
    ``myAUC < 1 - return_thresh``). 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).

    Within each cluster, rows follow Seurat 5.5.1's
    ``order(p_val, -abs(pct.1 - pct.2))``, which is how :func:`find_markers`
    returns them. 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. It used to be descending ``avg_log2FC`` here, older Seurat's
    ``-gde[, 2]``; on a full tie Seurat 5 keeps feature order instead.

    With ``test_use="roc"``, Seurat swaps the default ``return.thresh`` of 0.01,
    which means nothing for an AUC, for 0.7 before filtering, so only genes with
    ``myAUC`` above 0.7 or below 0.3 come back. It tests the value, not whether
    it was passed, so an explicit 0.01 is swapped too; truecell does the same.
    """
    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"]
        )

    # Clusters in the order visited, each already in Seurat's row order from
    # find_markers: what FindAllMarkers' rbind produces. Sorting again here used
    # to replace that order with a different tie-break.
    combined = pd.concat(all_results, axis=0)
    if return_thresh is not None:
        if test_use == "roc":
            # ROC has no p-value, so Seurat filters on myAUC, swapping its default
            # threshold of 0.01 for 0.7 first.
            thresh = 0.7 if return_thresh == 0.01 else return_thresh
            auc = combined["myAUC"]
            combined = combined[(auc > thresh) | (auc < 1 - 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"]
        return combined[[c for c in cols if c in combined.columns]]
    return combined[["cluster", "gene", "p_val", "avg_log2FC", "pct.1", "pct.2", "p_val_adj"]]

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.01, logfc_threshold: float = 0.1, 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.01,
    logfc_threshold: float = 0.1,
    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