# Example per-sample QC thresholds — copy to project utils/qc_thresholds.yaml
# Used by scripts/01-qc/03-filter_cells.*

default:
  min_genes: 200
  max_genes: 6000
  min_UMIs: 400
  max_UMIs: 100000
  max_mito_frac: 0.20
  max_hbb_score: 0.05
  max_ambient_frac: 0.25
  remove_predicted_doublet: true

# Optional per-sample overrides
samples:
  SAMPLE_A:
    max_mito_frac: 0.25
  SAMPLE_B:
    min_genes: 300

# MAD mode (alternative to fixed cutoffs)
mad:
  enabled: false
  nmads: 3
  metrics:
    mito_frac: upper
    n_genes: both
    n_UMIs: lower
  caps:
    max_mito_frac: 0.25
    max_genes: 8000
    min_UMIs: 300

# Filtering Strategies

**Human confirms thresholds before any filter script runs.** This file helps propose options; it does not authorize automatic filtering.

Choose explicitly **with the user** after reviewing diagnosis figures. Document agreed thresholds in `qc_thresholds.yaml` or a visible script block.

## Strategy Comparison

| Strategy | Pros | Cons | Best for |
|----------|------|------|----------|
| Fixed global | reproducible, comparable across studies | fails when depth varies by sample | homogeneous 10X cohort, reproducing published atlas |
| Per-sample fixed | handles known protocol differences | needs manual tuning per cohort | mixed sorting strategies |
| Per-sample MAD | adaptive to depth/outlier structure | less comparable across projects | public data, multi-lab merge |
| HQ gate + advanced | stable decontX/doublet inputs | removes low-count cells early | BM, blood-rich, high doublet rate |

## Recommended Default (new projects)

Per-sample **MAD-based** on core metrics + **per-sample Scrublet/DoubletFinder** on HQ cells — **propose first, apply only after sign-off**.

```text
HQ gate: n_genes >= 200 (or 500 for decontX)
mito: mito_frac < median + 3*MAD  (cap at 20-25% for lung)
n_genes high: n_genes < median + 3*MAD  (cap at 6000-8000)
n_UMIs low: n_UMIs >= median - 3*MAD  (floor at 300-500)
doublet: predicted_doublet == False (per sample)
optional: hbb_score < 0.05 for lung/parenchyma
```

Always compute loss rate per sample before applying.

## Fixed Cutoffs (legacy atlas style)

Example from lung merge-first pipeline:

```python
min_genes = 180
max_genes = 6000
min_counts = 400
max_counts = 100000
max_pct_mito = 20
```

Use only when reproducing that atlas or when cohort is pre-validated homogeneous.

## Per-Sample MAD (R sketch)

```r
mad_filter <- function(x, nmads = 3, type = c("lower", "upper", "both")) {
  med <- median(x, na.rm = TRUE)
  mad_val <- mad(x, na.rm = TRUE)
  if (mad_val == 0) mad_val <- sd(x, na.rm = TRUE) / 1.4826
  type <- match.arg(type)
  if (type == "lower") return(x >= med - nmads * mad_val)
  if (type == "upper") return(x <= med + nmads * mad_val)
  (x >= med - nmads * mad_val) & (x <= med + nmads * mad_val)
}
```

Apply within each `sample_id` group.

## Per-Sample MAD (Python sketch)

```python
import numpy as np
import pandas as pd

def mad_bounds(series, nmads=3):
    med = np.median(series)
    mad = np.median(np.abs(series - med))
    if mad == 0:
        mad = series.std(ddof=0) / 1.4826
    return med - nmads * mad, med + nmads * mad

def per_sample_mito_pass(df, sample_col="sample_id", mito_col="mito_frac", nmads=3, cap=0.25):
    out = []
    for sid, g in df.groupby(sample_col):
        lo, hi = mad_bounds(g[mito_col], nmads=nmads)
        hi = min(hi, cap)
        out.append(g[mito_col] <= hi)
    return pd.concat(out).sort_index()
```

## Filter Order

```text
1. optional coarse gate (n_genes >= 200)
2. mito_frac
3. n_genes / n_UMIs (low and high)
4. tissue contamination (hbb_score, etc.)
5. ambient_frac (if computed)
6. doublet removal (per sample)
7. filter_genes min_cells (on retained cells)
```

Do not normalize or scale before filtering.

## Doublet Calling

