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

---

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)

import harmonypy
sc.tl.pca(adata, n_comps=50)
ho = harmonypy.run_harmony(adata.obsm['X_pca'][:, :30], adata.obs, 'batch', random_state=0)
adata.obsm['X_pca_harmony'] = ho.Z_corr.T
sc.pp.neighbors(adata, use_rep='X_pca_harmony')
sc.tl.leiden(adata, resolution=0.5)
sc.tl.umap(adata)

---

ToolUniverse Integration

Data Discovery (before analysis)

• CxGDisc_search_datasets: Search CELLxGENE Discover for scRNA-seq datasets by disease, tissue, organism. Use broad disease terms (e.g., "breast cancer" not "triple-negative").
• GEO_search_rnaseq_datasets / geo_search_datasets: Search GEO for scRNA-seq studies
• NCBI_SRA_search_runs: Search SRA for sequencing runs (query="single cell RNA-seq [disease]")
• OmicsDI_search_datasets: Cross-repository dataset search

Cell Type Markers

• CellMarker_search_by_cell_type: Tissue-specific cell markers (use CellMarker_list_cell_types first — exact names required, e.g., "Regulatory T(Treg) cell" not "Regulatory T cell")
• CellMarker_search_cancer_markers: Cancer-context markers with experimental evidence
• CellMarker_search_by_gene: Reverse lookup — which cell types express a gene?
• HPA_search_genes_by_query: Cell-type marker gene search

Gene Annotation

• MyGene_query_genes / MyGene_batch_query: Gene ID conversion
• ensembl_lookup_gene: Ensembl gene details
• UniProt_get_function_by_accession: Protein function

Cell-Cell Communication

• OmniPath_get_ligand_receptor_interactions: L-R pairs (CellPhoneDB, CellChatDB)
• OmniPath_get_signaling_interactions: Downstream signaling
• OmniPath_get_complexes: Multi-subunit receptors

Enrichment (Post-DE)

• PANTHER_enrichment: GO enrichment (BP, MF, CC)
• STRING_functional_enrichment: Network-based enrichment
• ReactomeAnalysis_pathway_enrichment: Reactome pathways

Clinical Context (for tumor immunology)

• DGIdb_get_drug_gene_interactions: Drug interactions for immune checkpoint targets (genes=["CD274"] for PD-L1)
• civic_search_evidence_items: Clinical evidence for mutations/biomarkers
• TIMER2_immune_estimation: TCGA immune infiltration correlation
• search_clinical_trials: Clinical trial matching
• GTEx_get_expression_summary: Normal tissue baseline expression
• PubMed_search_articles: Literature context

---

Scanpy vs Seurat Equivalents

---

Reasoning Framework for Result Interpretation

Evidence Grading

Interpretation Guidance

• QC metric thresholds: Standard filters are nGenes > 200 (remove empty droplets), nGenes < 5000-6000 (remove doublets), pct_counts_mt < 20% (remove dying cells). These thresholds are tissue-dependent: immune cells tolerate stricter nGene filters; neurons may have higher mitochondrial content naturally. Always visualize distributions before setting cutoffs.
• Cluster resolution guidance: Leiden resolution 0.3-0.5 yields broad cell types (T cells, B cells, myeloid). Resolution 0.8-1.2 resolves subtypes (CD4 naive, CD4 memory, Treg). Resolution > 2.0 risks over-clustering (splitting biologically homogeneous populations). Validate by checking that each cluster has distinct marker genes.
• Marker gene confidence levels: A strong marker is highly specific (high pct_diff between cluster and rest) and highly expressed (high log2FC). Genes expressed in many clusters with small fold changes are poor markers. Cross-reference with known markers from CellMarker or HPA databases.
• Pseudo-bulk vs single-cell DE: For comparing conditions (treatment vs control), pseudo-bulk DE (aggregate by sample, then DESeq2) is more statistically valid than single-cell DE, which inflates significance due to non-independence of cells from the same sample.
• Batch effects: If samples cluster by batch rather than biology on UMAP, apply Harmony or other correction before biological interpretation.

Synthesis Questions

1. Do the identified clusters correspond to known cell types based on canonical markers, or do some clusters lack clear biological identity (potentially doublets or low-quality cells)? 2. At the chosen clustering resolution, are there clusters that merge when resolution is lowered, suggesting they may be a single cell type split by technical noise? 3. For differential expression between conditions, are the results consistent between single-cell and pseudo-bulk approaches, and do the top DE genes have known biological relevance? 4. Do QC-flagged cells (high mito, extreme gene counts) concentrate in specific clusters, and does removing them change the clustering structure? 5. If batch correction was applied, do post-correction clusters still maintain expected cell-type-specific marker expression?

---

Troubleshooting

---

Reference Documentation

Detailed Analysis Patterns: analysis_patterns.md (per-cell-type DE, correlation, PCA, ANOVA, cell communication)

Core Workflows:

• references/scrna_qc.md - scRNA-seq QC gating (mito%, doublets, ambient RNA, empty droplets, MAD thresholds)
• references/scanpy_workflow.md - Complete scanpy pipeline
• references/seurat_workflow.md - Seurat to Scanpy translation
• references/clustering_guide.md - Clustering methods
• references/marker_identification.md - Marker genes, annotation
• references/trajectory_analysis.md - Pseudotime
• references/cell_communication.md - OmniPath/CellPhoneDB workflow
• references/troubleshooting.md - Detailed error solutions

---

Analysis Conventions

DESeq2 library choice: match the authoritative script

If the data folder contains an authoritative script (run_*.py, analysis.R), use whichever DESeq2 library it uses (pydeseq2 or R DESeq2). The two libraries give slightly different DEG counts (~2-10% at the same thresholds), so matching matters. If no script exists, prefer R DESeq2 via the run_deseq2_analysis tool or Rscript:

tu run run_deseq2_analysis '{"operation":"deseq2","counts_file":"pseudo_bulk_counts.csv","metadata_file":"sample_meta.csv","design":"~ sex","contrast":"sex, M, F","lfc_shrinkage":true}'

Single-Cell Analysis Patterns

Detailed code patterns for common single-cell analysis types.

---

Pattern 1: Per-Cell-Type Differential Expression

Question: "Which immune cell type has the most DEGs after treatment?"

import scanpy as sc

adata = sc.read_h5ad("data.h5ad")
sc.pp.normalize_total(adata, target_sum=1e4)
sc.pp.log1p(adata)

cell_types = adata.obs['cell_type'].unique()
de_results = {}

for ct in cell_types:
    adata_ct = adata[adata.obs['cell_type'] == ct].copy()
    n_treat = (adata_ct.obs['condition'] == 'treatment').sum()
    n_ctrl = (adata_ct.obs['condition'] == 'control').sum()
    if n_treat < 3 or n_ctrl < 3:
        continue

    sc.tl.rank_genes_groups(adata_ct, groupby='condition',
                             groups=['treatment'], reference='control',
                             method='wilcoxon')
    df = sc.get.rank_genes_groups_df(adata_ct, group='treatment')
    sig = df[df['pvals_adj'] < 0.05]
    de_results[ct] = {'n_sig': len(sig), 'results': df}
    print(f"{ct}: {len(sig)} DEGs")

top_ct = max(de_results, key=lambda x: de_results[x]['n_sig'])
print(f"Answer: {top_ct} ({de_results[top_ct]['n_sig']} DEGs)")

---

Pattern 2: Gene Property vs Expression Correlation

Question: "What is the Pearson correlation between gene length and expression in CD4 T cells?"

import scanpy as sc
import pandas as pd
import numpy as np
from scipy import stats
from scipy.sparse import issparse

adata = sc.read_h5ad("data.h5ad")
gene_info = pd.read_csv("gene_info.tsv", sep='\t', index_col=0)
common = adata.var_names.intersection(gene_info.index)
adata.var['gene_length'] = gene_info.loc[common, 'gene_length'].reindex(adata.var_names)
adata.var['gene_type'] = gene_info.loc[common, 'gene_type'].reindex(adata.var_names)

mask = adata.var['gene_type'] == 'protein_coding'
adata_pc = adata[:, mask].copy()

cell_types = ['CD4 T cells', 'CD8 T cells', 'CD14 Monocytes']
for ct in cell_types:
    adata_ct = adata_pc[adata_pc.obs['cell_type'] == ct]
    X = adata_ct.X.toarray() if issparse(adata_ct.X) else adata_ct.X
    mean_expr = np.mean(X, axis=0)
    gene_lengths = adata_ct.var['gene_length'].values
    valid = ~np.isnan(gene_lengths) & ~np.isnan(mean_expr)
    r, p = stats.pearsonr(gene_lengths[valid], mean_expr[valid])
    print(f"{ct}: r = {r:.6f}, p = {p:.2e}, n = {valid.sum()} genes")

---

Pattern 3: PCA on Expression Matrix

Question: "What percentage of variance is explained by PC1 after log10 transform?"

import pandas as pd
import numpy as np
from sklearn.decomposition import PCA

df = pd.read_csv("expression.csv", index_col=0)
if df.shape[0] > df.shape[1] * 5:
    df = df.T  # Genes were rows, transpose

X = np.log10(df.values + 1)
n_components = min(X.shape[0], X.shape[1])
pca = PCA(n_components=n_components)
pca.fit(X)

print(f"PC1: {pca.explained_variance_ratio_[0]*100:.2f}% variance")
print(f"PC1+PC2: {sum(pca.explained_variance_ratio_[:2])*100:.2f}%")
print(f"Top 10 PCs: {sum(pca.explained_variance_ratio_[:10])*100:.2f}%")

---

Pattern 4: Statistical Comparison Between Cell Types

Question: "What is the t-statistic comparing LFCs between CD4/CD8 and other cell types?"

from scipy import stats
import numpy as np

# After running per-cell-type DE (Pattern 1):
cd4_lfc = de_results['CD4 T cells']['results']['log2FoldChange'].values
cd8_lfc = de_results['CD8 T cells']['results']['log2FoldChange'].values
cd4_cd8_lfc = np.concatenate([cd4_lfc, cd8_lfc])

other_lfc = []
for ct in ['CD14 Monocytes', 'NK cells', 'B cells']:
    other_lfc.append(de_results[ct]['results']['log2FoldChange'].values)
other_lfc = np.concatenate(other_lfc)

t_stat, p_val = stats.ttest_ind(cd4_cd8_lfc, other_lfc, equal_var=False)
print(f"t-statistic: {t_stat:.4f}")
print(f"p-value: {p_val:.4e}")

---

Pattern 5: ANOVA Across Cell Types

Question: "What is the F-statistic for miRNA expression across immune cell types?"

import pandas as pd
from scipy import stats

df = pd.read_csv("mirna_expr.csv", index_col=0)
meta = pd.read_csv("metadata.csv", index_col=0)

