Single-cell analysis workflow: Seurat v5 and Scanpy from QC to cell annotation

Step 1

From Cell Ranger output to an object: which file to read and what changed from v4 to v5

Read filtered_feature_bc_matrix (or the .h5 of the same name); empty droplets have already been removed. raw_feature_bc_matrix is only needed as input to SoupX or CellBender.

# Seurat v5 / SeuratObject 5.4: replacements for v4 code
GetAssayData(obj, slot = "counts")   # Error: The `slot` argument of `GetAssayData()` was deprecated in SeuratObject 5.0.0 and is now defunct.
obj@assays$RNA@counts                # Error: no slot of name "counts" for this object of class "Assay5"

obj[["RNA"]]$counts                  # v5 accessors
LayerData(obj, assay = "RNA", layer = "data")
GetAssayData(obj, layer = "data")    # works only while the assay has one "data" layer

# Old v4 object read into v5: convert the assay class when a tool needs Assay5 (e.g. BPCells, IntegrateLayers)
obj[["RNA"]] <- as(obj[["RNA"]], "Assay5")
# Tool that only accepts v4 assays: convert back (merge layers first)
obj <- JoinLayers(obj); obj[["RNA"]] <- as(obj[["RNA"]], "Assay")

Step 2

QC thresholds: derive them from the data instead of copying 200, 2500 and 5%

QC removes three kinds of barcodes: those with almost no signal (empty droplets or broken cells), those with unusually many genes and UMIs (mostly doublets), and those with a high mitochondrial fraction (damaged membranes that lost cytoplasmic mRNA). Compute thresholds separately for each sample.

Tested on this page: macOS arm64 (8 cores, 16 GB), Scanpy 1.12.4 / anndata 0.13.4 / Python 3.12, Seurat 5.5.1 / SeuratObject 5.4.0 / R 4.5, 10x PBMC 3k (2,700 cells). Other jobs were running on the same machine during the test, so run times indicate order of magnitude only.

In PBMC 3k the median mitochondrial fraction is 2.03%, the 95th percentile 4.01% and the 99th percentile 5.88%. The tutorial's fixed rule (200 < genes < 2500 and mitochondrial < 5%) keeps 2,638 cells. The MAD rule (5 MAD on UMIs, genes and top-20 share; mitochondrial 3 MAD and > 8%) removes 104 cells and keeps 2,596; the MAD interval is 366 to 1,820 for genes and 708 to 6,811 for UMIs.

On the same data the 3 MAD upper bound for the mitochondrial fraction is only 3.65%. Without the 8% floor, filtering on 3 MAD alone raises the number of removed cells from 104 to 276; most of the extra cells are normal PBMCs with 3.7%–5% mitochondrial reads. This is why the MAD rule needs a fixed floor.

The MAD rule also removes real populations with low RNA content. PBMC 3k contains 19 platelet-like cells with high PPBP and a median of only 397 genes; the MAD rule removes 10 of them, the rest are too few to form a cluster, and the Scanpy workflow ends without a platelet cluster, whereas the fixed-threshold workflow keeps a 13-cell platelet cluster. If platelets, erythroid precursors, neutrophils or other low-RNA cells matter for your question, check the markers of the removed cells before fixing the thresholds.

With several samples, compute MAD per sample. Pooled MAD lets a deeply sequenced sample turn normal cells from a shallow sample into outliers.

Step 3

Doublets and ambient RNA: when they must be handled and how

Run doublet detection separately for each capture channel (each 10x channel) and before strict QC. Ambient RNA correction is optional; it matters most for marker and differential expression analysis.

Tested on this page: on the same PBMC 3k data scDblFinder (1.24.10) calls 124 doublets (4.6%) and Scanpy's built-in Scrublet calls only 32 (1.2%). At 0.8% per 1,000 cells, about 58 doublets are expected among 2,700 cells, so the two tools fall on either side of the expectation. Report which tool you used and how many cells it removed, and check on the UMAP whether the called doublets sit between two clusters.

Called doublets can be flagged rather than removed. After clustering, if a small cluster expresses markers of two lineages (for example CD3E and CD14) and has high doublet scores, remove the whole cluster. Single-cell Best Practices also keeps flagged doublets at first.

Janssen et al. 2023 (Genome Biology) compared CellBender, DecontX and SoupX on genotype-mixed samples: background removal improved clustering and cell classification only slightly, improved marker gene detection the most, and CellBender estimated the background level most accurately. For clustering and coarse annotation the step can be skipped; for differential expression or claims that a gene is specific to a cluster, correct for ambient RNA.