| Method | Language | Scope |
|--------|----------|-------|
| Scrublet | Python | per sample on HQ cells |
| DoubletFinder | R | per sample on HQ cells |
| scDblFinder | R/Bioc | per sample SCE |

Never run doublet detection on merged multi-sample objects unless samples are technical replicates of the same library.

## Gene Filtering

After cell filtering:

```text
min_cells = 10   # common default
```

Recompute on per-sample filtered matrix before merge. Gene sets may differ slightly per sample — expected.

## Reporting

Save per sample:

```text
cells_before
cells_after
pct_removed
median_n_genes_before/after
median_n_UMIs_before/after
thresholds_used (JSON/YAML)
```

Aggregate to `QC_summary.tsv` at project root or `result/01-qc/QC_summary.tsv`.

# Human-in-the-Loop QC

This reference defines how agents and scripts should behave: **data informs, humans decide**.

## Core Principle

QC is not a button to press. It is an iterative conversation between:

```text
data inspection  ->  human judgment  ->  small reversible action  ->  re-inspect
```

Bundled scripts and project scripts compute and visualize. They do **not** replace the analyst's review of distributions, biology, and sample metadata.

## What the Human Owns

These decisions must come from the analyst, not from agent defaults:

| Decision | Why human |
|----------|-----------|
| which metrics to compute | tissue, species, available side data differ |
| whether a sample looks usable | metadata errors, failed libraries, batch disasters |
| filter thresholds per sample | depth and mito distributions vary |
| whether to remove doublets aggressively | tumor/high-RNA samples behave differently |
| whether ambient correction is needed | cannot infer from code alone |
| when to merge | only after per-sample sign-off |
| exceptions for outlier samples | biology, not statistics, may justify keeping cells |

The agent's job: **surface evidence and options**. The analyst's job: **confirm or override**.

## Iterative Loop (Default)

Do not plan a 4-script end-to-end run unless the user explicitly asks.

```text
Step A  Inspect inputs (matrix type, species, sample list, metadata)
   ↓
Step B  Compute metrics for 1–3 pilot samples OR summarize existing metadata
   ↓
Step C  Present findings table + suggested issues (stop here)
   ↓  ← human confirms metrics scope, pilot samples, next metrics
Step D  Generate diagnosis figures for confirmed samples
   ↓
Step E  Propose threshold *candidates* with expected cell loss (stop here)
   ↓  ← human confirms or edits thresholds per sample
Step F  Apply filter only after explicit approval; save checkpoint
   ↓
Step G  Re-summarize; human decides merge / doublet / ambient next steps
```

Each arrow labeled "stop here" is a **hard gate**. Do not skip gates to "save time".

## Agent Behavior Rules

When this skill is active:

1. **Inspect before coding** — read metadata, check one sample's distributions, report facts.
2. **Pilot first** — default to 1–3 representative samples (good / borderline / bad if known), not the full cohort.
3. **Propose, don't impose** — give threshold *ranges* and predicted loss; never silently pick cutoffs.
4. **Ask structured questions** — use `AskQuestion` or a short numbered list when choices are real forks.
5. **No silent heavy runs** — do not run filter, Scrublet, decontX, or merge on full data unless user says so.
6. **Show numbers** — median, MAD, percentiles, cells removed per sample; not just "looks fine".
7. **Flag conflicts** — e.g. chrY high but metadata says female; high hbb in lung parenchyma.
8. **Preserve reversibility** — keep pre-filter `metadata.tsv` and counts; filtering writes new artifacts.
9. **Script = one stage** — write or run `01-calculate_metrics` first; wait before `03-filter_cells`.
10. **Owner-editable output** — generated scripts must have visible threshold blocks the user can tweak.

## Presentation Template

After metrics computation, report in this shape:

```markdown
## QC 初检 — <sample_id>

**数据**: raw 10x | n_cells=... | median n_genes=... | median mito=...

**观察**:
- ...
- ...

**需你确认**:
1. 是否继续算 <metric>（如 hbb / ambient / doublet）？
2. 该样本是否纳入后续过滤？
3. mito 上限倾向：固定 20% / 按 MAD / 其他？

**若采用建议阈值** (仅估算，未执行):
| 规则 | 预计保留 | 预计剔除 |
```

After user confirms, run the next small step only.

## Threshold Proposal (Not Application)

When suggesting cutoffs, show per sample:

```text
sample_id | median_genes | p95_mito | MAD_mito_hi | fixed_20%_loss | mad_loss
```

Let the user pick column strategy per sample. Heterogeneous cohorts often need **different** rules per sample.

