mims-harvard/tooluniverse/tooluniverse-single-cell

Single-Cell Genomics and Expression Matrix Analysis

RULE ZERO — Check for pre-computed results FIRST

Before following any instruction below, scan the data folder for:

• *_executed.ipynb → read with tu run read_executed_notebook '{"data_folder":"<path>","search":"<keyword>"}' and cite its cell outputs as the authoritative answer
• Pre-computed result files (CSV/TSV with names like *results*, *deseq*, *enrich*, *stats*, *_simplified.csv) → read directly and report the requested value
• Canonical analysis scripts (analysis.R, run_*.py, find_*.R, *.Rmd) → execute as-is and read the output

Only follow this skill's re-analysis recipe below if none of the above exist. Re-running from raw data produces different numbers than the published answer and is much slower (often 5-10× turn count).

Comprehensive single-cell RNA-seq analysis and expression matrix processing using scanpy, anndata, scipy, and ToolUniverse.

LOOK UP, DON'T GUESS

When uncertain about any scientific fact, SEARCH databases first (PubMed, UniProt, ChEMBL, ClinVar, etc.) rather than reasoning from memory. A database-verified answer is always more reliable than a guess.

When to Use This Skill

Apply when users:

• Have scRNA-seq data (h5ad, 10X, CSV count matrices) and want analysis
• Need scRNA-seq quality control / QC gating: deciding cell filters by

mito % (pct_counts_mt), gene/UMI counts, doublets, ambient RNA, empty droplets

• Ask about cell type identification, clustering, or annotation
• Need differential expression analysis by cell type or condition
• Want gene-expression correlation analysis (e.g., gene length vs expression by cell type)
• Ask about PCA, UMAP, t-SNE for expression data
• Need Leiden/Louvain clustering on expression matrices
• Want statistical comparisons between cell types (t-test, ANOVA, fold change)
• Ask about marker genes, batch correction, trajectory, or cell-cell communication

NOT for (use other skills instead):

• Bulk RNA-seq DESeq2 only -> tooluniverse-rnaseq-deseq2
• Gene enrichment only -> tooluniverse-gene-enrichment
• VCF/variant analysis -> tooluniverse-variant-analysis
• Deep-learning integration/batch correction (scVI, scANVI, totalVI, PeakVI, MultiVI, DestVI, veloVI) -> tooluniverse-scvi-tools

Core Principles

1. Data-first - Load, inspect, validate before analysis
2. AnnData-centric - All data flows through anndata objects
3. Cell type awareness - Per-cell-type subsetting when needed
4. Statistical rigor - Normalization, multiple testing correction, effect sizes
5. Question-driven - Parse what the user is actually asking

Required Packages

import scanpy as sc, anndata as ad, pandas as pd, numpy as np
from scipy import stats
from scipy.cluster.hierarchy import linkage, fcluster
from sklearn.decomposition import PCA
from statsmodels.stats.multitest import multipletests
import gseapy as gp  # enrichment
import harmonypy     # batch correction (optional)

Install: pip install scanpy anndata leidenalg umap-learn harmonypy gseapy pandas numpy scipy scikit-learn statsmodels

Workflow Decision Tree

START: User question about scRNA-seq data
|
+-- FULL PIPELINE (raw counts -> annotated clusters)
|   Workflow: QC -> Normalize -> HVG -> PCA -> Cluster -> Annotate -> DE
|   See: references/scanpy_workflow.md
|
+-- DIFFERENTIAL EXPRESSION (per-cell-type comparison)
|   Most common pattern: per-cell-type DE
|   See: analysis_patterns.md "Pattern 1"
|
+-- CORRELATION ANALYSIS (gene property vs expression)
|   Pattern: Gene length vs expression correlation
|   See: analysis_patterns.md "Pattern 2"
|
+-- CLUSTERING & PCA (expression matrix analysis)
|   See: references/clustering_guide.md
|
+-- CELL COMMUNICATION (ligand-receptor interactions)
|   See: references/cell_communication.md
|
+-- TRAJECTORY ANALYSIS (pseudotime)
    See: references/trajectory_analysis.md

Data format handling:

• h5ad -> sc.read_h5ad()
• 10X -> sc.read_10x_mtx() or sc.read_10x_h5()
• CSV/TSV -> pd.read_csv() -> Convert to AnnData (check orientation!)

Data Loading

AnnData expects: cells/samples as rows (obs), genes as columns (var)

adata = sc.read_h5ad("data.h5ad")  # h5ad already oriented
# CSV/TSV: check orientation
df = pd.read_csv("counts.csv", index_col=0)
if df.shape[0] > df.shape[1] * 5:  # genes > samples by 5x => transpose
    df = df.T
adata = ad.AnnData(df)
# Load metadata
meta = pd.read_csv("metadata.csv", index_col=0)
common = adata.obs_names.intersection(meta.index)
adata = adata[common].copy()
for col in meta.columns:
    adata.obs[col] = meta.loc[common, col]

Quality Control and QC Gating (do this BEFORE downstream analysis)

QC gating decides which cells and genes are real before normalization, clustering, or DE. Skipping or rushing it propagates silently: doublets become fake "intermediate" states, ambient RNA smears markers across clusters, and empty droplets inflate cell counts. Never report a filtered cell count without the gates applied, and never report cutoffs you did not actually run.

HONEST EXECUTION: QC runs scanpy/AnnData via Bash/Python. If scanpy is not installed, do NOT fabricate metrics — print the install plan and stop: python scripts/scrna_qc.py --install-plan (exits 0, lists what is missing, suggests pip install scanpy anndata scrublet).

Per-cell QC metrics (what to compute)

