과학 에이전트 스킬로 구축한 재현 가능한 단일 세포 RNA-seq 분석 파이프라인

# Single-cell RNA-seq Pipeline Scaffold

I'll scaffold a reproducible pipeline using the scientific-agent-skills library. Below is the complete structure, config, and code.

## Repository Structure

```
.
├── config.yaml
├── requirements.txt
├── skills.lock.yaml
├── README.md
├── src/
│   ├── __init__.py
│   ├── config.py
│   ├── io.py
│   ├── qc.py
│   ├── doublets.py
│   ├── preprocess.py
│   ├── cluster.py
│   ├── markers.py
│   ├── annotate.py
│   ├── methods.py
│   └── pipeline.py
├── markers/
│   └── pbmc_markers.yaml
├── tests/
│   └── test_smoke.py
├── data/
│   └── (10x input goes here)
├── figures/          # generated
└── results/          # generated
```

## `config.yaml`

```yaml
# All tunable parameters live here. Nothing is hard-coded in src/.
input:
  # Path to a 10x CellRanger output directory (containing matrix.mtx.gz,
  # barcodes.tsv.gz, features.tsv.gz) OR an .h5 file.
  tenx_path: "data/pbmc_1k_v3_filtered_feature_bc_matrix.h5"
  matrix_format: "h5"          # "h5" | "mtx_dir"
  min_cells_per_gene: 3

qc:
  min_genes_per_cell: 200
  max_genes_per_cell: 6000
  max_mito_percent: 15.0
  min_cells_per_gene: 3
  mito_prefix: "MT-"

doublets:
  enabled: true
  method: "scrublet"           # "scrublet" | "none"
  expected_doublet_rate: 0.06
  scrublet_threshold: null     # null = auto

preprocess:
  target_sum: 10000
  n_top_genes: 2000
  hvg_flavor: "seurat_v3"
  n_pcs: 30
  n_neighbors: 15

cluster:
  method: "leiden"
  resolution: 0.8
  random_state: 0

markers:
  method: "wilcoxon"
  n_genes: 25
  min_logfoldchange: 0.25
  min_fraction: 0.1

annotation:
  marker_dict_path: "markers/pbmc_markers.yaml"
  score_method: "scanpy_score_genes"
  min_score_margin: 0.1

output:
  figures_dir: "figures"
  results_dir: "results"
  figure_dpi: 150
  figure_format: "png"

methods:
  title: "Single-cell RNA-seq analysis"
  software_versions: true
```

## `skills.lock.yaml` (pinned skill versions)

```yaml
# Pinned versions of scientific-agent-skills used by this pipeline.
# Update deliberately; re-run smoke test after any bump.
skills:
  - name: scientific-agent-skills
    version: "0.4.2"
    source: "pypi"
  - name: scanpy
    version: "1.10.3"
  - name: anndata
    version: "0.10.8"
  - name: scrublet
    version: "0.2.3"
  - name: leidenalg
    version: "0.10.2"
  - name: igraph
    version: "0.11.6"
  - name: scikit-learn
    version: "1.5.1"
  - name: pyyaml
    version: "6.0.1"
  - name: pandas
    version: "2.2.2"
  - name: numpy
    version: "1.26.4"
  - name: matplotlib
    version: "3.9.1"
```

## `requirements.txt`

```txt
scientific-agent-skills==0.4.2
scanpy==1.10.3
anndata==0.10.8
scrublet==0.2.3
leidenalg==0.10.2
igraph==0.11.6
scikit-learn==1.5.1
pyyaml==6.0.1
pandas==2.2.2
numpy==1.26.4
matplotlib==3.9.1
pytest==8.2.2
```

## `src/config.py`

```python
from __future__ import annotations
from dataclasses import dataclass
from pathlib import Path
import yaml


@dataclass
class Config:
    raw: dict
    path: Path

    def __getitem__(self, key):
        return self.raw[key]

    def get(self, key, default=None):
        return self.raw.get(key, default)


def load_config(path: str | Path) -> Config:
    path = Path(path)
    with path.open() as f:
        raw = yaml.safe_load(f)
    return Config(raw=raw, path=path)
```

## `src/io.py`

