Skills › Scientific Agent Skills › scanpy › score_genes.py
score_genes.py
86 lines · skills/scanpy/scripts/score_genes.py at 5f2bc03
1#!/usr/bin/env python3
2"""
3Score cells for one or more gene signatures.
5Computes ``sc.tl.score_genes`` for each named gene set in a JSON file and adds
6one obs column per signature (``<name>_score``). Also supports the built-in
7cell-cycle scoring with ``--cell-cycle``. Writes UMAPs colored by each score.
9Gene-set JSON format:
10 {"T_cell": ["CD3D","CD3E","CD3G"], "cytotoxic": ["GZMB","PRF1","NKG7"]}
12Examples:
13 python score_genes.py annotated.h5ad -o scored.h5ad --gene-sets signatures.json
14 python score_genes.py annotated.h5ad -o scored.h5ad --cell-cycle
15"""
17import argparse
18import json
20from _common import add_io_args, configure_scanpy, info, load_anndata, save_anndata, die
22# Tirosh et al. (2016) legacy human symbols. Map IDs/orthologs explicitly for other organisms.
23S_GENES = ["MCM5", "PCNA", "TYMS", "FEN1", "MCM2", "MCM4", "RRM1", "UNG", "GINS2",
24 "MCM6", "CDCA7", "DTL", "PRIM1", "UHRF1", "MLF1IP", "HELLS", "RFC2",
25 "RPA2", "NASP", "RAD51AP1", "GMNN", "WDR76", "SLBP", "CCNE2", "UBR7",
26 "POLD3", "MSH2", "ATAD2", "RAD51", "RRM2", "CDC45", "CDC6", "EXO1",
27 "TIPIN", "DSCC1", "BLM", "CASP8AP2", "USP1", "CLSPN", "POLA1", "CHAF1B",
28 "BRIP1", "E2F8"]
29G2M_GENES = ["HMGB2", "CDK1", "NUSAP1", "UBE2C", "BIRC5", "TPX2", "TOP2A", "NDC80",
30 "CKS2", "NUF2", "CKS1B", "MKI67", "TMPO", "CENPF", "TACC3", "FAM64A",
31 "SMC4", "CCNB2", "CKAP2L", "CKAP2", "AURKB", "BUB1", "KIF11", "ANP32E",
32 "TUBB4B", "GTSE1", "KIF20B", "HJURP", "CDCA3", "HN1", "CDC20", "TTK",
33 "CDC25C", "KIF2C", "RANGAP1", "NCAPD2", "DLGAP5", "CDCA2", "CDCA8",
34 "ECT2", "KIF23", "HMMR", "AURKA", "PSRC1", "ANLN", "LBR", "CKAP5",
35 "CENPE", "CTCF", "NEK2", "G2E3", "GAS2L3", "CBX5", "CENPA"]
38def main():
39 p = argparse.ArgumentParser(description=__doc__,
40 formatter_class=argparse.RawDescriptionHelpFormatter)
41 add_io_args(p, default_output="scored.h5ad")
42 p.add_argument("--gene-sets", default=None, help="JSON of {name: [genes]}")
43 p.add_argument("--cell-cycle", action="store_true",
44 help="Add S_score, G2M_score, and phase (Tirosh 2016 human genes)")
45 p.add_argument("--no-plots", action="store_true", help="Skip plots")
46 args = p.parse_args()
48 sc = configure_scanpy(figdir=args.figdir)
49 adata = load_anndata(args.input)
50 if not args.gene_sets and not args.cell_cycle:
51 die("provide --gene-sets and/or --cell-cycle")
52 var_names = adata.raw.var_names if adata.raw is not None else adata.var_names
53 plot_keys = []
55 if args.gene_sets:
56 with open(args.gene_sets) as fh:
57 sets = json.load(fh)
58 for name, genes in sets.items():
59 present = [g for g in genes if g in var_names]
60 if not present:
61 info(f"Skipping '{name}': no genes present")
62 continue
63 sc.tl.score_genes(adata, present, score_name=f"{name}_score",
64 use_raw=adata.raw is not None)
65 plot_keys.append(f"{name}_score")
66 info(f"Scored '{name}' ({len(present)}/{len(genes)} genes)")
68 if args.cell_cycle:
69 s = [g for g in S_GENES if g in var_names]
70 g2m = [g for g in G2M_GENES if g in var_names]
71 if not s or not g2m:
72 die("cell-cycle scoring needs both S and G2M genes; verify organism and gene identifiers")
73 sc.tl.score_genes_cell_cycle(adata, s_genes=s, g2m_genes=g2m,
74 use_raw=adata.raw is not None)
75 plot_keys += ["phase"]
76 info(f"Cell-cycle scored (S={len(s)}, G2M={len(g2m)} genes); phases: "
77 + ", ".join(f"{k}={v}" for k, v in adata.obs['phase'].value_counts().items()))
79 if not args.no_plots and plot_keys and "X_umap" in adata.obsm:
80 sc.pl.umap(adata, color=plot_keys, show=False, save="_scores.png")
82 save_anndata(adata, args.output)
85if __name__ == "__main__":
86 main()