Read `references/filtering-strategies.md` for methods; present as options with trade-offs.

## Scripts Are Tools, Not Pipelines

| Script | Agent may run without asking | Requires human confirmation first |
|--------|------------------------------|-----------------------------------|
| `inspect_qc_metadata.py` | yes (read-only summary) | — |
| `calculate_metrics.py/R` (core only) | pilot samples only | full cohort |
| `calculate_metrics` + `--run-scrublet` | no | yes |
| any `filter_cells` script | no | yes |
| merge QC-passed | no | yes |

## Red Flags — Stop and Ask

- Single sample with >50% predicted loss under proposed thresholds
- Mito median >15% for nominally healthy parenchyma
- chrY/metadata sex mismatch
- hbb_score high in non-blood tissue
- Orders-of-magnitude depth difference vs sibling samples
- User said "run QC" but did not specify tissue, species, or matrix type

## What "Done" Means

QC is not "done" when a filtered h5ad exists. It is done when the **analyst has reviewed** diagnosis outputs and **explicitly accepted** the filter summary for each sample (or documented exceptions).

Minimal sign-off artifact: `QC_summary.tsv` + user message or `QC_decisions.md` noting per-sample thresholds used.

# QC Metrics Catalog

Canonical column names for cross-language workflows. Map Seurat/scanpy native names to these when saving `metadata.tsv`.

## Core Metrics (always compute)

| Column | Definition | Human gene prefix | Mouse gene prefix |
|--------|------------|-------------------|-------------------|
| `n_genes` | genes with count > 0 per cell | — | — |
| `n_UMIs` | total UMI/counts per cell | — | — |
| `mito_frac` | mitochondrial UMIs / `n_UMIs` | `^MT-` | `^mt-` |
| `pct_counts_rb` | ribosomal UMIs / `n_UMIs` | `^RPS`, `^RPL` | same |

scanpy aliases: `n_genes_by_counts`, `total_counts`, `pct_counts_mt`, `pct_counts_rb`.

## Recommended Metrics

| Column | Definition | Notes |
|--------|------------|-------|
| `hbb_score` | hemoglobin gene UMIs / `n_UMIs` | blood contamination; use `assets/gene_sets/hbb_genes_*.txt` |
| `doublet_score` | continuous doublet score | DoubletFinder or Scrublet |
| `predicted_doublet` | logical/class doublet call | filter explicitly |
| `phase` | cell cycle phase | G1 / S / G2M |
| `s_score` | S phase score | Seurat CellCycleScoring |
| `g2m_score` | G2M phase score | Seurat CellCycleScoring |
| `is_HQ` | high-quality gate for advanced steps | common: `n_genes >= 500` |

## Extended Metrics

| Column | Definition | Requires |
|--------|------------|----------|
| `chrY_frac` | chrY gene UMIs / `n_UMIs` | `assets/gene_sets/chrY_genes_human.txt` or species list |
| `ambient_frac` | ambient RNA contamination | `celda::decontX` (R) or SoupX |
| `nuclear_frac` | intron / (intron + exon) | STARsolo Velocyto spliced/unspliced matrices |

## Derived / Workflow Columns

| Column | Definition |
|--------|------------|
| `pass_qc` | final retain flag after all filters |
| `fail_reason` | optional comma-separated reasons |
| `sample_id` | sample identifier |

## Tissue-Specific Optional Scores

Add when biologically relevant:

| Tissue | Suggested score | Genes |
|--------|-----------------|-------|
| Lung | epithelial ambient | `EPCAM`, `KRT8`, `KRT18` fraction |
| Lung | immune enrichment | `PTPRC` fraction |
| Liver | hepatocyte score | `ALB`, `APOA1` |
| PBMC/BM | RBC contamination | `hbb_score` (critical) |

Keep tissue scores in the project script or a project-local `utils/qc_gene_sets.yaml`, not in this skill, unless they become reusable.

## Metric Computation Order

```text
1. n_genes, n_UMIs
2. mito_frac, pct_counts_rb
3. hbb_score, chrY_frac (if species lists available)
4. is_HQ gate
5. ambient_frac (on HQ cells)
6. doublet_score (on HQ cells)
7. cell cycle (on HQ cells or all cells — document choice)
8. nuclear_frac (if STARsolo available)
```

## What Not To Treat As QC Metrics

- normalized expression, log1p values, scaled data;
- PCA/UMAP coordinates;
- cluster labels (downstream);
- batch-corrected embeddings.

