.catalog_stamp
1782943770000
스킬 불러오는 중
스킬 불러오는 중
aipoch/open-science · GitHub
Probabilistic single-cell RNA-seq with scvi-tools — scVI for a batch-corrected latent space, scANVI for semi-supervised label transfer, and Bayesian differential expression. Reach for this skill to integrate scRNA-seq batches, embed cells for clustering, transfer annotations from a reference onto a query, or score differentially expressed genes per cluster. For spatial deconvolution / mapping use the cell2location, DestVI, or Tangram methods instead.
프로젝트 폴더에서 아래 명령어를 실행하고, 설치할 에이전트를 선택하세요.
npx skills add aipoch/open-science --skill scvi-tools설치 명령을 직접 실행해야 적용됩니다. 지원 에이전트와 필요한 권한·라이선스는 제작자의 안내를 확인하세요.
.catalog_stamp1782943770000
kernel.py"""
scvi-tools kernel helpers. Load these definitions as described in SKILL.md:
h5ad_safe_obs
prepare_scvi_counts
Module top level is definition-only (functions + literal constants) and
stdlib-only so the sidecar AST gate accepts it and it loads under the
``python3 -I -S`` skeleton check; numpy/pandas/scipy are imported lazily inside
function bodies.
"""
def prepare_scvi_counts(adata):
"""Validate existing counts, using X only when the counts layer is absent.
Numerical checks cannot establish raw-count provenance. The caller must
identify the source from dataset documentation or preprocessing code.
Existing counts are never overwritten. No rounding, exponentiation, or
sparse-to-dense conversion is used to make an invalid input pass.
"""
import numpy as np
from scipy import sparse
if adata.isbacked or adata.is_view:
raise ValueError("Use an in-memory AnnData object, copying a view before preparation.")
source = "counts" if "counts" in adata.layers else "X"
matrix = adata.layers["counts"] if source == "counts" else adata.X
if matrix is None or matrix.shape != adata.shape or 0 in matrix.shape:
raise ValueError("Counts must be a nonempty cell-by-gene matrix aligned with AnnData.")
if sparse.issparse(matrix):
# CSR/CSC data can be checked directly; other SciPy formats remain sparse.
values = matrix.data if matrix.format in ("csr", "csc") else matrix.tocoo().data
elif isinstance(matrix, np.ndarray):
values = matrix
else:
raise ValueError("Load counts as a NumPy array or SciPy sparse matrix before validation.")
if values.dtype.kind not in "iuf":
raise ValueError("Raw counts must contain real numeric values, not booleans or strings.")
# Bound temporary allocations even for large dense expression matrices.
positive = False
for start in range(0, values.size, 65536):
block = values.flat[start:start + 65536]
if not np.isfinite(block).all() or (block < 0).any():
raise ValueError("Raw counts must be finite and nonnegative.")
if values.dtype.kind == "f" and (block != np.floor(block)).any():
raise ValueError("Fractional values are not raw integer counts; do not round or invert normalization.")
positive = positive or bool((block > 0).any())
if not positive:
raise ValueError("The count matrix contains no positive counts; check the input and QC.")
if "counts" not in adata.layers:
adata.layers["counts"] = matrix.copy()
return {"source": source, "checks": "finite, nonnegative, integer-valued"}
def h5ad_safe_obs(df):
"""Return a copy of `df` with index + string-like columns coerced so
anndata can `.write_h5ad()` without
``IORegistryError: No method registered for writing
<class 'pandas.arrays.ArrowStringArray'>`` (anndata #2377).
anndata 0.11.x's HDF5 writer has no registered method for
pyarrow-backed string columns (dtype ``string[pyarrow]`` /
``large_string[pyarrow]`` / dictionary-encoded), and
``anndata.settings.allow_write_nullable_strings`` gates only the
Python-backed ``pd.arrays.StringArray`` — it does not cover Arrow.
``.astype(str)`` alone is NOT enough: on a pyarrow-backed
Index/Series it returns another Arrow-backed array. Round-trip
through ``np.asarray(..., dtype=object)`` to force a plain object
array, then to ``Categorical``. Nulls are preserved (NA stays NA in
the categorical, not the literal string ``"<NA>"``/``"nan"``).
Typical use::
adata.obs = h5ad_safe_obs(adata.obs)
adata.var = h5ad_safe_obs(adata.var)
adata.write_h5ad("out.h5ad")
"""
import numpy as np
import pandas as pd
out = df.copy()
out.index = pd.Index(
np.asarray(out.index, dtype=object).astype(str), name=out.index.name
)
for c in out.columns:
dt = str(out[c].dtype)
# "string" substring covers string / string[python] / string[pyarrow]
# / large_string[pyarrow] / dictionary<values=string,...>[pyarrow]
# without sweeping in int64[pyarrow], double[pyarrow], bool[pyarrow],
# timestamp[...][pyarrow] etc.
if dt == "object" or "string" in dt:
mask = out[c].notna()
vals = np.asarray(out[c].astype(object).where(mask, None), dtype=object)
vals[mask.to_numpy()] = np.asarray(out[c][mask].astype(str), dtype=object)
out[c] = pd.Categorical(vals)
return out
SKILL.md---
name: scvi-tools
description: >
Probabilistic single-cell RNA-seq with scvi-tools — scVI for a
batch-corrected latent space, scANVI for semi-supervised label transfer,
and Bayesian differential expression. Reach for this skill to integrate
scRNA-seq batches, embed cells for clustering, transfer annotations from a
reference onto a query, or score differentially expressed genes per cluster.
For spatial deconvolution / mapping use the cell2location, DestVI, or
Tangram methods instead.
license: Apache-2.0
requirements: [gpu]
metadata:
display-name: scvi-tools
---
# scvi-tools — scVI / scANVI
scvi-tools (Gayoso et al. 2022, github.com/scverse/scvi-tools, BSD-3-Clause)
wraps a family
of deep generative models for single-cell omics. The scRNA-seq core is **scVI**
(unsupervised batch-corrected latent embedding) and **scANVI** (scVI + a
classifier head for semi-supervised cell-type label transfer). Both expect
**raw integer UMI counts** and emit a low-dimensional `X_scVI` / `X_scANVI`
that drops into the scanpy neighbors → leiden → umap pipeline.
## Setup (any agent, no API key)
This is a **pure skill** — `kernel.py` is deterministic Python and _you_ (the
base model) do all the reasoning. There is no `host` runtime and no LLM API.
The helpers are `prepare_scvi_counts` for input validation and `h5ad_safe_obs`
for obs/var frames that `anndata.write_h5ad()` can serialize. Load them once
per session in a Python cell:
```python
exec(open("scvi-tools/kernel.py", encoding="utf-8").read()) # path to this skill's kernel.py
```
Nothing auto-loads it outside Claude Science. Then call the helpers directly.
If a helper raises `NameError`, you haven't exec'd kernel.py.
Dependencies: `pip install scvi-tools scanpy anndata`. Training needs a
CUDA-capable GPU — see [Remote compute](#remote-compute-rent-a-gpu) to fall out
to a rented GPU when you don't have one locally.
## How to run
### scVI — batch-corrected latent space
`prepare_scvi_counts(adata)` checks and preserves existing `counts`. If absent,
it checks `.X` before copying it to `counts`. Values must be finite, nonnegative,
and integer-valued, with at least one positive count. Invalid existing counts
raise instead of being replaced by `.X`; sparse matrices stay sparse.
These numerical checks cannot prove raw-count provenance. Check the dataset
documentation or preprocessing history; do not round or exponentiate transformed
values to make them pass. If the history is unclear, resolve it before training.
The helper does not use `.raw.X`, which may be normalized or have a different
gene axis. Use an in-memory AnnData object; materialize backed data or copy views
only within the available memory budget.
```python
import scanpy as sc
import scvi
adata = sc.read_h5ad("dataset.h5ad")
# Verify count provenance from the input documentation or preprocessing history.
# Preserve counts if present; otherwise validate .X before copying it to counts.
counts_record = prepare_scvi_counts(adata)
print(counts_record) # numerical checks only; does not certify provenance
adata.X = adata.layers["counts"].copy() # derive plotting/HVG data from verified counts
sc.pp.normalize_total(adata); sc.pp.log1p(adata) # optional, for HVG / plotting only
sc.pp.highly_variable_genes(adata, n_top_genes=2000, batch_key="batch", subset=True)
scvi.model.SCVI.setup_anndata(adata, layer="counts", batch_key="batch")
model = scvi.model.SCVI(adata, n_latent=30)
model.train(max_epochs=200, early_stopping=True, accelerator="gpu", devices=1)
adata.obsm["X_scVI"] = model.get_latent_representation()
adata.layers["scvi_normalized"] = model.get_normalized_expression(library_size=1e4)
```
### scANVI — label transfer from a partially-annotated reference
```python
lvae = scvi.model.SCANVI.from_scvi_model(
model, labels_key="cell_type", unlabeled_category="Unknown",
)
lvae.train(max_epochs=20, n_samples_per_label=100, accelerator="gpu", devices=1)
adata.obsm["X_scANVI"] = lvae.get_latent_representation()
adata.obs["pred_cell_type"] = lvae.predict()
```
`accelerator="gpu", devices=1` is the PyTorch-Lightning spelling; the legacy
`use_gpu=` kwarg was **removed** in scvi-tools 1.x and now raises `TypeError`.
## Differential expression
```python
de = model.differential_expression(
groupby="leiden", group1="3", # group2=None → vs. all other cells
mode="change", delta=0.25,
)
top = de.sort_values("proba_de", ascending=False).head(50)
```
For one-vs-rest leave `group2` out — `"rest"` is scanpy's
`rank_genes_groups` convention, not scvi-tools'; here `group2` is a literal
category name and `"rest"` would match zero cells.
scvi-tools ≥1.4 defaults to `mode="vanilla"`, whose result columns are
exactly:
```
['proba_m1', 'proba_m2', 'bayes_factor', 'scale1', 'scale2', 'raw_mean1',
'raw_mean2', 'non_zeros_proportion1', 'non_zeros_proportion2',
'raw_normalized_mean1', 'raw_normalized_mean2', 'comparison', 'group1',
'group2']
```
— no `lfc_*`, no `proba_de`, no `is_de_fdr_*`. **Pass `mode="change"`** to
get `lfc_mean` / `lfc_median` / `proba_de` / `is_de_fdr_0.05`. Sort on
`proba_de` (or on `bayes_factor` if you deliberately stayed in vanilla
mode).
## Output format
| Key | What |
| --------------------------------- | ------------------------------------------------------ |
| `adata.obsm["X_scVI"]` | `n_cells × n_latent` batch-corrected embedding |
| `adata.obsm["X_scANVI"]` | label-aware embedding (better separates known classes) |
| `adata.obs["pred_cell_type"]` | scANVI predicted label per cell |
| `adata.layers["scvi_normalized"]` | decoded expression, library-size normalized |
| DE dataframe | per-gene `lfc_*` / `proba_de` (with `mode="change"`) |
## Remote compute (rent a GPU)
An A100-class GPU is recommended for >50k cells. Training is a plain Python
script (`pipeline.py`) that reads counts, trains scVI/scANVI, and writes the
output `.h5ad` — run it on whatever GPU you have (a local/cluster CUDA box, or a
serverless GPU host such as Modal). There is no Claude-Science compute broker
here; drive the GPU host directly.
**Modal** (serverless GPU) — wrap `pipeline.py` in a Modal app and run it with
the Modal CLI (`modal run pipeline.py`), which blocks until the job finishes, so
you read the result synchronously (no notification tool needed):
```python
# pipeline.py — run with: modal run pipeline.py
import modal
image = (modal.Image.debian_slim()
.pip_install("scvi-tools==1.4.2", "scanpy==1.11.5", "anndata==0.11.4"))
app = modal.App("scvi-run", image=image)
vol = modal.Volume.from_name("scvi-data", create_if_missing=True) # holds dataset.h5ad / out.h5ad
@app.function(gpu="A100", timeout=3600, volumes={"/data": vol})
def train():
import scanpy as sc, scvi # noqa
adata = sc.read_h5ad("/data/dataset.h5ad")
# ... setup_anndata / scVI / scANVI / DE — see the recipe above ...
adata.obs = h5ad_safe_obs(adata.obs) # paste the helper into THIS script (below)
adata.write_h5ad("/data/out.h5ad")
vol.commit()
@app.local_entrypoint()
def main():
train.remote() # blocks until done; then read /data/out.h5ad from the volume
```
Helpers loaded via `exec` in your **local** session (see **Setup**) are not
defined in a remote `pipeline.py`. Include `prepare_scvi_counts` from `kernel.py`
before the remote preprocessing/training recipe, and `h5ad_safe_obs` before
`.write_h5ad()`. Copy the helper definitions into that script or ship and load
`kernel.py` there; keep the documented count-source selection and validation.
For a local/cluster GPU, just run `pipeline.py` directly where CUDA is visible —
no wrapper needed. (For a fuller Modal workflow see the `remote-compute-modal`
skill.)
## Gotchas
| Gotcha | What happens / fix |
| -------------------------------------------------------------------------- | ---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- |
| `differential_expression()` defaults to `mode="vanilla"` (scvi-tools ≥1.4) | `KeyError: 'lfc_mean'` / `'proba_de'` when sorting — pass `mode="change"` to get `lfc_*`/`proba_de`/`is_de_fdr_*`; in vanilla mode sort on `bayes_factor`. |
| `adata.obs` index/columns are `string[pyarrow]` (`ArrowStringArray`) | `.write_h5ad()` dies with `IORegistryError: No method registered for writing <class 'pandas.arrays.ArrowStringArray'>` (anndata #2377). Coerce before writing: `adata.obs = h5ad_safe_obs(adata.obs)` (kernel helper — load it via `exec` locally, see Setup; inline the coercion in a remote `pipeline.py`). **`.astype(str)` alone is not enough** — on a pyarrow-backed Index/Series it returns another Arrow-backed array; round-trip through `np.asarray(..., dtype=object)`. `anndata.settings.allow_write_nullable_strings = True` does **not** cover Arrow-backed strings. |
| `use_gpu=` kwarg | Removed in 1.x → `TypeError: train() got an unexpected keyword argument 'use_gpu'`. Use `accelerator="gpu", devices=1`. |
| Log-normalized data fed to `setup_anndata` | Silent garbage — scVI's NB likelihood needs raw integer counts. Verify count provenance, call `prepare_scvi_counts`, and pass `layer="counts"`. Never overwrite counts from an unverified `.X`. |
## Troubleshooting
| Symptom | Fix |
| --------------------------------------------------------------------------------------------------------------- | ------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------ |
| `KeyError: 'lfc_mean'` (or `'proba_de'`, `'is_de_fdr_0.05'`) on DE result | Add `mode="change"` to `differential_expression()`; the default vanilla mode has no LFC columns. |
| `IORegistryError: No method registered for writing <class 'pandas.arrays.ArrowStringArray'>` on `.write_h5ad()` | `adata.obs = h5ad_safe_obs(adata.obs)` (and `adata.var` if needed) before writing. The `allow_write_nullable_strings` flag does not help here. |
| `TypeError: ... unexpected keyword argument 'use_gpu'` | Replace with `accelerator="gpu", devices=1`. |
| `ValueError: ... non-negative integers` / NB loss explodes | `layer="counts"` points at log/float data — restore raw counts. |
| `MisconfigurationException: No supported gpu backend found` | No CUDA visible — drop `accelerator`/`devices` to fall back to CPU, or dispatch via Remote compute. |
| `UnicodeEncodeError: 'ascii' codec can't encode character ...` writing a summary / printing | Container has no `LANG` so Python defaults to ASCII. Open files with `encoding="utf-8"` and/or `sys.stdout.reconfigure(encoding="utf-8")` at script top, or set `PYTHONIOENCODING=utf-8` in the image. |
---
**Next**: cluster on `X_scVI` with scanpy (`sc.pp.neighbors(use_rep="X_scVI")`
→ `sc.tl.leiden` → `sc.tl.umap`); for spatial deconvolution train
cell2location / DestVI / Tangram on the scRNA-seq reference.
test_kernel.py"""Run from the repository root: python -m unittest discover -s resources/skills/scvi-tools -v"""
import tempfile
import unittest
from pathlib import Path
import anndata as ad
import numpy as np
from scipy import sparse
from kernel import prepare_scvi_counts
class PrepareCountsTests(unittest.TestCase):
def test_existing_counts_survive_normalized_x_and_h5ad_roundtrip(self):
counts = np.array([[0, 2], [3, 1]])
data = ad.AnnData(np.log1p(counts), layers={"counts": counts.copy()})
with tempfile.TemporaryDirectory() as tmp:
path = Path(tmp) / "input.h5ad"
data.write_h5ad(path)
loaded = ad.read_h5ad(path)
existing = loaded.layers["counts"]
prepare_scvi_counts(loaded)
self.assertIs(loaded.layers["counts"], existing)
np.testing.assert_array_equal(loaded.X, np.log1p(counts))
np.testing.assert_array_equal(loaded.layers["counts"], counts)
loaded.X = loaded.layers["counts"].copy()
loaded.X[0, 0] = 99
np.testing.assert_array_equal(loaded.layers["counts"], counts)
# Preparation does not write back to the user's input file.
np.testing.assert_array_equal(ad.read_h5ad(path).X, np.log1p(counts))
def test_sparse_counts_remain_sparse_without_toarray(self):
from unittest.mock import patch
for cls in (sparse.csr_matrix, sparse.csc_matrix, sparse.csr_array, sparse.csc_array):
with self.subTest(format=cls.__name__):
matrix = cls(np.array([[0, 2], [3, 0]]))
data = ad.AnnData(matrix)
with patch.object(cls, "toarray", side_effect=AssertionError("densified")):
prepare_scvi_counts(data)
self.assertIsInstance(data.layers["counts"], cls)
self.assertEqual((data.layers["counts"] != matrix).nnz, 0)
def test_rejects_fractional_value_beyond_first_chunk(self):
matrix = np.ones((2, 40000))
matrix[-1, -1] = 0.5
data = ad.AnnData(matrix)
with self.assertRaisesRegex(ValueError, "Fractional"):
prepare_scvi_counts(data)
self.assertNotIn("counts", data.layers)
def test_invalid_existing_counts_do_not_fall_back_to_x(self):
data = ad.AnnData(np.ones((2, 2)), layers={"counts": np.full((2, 2), 0.5)})
with self.assertRaisesRegex(ValueError, "Fractional"):
prepare_scvi_counts(data)
np.testing.assert_array_equal(data.layers["counts"], np.full((2, 2), 0.5))
def test_absent_counts_defaults_to_x_without_using_raw(self):
for sparse_input in (False, True):
with self.subTest(sparse=sparse_input):
matrix = np.array([[1., 0.], [2., 3.]])
data = ad.AnnData(sparse.csr_matrix(matrix) if sparse_input else matrix.copy())
data.raw = ad.AnnData(np.full((2, 2), 99))
record = prepare_scvi_counts(data)
self.assertEqual(record["source"], "X")
if sparse_input:
self.assertTrue(sparse.issparse(data.layers["counts"]))
self.assertEqual((data.layers["counts"] != data.X).nnz, 0)
self.assertFalse(np.shares_memory(data.X.data, data.layers["counts"].data))
else:
np.testing.assert_array_equal(data.layers["counts"], matrix)
self.assertFalse(np.shares_memory(data.X, data.layers["counts"]))
def test_absent_counts_with_invalid_x_does_not_create_counts(self):
for value in (0.25, -1, float("nan"), float("inf")):
for sparse_input in (False, True):
with self.subTest(value=value, sparse=sparse_input):
matrix = np.array([[1., 0.], [2., value]])
data = ad.AnnData(sparse.csr_matrix(matrix) if sparse_input else matrix)
data.raw = ad.AnnData(np.ones((2, 2)))
with self.assertRaises(ValueError):
prepare_scvi_counts(data)
self.assertNotIn("counts", data.layers)
def test_rejects_empty_zero_and_nonnumeric_counts(self):
for matrix in (np.empty((0, 2)), np.zeros((2, 2)), np.ones((2, 2), dtype=bool),
np.ones((2, 2), dtype=complex), np.full((2, 2), "1")):
with self.subTest(dtype=matrix.dtype, shape=matrix.shape):
data = ad.AnnData(matrix)
with self.assertRaises(ValueError):
prepare_scvi_counts(data)
self.assertNotIn("counts", data.layers)
def test_rejects_views_and_backed_objects(self):
data = ad.AnnData(np.ones((2, 2)))
with self.assertRaisesRegex(ValueError, "in-memory"):
prepare_scvi_counts(data[:1])
with tempfile.TemporaryDirectory() as tmp:
path = Path(tmp) / "input.h5ad"
data.write_h5ad(path)
backed = ad.read_h5ad(path, backed="r")
try:
with self.assertRaisesRegex(ValueError, "in-memory"):
prepare_scvi_counts(backed)
finally:
backed.file.close()
if __name__ == "__main__":
unittest.main()