Step 4

Normalization, highly variable genes and the number of PCs

LogNormalize (normalize_total + log1p in Scanpy) is a safe default. Ahlmann-Eltze and Huber 2023 (Nature Methods) found that log(y/s + 1) followed by PCA performs as well as or better than SCTransform, Pearson residuals and other more complex transformations.

# Scanpy 1.12 (Python >= 3.12): PBMC 3k from Cell Ranger output to annotated clusters
import numpy as np, scanpy as sc
from scipy.stats import median_abs_deviation

adata = sc.read_10x_mtx("filtered_gene_bc_matrices/hg19", var_names="gene_symbols")
# Cell Ranger >= 3: sc.read_10x_h5("filtered_feature_bc_matrix.h5")
adata.var_names_make_unique()

# QC metrics: mitochondrial / ribosomal / hemoglobin genes (mouse: "mt-", ("Rps","Rpl"), "^Hb[^(p)]")
adata.var["mt"] = adata.var_names.str.startswith("MT-")
adata.var["ribo"] = adata.var_names.str.startswith(("RPS", "RPL"))
adata.var["hb"] = adata.var_names.str.contains("^HB[^(P)]")
sc.pp.calculate_qc_metrics(adata, qc_vars=["mt", "ribo", "hb"], percent_top=[20], log1p=True, inplace=True)

# Doublets on raw counts, per capture (batch_key="sample" when several samples are concatenated)
sc.pp.scrublet(adata, random_state=0)

# MAD outliers: 5 MAD on counts / genes / top-20 share, mito 3 MAD AND above a fixed floor
def outlier(x, nmads, upper_only=False):
    med, mad = np.median(x), median_abs_deviation(x)
    return (x > med + nmads * mad) if upper_only else ((x < med - nmads * mad) | (x > med + nmads * mad))
o = adata.obs
bad = (outlier(o.log1p_total_counts, 5) | outlier(o.log1p_n_genes_by_counts, 5)
       | outlier(o.pct_counts_in_top_20_genes, 5)
       | (outlier(o.pct_counts_mt, 3, upper_only=True) & (o.pct_counts_mt > 8)))
adata = adata[~bad & ~adata.obs.predicted_doublet].copy()
sc.pp.filter_genes(adata, min_cells=3)

# Normalize, HVG, PCA
adata.layers["counts"] = adata.X.copy()
sc.pp.normalize_total(adata, target_sum=1e4)
sc.pp.log1p(adata)
sc.pp.highly_variable_genes(adata, n_top_genes=2000, flavor="seurat_v3", layer="counts")  # add batch_key="sample" for multi-sample data
adata.raw = adata                                     # keep all genes (log-normalized) for markers and plots
hv = adata[:, adata.var.highly_variable].copy()       # scale only HVGs; scaling all genes makes X dense
sc.pp.scale(hv, max_value=10)
sc.tl.pca(hv, n_comps=50)

# Graph, clustering, UMAP
sc.pp.neighbors(hv, n_neighbors=15, n_pcs=30)
for r in (0.2, 0.4, 0.6, 0.8, 1.0):
    sc.tl.leiden(hv, resolution=r, flavor="igraph", n_iterations=2, key_added=f"leiden_{r}")
sc.tl.umap(hv)

# Markers on log-normalized expression of all genes
adata.obs["leiden"] = hv.obs["leiden_0.6"].values
sc.tl.rank_genes_groups(adata, "leiden", method="wilcoxon", use_raw=False)
sc.pl.rank_genes_groups_dotplot(adata, n_genes=5)

# Seurat v5: same workflow in R
library(Seurat)
counts <- Read10X("filtered_gene_bc_matrices/hg19")      # Cell Ranger >= 3: Read10X_h5("filtered_feature_bc_matrix.h5")
obj <- CreateSeuratObject(counts, project = "pbmc3k", min.cells = 3, min.features = 200)
obj[["percent.mt"]] <- PercentageFeatureSet(obj, pattern = "^MT-")   # mouse: "^mt-"

# Doublets: scDblFinder on the raw counts of each capture, before strict QC
library(scDblFinder); library(SingleCellExperiment)
set.seed(1)
sce <- scDblFinder(SingleCellExperiment(list(counts = LayerData(obj, layer = "counts"))))
obj$dbl_class <- sce$scDblFinder.class               # several captures: scDblFinder(sce, samples = "sample")