vn = adata.var_names.str.upper()
adata.var['mt']   = vn.str.startswith('MT-')            # mitochondrial
adata.var['ribo'] = vn.str.startswith(('RPS', 'RPL'))   # ribosomal protein
adata.var['hb']   = vn.str.contains(r'^HB[^P]', regex=True)  # hemoglobin (RBC)
sc.pp.calculate_qc_metrics(
    adata, qc_vars=['mt', 'ribo', 'hb'],
    percent_top=None,   # REQUIRED for small gene panels (<500) — else IndexError
    log1p=True, inplace=True)

Key metrics in adata.obs: n_genes_by_counts, total_counts, pct_counts_mt, pct_counts_ribo, pct_counts_hb.

WHY each metric flags a problem (reason from biology, not magic numbers)

• High `pct_counts_mt` -> dying / stressed cell. A ruptured membrane lets

cytoplasmic mRNA leak out while mito transcripts stay trapped, enriching the captured RNA for mito. Cutoff is tissue-dependent (cardiomyocytes/hepatocytes are mito-rich at baseline — a blanket 10% would discard healthy cells).

• Low `n_genes_by_counts` / `total_counts` -> empty droplet or debris

(only ambient RNA captured; few genes, low depth).

• Very high counts/genes -> doublet (two transcriptomes ~double depth and

diversity) — but high count alone is weak; use a doublet caller (below).

• High `pct_counts_hb` -> RBC/blood contamination in solid tissue.

Choosing thresholds — distribution-aware (MAD), not hardcoded

Hardcoded cutoffs (mt<5%, n_genes<2500) are a starting point only; they break on mito-rich tissues and on shallow vs deep libraries. Prefer a MAD-based rule (robust to the very outliers you are removing): flag cells

nmads median-absolute-deviations from the median. Use nmads=5 on

log1p counts/genes (both tails), nmads=3 upper-only on pct_counts_mt, and pair mito with a biological ceiling so a uniformly degraded sample doesn't pass. Always visualize distributions first (violin + total_counts vs pct_counts_mt scatter — dying cells sit in the low-count/high-mito corner).

Run the helper (computes metrics + MAD gating, reports per-step removals):

python scripts/scrna_qc.py data.h5ad --doublets        # or --install-plan first

Doublets, ambient RNA, empty droplets (per-cell metrics miss these)

• Doublets: Scrublet (sc.pp.scrublet, scanpy >=1.10) or scDblFinder (R).

Run per sample before merging; flag-cluster-drop (doublets form bridge clusters). expected_doublet_rate ~0.8%/1,000 cells recovered (10x).

• Ambient RNA: cell-free "soup" mRNA in every droplet — a count-correction

step (SoupX / DecontX, R), NOT a cell filter. Per-cell QC cannot detect it; suspect it when markers look implausibly ubiquitous.

• Empty droplets: upstream of per-cell QC. CellRanger's filtered matrix

already applies an EmptyDrops-style call; with only the raw matrix, run EmptyDrops (DropletUtils, R) or a barcode-rank knee before per-cell QC.

QC Interpretation Table (metric -> concern -> action)

Full reasoning, MAD code, Scrublet/SoupX/EmptyDrops recipes, and order of operations: references/scrna_qc.md. Inline pipeline QC: references/scanpy_workflow.md Phase 2.

Complete Pipeline (Quick Reference)

import scanpy as sc
adata = sc.read_10x_h5("filtered_feature_bc_matrix.h5")
# QC
adata.var['mt'] = adata.var_names.str.startswith('MT-')
sc.pp.calculate_qc_metrics(adata, qc_vars=['mt'], inplace=True)
adata = adata[adata.obs['pct_counts_mt'] < 20].copy()
sc.pp.filter_cells(adata, min_genes=200)
sc.pp.filter_genes(adata, min_cells=3)
# Normalize + HVG + PCA
sc.pp.normalize_total(adata, target_sum=1e4)
sc.pp.log1p(adata)
adata.raw = adata.copy()
sc.pp.highly_variable_genes(adata, n_top_genes=2000)
sc.tl.pca(adata, n_comps=50)
# Cluster + UMAP
sc.pp.neighbors(adata, n_pcs=30)
sc.tl.leiden(adata, resolution=0.5)
sc.tl.umap(adata)
# Find markers + Annotate + Per-cell-type DE
sc.tl.rank_genes_groups(adata, groupby='leiden', method='wilcoxon')

Differential Expression Decision Tree

Single-Cell DE (many cells per condition):
  Use: sc.tl.rank_genes_groups(), methods: wilcoxon, t-test, logreg
  Best for: Per-cell-type DE, marker gene finding
 
Pseudo-Bulk DE (aggregate counts by sample):
  Use: R DESeq2 via `tu run run_deseq2_analysis` or Rscript (NOT pydeseq2 — gives different DEG counts)
  Best for: Sample-level comparisons with replicates
 
Statistical Tests Only:
  Use: scipy.stats (ttest_ind, f_oneway, pearsonr)
  Best for: Correlation, ANOVA, t-tests on summaries

Statistical Tests (Quick Reference)

from scipy import stats
from statsmodels.stats.multitest import multipletests
# Pearson/Spearman correlation
r, p = stats.pearsonr(gene_lengths, mean_expression)
# Welch's t-test
t_stat, p_val = stats.ttest_ind(group1, group2, equal_var=False)
# ANOVA
f_stat, p_val = stats.f_oneway(group1, group2, group3)
# Multiple testing correction (BH)
reject, pvals_adj, _, _ = multipletests(pvals, method='fdr_bh')

Batch Correction (Harmony)