Scanpy

#!/usr/bin/env python3
"""
Complete Single-Cell Analysis Template

This template provides a complete workflow for single-cell RNA-seq analysis
using scanpy, from data loading through clustering and cell type annotation.

Customize the parameters and sections as needed for your specific dataset.
"""

import scanpy as sc
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt

# ============================================================================
# CONFIGURATION
# ============================================================================

# File paths
INPUT_FILE = 'data/raw_counts.h5ad'  # Change to your input file
OUTPUT_DIR = 'results/'
FIGURES_DIR = 'figures/'

# QC parameters
MIN_GENES = 200          # Minimum genes per cell
MIN_CELLS = 3            # Minimum cells per gene
MT_THRESHOLD = 5         # Maximum mitochondrial percentage

# Analysis parameters
N_TOP_GENES = 2000       # Number of highly variable genes
N_PCS = 40               # Number of principal components
N_NEIGHBORS = 10         # Number of neighbors for graph
LEIDEN_RESOLUTION = 0.5  # Clustering resolution

# Scanpy settings
sc.settings.verbosity = 3
sc.settings.set_figure_params(dpi=80, facecolor='white')
sc.settings.figdir = FIGURES_DIR
sc.settings.autosave = True

# ============================================================================
# 1. LOAD DATA
# ============================================================================

print("=" * 80)
print("LOADING DATA")
print("=" * 80)

# Load data (adjust based on your file format)
adata = sc.read_h5ad(INPUT_FILE)
# adata = sc.read_10x_mtx('data/filtered_gene_bc_matrices/')  # For 10X data
# adata = sc.read_csv('data/counts.csv')  # For CSV data

print(f"Loaded: {adata.n_obs} cells x {adata.n_vars} genes")

# ============================================================================
# 2. QUALITY CONTROL
# ============================================================================

print("\n" + "=" * 80)
print("QUALITY CONTROL")
print("=" * 80)

# Identify mitochondrial genes
adata.var['mt'] = adata.var_names.str.startswith('MT-')

# Calculate QC metrics
sc.pp.calculate_qc_metrics(adata, qc_vars=['mt'], percent_top=None,
                            log1p=False, inplace=True)

# Visualize QC metrics before filtering
sc.pl.violin(adata, ['n_genes_by_counts', 'total_counts', 'pct_counts_mt'],
             jitter=0.4, multi_panel=True, save='_qc_before_filtering')

sc.pl.scatter(adata, x='total_counts', y='pct_counts_mt', save='_qc_mt')
sc.pl.scatter(adata, x='total_counts', y='n_genes_by_counts', save='_qc_genes')

# Filter cells and genes
print(f"\nBefore filtering: {adata.n_obs} cells, {adata.n_vars} genes")

sc.pp.filter_cells(adata, min_genes=MIN_GENES)
sc.pp.filter_genes(adata, min_cells=MIN_CELLS)
adata = adata[adata.obs.pct_counts_mt < MT_THRESHOLD, :]

print(f"After filtering: {adata.n_obs} cells, {adata.n_vars} genes")

# Optional: doublet detection (uncomment; run before normalization)
# sc.pp.scrublet(adata)
# adata = adata[~adata.obs['predicted_doublet'], :].copy()
# print(f"After doublet removal: {adata.n_obs} cells")

# ============================================================================
# 3. NORMALIZATION
# ============================================================================

print("\n" + "=" * 80)
print("NORMALIZATION")
print("=" * 80)

# Normalize to 10,000 counts per cell
sc.pp.normalize_total(adata, target_sum=1e4)

# Log-transform
sc.pp.log1p(adata)

# Store normalized data
adata.raw = adata

# ============================================================================
# 4. FEATURE SELECTION
# ============================================================================

print("\n" + "=" * 80)
print("FEATURE SELECTION")
print("=" * 80)

# Identify highly variable genes
sc.pp.highly_variable_genes(adata, n_top_genes=N_TOP_GENES)

# Visualize
sc.pl.highly_variable_genes(adata, save='_hvg')

print(f"Selected {sum(adata.var.highly_variable)} highly variable genes")

# Subset to highly variable genes
adata = adata[:, adata.var.highly_variable]

# ============================================================================
# 5. SCALING AND REGRESSION
# ============================================================================

print("\n" + "=" * 80)
print("SCALING AND REGRESSION")
print("=" * 80)

# Regress out unwanted sources of variation
sc.pp.regress_out(adata, ['total_counts', 'pct_counts_mt'])

# Scale data
sc.pp.scale(adata, max_value=10)

# ============================================================================
# 6. DIMENSIONALITY REDUCTION
# ============================================================================

print("\n" + "=" * 80)
print("DIMENSIONALITY REDUCTION")
print("=" * 80)

# PCA
sc.tl.pca(adata, svd_solver='arpack')
sc.pl.pca_variance_ratio(adata, log=True, save='_pca_variance')

# Compute neighborhood graph
sc.pp.neighbors(adata, n_neighbors=N_NEIGHBORS, n_pcs=N_PCS)

# UMAP
sc.tl.umap(adata)

# ============================================================================
# 7. CLUSTERING
# ============================================================================

print("\n" + "=" * 80)
print("CLUSTERING")
print("=" * 80)

# Leiden clustering
sc.tl.leiden(adata, resolution=LEIDEN_RESOLUTION)

# Visualize
sc.pl.umap(adata, color='leiden', legend_loc='on data', save='_leiden')

print(f"Identified {len(adata.obs['leiden'].unique())} clusters")

# ============================================================================
# 8. MARKER GENE IDENTIFICATION
# ============================================================================

print("\n" + "=" * 80)
print("MARKER GENE IDENTIFICATION")
print("=" * 80)

# Find marker genes (exploratory — pseudobulk + pydeseq2 for rigorous DE)
sc.tl.rank_genes_groups(adata, 'leiden', method='wilcoxon')

# Visualize top markers
sc.pl.rank_genes_groups(adata, n_genes=25, sharey=False, save='_markers')
sc.pl.rank_genes_groups_heatmap(adata, n_genes=10, save='_markers_heatmap')
sc.pl.rank_genes_groups_dotplot(adata, n_genes=5, save='_markers_dotplot')

# Get top markers for each cluster
for cluster in adata.obs['leiden'].unique():
    print(f"\nCluster {cluster} top markers:")
    markers = sc.get.rank_genes_groups_df(adata, group=cluster).head(10)
    print(markers[['names', 'scores', 'pvals_adj']].to_string(index=False))

# ============================================================================
# 9. CELL TYPE ANNOTATION (CUSTOMIZE THIS SECTION)
# ============================================================================

print("\n" + "=" * 80)
print("CELL TYPE ANNOTATION")
print("=" * 80)

# Example marker genes for common cell types (customize for your data)
marker_genes = {
    'T cells': ['CD3D', 'CD3E', 'CD3G'],
    'B cells': ['MS4A1', 'CD79A', 'CD79B'],
    'Monocytes': ['CD14', 'LYZ', 'S100A8'],
    'NK cells': ['NKG7', 'GNLY', 'KLRD1'],
    'Dendritic cells': ['FCER1A', 'CST3'],
}

# Visualize marker genes
for cell_type, genes in marker_genes.items():
    available_genes = [g for g in genes if g in adata.raw.var_names]
    if available_genes:
        sc.pl.umap(adata, color=available_genes, use_raw=True,
                   save=f'_{cell_type.replace(" ", "_")}')

