Skills › Scientific Agent Skills › scanpy › preprocess.py
preprocess.py
79 lines · skills/scanpy/scripts/preprocess.py at 5f2bc03
1#!/usr/bin/env python3
2"""
3Normalize, log-transform, select highly variable genes, and optionally scale.
5Takes QC-filtered raw counts and produces a normalized, log1p-transformed
6object ready for dimensionality reduction. A copy of the raw counts is kept in
7``adata.layers['counts']`` and the normalized log values in ``adata.raw`` so
8downstream marker/expression plots can use ``use_raw=True``.
10Examples:
11 python preprocess.py filtered.h5ad -o normalized.h5ad
12 python preprocess.py filtered.h5ad -o normalized.h5ad --n-top-genes 3000 --scale
13 python preprocess.py filtered.h5ad -o normalized.h5ad --flavor seurat_v3 --batch-key sample
14 python preprocess.py filtered.h5ad -o normalized.h5ad --regress-out total_counts pct_counts_mt
15"""
17import argparse
19from _common import add_io_args, configure_scanpy, info, load_anndata, save_anndata, prepare_counts, normalize_hvg
22def main():
23 p = argparse.ArgumentParser(description=__doc__,
24 formatter_class=argparse.RawDescriptionHelpFormatter)
25 add_io_args(p, default_output="normalized.h5ad")
26 p.add_argument("--counts-layer", default=None, help="Explicit raw-count layer instead of X")
27 p.add_argument("--target-sum", type=float, default=1e4,
28 help="Counts per cell after normalization (default 1e4)")
29 p.add_argument("--n-top-genes", type=int, default=2000,
30 help="Number of highly variable genes (default 2000)")
31 p.add_argument("--flavor", default="seurat",
32 choices=["seurat", "cell_ranger", "seurat_v3", "seurat_v3_paper"],
33 help="HVG flavor. seurat_v3 expects raw counts (default seurat)")
34 p.add_argument("--batch-key", default=None,
35 help="obs column for batch-aware HVG selection")
36 p.add_argument("--subset-hvg", action="store_true",
37 help="Subset the matrix to HVGs (smaller, but drops other genes from X)")
38 p.add_argument("--regress-out", nargs="+", default=None,
39 help="obs columns to regress out (e.g. total_counts pct_counts_mt)")
40 p.add_argument("--scale", action="store_true",
41 help="Scale to unit variance and zero mean (max_value=10)")
42 p.add_argument("--no-plots", action="store_true", help="Skip HVG plot")
43 args = p.parse_args()
45 sc = configure_scanpy(figdir=args.figdir)
46 adata = load_anndata(args.input)
47 info(f"Loaded {adata.n_obs} cells x {adata.n_vars} genes")
49 prepare_counts(adata, args.counts_layer)
50 normalize_hvg(sc, adata, target_sum=args.target_sum, n_top_genes=args.n_top_genes,
51 flavor=args.flavor, batch_key=args.batch_key)
53 n_hvg = int(adata.var["highly_variable"].sum())
54 info(f"Selected {n_hvg} highly variable genes")
56 # Stash the full normalized log matrix so plots can use_raw=True.
57 # normalize_hvg saved an independent full log-normalized snapshot.
59 if not args.no_plots:
60 sc.pl.highly_variable_genes(adata, show=False, save="_hvg.png")
62 if args.subset_hvg:
63 info("HVG subsetting also subsets counts; retain the full-gene input for pseudobulk")
64 adata = adata[:, adata.var["highly_variable"]].copy()
65 info(f"Subset to {adata.n_vars} HVGs")
67 if args.regress_out:
68 info(f"Regressing out: {', '.join(args.regress_out)}")
69 sc.pp.regress_out(adata, args.regress_out)
71 if args.scale:
72 info("Scaling (max_value=10)")
73 sc.pp.scale(adata, max_value=10)
75 save_anndata(adata, args.output)
78if __name__ == "__main__":
79 main()