Skills › Scientific Agent Skills › scanpy › run_pipeline.py
run_pipeline.py
201 lines · skills/scanpy/scripts/run_pipeline.py at 5f2bc03
1#!/usr/bin/env python3
2"""
3End-to-end standard scRNA-seq pipeline in one command.
5Runs the full exploratory workflow on raw counts:
6 load -> QC + filter -> (optional doublets) -> normalize + log1p -> HVG
7 -> (optional scale/regress) -> PCA -> (optional batch correction)
8 -> neighbors -> UMAP -> Leiden -> marker genes -> save.
10Produces a processed .h5ad, marker CSVs, and figures. Tune the common knobs via
11flags, or pass a JSON config with ``--config`` (keys mirror the flag names with
12underscores). This is the fastest path from counts to a clustered, annotated-
13ready object; use the individual step scripts when you need to iterate on one stage.
15Examples:
16 python run_pipeline.py raw.h5ad -o processed.h5ad
17 python run_pipeline.py raw.h5ad -o processed.h5ad --resolution 0.8 --n-top-genes 3000
18 python run_pipeline.py raw.h5ad -o processed.h5ad --batch-key sample --batch-method harmony
19 python run_pipeline.py raw.h5ad -o processed.h5ad --config params.json
20"""
22import argparse
23import json
24import os
26from _common import configure_scanpy, info, load_anndata, save_anndata, prepare_counts, normalize_hvg, compute_pca, build_neighbors, die, ensure_categories, integrate_harmony
29def build_parser():
30 p = argparse.ArgumentParser(description=__doc__,
31 formatter_class=argparse.RawDescriptionHelpFormatter)
32 p.add_argument("input", help="Input raw-counts file (.h5ad, .h5, .csv, 10x dir, ...)")
33 p.add_argument("-o", "--output", default="processed.h5ad", help="Output .h5ad")
34 p.add_argument("--figdir", default="figures", help="Figure directory")
35 p.add_argument("--marker-dir", default="results/markers", help="Marker CSV directory")
36 p.add_argument("--config", default=None, help="JSON file overriding any of the options below")
37 # QC
38 p.add_argument("--min-genes", type=int, default=200)
39 p.add_argument("--max-genes", type=int, default=None)
40 p.add_argument("--min-cells", type=int, default=3)
41 p.add_argument("--mt-threshold", type=float, default=5)
42 p.add_argument("--scrublet", action="store_true")
43 # normalization / HVG
44 p.add_argument("--counts-layer", default=None, help="Explicit raw-count layer instead of X")
45 p.add_argument("--target-sum", type=float, default=1e4)
46 p.add_argument("--n-top-genes", type=int, default=2000)
47 p.add_argument("--hvg-flavor", default="seurat", choices=["seurat", "cell_ranger", "seurat_v3", "seurat_v3_paper"])
48 p.add_argument("--scale", action="store_true")
49 p.add_argument("--regress-out", nargs="+", default=None)
50 # dim reduction / clustering
51 p.add_argument("--n-pcs", type=int, default=40)
52 p.add_argument("--n-neighbors", type=int, default=15)
53 p.add_argument("--resolution", type=float, default=0.5)
54 # batch correction
55 p.add_argument("--batch-key", default=None)
56 p.add_argument("--batch-method", default="harmony", choices=["harmony", "combat"])
57 # markers
58 p.add_argument("--marker-method", default="wilcoxon", choices=["wilcoxon", "t-test", "t-test_overestim_var"])
59 p.add_argument("--skip-markers", action="store_true")
60 return p
63def apply_config(args):
64 if args.config:
65 with open(args.config) as fh:
66 cfg = json.load(fh)
67 if not isinstance(cfg, dict):
68 die("config must be a JSON object")
69 for k, v in cfg.items():
70 key = k.replace("-", "_")
71 if key not in vars(args) or key in {"input", "output", "config"}:
72 die(f"unknown or unsupported config option: {k}")
73 action = next(a for a in build_parser()._actions if a.dest == key)
74 if action.choices and v not in action.choices:
75 die(f"invalid config choice for {k}: {v}")
76 old = getattr(args, key)
77 if old is not None and not isinstance(v, type(old)):
78 if not (type(old) is float and type(v) is int):
79 die(f"invalid config value type for {k}")
80 setattr(args, key, v)
81 return args
84def main():
85 args = apply_config(build_parser().parse_args())
86 sc = configure_scanpy(figdir=args.figdir)
88 info("[1/8] Loading data")
89 adata = load_anndata(args.input)
90 prepare_counts(adata, args.counts_layer)
91 if args.batch_key:
92 ensure_categories(adata, args.batch_key)
93 adata.var_names_make_unique()
94 info(f" {adata.n_obs} cells x {adata.n_vars} genes")
96 info("[2/8] QC + filtering")
97 adata.var["mt"] = adata.var_names.str.startswith(("MT-", "mt-", "Mt-"))
98 qc_vars = ["mt"] if adata.var["mt"].any() else []
99 sc.pp.calculate_qc_metrics(adata, qc_vars=qc_vars, percent_top=None, log1p=False, inplace=True)
100 sc.pl.violin(adata, [k for k in ["n_genes_by_counts", "total_counts", "pct_counts_mt"]
101 if k in adata.obs.columns],
102 jitter=0.4, multi_panel=True, show=False, save="_qc.png")
103 n0 = adata.n_obs
104 sc.pp.filter_cells(adata, min_genes=args.min_genes)
105 sc.pp.filter_genes(adata, min_cells=args.min_cells)
106 if args.max_genes:
107 adata = adata[adata.obs["n_genes_by_counts"] < args.max_genes, :].copy()
108 if qc_vars and "pct_counts_mt" in adata.obs.columns:
109 adata = adata[adata.obs["pct_counts_mt"] < args.mt_threshold, :].copy()
110 if min(adata.shape) == 0:
111 die("QC removed all cells or genes; inspect thresholds")
112 if args.scrublet:
113 try:
114 sc.pp.scrublet(adata, batch_key=args.batch_key, random_state=0)
115 adata = adata[~adata.obs["predicted_doublet"], :].copy()
116 except (ImportError, ValueError) as e:
117 die(f"Requested Scrublet failed: {e}; check counts, capture/library groups and scikit-image")
118 info(f" {n0} -> {adata.n_obs} cells, {adata.n_vars} genes")
120 info("[3/8] Normalize + log1p + HVG")
121 normalize_hvg(sc, adata, target_sum=args.target_sum, n_top_genes=args.n_top_genes,
122 flavor=args.hvg_flavor, batch_key=args.batch_key)
123 info(f" {int(adata.var['highly_variable'].sum())} HVGs")
125 info("[4/8] Scale / regress (optional)")
126 work = adata[:, adata.var["highly_variable"]].copy()
127 if args.regress_out:
128 sc.pp.regress_out(work, args.regress_out)
129 if args.scale:
130 sc.pp.scale(work, max_value=10)
132 info("[5/8] PCA" + (f" + {args.batch_method} batch correction" if args.batch_key else ""))
133 compute_pca(sc, work, args.n_pcs)
134 sc.pl.pca_variance_ratio(work, n_pcs=work.obsm["X_pca"].shape[1], log=True, show=False, save="_variance.png")
135 use_rep = "X_pca"
136 if args.batch_key:
137 if args.batch_method == "harmony":
138 try:
139 integrate_harmony(work, args.batch_key)
140 use_rep = "X_pca_harmony"
141 except ImportError:
142 die("Requested Harmony requires harmonypy; install it before running integration")
143 else:
144 sc.pp.combat(work, key=args.batch_key)
145 compute_pca(sc, work, args.n_pcs)
147 info("[6/8] Neighbors + UMAP")
148 build_neighbors(sc, work, n_neighbors=args.n_neighbors, n_pcs=args.n_pcs, use_rep=use_rep)
149 sc.tl.umap(work)
151 info("[7/8] Leiden clustering")
152 sc.tl.leiden(work, resolution=args.resolution, flavor="igraph",
153 n_iterations=2, directed=False)
154 n_clusters = work.obs["leiden"].nunique()
155 color = ["leiden"] + ([args.batch_key] if args.batch_key else [])
156 sc.pl.umap(work, color=color, legend_loc="on data", show=False, save="_leiden.png")
157 info(f" {n_clusters} clusters at resolution {args.resolution}")
159 # Carry embeddings/clusters back onto the full-gene object so markers use all genes.
160 adata.obs["leiden"] = work.obs["leiden"].copy()
161 adata.obsm["X_pca"] = work.obsm["X_pca"]
162 adata.obsm["X_umap"] = work.obsm["X_umap"]
163 if use_rep in work.obsm:
164 adata.obsm[use_rep] = work.obsm[use_rep]
165 import numpy as np
166 adata.uns["pca"] = work.uns["pca"].copy()
167 adata.uns["pca"]["params"]["mask_var"] = "highly_variable"
168 loadings = np.zeros((adata.n_vars, work.obsm["X_pca"].shape[1]), dtype=work.varm["PCs"].dtype)
169 loadings[adata.var["highly_variable"].to_numpy()] = work.varm["PCs"]
170 adata.varm["PCs"] = loadings
171 adata.uns["umap"] = work.uns["umap"].copy()
172 adata.uns["leiden"] = work.uns["leiden"].copy()
173 adata.uns["neighbors"] = work.uns["neighbors"]
174 adata.obsp = work.obsp
176 if not args.skip_markers and (n_clusters < 2 or work.obs["leiden"].value_counts().min() < 2):
177 die("marker ranking requires at least two groups with at least two cells each; tune clustering or use --skip-markers")
178 if not args.skip_markers:
179 info("[8/8] Marker genes")
180 sc.tl.rank_genes_groups(adata, "leiden", method=args.marker_method, use_raw=True,
181 reference="rest", corr_method="benjamini-hochberg", tie_correct=True)
182 sc.pl.rank_genes_groups_dotplot(adata, n_genes=5, show=False, save="_markers_dotplot.png")
183 os.makedirs(args.marker_dir, exist_ok=True)
184 import pandas as pd
185 frames = []
186 for g in adata.obs["leiden"].cat.categories:
187 df = sc.get.rank_genes_groups_df(adata, group=g).head(25)
188 df.insert(0, "cluster", g)
189 frames.append(df)
190 pd.concat(frames, ignore_index=True).to_csv(
191 os.path.join(args.marker_dir, "markers_all.csv"), index=False)
192 info(f" marker tables in {args.marker_dir}/")
193 else:
194 info("[8/8] Skipping markers")
196 save_anndata(adata, args.output)
197 info("Pipeline complete. Next: inspect markers, then annotate.py with a cluster->cell_type mapping.")
200if __name__ == "__main__":
201 main()