# MAD thresholds instead of the 200 / 2500 / 5% values of the PBMC3k tutorial
mad_out <- function(x, n, upper = FALSE) {
  m <- median(x); d <- mad(x)
  if (upper) x > m + n * d else x < m - n * d | x > m + n * d
}
bad <- mad_out(log1p(obj$nCount_RNA), 5) | mad_out(log1p(obj$nFeature_RNA), 5) |
       (mad_out(obj$percent.mt, 3, upper = TRUE) & obj$percent.mt > 8)
obj <- subset(obj, cells = colnames(obj)[!bad & obj$dbl_class == "singlet"])

obj <- NormalizeData(obj)                            # LogNormalize, scale.factor = 1e4
obj <- FindVariableFeatures(obj, nfeatures = 2000)
obj <- ScaleData(obj)                                # HVGs only by default
obj <- RunPCA(obj, npcs = 50)
ElbowPlot(obj, ndims = 50)
obj <- FindNeighbors(obj, dims = 1:20)
obj <- FindClusters(obj, resolution = c(0.2, 0.4, 0.6, 0.8, 1.0))
obj <- RunUMAP(obj, dims = 1:20)

library(clustree)                                    # attach it; clustree::clustree() alone fails at ggsave
clustree(obj, prefix = "RNA_snn_res.")
Idents(obj) <- "RNA_snn_res.0.6"
markers <- FindAllMarkers(obj, only.pos = TRUE, min.pct = 0.25, logfc.threshold = 0.25)

Step 5

Multi-sample integration: look at the unintegrated result first, then choose Harmony, RPCA or scVI

Start without integration and color the UMAP by sample. If cells of the same type from different samples already mix, integration is unnecessary. Integration only changes the embedding (a harmony, integrated.rpca or similar reduction); the expression matrix is unchanged, and differential expression always uses the original RNA counts and data.

Tested on this page (Seurat 5.5.1, harmony 2.0.5, PBMC 3k randomly split into two "samples"): running IntegrateLayers(method = HarmonyIntegration) without split fails with attempt to set an attribute on NULL; after split Harmony completes. After integration, subsetting a subcluster leaves the layers joined, and running IntegrateLayers again fails with the same error; it succeeds after splitting the subset again. Calling split on a merged object that already has one layer per sample fails with The following layers are already split. CSDN authors report that CCAIntegration without split fails with no applicable method for 'Assays' applied to an object of class "NULL".

Tested on this page (Scanpy 1.12.4, harmonypy 2.0.2): sc.external.pp.harmony_integrate fails with ValueError: Value passed for key 'X_pca_harmony' is of incorrect shape. harmonypy 2.x already returns Z_corr as cells × PCs, and the Scanpy wrapper transposes it again. The code below calls harmonypy directly.

# Seurat v5 multi-sample integration: split -> per-layer preprocessing -> IntegrateLayers -> JoinLayers
obj <- merge(s1, y = list(s2, s3), add.cell.ids = c("s1", "s2", "s3"))   # already one layer per object: counts.1, counts.2, counts.3
# Object built from one matrix, or after JoinLayers: split by sample first
# obj[["RNA"]] <- split(obj[["RNA"]], f = obj$sample)
# Splitting a merged object again fails: The following layers are already split: 'counts.1', 'counts.2' Please join before splitting
obj <- NormalizeData(obj); obj <- FindVariableFeatures(obj); obj <- ScaleData(obj); obj <- RunPCA(obj)

# Unintegrated baseline first
obj <- FindNeighbors(obj, dims = 1:30, reduction = "pca")
obj <- FindClusters(obj, resolution = 0.5, cluster.name = "unintegrated_clusters")
obj <- RunUMAP(obj, dims = 1:30, reduction = "pca", reduction.name = "umap.unintegrated")

obj <- IntegrateLayers(obj, method = HarmonyIntegration, orig.reduction = "pca", new.reduction = "harmony")
# alternatives: method = RPCAIntegration (more conservative), CCAIntegration, FastMNNIntegration, scVIIntegration
obj <- FindNeighbors(obj, reduction = "harmony", dims = 1:30)
obj <- FindClusters(obj, resolution = 0.5, cluster.name = "harmony_clusters")
obj <- RunUMAP(obj, reduction = "harmony", dims = 1:30, reduction.name = "umap.harmony")