Compute or attach these only after QC filtering unless needed for diagnosis.

# QC Project Layout

Aligned with `analysis-workflow`. Scripts are **staged tools**; the analyst advances stage by stage after reviewing outputs. See `human-in-the-loop.md`.

## Principle

```text
script -> result/<stage>/<script-id>/ -> figure/<stage>/<script-id>/
```

QC is checkpoint-heavy: metric tables and filtered matrices are formal `result/` artifacts, not `tmp/`.

## Recommended Stage: `01-qc`

```text
scripts/01-qc/
  01-calculate_metrics.R|py    # metrics only, no hard filtering
  02-qc_diagnosis.R|py         # figures from metadata
  03-filter_cells.R|py           # apply thresholds + optional doublet removal
  04-merge_qc_passed.R|py      # optional; only after per-sample QC
utils/
  qc_io.R|py                   # optional shared loaders
  qc_thresholds.yaml           # optional explicit threshold table
result/01-qc/
  01-calculate_metrics/<sample_id>/
    metadata.tsv.gz
    counts.mtx.gz               # optional checkpoint
  03-filter_cells/<sample_id>/
    filtered_counts.mtx.gz
    filter_summary.json
figure/01-qc/
  02-qc_diagnosis/<sample_id>/
    qc_violin.pdf
    qc_scatter_mito_genes.pdf
```

## Per-Sample vs Project-Level

| Artifact | Granularity |
|----------|-------------|
| `metadata.tsv` | per sample |
| diagnosis figures | per sample + optional cohort overview |
| filtered matrix | per sample |
| `QC_summary.tsv` | project-level |
| merged h5ad/Seurat | project-level, after QC |

## Script Header Template

```python
# Script: 03-filter_cells.py
# Purpose: per-sample cell filtering from QC metadata
# Input:  result/01-qc/01-calculate_metrics/<sample_id>/metadata.tsv
# Output: result/01-qc/03-filter_cells/<sample_id>/
# Figure: figure/01-qc/03-filter_cells/<sample_id>/
# Status: draft | validated | production
```

## Git Boundary

Do not commit large matrices. Commit:

- scripts, utils, threshold YAML;
- small `QC_summary.tsv` if intentional;
- `workflow_map.md` / `data_lineage.md` entries.

## Workflow Map Entry Example

```markdown
| Step | Script | Result | Figure |
|------|--------|--------|--------|
| QC metrics | scripts/01-qc/01-calculate_metrics.R | result/01-qc/01-calculate_metrics/ | figure/01-qc/01-calculate_metrics/ |
| QC filter | scripts/01-qc/03-filter_cells.R | result/01-qc/03-filter_cells/ | figure/01-qc/03-filter_cells/ |
```

# Python / scanpy QC

## Environment

Prefer project-local pixi/conda env with: `scanpy`, `anndata`, `numpy`, `pandas`, `matplotlib`, `scrublet`.

Use bundled script first:

```bash
python <skill-root>/scripts/calculate_metrics.py --help
```

## Load 10x Matrix

```python
import scanpy as sc

# prefer raw for QC metric computation
adata = sc.read_10x_mtx(
    "data/cellranger/SAMPLE/outs/raw_feature_bc_matrix",
    var_names="gene_symbols",
    make_unique=True,
)
```

From existing h5ad: load and ensure `.X` is raw counts (not log-normalized).

## Core Metrics

```python
species = "human"  # or "mouse"
mt_prefix = "MT-" if species == "human" else "mt-"

adata.var["mt"] = adata.var_names.str.startswith(mt_prefix)
adata.var["rb"] = adata.var_names.str.startswith(("RPS", "RPL"))

sc.pp.calculate_qc_metrics(
    adata,
    qc_vars=["mt", "rb"],
    percent_top=None,
    log1p=False,
    inplace=True,
)

# map to canonical names when saving
adata.obs["n_genes"] = adata.obs["n_genes_by_counts"]
adata.obs["n_UMIs"] = adata.obs["total_counts"]
adata.obs["mito_frac"] = adata.obs["pct_counts_mt"] / 100.0
```

## Extended Metrics (hemoglobin)

```python
from pathlib import Path

def read_gene_list(path):
    return [g.strip() for g in Path(path).read_text().splitlines() if g.strip()]

skill_root = Path("<skill-root>")
hbb = read_gene_list(skill_root / "assets/gene_sets/hbb_genes_human.txt")
hbb = [g for g in hbb if g in adata.var_names]
adata.obs["hbb_score"] = adata[:, hbb].X.sum(axis=1).A1 / adata.obs["n_UMIs"]
```