meta_filtered = meta[meta['cell_type'] != 'PBMC']
df_filtered = df[meta_filtered.index]

cell_types = meta_filtered['cell_type'].unique()
groups = {}
for ct in cell_types:
    samples = meta_filtered[meta_filtered['cell_type'] == ct].index
    groups[ct] = df_filtered[samples].values.flatten()

f_stat, p_val = stats.f_oneway(*groups.values())
print(f"F-statistic: {f_stat:.4f}")
print(f"p-value: {p_val:.4e}")

---

Pattern 6: Cell-Cell Communication Analysis

Question: "Which ligand-receptor interactions are strongest between tumor and T cells?"

from tooluniverse import ToolUniverse
import pandas as pd

tu = ToolUniverse()
tu.load_tools()

# Step 1: Get ligand-receptor pairs from OmniPath
result = tu.run_tool(
    "OmniPath_get_ligand_receptor_interactions",
    databases="CellPhoneDB,CellChatDB"
)
lr_pairs = pd.DataFrame(result['data']['interactions'])

# Step 2: Filter to expressed pairs
expressed_lr = lr_pairs[
    lr_pairs['source_genesymbol'].isin(adata.var_names) &
    lr_pairs['target_genesymbol'].isin(adata.var_names)
]

# Step 3: Score communication between cell types
communication_scores = score_cell_communication(
    adata, expressed_lr, cell_type_col='cell_type'
)

# Step 4: Filter to tumor-T cell interactions
tumor_tcell = communication_scores[
    ((communication_scores['sender'] == 'Tumor') &
     (communication_scores['receiver'].str.contains('T cell'))) |
    ((communication_scores['receiver'] == 'Tumor') &
     (communication_scores['sender'].str.contains('T cell')))
]

top_interactions = tumor_tcell.nlargest(20, 'score')
print(top_interactions[['sender', 'receiver', 'ligand', 'receptor', 'score']])

See: references/cell_communication.md for complete helper functions and scoring workflow.

---

Report Generation

Always extract the specific answer to the user's question:

report = f"""
# Analysis Results

## Per-Cell-Type Differential Expression

| Cell Type | Significant DEGs (padj < 0.05) |
|-----------|-------------------------------|
{chr(10).join([f"| {ct} | {res['n_sig']} |" for ct, res in de_results.items()])}

## Answer

**{top_ct}** has the highest number of significantly differentially expressed
genes with **{de_results[top_ct]['n_sig']} DEGs** (Wilcoxon test, BH-corrected
p < 0.05).
"""

Cell-Cell Communication Analysis

Complete guide for analyzing ligand-receptor interactions and cell-cell communication using OmniPath database (integrates CellPhoneDB, CellChatDB, and 100+ other databases) via ToolUniverse.

---

Overview

Cell-cell communication analysis identifies which cell types are signaling to each other through ligand-receptor (L-R) pairs. This is critical for understanding:

• Immune cell interactions (T cell exhaustion, activation)
• Tumor-immune communication (checkpoint blockade targets)
• Development and differentiation (niche signals)
• Tissue homeostasis (stromal-epithelial crosstalk)

Data sources: OmniPath integrates CellPhoneDB, CellChatDB, ICELLNET, Kirouac2010, Ramilowski2015, and 100+ other curated databases.

---

Workflow Overview

1. Get L-R Pairs from OmniPath
   └─ Query databases (CellPhoneDB, CellChatDB)

2. Filter to Expressed Pairs
   └─ Check genes present in dataset
   └─ Filter by expression thresholds

3. Score Cell-Cell Communication
   └─ Sender-receiver matrix
   └─ Calculate communication scores

4. Identify Top Interactions
   └─ Rank by score
   └─ Filter by biology (e.g., tumor-immune)

5. Trace Signaling Cascades
   └─ Get downstream targets
   └─ Identify transcription factors

6. Validate and Report
   └─ Cross-check with literature
   └─ Generate communication network

---

Step 1: Get Ligand-Receptor Pairs from OmniPath

from tooluniverse import ToolUniverse
import pandas as pd

tu = ToolUniverse()
tu.load_tools()

def get_ligand_receptor_pairs(proteins=None, databases=None):
    """Get ligand-receptor pairs from OmniPath.

    Args:
        proteins: Comma-separated protein names (None = all)
        databases: Comma-separated database names (None = all)
                  Options: CellPhoneDB, CellChatDB, etc.

    Returns:
        DataFrame with L-R interactions
    """
    result = tu.run_tool(
        "OmniPath_get_ligand_receptor_interactions",
        proteins=proteins,
        databases=databases
    )

    if result['metadata']['success']:
        interactions = result['data']['interactions']
        df = pd.DataFrame(interactions)
        return df
    else:
        print(f"Error: {result.get('error', 'Unknown')}")
        return pd.DataFrame()

# Example: Get all CellPhoneDB pairs
lr_pairs = get_ligand_receptor_pairs(databases="CellPhoneDB")
print(f"Found {len(lr_pairs)} L-R pairs from CellPhoneDB")

# Example: Get specific immune checkpoints
immune_checkpoints = "CD274,PDCD1,CTLA4,CD80,CD86,HAVCR2,TIGIT"
checkpoint_lr = get_ligand_receptor_pairs(proteins=immune_checkpoints)
print(checkpoint_lr[['source_genesymbol', 'target_genesymbol', 'curation_effort']])

Important columns:

• source_genesymbol: Ligand gene name
• target_genesymbol: Receptor gene name
• is_directed: True for directed interactions
• sources: Database sources (comma-separated)
• references: PubMed IDs
• curation_effort: Number of supporting evidences

---

Step 2: Filter to Expressed L-R Pairs

def filter_expressed_lr_pairs(adata, lr_pairs, min_frac=0.1, min_mean=0.1):
    """Filter L-R pairs to only those expressed in the dataset.

    Args:
        adata: AnnData object with normalized expression
        lr_pairs: DataFrame from get_ligand_receptor_pairs()
        min_frac: Minimum fraction of cells expressing (default 0.1 = 10%)
        min_mean: Minimum mean expression level (default 0.1)

    Returns:
        DataFrame with expressed L-R pairs
    """
    from scipy.sparse import issparse

    # Get genes in dataset
    genes_in_data = set(adata.var_names)

    # Filter to genes in data
    expressed_lr = lr_pairs[
        lr_pairs['source_genesymbol'].isin(genes_in_data) &
        lr_pairs['target_genesymbol'].isin(genes_in_data)
    ].copy()

    # Calculate expression statistics
    def gene_stats(gene):
        if gene not in genes_in_data:
            return 0, 0
        X = adata[:, gene].X
        if issparse(X):
            X = X.toarray()
        mean_expr = X.mean()
        frac_expr = (X > 0).mean()
        return mean_expr, frac_expr

    expressed_lr['ligand_mean'] = expressed_lr['source_genesymbol'].apply(
        lambda g: gene_stats(g)[0]
    )
    expressed_lr['ligand_frac'] = expressed_lr['source_genesymbol'].apply(
        lambda g: gene_stats(g)[1]
    )
    expressed_lr['receptor_mean'] = expressed_lr['target_genesymbol'].apply(
        lambda g: gene_stats(g)[0]
    )
    expressed_lr['receptor_frac'] = expressed_lr['target_genesymbol'].apply(
        lambda g: gene_stats(g)[1]
    )

    # Filter by thresholds
    expressed_lr = expressed_lr[
        (expressed_lr['ligand_mean'] >= min_mean) &
        (expressed_lr['receptor_mean'] >= min_mean) &
        (expressed_lr['ligand_frac'] >= min_frac) &
        (expressed_lr['receptor_frac'] >= min_frac)
    ]

    return expressed_lr

# Example
expressed_lr = filter_expressed_lr_pairs(adata, lr_pairs, min_frac=0.05, min_mean=0.05)
print(f"Expressed: {len(expressed_lr)}/{len(lr_pairs)} L-R pairs ({100*len(expressed_lr)/len(lr_pairs):.1f}%)")

---

Step 3: Score Cell-Cell Communication

def score_cell_communication(adata, lr_pairs, cell_type_col='cell_type',
                             method='mean_product'):
    """Score cell-cell communication for each cell type pair.

    Args:
        adata: AnnData with cell type annotations
        lr_pairs: Expressed L-R pairs from filter_expressed_lr_pairs()
        cell_type_col: Column in adata.obs with cell type labels
        method: 'mean_product' or 'fraction_product'

    Returns:
        DataFrame with columns: sender, receiver, ligand, receptor, score
    """
    import numpy as np
    from scipy.sparse import issparse

    cell_types = adata.obs[cell_type_col].unique()
    results = []

    for _, row in lr_pairs.iterrows():
        ligand = row['source_genesymbol']
        receptor = row['target_genesymbol']

        if ligand not in adata.var_names or receptor not in adata.var_names:
            continue

        for sender_ct in cell_types:
            for receiver_ct in cell_types:
                # Sender cells
                sender_mask = adata.obs[cell_type_col] == sender_ct
                ligand_expr = adata[sender_mask, ligand].X
                if issparse(ligand_expr):
                    ligand_expr = ligand_expr.toarray().flatten()

                # Receiver cells
                receiver_mask = adata.obs[cell_type_col] == receiver_ct
                receptor_expr = adata[receiver_mask, receptor].X
                if issparse(receptor_expr):
                    receptor_expr = receptor_expr.toarray().flatten()

                # Calculate score
                if method == 'mean_product':
                    ligand_mean = np.mean(ligand_expr)
                    receptor_mean = np.mean(receptor_expr)
                    score = ligand_mean * receptor_mean

                elif method == 'fraction_product':
                    ligand_frac = np.mean(ligand_expr > 0)
                    receptor_frac = np.mean(receptor_expr > 0)
                    ligand_mean = np.mean(ligand_expr[ligand_expr > 0]) if ligand_frac > 0 else 0
                    receptor_mean = np.mean(receptor_expr[receptor_expr > 0]) if receptor_frac > 0 else 0
                    score = ligand_frac * receptor_frac * ligand_mean * receptor_mean

                if score > 0:
                    results.append({
                        'sender': sender_ct,
                        'receiver': receiver_ct,
                        'ligand': ligand,
                        'receptor': receptor,
                        'ligand_mean': ligand_mean,
                        'receptor_mean': receptor_mean,
                        'score': score,
                        'curation_effort': row.get('curation_effort', 0),
                        'databases': row.get('sources', 'Unknown')
                    })

    return pd.DataFrame(results)

# Example
communication_scores = score_cell_communication(
    adata, expressed_lr, cell_type_col='cell_type'
)

# Top 20 interactions
top_20 = communication_scores.nlargest(20, 'score')
print("\nTop 20 cell-cell interactions:")
print(top_20[['sender', 'receiver', 'ligand', 'receptor', 'score']])