obj <- JoinLayers(obj)          # before FindMarkers / FindAllMarkers / GetAssayData

# Re-clustering a subset: layers are joined now, so split again before IntegrateLayers
sub <- subset(obj, idents = c("0", "3"))
sub[["RNA"]] <- split(sub[["RNA"]], f = sub$sample)
sub <- NormalizeData(sub); sub <- FindVariableFeatures(sub); sub <- ScaleData(sub); sub <- RunPCA(sub)
sub <- IntegrateLayers(sub, method = HarmonyIntegration, orig.reduction = "pca", new.reduction = "harmony_sub")
sub <- FindNeighbors(sub, reduction = "harmony_sub", dims = 1:20); sub <- FindClusters(sub, resolution = 0.3)
sub <- JoinLayers(sub)

# Scanpy: Harmony (harmonypy) and scVI on the same AnnData
import scanpy as sc
sc.pp.highly_variable_genes(adata, n_top_genes=2000, flavor="seurat_v3", layer="counts", batch_key="sample")
hv = adata[:, adata.var.highly_variable].copy()
sc.pp.scale(hv, max_value=10); sc.tl.pca(hv, n_comps=50)
# sc.external.pp.harmony_integrate(hv, key="sample") fails with harmonypy 2.x:
#   ValueError: Value passed for key 'X_pca_harmony' is of incorrect shape ... (50,) while it should have had (n_cells,)
# harmonypy 2.x already returns Z_corr as cells x PCs, the Scanpy 1.12.4 wrapper transposes it again
import harmonypy
ho = harmonypy.run_harmony(hv.obsm["X_pca"], hv.obs, "sample")
Z = ho.Z_corr
hv.obsm["X_pca_harmony"] = Z if Z.shape[0] == hv.n_obs else Z.T
sc.pp.neighbors(hv, use_rep="X_pca_harmony", n_pcs=30)
sc.tl.leiden(hv, resolution=0.5, flavor="igraph", n_iterations=2)

# scVI: needs raw counts; a GPU shortens training, CPU works for tens of thousands of cells
import scvi
scvi.model.SCVI.setup_anndata(hv, layer="counts", batch_key="sample")
m = scvi.model.SCVI(hv, n_latent=30)
m.train()                                                      # max_epochs scales down automatically for large data
hv.obsm["X_scVI"] = m.get_latent_representation()
sc.pp.neighbors(hv, use_rep="X_scVI")

Step 6

Choosing the clustering resolution: clustree plus stability across seeds

Resolution sets how coarse or fine clusters are, and there is no single correct value. Scan a range of resolutions and pick the lowest one at which the cluster structure is stable and every cluster has interpretable markers; for finer subtypes, subset a large cluster and recluster it.

On PBMC 3k, at res 1.5 the number of clusters varies between 16 and 18 across seeds and ARI falls to 0.75, so the extra clusters are unstable. With 30 PCs and res 0.4 the five seeds give nearly identical results (ARI 0.99). At the same resolution, 10 PCs yield 1–3 more clusters than 30 PCs, so report the number of PCs together with the resolution.

In a clustree plot, a cluster receiving arrows from two clusters of the previous level means the resolution is high enough that cells are being reassigned. In this page's PBMC 3k clustree plot such crossings start between res 0.4 and 0.6, consistent with the 0.5 used in the Seurat tutorial. Seurat's rule of thumb is res 0.4–1.2 for about 3,000 cells, with higher values for more cells.

Step 7

Cell annotation: marker genes first, automatic tools as a cross-check

Automatic annotation tools can only assign labels present in their reference. Let an automatic tool propose a label for each cluster, check every cluster with a dot plot of canonical markers, and when they conflict, trust the markers and the tissue context.

Tested on this page (CellTypist 1.7.1, Immune_All_Low, majority_voting, over_clustering set to the 8 Leiden clusters at res 0.6): the voted labels agree with the marker-based calls for all 8 clusters: Tcm/Naive helper T (1,165 cells), Tem/Trm cytotoxic T (280), B cells (332), Classical monocytes (465), CD16+ NK (157), Non-classical monocytes (168), DC (27) and Megakaryocytes/platelets (13). Per-cell predictions produce 41 different labels, so reading only per-cell results splits one cluster into many small subtypes.

At this resolution naive and memory CD4 T cells are not separated (both are in the 1,165-cell cluster), and CCR7 and LEF1 are high only in part of it. To separate them, recluster this cluster on its own, or raise the resolution and check with CCR7 and S100A4.

