Skills › Scientific Agent Skills › scanpy › cluster.py
cluster.py
63 lines · skills/scanpy/scripts/cluster.py at 5f2bc03
1#!/usr/bin/env python3
2"""
3Leiden clustering on a precomputed neighborhood graph.
5Runs Leiden at one or more resolutions and writes a UMAP colored by each
6clustering. Requires that ``sc.pp.neighbors`` has already been run
7(use reduce_dimensions.py first).
9Examples:
10 python cluster.py reduced.h5ad -o clustered.h5ad --resolution 0.5
11 python cluster.py reduced.h5ad -o clustered.h5ad --resolution 0.3 0.5 0.8 1.0
12 python cluster.py reduced.h5ad -o clustered.h5ad --algorithm louvain
13"""
15import argparse
17from _common import add_io_args, configure_scanpy, die, info, load_anndata, save_anndata
20def main():
21 p = argparse.ArgumentParser(description=__doc__,
22 formatter_class=argparse.RawDescriptionHelpFormatter)
23 add_io_args(p, default_output="clustered.h5ad")
24 p.add_argument("--resolution", type=float, nargs="+", default=[0.5],
25 help="One or more resolutions (default 0.5). Higher = more clusters")
26 p.add_argument("--algorithm", default="leiden", choices=["leiden", "louvain"],
27 help="Clustering algorithm (default leiden)")
28 p.add_argument("--key", default=None,
29 help="obs key for the result (single resolution only; "
30 "default '<algorithm>'). Multiple resolutions use '<algorithm>_<res>'")
31 p.add_argument("--no-plots", action="store_true", help="Skip UMAP plots")
32 args = p.parse_args()
34 sc = configure_scanpy(figdir=args.figdir)
35 adata = load_anndata(args.input)
36 if "neighbors" not in adata.uns:
37 die("no neighborhood graph found. Run reduce_dimensions.py first.")
39 cluster_fn = sc.tl.leiden if args.algorithm == "leiden" else sc.tl.louvain
40 keys = []
41 for res in args.resolution:
42 if len(args.resolution) == 1:
43 key = args.key or args.algorithm
44 else:
45 key = f"{args.algorithm}_{res}"
46 # flavor='igraph' is the scanpy 1.12 default-recommended Leiden backend.
47 kwargs = {"resolution": res, "key_added": key}
48 if args.algorithm == "leiden":
49 kwargs.update(flavor="igraph", n_iterations=2, directed=False)
50 cluster_fn(adata, **kwargs)
51 n = adata.obs[key].nunique()
52 info(f"{key}: {n} clusters at resolution {res}")
53 keys.append(key)
55 if not args.no_plots:
56 sc.pl.umap(adata, color=keys, legend_loc="on data",
57 show=False, save="_clusters.png")
59 save_anndata(adata, args.output)
62if __name__ == "__main__":
63 main()