Scrnaseq Scanpy Core Analysis

---
id: scrnaseq-scanpy-core-analysis
name: Single-Cell RNA-seq Core Analysis (Scanpy)
category: transcriptomics
short-description: "Complete single-cell RNA-seq analysis using Scanpy from raw data to cell type annotation with clustering and visualization."
detailed-description: "Complete single-cell RNA-seq analysis using Scanpy from raw data to cell type annotation. Use when you have 10X Chromium, Drop-seq, or other scRNA-seq data requiring QC, normalization, clustering, and visualization. Implements current best practices including ambient RNA correction (CellBender), batch-aware adaptive QC (MAD), doublet detection (Scrublet), standard or Pearson residuals normalization, batch integration (scVI/Harmony), multi-resolution Leiden clustering, and pseudobulk differential expression for condition comparisons. Best for human or mouse data with 500+ cells per sample. Produces publication-ready plots and annotated AnnData objects."
starting-prompt: "Analyze single-cell RNA-seq data with Scanpy from QC through cell type annotation. Generate a PDF report with an intro, methods, results, conclusions and figures from all of the analyses you perform."
---

# 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](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`:
> ```python
> 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](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](knowledge_base/scripts/scrnaseq-scanpy-core-analysis/setup_and_import.py), [scripts/qc_metrics.py](knowledge_base/scripts/scrnaseq-scanpy-core-analysis/qc_metrics.py), [scripts/filter_cells.py](knowledge_base/scripts/scrnaseq-scanpy-core-analysis/filter_cells.py)

```python
# 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](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](knowledge_base/scripts/scrnaseq-scanpy-core-analysis/normalize_data.py), [scripts/scale_and_pca.py](knowledge_base/scripts/scrnaseq-scanpy-core-analysis/scale_and_pca.py), [scripts/integrate_scvi.py](knowledge_base/scripts/scrnaseq-scanpy-core-analysis/integrate_scvi.py)

```python
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 →](references/integration_methods.md)

**✅ 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](knowledge_base/scripts/scrnaseq-scanpy-core-analysis/cluster_cells.py), [scripts/find_markers.py](knowledge_base/scripts/scrnaseq-scanpy-core-analysis/find_markers.py), [scripts/annotate_celltypes.py](knowledge_base/scripts/scrnaseq-scanpy-core-analysis/annotate_celltypes.py)

```python
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 →](references/marker_gene_database.md)

**⚠️ 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](knowledge_base/scripts/scrnaseq-scanpy-core-analysis/pseudobulk_de.py). Script blocks with N=1. [Details →](references/pseudobulk_de_guide.md)

**✅ 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](knowledge_base/scripts/scrnaseq-scanpy-core-analysis/export_results.py)

```python
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](references/ambient_rna_correction.md) |
| **QC Strategy** | MAD (multi-batch). Fixed (single batch, tissue-specific) | [qc_guidelines.md](references/qc_guidelines.md) |
| **Normalization** | Standard (most data). Pearson (heteroscedastic) | [scanpy_best_practices.md](references/scanpy_best_practices.md) |
| **Integration** | scVI (complex batches). Harmony (fast, simple) | [integration_methods.md](references/integration_methods.md) |
| **Resolution** | Test 0.4, 0.6, 0.8, 1.0. Choose by biology and stability | [scanpy_best_practices.md](references/scanpy_best_practices.md) |
| **Annotation** | Manual (accurate). CellTypist (fast). Both (validate) | [marker_gene_database.md](references/marker_gene_database.md) |

**Complete working examples:** [references/common-patterns.md](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](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](references/workflow-details.md) | [common-patterns.md](references/common-patterns.md) | [scanpy_best_practices.md](references/scanpy_best_practices.md) | [qc_guidelines.md](references/qc_guidelines.md) | [integration_methods.md](references/integration_methods.md) | [pseudobulk_de_guide.md](references/pseudobulk_de_guide.md) | [marker_gene_database.md](references/marker_gene_database.md) | [troubleshooting_guide.md](references/troubleshooting_guide.md)

**Scripts:** [scripts/](knowledge_base/scripts/scrnaseq-scanpy-core-analysis/) | **Evaluation:** [assets/eval/complete_example_analysis.py](assets/eval/complete_example_analysis.py)