# CellTypist on log1p data normalized to 10,000 counts per cell (adata after normalize_total + log1p)
import celltypist
from celltypist import models
models.download_models(model=["Immune_All_Low.pkl"])
pred = celltypist.annotate(adata, model="Immune_All_Low.pkl",
                           majority_voting=True, over_clustering=adata.obs["leiden"])
adata.obs["celltypist"] = pred.predicted_labels["majority_voting"]
adata.obs["celltypist_conf"] = pred.probability_matrix.max(axis=1)   # low values = label not in the model

# SingleR + celldex (Bioconductor). On macOS arm64 bioconda has no SingleR build:
# install with BiocManager::install(c("SingleR", "celldex")) inside the R environment
library(SingleR); library(celldex)
ref <- celldex::MonacoImmuneData()                    # immune; HumanPrimaryCellAtlasData() for broad tissue types
pred <- SingleR(test = LayerData(obj, layer = "data"), ref = ref,
                labels = ref$label.fine, clusters = obj$seurat_clusters)   # per-cluster mode
obj$singler <- pred$pruned.labels[match(obj$seurat_clusters, rownames(pred))]
plotScoreHeatmap(pred)                                # check that the best score clearly beats the second

Two routes

Seurat and Scanpy side by side, and converting objects

Both routes can do the same analysis. Above about 100,000 cells, or with deep learning tools such as scVI, the Python route uses less memory; choose Seurat when you need R packages such as SingleR, CellChat or Monocle.

For conversion use anndataR from scverse (in Bioconductor; latest GitHub release v1.3.2 on 2026-10-07; this page tested version 1.2.2 from bioconda). The SeuratDisk repository has had no commits since November 2023 and has 156 open issues; common problems include Convert filling X from scale.data first so raw counts are not transferred, and failures reading LZF-compressed h5ad files.

Tested on this page (anndataR 1.2.2): write_h5ad on a v5 object without JoinLayers does not fail; it only warns Skipping Layer "counts.A" with unexpected dimensions, and the resulting h5ad contains no expression matrix. After JoinLayers the file contains counts and data layers with X empty, so adata.X has to be set in Python. Reading an h5ad written by Scanpy turns X and the layers into Seurat layers and keeps the obs columns as well as X_pca and X_umap. anndataR does not convert varp or Seurat's Neighbors and Images. zellkonverter, which goes through SingleCellExperiment, is another workable route.

# R: anndataR 1.2 (Bioconductor). rhdf5 is a separate install, otherwise:
#   Error: HDF5AnnData requires the rhdf5 package
# BiocManager::install(c("anndataR", "rhdf5"))
library(Seurat); library(anndataR)

obj <- read_h5ad("pbmc3k_processed.h5ad", as = "Seurat")   # h5ad -> Seurat: X and layers become layers, obsm -> reductions
obj <- JoinLayers(obj)                                       # split layers are skipped with a warning only
write_h5ad(obj, "from_seurat.h5ad")                          # Seurat -> h5ad: layers "counts" and "data", X left empty

# Python: X is empty after the Seurat -> h5ad direction, point it at a layer
import anndata as ad
adata = ad.read_h5ad("from_seurat.h5ad")
adata.X = adata.layers["counts"].copy()                      # or layers["data"] for log-normalized values

Resources

Memory and run time: estimating from cell numbers, and what to do when memory runs out

The memory bottleneck in single-cell analysis is dense matrices. Count matrices are sparse; once ScaleData or sc.pp.scale is applied to all genes, it creates a dense matrix of cells × genes.

Tested on this page: the filtered PBMC 3k sparse matrix takes 16.8 MB; the same matrix stored densely (2,607 × 13,607 × 8 bytes) takes 270.6 MB, 16 times more. At that ratio, scaling all genes for 100,000 cells × 20,000 genes needs about 16 GB for that one matrix alone.

When memory runs out, try these in order: run ScaleData on HVGs only (the Seurat default; do not pass features = rownames(obj)); delete scale.data and SCT assays you no longer need before saving; in Seurat keep counts on disk with BPCells, analyze a 50,000-cell sketch with SketchData and project back to all cells with ProjectData; in Python use anndata's read_lazy or backed mode and load only the subset you need.