# Manual annotation based on marker expression (customize this mapping)
cluster_to_celltype = {
    '0': 'CD4 T cells',
    '1': 'CD14+ Monocytes',
    '2': 'B cells',
    '3': 'CD8 T cells',
    '4': 'NK cells',
    # Add more mappings based on your marker analysis
}

# Apply annotations
adata.obs['cell_type'] = adata.obs['leiden'].map(cluster_to_celltype)
adata.obs['cell_type'] = adata.obs['cell_type'].fillna('Unknown')

# Visualize annotated cell types
sc.pl.umap(adata, color='cell_type', legend_loc='on data', save='_celltypes')

# ============================================================================
# 10. ADDITIONAL ANALYSES (OPTIONAL)
# ============================================================================

print("\n" + "=" * 80)
print("ADDITIONAL ANALYSES")
print("=" * 80)

# PAGA trajectory analysis (optional)
sc.tl.paga(adata, groups='leiden')
sc.pl.paga(adata, color='leiden', save='_paga')

# Gene set scoring (optional)
# example_gene_set = ['CD3D', 'CD3E', 'CD3G']
# sc.tl.score_genes(adata, example_gene_set, score_name='T_cell_score')
# sc.pl.umap(adata, color='T_cell_score', save='_gene_set_score')

# ============================================================================
# 11. SAVE RESULTS
# ============================================================================

print("\n" + "=" * 80)
print("SAVING RESULTS")
print("=" * 80)

import os
os.makedirs(OUTPUT_DIR, exist_ok=True)

# Save processed AnnData object
adata.write(f'{OUTPUT_DIR}/processed_data.h5ad')
print(f"Saved processed data to {OUTPUT_DIR}/processed_data.h5ad")

# Export metadata
adata.obs.to_csv(f'{OUTPUT_DIR}/cell_metadata.csv')
adata.var.to_csv(f'{OUTPUT_DIR}/gene_metadata.csv')
print(f"Saved metadata to {OUTPUT_DIR}/")

# Export marker genes
for cluster in adata.obs['leiden'].unique():
    markers = sc.get.rank_genes_groups_df(adata, group=cluster)
    markers.to_csv(f'{OUTPUT_DIR}/markers_cluster_{cluster}.csv', index=False)
print(f"Saved marker genes to {OUTPUT_DIR}/")

# ============================================================================
# 12. SUMMARY
# ============================================================================

print("\n" + "=" * 80)
print("ANALYSIS SUMMARY")
print("=" * 80)

print(f"\nFinal dataset:")
print(f"  Cells: {adata.n_obs}")
print(f"  Genes: {adata.n_vars}")
print(f"  Clusters: {len(adata.obs['leiden'].unique())}")

print(f"\nCell type distribution:")
print(adata.obs['cell_type'].value_counts())

print("\n" + "=" * 80)
print("ANALYSIS COMPLETE")
print("=" * 80)

{
  "0": "CD4 T cells",
  "1": "CD14+ Monocytes",
  "2": "B cells",
  "3": "CD8 T cells",
  "4": "NK cells",
  "5": "FCGR3A+ Monocytes",
  "6": "Dendritic cells",
  "7": "Platelets"
}

{
  "T_cell": ["CD3D", "CD3E", "CD3G", "CD2"],
  "B_cell": ["MS4A1", "CD79A", "CD79B", "CD19"],
  "monocyte": ["CD14", "LYZ", "S100A8", "S100A9"],
  "NK_cell": ["NKG7", "GNLY", "KLRD1", "NCAM1"],
  "cytotoxic": ["GZMB", "GZMK", "PRF1", "NKG7"],
  "dendritic": ["FCER1A", "CST3", "CLEC10A"],
  "exhaustion": ["PDCD1", "HAVCR2", "LAG3", "TIGIT", "CTLA4"]
}

{
  "min_genes": 200,
  "max_genes": 6000,
  "min_cells": 3,
  "mt_threshold": 10,
  "scrublet": true,
  "target_sum": 10000,
  "n_top_genes": 2000,
  "hvg_flavor": "seurat",
  "scale": false,
  "regress_out": null,
  "n_pcs": 40,
  "n_neighbors": 15,
  "resolution": 0.5,
  "batch_key": null,
  "batch_method": "harmony",
  "marker_method": "wilcoxon",
  "skip_markers": false
}

#!/usr/bin/env python3
"""
Shared helpers for the scanpy script toolkit.

Every CLI script in this directory imports from this module so that data
loading, saving, figure configuration, and logging behave consistently.
This file is NOT a CLI itself; import it:

    from _common import load_anndata, save_anndata, configure_scanpy, info
"""

import os
import sys


def info(msg):
    """Print a progress message to stderr-friendly stdout with a marker."""
    print(f"[scanpy] {msg}", flush=True)


def die(msg, code=1):
    """Print an error and exit."""
    print(f"Error: {msg}", file=sys.stderr, flush=True)
    sys.exit(code)


def _import_scanpy():
    try:
        import scanpy as sc  # noqa: F401
        return sc
    except ImportError:
        die('scanpy not installed. Install with: uv pip install "scanpy[leiden]"')


def configure_scanpy(figdir="figures", dpi=120, verbosity=1, autosave=False,
                     file_format="png"):
    """Apply consistent scanpy settings and return the scanpy module.

    Note: autosave is left off by default because the toolkit scripts pass
    explicit ``save=`` suffixes to plotting calls for predictable filenames.
    """
    sc = _import_scanpy()
    sc.settings.verbosity = verbosity
    sc.settings.set_figure_params(dpi=dpi, facecolor="white")
    sc.settings.figdir = figdir
    sc.settings.file_format_figs = file_format
    if autosave:
        sc.settings.autosave = True
    os.makedirs(figdir, exist_ok=True)
    return sc


def load_anndata(path, var_names="gene_symbols"):
    """Load an AnnData object, dispatching on the file extension / layout.

    Supported inputs:
      * ``.h5ad``                  -> sc.read_h5ad
      * ``.h5`` (10x CellRanger)   -> sc.read_10x_h5
      * ``.csv`` / ``.tsv`` / ``.txt`` -> sc.read_csv / read_text
      * ``.loom``                  -> sc.read_loom
      * ``.mtx``                   -> sc.read (matrix market)
      * a directory                -> sc.read_10x_mtx (10x mtx folder)
    """
    sc = _import_scanpy()
    if not os.path.exists(path):
        die(f"input not found: {path}")

    lower = path.lower()
    if os.path.isdir(path):
        info(f"Reading 10x mtx directory: {path}")
        return sc.read_10x_mtx(path, var_names=var_names)
    if lower.endswith(".h5ad"):
        return sc.read_h5ad(path)
    if lower.endswith(".h5"):
        info("Reading 10x HDF5 (.h5)")
        return sc.read_10x_h5(path)
    if lower.endswith(".loom"):
        return sc.read_loom(path)
    if lower.endswith(".csv"):
        return sc.read_csv(path)
    if lower.endswith((".tsv", ".txt")):
        return sc.read_text(path)
    if lower.endswith(".mtx") or lower.endswith(".mtx.gz"):
        return sc.read(path)
    die(f"unrecognized input format: {path}")


def save_anndata(adata, path):
    """Write an AnnData object to .h5ad, creating parent dirs as needed."""
    parent = os.path.dirname(os.path.abspath(path))
    os.makedirs(parent, exist_ok=True)
    adata.write_h5ad(path)
    info(f"Wrote {path}  ({adata.n_obs} cells x {adata.n_vars} genes)")


