Scrnaseq scanpy core analysis

Single-Cell RNA-seq Core Analysis (Scanpy)

Complete workflow for single-cell RNA-seq analysis using Scanpy and the scverse ecosystem. Process raw data through quality control, normalization, clustering, and cell type annotation with publication-ready visualizations.

When to Use This Skill

• Analyze 10X Chromium data (CellRanger output, H5 files, raw/filtered matrices)
• Process Drop-seq, Smart-seq2, or inDrop single-cell RNA-seq data
• Integrate multi-batch data using scVI, scANVI, or Harmony
• Annotate cell types manually or with automated reference-based methods
• Compare conditions using pseudobulk differential expression (multi-sample data)

Don't use for: Bulk RNA-seq (use bulk-rnaseq-counts-to-de-deseq2), R-based scRNA-seq (use scrnaseq-seurat-core-analysis), Spatial transcriptomics (coming soon)

Installation

Package	Version	License	Commercial Use	Installation
scanpy	≥1.9	BSD-3-Clause	Permitted	pip install scanpy
anndata	≥0.8	BSD-3-Clause	Permitted	pip install anndata
numpy	≥1.20	BSD-3-Clause	Permitted	pip install numpy
pandas	≥1.3	BSD-3-Clause	Permitted	pip install pandas
matplotlib	≥3.4	PSF	Permitted	pip install matplotlib
seaborn	≥0.12	BSD-3-Clause	Permitted	pip install seaborn
adjustText	≥0.8	MIT	Permitted	pip install adjustText
scrublet	≥0.2.3	MIT	Permitted	pip install scrublet
scvi-tools	≥1.0	BSD-3-Clause	Permitted	pip install scvi-tools
harmonypy	≥0.0.9	GPL-3	Permitted	pip install harmonypy
celltypist	≥1.0	MIT	Permitted	pip install celltypist
pydeseq2	≥0.4	MIT	Permitted	pip install pydeseq2

Install all: pip install scanpy anndata numpy pandas matplotlib seaborn adjustText scrublet

Minimum versions: Python ≥3.8, scanpy ≥1.9, anndata ≥0.8

Inputs

Required:

• Raw or filtered count matrix: CellRanger output (filtered_feature_bc_matrix/), H5 files (.h5), AnnData (.h5ad), or count matrices (CSV/TSV)

Optional: Sample metadata (CSV/TSV) with sample IDs, conditions, batches, donor IDs

Data requirements: Min 500 cells/sample (1000+ recommended), Human or Mouse, UMI-based or read counts. See references/qc_guidelines.md for tissue-specific thresholds.

Outputs

Analysis objects:

• adata_processed.h5ad - Complete annotated AnnData object
  • Load with: import scanpy as sc; adata = sc.read_h5ad('adata_processed.h5ad')
  • Contains: raw counts, normalized data, QC metrics, clusters, cell types, UMAP/PCA
  • Note: adata.X contains log-normalized (not scaled) data. Scaled data is used internally for PCA but not stored in .X. This is correct for downstream analysis (DE, visualization).
  • Required for: trajectory inference, cell-cell communication, downstream analyses

Reports:

• scrna_analysis_report.pdf - Agent-generated comprehensive PDF with Methods, Results, Figures, Conclusions
• analysis_summary.txt - Text summary of dataset, QC, clustering, integration (generated by export_anndata_results())

⚠️ PDF style rules:

• US Letter page size (8.5 × 11 in) — always set page dimensions explicitly; do not rely on library defaults
• No Unicode superscripts — use 3.36e-06 or 3.36 × 10^(-6), not Unicode superscript chars (they render as ■ in PDF fonts)
• No half-empty pages — group headings with their content; only page-break before major sections (Results, Conclusions)
• Figures ≥80% page width — multi-panel figures must be large enough to read; never embed below 50% width

Tables: cell_metadata.csv, expression_matrix_counts.csv, expression_matrix_normalized.csv, pca_coordinates.csv, umap_coordinates.csv, cluster_markers_all.csv, {celltype}_deseq2_results.csv

