Skip to content

Commit a9ef2e4

Browse files
mikesshclaude
andcommitted
feat(cli,build): convert/downsample/filter/pool commands, Parquet output, uv bootstrap
CLI: add the legacy-style data commands — `convert` (any supported format → canonical TSV/Parquet), `downsample`, `filter` (coding/frequency/segment), `pool` (flat or incidence --join). `-o` is now format-aware (.parquet/.pq → Parquet, else TSV). Packaging: switch the dev bootstrap to uv (repo-local .venv, no conda); rewrite setup.sh to be uv-first with a `python -m venv` fallback and portable across bash and zsh. environment.yml demoted to optional (mmseqs2 + slow arda tests only). Docs/skills: v3.0.0 status line; README + docs/index + usage gain the new CLI commands, Parquet output, the uv install path, a Longitudinal-dynamics usage section, and the AnnData/zarr bridge note. CHANGELOG updated; test_cli covers the four new commands. Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
1 parent 29f415f commit a9ef2e4

10 files changed

Lines changed: 365 additions & 75 deletions

File tree

CHANGELOG.md

Lines changed: 12 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -35,6 +35,18 @@ Notable changes to vdjtools v2. Releases before 3.0.0 are recorded in the git ta
3535
samples) and `--cohort DIR` (one streamed pass over a pre-ingested Parquet cohort); the `overlap`
3636
command now pre-aggregates each sample once.
3737

