Quickstart

Minimal API (recommended default path)

If you want the recommended default path without dozens of parameters, use the simple wrappers:

import scatrans as scat

# One-liner pipeline: score → filter → GO enrichment
result = scat.run_default_pipeline(
    adata,
    groupby="condition",
    target_group="Disease",
    reference_group="Control",
    sample_col="sample",   # optional; auto-selects pseudobulk when >=3 replicates/group
    organism="mouse",
)
print(result["candidates"].head())
print(result["enrichment"].head())

# Or just the core scoring step:
adata_res, significant, all_results = scat.active_score_simple(
    adata,
    groupby="condition",
    target_group="Disease",
    reference_group="Control",
    sample_col="sample",
)

active_score_simple / run_default_pipeline auto-attach gene features, pick Wilcoxon (single-cell) or pseudobulk+PyDESeq2 (when replicates allow), and keep permutation off by default. Use active_score(...) directly for advanced options (permutation, mixed models, Memento, etc.) — see the User Guide.

Complete end-to-end example

This is a complete, copy-paste friendly workflow for first-time users. It takes you from loaded data to differential results, filtering, enrichment, and visualization of enrichment results.

import scanpy as sc
import scatrans as scat

# 1. Load your data (must contain spliced/unspliced layers or use differential_expression instead)
adata = sc.read_h5ad("your_data.h5ad")

# 2. Store raw counts + original layers early (before HVG/normalization)
scat.store_raw_counts(adata, layer="counts", save_raw=False)

# 3. Standard preprocessing (adjust as needed for your analysis)
sc.pp.highly_variable_genes(adata, n_top_genes=3000)
sc.pp.normalize_total(adata, target_sum=1e4)
sc.pp.log1p(adata)
sc.pp.neighbors(adata)
sc.tl.umap(adata)
sc.tl.leiden(adata)

# 4. Attach gene features for bias correction (optional)
adata = scat.add_gene_features(adata, organism="mouse")  # or "human"

# 5. Run differential analysis (active transcription score)
adata_res, significant, all_results = scat.active_score(
    adata_input=adata,
    groupby="condition",
    target_group="Disease",
    reference_group="Control",
    show_plot=False,
)

print("Differential analysis results (top rows):")
print(all_results.head())

# 6. Gene filtering (use the full table; the built-in 'significant' is often empty)
candidates = scat.filter_active_genes(
    all_results,
    preset="heuristic",           # or "pseudobulk" / "permissive"
    # active_score_cutoff=30,
    # logfc_cutoff=0.3,
    # pval_cutoff=0.05,
)

print(f"\nFiltered candidate genes: {len(candidates)}")

# 7. Functional enrichment (GO)
enrich_res = scat.run_enrichment(
    gene_list=candidates.index.tolist(),
    gene_sets="GO_Biological_Process",   # or "GO_BP"
    organism="mouse",                    # or "human"
    adata=adata,                         # uses stored raw genes as background
    pval_cutoff=0.05,
)

print("\nTop GO enrichment terms:")
print(enrich_res.head())

# KEGG enrichment (alternative)
kegg_res = scat.run_kegg(
    gene_list=candidates.index.tolist(),
    organism="mouse",   # or "human"
    adata=adata,
)

# 8. Visualize enrichment results
scat.pl.enrich_dotplot(enrich_res, top_n=15, title="GO Enrichment")
scat.pl.enrich_dotplot(kegg_res, top_n=10, title="KEGG Pathways")

# Optional: save figures
# scat.pl.enrich_dotplot(enrich_res, top_n=12, save_path="enrich_go.pdf")

# Optional: main result plots
# scat.pl.comet_plot(all_results, top_n=12)
# scat.pl.volcano_plot(all_results, top_n=10)

You can now explore all_results, adjust filters in step 6, try different gene_sets, or run run_go / run_gsea (see Functional Enrichment).

For pure differential expression without spliced/unspliced layers, replace step 5 with scat.differential_expression(...) — see Standalone Differential Expression (no velocity data required).

Preserving raw counts and layers

Call store_raw_counts early (after loading and QC, before HVG or normalization). It writes the current .X to layers["counts"] and copies the original spliced/unspliced layers. These survive later subsetting and provide the correct background for enrichment and count-based DE.

The default save_raw=False avoids populating adata.raw.

After HVG-based visualization on a copy, restore or use the preserved layers for full-gene DE, active scoring, or enrichment (pass adata= to run_enrichment or run_kegg to use the stored gene list as background).

HVG subsetting also subsets the saved layers. This keeps velocity calculations consistent with .X. To analyze more genes than the HVG set, store before subsetting or operate on the unfiltered object for DE and enrichment steps.

To restore raw counts into .X for the current gene set:

adata_raw = scat.restore_raw_counts(adata, layer="counts", inplace=False)

See Standalone Differential Expression (no velocity data required) for the no-velocity use case.