## HQ Gate + Scrublet (per sample)

```python
import scrublet as scr

hq = adata.obs["n_genes"] >= 500
counts = adata[hq].X

scrub = scr.Scrublet(counts, expected_doublet_rate=0.06)
scores, preds = scrub.scrub_doublets(min_counts=2, min_cells=3, n_prin_comps=30)

adata.obs["doublet_score"] = float("nan")
adata.obs["predicted_doublet"] = False
adata.obs.loc[hq, "doublet_score"] = scores
adata.obs.loc[hq, "predicted_doublet"] = preds
```

## Filter Cells

```python
# example fixed thresholds — prefer per-sample MAD in production
keep = (
    (adata.obs["n_genes"] >= 200)
    & (adata.obs["n_genes"] <= 6000)
    & (adata.obs["n_UMIs"] >= 400)
    & (adata.obs["mito_frac"] < 0.20)
    & (~adata.obs["predicted_doublet"])
)
adata = adata[keep].copy()
sc.pp.filter_genes(adata, min_cells=10)
```

## Save Checkpoint

```python
out = "result/01-qc/03-filter_cells/SAMPLE"
adata.write_h5ad(f"{out}/filtered.h5ad", compression="gzip")
adata.obs.to_csv(f"{out}/metadata.tsv", sep="\t")
```

## Figures (publication defaults)

```python
import matplotlib as mpl
mpl.rcParams.update({
    "font.family": "Arial",
    "pdf.fonttype": 42,
    "ps.fonttype": 42,
})

import matplotlib.pyplot as plt
import seaborn as sns

fig, ax = plt.subplots(figsize=(4, 4))
ax.scatter(
    adata.obs["mito_frac"],
    adata.obs["n_genes"],
    s=1,
    alpha=0.3,
    rasterized=True,
)
fig.savefig("figure/01-qc/02-qc_diagnosis/SAMPLE/qc_scatter.pdf", dpi=300)
```

## Anti-Pattern: Merge-First QC

Do not copy this pattern for new work:

```python
merged = adatas[0].concatenate(adatas[1:])
sc.pp.filter_cells(merged, min_genes=180)  # global on merged
fc.scrublet_doublet_removal_10X(merged)    # global doublet
```

Split into per-sample metrics → filter → merge.

## Optional: SoupX (ambient RNA)

Requires raw + filtered matrices and cluster labels or auto-estimation. Implement in project `utils/` when user requests ambient correction; not in bundled minimal script.

# R / Seurat QC

## Environment

Typical packages: `Matrix`, `Seurat`, `ggplot2`, `celda` (decontX), `SingleCellExperiment`.

Use bundled script first:

```bash
Rscript <skill-root>/scripts/calculate_metrics.R --help
```

## Load 10x Raw Matrix

```r
library(Matrix)

read_10x_counts <- function(matrix_dir) {
  counts <- Matrix::readMM(file.path(matrix_dir, "matrix.mtx.gz"))
  counts <- as(counts, "CsparseMatrix")
  features <- read.table(file.path(matrix_dir, "features.tsv.gz"), header = FALSE)
  barcodes <- read.table(file.path(matrix_dir, "barcodes.tsv.gz"), header = FALSE)
  rownames(counts) <- make.unique(features$V2)
  colnames(counts) <- barcodes$V1
  counts
}
```

## Core Metrics

```r
metadata <- data.frame(
  row.names = colnames(counts),
  n_genes = Matrix::colSums(counts > 0),
  n_UMIs = Matrix::colSums(counts),
  stringsAsFactors = FALSE
)

mito_genes <- grep("^MT-", rownames(counts), value = TRUE)  # mouse: ^mt-
metadata$mito_frac <- Matrix::colSums(counts[mito_genes, , drop = FALSE]) / metadata$n_UMIs

rb_genes <- grep("^(RPS|RPL)", rownames(counts), value = TRUE)
metadata$pct_counts_rb <- 100 * Matrix::colSums(counts[rb_genes, , drop = FALSE]) / metadata$n_UMIs
```

## Extended Metrics

Reference implementation: `<project-root>/scripts/calculate_metrics_extended.R`

