Skills › Scientific Agent Skills › scanpy › reduce_dimensions.py
reduce_dimensions.py
64 lines · skills/scanpy/scripts/reduce_dimensions.py at 5f2bc03
1#!/usr/bin/env python3
2"""
3PCA, neighborhood graph, and UMAP / t-SNE embeddings.
5Takes a normalized object (optionally HVG-subset / scaled) and computes PCA,
6the kNN graph, and a UMAP (and optionally t-SNE) embedding. Writes a PCA
7variance-ratio elbow plot to help choose ``--n-pcs``.
9Examples:
10 python reduce_dimensions.py normalized.h5ad -o reduced.h5ad
11 python reduce_dimensions.py normalized.h5ad -o reduced.h5ad --n-pcs 50 --n-neighbors 15
12 python reduce_dimensions.py normalized.h5ad -o reduced.h5ad --tsne --use-rep X_pca_harmony
13"""
15import argparse
17from _common import add_io_args, configure_scanpy, info, load_anndata, save_anndata, compute_pca, build_neighbors, die
20def main():
21 p = argparse.ArgumentParser(description=__doc__,
22 formatter_class=argparse.RawDescriptionHelpFormatter)
23 add_io_args(p, default_output="reduced.h5ad")
24 p.add_argument("--n-comps", type=int, default=50, help="PCs to compute (default 50)")
25 p.add_argument("--n-pcs", type=int, default=40, help="PCs used for the kNN graph (default 40)")
26 p.add_argument("--n-neighbors", type=int, default=15, help="Neighbors for the graph (default 15)")
27 p.add_argument("--use-rep", default=None,
28 help="obsm key to build the graph from instead of PCA "
29 "(e.g. X_pca_harmony from batch_correct.py)")
30 p.add_argument("--tsne", action="store_true", help="Also compute t-SNE")
31 p.add_argument("--color", nargs="+", default=None,
32 help="obs/var keys to color the UMAP by (e.g. sample n_genes_by_counts)")
33 p.add_argument("--no-plots", action="store_true", help="Skip plots")
34 args = p.parse_args()
36 sc = configure_scanpy(figdir=args.figdir)
37 adata = load_anndata(args.input)
38 info(f"Loaded {adata.n_obs} cells x {adata.n_vars} genes")
40 if args.use_rep is None or args.use_rep == "X_pca":
41 n_comps = compute_pca(sc, adata, args.n_comps)
42 if not args.no_plots:
43 sc.pl.pca_variance_ratio(adata, n_pcs=n_comps, log=True,
44 show=False, save="_variance.png")
45 build_neighbors(sc, adata, n_neighbors=args.n_neighbors, n_pcs=args.n_pcs,
46 use_rep=args.use_rep or "X_pca")
47 info("Computing UMAP...")
48 sc.tl.umap(adata)
50 if args.tsne:
51 info("Computing t-SNE...")
52 sc.tl.tsne(adata, use_rep=args.use_rep or "X_pca",
53 perplexity=min(30, max(1, (adata.n_obs - 1) / 3)), random_state=0)
55 if not args.no_plots and args.color:
56 color = [c for c in args.color if c in adata.obs.columns or c in adata.var_names]
57 if color:
58 sc.pl.umap(adata, color=color, show=False, save="_colored.png")
60 save_anndata(adata, args.output)
63if __name__ == "__main__":
64 main()