PBMC 3k Tutorial — R Seurat vs Truecell (Python)¶

A step-by-step translation of the official Seurat PBMC 3k tutorial into Python using the Truecell package. Every R code block is paired with the equivalent Python code and both plots are shown side by side so R users can follow along directly.

Dataset: 3k PBMCs from a Healthy Donor — 10x Genomics (2016)
R reference: Seurat v5 · Hao et al. 2024
Python: Truecell

Setup¶

R (Seurat) Python (Truecell)

library(Seurat)
library(dplyr)
library(patchwork)   # combine plots
library(ggplot2)     # theming

# Data structures + analysis
from truecell.datasets import pbmc3k
from truecell.truecell import create_truecell_object
from truecell.preprocessing import (
    normalize_data, find_variable_features,
    scale_data, percentage_feature_set,
)
from truecell.reduction import run_pca
from truecell.neighbors import find_neighbors
from truecell.clustering import find_clusters
from truecell.umap import run_umap
from truecell.markers import find_markers, find_all_markers

# Plotting (mirrors Seurat's plotting API)
from truecell.plotting import (
    dim_plot, feature_plot, vln_plot,
    elbow_plot, feature_scatter,
    variable_feature_plot, dim_heatmap,
    do_heatmap, ridge_plot,
)

import matplotlib.pyplot as plt   # save / display figures

Key difference: R plots render to the graphics device automatically. Truecell functions return a matplotlib.Figure — display it in a Jupyter notebook (it shows inline), or save it with fig.savefig("out.png").

Step 1 · Load Data¶

R (Seurat) Python (Truecell)

# Read 10x sparse matrix from disk
pbmc.data <- Read10X(
  data.dir = "pbmc3k/filtered_gene_bc_matrices/hg19/"
)
# 32,738 genes × 2,700 cells

# Downloads automatically to ~/.truecell_data/pbmc3k
# (or pass data_dir= to point at an existing folder)
counts, genes, cells = pbmc3k()
# counts : scipy.sparse.csc_matrix  (32,738 × 2,700)
# genes  : list[str]  – feature names
# cells  : list[str]  – barcode names

Step 2 · Create Object¶

R (Seurat) Python (Truecell)

pbmc <- CreateSeuratObject(
  counts     = pbmc.data,
  project    = "pbmc3k",
  min.cells  = 3,     # keep genes in ≥3 cells
  min.features = 200  # keep cells with ≥200 genes
)
# 13,714 features × 2,700 cells

pbmc = create_truecell_object(
    counts       = counts,
    project      = "pbmc3k",
    min_cells    = 3,    # keep genes in ≥3 cells
    min_features = 200,  # keep cells with ≥200 genes
    feature_names = genes,
    cell_names    = cells,
)
# 13,714 features × 2,700 cells
print(pbmc)
# Truecell object — pbmc3k
#   2700 cells × 13714 features
#   Active assay: 'RNA'

Naming: CreateSeuratObject → create_truecell_object. Arguments use _ instead of . (min.cells → min_cells). The returned object exposes the same slots: pbmc.meta_data, pbmc.assays, etc.

Step 3 · QC Metrics & Violin Plot¶

R (Seurat) Python (Truecell)

# Percent mitochondrial reads
pbmc[["percent.mt"]] <- PercentageFeatureSet(
  pbmc, pattern = "^MT-"
)

# Violin plot of three QC metrics
VlnPlot(
  pbmc,
  features = c("nFeature_RNA",
               "nCount_RNA",
               "percent.mt"),
  ncol = 3
)

# Percent mitochondrial reads
percentage_feature_set(
    pbmc, pattern=r"^MT-", col_name="percent.mt"
)

# Violin plot of three QC metrics
fig = vln_plot(
    pbmc,
    features=["nFeature_RNA", "nCount_RNA", "percent.mt"],
    ncol=3,
    figsize=(12, 4),
)
fig.savefig("qc_violin.png", dpi=150, bbox_inches="tight")