---

Step 4: Identify Top Interactions

Filter by Cell Type Pair

# Example: Tumor → T cell interactions
tumor_to_tcell = communication_scores[
    (communication_scores['sender'] == 'Tumor') &
    (communication_scores['receiver'].str.contains('T cell|CD4|CD8'))
]

print(f"\nTumor → T cell interactions: {len(tumor_to_tcell)}")
print(tumor_to_tcell.nlargest(10, 'score'))

Filter by Pathway

# Immune checkpoints
checkpoints = ['CD274', 'PDCD1', 'CTLA4', 'CD80', 'CD86', 'HAVCR2', 'TIGIT']
checkpoint_interactions = communication_scores[
    communication_scores['ligand'].isin(checkpoints) |
    communication_scores['receptor'].isin(checkpoints)
].sort_values('score', ascending=False)

print(f"\nCheckpoint interactions: {len(checkpoint_interactions)}")
print(checkpoint_interactions.head(10))

---

Step 5: Trace Downstream Signaling

def get_downstream_signaling(receptor_gene):
    """Get downstream signaling from a receptor.

    Args:
        receptor_gene: Receptor gene symbol

    Returns:
        DataFrame with signaling interactions
    """
    result = tu.run_tool(
        "OmniPath_get_signaling_interactions",
        proteins=receptor_gene,
        is_directed=True
    )

    if result['metadata']['success']:
        interactions = result['data']['interactions']
        df = pd.DataFrame(interactions)

        # Filter to receptor as source
        df = df[df['source_genesymbol'] == receptor_gene]

        return df[[
            'source_genesymbol', 'target_genesymbol',
            'is_stimulation', 'is_inhibition',
            'sources', 'references'
        ]]

    return pd.DataFrame()

# Example: PDCD1 (PD-1) signaling
pdcd1_signaling = get_downstream_signaling('PDCD1')
print(f"\nPDCD1 signals to {len(pdcd1_signaling)} targets")
print(pdcd1_signaling.head(10))

# Find transcription factors
tfs = pdcd1_signaling[pdcd1_signaling['target_genesymbol'].str.contains('NFAT|FOS|JUN|STAT')]
print(f"\nTranscription factors: {list(tfs['target_genesymbol'])}")

---

Step 6: Handle Protein Complexes

Some receptors are multi-subunit complexes (e.g., TGF-beta receptors):

def check_complex_expression(adata, complex_name, cell_type=None):
    """Check if all subunits of a protein complex are expressed.

    Args:
        adata: AnnData object
        complex_name: Complex name (e.g., "TGFBR2")
        cell_type: Optional cell type to subset

    Returns:
        List of dicts with complex info
    """
    import numpy as np

    # Get complex composition from OmniPath
    result = tu.run_tool("OmniPath_get_complexes", proteins=complex_name)

    if not result['metadata']['success'] or not result['data']['complexes']:
        return {'complex_found': False}

    complexes = result['data']['complexes']

    # Subset to cell type
    if cell_type:
        adata_subset = adata[adata.obs.cell_type == cell_type, :]
    else:
        adata_subset = adata

    results = []
    for complex_info in complexes:
        components = complex_info.get('components_genesymbols', '').split('_')

        # Check expression of each component
        component_expr = {}
        for comp in components:
            if comp in adata_subset.var_names:
                X = adata_subset[:, comp].X
                if issparse(X):
                    X = X.toarray().flatten()
                mean_expr = np.mean(X)
                frac_expr = np.mean(X > 0)
                component_expr[comp] = {'mean': mean_expr, 'fraction': frac_expr}
            else:
                component_expr[comp] = {'mean': 0, 'fraction': 0}

        # Complex score = minimum of subunit expressions
        min_mean = min([v['mean'] for v in component_expr.values()])
        min_frac = min([v['fraction'] for v in component_expr.values()])

        results.append({
            'complex_name': complex_info.get('name', 'Unknown'),
            'components': components,
            'component_expression': component_expr,
            'complex_score': min_mean * min_frac,
            'all_subunits_expressed': all([v['fraction'] > 0.1 for v in component_expr.values()])
        })

    return results

# Example: TGF-beta receptor
tgfb_receptor = check_complex_expression(adata, "TGFBR2", cell_type="Fibroblast")
for comp_result in tgfb_receptor:
    print(f"\nComplex: {comp_result['complex_name']}")
    print(f"Components: {comp_result['components']}")
    print(f"All expressed: {comp_result['all_subunits_expressed']}")
    print(f"Score: {comp_result['complex_score']:.4f}")

---

Step 7: Generate Communication Report

def generate_communication_report(adata, cell_type_col='cell_type',
                                 databases="CellPhoneDB,CellChatDB",
                                 min_score=0.01):
    """Generate complete cell-cell communication report."""

    report = []
    report.append("# Cell-Cell Communication Analysis Report\n")

    # Step 1: Get L-R pairs
    report.append("## 1. Ligand-Receptor Database Query")
    lr_pairs = get_ligand_receptor_pairs(databases=databases)
    report.append(f"- Total L-R pairs: {len(lr_pairs)}")
    report.append(f"- Databases: {databases}\n")

    # Step 2: Filter expressed
    report.append("## 2. Expressed Ligand-Receptor Pairs")
    expressed_lr = filter_expressed_lr_pairs(adata, lr_pairs, min_frac=0.05, min_mean=0.05)
    report.append(f"- Expressed: {len(expressed_lr)}/{len(lr_pairs)} ({100*len(expressed_lr)/len(lr_pairs):.1f}%)")
    report.append(f"- Thresholds: >5% cells, mean >0.05\n")

    # Step 3: Score communication
    report.append("## 3. Cell-Cell Communication Scores")
    communication_scores = score_cell_communication(adata, expressed_lr, cell_type_col=cell_type_col)
    communication_scores = communication_scores[communication_scores['score'] >= min_score]
    report.append(f"- Total interactions: {len(communication_scores)}")
    report.append(f"- Min score: {min_score}\n")

    # Top 20 table
    report.append("### Top 20 Interactions")
    report.append("| Sender | Receiver | Ligand | Receptor | Score | Curation |")
    report.append("|--------|----------|--------|----------|-------|----------|")
    for _, row in communication_scores.nlargest(20, 'score').iterrows():
        report.append(
            f"| {row['sender']} | {row['receiver']} | {row['ligand']} | "
            f"{row['receptor']} | {row['score']:.4f} | {row['curation_effort']} |"
        )
    report.append("")

    # Communication by cell type
    report.append("## 4. Communication Summary by Cell Type")
    sender_counts = communication_scores.groupby('sender').size().sort_values(ascending=False)
    receiver_counts = communication_scores.groupby('receiver').size().sort_values(ascending=False)

    report.append("\n### Top Sender Cell Types")
    report.append("| Cell Type | Outgoing Interactions |")
    report.append("|-----------|----------------------|")
    for ct, count in sender_counts.head(10).items():
        report.append(f"| {ct} | {count} |")

    report.append("\n### Top Receiver Cell Types")
    report.append("| Cell Type | Incoming Interactions |")
    report.append("|-----------|----------------------|")
    for ct, count in receiver_counts.head(10).items():
        report.append(f"| {ct} | {count} |")

    return "\n".join(report)

# Generate and save report
report = generate_communication_report(adata, cell_type_col='cell_type')
print(report)

with open('cell_communication_report.md', 'w') as f:
    f.write(report)

---

Example: Tumor-Immune Cell Communication

Complete workflow for T cell exhaustion analysis:

# Step 1: Get immune checkpoint L-R pairs
checkpoint_proteins = "CD274,PDCD1,CTLA4,CD80,CD86,HAVCR2,TIGIT,CD96,NECTIN2,LAG3"
checkpoint_lr = get_ligand_receptor_pairs(proteins=checkpoint_proteins)

# Step 2: Filter to expressed
expressed_checkpoints = filter_expressed_lr_pairs(adata, checkpoint_lr, min_frac=0.05)

# Step 3: Score communication
communication_scores = score_cell_communication(adata, expressed_checkpoints)

# Step 4: Tumor-T cell interactions
tumor_tcell = communication_scores[
    ((communication_scores['sender'] == 'Tumor') &
     (communication_scores['receiver'].str.contains('T cell|CD4|CD8'))) |
    ((communication_scores['receiver'] == 'Tumor') &
     (communication_scores['sender'].str.contains('T cell|CD4|CD8')))
]

# Step 5: Find exhaustion signals
exhaustion_pairs = tumor_tcell[
    tumor_tcell['ligand'].isin(['CD274', 'HAVCR2']) |  # PD-L1, TIM-3
    tumor_tcell['receptor'].isin(['PDCD1', 'HAVCR2', 'CTLA4', 'LAG3'])
].sort_values('score', ascending=False)

print("\nTop tumor-T cell exhaustion signals:")
print(exhaustion_pairs[['sender', 'receiver', 'ligand', 'receptor', 'score']])

# Step 6: Expression levels
print(f"\nCD274 (PD-L1) in tumor: {adata[adata.obs.cell_type=='Tumor', 'CD274'].X.mean():.3f}")
print(f"PDCD1 (PD-1) in T cells: {adata[adata.obs.cell_type.str.contains('T cell'), 'PDCD1'].X.mean():.3f}")

# Step 7: Downstream signaling
pdcd1_cascade = get_downstream_signaling('PDCD1')
print(f"\nPDCD1 downstream targets: {len(pdcd1_cascade)}")

---

Visualization

def plot_communication_network(communication_scores, top_n=50, min_score=0.01):
    """Plot cell-cell communication network."""
    import matplotlib.pyplot as plt
    import networkx as nx

    # Filter
    comm_filtered = communication_scores[communication_scores['score'] >= min_score]
    comm_top = comm_filtered.nlargest(top_n, 'score')

    # Build network
    G = nx.DiGraph()
    for _, row in comm_top.iterrows():
        edge_label = f"{row['ligand']}→{row['receptor']}"
        G.add_edge(row['sender'], row['receiver'],
                  weight=row['score'], label=edge_label)

    # Plot
    plt.figure(figsize=(12, 10))
    pos = nx.spring_layout(G, k=2, iterations=50)

    # Nodes
    nx.draw_networkx_nodes(G, pos, node_size=3000, node_color='lightblue', alpha=0.9)
    nx.draw_networkx_labels(G, pos, font_size=10, font_weight='bold')

    # Edges
    edges = G.edges()
    weights = [G[u][v]['weight'] for u, v in edges]
    max_weight = max(weights)
    widths = [5 * w / max_weight for w in weights]

    nx.draw_networkx_edges(G, pos, width=widths, alpha=0.6,
                          edge_color='gray', arrows=True, arrowsize=20)

    plt.title(f"Cell-Cell Communication Network (Top {top_n})", fontsize=14)
    plt.axis('off')
    plt.tight_layout()

    return plt

