Prompt Engineering for Seurat and Scanpy: A Reproducible Workflow

Learn to specify single-cell analysis tasks for optional AI-assisted code drafting. State the data, biological question, method, and outputs, then review any proposed code and inspect results if you choose to run it.

What this lesson provides

This lesson provides a task-specific prompt, a primary Seurat marker-analysis template in R, and a separate optional Scanpy QC/clustering template in Python. The examples are not equivalent implementations: they require different prepared inputs and answer different questions. A generated script and its outputs are not supplied; using an assistant and running either example are optional.

Learning objectives

• Learn how to specify reproducible Seurat (R) and Scanpy (Python) analysis tasks for optional AI-assisted code drafting.
• Be able to craft a prompt that requests code with dataset context, a biological goal, practical constraints, and concrete output checks, then verify any proposed code against those requirements.
• Understand how to validate any proposed code against explicit checks and inspect outputs if you choose to run it.

Prerequisites

• Basic familiarity with single-cell RNA-seq (scRNA-seq) data concepts: cells, genes, counts, metadata, clustering, and differential expression.
• An AI coding environment is needed only if you choose to request generated code. Running a supplied template requires local R or Python and the corresponding Seurat or Scanpy dependencies.
• Command-line literacy is useful if you choose to run scripts and save outputs to explicit file paths.

Why intent-based prompting matters for Seurat & Scanpy workflows

A prompt should specify the question and the constraints that affect its answer, then ask the assistant to state assumptions and produce checks you can run. Intent-based prompting (a.k.a. "vibe coding") focuses the assistant on what you want to achieve (for example, "identify markers distinguishing cluster 2") and leaves the assistant to choose reasonable API calls within explicit constraints you provide. A prompt does not remove the need to verify API calls against the installed package version.

In single-cell analysis, small changes in preprocessing (filtering thresholds, normalization, variable-gene selection) materially affect downstream results. Therefore, prompt templates should require the assistant to state assumptions, include checks, and write files that let you re-run analyses deterministically.

Reproducible prompting loop

Before asking an assistant to write code, establish a reproducible loop:

• Provide a concise dataset summary and file paths.
• State the biological intent (what question you want answered).
• Define constraints (statistical test, thresholds, plotting formats, language).
• Require output files and programmatic checks (dimensions, column names, saved plots).
• Run the script, inspect outputs, and iterate with small changes.

Use this prompt as a specification, not as evidence that an assistant will return correct code. The prompt below uses the same input and output locations as the Seurat template: ./data/seurat_obj.rds and ./results/cluster2_markers.csv, plus the cluster-colored UMAP at ./results/umap_seurat_clusters.png. Replace assumptions with facts from your object; biological results cannot be guaranteed by a minimum row count.

Goal: Identify marker genes for cluster 2 vs all in a PBMC dataset using the Seurat example below.
Input: RDS Seurat object at ./data/seurat_obj.rds; inspect and report its actual dimensions; metadata column 'seurat_clusters' must be present.
Constraints: Use R and Seurat >=4; test = MAST; logfc.threshold = 0.25; return genes with p_val_adj < 0.05; save a CSV report.
Output: ./results/cluster2_markers.csv and UMAP plot ./results/umap_seurat_clusters.png
Checks: verify the file exists, inspect actual column names and row count (zero rows can be valid), confirm the intended identities were compared, and open the UMAP PNG.
Ask: show all package imports and version prints; explain any assumptions and cite Seurat functions used.

Reproducible, sequential workflow examples

The R Seurat workflow below is the main route and implements the cluster-2 marker prompt. The optional Python Scanpy workflow is a separate QC → clustering → exploratory marker-ranking task; it is not an equivalent implementation of the Seurat prompt or a shared sequential pipeline. Both are templates for prepared objects.

Main workflow: Seurat cluster markers (R)

Prepare a serialized Seurat object at ./data/seurat_obj.rds, with an RNA expression assay and a seurat_clusters metadata column. Check its structure from the project root before running the workflow:

Rscript -e 'x <- readRDS("./data/seurat_obj.rds"); print(dim(x)); print(colnames(x@meta.data)); stopifnot("seurat_clusters" %in% colnames(x@meta.data))'