vln_plot accepts both gene names and metadata column names as features. The group_by argument (default: active idents) controls the x-axis grouping — same as R's group.by.

Step 4 · QC Scatter Plots¶

R (Seurat) Python (Truecell)

plot1 <- FeatureScatter(
  pbmc,
  feature1 = "nCount_RNA",
  feature2 = "percent.mt"
)
plot2 <- FeatureScatter(
  pbmc,
  feature1 = "nCount_RNA",
  feature2 = "nFeature_RNA"
)
plot1 + plot2   # patchwork combines them

fig1 = feature_scatter(
    pbmc, "nCount_RNA", "percent.mt"
)
fig2 = feature_scatter(
    pbmc, "nCount_RNA", "nFeature_RNA"
)
# Save individually or combine with matplotlib subplots

R uses patchwork's + operator to combine plots. In Python, pass figsize or use matplotlib.gridspec to arrange figures manually.

Step 5 · Filter Cells¶

R (Seurat) Python (Truecell)

pbmc <- subset(
  pbmc,
  subset = nFeature_RNA > 200 &
           nFeature_RNA < 2500 &
           percent.mt < 5
)
# 2,638 cells retained

md   = pbmc.meta_data
keep = (
    (md["nFeature_RNA"] > 200) &
    (md["nFeature_RNA"] < 2500) &
    (md["percent.mt"]   < 5)
)
pbmc = pbmc.subset(cells=list(md.index[keep]))
# 2,638 cells retained

subset() works in both languages. Python uses boolean pandas indexing; the result is identical (2,638 cells).

Step 6 · Normalise Data¶

R (Seurat) Python (Truecell)

pbmc <- NormalizeData(
  pbmc,
  normalization.method = "LogNormalize",
  scale.factor = 10000
)

normalize_data(
    pbmc,
    normalization_method = "LogNormalize",
    scale_factor         = 10000,
)
# Stored in pbmc.assays["RNA"].layers["data"]

Both apply log1p(counts / total_counts × 10,000). Result is stored in the data layer of the active assay.

Step 7 · Highly Variable Features¶

R (Seurat) Python (Truecell)

pbmc <- FindVariableFeatures(
  pbmc,
  selection.method = "vst",
  nfeatures = 2000
)

# Visualise
top10 <- head(VariableFeatures(pbmc), 10)
plot1 <- VariableFeaturePlot(pbmc)
plot2 <- LabelPoints(
  plot = plot1, points = top10, repel = TRUE
)
plot1 + plot2

find_variable_features(
    pbmc,
    selection_method = "vst",
    nfeatures        = 2000,
)

# Visualise (labels top 10 automatically)
fig = variable_feature_plot(
    pbmc, label=True, n_label=10, repel=True
)
# Access the list:
hvg = pbmc.assays["RNA"].variable_features
top10 = hvg[:10]
print(top10)

Both plots use Standardized Variance on the y-axis, matching R's VariableFeaturePlot. Truecell reproduces Seurat's vst algorithm faithfully: it fits the mean–variance LOESS on the raw counts, then ranks genes by the variance of the standardized values after clipping each to sqrt(n_cells) (the clip step that stops single-cell outliers from dominating). Both move the ten names apart: R's LabelPoints(repel = TRUE) with ggrepel, which starts from random jitter, and Truecell's repel=True with a deterministic placement worked out whenever the figure is drawn. Without repel, names in the crowded top-right corner print over one another. Top-10 gene overlap: 9/10 (90%) — the same PPBP, LYZ, S100A9, IGLL5, GNLY, FTL, PF4, FTH1, S100A8 HVGs as the R tutorial; GNG11 sits just outside the top 10 (rank 11). Minor rank differences come from the LOESS implementation (R: Fortran; Python: local quadratic fit).

Step 8 · Scale Data¶

R (Seurat) Python (Truecell)

