Skip to content

Commit 3be37bc

Browse files
Merge-in the dev-1.0 branch
2 parents fd21228 + 3ded6b9 commit 3be37bc

52 files changed

Lines changed: 18641 additions & 19654 deletions

Some content is hidden

Large Commits have some content hidden by default. Use the searchbox below for content that may be hidden.

.github/workflows/docs.yml

Lines changed: 9 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -6,6 +6,14 @@ on:
66
tags:
77
- '*'
88
merge_group:
9+
910
jobs:
1011
Docs:
11-
uses: tskit-dev/.github/.github/workflows/docs.yml@v15
12+
#Disabled until major code changes completed.
13+
# uses: tskit-dev/.github/.github/workflows/docs.yml@v15
14+
15+
# Placeholder
16+
runs-on: ubuntu-latest
17+
steps:
18+
- name: Placeholder
19+
run: echo "Docs broken, fix when code is ready!"

.github/workflows/tests.yml

Lines changed: 0 additions & 28 deletions
Original file line numberDiff line numberDiff line change
@@ -38,31 +38,3 @@ jobs:
3838
matrix:
3939
python: [ 3.11, 3.13 ]
4040
os: [ macos-latest, ubuntu-24.04, windows-latest ]
41-
42-
zarr-v3-compatibility:
43-
name: Test zarr v3 error
44-
runs-on: ubuntu-latest
45-
steps:
46-
- name: Checkout
47-
uses: actions/checkout@v4.2.2
48-
with:
49-
submodules: true
50-
51-
- name: Install uv and set the python version
52-
uses: astral-sh/setup-uv@v6
53-
with:
54-
python-version: "3.12"
55-
version: "0.10.0"
56-
57-
- name: Setup venv
58-
run: |
59-
uv venv
60-
61-
- name: Install zarr v3 and tsinfer
62-
run: |
63-
uv pip install -e .
64-
uv pip install --upgrade zarr
65-
66-
- name: Test zarr v3 error
67-
run: |
68-
uv run --no-sync python -c "import tsinfer" 2>&1 | grep -q "zarr version.*is not supported"

CLAUDE.md