# Plot
plot = plot_communication_network(communication_scores, top_n=30)
plot.savefig('communication_network.png', dpi=300, bbox_inches='tight')

---

Tips and Best Practices

1. Expression thresholds: Balance sensitivity vs specificity

• Stringent: min_frac=0.1, min_mean=0.1 (fewer false positives)
• Permissive: min_frac=0.05, min_mean=0.05 (capture rare interactions)

2. Communication score: Mean product method is simple and interpretable

• Fraction product accounts for expression breadth
• Can add log-transform for very skewed distributions

3. Database selection:

• CellPhoneDB: Well-curated, human-focused
• CellChatDB: Broader coverage, mouse + human
• Use both for comprehensive analysis

4. Validation: Always cross-validate top hits:

• Check expression in UMAP/violin plots
• Verify with literature (PubMed IDs in references column)
• Compare across replicates/cohorts

5. Protein complexes: Check multi-subunit receptors

• Use check_complex_expression() for receptors like TGFBR1/2
• All subunits must be expressed for functional complex

6. Statistical testing: For rigorous analysis:

• Permutation test (shuffle cell labels)
• Compare to random cell type assignments
• Correct for multiple testing

---

See Also

• scanpy_workflow.md - Load and normalize data before communication analysis
• marker_identification.md - Cell type annotation for communication analysis
• troubleshooting.md - Common OmniPath API issues

Clustering Methods for Single-Cell Data

Guide to different clustering methods: Leiden, Louvain, hierarchical, and bootstrap consensus clustering.

---

Leiden Clustering (Recommended)

Best all-around method for single-cell data.

import scanpy as sc

# Build neighbor graph
sc.pp.neighbors(adata, n_neighbors=15, n_pcs=30)

# Run Leiden
sc.tl.leiden(adata, resolution=0.5, random_state=0)

print(f"Clusters: {adata.obs['leiden'].nunique()}")

Resolution parameter:

• Higher = More clusters
• Typical range: 0.3 - 1.5
• Start with 0.5, adjust based on biology

---

Louvain Clustering

Alternative to Leiden (older algorithm).

sc.tl.louvain(adata, resolution=0.5, random_state=0)

Leiden vs Louvain:

• Leiden: Better optimization, guaranteed connected communities
• Louvain: Faster, may produce disconnected communities
• For publication: Use Leiden

---

Hierarchical Clustering

For expression matrices (not single-cell level).

from scipy.cluster.hierarchy import linkage, fcluster, dendrogram
from scipy.spatial.distance import pdist
import numpy as np

def hierarchical_clustering(expression_df, n_clusters=3, method='ward', metric='euclidean'):
    """Hierarchical clustering on expression matrix.

    Args:
        expression_df: DataFrame (samples as rows, genes as columns)
        n_clusters: Number of clusters
        method: 'ward', 'complete', 'average', 'single'
        metric: Distance metric

    Returns:
        dict with labels, linkage_matrix
    """
    # Compute linkage
    if method == 'ward':
        Z = linkage(expression_df.values, method='ward')
    else:
        dist = pdist(expression_df.values, metric=metric)
        Z = linkage(dist, method=method)

    # Cut tree
    labels = fcluster(Z, t=n_clusters, criterion='maxclust')

    # Cluster sizes
    unique, counts = np.unique(labels, return_counts=True)
    for c, n in zip(unique, counts):
        print(f"  Cluster {c}: {n} samples")

    return {
        'labels': labels,
        'linkage_matrix': Z,
        'n_clusters': n_clusters
    }

# Example
result = hierarchical_clustering(expr_df, n_clusters=3, method='ward')

Linkage methods:

• ward: Minimizes within-cluster variance (best for most data)
• complete: Maximum distance between clusters (compact clusters)
• average: Average distance (balanced)
• single: Minimum distance (can create chains)

---

Bootstrap Consensus Clustering

Robust clustering with logistic regression prediction.

from sklearn.linear_model import LogisticRegression
from sklearn.preprocessing import StandardScaler

def bootstrap_consensus_clustering(expression_df, n_clusters=3, n_iterations=50,
                                    train_fraction=0.7, random_state=42):
    """Bootstrap consensus clustering.

    Args:
        expression_df: DataFrame (samples as rows, genes as columns)
        n_clusters: Number of clusters
        n_iterations: Bootstrap iterations
        train_fraction: Fraction for training
        random_state: Random seed

    Returns:
        dict with labels, consensus_matrix, consistent_count
    """
    np.random.seed(random_state)
    n_samples = len(expression_df)
    n_train = int(n_samples * train_fraction)

    # Consensus matrices
    train_consensus = np.zeros((n_samples, n_samples))
    test_consensus = np.zeros((n_samples, n_samples))
    train_count = np.zeros((n_samples, n_samples))
    test_count = np.zeros((n_samples, n_samples))

    for i in range(n_iterations):
        # Random split
        indices = np.random.permutation(n_samples)
        train_idx = indices[:n_train]
        test_idx = indices[n_train:]

        # Cluster training data
        train_data = expression_df.iloc[train_idx]
        Z = linkage(train_data.values, method='ward')
        train_labels = fcluster(Z, t=n_clusters, criterion='maxclust')

        # Update train consensus
        for a_i, a in enumerate(train_idx):
            for b_i, b in enumerate(train_idx):
                train_count[a, b] += 1
                if train_labels[a_i] == train_labels[b_i]:
                    train_consensus[a, b] += 1

        # Predict test labels
        scaler = StandardScaler()
        X_train = scaler.fit_transform(train_data.values)
        X_test = scaler.transform(expression_df.iloc[test_idx].values)

        lr = LogisticRegression(max_iter=1000, random_state=random_state)
        lr.fit(X_train, train_labels)
        test_labels = lr.predict(X_test)

        # Update test consensus
        for a_i, a in enumerate(test_idx):
            for b_i, b in enumerate(test_idx):
                test_count[a, b] += 1
                if test_labels[a_i] == test_labels[b_i]:
                    test_consensus[a, b] += 1

    # Normalize
    with np.errstate(divide='ignore', invalid='ignore'):
        train_consensus_norm = np.where(train_count > 0, train_consensus / train_count, 0)
        test_consensus_norm = np.where(test_count > 0, test_consensus / test_count, 0)

    # Final clustering
    combined_consensus = (train_consensus_norm + test_consensus_norm) / 2
    np.fill_diagonal(combined_consensus, 1.0)

    from scipy.spatial.distance import squareform
    dist = 1 - combined_consensus
    np.fill_diagonal(dist, 0)
    dist = np.maximum(dist, 0)
    dist = (dist + dist.T) / 2
    condensed = squareform(dist)
    Z_final = linkage(condensed, method='average')
    final_labels = fcluster(Z_final, t=n_clusters, criterion='maxclust')

    # Count consistent samples
    consistent_count = 0
    for s in range(n_samples):
        cluster = final_labels[s]
        same_cluster = np.where(final_labels == cluster)[0]
        same_cluster = same_cluster[same_cluster != s]

        if len(same_cluster) > 0:
            train_scores = [train_consensus_norm[s, j] for j in same_cluster if train_count[s, j] > 0]
            test_scores = [test_consensus_norm[s, j] for j in same_cluster if test_count[s, j] > 0]

            if train_scores and test_scores:
                if np.mean(train_scores) > 0.7 and np.mean(test_scores) > 0.7:
                    consistent_count += 1

    print(f"Consistently classified: {consistent_count}/{n_samples}")

    return {
        'labels': final_labels,
        'train_consensus': train_consensus_norm,
        'test_consensus': test_consensus_norm,
        'combined_consensus': combined_consensus,
        'consistent_count': consistent_count
    }

# Example
result = bootstrap_consensus_clustering(expr_df, n_clusters=3, n_iterations=50)
print(f"Answer: {result['consistent_count']} samples consistently classified")

---

PCA for Clustering

Perform PCA on expression matrix.

from sklearn.decomposition import PCA
import pandas as pd
import numpy as np

def manual_pca(expression_df, log_transform='log10', pseudocount=1):
    """Run PCA with specific transforms.

    Args:
        expression_df: DataFrame (samples as rows, genes as columns)
        log_transform: 'log10', 'log2', 'log1p', or None
        pseudocount: Pseudocount for log

    Returns:
        dict with variance_ratio, pc_coords, loadings
    """
    X = expression_df.values.astype(float)

    # Transform
    if log_transform == 'log10':
        X = np.log10(X + pseudocount)
    elif log_transform == 'log2':
        X = np.log2(X + pseudocount)
    elif log_transform == 'log1p':
        X = np.log1p(X)

    # Run PCA
    n_components = min(X.shape[0], X.shape[1])
    pca = PCA(n_components=n_components)
    pc_coords = pca.fit_transform(X)

    # Results
    result = {
        'variance_ratio': pca.explained_variance_ratio_,
        'variance_explained': pca.explained_variance_,
        'pc_coords': pd.DataFrame(
            pc_coords,
            index=expression_df.index,
            columns=[f'PC{i+1}' for i in range(n_components)]
        ),
        'loadings': pd.DataFrame(
            pca.components_.T,
            index=expression_df.columns,
            columns=[f'PC{i+1}' for i in range(n_components)]
        ),
        'cumulative_variance': np.cumsum(pca.explained_variance_ratio_)
    }

    print(f"PC1: {result['variance_ratio'][0]*100:.2f}% variance")
    print(f"PC1-10: {result['cumulative_variance'][9]*100:.2f}% variance")

    return result

# Example
pca_result = manual_pca(expr_df, log_transform='log10', pseudocount=1)

---

Choosing the Right Method

---

Validation

Silhouette Score

from sklearn.metrics import silhouette_score

# After clustering
silhouette_avg = silhouette_score(X, labels)
print(f"Silhouette score: {silhouette_avg:.3f}")
# Range: -1 to 1. >0.5 is good.

Cluster Stability

# Run clustering multiple times with different random seeds
from collections import Counter

all_labels = []
for seed in range(10):
    sc.tl.leiden(adata, resolution=0.5, random_state=seed, key_added=f'leiden_{seed}')
    all_labels.append(adata.obs[f'leiden_{seed}'])

# Check consistency
# (Implementation depends on label alignment)

---

See Also

• scanpy_workflow.md - Prepare data for clustering
• marker_identification.md - Annotate clusters after clustering

Marker Gene Identification and Cell Type Annotation

Guide to finding marker genes and annotating cell types in single-cell data.

---

Find Marker Genes for Clusters

import scanpy as sc

# Run DE for all clusters
sc.tl.rank_genes_groups(
    adata,
    groupby='leiden',
    method='wilcoxon',  # or 't-test', 'logreg'
    n_genes=100,
    corr_method='benjamini-hochberg'
)