all.genes <- rownames(pbmc)
pbmc <- ScaleData(pbmc, features = all.genes)
# Stores z-scores in pbmc[["RNA"]]@scale.data

all_genes = pbmc.assays["RNA"]._all_feature_names
scale_data(pbmc, features=all_genes)
# Stored in pbmc.assays["RNA"].layers["scale.data"]

Step 9 · PCA¶

R (Seurat) Python (Truecell)

pbmc <- RunPCA(
  pbmc,
  features = VariableFeatures(pbmc),
  npcs     = 50
)

# Top positive/negative loading genes per PC
VizDimLoadings(pbmc, dims = 1:2, reduction = "pca")

# PCA scatter
DimPlot(pbmc, reduction = "pca")

hvg = pbmc.assays["RNA"].variable_features
run_pca(pbmc, n_pcs=50, features=hvg,
        reduction_name="pca")

# Top positive/negative loading genes per PC
# viz_dim_loadings mirrors VizDimLoadings exactly
fig = viz_dim_loadings(
    pbmc, reduction="pca", dims=[1, 2], n_features=15
)

# PCA scatter
fig = dim_plot(pbmc, reduction="pca", label=False)

Both show horizontal bar charts of top loading genes per PC. The same myeloid genes dominate PC1 positively (CST3, TYROBP, LST1, AIF1) and T-cell genes dominate PC1 negatively (IL7R, LTB, IL32).

Step 10 · Elbow Plot — Choose Dimensionality¶

R (Seurat)	Python (Truecell)
`ElbowPlot(pbmc) # Elbow visible around PC 9–10`	`fig = elbow_plot(pbmc, ndims=20) # Elbow visible around PC 9–10`

PC	R stdev	Truecell stdev
1	~6.8	6.766
2	~4.8	4.808
10	~1.7	1.684

Step 11 · Find Neighbors & Cluster¶

R (Seurat) Python (Truecell)

pbmc <- FindNeighbors(pbmc, dims = 1:10)
pbmc <- FindClusters(
  pbmc,
  resolution = 0.5   # → 9 clusters
)
# pbmc$seurat_clusters holds the labels

find_neighbors(pbmc, dims=range(10), k_param=20)
find_clusters(
    pbmc,
    resolution  = 0.5,   # → 9 clusters
    algorithm   = 1,     # 1 = Louvain, 4 = Leiden
    random_seed = 0,
)
# pbmc.meta_data["seurat_clusters"] holds labels

Both run Seurat's modularity optimiser: truecell's is a translation of Seurat's C++, with the same ten restarts and the same random stream. At resolution = 0.5 each returns 9 clusters, numbered alike, and the two runs agree about 2,562 of the 2,638 cells (ARI 0.928). Given the same graph they would agree on every cell, so what is left comes from the graph, traced in the note below. The largest difference is the CD8 T cluster, 303 cells here and 270 in Seurat.

The whole sweep, not one point on it¶

Picking a resolution means running a few and comparing them, so a fidelity claim pinned to a single setting says less than it looks like it does. The tutorial runs find_clusters(pbmc, resolution=[0.4, 0.8, 1.2, 0.5]) — Seurat's own vector idiom — and scores every one against R:

resolution	truecell	Seurat	ARI	concordance
0.4	9	9	0.9185	0.9701
0.5	9	9	0.9284	0.9712
0.8	11	11	0.9434	0.9746
1.2	11	12	0.9253	0.9632

The cluster count matches at 0.4, 0.5 and 0.8, and Seurat finds one more at 1.2. Agreement stays between 0.92 and 0.94 across the sweep. It used to fall from 0.90 to 0.80 as the resolution rose, while the two tools ran different optimisers.

0.5 is given last on purpose. Seurat leaves the object on the last resolution in the sequence, so every step below — UMAP, markers, annotation, the handoff — sees exactly the partition it saw before this sweep existed. The four outputs this tutorial writes are byte-identical to what it produced without it.