Expect nonzero feature-by-cell dimensions and the required metadata field. This structural check does not establish that counts, assay choice, identities, or biological annotations are appropriate. If the object or required field is missing, stop and prepare a compatible object first. The existing seurat_clusters labels define the marker comparison; if UMAP is absent, the script computes a visualization without reclustering or replacing those labels.

• Input Seurat object: ./data/seurat_obj.rds (serialized Seurat object from saveRDS)
• Output directory: ./results/

# Save this as ./scripts/seurat_marker_workflow.R
library(Seurat)
library(ggplot2)
library(MAST) # used for statistical test
library(dplyr)
library(tibble) # provides rownames_to_column

# Paths
infile <- "./data/seurat_obj.rds"
outdir <- "./results"
dir.create(outdir, showWarnings = FALSE, recursive = TRUE)

# Load object and report
seurat_obj <- readRDS(infile)
cat("Seurat object dimensions (cells x genes):", dim(seurat_obj), "\n")
cat("Available metadata columns:", paste(colnames(seurat_obj@meta.data), collapse=", "), "\n")

# Validate metadata before plotting or using the cluster identities.
if (!"seurat_clusters" %in% colnames(seurat_obj@meta.data)) {
  stop("Required metadata column 'seurat_clusters' is missing; prepare or select the intended identities first.")
}

# Record the incoming biological labels and intended marker comparison.
incoming_cluster_labels <- unname(as.character(seurat_obj[[]][["seurat_clusters"]]))
if (anyNA(incoming_cluster_labels)) {
  stop("Required metadata column 'seurat_clusters' contains missing labels.")
}
cluster_of_interest <- "2"
comparison_metadata <- list(
  identity_field = "seurat_clusters",
  cluster_of_interest = cluster_of_interest,
  ident.1 = cluster_of_interest,
  ident.2 = NULL,
  reference = "all other cells",
  test = "MAST",
  assay = DefaultAssay(seurat_obj),
  logfc.threshold = 0.25,
  min.pct = 0.1
)

# Compute a visualization only when missing. Do not call FindClusters here:
# a new clustering would silently replace the required marker-comparison labels.
if (!"umap" %in% names(seurat_obj@reductions)) {
  seurat_obj <- NormalizeData(seurat_obj, verbose = FALSE)
  seurat_obj <- FindVariableFeatures(seurat_obj, selection.method = "vst", nfeatures = 2000, verbose = FALSE)
  seurat_obj <- ScaleData(seurat_obj, verbose = FALSE)
  seurat_obj <- RunPCA(seurat_obj, verbose = FALSE)
  available_pcs <- ncol(Embeddings(seurat_obj, reduction = "pca"))
  if (available_pcs < 2L) {
    stop("At least two PCA dimensions are needed to compute the UMAP visualization.")
  }
  umap_dims <- seq_len(min(20L, available_pcs))
  seurat_obj <- FindNeighbors(seurat_obj, dims = umap_dims, verbose = FALSE)
  seurat_obj <- RunUMAP(seurat_obj, dims = umap_dims, verbose = FALSE)
}

# Visualization work must leave the input biological labels unchanged.
if (!identical(unname(as.character(seurat_obj[[]][["seurat_clusters"]])), incoming_cluster_labels)) {
  stop("The incoming seurat_clusters labels changed while preparing the UMAP.")
}

# Save UMAP plot (colored by clusters)
p <- DimPlot(seurat_obj, reduction = "umap", group.by = "seurat_clusters") + ggtitle("UMAP by seurat_clusters")
ggsave(filename = file.path(outdir, "umap_seurat_clusters.png"), plot = p, width = 6, height = 5)

# Select and persist the explicitly intended comparison metadata.
Idents(seurat_obj) <- comparison_metadata$identity_field
comparison_metadata$observed_identity_labels <- sort(unique(as.character(seurat_obj[[]][["seurat_clusters"]])))
seurat_obj@misc$cluster2_marker_comparison <- comparison_metadata
print(comparison_metadata)

# Exploratory cell-level marker test. For condition-level inference across samples,
# use a sample-aware design rather than treating cells as replicates.
markers <- FindMarkers(seurat_obj,
                       ident.1 = comparison_metadata$ident.1,
                       ident.2 = comparison_metadata$ident.2,
                       assay = comparison_metadata$assay,
                       test.use = comparison_metadata$test,
                       logfc.threshold = comparison_metadata$logfc.threshold,
                       min.pct = comparison_metadata$min.pct)