def add_io_args(parser, default_output=None):
    """Attach the standard input/output/figdir arguments to an argparse parser."""
    parser.add_argument("input", help="Input file (.h5ad, .h5, .csv, .loom, or 10x mtx dir)")
    parser.add_argument("-o", "--output", default=default_output,
                        help="Output .h5ad path" +
                             (f" (default: {default_output})" if default_output else ""))
    parser.add_argument("--figdir", default="figures",
                        help="Directory for saved figures (default: figures)")
    return parser


def summarize(adata):
    """Return a short human-readable summary string of an AnnData object."""
    lines = [f"{adata.n_obs} cells x {adata.n_vars} genes"]
    if len(adata.obs.columns):
        lines.append("obs: " + ", ".join(adata.obs.columns[:20]))
    if list(adata.obsm.keys()):
        lines.append("obsm: " + ", ".join(adata.obsm.keys()))
    if list(adata.layers.keys()):
        lines.append("layers: " + ", ".join(adata.layers.keys()))
    return "\n".join(lines)

#!/usr/bin/env python3
"""
Annotate clusters with cell-type labels from a mapping file.

Maps a cluster column (e.g. leiden) to cell-type names using a JSON or CSV
mapping, writes the labels into a new obs column, and saves a UMAP and dotplot.

The mapping file is one of:
  * JSON  : {"0": "CD4 T cells", "1": "B cells", ...}
  * CSV    : two columns ``cluster,cell_type``

Optionally provide a marker-gene JSON to draw a reference dotplot that helps
decide the mapping:
  {"T cells": ["CD3D","CD3E"], "B cells": ["MS4A1","CD79A"]}

Examples:
    python annotate.py clustered.h5ad -o annotated.h5ad --mapping celltypes.json
    python annotate.py clustered.h5ad -o annotated.h5ad --mapping map.csv --cluster-key leiden
    python annotate.py clustered.h5ad --markers markers.json --cluster-key leiden   # dotplot only
"""

import argparse
import json
import os

from _common import add_io_args, configure_scanpy, die, info, load_anndata, save_anndata


def load_mapping(path):
    if path.lower().endswith(".json"):
        with open(path) as fh:
            return {str(k): v for k, v in json.load(fh).items()}
    import pandas as pd
    df = pd.read_csv(path)
    if df.shape[1] < 2:
        die("CSV mapping needs at least two columns: cluster,cell_type")
    return {str(k): v for k, v in zip(df.iloc[:, 0], df.iloc[:, 1])}


def main():
    p = argparse.ArgumentParser(description=__doc__,
                                formatter_class=argparse.RawDescriptionHelpFormatter)
    add_io_args(p, default_output="annotated.h5ad")
    p.add_argument("--cluster-key", default="leiden", help="obs column with clusters (default leiden)")
    p.add_argument("--mapping", default=None, help="JSON or CSV cluster->cell_type mapping")
    p.add_argument("--label-key", default="cell_type", help="New obs column name (default cell_type)")
    p.add_argument("--markers", default=None,
                   help="JSON of {cell_type: [genes]} to draw a reference dotplot")
    p.add_argument("--no-plots", action="store_true", help="Skip plots")
    args = p.parse_args()

    sc = configure_scanpy(figdir=args.figdir)
    adata = load_anndata(args.input)
    if args.cluster_key not in adata.obs.columns:
        die(f"cluster key '{args.cluster_key}' not in obs: {list(adata.obs.columns)}")

    if args.markers:
        with open(args.markers) as fh:
            marker_dict = json.load(fh)
        present = {ct: [g for g in genes if g in (adata.raw.var_names if adata.raw is not None else adata.var_names)]
                   for ct, genes in marker_dict.items()}
        present = {ct: g for ct, g in present.items() if g}
        if present and not args.no_plots:
            sc.pl.dotplot(adata, present, groupby=args.cluster_key,
                          use_raw=adata.raw is not None, show=False, save="_marker_reference.png")
            info("Wrote marker reference dotplot")

    if args.mapping:
        mapping = load_mapping(args.mapping)
        adata.obs[args.label_key] = (
            adata.obs[args.cluster_key].astype(str).map(mapping).fillna("Unknown").astype("category")
        )
        info(f"Annotated '{args.label_key}': "
             + ", ".join(f"{k}={v}" for k, v in adata.obs[args.label_key].value_counts().items()))
        if not args.no_plots:
            sc.pl.umap(adata, color=args.label_key, legend_loc="on data",
                       show=False, save="_celltypes.png")
        save_anndata(adata, args.output)
    elif not args.markers:
        die("provide --mapping (to annotate) and/or --markers (for a reference dotplot)")


if __name__ == "__main__":
    main()

#!/usr/bin/env python3
"""
Batch correction / integration across samples.

Supports three methods:
  * harmony  : corrects the PCA embedding -> writes obsm['X_pca_harmony'].
               Fast, recommended default. Needs harmonypy (uv pip install harmonypy).
               Follow with: reduce_dimensions.py --use-rep X_pca_harmony
  * bbknn     : batch-balanced kNN graph (replaces sc.pp.neighbors). Then cluster directly.
               Needs bbknn (uv pip install bbknn).
  * combat    : corrects the expression matrix in place (sc.pp.combat). Built into scanpy.

Run on a normalized object that already has PCA (harmony/bbknn) computed.

Examples:
    python batch_correct.py reduced.h5ad -o integrated.h5ad --method harmony --batch-key sample
    python batch_correct.py reduced.h5ad -o integrated.h5ad --method bbknn --batch-key sample
    python batch_correct.py normalized.h5ad -o integrated.h5ad --method combat --batch-key batch
"""

import argparse

from _common import add_io_args, configure_scanpy, die, info, load_anndata, save_anndata


def main():
    p = argparse.ArgumentParser(description=__doc__,
                                formatter_class=argparse.RawDescriptionHelpFormatter)
    add_io_args(p, default_output="integrated.h5ad")
    p.add_argument("--method", default="harmony", choices=["harmony", "bbknn", "combat"])
    p.add_argument("--batch-key", required=True, help="obs column identifying batches")
    args = p.parse_args()

    sc = configure_scanpy(figdir=args.figdir)
    adata = load_anndata(args.input)
    if args.batch_key not in adata.obs.columns:
        die(f"batch key '{args.batch_key}' not in obs: {list(adata.obs.columns)}")

    if args.method == "harmony":
        if "X_pca" not in adata.obsm:
            sc.tl.pca(adata, svd_solver="arpack")
        try:
            sc.external.pp.harmony_integrate(adata, args.batch_key)
        except ImportError:
            die("harmonypy not installed. Install with: uv pip install harmonypy")
        info("Wrote obsm['X_pca_harmony']. Next: "
             "reduce_dimensions.py --use-rep X_pca_harmony")
    elif args.method == "bbknn":
        if "X_pca" not in adata.obsm:
            sc.tl.pca(adata, svd_solver="arpack")
        try:
            sc.external.pp.bbknn(adata, batch_key=args.batch_key)
        except ImportError:
            die("bbknn not installed. Install with: uv pip install bbknn")
        sc.tl.umap(adata)
        info("Built batch-balanced graph + UMAP. Next: cluster.py")
    elif args.method == "combat":
        sc.pp.combat(adata, key=args.batch_key)
        info("Corrected expression matrix with ComBat. Re-run reduce_dimensions.py.")

    save_anndata(adata, args.output)


if __name__ == "__main__":
    main()

#!/usr/bin/env python3
"""
Leiden clustering on a precomputed neighborhood graph.

Runs Leiden at one or more resolutions and writes a UMAP colored by each
clustering. Requires that ``sc.pp.neighbors`` has already been run
(use reduce_dimensions.py first).

Examples:
    python cluster.py reduced.h5ad -o clustered.h5ad --resolution 0.5
    python cluster.py reduced.h5ad -o clustered.h5ad --resolution 0.3 0.5 0.8 1.0
    python cluster.py reduced.h5ad -o clustered.h5ad --algorithm louvain
"""