Each ARI now carries a declared band (CLUSTER_BANDS), for a reason: this file documented ARI 0.938 at resolution 0.5 while measuring 0.899, across six documents. The number had drifted — most likely with the graph fixes in #67–#71, which moved cells between clusters — and with no band on it, nothing failed. That is the same drift the DE tutorial's deseq2 top50 band did catch at the time (25 → 22). The bands were re-measured when find_clusters became Seurat's optimiser, which raised every ARI in the table above.

Where it comes from is traceable: the two runs keep the same 2,638 barcodes and agree on the per-gene VST means to 4.8e-14 and observed variances to 5.1e-14, but variance.standardized depends on a degree-2 LOESS fit that differs by up to 2.6e-2 relative — and variance.expected, the fit itself, differs by 2.5e-2, so the disagreement is the fit and only the fit. That flips 2 of the 2,000 variable features (both at ranks 1,982–2,000, and 0.03 apart in the fit), which moves the PCA slightly — matched |r| 0.9986 over the 10 dims clustering uses — which leaves the two SNN graphs 328 entries apart (194,226 against 194,554), which moves the cells the two partitions disagree about.

Run Rscript tutorials/pbmc3k_verify.R then python tutorials/pbmc3k_tutorial.py --report to reproduce every number in this note.

Step 12 · UMAP¶

R (Seurat) Python (Truecell)

pbmc <- RunUMAP(pbmc, dims = 1:10)
DimPlot(pbmc, reduction = "umap")

run_umap(pbmc, dims=range(10),
         reduction_name="umap", seed=42)
# label=False matches R's default (no centroid labels)
fig = dim_plot(pbmc, reduction="umap", label=False)

Note — UMAP layout difference is expected, not a bug. R Seurat uses the uwot UMAP implementation while Python uses umap-learn. Both libraries apply the same UMAP algorithm but differ in initialisation strategy, nearest-neighbour graph construction, and optimisation details. As a result, the spatial arrangement of clusters will look different even with the same random seed — clusters may be rotated, reflected, or positioned further apart. This is a known and well-documented difference between the two libraries and does not indicate an error in the analysis. The important thing is that the same 9 biologically meaningful clusters are recovered in both, numbered alike: both tools number clusters by size.

Step 13 · Feature Plots — Canonical Markers¶

R (Seurat) Python (Truecell)

FeaturePlot(
  pbmc,
  features = c(
    "MS4A1", "CD79A",
    "NKG7",  "GNLY",
    "FCGR3A","LYZ",
    "PPBP",  "CD8A",
    "IL7R"
  )
)

fig = feature_plot(
    pbmc,
    features=[
        "MS4A1", "CD79A",
        "NKG7",  "GNLY",
        "FCGR3A","LYZ",
        "PPBP",  "CD8A",
        "IL7R",
    ],
    reduction = "umap",
    ncol      = 3,
    colormap  = "YlOrRd",
)

Note — minor visual differences are expected, not bugs. Both plots use the same yellow→red expression gradient and the order=True behaviour (high-expression cells drawn on top). Cells with zero expression are rendered in light gray in both R and Truecell, making the expressing cells stand out clearly. The spatial layout of the underlying UMAP differs for the same reason described in Step 12 (different UMAP implementations), but the gene expression patterns are biologically identical: MS4A1/CD79A in B cells, NKG7/GNLY in NK cells, LYZ/FCGR3A in monocytes, PPBP in platelets, CD8A in CD8 T cells, and IL7R in T cells.

Step 14 · Violin Plots — Marker Expression¶

R (Seurat) Python (Truecell)

VlnPlot(
  pbmc,
  features = c("MS4A1", "CD79A")
)
VlnPlot(
  pbmc,
  features = c("NKG7", "PF4"),
  slot = "counts",
  log  = TRUE
)

# Matches R's first VlnPlot call (MS4A1 + CD79A)
fig = vln_plot(
    pbmc,
    features=["MS4A1", "CD79A"],
    figsize=(12, 4),
)

