Skills › Scientific Agent Skills › scanpy › qc_analysis.py
qc_analysis.py
112 lines · skills/scanpy/scripts/qc_analysis.py at 5f2bc03
1#!/usr/bin/env python3
2"""
3Quality control and filtering for single-cell RNA-seq data.
5Calculates QC metrics (genes/counts per cell, mitochondrial / ribosomal /
6hemoglobin fractions), writes before/after QC plots, optionally runs Scrublet
7doublet detection, and filters cells and genes by the given thresholds.
9Run this FIRST on raw counts, before normalization.
11Examples:
12 python qc_analysis.py raw.h5ad -o filtered.h5ad
13 python qc_analysis.py raw.h5ad -o filtered.h5ad --mt-threshold 10 --min-genes 500
14 python qc_analysis.py 10x_dir/ -o filtered.h5ad --max-genes 6000 --scrublet
15"""
17import argparse
19from _common import add_io_args, configure_scanpy, info, load_anndata, save_anndata, die, prepare_counts, ensure_categories
22def annotate_gene_classes(adata):
23 """Flag mitochondrial, ribosomal, and hemoglobin genes (human or mouse names)."""
24 names = adata.var_names
25 adata.var["mt"] = names.str.startswith(("MT-", "mt-", "Mt-"))
26 adata.var["ribo"] = names.str.startswith(("RPS", "RPL", "Rps", "Rpl"))
27 adata.var["hb"] = names.str.contains(r"^HB[^P]|^Hb[^p]", regex=True)
28 return [v for v in ["mt", "ribo", "hb"] if adata.var[v].any()]
31def make_qc_plots(sc, adata, prefix):
32 qc_keys = ["n_genes_by_counts", "total_counts", "pct_counts_mt"]
33 qc_keys = [k for k in qc_keys if k in adata.obs.columns]
34 sc.pl.violin(adata, qc_keys, jitter=0.4, multi_panel=True,
35 show=False, save=f"_{prefix}_violin.png")
36 if "pct_counts_mt" in adata.obs.columns:
37 sc.pl.scatter(adata, x="total_counts", y="pct_counts_mt",
38 show=False, save=f"_{prefix}_mt.png")
39 sc.pl.scatter(adata, x="total_counts", y="n_genes_by_counts",
40 show=False, save=f"_{prefix}_counts.png")
43def main():
44 p = argparse.ArgumentParser(description=__doc__,
45 formatter_class=argparse.RawDescriptionHelpFormatter)
46 add_io_args(p, default_output="qc_filtered.h5ad")
47 p.add_argument("--counts-layer", default=None, help="Explicit raw-count layer instead of X")
48 p.add_argument("--batch-key", default=None, help="Capture/library ID for separate Scrublet runs")
49 p.add_argument("--min-genes", type=int, default=200, help="Min genes per cell (default 200)")
50 p.add_argument("--max-genes", type=int, default=None, help="Max genes per cell (upper outliers)")
51 p.add_argument("--min-counts", type=int, default=None, help="Min total counts per cell")
52 p.add_argument("--max-counts", type=int, default=None, help="Max total counts per cell")
53 p.add_argument("--min-cells", type=int, default=3, help="Min cells per gene (default 3)")
54 p.add_argument("--mt-threshold", type=float, default=5, help="Max pct mitochondrial counts (default 5)")
55 p.add_argument("--scrublet", action="store_true", help="Run Scrublet doublet detection and drop doublets")
56 p.add_argument("--no-plots", action="store_true", help="Skip QC plots")
57 args = p.parse_args()
59 sc = configure_scanpy(figdir=args.figdir)
60 adata = load_anndata(args.input)
61 prepare_counts(adata, args.counts_layer)
62 if args.batch_key:
63 ensure_categories(adata, args.batch_key)
64 adata.var_names_make_unique()
65 info(f"Loaded {adata.n_obs} cells x {adata.n_vars} genes")
67 qc_vars = annotate_gene_classes(adata)
68 sc.pp.calculate_qc_metrics(adata, qc_vars=qc_vars, percent_top=None,
69 log1p=False, inplace=True)
70 info(f"Mean genes/cell={adata.obs['n_genes_by_counts'].mean():.0f} "
71 f"mean counts/cell={adata.obs['total_counts'].mean():.0f} "
72 f"mean pct_mt={(adata.obs['pct_counts_mt'].mean() if 'pct_counts_mt' in adata.obs else 0):.1f}")
74 if not args.no_plots:
75 make_qc_plots(sc, adata, "qc_before")
77 n0, g0 = adata.n_obs, adata.n_vars
78 sc.pp.filter_cells(adata, min_genes=args.min_genes)
79 if args.min_counts:
80 sc.pp.filter_cells(adata, min_counts=args.min_counts)
81 if args.max_genes:
82 adata = adata[adata.obs["n_genes_by_counts"] < args.max_genes, :].copy()
83 if args.max_counts:
84 adata = adata[adata.obs["total_counts"] < args.max_counts, :].copy()
85 if "pct_counts_mt" in adata.obs.columns:
86 adata = adata[adata.obs["pct_counts_mt"] < args.mt_threshold, :].copy()
87 sc.pp.filter_genes(adata, min_cells=args.min_cells)
89 if min(adata.shape) == 0:
90 die("QC removed all cells or genes; inspect thresholds")
92 if args.scrublet:
93 info("Running Scrublet doublet detection...")
94 try:
95 sc.pp.scrublet(adata, batch_key=args.batch_key, random_state=0)
96 except (ImportError, ValueError) as e:
97 die(f"Scrublet failed ({e}). Check raw counts, library size and scikit-image installation")
98 n_dbl = int(adata.obs["predicted_doublet"].sum())
99 adata = adata[~adata.obs["predicted_doublet"], :].copy()
100 info(f"Removed {n_dbl} predicted doublets")
102 info(f"Cells {n0} -> {adata.n_obs} ({adata.n_obs / n0 * 100:.1f}% kept) "
103 f"Genes {g0} -> {adata.n_vars}")
105 if not args.no_plots:
106 make_qc_plots(sc, adata, "qc_after")
108 save_anndata(adata, args.output)
111if __name__ == "__main__":
112 main()