# Get results for cluster 0
markers_0 = sc.get.rank_genes_groups_df(adata, group='0')
print(markers_0.head(10))

Available methods:

• wilcoxon: Non-parametric, robust (default)
• t-test: Fast, parametric
• logreg: Logistic regression, good for classification

---

Cell Type Annotation Strategies

1. Known Marker Genes

# Define markers for cell types
marker_genes = {
    'T cells': ['CD3D', 'CD3E', 'CD8A', 'CD4'],
    'CD4 T cells': ['CD3D', 'CD4', 'IL7R'],
    'CD8 T cells': ['CD3D', 'CD8A', 'CD8B'],
    'B cells': ['CD19', 'MS4A1', 'CD79A'],
    'Monocytes': ['CD14', 'LYZ', 'S100A9'],
    'CD14 Monocytes': ['CD14', 'LYZ'],
    'CD16 Monocytes': ['FCGR3A', 'MS4A7'],
    'NK cells': ['NKG7', 'GNLY', 'KLRB1'],
    'Dendritic cells': ['FCER1A', 'CD1C'],
}

# Score clusters
from scipy.sparse import issparse
X = adata.X.toarray() if issparse(adata.X) else adata.X
expr_df = pd.DataFrame(X, index=adata.obs_names, columns=adata.var_names)

cluster_scores = {}
for ct, markers in marker_genes.items():
    available = [m for m in markers if m in adata.var_names]
    if available:
        scores = expr_df[available].mean(axis=1)
        cluster_scores[ct] = scores.groupby(adata.obs['leiden']).mean()

# Assign cell types
score_df = pd.DataFrame(cluster_scores)
assignments = score_df.idxmax(axis=1)
adata.obs['cell_type'] = adata.obs['leiden'].map(assignments)

2. ToolUniverse HPA Database

from tooluniverse import ToolUniverse
tu = ToolUniverse()
tu.load_tools()

# Search for cell-type markers
result = tu.tools.HPA_search_genes_by_query(query="T cell marker blood")
if result:
    print("HPA T cell markers:", [g.get('gene_name') for g in result[:10]])

# Get tissue-specific expression
result = tu.tools.HPA_get_rna_expression_in_specific_tissues(
    ensembl_id="ENSG00000167286",  # CD3D
    tissue_name="blood"
)

3. Automated Annotation (scanpy)

# Using marker gene scores
sc.tl.marker_gene_overlap(adata, marker_genes)

# Using cell type prediction (if reference available)
import scanpy.external as sce
# sce.tl.ingest(adata, adata_ref, obs='cell_type')

---

Visualize Marker Expression

# Dot plot
sc.pl.dotplot(adata, marker_genes, groupby='leiden')

# Stacked violin
sc.pl.stacked_violin(adata, marker_genes, groupby='leiden')

# UMAP with marker overlay
sc.pl.umap(adata, color=['leiden', 'CD3D', 'CD79A', 'CD14'])

---

Validate Annotations

Check Marker Specificity

# Calculate marker specificity
for ct in adata.obs['cell_type'].unique():
    mask = adata.obs['cell_type'] == ct
    markers_for_ct = marker_genes.get(ct, [])

    for marker in markers_for_ct:
        if marker in adata.var_names:
            expr_in_ct = adata[mask, marker].X.mean()
            expr_out_ct = adata[~mask, marker].X.mean()
            fold_enrichment = expr_in_ct / (expr_out_ct + 0.001)
            print(f"{ct} - {marker}: {fold_enrichment:.2f}x enriched")

Cross-Reference with Literature

# Get top expressed genes in cell type
ct_data = adata[adata.obs['cell_type'] == 'CD4 T cells']
X = ct_data.X.toarray() if issparse(ct_data.X) else ct_data.X
mean_expr = np.mean(X, axis=0)
top_genes = ct_data.var_names[np.argsort(mean_expr)[::-1][:20]]
print(f"Top genes in CD4 T cells: {list(top_genes)}")

---

Common Cell Type Markers

Immune Cells (PBMC)

immune_markers = {
    'T cells': ['CD3D', 'CD3E', 'CD3G'],
    'CD4 T cells': ['CD3D', 'CD4', 'IL7R', 'TCF7'],
    'CD8 T cells': ['CD3D', 'CD8A', 'CD8B', 'GZMK'],
    'Regulatory T cells': ['FOXP3', 'IL2RA', 'CTLA4'],
    'NK cells': ['NKG7', 'GNLY', 'KLRD1', 'KLRB1'],
    'B cells': ['CD19', 'MS4A1', 'CD79A', 'CD79B'],
    'Plasma cells': ['IGHG1', 'MZB1', 'SDC1', 'XBP1'],
    'Monocytes': ['CD14', 'LYZ', 'S100A9', 'S100A8'],
    'CD14 Monocytes': ['CD14', 'LYZ', 'FCN1'],
    'CD16 Monocytes': ['FCGR3A', 'MS4A7', 'CDKN1C'],
    'Dendritic cells': ['FCER1A', 'CD1C', 'CLEC10A'],
    'pDC': ['IL3RA', 'GZMB', 'SERPINF1', 'ITM2C'],
    'Platelets': ['PPBP', 'PF4', 'TUBB1'],
}

Tumor Microenvironment

tme_markers = {
    'Cancer cells': ['EPCAM', 'KRT8', 'KRT18', 'KRT19'],
    'CAFs': ['COL1A1', 'COL1A2', 'DCN', 'ACTA2'],
    'Endothelial': ['PECAM1', 'VWF', 'CDH5'],
    'TAMs': ['CD68', 'CD163', 'MSR1'],
    'Exhausted T cells': ['PDCD1', 'HAVCR2', 'LAG3', 'TIGIT'],
}

---

Tips

1. Multiple markers: Use 3-5 markers per cell type for robustness 2. Negative markers: Also check genes NOT expressed 3. Expression levels: Mean expression matters, not just presence 4. Subclustering: For ambiguous clusters, re-cluster at higher resolution 5. Manual curation: Always manually review top markers

---

See Also

• scanpy_workflow.md - Clustering before annotation
• cell_communication.md - Use annotated cell types for communication

Complete Scanpy Workflow

Complete reference for single-cell RNA-seq analysis using scanpy, from raw counts to annotated cell types.

---

Phase 1: Data Loading and Validation

Load h5ad Files

import scanpy as sc
adata = sc.read_h5ad("data.h5ad")
print(f"Shape: {adata.n_obs} cells x {adata.n_vars} genes")
print(f"Obs columns: {list(adata.obs.columns)}")
print(f"Var columns: {list(adata.var.columns)}")

Load 10X Files

# From directory
adata = sc.read_10x_mtx("filtered_gene_bc_matrices/hg19/")

# From HDF5
adata = sc.read_10x_h5("filtered_feature_bc_matrix.h5")

Load CSV/TSV and Convert to AnnData

import anndata as ad
import pandas as pd

df = pd.read_csv("counts.csv", index_col=0)

# Check orientation (genes vs cells)
if df.shape[0] > df.shape[1] * 5:
    print("Transposing: genes were rows")
    df = df.T

adata = ad.AnnData(df)

Attach Metadata

meta = pd.read_csv("metadata.csv", index_col=0, sep='\t')

# Align indices
common = adata.obs_names.intersection(meta.index)
adata = adata[common].copy()

for col in meta.columns:
    adata.obs[col] = meta.loc[common, col]

Attach Gene Annotations

gene_info = pd.read_csv("gene_info.tsv", sep='\t', index_col=0)

common_genes = adata.var_names.intersection(gene_info.index)

for col in ['gene_length', 'gene_type', 'chromosome']:
    if col in gene_info.columns:
        adata.var[col] = gene_info.loc[adata.var_names, col]

---

Phase 2: Quality Control

Calculate QC Metrics

# Identify mitochondrial genes
adata.var['mt'] = adata.var_names.str.startswith(('MT-', 'mt-'))

# Calculate metrics
sc.pp.calculate_qc_metrics(
    adata, qc_vars=['mt'],
    percent_top=None,
    log1p=False,
    inplace=True
)

# Available metrics:
# - total_counts: Total UMI counts per cell
# - n_genes_by_counts: Number of genes expressed per cell
# - pct_counts_mt: Percentage of counts in mitochondrial genes

Filter Cells

n_before = adata.n_obs

# Minimum genes per cell
sc.pp.filter_cells(adata, min_genes=200)

# Maximum mitochondrial percentage
adata = adata[adata.obs['pct_counts_mt'] < 20].copy()

# Optional: Remove doublets by gene count
# adata = adata[adata.obs['n_genes_by_counts'] < 5000].copy()

# Optional: Minimum UMI counts
# sc.pp.filter_cells(adata, min_counts=500)

n_after = adata.n_obs
print(f"Filtered: {n_before} → {n_after} cells ({n_before - n_after} removed)")

Filter Genes

# Minimum cells per gene (remove rare genes)
sc.pp.filter_genes(adata, min_cells=3)

print(f"After gene filtering: {adata.n_vars} genes")

Doublet Detection (Optional)

# Using scrublet via scanpy
sc.external.pp.scrublet(adata, expected_doublet_rate=0.06)

n_doublets = adata.obs['predicted_doublet'].sum()
print(f"Detected {n_doublets} doublets ({n_doublets/adata.n_obs*100:.1f}%)")

# Remove doublets
adata = adata[~adata.obs['predicted_doublet']].copy()

---

Phase 3: Normalization and Scaling

Store Raw Counts

# Important: Store raw counts before normalization
adata.raw = adata.copy()

Library-Size Normalization

# Normalize each cell to 10,000 total counts
sc.pp.normalize_total(adata, target_sum=1e4)

Log Transformation

# Natural log (log1p = log(x + 1))
sc.pp.log1p(adata)

Highly Variable Genes

# Find top variable genes
sc.pp.highly_variable_genes(
    adata,
    n_top_genes=2000,
    flavor='seurat_v3'  # Use 'seurat' if already log-transformed
)

print(f"Highly variable genes: {adata.var['highly_variable'].sum()}")

Scale Data

# Z-score scaling
sc.pp.scale(adata, max_value=10)

---

Phase 4: Dimensionality Reduction

PCA

# Run PCA on highly variable genes
sc.tl.pca(adata, n_comps=50, use_highly_variable=True)

# Variance explained
var_ratio = adata.uns['pca']['variance_ratio']
print(f"PC1: {var_ratio[0]*100:.2f}% variance")
print(f"Top 10 PCs: {sum(var_ratio[:10])*100:.2f}% variance")

# PC coordinates: adata.obsm['X_pca']
# PC loadings: adata.varm['PCs']

UMAP

