|
| 1 | +#!/usr/bin/env python3 |
| 2 | +"""Create a tiny synthetic scRNA-seq AnnData file for smoke tests and tutorials. |
| 3 | +
|
| 4 | +The dataset is intentionally simple: three PBMC-like populations, sparse integer |
| 5 | +counts, standard QC gene prefixes, and basic sample metadata. It is not intended |
| 6 | +for biological benchmarking. |
| 7 | +""" |
| 8 | + |
| 9 | +import argparse |
| 10 | +import os |
| 11 | + |
| 12 | +import anndata as ad |
| 13 | +import numpy as np |
| 14 | +import pandas as pd |
| 15 | +from scipy import sparse |
| 16 | + |
| 17 | + |
| 18 | +def make_gene_names(n_background: int = 560) -> list[str]: |
| 19 | + marker_genes = [ |
| 20 | + "CD3D", "CD3E", "CD3G", "CD4", "IL7R", "CCR7", "SELL", |
| 21 | + "CD8A", "CD8B", "GZMB", "GZMK", "NKG7", "GNLY", "MS4A1", |
| 22 | + "CD79A", "CD79B", "CD14", "LYZ", "S100A8", "S100A9", |
| 23 | + "FCGR3A", "MS4A7", "FCER1A", "CST3", "LILRA4", "IL3RA", |
| 24 | + "PPBP", "PF4", "HBB", "HBA1", "HBA2", |
| 25 | + ] |
| 26 | + mt_genes = [f"MT-{name}" for name in ["ND1", "ND2", "CO1", "CO2", "ATP6", "CYB"]] |
| 27 | + ribo_genes = [f"RPL{i}" for i in range(3, 13)] + [f"RPS{i}" for i in range(3, 13)] |
| 28 | + background = [f"GENE{i:04d}" for i in range(n_background)] |
| 29 | + return marker_genes + mt_genes + ribo_genes + background |
| 30 | + |
| 31 | + |
| 32 | +def create_toy_adata(n_cells: int, seed: int) -> ad.AnnData: |
| 33 | + rng = np.random.default_rng(seed) |
| 34 | + genes = make_gene_names() |
| 35 | + n_genes = len(genes) |
| 36 | + |
| 37 | + cell_types = np.array(["CD4 T cell", "B cell", "Classical Mono"]) |
| 38 | + probs = np.array([0.45, 0.25, 0.30]) |
| 39 | + assigned = rng.choice(cell_types, size=n_cells, p=probs) |
| 40 | + |
| 41 | + base = rng.negative_binomial(n=2, p=0.92, size=(n_cells, n_genes)).astype(np.int32) |
| 42 | + gene_index = {g: i for i, g in enumerate(genes)} |
| 43 | + |
| 44 | + signatures = { |
| 45 | + "CD4 T cell": ["CD3D", "CD3E", "CD3G", "CD4", "IL7R", "CCR7"], |
| 46 | + "B cell": ["MS4A1", "CD79A", "CD79B"], |
| 47 | + "Classical Mono": ["CD14", "LYZ", "S100A8", "S100A9"], |
| 48 | + } |
| 49 | + for ct, markers in signatures.items(): |
| 50 | + rows = np.where(assigned == ct)[0] |
| 51 | + cols = [gene_index[g] for g in markers] |
| 52 | + base[np.ix_(rows, cols)] += rng.poisson(lam=4, size=(len(rows), len(cols))).astype(np.int32) |
| 53 | + |
| 54 | + # Add a few mitochondrial and ribosomal counts so QC plots are meaningful. |
| 55 | + mt_cols = [i for i, g in enumerate(genes) if g.startswith("MT-")] |
| 56 | + ribo_cols = [i for i, g in enumerate(genes) if g.startswith(("RPL", "RPS"))] |
| 57 | + base[:, mt_cols] += rng.poisson(lam=1, size=(n_cells, len(mt_cols))).astype(np.int32) |
| 58 | + base[:, ribo_cols] += rng.poisson(lam=1, size=(n_cells, len(ribo_cols))).astype(np.int32) |
| 59 | + |
| 60 | + obs = pd.DataFrame( |
| 61 | + { |
| 62 | + "sample": ["toy_pbmc"] * n_cells, |
| 63 | + "batch": ["batch1"] * n_cells, |
| 64 | + "condition": ["control"] * n_cells, |
| 65 | + "cell_type_truth": assigned, |
| 66 | + }, |
| 67 | + index=[f"toy_cell_{i:04d}" for i in range(n_cells)], |
| 68 | + ) |
| 69 | + var = pd.DataFrame(index=genes) |
| 70 | + adata = ad.AnnData(X=sparse.csr_matrix(base), obs=obs, var=var) |
| 71 | + adata.var_names_make_unique() |
| 72 | + return adata |
| 73 | + |
| 74 | + |
| 75 | +def main() -> None: |
| 76 | + parser = argparse.ArgumentParser(description="Create toy scRNA-seq .h5ad data") |
| 77 | + parser.add_argument("--out", default="tests/data/toy_pbmc_raw.h5ad") |
| 78 | + parser.add_argument("--n-cells", type=int, default=300) |
| 79 | + parser.add_argument("--seed", type=int, default=7) |
| 80 | + args = parser.parse_args() |
| 81 | + |
| 82 | + os.makedirs(os.path.dirname(args.out) or ".", exist_ok=True) |
| 83 | + adata = create_toy_adata(n_cells=args.n_cells, seed=args.seed) |
| 84 | + adata.write_h5ad(args.out) |
| 85 | + print(f"Wrote {args.out}: {adata.n_obs} cells x {adata.n_vars} genes") |
| 86 | + |
| 87 | + |
| 88 | +if __name__ == "__main__": |
| 89 | + main() |
0 commit comments