Visualizations (PNG + SVG at 300 DPI): QC violins, UMAP plots, marker heatmaps, dot plots, volcano/MA plots

Clarification Questions

Default settings (use unless user specifies otherwise):

• Format: Filtered 10X CellRanger output | Species: Human | Tissue: PBMC
• Normalization: Standard (target sum + log1p) | Clustering: Test 0.4, 0.6, 0.8, 1.0

1. Input Files (ASK THIS FIRST):

• Do you have specific single-cell data file(s) to analyze?
• Expected: CellRanger output, H5, H5AD, or count matrix
• Or use example/demo data? (PBMC 3k dataset available)

2. Data format and species:

• (If using your own data) What format? Filtered 10X (default), Raw 10X, H5, H5AD? Human or Mouse? Tissue type?
• (If using example data) Defaults apply: filtered 10X, human, PBMC — skip to Q3

3. Batch structure:

• a) Single sample (no integration needed)
• b) Multiple batches (scVI recommended, Harmony for speed)

4. Analysis scope:

• a) Standard: QC + normalize + cluster + annotate (recommended)
• b) Standard + pseudobulk DE (requires ≥2 samples per condition)
• c) Custom: specify which steps

5. Clustering granularity:

• a) Coarse (0.3-0.5) — major cell types only
• b) Standard (0.6-0.8) — recommended
• c) Fine (1.0-1.5) — subtypes
• d) Test multiple: 0.4, 0.6, 0.8, 1.0 (recommended)

Standard Workflow

Note: Run from the OmicsClaw root directory and add the workflow scripts to sys.path:

import sys; import os; sys.path.insert(0, os.path.abspath('knowledge_base/scripts/scrnaseq-scanpy-core-analysis'))

🚨 MANDATORY: USE SCRIPTS EXACTLY AS SHOWN - DO NOT WRITE INLINE CODE 🚨

Detailed step-by-step code: references/workflow-details.md

CRITICAL - DO NOT:

• Write inline analysis code → STOP: Use the script functions
• Write custom export code → STOP: Use export_anndata_results()
• Skip verification messages → STOP: Check for "✓" messages after each step

IF SCRIPTS FAIL - Script Failure Hierarchy:

1. Fix and Retry (90%) - Install missing package, re-run script
2. Modify Script (5%) - Edit the script file itself, document changes
3. Use as Reference (4%) - Read script, adapt approach, cite source
4. Write from Scratch (1%) - Only if genuinely impossible, explain why

NEVER skip directly to writing inline code without trying the script first.

---

Step 1 — Load and QC | knowledge_base/scripts/scrnaseq-scanpy-core-analysis/setup_and_import.py, scripts/qc_metrics.py, scripts/filter_cells.py

# Load data (Option A: example, Option B: your data)
from load_example_data import load_example_data
adata = load_example_data("pbmc3k")

# QC metrics + adaptive filtering + doublet detection
from qc_metrics import calculate_qc_metrics, batch_mad_outlier_detection
from filter_cells import run_scrublet_detection, filter_by_mad_outliers
adata = calculate_qc_metrics(adata, species="human")
adata = batch_mad_outlier_detection(adata, batch_key="batch")  # creates 'outlier' column
adata = run_scrublet_detection(adata, batch_key="batch")
adata = filter_by_mad_outliers(adata, remove_doublets=True)

DO NOT write inline QC code. Doublet rate auto-scales per batch (~0.8% per 1,000 cells). Aim for >70% cell retention. For raw data, prepend ambient RNA correction: references/ambient_rna_correction.md

✅ VERIFICATION: "✓ Data loaded successfully!" → QC metrics added → filtering summary with retention %.

Step 2 — Normalize, reduce, integrate | scripts/normalize_data.py, scripts/scale_and_pca.py, scripts/integrate_scvi.py

