Skills › Scientific Agent Skills › scanpy › find_markers.py
find_markers.py
90 lines · skills/scanpy/scripts/find_markers.py at 5f2bc03
1#!/usr/bin/env python3
2"""
3Rank marker genes per group and export tables + plots.
5Runs ``sc.tl.rank_genes_groups`` for a grouping (e.g. leiden clusters), writes
6a combined CSV of the top markers per group, per-group CSVs, and the standard
7marker plots (rank panel, heatmap, dotplot).
9NOTE: per-cell tests inflate significance because cells are not independent.
10Use this for EXPLORATORY cluster markers. For rigorous DE between conditions,
11use pseudobulk.py + pydeseq2.
13Examples:
14 python find_markers.py clustered.h5ad -o markers.h5ad
15 python find_markers.py clustered.h5ad --groupby leiden --method wilcoxon --n-genes 50
16 python find_markers.py clustered.h5ad --csv-dir results/markers --top 10
17"""
19import argparse
20import os
22from _common import add_io_args, configure_scanpy, die, info, load_anndata, save_anndata, ensure_categories
25def safe_name(value):
26 from urllib.parse import quote
27 return quote(str(value), safe="")
30def main():
31 p = argparse.ArgumentParser(description=__doc__,
32 formatter_class=argparse.RawDescriptionHelpFormatter)
33 add_io_args(p, default_output=None)
34 p.add_argument("--groupby", default="leiden", help="obs column to group by (default leiden)")
35 p.add_argument("--method", default="wilcoxon",
36 choices=["wilcoxon", "t-test", "t-test_overestim_var", "logreg"],
37 help="Test method (default wilcoxon)")
38 p.add_argument("--n-genes", type=int, default=25, help="Genes shown in plots (default 25)")
39 p.add_argument("--top", type=int, default=25, help="Top markers per group in CSV (default 25)")
40 p.add_argument("--use-raw", action=argparse.BooleanOptionalAction, default=None,
41 help="Rank on adata.raw (normalized log values) instead of X")
42 p.add_argument("--layer", default=None, help="Explicit log-normalized layer (disables raw)")
43 p.add_argument("--groups", nargs="+", default=None, help="Groups to compare (default all)")
44 p.add_argument("--reference", default="rest", help="Reference label or rest")
45 p.add_argument("--csv-dir", default="results/markers", help="Directory for marker CSVs")
46 p.add_argument("--no-plots", action="store_true", help="Skip marker plots")
47 args = p.parse_args()
49 sc = configure_scanpy(figdir=args.figdir)
50 adata = load_anndata(args.input)
51 if args.groupby not in adata.obs.columns:
52 die(f"groupby column '{args.groupby}' not found in obs: {list(adata.obs.columns)}")
54 ensure_categories(adata, args.groupby)
55 if args.use_raw is True and adata.raw is None:
56 die("--use-raw requested but .raw is missing")
57 if args.layer and args.use_raw is True:
58 die("choose --layer or --use-raw, not both")
59 use_raw = False if args.layer else (adata.raw is not None if args.use_raw is None else args.use_raw)
60 sc.tl.rank_genes_groups(adata, args.groupby, method=args.method,
61 use_raw=use_raw, layer=args.layer, groups=args.groups or "all",
62 reference=args.reference, corr_method="benjamini-hochberg", tie_correct=True)
64 os.makedirs(args.csv_dir, exist_ok=True)
65 # logreg, selected groups and explicit references can return fewer groups.
66 groups = list(adata.uns["rank_genes_groups"]["names"].dtype.names)
67 combined = []
68 for g in groups:
69 df = sc.get.rank_genes_groups_df(adata, group=str(g)).head(args.top)
70 df.insert(0, "group", g)
71 df.to_csv(os.path.join(args.csv_dir, f"markers_{safe_name(args.groupby)}_{safe_name(g)}.csv"), index=False)
72 combined.append(df)
73 if combined:
74 import pandas as pd
75 all_path = os.path.join(args.csv_dir, f"markers_{safe_name(args.groupby)}_all.csv")
76 pd.concat(combined, ignore_index=True).to_csv(all_path, index=False)
77 info(f"Wrote marker tables to {args.csv_dir}/ ({len(groups)} groups)")
79 if not args.no_plots:
80 sc.pl.rank_genes_groups(adata, n_genes=args.n_genes, sharey=False,
81 show=False, save="_markers.png")
82 sc.pl.rank_genes_groups_dotplot(adata, n_genes=5, use_raw=use_raw, layer=args.layer, show=False, save="_markers_dotplot.png")
83 sc.pl.rank_genes_groups_heatmap(adata, n_genes=10, use_raw=use_raw, layer=args.layer, show=False, save="_markers_heatmap.png")
85 if args.output:
86 save_anndata(adata, args.output)
89if __name__ == "__main__":
90 main()