# Filter by adjusted p-value
significant_markers <- markers %>% rownames_to_column(var = "gene") %>% filter(p_val_adj < 0.05)
write.csv(significant_markers, file = file.path(outdir, paste0("cluster", cluster_of_interest, "_markers.csv")), row.names = FALSE)

# Save computed reductions, unchanged incoming labels, and comparison metadata.
saveRDS(seurat_obj, file = file.path(outdir, "seurat_obj_processed.rds"))
message("Wrote results to ", outdir)

What each command produces:

• ./results/umap_seurat_clusters.png: UMAP colored by 'seurat_clusters'.
• ./results/cluster2_markers.csv: exploratory differential-expression results for cluster 2 versus the remaining cells, filtered by adjusted p-value.
• ./results/seurat_obj_processed.rds: Seurat object with any computed reductions, unchanged incoming cluster labels, and the recorded cluster-2 comparison metadata.

Optional Python route: the following uses a separate AnnData input and performs its own QC and clustering; do not compare its outputs as if they came from the same preprocessing or marker-test contract as the Seurat analysis.

Optional, distinct workflow: Scanpy QC and exploratory markers (Python)

Prepare an AnnData file at ./data/adata.h5ad with gene symbols in var_names and count-like values in .X; this template computes QC metrics but does not create the object. From the project root, check the input before running the workflow:

python -c "import scanpy as sc; a=sc.read_h5ad('./data/adata.h5ad'); print(a.shape); print(list(a.obs.columns)); print(a.var_names[:5].tolist())"

Confirm that the printed dimensions are nonzero, feature names are the expected organism-specific identifiers, and .X contains the intended count-like data. If the file is absent or the structure is different, stop and prepare or locate a compatible object rather than proceeding.

• Input AnnData file: ./data/adata.h5ad. This example expects gene symbols in var_names and count-like values in .X; it calculates n_genes_by_counts and pct_counts_mt itself.
• Output directory: ./results/ (create if missing)

# Save this as ./scripts/scanpy_quick_workflow.py
import scanpy as sc
import matplotlib.pyplot as plt
import anndata
import os
import sys

# Input/outputs (change paths as needed)
adata_path = "./data/adata.h5ad"
outdir = "./results"
os.makedirs(outdir, exist_ok=True)

# Load and print basic info
adata = sc.read_h5ad(adata_path)
print("AnnData shape:", adata.shape)
print("obs columns:", list(adata.obs.columns)[:20])

# Calculate QC metrics (adds n_genes_by_counts, total_counts, pct_counts_mt if mt genes detected)
# Example convention for human symbols only; adapt for organism and annotation.
adata.var['mt'] = [name.upper().startswith('MT-') for name in adata.var_names]
sc.pp.calculate_qc_metrics(adata, qc_vars=['mt'], inplace=True)

# Plot QC violin (before filtering)
sc.pl.violin(adata, ['n_genes_by_counts', 'total_counts', 'pct_counts_mt'],
             jitter=0.4, multi_panel=True, show=False)
plt.savefig(os.path.join(outdir, "qc_violins_before.png"))
plt.clf()

# Illustrative thresholds only. Inspect distributions and choose per dataset.
min_genes = 200
max_genes = 8000
max_mt_pct = 5.0

sc.pp.filter_cells(adata, min_genes=min_genes)
sc.pp.filter_cells(adata, max_genes=max_genes)
adata = adata[adata.obs['pct_counts_mt'] < max_mt_pct, :]

sc.pp.filter_genes(adata, min_cells=3)

# Normalize counts, then preserve an unscaled log-normalized source before HVG selection.
sc.pp.normalize_total(adata, target_sum=1e4)
sc.pp.log1p(adata)
adata.layers["log1p_norm"] = adata.X.copy()
# Keep adata full-gene for ranking; use an HVG-only copy for scaled embeddings.
# Any input .raw is deliberately ignored by the explicit layer selection below.
sc.pp.highly_variable_genes(adata, n_top_genes=min(2000, adata.n_vars), subset=False)
adata_hvg = adata[:, adata.var["highly_variable"]].copy()
sc.pp.scale(adata_hvg, max_value=10)
available_pcs = min(adata_hvg.n_obs - 1, adata_hvg.n_vars - 1)
if available_pcs < 2:
    raise ValueError("Need at least two available PCs after cell and HVG filtering.")