Lines changed: 100 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,100 @@
1+
# CLAUDE.md
2+
3+
This file provides guidance to Claude Code (claude.ai/code) when working with code in this repository.
4+
5+
## Build & Development Setup
6+
7+
```bash
8+
uv sync # Install all dependencies (including dev groups)
9+
```
10+
11+
The package includes a C extension (`_tsinfer`) built from `lib/` sources via setuptools.
12+
`uv sync` compiles it automatically. If you modify C code, re-run `uv sync` to rebuild.
13+
14+
## Common Commands
15+
16+
```bash
17+
uv run pytest tests/ -v # Run all tests
18+
uv run pytest tests/test_matching.py # Run a single test file
19+
uv run pytest tests/test_matching.py::TestFoo::test_bar -v # Run a single test
20+
uv run pytest --skip-slow # Skip slow tests
21+
22+
uv run ruff check --fix # Lint Python code (auto-fix)
23+
uv run ruff format # Format Python code
24+
```
25+
26+
Line length is 89 characters (configured in pyproject.toml for both ruff and clang-format).
27+
28+
## Architecture
29+
30+
**tsinfer** infers tree sequences from genetic variation data stored in VCZ (Variant Call Zarr) format.
31+
32+
### Pipeline (three-stage inference)
33+
34+
The public API is in `tsinfer/__init__.py`, exposing three main functions from `pipeline.py`:
35+
36+
1. **`infer_ancestors`** (`ancestors.py`) — Generates ancestral haplotypes from sample genotypes. Two-pass chunk-aware approach: first computes per-site statistics, then builds ancestors using the C `_tsinfer.AncestorBuilder`. Output is an ancestor VCZ store.
37+
38+
2. **`match`** (`pipeline.py``matching.py`) — Matches ancestors against each other, then matches samples against the ancestor tree sequence. Uses the C `_tsinfer.AncestorMatcher` (Li & Stephens HMM) and `_tsinfer.TreeSequenceBuilder`. Ancestors are grouped for parallel matching via `grouping.py`.
39+
40+
3. **`post_process`** (`pipeline.py`) — Cleans up the raw inferred tree sequence (edge extension, parsimony-based refinement).
41+
42+
`run()` in `pipeline.py` chains all three stages.
43+
44+
### Key modules
45+
46+
- **`config.py`** — Dataclass-based configuration (`Config`, `AncestorsConfig`, `MatchConfig`, `PostProcessConfig`, `Source`)
47+
- **`vcz.py`** — VCZ/Zarr I/O layer; chunk-aware genotype loading (`get_genotypes_for_sites`)
48+
- **`grouping.py`** — Ancestor grouping and match job computation for parallel processing
49+
- **`matching.py`** — Core matching logic; `_ts_from_tsb` converts `TreeSequenceBuilder` to tskit tree sequence
50+
- **`ancestors.py`** — Ancestor generation with `InferenceSites` and `AncestorWriter`
51+
- **`tests/algorithm.py`** — Pure Python reference implementations of `AncestorBuilder`, `AncestorMatcher`, `TreeSequenceBuilder` (used for testing correctness against C)
52+
53+
### C extension (`_tsinfer`)
54+
55+
Source in `lib/`. Three main classes exposed to Python:
56+
- `AncestorBuilder` — builds inferred ancestors from genotype data
57+
- `AncestorMatcher` — Li & Stephens HMM matching algorithm
58+
- `TreeSequenceBuilder` — constructs tree sequences incrementally
59+
60+
Vendored dependencies in `lib/subprojects/`: tskit C library and kastore.
61+
62+
### Data flow
63+
64+
Sample VCZ → `infer_ancestors` → Ancestor VCZ → `match` → raw `tskit.TreeSequence``post_process` → final tree sequence
65+
66+
## Git usage
67+
68+
- Don't include Co-authored-By lines in git commits.
69+
70+
## Code Style
71+
72+
- Do not be overly defensive - defend only against circumstances that can
73+
occur within the current codebase.
74+
- Prefer dataclasses over tuples when returning multiple values.
75+
- Use explicit `None` comparisons: `if x is not None` not `if x`.
76+
- Zarr v3 is now used (dependency: `zarr>=3`).
77+
- Import all modules at the top of the file, not inside functions or methods.
78+
- Prefer importing a module and using module.function instead of
79+
using ``from module import function``. This applies to intra-package
80+
imports too: use ``from . import config`` then ``config.X``, not
81+
``from .config import X``. Exceptions: ``from typing import ...`` is
82+
acceptable; ``from .X import Y`` is acceptable in ``__init__.py`` for
83+
defining the public API. Use ``import concurrent.futures as cf``.
84+
- Use idiomatic pathlib.Path operations instead of os.path operations.
85+
- When a parameter has a computed default derived from another parameter,
86+
compute it once at the point of use (the leaf function), not at every
87+
layer in the call chain. Pass `None` through intermediate layers.
88+
- Use PEP 604 union syntax: `int | None`, not `Optional[int]`.
89+
- One `logger = logging.getLogger(__name__)` per module at top level.
90+
91+
## Tool use
92+
93+
- Use `uv run` for all Python tooling (never bare `python -m`)
94+
95+
## Testing
96+
97+
- Organise tests in classes, not flat functions. Use pytest fixtures for setup.
98+
- Test helpers are in `tests/helpers.py` (e.g., `make_sample_vcz`, `make_ancestor_vcz`)
99+
- `tests/algorithm.py` contains Python reference implementations used to verify C code
100+
- `msprime` is used to simulate test data

README.md