| Step | Metric | Function |
|------|--------|----------|
| chrY | `chrY_frac` | gene list intersection |
| nuclear | `nuclear_frac` | STARsolo Velocyto spliced/unspliced |
| ambient | `ambient_frac_decontX` | `celda::decontX` on HQ cells |
| doublet | `doublet_score` | DoubletFinder on HQ cells |
| cycle | `phase`, `s_score`, `g2m_score` | `Seurat::CellCycleScoring` |
| blood | `hbb_score` | hemoglobin gene set |

```r
skill_root <- "<skill-root>"
hbb_genes <- intersect(
  readLines(file.path(skill_root, "assets/gene_sets/hbb_genes_human.txt")),
  rownames(counts)
)
metadata$hbb_score <- Matrix::colSums(counts[hbb_genes, , drop = FALSE]) / metadata$n_UMIs
```

## HQ Gate

```r
metadata$is_HQ <- metadata$n_genes >= 500
```

Use HQ cells for decontX and DoubletFinder to stabilize estimation.

## decontX (ambient RNA)

```r
library(celda)
library(SingleCellExperiment)

sce <- SingleCellExperiment(
  assays = list(counts = counts[, metadata$is_HQ]),
  colData = metadata[metadata$is_HQ, , drop = FALSE]
)
sce <- decontX(sce)
metadata$ambient_frac <- NA_real_
metadata[rownames(colData(sce)), "ambient_frac"] <- colData(sce)$decontX_contamination
```

## DoubletFinder

Project-local `R/doubletFinder.R` may exist (GZL course). Pattern:

```r
seu <- CreateSeuratObject(counts = counts[, metadata$is_HQ], meta.data = metadata[metadata$is_HQ, ])
seu <- FindVariableFeatures(seu, selection.method = "vst", nfeatures = 2000)
# source project doubletFinder.R
# df <- doubletFinder(counts = counts[, metadata$is_HQ], select.genes = VariableFeatures(seu))
```

Prefer `scDblFinder` for new projects if DoubletFinder is not already available.

## Cell Cycle

```r
seu <- NormalizeData(seu)
seu <- CellCycleScoring(
  seu,
  s.features = intersect(cc.genes.updated.2019$s.genes, rownames(counts)),
  g2m.features = intersect(cc.genes.updated.2019$g2m.genes, rownames(counts))
)
metadata$phase <- NA_character_
metadata[rownames(seu@meta.data), "phase"] <- seu@meta.data$Phase
```

## Filter Cells

```r
keep <- with(metadata,
  n_genes >= 200 &
  n_genes <= 6000 &
  n_UMIs >= 400 &
  mito_frac < 0.20 &
  (is.na(ambient_frac) | ambient_frac < 0.25) &
  (is.na(doublet_score) | doublet_score < 0.25)  # project-specific
)
counts_filt <- counts[, keep]
```

## Save Checkpoint

```r
out_dir <- "result/01-qc/01-calculate_metrics/SAMPLE"
dir.create(out_dir, recursive = TRUE, showWarnings = FALSE)
write.table(metadata, file.path(out_dir, "metadata.tsv"),
  sep = "\t", quote = FALSE, row.names = TRUE, col.names = TRUE)
```

## Figures

```r
plot_font_family <- "Arial"
theme_set(theme_linedraw(base_family = plot_font_family))

p <- ggplot(metadata, aes(x = mito_frac, y = n_genes)) +
  geom_point(size = 0.1, alpha = 0.2, raster = TRUE) +
  scale_y_log10()

ggsave("figure/01-qc/02-qc_diagnosis/SAMPLE/qc_scatter.pdf",
  p, width = 4, height = 4, device = cairo_pdf)
```

## Seurat Object Path

If the project is Seurat-native, save filtered object:

```r
seu_filt <- CreateSeuratObject(counts = counts_filt, meta.data = metadata[keep, , drop = FALSE])
qs::qsave(seu_filt, "result/01-qc/03-filter_cells/SAMPLE/filtered.seurat.qs")
```

#!/usr/bin/env python3
"""Per-sample scRNA QC metric calculation (scanpy).

Human-in-the-loop tool: computes metrics and writes metadata.tsv.
Does NOT apply filtering thresholds. Run filter scripts only after
the analyst reviews diagnosis outputs and confirms cutoffs.

Bundled with singlecell-qc skill.

Example:
  python calculate_metrics.py \\
    --matrix-dir data/cellranger/SAMPLE/outs/raw_feature_bc_matrix \\
    --sample-id SAMPLE \\
    --species human \\
    --output-dir result/01-qc/01-calculate_metrics/SAMPLE
"""

from __future__ import annotations