n_pcs = min(40, available_pcs)
sc.tl.pca(adata_hvg, n_comps=n_pcs, svd_solver='arpack')
sc.pp.neighbors(adata_hvg, n_neighbors=10, n_pcs=n_pcs)
sc.tl.umap(adata_hvg)
sc.tl.leiden(adata_hvg, key_added='leiden_r10')
# Carry embedding labels and coordinates back to the full-gene object.
adata.obs['leiden_r10'] = adata_hvg.obs['leiden_r10'].copy()
adata.obsm['X_umap'] = adata_hvg.obsm['X_umap'].copy()

# Save UMAP plot directly to the requested results path.
sc.pl.umap(adata, color='leiden_r10', show=False)
plt.savefig(os.path.join(outdir, "umap_leiden_r10.png"), bbox_inches='tight')
plt.close()

# Exploratory cell-level ranking for each Leiden group versus the remaining cells.
# For replicated condition-level inference use a sample-aware design, often pseudobulk.
# Rank all genes retained after filtering from the saved, unscaled log1p-normalized layer.
# Explicit use_raw=False prevents a pre-existing .raw from silently changing the source.
sc.tl.rank_genes_groups(adata, 'leiden_r10', method='t-test_overestim_var',
                        reference='rest', layer='log1p_norm', use_raw=False)
# Convert to a readable table and save
import pandas as pd
r = adata.uns['rank_genes_groups']
groups = r['names'].dtype.names
records = []
for group in groups:
    names = r['names'][group].tolist()
    pvals_adj = r['pvals_adj'][group].tolist()
    logfoldchanges = r['logfoldchanges'][group].tolist() if 'logfoldchanges' in r else [None]*len(names)
    for name, p, lf in zip(names, pvals_adj, logfoldchanges):
        records.append({'group': group, 'gene': name, 'p_val_adj': p, 'logfoldchange': lf})
df = pd.DataFrame.from_records(records)
df.to_csv(os.path.join(outdir, "scanpy_rank_genes.csv"), index=False)

# Save AnnData with processed fields
adata.write_h5ad(os.path.join(outdir, "adata_processed.h5ad"))
print("Wrote results to", outdir)

What each command produces:

• ./results/qc_violins_before.png: violin plots of QC metrics before filtering.
• Filtered in-memory AnnData (then saved): ./results/adata_processed.h5ad
• UMAP image: ./results/umap_leiden_r10.png (the script explicitly saves this file).
• Cluster ranking table: ./results/scanpy_rank_genes.csv (all genes retained after filtering, not only HVGs; computed from the log1p_norm layer).

Check your reasoning: labels and expression scale

Task. Your Seurat input already has seurat_clusters values 0, 1, 2, but no UMAP. Which labels define cluster 2 versus the reference group after the visualization step? For Scanpy, which values feed marker ranking if the input already has a .raw field?

Validation boundary: focused synthetic checks exercised label preservation with R 4.3.3/Seurat 5.0.1 and expression-source selection/full-retained-gene ranking with Python 3.12.3/Scanpy 1.11.5. The Seurat check substituted Wilcoxon for MAST because MAST was unavailable. Neither complete script nor the MAST path was executed; these remain advanced templates, not biologically validated analyses.

Output inspection and interpretation

After running a generated script, inspect outputs with these checks and interpret them cautiously.

• QC plots Check whether the distributions of gene counts, total counts, and mitochondrial percentage look sensible for your sample type. Very high mitochondrial percentages often indicate damaged cells, but thresholds are dataset-specific; 5% is typical for fresh PBMCs and not a universal rule.
• Filtering results: Note the number of cells and genes retained. If you lose an unexpectedly large fraction (>30–50%) of cells, revisit thresholds and confirm metadata quality (e.g., correct gene name capitalization for mitochondrial genes).
• Clustering and UMAP: UMAP and cluster labels are descriptive summaries. A clear separation of known cell-type markers across clusters increases confidence, but cluster boundaries can change with preprocessing or neighbor parameter changes.
• Differential expression tables: Check the installed version's actual columns, group definition, assay or layer, and test. Adjusted p-values describe evidence under that test and correction; they do not establish causality. Cell-level comparisons can overstate evidence when cells from one sample are treated as independent replicates.