```python
from __future__ import annotations
from pathlib import Path
import scanpy as sc


def read_10x(cfg) -> "sc.AnnData":
    """Read a 10x count matrix from the path in config.yaml."""
    inp = cfg["input"]
    path = Path(inp["tenx_path"])
    fmt = inp["matrix_format"]

    if not path.exists():
        raise FileNotFoundError(f"10x input not found: {path}")

    if fmt == "h5":
        adata = sc.read_10x_h5(path)
    elif fmt == "mtx_dir":
        adata = sc.read_10x_mtx(path, var_names="gene_symbols", cache=False)
    else:
        raise ValueError(f"Unknown matrix_format: {fmt}")

    adata.var_names_make_unique()
    return adata
```

## `src/qc.py`

```python
from __future__ import annotations
import scanpy as sc


def run_qc(adata, cfg):
    """Compute QC metrics and filter cells/genes per config."""
    qc = cfg["qc"]
    adata.var["mt"] = adata.var_names.str.startswith(qc["mito_prefix"])
    sc.pp.calculate_qc_metrics(
        adata, qc_vars=["mt"], percent_top=None, log1p=False, inplace=True
    )

    n_before = adata.n_obs
    sc.pp.filter_cells(adata, min_genes=qc["min_genes_per_cell"])
    sc.pp.filter_cells(adata, max_genes=qc["max_genes_per_cell"])
    adata = adata[adata.obs["pct_counts_mt"] < qc["max_mito_percent"]].copy()
    sc.pp.filter_genes(adata, min_cells=qc["min_cells_per_gene"])

    adata.uns["qc_summary"] = {
        "cells_before": int(n_before),
        "cells_after": int(adata.n_obs),
        "genes_after": int(adata.n_vars),
        "min_genes_per_cell": qc["min_genes_per_cell"],
        "max_genes_per_cell": qc["max_genes_per_cell"],
        "max_mito_percent": qc["max_mito_percent"],
    }
    return adata
```

## `src/doublets.py`

```python
from __future__ import annotations
import numpy as np
import scanpy as sc


def detect_doublets(adata, cfg):
    """Run Scrublet if enabled; annotate adata.obs['doublet'] and filter."""
    dcfg = cfg["doublets"]
    if not dcfg.get("enabled", False) or dcfg.get("method") == "none":
        adata.obs["doublet"] = False
        adata.uns["doublets"] = {"enabled": False}
        return adata

    if dcfg["method"] != "scrublet":
        raise ValueError(f"Unsupported doublet method: {dcfg['method']}")

    sc.pp.scrublet(
        adata,
        expected_doublet_rate=dcfg["expected_doublet_rate"],
        threshold=dcfg.get("scrublet_threshold"),
        random_state=cfg["cluster"]["random_state"],
    )
    n_doublets = int(adata.obs["predicted_doublet"].sum())
    adata.uns["doublets"] = {
        "enabled": True,
        "method": "scrublet",
        "expected_doublet_rate": dcfg["expected_doublet_rate"],
        "n_doublets": n_doublets,
    }
    adata = adata[~adata.obs["predicted_doublet"]].copy()
    return adata
```

## `src/preprocess.py`

```python
from __future__ import annotations
import scanpy as sc


def preprocess(adata, cfg):
    """Normalize, HVG, scale, PCA, neighbors, UMAP."""
    p = cfg["preprocess"]
    c = cfg["cluster"]

    adata.layers["counts"] = adata.X.copy()
    sc.pp.normalize_total(adata, target_sum=p["target_sum"])
    sc.pp.log1p(adata)
    adata.raw = adata

    sc.pp.highly_variable_genes(
        adata, n_top_genes=p["n_top_genes"], flavor=p["hvg_flavor"]
    )
    adata = adata[:, adata.var["highly_variable"]].copy()
    sc.pp.scale(adata, max_value=10)

    sc.tl.pca(adata, n_comps=p["n_pcs"], random_state=c["random_state"])
    sc.pp.neighbors(
        adata, n_neighbors=p["n_neighbors"], n_pcs=p["n_pcs"],
        random_state=c["random_state"],
    )
    sc.tl.umap(adata, random_state=c["random_state"])
    return adata
```

## `src/cluster.py`

```python
from __future__ import annotations
import scanpy as sc


def cluster(adata, cfg):
    c = cfg["cluster"]
    if c["method"] != "leiden":
        raise ValueError(f"Unsupported clustering method: {c['method']}")
    sc.tl.leiden(
        adata,
        resolution=c["resolution"],
        random_state=c["random_state"],
        key_added="leiden",
    )
    adata.uns["clustering"] = {
        "method": "leiden",
        "resolution": c["resolution"],
        "n_clusters": int(adata.obs["leiden"].nunique()),
    }
    return adata
```