import argparse

from _common import add_io_args, configure_scanpy, die, info, load_anndata, save_anndata


def main():
    p = argparse.ArgumentParser(description=__doc__,
                                formatter_class=argparse.RawDescriptionHelpFormatter)
    add_io_args(p, default_output="clustered.h5ad")
    p.add_argument("--resolution", type=float, nargs="+", default=[0.5],
                   help="One or more resolutions (default 0.5). Higher = more clusters")
    p.add_argument("--algorithm", default="leiden", choices=["leiden", "louvain"],
                   help="Clustering algorithm (default leiden)")
    p.add_argument("--key", default=None,
                   help="obs key for the result (single resolution only; "
                        "default '<algorithm>'). Multiple resolutions use '<algorithm>_<res>'")
    p.add_argument("--no-plots", action="store_true", help="Skip UMAP plots")
    args = p.parse_args()

    sc = configure_scanpy(figdir=args.figdir)
    adata = load_anndata(args.input)
    if "neighbors" not in adata.uns:
        die("no neighborhood graph found. Run reduce_dimensions.py first.")

    cluster_fn = sc.tl.leiden if args.algorithm == "leiden" else sc.tl.louvain
    keys = []
    for res in args.resolution:
        if len(args.resolution) == 1:
            key = args.key or args.algorithm
        else:
            key = f"{args.algorithm}_{res}"
        # flavor='igraph' is the scanpy 1.12 default-recommended Leiden backend.
        kwargs = {"resolution": res, "key_added": key}
        if args.algorithm == "leiden":
            kwargs.update(flavor="igraph", n_iterations=2, directed=False)
        cluster_fn(adata, **kwargs)
        n = adata.obs[key].nunique()
        info(f"{key}: {n} clusters at resolution {res}")
        keys.append(key)

    if not args.no_plots:
        sc.pl.umap(adata, color=keys, legend_loc="on data",
                   show=False, save="_clusters.png")

    save_anndata(adata, args.output)


if __name__ == "__main__":
    main()

#!/usr/bin/env python3
"""
Load any supported single-cell format and write it as .h5ad.

Convenience wrapper to get 10x mtx folders, 10x .h5, CSV/TSV, loom, or mtx
files into AnnData .h5ad once, so later steps all read a single fast format.
For R-native files (.rds / Seurat / SingleCellExperiment), see
references/r_interop.md — those must be converted with R first.

Examples:
    python convert.py filtered_feature_bc_matrix/ -o data.h5ad
    python convert.py raw_counts.csv -o data.h5ad --transpose
    python convert.py matrix.h5 -o data.h5ad
"""

import argparse

from _common import add_io_args, configure_scanpy, info, load_anndata, save_anndata, summarize


def main():
    p = argparse.ArgumentParser(description=__doc__,
                                formatter_class=argparse.RawDescriptionHelpFormatter)
    add_io_args(p, default_output="data.h5ad")
    p.add_argument("--transpose", action="store_true",
                   help="Transpose after loading (use if matrix is genes x cells)")
    p.add_argument("--make-unique", action="store_true",
                   help="Make var (gene) names unique")
    args = p.parse_args()

    configure_scanpy(figdir=args.figdir)
    adata = load_anndata(args.input)
    if args.transpose:
        adata = adata.T
        info(f"Transposed -> {adata.n_obs} cells x {adata.n_vars} genes")
    if args.make_unique:
        adata.var_names_make_unique()
    print(summarize(adata))
    save_anndata(adata, args.output)


if __name__ == "__main__":
    main()

#!/usr/bin/env python3
"""
Rank marker genes per group and export tables + plots.

Runs ``sc.tl.rank_genes_groups`` for a grouping (e.g. leiden clusters), writes
a combined CSV of the top markers per group, per-group CSVs, and the standard
marker plots (rank panel, heatmap, dotplot).

NOTE: per-cell tests inflate significance because cells are not independent.
Use this for EXPLORATORY cluster markers. For rigorous DE between conditions,
use pseudobulk.py + pydeseq2.

Examples:
    python find_markers.py clustered.h5ad -o markers.h5ad
    python find_markers.py clustered.h5ad --groupby leiden --method wilcoxon --n-genes 50
    python find_markers.py clustered.h5ad --csv-dir results/markers --top 10
"""

import argparse
import os

from _common import add_io_args, configure_scanpy, die, info, load_anndata, save_anndata


def main():
    p = argparse.ArgumentParser(description=__doc__,
                                formatter_class=argparse.RawDescriptionHelpFormatter)
    add_io_args(p, default_output=None)
    p.add_argument("--groupby", default="leiden", help="obs column to group by (default leiden)")
    p.add_argument("--method", default="wilcoxon",
                   choices=["wilcoxon", "t-test", "t-test_overestim_var", "logreg"],
                   help="Test method (default wilcoxon)")
    p.add_argument("--n-genes", type=int, default=25, help="Genes shown in plots (default 25)")
    p.add_argument("--top", type=int, default=25, help="Top markers per group in CSV (default 25)")
    p.add_argument("--use-raw", action="store_true",
                   help="Rank on adata.raw (normalized log values) instead of X")
    p.add_argument("--csv-dir", default="results/markers", help="Directory for marker CSVs")
    p.add_argument("--no-plots", action="store_true", help="Skip marker plots")
    args = p.parse_args()

    sc = configure_scanpy(figdir=args.figdir)
    adata = load_anndata(args.input)
    if args.groupby not in adata.obs.columns:
        die(f"groupby column '{args.groupby}' not found in obs: {list(adata.obs.columns)}")

    sc.tl.rank_genes_groups(adata, args.groupby, method=args.method,
                            use_raw=args.use_raw if adata.raw is not None else False)

    os.makedirs(args.csv_dir, exist_ok=True)
    groups = list(adata.obs[args.groupby].cat.categories) \
        if hasattr(adata.obs[args.groupby], "cat") else sorted(adata.obs[args.groupby].unique())
    combined = []
    for g in groups:
        df = sc.get.rank_genes_groups_df(adata, group=str(g)).head(args.top)
        df.insert(0, "group", g)
        df.to_csv(os.path.join(args.csv_dir, f"markers_{args.groupby}_{g}.csv"), index=False)
        combined.append(df)
    if combined:
        import pandas as pd
        all_path = os.path.join(args.csv_dir, f"markers_{args.groupby}_all.csv")
        pd.concat(combined, ignore_index=True).to_csv(all_path, index=False)
        info(f"Wrote marker tables to {args.csv_dir}/ ({len(groups)} groups)")

    if not args.no_plots:
        sc.pl.rank_genes_groups(adata, n_genes=args.n_genes, sharey=False,
                                show=False, save="_markers.png")
        sc.pl.rank_genes_groups_dotplot(adata, n_genes=5, show=False, save="_markers_dotplot.png")
        sc.pl.rank_genes_groups_heatmap(adata, n_genes=10, show=False, save="_markers_heatmap.png")

    if args.output:
        save_anndata(adata, args.output)


if __name__ == "__main__":
    main()

#!/usr/bin/env python3
"""
Inspect an AnnData / single-cell file and print a structured summary.

Reports shape, obs/var columns (with dtypes and category counts), layers,
obsm/varm/uns keys, X dtype and value range, and whether the data looks like
raw counts or normalized values. Use this before analysis to understand an
unfamiliar dataset and decide which pipeline steps still need to run.

Examples:
    python inspect_data.py data.h5ad
    python inspect_data.py 10x_dir/
"""