from normalize_data import run_standard_normalization
from find_variable_genes import find_highly_variable_genes
from scale_and_pca import scale_data, run_pca_analysis
adata = run_standard_normalization(adata, target_sum=1e4)
adata = find_highly_variable_genes(adata, n_top_genes=2000)
adata = scale_data(adata, vars_to_regress=["total_counts", "pct_counts_mt"])
adata = run_pca_analysis(adata, n_pcs=50)

# Multi-batch only: integration + diagnostics
from integrate_scvi import run_scvi_integration
from integration_diagnostics import compute_lisi_scores
adata = run_scvi_integration(adata, batch_key="batch", condition_key="condition")
lisi = compute_lisi_scores(adata, batch_key="batch", use_rep="X_scVI")

DO NOT write inline normalization or integration code. The integration script auto-detects batch-condition confounding. Details →

✅ VERIFICATION:

• Normalization: "✓ Normalization complete"
• PCA: "PCA loadings verified: N HVG rows have non-zero loadings" — if you see a WARNING about zero loadings, re-run PCA with use_highly_variable=True
• After PCA, call suggest_n_pcs(adata) to get recommended PC count for Step 3
• Integration (multi-batch): LISI scores printed

Step 3 — Cluster, annotate, visualize | scripts/cluster_cells.py, scripts/find_markers.py, scripts/annotate_celltypes.py

from cluster_cells import build_neighbor_graph, cluster_leiden_multiple_resolutions
from run_umap import run_umap_reduction
from find_markers import find_all_cluster_markers
from plot_dimreduction import plot_umap_clusters
use_rep = "X_scVI" if "X_scVI" in adata.obsm else "X_pca"
# Default n_pcs=30 is standard. NEVER use <15 PCs.
adata = build_neighbor_graph(adata, use_rep=use_rep, n_neighbors=10, n_pcs=30)
adata = cluster_leiden_multiple_resolutions(adata, resolutions=[0.4, 0.6, 0.8, 1.0])
adata = run_umap_reduction(adata)
markers = find_all_cluster_markers(adata, cluster_key="leiden_0.8")
plot_umap_clusters(adata, cluster_key="leiden_0.8", output_dir="results/umap")

# Annotate (manual or CellTypist)
from annotate_celltypes import annotate_clusters_manual
annotations = {"0": "CD4 T cells", "1": "CD14+ Monocytes", ...}
adata = annotate_clusters_manual(adata, annotations, cluster_key="leiden_0.8")

DO NOT write inline clustering or annotation code. CellTypist validates labels post-hoc (flags suspect ILC/HSC, contamination, low-complexity). Markers →

⚠️ n_pcs for neighbor graph: Default is 30 PCs (standard). Using <15 PCs risks collapsing distinct populations. If you used suggest_n_pcs() in Step 2, pass that value here.

For pseudobulk DE (multi-sample, ≥2 replicates/condition): scripts/pseudobulk_de.py. Script blocks with N=1. Details →

✅ VERIFICATION:

• Cluster counts → UMAP plots saved → marker genes identified → cell type annotations added
• After find_all_cluster_markers(): verify the returned DataFrame columns and print markers.head() to confirm values are sensible before saving to CSV.

Step 4 — Export results | scripts/export_results.py

from export_results import export_anndata_results

export_anndata_results(adata, output_dir="results", cluster_key="cell_type")

DO NOT write custom export code. Use export_anndata_results().

Exports: H5AD, expression matrices (raw + normalized CSV), cell metadata, UMAP/PCA coordinates, text summary.

✅ VERIFICATION: You MUST see:

=== Export Complete ===

If you don't see "Export Complete": The export did not complete. Re-run the export function.

⚠️ Large datasets (>20k cells): H5AD export may take 10-60 seconds depending on dataset size. The script prints progress updates (size estimate, compression method, elapsed time). If export appears stuck for >2 minutes, interrupt and retry with save_h5ad(adata, path, compression=None) to skip compression.

Decision Guide