# Compute neighbors first
sc.pp.neighbors(adata, n_neighbors=15, n_pcs=30)

# Run UMAP
sc.tl.umap(adata)

# UMAP coordinates: adata.obsm['X_umap']

t-SNE (Alternative)

sc.tl.tsne(adata, n_pcs=30)
# t-SNE coordinates: adata.obsm['X_tsne']

---

Phase 5: Clustering

Leiden Clustering (Recommended)

# Build neighbor graph (if not done for UMAP)
sc.pp.neighbors(adata, n_neighbors=15, n_pcs=30)

# Leiden clustering
sc.tl.leiden(adata, resolution=0.5, random_state=0)

n_clusters = adata.obs['leiden'].nunique()
print(f"Leiden clustering: {n_clusters} clusters")

Louvain Clustering (Alternative)

sc.tl.louvain(adata, resolution=0.5, random_state=0)

Resolution Parameter

• Higher resolution = More clusters
• Typical range: 0.3 - 1.5
• Start with 0.5, adjust based on biological expectations

---

Phase 6: Marker Gene Identification

Find Marker Genes for Each Cluster

# Run DE test for all clusters
sc.tl.rank_genes_groups(
    adata,
    groupby='leiden',
    method='wilcoxon',  # or 't-test', 'logreg'
    n_genes=100,
    corr_method='benjamini-hochberg'
)

# Get results for cluster 0
markers_0 = sc.get.rank_genes_groups_df(adata, group='0')
print(markers_0.head(10))

# Top marker genes per cluster
sc.pl.rank_genes_groups(adata, n_genes=5, sharey=False)

---

Phase 7: Cell Type Annotation

Marker-Based Annotation

# Known markers
marker_genes = {
    'T cells': ['CD3D', 'CD3E', 'CD8A', 'CD4'],
    'B cells': ['CD19', 'MS4A1', 'CD79A'],
    'Monocytes': ['CD14', 'LYZ', 'S100A9'],
    'NK cells': ['NKG7', 'GNLY', 'KLRB1'],
    'Dendritic cells': ['FCER1A', 'CD1C'],
}

# Score each cluster
from scipy.sparse import issparse
X = adata.X.toarray() if issparse(adata.X) else adata.X
expr_df = pd.DataFrame(X, index=adata.obs_names, columns=adata.var_names)

cluster_scores = {}
for ct, markers in marker_genes.items():
    available_markers = [m for m in markers if m in adata.var_names]
    if available_markers:
        scores = expr_df[available_markers].mean(axis=1)
        cluster_scores[ct] = scores.groupby(adata.obs['leiden']).mean()

# Assign cell types
score_df = pd.DataFrame(cluster_scores)
assignments = score_df.idxmax(axis=1)
adata.obs['cell_type'] = adata.obs['leiden'].map(assignments)

Use ToolUniverse for Marker Discovery

from tooluniverse import ToolUniverse
tu = ToolUniverse()
tu.load_tools()

# Search HPA for tissue-specific markers
result = tu.tools.HPA_search_genes_by_query(
    query="T cell marker blood"
)

---

Phase 8: Differential Expression Analysis

Per-Cell-Type DE Between Conditions

cell_types = adata.obs['cell_type'].unique()
de_results = {}

for ct in cell_types:
    # Subset to cell type
    adata_ct = adata[adata.obs['cell_type'] == ct].copy()

    # Check sufficient cells
    n_treat = (adata_ct.obs['condition'] == 'treatment').sum()
    n_ctrl = (adata_ct.obs['condition'] == 'control').sum()

    if n_treat < 3 or n_ctrl < 3:
        print(f"{ct}: Skipped (insufficient cells)")
        continue

    # Run DE
    sc.tl.rank_genes_groups(
        adata_ct,
        groupby='condition',
        groups=['treatment'],
        reference='control',
        method='wilcoxon',
        n_genes=adata_ct.n_vars
    )

    # Get results
    df = sc.get.rank_genes_groups_df(adata_ct, group='treatment')

    # Filter significant
    sig = df[(df['pvals_adj'] < 0.05) & (df['logfoldchanges'].abs() > 0.5)]

    de_results[ct] = {
        'all': df,
        'significant': sig,
        'n_sig': len(sig),
        'n_up': (sig['logfoldchanges'] > 0).sum(),
        'n_down': (sig['logfoldchanges'] < 0).sum(),
    }

    print(f"{ct}: {len(sig)} DEGs ({de_results[ct]['n_up']} up, {de_results[ct]['n_down']} down)")

---

Phase 9: Batch Correction with Harmony

import harmonypy

# After PCA
sc.tl.pca(adata, n_comps=50)

# Run Harmony
ho = harmonypy.run_harmony(
    adata.obsm['X_pca'][:, :30],  # Use first 30 PCs
    adata.obs,
    'batch',  # Batch column name
    random_state=0
)

# Store corrected PCs
adata.obsm['X_pca_harmony'] = ho.Z_corr.T

# Re-compute neighbors and cluster on corrected PCs
sc.pp.neighbors(adata, use_rep='X_pca_harmony', n_pcs=30)
sc.tl.leiden(adata, resolution=0.5)
sc.tl.umap(adata)

---

Complete Pipeline Example

import scanpy as sc

# 1. Load
adata = sc.read_10x_h5("data.h5")

# 2. 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)

# 3. Normalize
adata.raw = adata.copy()
sc.pp.normalize_total(adata, target_sum=1e4)
sc.pp.log1p(adata)

# 4. HVG + Scale
sc.pp.highly_variable_genes(adata, n_top_genes=2000)
sc.pp.scale(adata, max_value=10)

# 5. PCA
sc.tl.pca(adata, n_comps=50)

# 6. Cluster
sc.pp.neighbors(adata, n_pcs=30)
sc.tl.leiden(adata, resolution=0.5)
sc.tl.umap(adata)

# 7. Markers
sc.tl.rank_genes_groups(adata, groupby='leiden', method='wilcoxon')

# 8. Annotate (manual step)
# 9. DE analysis (per cell type, if conditions present)

# Save
adata.write_h5ad("processed.h5ad")

---

Tips and Best Practices

1. Always store raw counts before normalization (adata.raw = adata.copy()) 2. QC thresholds depend on dataset:

• min_genes: 200-500
• pct_counts_mt: 10-20%
• max_genes: 5000-7000 (doublet filter)

3. Highly variable genes: 2000-3000 for most datasets 4. PCA components: 30-50 sufficient for most analyses 5. Resolution tuning: Start with 0.5, increase for finer clusters 6. Batch correction: Use Harmony for multiple batches/samples 7. DE method: Wilcoxon (default) good for most cases; t-test faster 8. Statistical power: Need >= 3 cells per condition per cell type

---

See Also

• clustering_guide.md - Advanced clustering methods
• marker_identification.md - Cell type annotation strategies
• troubleshooting.md - Common errors and solutions

scRNA-seq Quality Control: Gating Before Downstream Analysis

The QC-gating step decides which cells (and genes) are real before any normalization, clustering, or DE. Bad gating propagates silently: doublets become fake "intermediate" cell states, ambient RNA makes every cluster express every marker, and empty droplets inflate cell counts. This doc covers the per-cell QC metrics, why each one flags a problem, how to choose thresholds from the data (not magic numbers), and the technical artifacts (doublets, ambient RNA, empty droplets) that per-cell metrics alone do not catch.

Honesty note: every command here runs scanpy/AnnData via Bash/Python.

If scanpy is not installed, do NOT fabricate numbers. Emit the install plan

(pip install scanpy scrublet) and run scripts/scrna_qc.py --install-plan

to confirm the environment, then stop and report what is missing.

---

1. The standard per-cell QC metrics

After sc.pp.calculate_qc_metrics(adata, qc_vars=['mt','ribo','hb'], ...), each cell (row of adata.obs) carries:

Gene flags are set on adata.var BEFORE calling calculate_qc_metrics:

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

---

2. WHY each metric flags a problem (the biology)

• High `pct_counts_mt` -> dying / stressed cell. When a cell's membrane

ruptures during dissociation, cytoplasmic mRNA leaks out but mitochondrial transcripts stay trapped inside mitochondria. The surviving captured RNA is therefore enriched for mito transcripts. A high mito fraction is a hallmark of a broken/apoptotic cell. Typical cutoff 5-20% but tissue-dependent: cardiomyocytes, hepatocytes, and brown fat are mito-rich at baseline, so a 10% blanket cutoff would discard healthy cells.

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

An empty droplet captures only ambient RNA, so it has few distinct genes and low depth. Cells below ~200 genes are usually not real cells.

• Very high `n_genes_by_counts` / `total_counts` -> doublet. Two cells in

one droplet contribute two transcriptomes, roughly doubling both depth and gene diversity. Extreme upper-tail cells are doublet-suspicious — but a high count alone is weak evidence (a large, transcriptionally active cell also has high counts). Use a dedicated doublet caller (Section 4) rather than a hard upper gene cap.

• High `pct_counts_hb` -> red-blood-cell / blood contamination in a solid

tissue dissociation. Often filtered in non-blood tissues.

• `pct_counts_ribo` is informative but rarely a hard filter: very high ribo

fraction can indicate low-complexity / stressed cells, and ribo content is strongly cell-type-specific (proliferating cells are ribo-high), so flag and investigate rather than blindly filter.

---

3. Choosing thresholds — distribution-aware, not magic numbers

Hardcoded cutoffs (mt < 5%, n_genes < 2500) are a starting point, not an answer. They fail on mito-rich tissues and on shallow vs deep libraries. Prefer a data-driven, MAD-based outlier rule that adapts to the dataset.

MAD-based outlier detection (recommended)

Flag a cell as an outlier on a metric when it lies more than nmads median absolute deviations from the median. MAD is robust to the very outliers we are trying to find (unlike mean/SD).

import numpy as np

def is_outlier(adata, metric, nmads=5, upper_only=False):
    M = adata.obs[metric].astype(float)
    med = np.median(M)
    mad = np.median(np.abs(M - med))
    if mad == 0:
        return np.zeros(len(M), dtype=bool)
    lower = med - nmads * mad
    upper = med + nmads * mad
    if upper_only:
        return M > upper
    return (M < lower) | (M > upper)

# Apply on log1p-scaled count metrics (counts are right-skewed)
import scanpy as sc
sc.pp.calculate_qc_metrics(adata, qc_vars=['mt'], percent_top=None,
                           log1p=True, inplace=True)
adata.obs['outlier'] = (
    is_outlier(adata, 'log1p_total_counts', 5)
    | is_outlier(adata, 'log1p_n_genes_by_counts', 5)
)
# Mito is asymmetric: only the HIGH side is bad. Use a tighter nmads (3) AND a
# biological ceiling so a mito-rich tissue doesn't pass everything.
adata.obs['mt_outlier'] = (
    is_outlier(adata, 'pct_counts_mt', 3, upper_only=True)
    | (adata.obs['pct_counts_mt'] > 20)
)
keep = ~(adata.obs['outlier'] | adata.obs['mt_outlier'])
print(f"Keeping {keep.sum()}/{adata.n_obs} cells")
adata = adata[keep].copy()