import argparse

import numpy as np

from _common import configure_scanpy, info, load_anndata


def _is_integer_matrix(X, n=10000):
    sub = X[:50] if X.shape[0] > 50 else X
    arr = sub.toarray() if hasattr(sub, "toarray") else np.asarray(sub)
    arr = arr.ravel()[:n]
    if arr.size == 0:
        return False
    return bool(np.all(np.equal(np.mod(arr, 1), 0)))


def main():
    p = argparse.ArgumentParser(description=__doc__,
                                formatter_class=argparse.RawDescriptionHelpFormatter)
    p.add_argument("input", help="Input file (.h5ad, .h5, .csv, .loom, or 10x mtx dir)")
    p.add_argument("--max-cols", type=int, default=40, help="Max obs/var columns to list")
    args = p.parse_args()

    configure_scanpy()
    adata = load_anndata(args.input)

    print("=" * 70)
    print(f"AnnData: {adata.n_obs} cells x {adata.n_vars} genes")
    print("=" * 70)

    X = adata.X
    looks_int = _is_integer_matrix(X)
    xmax = X.max() if not hasattr(X, "toarray") else X.max()
    print(f"\nX dtype={X.dtype}  max={float(xmax):.2f}  "
          f"=> {'raw counts (likely)' if looks_int else 'normalized/log (likely)'}")

    print(f"\nobs columns ({len(adata.obs.columns)}):")
    for c in adata.obs.columns[:args.max_cols]:
        col = adata.obs[c]
        if str(col.dtype) in ("category", "object"):
            nuniq = col.nunique()
            extra = f"  {nuniq} categories" + (f": {list(col.unique()[:8])}" if nuniq <= 8 else "")
        else:
            extra = f"  range=[{col.min():.2f}, {col.max():.2f}]"
        print(f"  - {c} ({col.dtype}){extra}")

    print(f"\nvar columns ({len(adata.var.columns)}): {list(adata.var.columns[:args.max_cols])}")
    print(f"\nlayers: {list(adata.layers.keys())}")
    print(f"obsm: {list(adata.obsm.keys())}")
    print(f"varm: {list(adata.varm.keys())}")
    print(f"uns:  {list(adata.uns.keys())}")
    print(f"raw:  {'present' if adata.raw is not None else 'none'}")

    done = []
    if any(k in adata.obsm for k in ("X_pca",)):
        done.append("PCA")
    if "neighbors" in adata.uns:
        done.append("neighbors")
    if "X_umap" in adata.obsm:
        done.append("UMAP")
    if any(k in adata.obs.columns for k in ("leiden", "louvain")):
        done.append("clustering")
    print("\nPipeline steps already present: " + (", ".join(done) if done else "none"))


if __name__ == "__main__":
    main()

#!/usr/bin/env python3
"""
Generate common single-cell plots from a processed AnnData object.

A flexible plotting front-end so the agent doesn't hand-write matplotlib /
scanpy plotting calls. Pick a plot type and the keys/genes to show.

Plot types:
  umap / tsne / pca   : embedding colored by --color obs columns and/or --genes
  violin               : --genes (or QC metrics) split by --groupby
  dotplot / matrixplot / heatmap / tracksplot / stacked_violin : --genes by --groupby

Examples:
    python plot.py annotated.h5ad --kind umap --color leiden cell_type
    python plot.py annotated.h5ad --kind umap --genes CD3D MS4A1 NKG7 --use-raw
    python plot.py annotated.h5ad --kind dotplot --genes CD3D CD14 MS4A1 --groupby cell_type
    python plot.py annotated.h5ad --kind violin --genes CD3D --groupby leiden
"""

import argparse

from _common import configure_scanpy, die, info, load_anndata


def main():
    p = argparse.ArgumentParser(description=__doc__,
                                formatter_class=argparse.RawDescriptionHelpFormatter)
    p.add_argument("input", help="Input .h5ad")
    p.add_argument("--kind", required=True,
                   choices=["umap", "tsne", "pca", "violin", "dotplot",
                            "matrixplot", "heatmap", "tracksplot", "stacked_violin"])
    p.add_argument("--color", nargs="+", default=None, help="obs columns to color by")
    p.add_argument("--genes", nargs="+", default=None, help="genes to display")
    p.add_argument("--groupby", default=None, help="obs column to group by (violin/dotplot/...)")
    p.add_argument("--use-raw", action="store_true", help="Use adata.raw (normalized log values)")
    p.add_argument("--figdir", default="figures", help="Figure output directory")
    p.add_argument("--save", default=None, help="Filename suffix (default derived from --kind)")
    args = p.parse_args()

    sc = configure_scanpy(figdir=args.figdir)
    adata = load_anndata(args.input)
    # scanpy prepends the plot-function name, so ".png" yields e.g. "umap.png".
    save = args.save or ".png"
    use_raw = args.use_raw and adata.raw is not None

    var_names = adata.raw.var_names if use_raw else adata.var_names
    genes = [g for g in (args.genes or []) if g in var_names]
    if args.genes and len(genes) < len(args.genes):
        missing = set(args.genes) - set(genes)
        info(f"Genes not found (skipped): {sorted(missing)}")

    if args.kind in ("umap", "tsne", "pca"):
        color = (args.color or []) + genes
        if not color:
            die("provide --color and/or --genes for embedding plots")
        fn = {"umap": sc.pl.umap, "tsne": sc.pl.tsne, "pca": sc.pl.pca}[args.kind]
        fn(adata, color=color, use_raw=use_raw, show=False, save=save)
    elif args.kind == "violin":
        keys = genes or [c for c in (args.color or []) if c in adata.obs.columns]
        if not keys:
            die("provide --genes (or obs keys via --color) for violin")
        sc.pl.violin(adata, keys, groupby=args.groupby, use_raw=use_raw,
                     show=False, save=save)
    else:
        if not genes or not args.groupby:
            die(f"{args.kind} requires --genes and --groupby")
        fn = {
            "dotplot": sc.pl.dotplot, "matrixplot": sc.pl.matrixplot,
            "heatmap": sc.pl.heatmap, "tracksplot": sc.pl.tracksplot,
            "stacked_violin": sc.pl.stacked_violin,
        }[args.kind]
        fn(adata, genes, groupby=args.groupby, use_raw=use_raw, show=False, save=save)

    info(f"Saved {args.figdir}/{args.kind}{save}")


if __name__ == "__main__":
    main()

#!/usr/bin/env python3
"""
Normalize, log-transform, select highly variable genes, and optionally scale.

Takes QC-filtered raw counts and produces a normalized, log1p-transformed
object ready for dimensionality reduction. A copy of the raw counts is kept in
``adata.layers['counts']`` and the normalized log values in ``adata.raw`` so
downstream marker/expression plots can use ``use_raw=True``.

Examples:
    python preprocess.py filtered.h5ad -o normalized.h5ad
    python preprocess.py filtered.h5ad -o normalized.h5ad --n-top-genes 3000 --scale
    python preprocess.py filtered.h5ad -o normalized.h5ad --flavor seurat_v3 --batch-key sample
    python preprocess.py filtered.h5ad -o normalized.h5ad --regress-out total_counts pct_counts_mt
"""

import argparse

from _common import add_io_args, configure_scanpy, info, load_anndata, save_anndata