# Matches R's second VlnPlot call (raw counts)
fig = vln_plot(
    pbmc,
    features=["NKG7", "PF4"],
    layer="counts",   # slot="counts" in R
    figsize=(12, 4),
)

Both show cluster-specific marker expression. The layer argument in Truecell maps directly to Seurat's slot argument.

Step 15 · Find Markers¶

R (Seurat) Python (Truecell)

# Markers for one cluster
cluster2.markers <- FindMarkers(
  pbmc, ident.1 = 2
)
head(cluster2.markers, 5)

# Markers for all clusters
pbmc.markers <- FindAllMarkers(
  pbmc,
  only.pos        = TRUE,
  min.pct         = 0.25,
  logfc.threshold = 0.25
)
pbmc.markers %>%
  group_by(cluster) %>%
  slice_max(n = 2, order_by = avg_log2FC)

# Markers for one cluster
c2_markers = find_markers(
    pbmc, ident_1="2", only_pos=True
)
print(c2_markers.head(5))

# Markers for all clusters
all_markers = find_all_markers(
    pbmc,
    only_pos        = True,
    min_pct         = 0.25,
    logfc_threshold = 0.25,
)
top2 = (
    all_markers
    .groupby("cluster", group_keys=False)
    .apply(lambda x: x.nlargest(2, "avg_log2FC"))
)
print(top2)

Argument names use _ instead of . (ident.1 → ident_1, only.pos → only_pos). The returned DataFrame has the same columns: p_val, avg_log2FC, pct.1, pct.2, p_val_adj, cluster, gene.

Step 16 · Expression Heatmap¶

R (Seurat) Python (Truecell)

top10 <- pbmc.markers %>%
  group_by(cluster) %>%
  top_n(n = 10, wt = avg_log2FC)

DoHeatmap(pbmc, features = top10$gene) +
  NoLegend()

top10_genes = (
    all_markers
    .groupby("cluster", group_keys=False)
    .apply(lambda x: x.nlargest(10, "avg_log2FC"))
    ["gene"].tolist()
)
top10_genes = list(dict.fromkeys(top10_genes))

fig = do_heatmap(pbmc, features=top10_genes)

do_heatmap automatically sorts cells by cluster and draws a coloured cluster bar at the top — matching R's DoHeatmap layout. NoLegend() in R is handled by do_heatmap's internal layout.

Step 17 · Cell Type Annotation¶

R (Seurat) Python (Truecell)

new.cluster.ids <- c(
  "Naive CD4 T", "CD14+ Mono",
  "Memory CD4 T", "B",
  "CD8 T", "FCGR3A+ Mono",
  "NK", "DC", "Platelet"
)
names(new.cluster.ids) <- levels(pbmc)
pbmc <- RenameIdents(pbmc, new.cluster.ids)

DimPlot(
  pbmc,
  reduction = "umap",
  label     = TRUE,
  pt.size   = 0.5
) + NoLegend()

cell_type_map = {
    "0": "Naive CD4 T",
    "1": "CD14+ Mono",
    "2": "Memory CD4 T",
    "3": "B",
    "4": "CD8 T",
    "5": "FCGR3A+ Mono",
    "6": "NK",
    "7": "DC",
    "8": "Platelet",
}
pbmc.rename_idents(cell_type_map)

fig = dim_plot(
    pbmc,
    reduction = "umap",
    label     = True,
    repel     = True,
    pt_size   = 0.5,
    title     = "UMAP — Cell Type Annotations",
)

Both tools number clusters by size, and here the nine come out in the same order on both sides, so the map is Seurat's own, key for key. The two partitions agree at ARI 0.928 — see Step 11.

Both put each label at the median of its cluster's cells. repel=True then moves a label only where it would overlap another, run off the panel, or cover more than half of a cluster's cells. Here that moves the labels of the small clusters off the cells they name: at its median, the "Platelet" label's box covered all 14 platelets.