Guidance:

• nmads=5 for total_counts / n_genes_by_counts (catches extreme tails on

both ends — empty droplets and gross doublets).

• nmads=3, upper-only for pct_counts_mt (dying cells are one-sided).
• Always pair the MAD rule with a biological sanity ceiling on mito

(e.g. 20% general, higher for mito-rich tissue) so a uniformly degraded sample doesn't "pass" just because everything is equally bad.

• Visualize before committing: violin/scatter of the three metrics, and a

total_counts vs pct_counts_mt scatter (the classic L-shape — dying cells sit in the low-count / high-mito corner).

Always inspect distributions first

sc.pl.violin(adata, ['n_genes_by_counts', 'total_counts', 'pct_counts_mt'],
             jitter=0.4, multi_panel=True)
sc.pl.scatter(adata, x='total_counts', y='pct_counts_mt')
sc.pl.scatter(adata, x='total_counts', y='n_genes_by_counts')

---

4. Doublet detection (Scrublet / scDblFinder)

Per-cell counts cannot reliably separate a doublet from a large active cell. Use a simulation-based caller:

# Scrublet via scanpy — run on RAW counts, per-sample (not on merged batches)
sc.pp.scrublet(adata, expected_doublet_rate=0.06)   # scanpy >=1.10
# older scanpy: sc.external.pp.scrublet(adata, expected_doublet_rate=0.06)
n_doub = int(adata.obs['predicted_doublet'].sum())
print(f"{n_doub} predicted doublets ({100*n_doub/adata.n_obs:.1f}%)")
# Often better to FLAG, cluster, then drop — doublets form bridge clusters
adata = adata[~adata.obs['predicted_doublet']].copy()

• expected_doublet_rate scales with loading: ~0.8% per 1,000 cells recovered

(10x). 5,000 cells -> ~4%; 10,000 -> ~8%.

• Run per sample/lane, before merging — cross-sample "doublets" aren't real.
• scDblFinder (R/Bioconductor) is the alternative and often more accurate;

run via Rscript if R is available. Same principle: simulate doublets, score each cell, threshold.

• Prefer flag-cluster-drop over hard pre-clustering removal: real doublets

collapse into recognizable bridge clusters between two parent types.

---

5. Ambient RNA awareness (SoupX / DecontX)

Ambient ("soup") RNA is cell-free mRNA released by lysed cells that gets co-encapsulated into every droplet. Effect: highly expressed genes from one population bleed into all clusters, smearing marker specificity.

• Per-cell QC metrics do NOT detect ambient contamination — it is a

count-correction step, not a cell-filtering step.

• SoupX (R) estimates the soup profile from empty droplets and subtracts it

from cell counts. Needs the raw (unfiltered) + filtered matrices.

• DecontX (celda, R) is a Python-callable-via-R alternative.
• Awareness rule for this skill: if downstream markers look implausibly