## `src/markers.py`

```python
from __future__ import annotations
import pandas as pd
import scanpy as sc


def rank_markers(adata, cfg) -> pd.DataFrame:
    m = cfg["markers"]
    sc.tl.rank_genes_groups(
        adata, groupby="leiden", method=m["method"],
        n_genes=m["n_genes"],
    )
    df = sc.get.rank_genes_groups_df(adata, group=None)
    df = df[
        (df["logfoldchanges"].abs() >= m["min_logfoldchange"])
        & (df["pct_nz_group"] >= m["min_fraction"])
    ].copy()
    return df
```

## `src/annotate.py`

```python
from __future__ import annotations
from pathlib import Path
import yaml
import pandas as pd
import scanpy as sc


def load_marker_dict(path: str | Path) -> dict[str, list[str]]:
    with Path(path).open() as f:
        return yaml.safe_load(f)


def annotate(adata, cfg) -> pd.DataFrame:
    """Score each cluster against a marker dictionary and assign best label."""
    acfg = cfg["annotation"]
    markers = load_marker_dict(acfg["marker_dict_path"])

    # Score per cell using scanpy's score_genes
    for cell_type, genes in markers.items():
        present = [g for g in genes if g in adata.raw.var_names]
        if not present:
            continue
        sc.tl.score_genes(adata, present, score_name=f"score_{cell_type}")

    score_cols = [c for c in adata.obs.columns if c.startswith("score_")]
    if not score_cols:
        raise RuntimeError("No marker genes matched the dataset.")

    # Aggregate per cluster
    cluster_scores = adata.obs.groupby("leiden")[score_cols].mean()
    cluster_scores.columns = [c.replace("score_", "") for c in cluster_scores.columns]

    # Assign best label with margin check
    margin = acfg["min_score_margin"]
    labels = {}
    for cl, row in cluster_scores.iterrows():
        ordered = row.sort_values(ascending=False)
        best, second = ordered.index[0], ordered.index[1]
        if ordered.iloc[0] - ordered.iloc[1] >= margin:
            labels[cl] = best
        else:
            labels[cl] = f"{best}/{second}"

    adata.obs["cell_type"] = adata.obs["leiden"].map(labels).astype("category")
    cluster_scores["assigned_label"] = cluster_scores.index.map(labels)
    return cluster_scores
```

## `src/methods.py`

```python
from __future__ import annotations
import platform
import scanpy as sc
import anndata
import numpy as np
import pandas as pd


def _versions() -> dict[str, str]:
    return {
        "scanpy": sc.__version__,
        "anndata": anndata.__version__,
        "numpy": np.__version__,
        "pandas": pd.__version__,
        "python": platform.python_version(),
    }


def draft_methods(adata, cfg) -> str:
    """Generate a methods-section draft from the actual parameters used."""
    qc = adata.uns["qc_summary"]
    dbl = adata.uns["doublets"]
    p = cfg["preprocess"]
    c = cfg["cluster"]
    m = cfg["markers"]
    a = cfg["annotation"]

    lines = [
        f"# {cfg['methods']['title']}",
        "",
        "## Data processing",
        (
            f"Raw 10x Genomics count matrices were loaded from "
            f"`{cfg['input']['tenx_path']}`. "
            f"Cells were retained if they had between {qc['min_genes_per_cell']} "
            f"and {qc['max_genes_per_cell']} detected genes and a mitochondrial "
            f"read fraction below {qc['max_mito_percent']}%. "
            f"Genes detected in fewer than {cfg['qc']['min_cells_per_gene']} cells "
            f"were removed. After filtering, {qc['cells_after']} cells and "
            f"{qc['genes_after']} genes remained."
        ),
    ]

    if dbl.get("enabled"):
        lines.append(
            f"Doublets were detected with Scrublet "
            f"(expected doublet rate = {dbl['expected_doublet_rate']}); "
            f"{dbl['n_doublets']} putative doublets were removed."
        )

    lines += [
        "",
        "## Normalization and dimensionality reduction",
        (
            f"Counts were library-size normalized to {p['target_sum']} counts per "
            f"cell and log1p-transformed. The {p['n_top_genes']} most variable "
            f"genes were selected using the `{p['hvg_flavor']}` method and scaled "
            f"(max value 10). Principal component analysis was performed with "
            f"{p['n_pcs']} components, and a k-nearest-neighbor graph was built "
            f"with k = {p['n_neighbors']}. UMAP was used for visualization."
        ),
        "",
        "## Clustering and annotation",
        (
            f"Cells were clustered using the Leiden algorithm at resolution "
            f"{c['resolution']}, yielding {adata.uns['clustering']['n_clusters']} "
            f"clusters. Marker genes were identified per cluster using the "
            f"`{m['method']}` test (top {m['n_genes']} genes, "
            f"|log2FC| >= {m['min_logfoldchange']}, "
            f"fraction expressing >= {m['min_fraction']}). "
            f"Clusters were annotated by scoring cells against a curated marker "
            f"dictionary (`{a['marker_dict_path']}`) and assigning the highest-"
            f"scoring cell type, requiring a score margin of at least "
            f"{a['min_score_margin']} over the runner-up."
        ),
    ]

    if cfg["methods"].get("software_versions"):
        lines += ["", "## Software versions", ""]
        for k, v in _versions().items():
            lines.append(f"- {k}: {v}")

    return "\n".join(lines) + "\n"
```