38+
### Added — CLI & packaging
39+
40+
- New `vdjtools` subcommands: **`convert`** (read any supported format — native / AIRR / Parquet /
41+
MiXcr / MiGec / MiTCR / immunoSEQ / IMGT / Vidjil / RTCR / TRUST4 / arda — and write the canonical
42+
table), **`downsample`**, **`filter`** (coding / non-coding / frequency / V-J segment), and
43+
**`pool`** (flat pool or incidence `--join`).
44+
- Every command's `-o` is now **format-aware**: a `.parquet` / `.pq` path writes Parquet, anything
45+
else (or stdout) writes TSV.
46+
- Development switched to **uv** — one repo-local `.venv`, no conda. `setup.sh` is rewritten to be
47+
uv-first (with a `python -m venv` fallback) and **portable across bash and zsh**. `environment.yml`
48+
is now optional, needed only for MMseqs2 (arda's aligner) + the slow arda round-trip tests.
49+
3850
### Added — notebooks (marimo, `[examples]` extra)
3951

4052
- `notebooks/vaccination_tracking.py` — clonotype tracking + recapture model across YFV / influenza

CLAUDE.md

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -32,8 +32,8 @@ and are run ad hoc; they were never importable from here, so nothing in the libr
3232

3333
## Build / test / run
3434
```bash
35-
conda env create -f environment.yml && conda activate vdjtools # or reuse .venv
36-
pip install -e ".[dev,test]" # builds _core
35+
uv venv && source .venv/bin/activate && uv pip install -e ".[dev,test]" # builds _core (default)
36+
# or `bash setup.sh` (uv-first, portable bash/zsh); conda env.yml only for mmseqs2 + slow arda tests
3737
pytest tests/python -q
3838
cmake -S . -B build -DVDJTOOLS_TESTS=ON && cmake --build build && ctest --test-dir build
3939
sphinx-build -W --keep-going -b html docs docs/_build/html # docs gate (zero warnings)

README.md

Lines changed: 57 additions & 15 deletions
Original file line numberDiff line numberDiff line change
@@ -27,12 +27,13 @@ Built on the antigenomics ecosystem:
2727
[vdjmatch](https://github.com/antigenomics/vdjmatch) (overlap + TCRnet),
2828
[arda](https://github.com/antigenomics/arda) (AIRR annotation + markup repair).
2929

30-
> **Status: `v2.2.0`** — the native V(D)J model engine plus the full analytics suite (diversity,
31-
> overlap/TCRnet, preprocessing, biomarkers, single-cell), CDR features, and legacy-format ingestion
32-
> (MiXcr, MiGec, immunoSEQ, IMGT/HighV-QUEST, Vidjil, RTCR, TRUST4, arda). Clonotype columns follow the AIRR
33-
> **junction** convention (`junction_nt` / `junction_aa`). The legacy v1.x tool lives on the
34-
> [`legacy-1.x`](https://github.com/antigenomics/vdjtools/tree/legacy-1.x) branch and its releases
35-
> remain available under the repository tags (`v0.0.1``1.2.1`).
30+
> **Status: `v3.0.0`** — the native V(D)J model engine plus the full analytics suite (diversity,
31+
> overlap/TCRnet, preprocessing, biomarkers, single-cell), **longitudinal clonotype dynamics**
32+
> (paired expansion testing + the VDJtrack recapture model), CDR features, and legacy-format
33+
> ingestion (MiXcr, MiGec, immunoSEQ, IMGT/HighV-QUEST, Vidjil, RTCR, TRUST4, arda). Clonotype
34+
> columns follow the AIRR **junction** convention (`junction_nt` / `junction_aa`). The legacy v1.x
35+
> tool lives on the [`legacy-1.x`](https://github.com/antigenomics/vdjtools/tree/legacy-1.x) branch
36+
> and its releases remain available under the repository tags (`v0.0.1``1.2.1`).
3637
3738
## Install
3839

@@ -71,16 +72,23 @@ for germline lookup, Pgen, generation, or the analytics. Install it via conda/br
7172

7273
### Development
7374

75+
Uses [uv](https://docs.astral.sh/uv/) — one repo-local `.venv`, no conda:
76+
7477
```bash
75-
conda env create -f environment.yml # python + C++ toolchain + mmseqs2 (arda's aligner)
76-
conda activate vdjtools
77-
pip install -e ".[dev,test]" # builds the _core C++ extension
78+
uv venv && source .venv/bin/activate
79+
uv pip install -e ".[dev,test]" # builds the _core C++ extension (scikit-build-core)
7880
```
7981

80-
The conda env is a convenience, not a requirement — `pip install -e ".[dev,test]"` in any venv
81-
works. It supplies MMseqs2 so the slow arda annotation round-trips in the test suite run too.
82+
Or run the bootstrap script (portable across bash/zsh, uv-first with a `python -m venv` fallback):
83+
84+
```bash
85+
bash setup.sh --dev-parents --tests # or: zsh setup.sh
86+
```
8287

83-
Or run the bootstrap script: `bash setup.sh --dev-parents --tests`.
88+
You need a C++ toolchain (Xcode CLT on macOS, build-essential on Linux) for the native `_core`
89+
extension. MMseqs2 is arda's aligner, needed **only** for the annotation path and the slow arda
90+
round-trip tests — `brew install mmseqs2`, or use the optional `environment.yml` conda env which
91+
bundles it.
8492

8593
## Quickstart — recombination model engine
8694

@@ -130,15 +138,28 @@ vdjtools generate -m TRB -n 1000 -o gen.tsv # sample sequences (cf. olga-ge
130138
vdjtools pgen seqs.tsv -m TRB -o pgen.tsv # Pgen per CDR3 (cf. olga-compute_pgen)
131139
vdjtools pgen seqs.tsv -m TRB --mismatches 1 # + the Hamming-1 ball; --v-col/--j-col to condition
132140

141+
# data — convert any format to the canonical table (TSV, or Parquet by extension), preprocess
142+
vdjtools convert mixcr.txt.gz -o clones.parquet # MiXcr/immunoSEQ/AIRR/… → canonical Parquet
143+
vdjtools downsample clones.parquet 100000 -o ds.tsv
144+
vdjtools filter clones.parquet --coding --min-freq 1e-4 -o coding.tsv
145+
vdjtools pool s1.tsv s2.tsv s3.tsv --join --min-samples 2 -o joint.tsv
146+
133147
# repertoire analytics — sample files, or a cohort via -m/--metadata + --base-dir
134148
vdjtools diversity sampleA.tsv sampleB.tsv -o diversity.tsv
135149
vdjtools overlap *.tsv -o overlap.tsv
136150
vdjtools segment-usage *.tsv --segment v -o usage.tsv
137151
vdjtools spectratype *.tsv -o spectra.tsv
152+
vdjtools diversity -m metadata.txt --base-dir samples/ --threads 8 -o div.tsv # parallel cohort
153+
vdjtools spectratype --cohort cohort_parquet/ -o spectra.tsv # one streamed pass
154+
155+
# longitudinal — paired within-donor expansion test between two timepoints
156+
vdjtools dynamics day0.tsv day15.tsv -o tracked.tsv
138157
```
139158

140-
Native vdjtools and AIRR Rearrangement inputs are auto-detected; every command writes TSV to `-o`
141-
(or stdout, so it pipes). Run `vdjtools <command> --help` for options.
159+
Native vdjtools, AIRR Rearrangement, Parquet, and third-party inputs are auto-detected; every
160+
command writes to `-o`**TSV, or Parquet when the path ends in `.parquet` / `.pq`** — or to stdout
161+
(so it pipes). Cohort commands parallelise over samples with `-t/--threads` or stream a pre-ingested
162+
Parquet cohort with `--cohort`. Run `vdjtools <command> --help` for options.
142163

143164
## Analytics (Python API)
144165

@@ -175,6 +196,22 @@ usage = preprocess.correct_vj_usage(cohort, batch_col="batch", transform="sigmoi
175196
fixed = preprocess.apply_vj_correction(sampleA, usage, sample_id="A0")
176197
```
177198

199+
Longitudinal tracking — which clonotypes changed between two timepoints, and the VDJtrack recapture
200+
model (Pavlova, Zvyagin & Shugay 2024):
201+
202+
```python
203+
from vdjtools import dynamics
204+
205+
# paired within-donor test: emergent / expanded / persistent / contracted / vanishing
206+
tracked = dynamics.test_pair(day0, day15) # depth handled per-pair (effective N)
207+
grouped = dynamics.test_metaclonotypes(day0, day15, scope="1,0,0,1") # 1-Hamming CDR3 ball first
208+
called = dynamics.expansion_test(day0, day15) # edgeR NB-exact caller (log2FC + p)
209+
210+
# VDJtrack size-bucket recapture model — recapture fraction per clone-size class (Beta bands);
211+
# split by a group column + capture_test() for the group effect (see notebooks/vaccination_tracking.py)
212+
rates = dynamics.capture_rates(pre, post)
213+
```
214+
178215
Incidence-based clonotype association (Emerson 2017 / Howie 2015 / De Witt 2018 / Vlasova 2026)
179216
— a choice of test, condition, and co-occurrence — and single-cell paired-chain Pgen:
180217

@@ -232,7 +269,12 @@ suites (`RUN_BENCHMARK=1`).
232269
- **Overlap** — sample overlap and TCRnet (via vdjmatch/seqtree), similarity-aware overlap, clustering.
233270
- **Preprocess** — downsampling, error-correction, VJ-usage batch-effect correction, pooling/joining.
234271
- **Biomarker** — incidence association (Fisher / χ² / Bayesian / permutation) vs binary / HLA-allele / CMH-stratified conditions; α-β & same-chain co-occurrence pairing; metaclonotypes.
235-
- **Single-cell** — AIRR Cell / 10x interoperability, chain pairing + QC, and paired α/β Pgen.
272+
- **Dynamics** — longitudinal clonotype tracking between timepoints: the paired within-donor
273+
expansion test (emergent / expanded / persistent / contracted / vanishing), the VDJtrack
274+
size-bucket **recapture model**, metaclonotype-grouped testing, and an edgeR NB-exact caller
275+
([`vdjtools.dynamics`](python/vdjtools/dynamics)).
276+
- **Single-cell** — AIRR Cell / 10x interoperability, chain pairing + QC, paired α/β Pgen, and a
277+
`to_anndata` bridge into the scverse ecosystem (writes `.h5ad` / `.zarr` via AnnData).
236278

237279
## License
238280

docs/index.rst

Lines changed: 17 additions & 7 deletions
Original file line numberDiff line numberDiff line change
@@ -10,8 +10,9 @@ built on the antigenomics ecosystem (`seqtree <https://github.com/antigenomics/s
1010

1111
.. note::
1212

13-
**v2.2.0** — the native V(D)J model engine plus the full analytics suite (diversity, overlap/TCRnet,
14-
preprocessing, biomarkers, single-cell), CDR features, and legacy-format ingestion (MiXcr, MiGec,
13+
**v3.0.0** — the native V(D)J model engine plus the full analytics suite (diversity, overlap/TCRnet,
14+
preprocessing, biomarkers, single-cell), longitudinal clonotype **dynamics** (paired expansion
15+
testing + the VDJtrack recapture model), CDR features, and legacy-format ingestion (MiXcr, MiGec,
1516
immunoSEQ, IMGT/HighV-QUEST, Vidjil, RTCR, TRUST4, arda). Clonotype columns follow the AIRR **junction**
1617
convention (``junction_nt`` / ``junction_aa``). The legacy Groovy/Java vdjtools (v1.x) lives on the
1718
``legacy-1.x`` branch and its releases remain available under the repository's tags.
@@ -51,15 +52,21 @@ the repertoire analytics (sample files or a metadata table, like the legacy tool
5152
vdjtools models # list the bundled models
5253
vdjtools generate -m TRB -n 1000 -o gen.tsv # sample sequences (cf. olga-generate_sequences)
5354
vdjtools pgen seqs.tsv -m TRB -o pgen.tsv # Pgen per CDR3 (cf. olga-compute_pgen)
54-
vdjtools pgen seqs.tsv -m TRB --mismatches 1 # + the Hamming-1 ball; --v-col/--j-col to condition
55+
56+
vdjtools convert mixcr.txt.gz -o clones.parquet # any format → canonical TSV / Parquet
57+
vdjtools filter clones.parquet --coding --min-freq 1e-4 -o coding.tsv
58+
vdjtools downsample clones.parquet 100000 -o ds.tsv
5559
5660
vdjtools diversity sampleA.tsv sampleB.tsv -o diversity.tsv
5761
vdjtools overlap *.tsv -o overlap.tsv
5862
vdjtools segment-usage *.tsv --segment v -o usage.tsv
59-
vdjtools spectratype *.tsv -o spectra.tsv
63+
vdjtools spectratype --cohort cohort_parquet/ -o spectra.tsv # one streamed pass
64+
vdjtools dynamics day0.tsv day15.tsv -o tracked.tsv # paired expansion test
6065
61-
Native vdjtools and AIRR inputs are auto-detected; every command writes TSV to ``-o`` (or stdout).
62-
Run ``vdjtools <command> --help`` for options.
66+
Native vdjtools, AIRR, Parquet, and third-party inputs are auto-detected; every command writes to
67+
``-o`` — TSV, or Parquet when the path ends in ``.parquet`` / ``.pq`` — or to stdout. Cohort commands
68+
parallelise over samples with ``-t/--threads`` or stream a Parquet cohort with ``--cohort``. Run
69+
``vdjtools <command> --help`` for options.
6370

6471
Performance
6572
-----------
@@ -89,7 +96,10 @@ Capabilities (see the :doc:`API reference <api>` and the project ROADMAP):
8996
- **Overlap** — sample overlap and TCRnet (via vdjmatch/seqtree), similarity-aware overlap, clustering.
9097
- **Preprocess** — downsampling, error-correction, batch-effect correction, pooling/joining.
9198
- **Biomarker** — incidence-based association (Fisher) and metaclonotype grouping.
92-
- **Single-cell** — AIRR Cell / 10x interoperability, chain pairing + QC, and paired α/β Pgen.
99+
- **Dynamics** — longitudinal clonotype tracking: paired within-donor expansion testing, the VDJtrack
100+
size-bucket recapture model, metaclonotype-grouped testing, and an edgeR NB-exact caller.
101+
- **Single-cell** — AIRR Cell / 10x interoperability, chain pairing + QC, paired α/β Pgen, and an
102+
AnnData / scverse bridge (``.h5ad`` / ``.zarr``).
93103

94104
.. toctree::
95105
:hidden:

docs/usage.rst

Lines changed: 66 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -167,6 +167,56 @@ the paper's own FMBA covid TCRβ cohort (deep repertoires, ~3.3M reads/sample):
167167
explained by batch drops from η²≈0.11 to ≈0.002 while the grand-mean usage and per-sample read
168168
depth are preserved.
169169

170+
Longitudinal dynamics
171+
---------------------
172+
173+
:mod:`vdjtools.dynamics` compares two samples of the **same donor** — which clonotypes changed
174+
between timepoints. Depth is handled *per pair* via a two-step effective sample size (never by
175+
normalising a whole cohort to a common depth, which is not defined here):
176+
177+
.. code-block:: python
178+
179+
from vdjtools import dynamics, preprocess
180+
181+
day0 = demo_sample(seed=1)
182+
day15 = preprocess.downsample(day0, day0["duplicate_count"].sum() // 2, seed=7) # a later draw
183+
184+
# paired within-donor test: every clonotype -> emergent / expanded / persistent /
185+
# contracted / vanishing (or `untested` when the pair is too shallow to call it)
186+
tracked = dynamics.test_pair(day0, day15)
187+
tracked["dynamics"].value_counts()
188+
189+
# group near-identical CDR3s (a 1-Hamming ball, V/J-pinned) BEFORE testing, for power on
190+
# convergent expansions — scope "1,0,0,1" = 1 substitution, "1,1,1,1" = 1 edit (Levenshtein)
191+
grouped = dynamics.test_metaclonotypes(day0, day15, scope="1,0,0,1")
192+
193+
# complementary edgeR NB-exact caller (TMM + qCML common dispersion) -> log2FC + p
194+
called = dynamics.expansion_test(day0, day15)
195+
196+
The **recapture model** (Pavlova, Zvyagin & Shugay, *Front Immunol* 2024) asks a different
197+
question: do *particular* clonotypes reappear more than the background rate, and does that depend
198+
on clone size? Clonotypes are binned by size (singleton / doubleton / tripleton / large), each
199+
bin's recapture fraction gets a ``Beta`` credible interval, and a log-linear model tests the
200+
**group** effect:
201+
202+
.. code-block:: python
203+
204+
import polars as pl
205+
206+
# tag the clonotypes of interest (e.g. antigen-specific) vs the rest, in the PRE sample
207+
specific = set(day0["junction_aa"].head(20))
208+
pre = day0.with_columns(
209+
pl.when(pl.col("junction_aa").is_in(specific)).then(pl.lit("specific"))
210+
.otherwise(pl.lit("background")).alias("group"))
211+
212+
rates = dynamics.capture_rates(pre, day15, group_col="group") # recapture per (group, size)
213+
dynamics.capture_test(rates) # log-linear group effect + p
214+
215+
Across several donors, concatenate the per-donor ``capture_rates(..., donor=<id>)`` frames and use
216+
:func:`~vdjtools.dynamics.capture_paired_test` for the per-size-class paired t-test. The full
217+
worked example — sunken/alluvial tracking plots, a VDJdb overlay, and the capture ribbons — is
218+
``marimo edit notebooks/vaccination_tracking.py``.
219+
170220
Biomarker association
171221
---------------------
172222

@@ -346,17 +396,30 @@ surface.
346396
Command line
347397
------------
348398

349-
Every workflow above has a CLI counterpart; inputs are auto-detected and results are written
350-
as TSV to ``-o`` (or stdout):
399+
Every workflow above has a CLI counterpart; inputs are auto-detected and results are written to
400+
``-o`` — TSV, or **Parquet when the path ends in** ``.parquet`` / ``.pq`` — or to stdout:
351401

352402
.. code-block:: bash
353403
354404
vdjtools models # list the bundled models
355405
vdjtools generate -m TRB -n 1000 -o gen.tsv # (cf. olga-generate_sequences)
356406
vdjtools pgen seqs.tsv -m TRB -o pgen.tsv # (cf. olga-compute_pgen)
407+
408+
# data: convert any format to the canonical table, and preprocess
409+
vdjtools convert mixcr.txt.gz -o clones.parquet # → canonical Parquet (or .tsv)
410+
vdjtools filter clones.parquet --coding --min-freq 1e-4 -o coding.tsv
411+
vdjtools downsample clones.parquet 100000 -o ds.tsv
412+
vdjtools pool s1.tsv s2.tsv --join --min-samples 2 -o joint.tsv
413+
414+
# analytics (sample files, or -m metadata + --base-dir; parallel with -t, or stream --cohort)
357415
vdjtools diversity sampleA.tsv sampleB.tsv -o diversity.tsv
358416
vdjtools overlap *.tsv -o overlap.tsv
359417
vdjtools segment-usage *.tsv --segment v -o usage.tsv
360-
vdjtools spectratype *.tsv -o spectra.tsv
418+
vdjtools spectratype -m metadata.txt --base-dir samples/ -t 8 -o spectra.tsv
419+
vdjtools diversity --cohort cohort_parquet/ -o div.tsv
420+
421+
# longitudinal + enrichment
422+
vdjtools dynamics day0.tsv day15.tsv -o tracked.tsv # paired within-donor expansion test
423+
vdjtools tcrnet sample.tsv -o net.tsv # neighbourhood enrichment (control cohort)
361424
362425
Run ``vdjtools <command> --help`` for options.

environment.yml

Lines changed: 9 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -1,9 +1,15 @@
1-
# Dedicated conda environment for vdjtools.
2-
# Provides Python, a C++ toolchain (to build the _core pybind11 extension), and MMseqs2.
1+
# OPTIONAL conda environment for vdjtools — only if you want MMseqs2 (arda's aligner) and a
2+
# bundled C++ toolchain from conda. The default/recommended dev setup is uv (no conda needed):
3+
#
4+
# bash setup.sh # uv venv + editable install (builds _core); see setup.sh
5+
# # or by hand:
6+
# uv venv && source .venv/bin/activate && uv pip install -e ".[dev,test]"
7+
#
8+
# Use this file only when you need MMseqs2 for the annotation / slow arda round-trip tests:
39
#
410
# conda env create -f environment.yml
511
# conda activate vdjtools
6-
# pip install -e ".[dev,test]"
12+
# uv pip install -e ".[dev,test]" # (or: pip install -e ".[dev,test]")
713
#
814
# NB: this conda env is a convenience for DEVELOPMENT, not a requirement. arda (the germline
915
# reference + VDJ-model backend) is a plain pip base dependency, so `pip install vdjtools`

0 commit comments

Comments
 (0)