What plausible results mean and what they don't:

• Marker genes with strong log fold-change and low p_adj suggest genes enriched in that cluster under analysis choices (normalization, test). They do not by themselves validate cell identity; you should cross-reference with known marker gene sets and, if possible, orthogonal assays.
• Absence of significant markers could mean clusters are not biologically distinct, the cluster is too small, or the test lacks power (low cell counts or high variability).

Troubleshooting: symptoms, likely causes, and fixes

• Symptom: Script errors on loading files (FileNotFoundError / cannot open file).
  Cause: Wrong relative path or missing files.
  Fix: Verify file existence with ls or dir; use absolute paths or check working directory in script (e.g., print getwd() in R / os.getcwd() in Python).
• Symptom: QC metrics like pct_counts_mt are all zeros or missing.
  Cause: Mitochondrial gene names do not match detection pattern (e.g., 'MT-' vs 'mt-').
  Fix: Inspect var/feature names and adjust the pattern detection (toupper/lower) or provide a list of mitochondrial gene symbols.
• Symptom: Very few or no variable genes selected.
  Cause: Raw counts not provided or log-normalized values present; params for highly variable gene selection may be inappropriate.
  Fix: Ensure you run appropriate normalization steps before HVG selection or set selection.method/nfeatures to reasonable values.
• Symptom: DE test fails with NA p-values or warnings about zeros.
  Cause: Low-expression genes or small cluster sizes make statistical tests unstable; some tests require non-zero variance.
  Fix: Filter genes (min.pct, min.cells), or switch to a test more appropriate for small groups (e.g., Wilcoxon, or aggregate pseudo-bulk approaches).
• Symptom: UMAP layout looks collapsed or clusters overlap completely.
  Cause: PCA not computed or insufficient PCs used; scaling issues.
  Fix: Recompute PCA, increase n_pcs, or re-run neighbor graph construction; check for batch effects that may need correction.

Limitations, version and reference-data caveats, and responsible use

Limitations and practical caveats:

• AI-generated code is only as good as the prompt and the training distribution of the model. Always read, run, and vet generated scripts before applying them to sensitive data.
• Different versions of Seurat, Scanpy, and dependent packages change defaults and returned structures. Record package versions, consult installed documentation, and check actual column names instead of assuming a fixed output schema.
• Single-cell pipelines are sensitive to preprocessing. Cluster marker tests compare cells and may inflate significance when cells are not independent biological replicates; compare conditions with a sample-aware design. Exploratory results are hypothesis-generating, not confirmatory.
• Public gene name conventions (e.g., capitalization in mitochondrial genes) and reference annotations differ by organism and build. Do not assume gene naming conventions without checking your dataset.

Responsible-use notes:

• Never execute code that performs destructive file operations (rm, unlink) without explicit manual review. If an assistant proposes file deletion, decline or require safe-guards (e.g., move to backup directories instead).
• Check network calls in generated code. If code downloads reference data, ensure the source URL is trustworthy and that you are allowed to access it in your environment.
• Do not paste identifiable or otherwise protected human data into an external AI service unless institutional policies, approvals, and service terms permit it. Share schema or synthetic examples instead when possible.

Concise next-step checklist

1. Before requesting code, record each prepared object's dimensions, metadata columns, feature identifiers, assay or layer, and provenance; stop if required fields or intended data are unclear.
2. Use the prompt template at the top of this tutorial to request an initial script from your AI assistant. Require that the assistant prints package versions and lists its assumptions.
3. Run the generated script locally, inspect the outputs listed in the prompt, and record failures or warnings.
4. Iterate with small, targeted prompts (e.g., "change logfc.threshold to 0.25" or "use Wilcoxon instead of MAST") and rerun checks.
5. Once stable, save the final prompt and script together in version control (e.g., ./prompts/cluster2_prompt.txt and ./scripts/cluster2_analysis.R).

References

• Seurat official documentation (Satija Lab)
• Scanpy documentation
• MAST package (Bioconductor)
• Stuart et al., Comprehensive integration of single-cell data (Seurat v3 paper)