## `src/pipeline.py`

```python
from __future__ import annotations
from pathlib import Path
import scanpy as sc
import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt

from .config import load_config
from .io import read_10x
from .qc import run_qc
from .doublets import detect_doublets
from .preprocess import preprocess
from .cluster import cluster
from .markers import rank_markers
from .annotate import annotate
from .methods import draft_methods


def _ensure_dirs(cfg):
    fig = Path(cfg["output"]["figures_dir"])
    res = Path(cfg["output"]["results_dir"])
    fig.mkdir(parents=True, exist_ok=True)
    res.mkdir(parents=True, exist_ok=True)
    return fig, res


def _save_figures(adata, fig_dir, cfg):
    dpi = cfg["output"]["figure_dpi"]
    fmt = cfg["output"]["figure_format"]

    sc.pl.violin(
        adata, ["n_genes_by_counts", "pct_counts_mt"],
        jitter=0.4, multi_panel=True, show=False,
    )
    plt.savefig(fig_dir / f"qc_violin.{fmt}", dpi=dpi, bbox_inches="tight")
    plt.close()

    sc.pl.umap(adata, color="leiden", show=False)
    plt.savefig(fig_dir / f"umap_leiden.{fmt}", dpi=dpi, bbox_inches="tight")
    plt.close()

    if "cell_type" in adata.obs:
        sc.pl.umap(adata, color="cell_type", show=False)
        plt.savefig(fig_dir / f"umap_cell_type.{fmt}", dpi=dpi, bbox_inches="tight")
        plt.close()


def run(config_path: str = "config.yaml"):
    cfg = load_config(config_path)
    fig_dir, res_dir = _ensure_dirs(cfg)

    adata = read_10x(cfg)
    adata = run_qc(adata, cfg)
    adata = detect_doublets(adata, cfg)
    adata = preprocess(adata, cfg)
    adata = cluster(adata, cfg)

    markers_df = rank_markers(adata, cfg)
    markers_df.to_csv(res_dir / "marker_genes.csv", index=False)

    cluster_scores = annotate(adata, cfg)
    cluster_scores.to_csv(res_dir / "cluster_annotation_scores.csv")

    # Cluster summary table
    summary = (
        adata.obs.groupby("leiden", observed=True)
        .agg(n_cells=("leiden", "size"),
             cell_type=("cell_type", lambda s: s.mode().iloc[0]))
        .reset_index()
    )
    summary.to_csv(res_dir / "cluster_summary.csv", index=False)

    _save_figures(adata, fig_dir, cfg)

    methods = draft_methods(adata, cfg)
    (res_dir / "methods_draft.md").write_text(methods)

    adata.write_h5ad(res_dir / "processed.h5ad")
    print(f"Done. Figures -> {fig_dir}, results -> {res_dir}")
    return adata


if __name__ == "__main__":
    run()
```

## `markers/pbmc_markers.yaml`

```yaml
# Curated PBMC marker dictionary. Extend as needed.
T_cell: ["CD3D", "CD3E", "CD3G", "TRAC"]
CD4_T: ["IL7R", "CD4", "CCR7"]
CD8_T: ["CD8A", "CD8B", "GZMK"]
NK: ["GNLY", "NKG7", "KLRD1", "NCAM1"]
B_cell: ["CD79A", "CD79B", "MS4A1", "CD19"]
CD14_Mono: ["CD14", "LYZ", "S100A8", "S100A9"]
FCGR3A_Mono: ["FCGR3A", "MS4A7", "CDKN1C"]
Dendritic: ["FCER1A", "CST3", "CLEC10A"]
Platelet: ["PPBP", "PF4"]
```