Decision	Quick Guide	Reference
Ambient RNA	Skip for filtered/PBMC. CellBender for raw/high-soup (brain, lung, tumor)	ambient_rna_correction.md
QC Strategy	MAD (multi-batch). Fixed (single batch, tissue-specific)	qc_guidelines.md
Normalization	Standard (most data). Pearson (heteroscedastic)	scanpy_best_practices.md
Integration	scVI (complex batches). Harmony (fast, simple)	integration_methods.md
Resolution	Test 0.4, 0.6, 0.8, 1.0. Choose by biology and stability	scanpy_best_practices.md
Annotation	Manual (accurate). CellTypist (fast). Both (validate)	marker_gene_database.md

Complete working examples: references/common-patterns.md

Common Issues

Issue	Cause	Solution
ImportError: No module named 'scanpy'	Not installed	pip install scanpy anndata numpy pandas matplotlib seaborn
Low cell retention (<70%)	Strict QC thresholds	Use MAD (nmads=5→7) or tissue-specific thresholds
Out of memory	Large dataset (>50k cells)	Use backed mode or subsample
Clusters driven by batch	Insufficient integration	Use scVI, increase n_latent, check confounding
Poor UMAP separation	Wrong parameters	Check PCA elbow, use 20-40 PCs, adjust n_neighbors
High MT% in all cells	Degradation or tissue-specific	Check distribution — bimodal: stricter filter; uniform: may be biological
FileNotFoundError: barcodes.tsv.gz	Wrong directory	Verify 10X output files present. Use import_h5_data() for .h5
H5AD export hangs or is very slow	Large file write with compression	Normal for >50 MB files. Script uses fast lzf compression. If still slow, pass compression=None to save_h5ad().
NameError after export interruption	Kernel restart lost variables	Re-run from Step 1 to restore adata. Export is idempotent — safe to re-run.

Expected warnings (not errors):

Warning	Meaning	Action
SVG export failed	Optional SVG dependency unavailable	Normal — PNG always generated. Both created in most environments.
Detected doublet rate differs from expected	Scrublet threshold or pre-filtering	Inspect adata.obs['doublet_score'].hist(). Adjust threshold if needed.
Pseudobulk DE blocked: N=1	Insufficient replicates	Need ≥2 per condition. Use cell-level Wilcoxon for exploratory only.
Batch-condition confounding	Condition has 1 sample	Clustering valid. Composition comparisons need caveats.
CellTypist labels suspect (ILC, HSC)	Automated misclassification	Cross-check markers. Relabel ILC→NK or HSC→Unassigned.
Multimodal data detected	RNA-only workflow on CITE-seq	Note in reports: "RNA modality only; ADT not analyzed."

Detailed troubleshooting: references/troubleshooting_guide.md

Suggested Next Steps

1. Functional Enrichment — functional-enrichment-from-degs for pathway analysis of DE results
2. Trajectory Analysis — PAGA, Palantir, or scVelo for developmental datasets
3. Cell-Cell Communication — CellPhoneDB, LIANA, or NicheNet for ligand-receptor interactions

Related Skills

Alternative: scrnaseq-seurat-core-analysis (R-based) | Downstream: functional-enrichment-from-degs, de-results-to-plots, de-results-to-gene-lists | Complementary: bulk-omics-clustering, experimental-design-statistics

References

1. Scanpy: Wolf FA, et al. (2018) Genome Biol. 19:15.
2. Best Practices: Luecken MD, Theis FJ. (2019) Mol Syst Biol. 15:e8746.
3. Pseudobulk DE: Squair JW, et al. (2021) Nat Commun. 12:5692.
4. scVI: Lopez R, et al. (2018) Nat Methods. 15:1053-1058.

Detailed guides: workflow-details.md | common-patterns.md | scanpy_best_practices.md | qc_guidelines.md | integration_methods.md | pseudobulk_de_guide.md | marker_gene_database.md | troubleshooting_guide.md

Scripts: scripts/ | Evaluation: assets/eval/complete_example_analysis.py