import argparse
import gzip
import json
from pathlib import Path

import numpy as np
import pandas as pd
import scanpy as sc
from scipy import sparse


SKILL_ROOT = Path(__file__).resolve().parents[1]
GENE_SETS = SKILL_ROOT / "assets" / "gene_sets"


def read_gene_list(path: Path) -> list[str]:
    genes = []
    for line in path.read_text().splitlines():
        g = line.strip()
        if g and not g.startswith("#"):
            genes.append(g)
    return genes


def fraction_score(adata, genes: list[str], denom: np.ndarray) -> np.ndarray:
    if not genes:
        return np.zeros(adata.n_obs, dtype=float)
    sub = adata[:, genes].X
    if sparse.issparse(sub):
        num = np.asarray(sub.sum(axis=1)).ravel()
    else:
        num = sub.sum(axis=1)
    with np.errstate(divide="ignore", invalid="ignore"):
        out = num / denom
    out[~np.isfinite(out)] = 0.0
    return out


def load_10x(matrix_dir: Path, min_genes: int = 50) -> sc.AnnData:
    adata = sc.read_10x_mtx(
        str(matrix_dir),
        var_names="gene_symbols",
        make_unique=True,
    )
    if min_genes > 0:
        sc.pp.filter_cells(adata, min_genes=min_genes)
    return adata


def compute_metrics(
    adata: sc.AnnData,
    species: str,
    hq_n_genes: int,
    run_scrublet: bool,
    expected_doublet_rate: float,
) -> pd.DataFrame:
    mt_prefix = "MT-" if species == "human" else "mt-"
    adata.var["mt"] = adata.var_names.str.startswith(mt_prefix)
    adata.var["rb"] = adata.var_names.str.startswith(("RPS", "RPL"))

    sc.pp.calculate_qc_metrics(
        adata,
        qc_vars=["mt", "rb"],
        percent_top=None,
        log1p=False,
        inplace=True,
    )

    obs = adata.obs.copy()
    obs["n_genes"] = obs["n_genes_by_counts"].astype(int)
    obs["n_UMIs"] = obs["total_counts"].astype(float)
    obs["mito_frac"] = obs["pct_counts_mt"].astype(float) / 100.0
    obs["pct_counts_rb"] = obs["pct_counts_rb"].astype(float)
    obs["is_HQ"] = obs["n_genes"] >= hq_n_genes

    hbb_file = GENE_SETS / f"hbb_genes_{species}.txt"
    if hbb_file.exists():
        hbb = [g for g in read_gene_list(hbb_file) if g in adata.var_names]
        obs["hbb_score"] = fraction_score(adata, hbb, obs["n_UMIs"].to_numpy())

    if species == "human":
        chry_file = GENE_SETS / "chrY_genes_human.txt"
        if chry_file.exists():
            chry = [g for g in read_gene_list(chry_file) if g in adata.var_names]
            obs["chrY_frac"] = fraction_score(adata, chry, obs["n_UMIs"].to_numpy())

    obs["doublet_score"] = np.nan
    obs["predicted_doublet"] = False

    if run_scrublet and obs["is_HQ"].sum() >= 100:
        import scrublet as scr

        hq_idx = obs.index[obs["is_HQ"]]
        hq = adata[hq_idx]
        x = hq.X.toarray() if hasattr(hq.X, "toarray") else np.asarray(hq.X)
        scrub = scr.Scrublet(x, expected_doublet_rate=expected_doublet_rate)
        scores, preds = scrub.scrub_doublets(
            min_counts=2,
            min_cells=3,
            min_gene_variability_pctl=85,
            n_prin_comps=min(30, hq.n_obs - 1, hq.n_vars - 1),
        )
        obs.loc[hq_idx, "doublet_score"] = scores
        obs.loc[hq_idx, "predicted_doublet"] = preds

    return obs


def write_metadata(obs: pd.DataFrame, out_dir: Path, sample_id: str) -> None:
    out_dir.mkdir(parents=True, exist_ok=True)
    obs = obs.copy()
    obs.insert(0, "barcode", obs.index.astype(str))
    obs.insert(0, "sample_id", sample_id)

    tsv = out_dir / "metadata.tsv"
    obs.to_csv(tsv, sep="\t")

    with gzip.open(out_dir / "metadata.tsv.gz", "wt") as fh:
        obs.to_csv(fh, sep="\t")

    summary = {
        "sample_id": sample_id,
        "n_cells": int(obs.shape[0]),
        "median_n_genes": float(np.median(obs["n_genes"])),
        "median_n_UMIs": float(np.median(obs["n_UMIs"])),
        "median_mito_frac": float(np.median(obs["mito_frac"])),
        "n_HQ": int(obs["is_HQ"].sum()),
        "n_predicted_doublet": int(obs["predicted_doublet"].sum()),
    }
    (out_dir / "metrics_summary.json").write_text(json.dumps(summary, indent=2))