Step 18 · Ridge Plot (Bonus)¶

R (Seurat)	Python (Truecell)
`RidgePlot( pbmc, features = c("LYZ", "NKG7", "MS4A1", "CD8A"), ncol = 2 )`	`fig = ridge_plot( pbmc, features=["LYZ", "NKG7", "MS4A1", "CD8A"], ncol=2, figsize=(12, 8), )`

The published Seurat vignette prints no RidgePlot, so the R panel here is generated by pbmc3k_verify.R (standard PBMC 3k workflow). R's RidgePlot uses the ggridges package; Truecell implements equivalent KDE ridgelines using scipy.stats.gaussian_kde.

Quick Reference — API Translation¶

Task	R (Seurat)	Python (Truecell)
Create object	`CreateSeuratObject(counts, min.cells, min.features)`	`create_truecell_object(counts, min_cells, min_features)`
% mito genes	`PercentageFeatureSet(pbmc, pattern="^MT-")`	`percentage_feature_set(pbmc, pattern=r"^MT-", col_name=...)`
Normalise	`NormalizeData(pbmc, normalization.method, scale.factor)`	`normalize_data(pbmc, normalization_method, scale_factor)`
HVGs	`FindVariableFeatures(pbmc, selection.method, nfeatures)`	`find_variable_features(pbmc, selection_method, nfeatures)`
Scale	`ScaleData(pbmc, features)`	`scale_data(pbmc, features)`
PCA	`RunPCA(pbmc, features, npcs)`	`run_pca(pbmc, features, n_pcs)`
Neighbors	`FindNeighbors(pbmc, dims)`	`find_neighbors(pbmc, dims, k_param)`
Cluster	`FindClusters(pbmc, resolution)`	`find_clusters(pbmc, resolution, algorithm)`
UMAP	`RunUMAP(pbmc, dims)`	`run_umap(pbmc, dims)`
Markers	`FindMarkers(pbmc, ident.1)`	`find_markers(pbmc, ident_1)`
All markers	`FindAllMarkers(pbmc, only.pos, logfc.threshold)`	`find_all_markers(pbmc, only_pos, logfc_threshold)`
Rename idents	`RenameIdents(pbmc, new.ids)`	`pbmc.rename_idents(mapping_dict)`
Subset	`subset(pbmc, subset = condition)`	`pbmc.subset(cells=keep_list)`
Get expression	`FetchData(pbmc, vars)`	`pbmc.fetch_data(vars)`
Access metadata	`pbmc@meta.data`	`pbmc.meta_data`
Access assay	`pbmc[["RNA"]]`	`pbmc.assays["RNA"]`
Active idents	`Idents(pbmc)`	`pbmc.idents`

Plotting API Translation¶

R (Seurat)	Python (Truecell)	Key argument changes
`VlnPlot(pbmc, features, ncol, slot)`	`vln_plot(pbmc, features, ncol, layer)`	`slot` → `layer`
`FeaturePlot(pbmc, features, order)`	`feature_plot(pbmc, features, order)`	same
`DimPlot(pbmc, reduction, label, pt.size)`	`dim_plot(pbmc, reduction, label, pt_size)`	`.` → `_`
`ElbowPlot(pbmc)`	`elbow_plot(pbmc)`	same
`FeatureScatter(pbmc, feature1, feature2)`	`feature_scatter(pbmc, feature1, feature2)`	same
`VariableFeaturePlot(pbmc)`	`variable_feature_plot(pbmc)`	same
`VizDimLoadings(pbmc, dims, reduction)`	`viz_dim_loadings(pbmc, dims, reduction)`	same
`DimHeatmap(pbmc, dims, cells, balanced)`	`dim_heatmap(pbmc, dims, cells, balanced)`	same
`DoHeatmap(pbmc, features)`	`do_heatmap(pbmc, features)`	same
`RidgePlot(pbmc, features, ncol)`	`ridge_plot(pbmc, features, ncol)`	same

Plot output difference¶

