scRNA-seq Quality Control & Preprocessing

When to Use

• Starting from raw count data that needs QC and filtering
• Removing low-quality cells (high mito%, low gene count, doublets)
• Normalizing and log-transforming expression data
• Selecting highly variable genes for downstream analysis
• Preparing data for clustering, trajectory, or communication analysis

Do NOT use when you need to:

• 数据已预处理且有 PCA/embeddings → (use scrna-clustering-umap directly)
• 需要批次校正多个样本 → (use scrna-batch-correction)
• 分析空间转录组数据 → (use spatial-neighborhood-analysis)

Prerequisites

• Python >= 3.10
• Required packages: scanpy, anndata, numpy, scipy, pandas
• Optional: scrublet (doublet detection)
• Input: .h5ad file with raw count matrix (cells x genes). Expected adata.X to contain unnormalized integer counts.

Core Rules

1. Always check if adata.X is sparse before operations; prefer sparse-compatible code (sc.pp.* handles this automatically)
2. Save raw counts in adata.layers["counts"] before normalization — never lose the originals
3. Set batch_key in sc.pp.highly_variable_genes() when multiple samples/batches are present
4. For HVG selection, try seurat_v3 first (requires raw counts), fall back to seurat or cell_ranger with logged data
5. Save intermediate results (e.g., after filtering, after normalization) to avoid re-computation on failures

Common Mistakes

• Re-running on already preprocessed data: Check adata.raw or whether adata.X contains log-normalized values; if so, skip normalization and HVG steps
• Using dense matrix for large datasets: Convert to sparse with adata.X = sparse.csr_matrix(adata.X) to save memory
• Forgetting batch_key: When multiple samples exist, omitting batch_key in HVG selection introduces batch-driven gene selection bias
• Over-filtering: Setting min_genes too high removes rare but biologically important cell types

Troubleshooting

Output

写入 --output 指定的 JSON 文件：

{
  "skill": "scrna-qc-preprocess",
  "status": "success",
  "steps": ["load", "qc_metrics", "doublet_detection", "filtering", "normalization", "hvg_selection", "scaling", "pca"],
  "result": {
    "output_h5ad": "path/to/preprocessed.h5ad",
    "n_cells": 5000,
    "n_genes": 2000,
    "n_hvg": 2000,
    "pct_mito_filtered": 5.2,
    "n_doublets_detected": 120
  }
}

Workflow

Step 1: Load and Inspect Data

import scanpy as sc
import numpy as np
import pandas as pd
from scipy import sparse

adata = sc.read_h5ad("input.h5ad")
print(f"Shape: {adata.shape}")
print(f"X type: {type(adata.X).__name__}")
print(f"Is sparse: {sparse.issparse(adata.X)}")
print(f"obs columns: {list(adata.obs.columns)}")
print(f"var columns: {list(adata.var.columns)}")

Step 2: QC Metrics Calculation

# Identify mitochondrial and ribosomal genes
adata.var["mt"] = adata.var_names.str.startswith(("MT-", "mt-", "Mt-"))
adata.var["ribo"] = adata.var_names.str.startswith(("RPS", "RPL", "rps", "rpl"))

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

# Key QC columns now in adata.obs:
# - n_genes_by_counts: number of genes expressed per cell
# - total_counts: total UMI counts per cell
# - pct_counts_mt: percentage of mitochondrial reads
# - pct_counts_ribo: percentage of ribosomal reads

Step 3: Doublet Detection (Optional but Recommended)

try:
    import scrublet as scr
    scrub = scr.Scrublet(adata.X)
    doublet_scores, predicted_doublets = scrub.scrub_doublets(
        min_counts=2, min_cells=3, n_prin_comps=30
    )
    adata.obs["doublet_score"] = doublet_scores
    adata.obs["predicted_doublet"] = predicted_doublets
    print(f"Detected {predicted_doublets.sum()} doublets ({predicted_doublets.mean():.1%})")
except ImportError:
    print("scrublet not available, skipping doublet detection")

Step 4: Filtering

# Standard filtering thresholds (adjust per dataset)
sc.pp.filter_cells(adata, min_genes=200)
sc.pp.filter_genes(adata, min_cells=3)

# Mito filter: keep cells with < 20% mito reads
adata = adata[adata.obs.pct_counts_mt < 20].copy()

# Gene count filter: remove cells with too many genes (potential doublets)
upper_limit = np.percentile(adata.obs.n_genes_by_counts, 98)
adata = adata[adata.obs.n_genes_by_counts < upper_limit].copy()

print(f"After filtering: {adata.shape}")

Step 5: Normalization and Log Transform

# Save raw counts
adata.layers["counts"] = adata.X.copy()

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

# Log transform
sc.pp.log1p(adata)

# Store normalized data as raw for DEG testing
adata.raw = adata

Step 6: Highly Variable Gene Selection

# Find highly variable genes
# Try multiple flavors for robustness across data types
n_hvg = 0
for flavor in ["seurat_v3", "cell_ranger", "seurat"]:
    try:
        sc.pp.highly_variable_genes(
            adata,
            n_top_genes=2000,
            flavor=flavor,
            layer="counts" if flavor == "seurat_v3" else None,
            batch_key="batch" if "batch" in adata.obs.columns else None,
        )
        n_hvg = adata.var.highly_variable.sum()
        if n_hvg > 0:
            break
    except (ValueError, TypeError, ImportError):
        continue

# Ultimate fallback: variance-based selection
if n_hvg == 0:
    from scipy import sparse as sp
    if sp.issparse(adata.X):
        variances = np.array(adata.X.toarray()).var(axis=0)
    else:
        variances = adata.X.var(axis=0)
    top_idx = np.argsort(variances)[-2000:]
    adata.var["highly_variable"] = False
    adata.var.iloc[top_idx, adata.var.columns.get_loc("highly_variable")] = True
    print("Used variance-based HVG selection (fallback)")

print(f"HVGs: {adata.var.highly_variable.sum()}")

# Subset to HVGs
adata = adata[:, adata.var.highly_variable].copy()

Step 7: Scaling and PCA

# Handle NaN values before scaling
import numpy as np
if sparse.issparse(adata.X):
    adata.X = np.nan_to_num(adata.X.toarray(), nan=0.0)
else:
    adata.X = np.nan_to_num(adata.X, nan=0.0)

# Scale to unit variance
sc.pp.scale(adata, max_value=10)

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

# Store PCA result
print(f"PCA shape: {adata.obsm["X_pca"].shape}")

Step 8: Save Output

adata.write_h5ad("preprocessed.h5ad")
print(f"Saved: {adata.shape}")

Key Parameters to Adjust