def parse_args() -> argparse.Namespace:
    p = argparse.ArgumentParser(description="Per-sample scRNA QC metrics (scanpy)")
    p.add_argument("--matrix-dir", required=True, type=Path, help="10x matrix directory")
    p.add_argument("--sample-id", required=True, help="Sample identifier")
    p.add_argument("--species", choices=["human", "mouse"], default="human")
    p.add_argument("--output-dir", required=True, type=Path)
    p.add_argument("--min-genes", type=int, default=50, help="Pre-filter min genes")
    p.add_argument("--hq-n-genes", type=int, default=500, help="HQ gate for Scrublet")
    p.add_argument("--run-scrublet", action="store_true", help="Run Scrublet on HQ cells")
    p.add_argument("--expected-doublet-rate", type=float, default=0.06)
    p.add_argument("--save-h5ad", action="store_true", help="Save AnnData with obs metrics")
    return p.parse_args()


def main() -> None:
    args = parse_args()
    adata = load_10x(args.matrix_dir, min_genes=args.min_genes)
    obs = compute_metrics(
        adata,
        species=args.species,
        hq_n_genes=args.hq_n_genes,
        run_scrublet=args.run_scrublet,
        expected_doublet_rate=args.expected_doublet_rate,
    )
    adata.obs = obs
    write_metadata(obs, args.output_dir, args.sample_id)
    if args.save_h5ad:
        adata.write_h5ad(args.output_dir / "with_metrics.h5ad", compression="gzip")
    print(f"Wrote metrics for {args.sample_id}: {obs.shape[0]} cells -> {args.output_dir}")


if __name__ == "__main__":
    main()

#!/usr/bin/env python3
"""Summarize per-sample QC metadata.tsv files."""

from __future__ import annotations

import argparse
import json
from pathlib import Path

import numpy as np
import pandas as pd


CORE_COLS = ["n_genes", "n_UMIs", "mito_frac", "pct_counts_rb", "hbb_score", "is_HQ", "predicted_doublet"]


def load_metadata(path: Path) -> pd.DataFrame:
    if path.suffix == ".gz":
        df = pd.read_csv(path, sep="\t", index_col=0)
    else:
        df = pd.read_csv(path, sep="\t", index_col=0)
    return df


def summarize_sample(df: pd.DataFrame, sample_id: str) -> dict:
    out = {"sample_id": sample_id, "n_cells": int(df.shape[0])}
    for col in CORE_COLS:
        if col in df.columns:
            if df[col].dtype == bool or col in {"is_HQ", "predicted_doublet"}:
                out[f"n_{col}"] = int(df[col].sum())
            else:
                out[f"median_{col}"] = float(np.median(df[col]))
                out[f"p95_{col}"] = float(np.percentile(df[col], 95))
    return out


def find_metadata_files(root: Path) -> list[Path]:
    files = list(root.rglob("metadata.tsv")) + list(root.rglob("metadata.tsv.gz"))
    # prefer deepest sample-level paths
    return sorted(set(files))


def parse_args() -> argparse.Namespace:
    p = argparse.ArgumentParser(description="Summarize QC metadata.tsv files")
    p.add_argument("--input-dir", required=True, type=Path, help="Root with per-sample metadata")
    p.add_argument("--output", type=Path, default=None, help="QC_summary.tsv path")
    return p.parse_args()


def main() -> None:
    args = parse_args()
    files = find_metadata_files(args.input_dir)
    if not files:
        raise SystemExit(f"No metadata.tsv found under {args.input_dir}")

    rows = []
    for f in files:
        df = load_metadata(f)
        sample_id = df["sample_id"].iloc[0] if "sample_id" in df.columns else f.parent.name
        rows.append(summarize_sample(df, str(sample_id)))

    summary = pd.DataFrame(rows).sort_values("sample_id")
    out = args.output or (args.input_dir / "QC_summary.tsv")
    summary.to_csv(out, sep="\t", index=False)
    print(summary.to_string(index=False))
    print(f"\nWrote {out}")


if __name__ == "__main__":
    main()