ubiquitous (e.g. hemoglobin in every cluster of a non-blood tissue, or a dominant cell type's markers everywhere), suspect ambient RNA and recommend SoupX/DecontX correction. Do not silently proceed.

---

6. Empty-droplet filtering (knee point / EmptyDrops)

Distinguishing real cells from empty droplets is upstream of per-cell QC.

• Knee/inflection on the barcode-rank plot: rank barcodes by total UMI,

plot rank vs total counts (log-log). The "knee" separates real cells (high counts) from the ambient plateau. CellRanger's filtered matrix already applies an EmptyDrops-style call.

• EmptyDrops (DropletUtils, R) tests each low-count barcode against the

ambient profile — recovers real small cells that a hard knee cutoff drops.

• If you only have the filtered matrix, empty-droplet removal is largely

already done; still apply min_genes (~200) and min_counts as a backstop.

• If you have the raw matrix, run EmptyDrops (or the knee heuristic) before

per-cell QC.

# Knee heuristic (when only raw counts available, no DropletUtils)
import numpy as np
tot = np.asarray(adata.X.sum(1)).ravel()
order = np.argsort(tot)[::-1]
ranked = tot[order]
# crude knee: largest drop in log-counts among the top barcodes
log_counts = np.log10(ranked + 1)
knee = np.argmax(np.diff(log_counts[:5000]) * -1) if len(ranked) > 1 else 0
print(f"Approx knee at rank {knee}, threshold ~{ranked[knee]:.0f} UMIs")

---

7. Recommended order of operations

1. (raw matrix only) Empty-droplet call — EmptyDrops or knee. Skip if you have the CellRanger filtered matrix. 2. calculate_qc_metrics with mt/ribo/hb flags (percent_top=None). 3. Inspect distributions (violin + count-vs-mito scatter). 4. Gate cells: MAD-based outliers on counts/genes + mito (3 MAD upper, capped), plus min_genes ~200, min_cells ~3 for genes. 5. Doublet detection (Scrublet/scDblFinder) per sample — flag, optionally cluster, then drop. 6. (optional) Ambient correction (SoupX/DecontX) if markers look smeared. 7. Proceed to normalization (references/scanpy_workflow.md Phase 3+).

Thresholds are dataset-specific. Report the cutoffs used and how many cells each step removed. Never report a filtered cell count without the gates applied.

---

See also

• references/scanpy_workflow.md — Phase 2 inline QC + full pipeline
• scripts/scrna_qc.py — run-if-available helper (computes metrics + MAD

gating from an .h5ad; prints an install plan if scanpy is absent)

• scripts/qc_metrics.py — existing filter/scrublet helpers

Seurat to Scanpy Translation

Quick reference for Seurat (R) users transitioning to Scanpy (Python).

---

Seurat vs Scanpy Equivalents

---

Side-by-Side Workflows

Seurat (R)

library(Seurat)

# Load and create object
data <- Read10X("filtered_gene_bc_matrices/hg19/")
seurat <- CreateSeuratObject(counts = data, min.cells = 3, min.features = 200)

# QC
seurat[["percent.mt"]] <- PercentageFeatureSet(seurat, pattern = "^MT-")
seurat <- subset(seurat, subset = percent.mt < 20)

# Normalize
seurat <- NormalizeData(seurat)
seurat <- FindVariableFeatures(seurat, nfeatures = 2000)
seurat <- ScaleData(seurat)

# PCA
seurat <- RunPCA(seurat, npcs = 50)

# Cluster
seurat <- FindNeighbors(seurat, dims = 1:30)
seurat <- FindClusters(seurat, resolution = 0.5)
seurat <- RunUMAP(seurat, dims = 1:30)

# Find markers
markers <- FindMarkers(seurat, ident.1 = 0, ident.2 = 1)

Scanpy (Python)

import scanpy as sc

# Load
adata = sc.read_10x_mtx("filtered_gene_bc_matrices/hg19/")

# 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
sc.pp.normalize_total(adata, target_sum=1e4)
sc.pp.log1p(adata)

# HVG + Scale
sc.pp.highly_variable_genes(adata, n_top_genes=2000)
sc.pp.scale(adata, max_value=10)

# PCA
sc.tl.pca(adata, n_comps=50)

# Cluster
sc.pp.neighbors(adata, n_pcs=30)
sc.tl.leiden(adata, resolution=0.5)
sc.tl.umap(adata)

# Find markers
sc.tl.rank_genes_groups(adata, groupby='leiden', method='wilcoxon')
markers_0 = sc.get.rank_genes_groups_df(adata, group='0')

---

Key Differences

1. Data Structure

• Seurat: Seurat object with @data, @meta.data, @reductions
• Scanpy: AnnData object with .X, .obs, .var, .obsm, .varm

2. Normalization

• Seurat: NormalizeData() does log-normalization in one step
• Scanpy: Two steps: normalize_total() then log1p()

3. Clustering

• Seurat: Louvain algorithm via FindClusters()
• Scanpy: Leiden (recommended) or Louvain via sc.tl.leiden() or sc.tl.louvain()

4. Subsetting

• Seurat: subset(seurat, subset = cell_type == "T cells")
• Scanpy: adata[adata.obs['cell_type'] == 'T cells']

5. Metadata

• Seurat: seurat@meta.data$new_column <- values
• Scanpy: adata.obs['new_column'] = values

---

Converting Between Formats

Seurat to AnnData

import scanpy as sc

# Save from R
# saveRDS(seurat, "seurat.rds")

# Load in Python (requires rpy2)
import anndata2ri
anndata2ri.activate()

from rpy2.robjects import r
r('library(Seurat)')
r('seurat <- readRDS("seurat.rds")')
r('SaveH5Seurat(seurat, "seurat.h5seurat")')
r('Convert("seurat.h5seurat", "h5ad")')

adata = sc.read_h5ad("seurat.h5ad")

AnnData to Seurat

# Save from Python
adata.write_h5ad("adata.h5ad")

# Load in R
library(SeuratDisk)
Convert("adata.h5ad", "seurat.h5seurat")
seurat <- LoadH5Seurat("seurat.h5seurat")

---

When to Use Each

Use Seurat (R) when:

• Working in R ecosystem (Bioconductor, ggplot2)
• Using Seurat-specific methods (SCTransform, CCA integration)
• Team uses R

Use Scanpy (Python) when:

• Working in Python ecosystem (pandas, scikit-learn, PyTorch)
• Need integration with ML pipelines
• Large datasets (Scanpy generally faster)
• Prefer Leiden clustering

---

See Also

• scanpy_workflow.md - Full Scanpy pipeline
• clustering_guide.md - Leiden vs Louvain

Trajectory Inference and Pseudotime Analysis

Guide to trajectory analysis and pseudotime ordering in single-cell data.

---

Overview

Trajectory analysis orders cells along a developmental or temporal continuum (pseudotime). Use for:

• Differentiation studies (stem cells → mature cells)
• Cell cycle analysis
• Response to stimulation over time
• Disease progression

---

Diffusion Pseudotime (DPT)

Built into scanpy, based on diffusion maps.

import scanpy as sc

# After standard preprocessing and clustering
# 1. Identify root cell (start of trajectory)
adata.uns['iroot'] = np.flatnonzero(
    adata.obs['leiden'] == '0'  # Choose starting cluster
)[0]

# 2. Compute diffusion map
sc.tl.diffmap(adata)

# 3. Compute DPT
sc.tl.dpt(adata)

# 4. Visualize
sc.pl.umap(adata, color=['dpt_pseudotime', 'leiden'])

# Pseudotime values in: adata.obs['dpt_pseudotime']

---

PAGA (Partition-based Graph Abstraction)

Models trajectories as cluster connectivity graph.

# After clustering
sc.tl.paga(adata, groups='leiden')

# Plot PAGA graph
sc.pl.paga(adata, color='leiden')

# Initialize UMAP positions from PAGA
sc.tl.umap(adata, init_pos='paga')

# Plot trajectory
sc.pl.umap(adata, color=['leiden', 'CD34'])  # Stem cell marker

---

Gene Expression Along Pseudotime

# Genes that change along pseudotime
sc.tl.rank_genes_groups(adata, groupby='leiden', method='wilcoxon')

# Plot gene expression vs pseudotime
import matplotlib.pyplot as plt

genes_of_interest = ['CD34', 'CD38', 'CD14']
for gene in genes_of_interest:
    if gene in adata.var_names:
        plt.figure()
        plt.scatter(
            adata.obs['dpt_pseudotime'],
            adata[:, gene].X.toarray().flatten(),
            alpha=0.3, s=5
        )
        plt.xlabel('Pseudotime')
        plt.ylabel(f'{gene} expression')
        plt.title(f'{gene} along trajectory')
        plt.show()

---

External Tools

PAGA with RNA velocity (scVelo)

import scvelo as scv

# Compute RNA velocity
scv.pp.filter_and_normalize(adata)
scv.pp.moments(adata)
scv.tl.velocity(adata)
scv.tl.velocity_graph(adata)

# Visualize
scv.pl.velocity_embedding_stream(adata, basis='umap')

Monocle-style (via Python wrapper)

# Not native to scanpy, requires specialized packages
# For Monocle analysis, consider using R/Seurat integration

---

Tips

1. Root cell selection: Critical for pseudotime. Choose starting cell type. 2. Linear vs branched: DPT handles linear trajectories, PAGA handles branching. 3. Validation: Check marker gene expression matches biological expectation. 4. Multiple trajectories: For complex differentiation, use PAGA.

---

See Also

• scanpy_workflow.md - Preprocessing before trajectory
• clustering_guide.md - Clustering for PAGA

Troubleshooting Guide

Common errors and solutions for single-cell analysis.

---

Installation Issues

ModuleNotFoundError: scanpy

pip install scanpy anndata

ModuleNotFoundError: leidenalg

Leiden clustering requires leidenalg package:

pip install leidenalg

ModuleNotFoundError: umap

UMAP requires umap-learn:

pip install umap-learn

ModuleNotFoundError: harmonypy

Batch correction requires Harmony:

pip install harmonypy

---

Data Loading Issues

Wrong Matrix Orientation

Problem: Genes and samples/cells swapped

Solution: Check and transpose

# AnnData expects: cells as rows, genes as columns
if df.shape[0] > df.shape[1] * 5:
    print("Transposing: genes were rows")
    df = df.T
adata = ad.AnnData(df)

Index Mismatch

Problem: Metadata doesn't align with expression

Solution: Find common indices

common = adata.obs_names.intersection(meta.index)
adata = adata[common].copy()
for col in meta.columns:
    adata.obs[col] = meta.loc[common, col]

Gene Name Mismatch

Problem: Gene names in different formats (Ensembl vs Symbol)

Solution: Convert IDs

from tooluniverse import ToolUniverse
tu = ToolUniverse()
tu.load_tools()

result = tu.tools.MyGene_batch_query(
    gene_ids=['ENSG00000141510', 'ENSG00000139618'],
    fields='symbol,ensembl.gene'
)
# Extract symbols from result

---

Sparse Matrix Issues

TypeError: can't multiply sequence

Problem: Operating on sparse matrix incorrectly

Solution: Convert to dense

from scipy.sparse import issparse

X = adata.X.toarray() if issparse(adata.X) else adata.X

Memory Error

Problem: Dataset too large to convert to dense

Solution: Use sparse operations

# Good: Works on sparse
mean_expr = np.array(adata.X.mean(axis=0)).flatten()

# Bad: Converts entire matrix to dense
# mean_expr = adata.X.toarray().mean(axis=0)

---

QC and Filtering Issues

No mitochondrial genes found

Problem: Gene names don't start with "MT-"

Solution: Check prefix

# Check gene name format
print(adata.var_names[:10])

# Try different prefixes
adata.var['mt'] = adata.var_names.str.startswith(('MT-', 'mt-', 'Mt-'))

# Or manually specify
mt_genes = ['MT-CO1', 'MT-CO2', 'MT-ND1', ...]  # Add all MT genes
adata.var['mt'] = adata.var_names.isin(mt_genes)

All cells filtered out

Problem: QC thresholds too stringent

Solution: Adjust thresholds

# Check distributions first
print(adata.obs['n_genes_by_counts'].describe())
print(adata.obs['pct_counts_mt'].describe())

# Adjust thresholds
sc.pp.filter_cells(adata, min_genes=100)  # Lower from 200
adata = adata[adata.obs['pct_counts_mt'] < 25].copy()  # Higher from 20

---

Clustering Issues

ValueError: n_neighbors must be less than n_samples

Problem: Too few cells for neighbor graph

Solution: Reduce n_neighbors

# Default is 15, reduce for small datasets
sc.pp.neighbors(adata, n_neighbors=min(15, adata.n_obs - 1))

Only one cluster

Problem: Resolution too low

Solution: Increase resolution

# Default is 1.0, increase for more clusters
sc.tl.leiden(adata, resolution=1.5)  # or 2.0

Too many clusters

Problem: Resolution too high

Solution: Decrease resolution

sc.tl.leiden(adata, resolution=0.3)  # or 0.5

---

Differential Expression Issues

ValueError: Not enough cells

Problem: < 3 cells per condition

Solution: Skip cell types with insufficient cells

n_treat = (adata_ct.obs['condition'] == 'treatment').sum()
n_ctrl = (adata_ct.obs['condition'] == 'control').sum()

if n_treat < 3 or n_ctrl < 3:
    print(f"Skipping {cell_type}: insufficient cells")
    continue

NaN values in results

Problem: Gene not expressed or zero variance

Solution: Filter before analysis

# Filter lowly expressed genes
sc.pp.filter_genes(adata, min_cells=3)

# Remove genes with zero variance
from sklearn.feature_selection import VarianceThreshold
selector = VarianceThreshold(threshold=0)
X_filtered = selector.fit_transform(X)

---

Statistical Issues

NaN in correlation

Problem: NaN values in data

Solution: Filter NaN before computing

valid = ~np.isnan(gene_lengths) & ~np.isnan(mean_expr)
r, p = stats.pearsonr(gene_lengths[valid], mean_expr[valid])

Inf values in log transform

Problem: log(0) = -inf

Solution: Use log1p or pseudocount

# Good: log1p(x) = log(x + 1)
sc.pp.log1p(adata)

# Good: log(x + pseudocount)
X_log = np.log10(X + 1)

# Bad: log(x) directly
# X_log = np.log10(X)  # Will produce -inf for zeros

---

Performance Issues

Memory error on large datasets

Problem: Dataset too large to fit in memory

Solution: Subsample or use HVG

# Option 1: Subsample cells
sc.pp.subsample(adata, fraction=0.5)

# Option 2: Use only highly variable genes
sc.pp.highly_variable_genes(adata, n_top_genes=2000)
adata = adata[:, adata.var['highly_variable']].copy()

Slow clustering

Problem: Large dataset (>100k cells)

Solution: Reduce PCs or subsample

# Use fewer PCs
sc.pp.neighbors(adata, n_pcs=20)  # instead of 30

# Or subsample for exploration
adata_sample = sc.pp.subsample(adata, fraction=0.1, copy=True)
# Analyze adata_sample first

---

ToolUniverse Integration Issues

API timeout

Problem: ToolUniverse API call times out

Solution: Reduce query size or retry

# Split large queries
batch_size = 50
for i in range(0, len(gene_list), batch_size):
    batch = gene_list[i:i+batch_size]
    result = tu.tools.MyGene_batch_query(gene_ids=batch)

ID conversion fails

Problem: Gene IDs not recognized

Solution: Try different ID types

# Try Ensembl ID
result = tu.tools.MyGene_query_genes(query="ENSG00000141510")

# Try gene symbol
result = tu.tools.MyGene_query_genes(query="TP53")

# Try with species
result = tu.tools.ensembl_lookup_gene(
    gene_id="ENSG00000141510",
    species="homo_sapiens"
)

---

OmniPath / Cell Communication Issues

No interactions found

Problem: Gene names don't match database

Solution: Check gene name format

# OmniPath uses gene symbols (not Ensembl IDs)
# Convert first if needed
result = tu.tools.STRING_map_identifiers(
    protein_ids=['ENSG00000141510'],
    species=9606
)
# Use preferredName from result

Empty communication matrix

Problem: Expression thresholds too stringent

Solution: Reduce thresholds

expressed_lr = filter_expressed_lr_pairs(
    adata, lr_pairs,
    min_frac=0.02,  # Lower from 0.05
    min_mean=0.01   # Lower from 0.05
)

---

Batch Correction Issues

Harmony not converging

Problem: Batch effects too strong

Solution: Adjust Harmony parameters

# Increase max iterations
ho = harmonypy.run_harmony(
    adata.obsm['X_pca'],
    adata.obs,
    'batch',
    max_iter_harmony=20  # Increase from default 10
)

Overcorrection

Problem: Biological variation removed

Solution: Use fewer PCs or weaker correction

# Use fewer PCs
ho = harmonypy.run_harmony(
    adata.obsm['X_pca'][:, :20],  # Instead of :30
    adata.obs,
    'batch'
)

---

File Format Issues

h5ad version mismatch

Problem: "Cannot read h5ad file"

Solution: Update anndata

pip install --upgrade anndata

10X format changed

Problem: Features file instead of genes file

Solution: Specify file names

# CellRanger v3+
adata = sc.read_10x_mtx(
    "filtered_feature_bc_matrix/",
    var_names='gene_symbols',
    cache=True
)

---

Tips

1. Always check data shape: Cells vs genes orientation 2. Check for NaN/Inf: Before statistical tests 3. Visualize before filtering: Check QC metric distributions 4. Save intermediate results: After time-consuming steps 5. Use try/except: For robust per-cell-type analysis

---

Debugging Checklist

• [ ] Data loaded correctly (check shape, obs, var)
• [ ] Matrix oriented correctly (cells as rows)
• [ ] Metadata aligned with expression data
• [ ] No NaN/Inf values in expression
• [ ] QC thresholds appropriate for dataset
• [ ] Sufficient cells per condition (>= 3)
• [ ] Gene names match between data and annotations
• [ ] Sparse matrix handled correctly
• [ ] Random seed set for reproducibility

---

See Also

• scanpy_workflow.md - Standard pipeline
• clustering_guide.md - Clustering troubleshooting
• cell_communication.md - OmniPath API issues