def main():
    p = argparse.ArgumentParser(description=__doc__,
                                formatter_class=argparse.RawDescriptionHelpFormatter)
    add_io_args(p, default_output="normalized.h5ad")
    p.add_argument("--target-sum", type=float, default=1e4,
                   help="Counts per cell after normalization (default 1e4)")
    p.add_argument("--n-top-genes", type=int, default=2000,
                   help="Number of highly variable genes (default 2000)")
    p.add_argument("--flavor", default="seurat",
                   choices=["seurat", "cell_ranger", "seurat_v3"],
                   help="HVG flavor. seurat_v3 expects raw counts (default seurat)")
    p.add_argument("--batch-key", default=None,
                   help="obs column for batch-aware HVG selection")
    p.add_argument("--subset-hvg", action="store_true",
                   help="Subset the matrix to HVGs (smaller, but drops other genes from X)")
    p.add_argument("--regress-out", nargs="+", default=None,
                   help="obs columns to regress out (e.g. total_counts pct_counts_mt)")
    p.add_argument("--scale", action="store_true",
                   help="Scale to unit variance and zero mean (max_value=10)")
    p.add_argument("--no-plots", action="store_true", help="Skip HVG plot")
    args = p.parse_args()

    sc = configure_scanpy(figdir=args.figdir)
    adata = load_anndata(args.input)
    info(f"Loaded {adata.n_obs} cells x {adata.n_vars} genes")

    # Preserve raw counts in a dedicated layer for pseudobulk / DE later.
    adata.layers["counts"] = adata.X.copy()

    if args.flavor == "seurat_v3":
        # seurat_v3 selects HVGs on raw counts, before normalization.
        sc.pp.highly_variable_genes(adata, n_top_genes=args.n_top_genes,
                                    flavor="seurat_v3", batch_key=args.batch_key)
        sc.pp.normalize_total(adata, target_sum=args.target_sum)
        sc.pp.log1p(adata)
    else:
        sc.pp.normalize_total(adata, target_sum=args.target_sum)
        sc.pp.log1p(adata)
        sc.pp.highly_variable_genes(adata, n_top_genes=args.n_top_genes,
                                    flavor=args.flavor, batch_key=args.batch_key)

    n_hvg = int(adata.var["highly_variable"].sum())
    info(f"Selected {n_hvg} highly variable genes")

    # Stash the full normalized log matrix so plots can use_raw=True.
    adata.raw = adata

    if not args.no_plots:
        sc.pl.highly_variable_genes(adata, show=False, save="_hvg.png")

    if args.subset_hvg:
        adata = adata[:, adata.var["highly_variable"]].copy()
        info(f"Subset to {adata.n_vars} HVGs")

    if args.regress_out:
        info(f"Regressing out: {', '.join(args.regress_out)}")
        sc.pp.regress_out(adata, args.regress_out)

    if args.scale:
        info("Scaling (max_value=10)")
        sc.pp.scale(adata, max_value=10)

    save_anndata(adata, args.output)


if __name__ == "__main__":
    main()

#!/usr/bin/env python3
"""
Aggregate single cells into pseudobulk profiles for rigorous DE.

Sums raw counts within each combination of grouping columns (e.g. sample x
cell_type) using ``sc.get.aggregate`` and exports a genes x pseudobulk-sample
count matrix plus a sample-metadata table. Feed these into pydeseq2 / edgeR /
limma for condition comparisons — this is the statistically correct route,
unlike per-cell rank_genes_groups.

Examples:
    python pseudobulk.py annotated.h5ad --by sample cell_type --out-prefix results/pb
    python pseudobulk.py annotated.h5ad --by sample cell_type --layer counts --out-prefix pb
"""

import argparse
import os

from _common import configure_scanpy, die, info, load_anndata


def main():
    p = argparse.ArgumentParser(description=__doc__,
                                formatter_class=argparse.RawDescriptionHelpFormatter)
    p.add_argument("input", help="Input .h5ad")
    p.add_argument("--by", nargs="+", required=True,
                   help="obs columns to aggregate by (e.g. sample cell_type)")
    p.add_argument("--layer", default="counts",
                   help="Layer with raw counts to sum (default 'counts'; falls back to X)")
    p.add_argument("--func", default="sum", choices=["sum", "mean", "count_nonzero"],
                   help="Aggregation function (default sum)")
    p.add_argument("--out-prefix", default="pseudobulk",
                   help="Output prefix; writes <prefix>_counts.csv and <prefix>_samples.csv")
    args = p.parse_args()

    sc = configure_scanpy()
    adata = load_anndata(args.input)
    for col in args.by:
        if col not in adata.obs.columns:
            die(f"grouping column '{col}' not in obs: {list(adata.obs.columns)}")

    layer = args.layer if args.layer in adata.layers else None
    if args.layer and layer is None:
        info(f"Layer '{args.layer}' not found; aggregating adata.X instead")

    pb = sc.get.aggregate(adata, by=args.by, func=args.func, layer=layer)
    # aggregate stores the result in a layer named after func.
    mat = pb.layers[args.func]

    import pandas as pd
    sample_ids = ["_".join(str(pb.obs.iloc[i][c]) for c in args.by) for i in range(pb.n_obs)]
    counts_df = pd.DataFrame(
        mat.T.toarray() if hasattr(mat, "toarray") else mat.T,
        index=pb.var_names, columns=sample_ids,
    )

    parent = os.path.dirname(os.path.abspath(args.out_prefix))
    if parent:
        os.makedirs(parent, exist_ok=True)
    counts_path = f"{args.out_prefix}_counts.csv"
    samples_path = f"{args.out_prefix}_samples.csv"
    counts_df.to_csv(counts_path)
    meta = pb.obs[args.by].copy()
    meta.index = sample_ids
    meta.to_csv(samples_path)

    info(f"Wrote {counts_df.shape[0]} genes x {counts_df.shape[1]} pseudobulk samples")
    info(f"  counts:  {counts_path}")
    info(f"  samples: {samples_path}")
    info("Next: load these into pydeseq2 for differential expression.")


if __name__ == "__main__":
    main()

#!/usr/bin/env python3
"""
Quality control and filtering for single-cell RNA-seq data.

Calculates QC metrics (genes/counts per cell, mitochondrial / ribosomal /
hemoglobin fractions), writes before/after QC plots, optionally runs Scrublet
doublet detection, and filters cells and genes by the given thresholds.

Run this FIRST on raw counts, before normalization.

Examples:
    python qc_analysis.py raw.h5ad -o filtered.h5ad
    python qc_analysis.py raw.h5ad -o filtered.h5ad --mt-threshold 10 --min-genes 500
    python qc_analysis.py 10x_dir/ -o filtered.h5ad --max-genes 6000 --scrublet
"""

import argparse

from _common import add_io_args, configure_scanpy, info, load_anndata, save_anndata


def annotate_gene_classes(adata):
    """Flag mitochondrial, ribosomal, and hemoglobin genes (human or mouse names)."""
    names = adata.var_names
    adata.var["mt"] = names.str.startswith(("MT-", "mt-", "Mt-"))
    adata.var["ribo"] = names.str.startswith(("RPS", "RPL", "Rps", "Rpl"))
    adata.var["hb"] = names.str.contains(r"^HB[^P]|^Hb[^p]", regex=True)
    return [v for v in ["mt", "ribo", "hb"] if adata.var[v].any()]


def make_qc_plots(sc, adata, prefix):
    qc_keys = ["n_genes_by_counts", "total_counts", "pct_counts_mt"]
    qc_keys = [k for k in qc_keys if k in adata.obs.columns]
    sc.pl.violin(adata, qc_keys, jitter=0.4, multi_panel=True,
                 show=False, save=f"_{prefix}_violin.png")
    if "pct_counts_mt" in adata.obs.columns:
        sc.pl.scatter(adata, x="total_counts", y="pct_counts_mt",
                      show=False, save=f"_{prefix}_mt.png")
    sc.pl.scatter(adata, x="total_counts", y="n_genes_by_counts",
                  show=False, save=f"_{prefix}_counts.png")