# R — plots render to graphics device automatically
VlnPlot(pbmc, features = "LYZ")

# Python — functions return a Figure; display or save explicitly
fig = vln_plot(pbmc, features="LYZ")

# In a Jupyter notebook: just call the function — it displays inline
# To save:
fig.savefig("lyz_violin.png", dpi=150, bbox_inches="tight")
# To display interactively:
plt.show()

Reproducing All Plots¶

git clone https://github.com/GenomicAI/truecell.git
cd truecell
uv venv && source .venv/bin/activate   # Windows: .venv\Scripts\activate
uv pip install -e ".[analysis]"
python tutorials/generate_plots.py     # Truecell plots -> tutorials/figures/
Rscript tutorials/pbmc3k_verify.R      # R Seurat RidgePlot + the numeric handoff

Every other R panel in this tutorial links the canonical Seurat vignette image; only the RidgePlot (which the vignette omits) is generated locally by pbmc3k_verify.R.

Checking the numbers, not the pictures¶

Because the R panels are canonical images, every step above is compared by eye. pbmc3k_verify.R also writes a numeric handoff — per-cell QC and PCA keyed by barcode, per-gene VST statistics, the marker table, and a set of scalar anchors — which --report compares:

python tutorials/pbmc3k_tutorial.py    # writes the Python side
Rscript tutorials/pbmc3k_verify.R      # writes the R side
python tutorials/pbmc3k_tutorial.py --report

Neither side is pinned to the other: both run their own pipeline from the same 10x bytes, which is what makes the comparison worth anything. What it reports on this dataset:

	truecell vs Seurat 5.5.1
Cells surviving QC	the same 2,638 barcodes; nCount and nFeature exact, percent.mt to 5.3e-15
VST per gene (13,714)	mean 4.8e-14 · variance 1.6e-11 · `variance.expected` 2.5e-2 · `variance.standardized` 2.6e-2, all relative
Variable features	1,998 of 2,000 shared, rank Spearman 0.9999
PCA (10 dims)	matched \|r\| mean 0.9986, min 0.9905, no reordering
kNN graph	52,760 on both (2,638 × 20)
Clusters	9 vs 9 — ARI 0.928, concordance 0.971 (see Step 11)
Markers, where the clusters hold identical cells	identical gene sets (151/151, 927/927 and 242/242), `avg_log2FC` to 4.9e-15, 5.3e-15 and 4.6e-14 respectively

That last row is the one to read carefully. Three clusters — B cells, dendritic cells and platelets — came out with exactly the same membership on both sides, and on those three the differential expression agrees gene for gene and to machine precision. So the marker differences elsewhere are downstream of the clustering, not of the DE code.

Two notes on reading the report's marker table:

top10 by p is not the right column to judge. R's Wilcoxon p-values underflow to exactly 0 for the strongest markers — between 40 and 302 genes per cluster on this dataset — so Seurat's own top-10 ordering is decided by its fold-change tie-break, not by p. top10 by FC, which both tools resolve exactly, runs 8–10 of 10.
Building this handoff found a real defect in find_all_markers: it was missing Seurat's return.thresh = 0.01, so it returned rows Seurat does not. It now applies the same filter and breaks p-value ties on descending avg_log2FC the way FindAllMarkers does.

References¶

Hao Y, Stuart T, Kowalski MH, et al. (2024). Dictionary learning for integrative, multimodal and scalable single-cell analysis. Nature Biotechnology, 42, 293–304. https://doi.org/10.1038/s41587-023-01767-y

Stuart T, Butler A, Hoffman P, et al. (2019). Comprehensive Integration of Single-Cell Data. Cell, 177(7), 1888–1902. https://doi.org/10.1016/j.cell.2019.05.031

10x Genomics (2016). 3k PBMCs from a Healthy Donor. https://www.10xgenomics.com/resources/datasets/3-k-pb-mcs-from-a-healthy-donor-1-standard-1-1-0