If R reports vector memory limit of 16.0 Gb reached on macOS, set R_MAX_VSIZE=32Gb in ~/.Renviron. Multithreaded Seurat also hits the future.globals.maxSize limit; for large datasets set options(future.globals.maxSize = 4e9) or higher as in the official vignette.

# Seurat v5 + BPCells: counts stay on disk, analyze a 50,000-cell sketch, project back
library(Seurat); library(BPCells)
options(future.globals.maxSize = 4e9)
mat <- open_matrix_10x_hdf5("filtered_feature_bc_matrix.h5")
write_matrix_dir(mat, dir = "bpcells_counts")              # one-time conversion to the on-disk format
obj <- CreateSeuratObject(open_matrix_dir("bpcells_counts"))
obj <- NormalizeData(obj); obj <- FindVariableFeatures(obj)
obj <- SketchData(obj, ncells = 50000, method = "LeverageScore", sketched.assay = "sketch")
DefaultAssay(obj) <- "sketch"
obj <- FindVariableFeatures(obj); obj <- ScaleData(obj); obj <- RunPCA(obj)
obj <- FindNeighbors(obj, dims = 1:50); obj <- FindClusters(obj, resolution = 1)
obj <- RunUMAP(obj, dims = 1:50, return.model = TRUE)
obj <- ProjectData(obj, assay = "RNA", full.reduction = "pca.full", sketched.assay = "sketch",
                   sketched.reduction = "pca", umap.model = "umap", dims = 1:50,
                   refdata = list(cluster_full = "seurat_clusters"))

# Python side: open .h5ad without loading X (anndata >= 0.12)
import anndata as ad
a = ad.experimental.read_lazy("big.h5ad")   # lazy obs/var/X; subset first, then .to_memory()
a = ad.read_h5ad("big.h5ad", backed="r")    # older backed mode: read-only X on disk

Tested on this page

Full results on PBMC 3k

Tested on this page: macOS arm64 (8 cores, 16 GB), Scanpy 1.12.4 / anndata 0.13.4 / Python 3.12, Seurat 5.5.1 / SeuratObject 5.4.0 / R 4.5, 10x PBMC 3k (2,700 cells). Other jobs were running on the same machine during the test, so run times indicate order of magnitude only.

• Errors reproduced on Seurat 5.5.1: GetAssayData(slot = ) is defunct; @counts does not exist; GetAssayData and FindMarkers fail with multiple layers; IntegrateLayers without split fails with attempt to set an attribute on NULL; re-integrating a subcluster fails with the same error.
• RunUMAP in Seurat 5.5.1 uses R's uwot with the cosine metric by default, so results will not match old tutorials that called Python umap-learn through reticulate.
• SingleR was not run on this machine: bioconda has no bioconductor-singler build for macOS arm64; only function and argument names were checked.
• CellBender and scVI need a GPU or long CPU time, and BPCells was not installed on this machine; for these three only argument names and documentation were checked.
• anndataR 1.2.2 was run in both directions: rhdf5 has to be installed separately, and h5ad files written without JoinLayers contain no expression matrix.

Chinese community

Verified experience from the Chinese community

The items below come from CSDN and Zhihu columns. Only posts with original error messages or the author's own measurements, consistent with official documentation or this page's tests, are included.

Hand it to an agent

Hand this workflow to an agent

You can describe the analysis in one sentence and let Scientify's research agent run it in a cloud computer.

Example instruction: "Analyze the six 10x samples in data/ (3 control, 3 treated): QC each sample with the MAD rule and run scDblFinder, show the unintegrated UMAP, integrate with Harmony and scVI and compare them; choose the resolution with clustree and multi-seed ARI; annotate with CellTypist and marker dot plots, and output cell-type proportions and treated-versus-control differential genes for every cluster."

In a cloud computer with Scanpy preinstalled, the agent installs the other packages it needs, runs QC and doublet detection per sample and records how many cells each step removes; after integration, clustering and annotation it outputs UMAPs, the clustree plot, marker dot plots, a proportion table and differential gene tables. When scVI needs a GPU it rents one for the task. The workspace keeps scripts, parameters, logs and h5ad files, so the run can be reproduced; the agent reviews its results adversarially, for example checking whether a treatment-specific cluster disappeared after integration and whether each annotated cluster has consistent markers. The task keeps running after you shut down your computer.

You still need to check: whether the QC thresholds suit your tissue, whether integration erased the biological differences you care about, whether annotation labels agree with the tissue context and the literature, and whether differential analysis uses samples, not cells, as replicates.