def main():
    p = argparse.ArgumentParser(description=__doc__,
                                formatter_class=argparse.RawDescriptionHelpFormatter)
    add_io_args(p, default_output="qc_filtered.h5ad")
    p.add_argument("--min-genes", type=int, default=200, help="Min genes per cell (default 200)")
    p.add_argument("--max-genes", type=int, default=None, help="Max genes per cell (upper outliers)")
    p.add_argument("--min-counts", type=int, default=None, help="Min total counts per cell")
    p.add_argument("--max-counts", type=int, default=None, help="Max total counts per cell")
    p.add_argument("--min-cells", type=int, default=3, help="Min cells per gene (default 3)")
    p.add_argument("--mt-threshold", type=float, default=5, help="Max pct mitochondrial counts (default 5)")
    p.add_argument("--scrublet", action="store_true", help="Run Scrublet doublet detection and drop doublets")
    p.add_argument("--no-plots", action="store_true", help="Skip QC plots")
    args = p.parse_args()

    sc = configure_scanpy(figdir=args.figdir)
    adata = load_anndata(args.input)
    adata.var_names_make_unique()
    info(f"Loaded {adata.n_obs} cells x {adata.n_vars} genes")

    qc_vars = annotate_gene_classes(adata)
    sc.pp.calculate_qc_metrics(adata, qc_vars=qc_vars, percent_top=None,
                               log1p=False, inplace=True)
    info(f"Mean genes/cell={adata.obs['n_genes_by_counts'].mean():.0f}  "
         f"mean counts/cell={adata.obs['total_counts'].mean():.0f}  "
         f"mean pct_mt={adata.obs.get('pct_counts_mt', 0).mean():.1f}")

    if not args.no_plots:
        make_qc_plots(sc, adata, "qc_before")

    n0, g0 = adata.n_obs, adata.n_vars
    sc.pp.filter_cells(adata, min_genes=args.min_genes)
    if args.min_counts:
        sc.pp.filter_cells(adata, min_counts=args.min_counts)
    if args.max_genes:
        adata = adata[adata.obs["n_genes_by_counts"] < args.max_genes, :].copy()
    if args.max_counts:
        adata = adata[adata.obs["total_counts"] < args.max_counts, :].copy()
    if "pct_counts_mt" in adata.obs.columns:
        adata = adata[adata.obs["pct_counts_mt"] < args.mt_threshold, :].copy()
    sc.pp.filter_genes(adata, min_cells=args.min_cells)

    if args.scrublet:
        info("Running Scrublet doublet detection...")
        try:
            sc.pp.scrublet(adata)
        except (ImportError, ValueError) as e:
            die(f"Scrublet failed ({e}). Install with: uv pip install scikit-image")
        n_dbl = int(adata.obs["predicted_doublet"].sum())
        adata = adata[~adata.obs["predicted_doublet"], :].copy()
        info(f"Removed {n_dbl} predicted doublets")

    info(f"Cells {n0} -> {adata.n_obs} ({adata.n_obs / n0 * 100:.1f}% kept)  "
         f"Genes {g0} -> {adata.n_vars}")

    if not args.no_plots:
        make_qc_plots(sc, adata, "qc_after")

    save_anndata(adata, args.output)


if __name__ == "__main__":
    main()

#!/usr/bin/env python3
"""
PCA, neighborhood graph, and UMAP / t-SNE embeddings.

Takes a normalized object (optionally HVG-subset / scaled) and computes PCA,
the kNN graph, and a UMAP (and optionally t-SNE) embedding. Writes a PCA
variance-ratio elbow plot to help choose ``--n-pcs``.

Examples:
    python reduce_dimensions.py normalized.h5ad -o reduced.h5ad
    python reduce_dimensions.py normalized.h5ad -o reduced.h5ad --n-pcs 50 --n-neighbors 15
    python reduce_dimensions.py normalized.h5ad -o reduced.h5ad --tsne --use-rep X_pca_harmony
"""

import argparse

from _common import add_io_args, configure_scanpy, info, load_anndata, save_anndata


def main():
    p = argparse.ArgumentParser(description=__doc__,
                                formatter_class=argparse.RawDescriptionHelpFormatter)
    add_io_args(p, default_output="reduced.h5ad")
    p.add_argument("--n-comps", type=int, default=50, help="PCs to compute (default 50)")
    p.add_argument("--n-pcs", type=int, default=40, help="PCs used for the kNN graph (default 40)")
    p.add_argument("--n-neighbors", type=int, default=15, help="Neighbors for the graph (default 15)")
    p.add_argument("--use-rep", default=None,
                   help="obsm key to build the graph from instead of PCA "
                        "(e.g. X_pca_harmony from batch_correct.py)")
    p.add_argument("--tsne", action="store_true", help="Also compute t-SNE")
    p.add_argument("--color", nargs="+", default=None,
                   help="obs/var keys to color the UMAP by (e.g. sample n_genes_by_counts)")
    p.add_argument("--no-plots", action="store_true", help="Skip plots")
    args = p.parse_args()

    sc = configure_scanpy(figdir=args.figdir)
    adata = load_anndata(args.input)
    info(f"Loaded {adata.n_obs} cells x {adata.n_vars} genes")

    n_comps = min(args.n_comps, adata.n_vars - 1, adata.n_obs - 1)
    sc.tl.pca(adata, n_comps=n_comps, svd_solver="arpack")
    if not args.no_plots:
        sc.pl.pca_variance_ratio(adata, n_pcs=n_comps, log=True,
                                 show=False, save="_variance.png")

    sc.pp.neighbors(adata, n_neighbors=args.n_neighbors, n_pcs=args.n_pcs,
                    use_rep=args.use_rep)
    info("Computing UMAP...")
    sc.tl.umap(adata)

    if args.tsne:
        info("Computing t-SNE...")
        sc.tl.tsne(adata, use_rep=args.use_rep or "X_pca")

    if not args.no_plots and args.color:
        color = [c for c in args.color if c in adata.obs.columns or c in adata.var_names]
        if color:
            sc.pl.umap(adata, color=color, show=False, save="_colored.png")

    save_anndata(adata, args.output)


if __name__ == "__main__":
    main()

#!/usr/bin/env python3
"""
End-to-end standard scRNA-seq pipeline in one command.

Runs the full exploratory workflow on raw counts:
    load -> QC + filter -> (optional doublets) -> normalize + log1p -> HVG
    -> (optional scale/regress) -> PCA -> (optional batch correction)
    -> neighbors -> UMAP -> Leiden -> marker genes -> save.

Produces a processed .h5ad, marker CSVs, and figures. Tune the common knobs via
flags, or pass a JSON config with ``--config`` (keys mirror the flag names with
underscores). This is the fastest path from counts to a clustered, annotated-
ready object; use the individual step scripts when you need to iterate on one stage.

Examples:
    python run_pipeline.py raw.h5ad -o processed.h5ad
    python run_pipeline.py raw.h5ad -o processed.h5ad --resolution 0.8 --n-top-genes 3000
    python run_pipeline.py raw.h5ad -o processed.h5ad --batch-key sample --batch-method harmony
    python run_pipeline.py raw.h5ad -o processed.h5ad --config params.json
"""

import argparse
import json
import os

from _common import configure_scanpy, info, load_anndata, save_anndata