Lines changed: 10 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -3,6 +3,16 @@
33
[![License](https://img.shields.io/github/license/tskit-dev/tsinfer)](https://github.com/tskit-dev/tsinfer/blob/main/LICENSE) [![PyPI version](https://img.shields.io/pypi/v/tsinfer.svg)](https://pypi.org/project/tsinfer/) [![Supported Python Versions](https://img.shields.io/pypi/pyversions/tsinfer.svg)](https://pypi.org/project/tsinfer/) [![Docs Build](https://github.com/tskit-dev/tsinfer/actions/workflows/docs.yml/badge.svg)](https://github.com/tskit-dev/tsinfer/actions/workflows/docs.yml) [![Tests](https://github.com/tskit-dev/tsinfer/actions/workflows/tests.yml/badge.svg)](https://github.com/tskit-dev/tsinfer/actions/workflows/tests.yml) [![codecov](https://codecov.io/gh/tskit-dev/tsinfer/branch/main/graph/badge.svg)](https://codecov.io/gh/tskit-dev/tsinfer)
44

55

6+
## NOTICE
7+
8+
**WARNING** The main branch on GitHub is a development branch for the 1.0 release.
9+
It has many changes from the current 0.x series, and is largely incompatible.
10+
**DO NOT USE THIS CODE** for your research - it is not ready. Please use the
11+
released version of tsinfer available from PyPI and conda-forge.
12+
13+
14+
## README (to updated once 1.0 stabilises)
15+
616
Infer whole-genome tree sequences from genetic variation data. Tsinfer implements efficient algorithms to reconstruct ancestral haplotypes and recombination breakpoints, producing succinct tree sequences that capture shared ancestry across the genome. It scales to large cohorts and integrates cleanly with the broader tskit ecosystem for downstream statistics and analysis.
717

818
The documentation ([stable](https://tskit.dev/tsinfer/docs/stable/)[latest](https://tskit.dev/tsinfer/docs/latest/)) contains details of how to use this software, including [installation instructions](https://tskit.dev/tsinfer/docs/stable/installation.html).

_tsinfermodule.c

Lines changed: 2 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -104,7 +104,8 @@ AncestorBuilder_init(AncestorBuilder *self, PyObject *args, PyObject *kwds)
104104
int err;
105105
static char *kwlist[]
106106
= { "num_samples", "max_sites", "genotype_encoding", "mmap_fd", NULL };
107-
int num_samples, max_sites, genotype_encoding;
107+
int num_samples, max_sites;
108+
int genotype_encoding = 0;
108109
int flags = 0;
109110
int mmap_fd = -1;
110111

convert_hdf5.py

Lines changed: 0 additions & 23 deletions
This file was deleted.

example_config.toml

Lines changed: 153 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,153 @@
1+
# Example tsinfer config for 1000 Genomes Project chromosome 20.
2+
#
3+
# Assumes input data has been converted to VCZ format using bio2zarr:
4+
# bio2zarr vcf2zarr convert 1kgp_chr20.vcf.gz data/1kgp_chr20.vcz
5+
#
6+
# Run the full pipeline:
7+
# tsinfer run example_config.toml --threads 4 -v
8+
#
9+
# Or run steps individually:
10+
# tsinfer infer-ancestors example_config.toml --threads 4 -v
11+
# tsinfer match example_config.toml --threads 4 -v
12+
#
13+
# Paths are resolved relative to this file's directory.
14+
15+
# ============================================================================
16+
# Sources
17+
# ============================================================================
18+
# Each [[source]] block defines a named view over a VCZ store. The same store
19+
# can appear multiple times with different filters.
20+
21+
[[source]]
22+
name = "1kgp_chr20"
23+
path = "data/1kgp_chr20.vcz"
24+
25+
# Variant filters use bcftools-style expressions (via vcztools).
26+
# Uncomment to restrict to biallelic SNPs with QUAL >= 30:
27+
# include = "TYPE='snp' && N_ALT=1 && QUAL >= 30"
28+
29+
# Exclude specific variants:
30+
# exclude = "FILTER != 'PASS'"
31+
32+
# Restrict to a genomic region (half-open coordinates):
33+
# regions = "chr20:1000000-50000000"
34+
35+
# Restrict to exact target positions (useful for known-site lists):
36+
# targets = "chr20:1000000-50000000"
37+
38+
# Subset samples by name (comma-separated). Prefix with ^ to exclude:
39+
# samples = "HG00096,HG00097,HG00099"
40+
# samples = "^NA12878,NA12891" # all samples except these two
41+
42+
# Per-sample times. Use a field name in the VCZ store, a constant, or a
43+
# path to an external VCZ with a matching sample dimension:
44+
# sample_time = 0 # constant: all contemporary
45+
# sample_time = "variant_sample_time" # field in the source store
46+
# sample_time = { path = "ages.vcz", field = "sample_time" } # external
47+
48+
49+
# ============================================================================
50+
# Ancestral state
51+
# ============================================================================
52+
# Where to find the ancestral allele for each variant. Required when the
53+
# source VCZ does not contain a "variant_ancestral_allele" array.
54+
#
55+
# For 1000 Genomes, Ensembl provides ancestral allele annotations in a
56+
# separate VCF/VCZ. The "field" is the array name holding the allele string.
57+
58+
[ancestral_state]
59+
path = "data/homo_sapiens-chr20.vcz"
60+
field = "variant_AA"
61+
62+
63+
# ============================================================================
64+
# Ancestors
65+
# ============================================================================
66+
# Controls the ancestor-generation step (infer-ancestors).
67+
68+
[[ancestors]]
69+
name = "ancestors"
70+
path = "data/ancestors.vcz" # output path for the ancestor store
71+
sources = ["1kgp_chr20"] # which source(s) to build from
72+
73+
# Maximum gap (in base pairs) between adjacent inference sites before
74+
# splitting into separate intervals. Sites further apart than this are
75+
# processed independently, reducing memory for sparse regions.
76+
# Human chromosomes with uniform coverage: 500 kb works well.
77+
# For sparser data or non-contiguous targets, try a smaller value.
78+
max_gap_length = 500_000 # default: 500,000
79+
80+
# Genotype encoding for the C ancestor builder.
81+
# "eight_bit" — one byte per haplotype per site (default, broad compatibility)
82+
# "one_bit" — one bit per haplotype per site (~8x less memory, biallelic only)
83+
# Use one_bit for large cohorts (>10k samples) to reduce memory pressure.
84+
genotype_encoding = "eight_bit"
85+
86+
# Zarr chunk sizes for the output ancestor store. Smaller chunks reduce
87+
# peak memory when the match step streams ancestors, but add I/O overhead.
88+
# samples_chunk_size = 1000 # default: 1000 (ancestor dimension)
89+
# variants_chunk_size = 1000 # default: 1000 (site dimension)
90+
91+
# Zarr compression for the output ancestor store. Uses Blosc with bitshuffle.
92+
# compressor = "zstd" # default: "zstd" (blosc cname)
93+
# compression_level = 7 # default: 7 (0-9)
94+
95+
96+
# ============================================================================
97+
# Match
98+
# ============================================================================
99+
# Controls the HMM matching step: copies ancestors against the tree, then
100+
# copies sample haplotypes against the ancestor tree.
101+
102+
[match]
103+
output = "data/output.trees" # output tree sequence path
104+
105+
# Enable Viterbi path compression. Reduces memory and speeds up matching
106+
# at the cost of a slightly less optimal Viterbi path. Almost always
107+
# beneficial; disable only for debugging.
108+
# path_compression = true # default: true
109+
110+
# Number of worker threads for the match step.
111+
# num_threads = 1 # default: 1
112+
113+
# Per-source match parameters. Each [match.sources.<name>] block configures
114+
# how nodes from that source are represented in the output tree sequence.
115+
#
116+
# node_flags: tskit node flags (default: 1 = NODE_IS_SAMPLE)
117+
# create_individuals: group nodes into tskit individuals (default: true)
118+
119+
[match.sources.ancestors]
120+
node_flags = 0
121+
create_individuals = false
122+
123+
[match.sources.1kgp_chr20]
124+
# node_flags = 1 # default: 1 (NODE_IS_SAMPLE)
125+
# create_individuals = true # default: true
126+
127+
128+
# ============================================================================
129+
# Post-processing
130+
# ============================================================================
131+
# Optional cleanup applied after matching.
132+
133+
[post_process]
134+
# Split the ultimate ancestor (virtual root) into per-tree roots.
135+
# split_ultimate = true # default: true
136+
137+
# Erase flanking material outside each sample's first/last informative site.
138+
# erase_flanks = true # default: true
139+
140+
141+
# ============================================================================
142+
# Individual metadata
143+
# ============================================================================
144+
# Map VCZ sample-dimensioned arrays into tskit individual metadata.
145+
# "fields" maps tskit metadata keys to VCZ array names.
146+
# "population" names a VCZ array whose unique values become tskit populations.
147+
148+
# [individual_metadata]
149+
# population = "sample_population"
150+
# [individual_metadata.fields]
151+
# name = "sample_id"
152+
# sex = "sample_sex"
153+
# population_name = "sample_population"

examples/1kg_conf.toml

Lines changed: 41 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,41 @@
1+
2+
# wget https://ftp.1000genomes.ebi.ac.uk/vol1/ftp/data_collections/1000G_2504_high_coverage/working/20220422_3202_phased_SNV_INDEL_SV/1kGP_high_coverage_Illumina.chr20.filtered.SNV_INDEL_SV_phased_panel.vcf.gz
3+
# wget https://ftp.1000genomes.ebi.ac.uk/vol1/ftp/data_collections/1000G_2504_high_coverage/working/20220422_3202_phased_SNV_INDEL_SV/1kGP_high_coverage_Illumina.chr20.filtered.SNV_INDEL_SV_phased_panel.vcf.gz.tbi
4+
# vcf2zarr explode 1kGP_high_coverage_Illumina.chr20.filtered.SNV_INDEL_SV_phased_panel.vcf.gz 1kgp_chr20.icf -p 4
5+
# vcf2zarr encode 1kgp_chr20.icf 1kgp_chr20.vcz -p 4
6+
[[source]]
7+
name = "1kgp_chr20"
8+
path = "data/1kgp_chr20.vcz"
9+
# First 5Mb for development
10+
regions = "chr20:1-10000000"
11+
include = 'TYPE="snp"'
12+
13+
# wget https://ftp.ensembl.org/pub/release-115/variation/vcf/homo_sapiens/homo_sapiens-chr20.vcf.gz
14+
# wget https://ftp.ensembl.org/pub/release-115/variation/vcf/homo_sapiens/homo_sapiens-chr20.vcf.gz.csi
15+
# vcf2zarr explode homo_sapiens-chr20.vcf.gz homo_sapiens-chr20.icf -p 4
16+
# vcf2zarr encode homo_sapiens-chr20.icf homo_sapiens-chr20.vcz -p 4
17+
18+
[ancestral_state]
19+
path = "data/homo_sapiens-chr20.vcz"
20+
field = "variant_AA"
21+
22+
[[ancestors]]
23+
name = "ancestors"
24+
path = "data/ancestors.vcz" # output path for the ancestor store
25+
sources = ["1kgp_chr20"] # which source(s) to build from
26+
max_gap_length = 500_000
27+
genotype_encoding = "one_bit"
28+
# samples_chunk_size = 1000 # default: 1000 (ancestor dimension)
29+
# variants_chunk_size = 1000 # default: 1000 (site dimension)
30+
31+
[match]
32+
output = "data/output.trees" # output tree sequence path
33+
34+
[match.sources.ancestors]
35+
node_flags = 0
36+
create_individuals = false
37+
38+
[match.sources.1kgp_chr20]
39+
node_flags = 1 # default: 1 (NODE_IS_SAMPLE)
40+
create_individuals = true # default: true
41+

0 commit comments

Comments
 (0)