Orchestrates the end-to-end single-cell RNA-seq pipeline from 10x Cell Ranger output to annotated cell types, chaining ambient-RNA removal, doublet detection, MAD-adaptive QC, normalization, integration, clustering, marker annotation, and (separately) pseudobulk DE + differential abundance. Use when honoring the made-once counting commitments (reference build/Ensembl vintage, --include-introns, cell-calling, feature namespace), ordering correct-then-detect-then-normalize (ambient before doublet before normalize), running per-sample QC before merge, integrating-then-clustering (never testing on integrated values), aggregating to pseudobulk for condition DE instead of cells-as-replicates, or pairing DE with differential abundance. Hands mechanism to the single-cell component skills; not a re-teach of any single step.
Install with Codex or Claude Copy this prompt, paste it into Codex, Claude, or another assistant, and let it review the skill page and install it for you.
A direct command skips the review prompt. Inspect the source before running it.
Orchestrates the end-to-end single-cell RNA-seq pipeline from 10x Cell Ranger output to annotated cell types, chaining ambient-RNA removal, doublet detection, MAD-adaptive QC, normalization, integration, clustering, marker annotation, and (separately) pseudobulk DE + differential abundance. Use when honoring the made-once counting commitments (reference build/Ensembl vintage, --include-introns, cell-calling, feature namespace), ordering correct-then-detect-then-normalize (ambient before doublet before normalize), running per-sample QC before merge, integrating-then-clustering (never testing on integrated values), aggregating to pseudobulk for condition DE instead of cells-as-replicates, or pairing DE with differential abundance. Hands mechanism to the single-cell component skills; not a re-teach of any single step.
Before using code patterns, verify installed versions match. If versions differ:
Python: pip show <package> then help(module.function) to check signatures
R: packageVersion('<pkg>') then ?function_name to verify parameters
If code throws ImportError, AttributeError, or TypeError, introspect the installed
package and adapt the example to match the actual API rather than retrying.
Note: the Cell Ranger filtered_feature_bc_matrix is cell-CALLING only, NOT ambient-corrected; recovering low-RNA cells or running SoupX/CellBender needs the RAW matrix. --include-introns is default TRUE since Cell Ranger 7 (essential for nuclei). gex_only=True/default silently drops Antibody Capture/CRISPR/HTO features. Confirm in-tool before quoting.
Single-Cell RNA-seq Pipeline
"Analyze my single-cell RNA-seq data from counts to cell types" -> Orchestrate QC filtering, normalization (scanpy/Seurat), batch integration (scVI/Harmony), clustering, marker detection, cell type annotation, and trajectory inference.
This is a workflow skill: it owns the chaining decisions and hand-offs, not the internals of any one step. Every step below cross-references the component skill that teaches its mechanism.
Made-once commitments (decided at cellranger count, inherited downstream)
Every gene ID/marker lookup/cross-dataset merge; join on gene_ids (stable), not symbols (lossy across releases)
--include-introns (default TRUE since CR 7)
Total UMI/cell and snRNA-vs-scRNA comparability; nuclei need introns; mixing intron-on/off samples is a hidden batch axis
Cell-vs-empty-droplet calling
The filtered matrix is cell-CALLING only, NOT ambient-corrected; SoupX/CellBender and low-RNA-cell rescue need the RAW matrix (filtered-only is irreversible loss)
Feature/barcode namespace (gex_only)
gex_only=True silently drops ADT/HTO/CRISPR features; set False and split by for CITE-seq/hashed pools
feature_types
Pipeline orchestration: the ordering decisions that make or break the result
Stage order is not arbitrary; getting it wrong silently corrupts every downstream step. The cross-cutting decisions this pipeline must get right:
Correct-then-detect-then-normalize, never reorder. Ambient-RNA removal (SoupX/CellBender) reads the RAW droplet matrix and MUST precede doublet detection and the final QC thresholds (ambient inflation distorts both doublet scores and per-cell QC); doublet detection precedes integration (doublets seed fake intermediate clusters). See single-cell/preprocessing.
Multiplexed (hashed) pools are demultiplexed FIRST. When samples are pooled with cell hashing (HTO/MULTI-seq) or genotype multiplexing, assign each cell to its sample of origin and remove cross-sample doublets before any per-sample step; the per-sample QC and doublet operations below assume the pool is already split by sample. See single-cell/hashing-demultiplexing.
Per-sample QC and doublet detection come BEFORE any merge or integration. Ambient-RNA cleanup (SoupX/CellBender), mito/gene-count filtering, and doublet calling (expected rate ~0.8% per 1,000 recovered cells) are per-capture operations; running them on a pooled object lets one lane's artifacts contaminate the shared null and lets surviving doublets form fake intermediate clusters. See single-cell/preprocessing and single-cell/doublet-detection.
Integrate, THEN cluster, THEN annotate - never reorder. For multi-sample designs, batch-correct on a shared embedding (Harmony/scVI/RPCA) and cluster on the corrected graph; clustering before correction makes clusters track samples or lanes instead of biology, and annotation has nothing to label until clusters exist. See single-cell/batch-integration, single-cell/clustering, single-cell/cell-annotation.
Over-integration erases biology. Batch-mixing metrics (kBET, iLISI) are maximized by destroying structure, so a method that flattens real cell states scores perfectly; pair every batch metric with a bio-conservation metric and re-inspect rare populations after correction. Batch-corrected expression is for embedding and clustering only - never feed it to differential expression. See single-cell/batch-integration.
Condition-level DE must use pseudobulk, not cells-as-replicates. Testing thousands of cells per donor as independent replicates is pseudoreplication and inflates false positives by orders of magnitude (Squair 2021); aggregate RAW counts per sample x cell type and hand them to DESeq2/edgeR/limma-voom. Step 7 marker p-values are descriptive ranking for labeling, not condition inference. See single-cell/markers-annotation and differential-expression/deseq2-basics.
Composition shifts are tested separately and masquerade as DE. A cluster-level expression change between conditions can be a pure proportion shift (a mixed cluster's substates re-balance with no gene changing per cell), invisible if only DE is run; always pair condition DE with a differential-abundance test (Milo cluster-free, or scCODA/sccomp/propeller cluster-based). See single-cell/differential-abundance.
The Seurat and Scanpy paths below are written single-sample for clarity. A multiplexed design demultiplexes the pool first (single-cell/hashing-demultiplexing); a multi-sample design then inserts per-sample QC and doublet removal, a merge, and an integration step between normalization (Step 4) and dimensionality reduction (Step 5); the rest of the order is unchanged.
Workflow Overview
10X data (RAW + filtered_feature_bc_matrix)
|
v
[0. Ambient removal] ---> SoupX/CellBender on RAW (per sample) (single-cell/preprocessing)
|
v
[1. Load Data] ---------> Read10X / read_10x_h5
|
v
[2. QC + Doublets] -----> MAD-adaptive nFeature/percent.mt; scDblFinder/Scrublet (per sample, before merge)
|
v
[3. Normalization] -----> SCTransform or LogNormalize
|
v
[4. HVG Selection] -----> FindVariableFeatures
|
v
[5. Dim Reduction] -----> PCA -> UMAP
|
v
[6. Clustering] --------> FindNeighbors -> FindClusters
|
v
[7. Markers] -----------> FindAllMarkers
|
v
[8. Annotation] --------> Manual or automated
|
v
Annotated Seurat/AnnData object
Primary Path: Seurat (R)
The Seurat/Scanpy paths below start from the filtered matrix for clarity. In practice, run ambient-RNA removal FIRST on the RAW droplet matrix per sample (SoupX needs raw+filtered+clusters; CellBender folds calling+denoise and is best for snRNA/high-ambient) — pick ONE (double-correction over-strips). See single-cell/preprocessing.
Regressing out percent.mt (or nCount/cell-cycle) is NOT reflexive: those covariates are confounded with real cell state and regressing them can erase biology, so only pass vars.to.regress for a covariate verified not confounded with the signal of interest. See single-cell/preprocessing.
# Set identity to chosen resolution
Idents(seurat_obj)<-'SCT_snn_res.0.4'# Find markers for all clusters
markers <- FindAllMarkers(seurat_obj, only.pos =TRUE, min.pct =0.25, logfc.threshold =0.25)# Top markers per cluster
top_markers <- markers %>%
group_by(cluster)%>%
slice_max(n =10, order_by = avg_log2FC)# Visualize top markers
DoHeatmap(seurat_obj, features = top_markers$gene)+ NoLegend()
Step 8: Cell Type Annotation
# Manual annotation based on known markers# Example for PBMC data:
cluster_annotations <-c('0'='CD4 T cells','1'='CD14 Monocytes','2'='B cells','3'='CD8 T cells','4'='NK cells','5'='CD16 Monocytes','6'='Dendritic cells')
seurat_obj$cell_type <- cluster_annotations[as.character(Idents(seurat_obj))]# Final UMAP
DimPlot(seurat_obj, reduction ='umap', group.by ='cell_type', label =TRUE)# Save object
saveRDS(seurat_obj,'seurat_annotated.rds')
Alternative Path: Scanpy (Python)
import scanpy as sc
import numpy as np
# Load 10X data
adata = sc.read_10x_h5('filtered_feature_bc_matrix.h5')
adata.var_names_make_unique()
# QC metrics
adata.var['mt'] = adata.var_names.str.startswith('MT-')
sc.pp.calculate_qc_metrics(adata, qc_vars=['mt'], percent_top=None, log1p=False, inplace=True)
# Filter (flat cutoffs are illustrative; prefer MAD-adaptive, tissue-aware thresholds, see single-cell/preprocessing)
sc.pp.filter_cells(adata, min_genes=200)
sc.pp.filter_genes(adata, min_cells=3)
adata = adata[adata.obs.n_genes_by_counts < 5000, :]
adata = adata[adata.obs.pct_counts_mt < 20, :]
# Doublet detection (rate from recovered cells, not the 0.05 placeholder, see single-cell/doublet-detection)
expected_rate = 0.008 * adata.n_obs / 1000
sc.pp.scrublet(adata, expected_doublet_rate=expected_rate)
adata = adata[~adata.obs['predicted_doublet'], :]
# Normalize and HVGs
sc.pp.normalize_total(adata, target_sum=1e4)
sc.pp.log1p(adata)
sc.pp.highly_variable_genes(adata, n_top_genes=2000)
# PCA, neighbors, UMAP -- stash log-normalized values in .raw BEFORE scaling, or rank_genes_groups# later computes logfoldchanges from z-scores (negative means -> NaN/garbage LFCs)
adata.raw = adata
adata = adata[:, adata.var.highly_variable]
sc.pp.scale(adata, max_value=10)
sc.tl.pca(adata, n_comps=50)
sc.pp.neighbors(adata, n_neighbors=15, n_pcs=30)
sc.tl.umap(adata)
# Clustering (pin the backend for reproducibility; default flips to igraph)
sc.tl.leiden(adata, resolution=0.5, flavor='igraph', n_iterations=2, directed=False)
# Markers
sc.tl.rank_genes_groups(adata, 'leiden', method='wilcoxon')
sc.pl.rank_genes_groups(adata, n_genes=10, sharey=False)
# Save
adata.write('scanpy_annotated.h5ad')
Parameter Recommendations
Step
Parameter
Recommendation
QC
min.features
200-500
QC
max.features
2500-5000 (depends on data)
QC
percent.mt
<10-20% (tissue-dependent)
SCTransform
vars.to.regress
none by default; only a validated non-confounded covariate
PCA
npcs
30-50
UMAP
dims
15-30 (check elbow plot)
Clustering
resolution
0.4-0.8 (start with 0.5)
QC cutoffs above are illustrative starting points, not universal thresholds: set min/max features and mito % per dataset from the data (MAD-adaptive, tissue-aware) rather than porting flat values, since healthy mito fraction varies by tissue. See single-cell/preprocessing.
Ambient markers everywhere (e.g. Hb across all PBMCs)
Ran SoupX/CellBender on the FILTERED matrix (soup already removed)
Run ambient removal on the RAW droplet matrix, before QC
Fake intermediate cell states
Doublets removed AFTER clustering/integration
Call doublets per sample before merge/normalize
Clusters track sample/lane, not biology
Clustered before integration
Integrate -> cluster -> annotate
Inflated DE, thousands of false positives
Tested cells as replicates (on integrated values)
Pseudobulk RAW counts per sample x cell-type -> DESeq2/edgeR/limma (Squair 2021)
"DE change" is really a proportion shift
Only ran DE, not abundance
Pair with differential abundance (Milo/scCODA/propeller)
Biology collapses to one blob
Reflexively regressed total_counts/cell-cycle
Regress only a validated, non-confounded covariate
CITE-seq/hashtag features vanish
gex_only=True default
Set False; split by feature_types
References
Heumos L, Schaar AC, Lance C, et al (2023) Best practices for single-cell analysis across modalities. Nature Reviews Genetics 24:550-572. DOI 10.1038/s41576-023-00586-w. (canonical pipeline order.)
Squair JW, Gautier M, Kathe C, et al (2021) Confronting false discoveries in single-cell differential expression. Nature Communications 12:5692. DOI 10.1038/s41467-021-25960-2. (cells-as-replicates is pseudoreplication; pseudobulk.)
Fleming SJ, Chaffin MD, Arduini A, et al (2023) Unsupervised removal of systematic background noise from droplet-based single-cell experiments using CellBender. Nature Methods 20:1323-1335. DOI 10.1038/s41592-023-01943-7. (ambient removal on the RAW matrix.)
Related Skills
database-access/geo-data - Resolve GSE to SRA; detect SuperSeries before processing
database-access/sra-data - Download 10x records with --include-technical for barcodes/UMIs
single-cell/data-io - Loading 10X, h5ad, RDS, and h5mu formats