def build_parser():
    p = argparse.ArgumentParser(description=__doc__,
                                formatter_class=argparse.RawDescriptionHelpFormatter)
    p.add_argument("input", help="Input raw-counts file (.h5ad, .h5, .csv, 10x dir, ...)")
    p.add_argument("-o", "--output", default="processed.h5ad", help="Output .h5ad")
    p.add_argument("--figdir", default="figures", help="Figure directory")
    p.add_argument("--marker-dir", default="results/markers", help="Marker CSV directory")
    p.add_argument("--config", default=None, help="JSON file overriding any of the options below")
    # QC
    p.add_argument("--min-genes", type=int, default=200)
    p.add_argument("--max-genes", type=int, default=None)
    p.add_argument("--min-cells", type=int, default=3)
    p.add_argument("--mt-threshold", type=float, default=5)
    p.add_argument("--scrublet", action="store_true")
    # normalization / HVG
    p.add_argument("--target-sum", type=float, default=1e4)
    p.add_argument("--n-top-genes", type=int, default=2000)
    p.add_argument("--hvg-flavor", default="seurat", choices=["seurat", "cell_ranger", "seurat_v3"])
    p.add_argument("--scale", action="store_true")
    p.add_argument("--regress-out", nargs="+", default=None)
    # dim reduction / clustering
    p.add_argument("--n-pcs", type=int, default=40)
    p.add_argument("--n-neighbors", type=int, default=15)
    p.add_argument("--resolution", type=float, default=0.5)
    # batch correction
    p.add_argument("--batch-key", default=None)
    p.add_argument("--batch-method", default="harmony", choices=["harmony", "combat"])
    # markers
    p.add_argument("--marker-method", default="wilcoxon")
    p.add_argument("--skip-markers", action="store_true")
    return p


def apply_config(args):
    if args.config:
        with open(args.config) as fh:
            cfg = json.load(fh)
        for k, v in cfg.items():
            setattr(args, k.replace("-", "_"), v)
    return args


def main():
    args = apply_config(build_parser().parse_args())
    sc = configure_scanpy(figdir=args.figdir)

    info("[1/8] Loading data")
    adata = load_anndata(args.input)
    adata.var_names_make_unique()
    info(f"      {adata.n_obs} cells x {adata.n_vars} genes")

    info("[2/8] QC + filtering")
    adata.var["mt"] = adata.var_names.str.startswith(("MT-", "mt-", "Mt-"))
    qc_vars = ["mt"] if adata.var["mt"].any() else []
    sc.pp.calculate_qc_metrics(adata, qc_vars=qc_vars, percent_top=None, log1p=False, inplace=True)
    sc.pl.violin(adata, [k for k in ["n_genes_by_counts", "total_counts", "pct_counts_mt"]
                         if k in adata.obs.columns],
                 jitter=0.4, multi_panel=True, show=False, save="_qc.png")
    n0 = adata.n_obs
    sc.pp.filter_cells(adata, min_genes=args.min_genes)
    sc.pp.filter_genes(adata, min_cells=args.min_cells)
    if args.max_genes:
        adata = adata[adata.obs["n_genes_by_counts"] < args.max_genes, :].copy()
    if qc_vars and "pct_counts_mt" in adata.obs.columns:
        adata = adata[adata.obs["pct_counts_mt"] < args.mt_threshold, :].copy()
    if args.scrublet:
        try:
            sc.pp.scrublet(adata)
            adata = adata[~adata.obs["predicted_doublet"], :].copy()
        except (ImportError, ValueError) as e:
            info(f"      Scrublet skipped ({e}); install scikit-image to enable")
    info(f"      {n0} -> {adata.n_obs} cells, {adata.n_vars} genes")

    info("[3/8] Normalize + log1p + HVG")
    adata.layers["counts"] = adata.X.copy()
    if args.hvg_flavor == "seurat_v3":
        sc.pp.highly_variable_genes(adata, n_top_genes=args.n_top_genes,
                                    flavor="seurat_v3", batch_key=args.batch_key)
        sc.pp.normalize_total(adata, target_sum=args.target_sum)
        sc.pp.log1p(adata)
    else:
        sc.pp.normalize_total(adata, target_sum=args.target_sum)
        sc.pp.log1p(adata)
        sc.pp.highly_variable_genes(adata, n_top_genes=args.n_top_genes,
                                    flavor=args.hvg_flavor, batch_key=args.batch_key)
    adata.raw = adata
    info(f"      {int(adata.var['highly_variable'].sum())} HVGs")

    info("[4/8] Scale / regress (optional)")
    work = adata[:, adata.var["highly_variable"]].copy()
    if args.regress_out:
        sc.pp.regress_out(work, args.regress_out)
    if args.scale:
        sc.pp.scale(work, max_value=10)

    info("[5/8] PCA" + (f" + {args.batch_method} batch correction" if args.batch_key else ""))
    sc.tl.pca(work, svd_solver="arpack")
    sc.pl.pca_variance_ratio(work, log=True, show=False, save="_variance.png")
    use_rep = "X_pca"
    if args.batch_key:
        if args.batch_method == "harmony":
            try:
                sc.external.pp.harmony_integrate(work, args.batch_key)
                use_rep = "X_pca_harmony"
            except ImportError:
                info("      harmonypy not installed; skipping (uv pip install harmonypy)")
        else:
            sc.pp.combat(work, key=args.batch_key)
            sc.tl.pca(work, svd_solver="arpack")

    info("[6/8] Neighbors + UMAP")
    sc.pp.neighbors(work, n_neighbors=args.n_neighbors, n_pcs=args.n_pcs, use_rep=use_rep)
    sc.tl.umap(work)

    info("[7/8] Leiden clustering")
    sc.tl.leiden(work, resolution=args.resolution, flavor="igraph",
                 n_iterations=2, directed=False)
    n_clusters = work.obs["leiden"].nunique()
    color = ["leiden"] + ([args.batch_key] if args.batch_key else [])
    sc.pl.umap(work, color=color, legend_loc="on data", show=False, save="_leiden.png")
    info(f"      {n_clusters} clusters at resolution {args.resolution}")

    # Carry embeddings/clusters back onto the full-gene object so markers use all genes.
    adata.obs["leiden"] = work.obs["leiden"].values
    adata.obsm["X_pca"] = work.obsm["X_pca"]
    adata.obsm["X_umap"] = work.obsm["X_umap"]
    if use_rep in work.obsm:
        adata.obsm[use_rep] = work.obsm[use_rep]
    adata.uns["neighbors"] = work.uns["neighbors"]
    adata.obsp = work.obsp

    if not args.skip_markers:
        info("[8/8] Marker genes")
        sc.tl.rank_genes_groups(adata, "leiden", method=args.marker_method, use_raw=True)
        sc.pl.rank_genes_groups_dotplot(adata, n_genes=5, show=False, save="_markers_dotplot.png")
        os.makedirs(args.marker_dir, exist_ok=True)
        import pandas as pd
        frames = []
        for g in adata.obs["leiden"].cat.categories:
            df = sc.get.rank_genes_groups_df(adata, group=g).head(25)
            df.insert(0, "cluster", g)
            frames.append(df)
        pd.concat(frames, ignore_index=True).to_csv(
            os.path.join(args.marker_dir, "markers_all.csv"), index=False)
        info(f"      marker tables in {args.marker_dir}/")
    else:
        info("[8/8] Skipping markers")

    save_anndata(adata, args.output)
    info("Pipeline complete. Next: inspect markers, then annotate.py with a cluster->cell_type mapping.")


if __name__ == "__main__":
    main()