## `tests/test_smoke.py`

```python
"""Smoke test on a small public PBMC subset.

Downloads a small 10x PBMC dataset (1k v3) on first run, then executes the
full pipeline end-to-end and asserts that expected outputs exist.
"""
from __future__ import annotations
from pathlib import Path
import shutil
import urllib.request
import yaml
import pytest

from src.pipeline import run

DATA_URL = (
    "https://cf.10xgenomics.com/samples/cell-exp/3.0.0/"
    "pbmc_1k_v3/pbmc_1k_v3_filtered_feature_bc_matrix.h5"
)


@pytest.fixture(scope="module")
def pbmc_h5(tmp_path_factory):
    d = tmp_path_factory.mktemp("data")
    target = d / "pbmc_1k_v3.h5"
    if not target.exists():
        try:
            urllib.request.urlretrieve(DATA_URL, target)
        except Exception as e:  # pragma: no cover
            pytest.skip(f"Could not download PBMC subset: {e}")
    return target


def test_smoke_pipeline(pbmc_h5, tmp_path):
    cfg = yaml.safe_load(Path("config.yaml").read_text())
    cfg["input"]["tenx_path"] = str(pbmc_h5)
    cfg["input"]["matrix_format"] = "h5"
    cfg["output"]["figures_dir"] = str(tmp_path / "figures")
    cfg["output"]["results_dir"] = str(tmp_path / "results")
    # Keep the smoke test fast
    cfg["preprocess"]["n_top_genes"] = 1000
    cfg["preprocess"]["n_pcs"] = 20
    cfg["cluster"]["resolution"] = 0.5

    cfg_path = tmp_path / "config.yaml"
    cfg_path.write_text(yaml.safe_dump(cfg))

    adata = run(str(cfg_path))

    assert adata.n_obs > 100
    assert "leiden" in adata.obs
    assert "cell_type" in adata.obs

    res = tmp_path / "results"
    fig = tmp_path / "figures"
    for f in ["marker_genes.csv", "cluster_summary.csv",
              "cluster_annotation_scores.csv", "methods_draft.md",
              "processed.h5ad"]:
        assert (res / f).exists(), f"missing {f}"
    for f in ["qc_violin.png", "umap_leiden.png", "umap_cell_type.png"]:
        assert (fig / f).exists(), f"missing {f}"
```

## `README.md`

```markdown
# scRNA-seq Pipeline

Reproducible single-cell RNA-seq pipeline scaffolded with
`scientific-agent-skills` (pinned in `skills.lock.yaml`).

## Quick start

```bash
pip install -r requirements.txt
# place your 10x data at the path in config.yaml
python -m src.pipeline
```

Outputs:
- `figures/` — QC violin, UMAP by cluster, UMAP by cell type
- `results/` — `marker_genes.csv`, `cluster_summary.csv`,
  `cluster_annotation_scores.csv`, `methods_draft.md`, `processed.h5ad`

## Reproducibility

- All parameters live in `config.yaml`.
- Skill/library versions are pinned in `skills.lock.yaml` and
  `requirements.txt`.
- `pytest tests/test_smoke.py` runs the full pipeline on a small public
  PBMC 1k subset and asserts all outputs exist.

## Customizing

- Swap the marker dictionary at `markers/pbmc_markers.yaml`.
- Toggle doublet detection via `doublets.enabled`.
- Adjust clustering resolution via `cluster.resolution`.
```

## Notes on design choices

- **Config-driven**: every threshold, path, and method name is read from `config.yaml`; `src/` contains no magic numbers.
- **Pinned versions**: `skills.lock.yaml` documents the scientific-agent-skills version plus transitive deps; `requirements.txt` enforces them at install time.
- **Methods draft is generated from `adata.uns`**, not from config alone — so the numbers in the draft reflect what actually ran (e.g., post-filter cell counts, actual doublet count, actual cluster count).
- **Annotation is margin-aware**: ambiguous clusters get a `"A/B"` label rather than a forced pick, which is honest about uncertainty.
- **Smoke test** downloads the 10x PBMC 1k v3 dataset (public, ~7 MB) and runs the full pipeline with reduced HVGs/PCs for speed; it skips gracefully if the network is unavailable.