From 59d1ec93908c27b7b4b95fcac2ae0e4c07e11f23 Mon Sep 17 00:00:00 2001 From: Jerome Kelleher Date: Mon, 30 Mar 2026 15:57:12 +0100 Subject: [PATCH 01/11] Add CLI + TOML config quickstart docs MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Replace the Python API-focused usage page with a quickstart guide covering the VCF → VCZ → TOML config → CLI pipeline workflow. New files: - docs/quickstart.md: walkthrough of data prep, config, CLI commands, config reference, and CLI reference - docs/_static/example_data.vcf.gz: hand-crafted bgzipped VCF with 8 sites, 3 diploid samples, and AA INFO field - docs/_static/example_config.toml: minimal TOML config for the example --- docs/_static/example_config.toml | 27 +++ docs/_static/example_data.vcf.gz | Bin 0 -> 398 bytes docs/_toc.yml | 2 +- docs/quickstart.md | 337 +++++++++++++++++++++++++++++++ 4 files changed, 365 insertions(+), 1 deletion(-) create mode 100644 docs/_static/example_config.toml create mode 100644 docs/_static/example_data.vcf.gz create mode 100644 docs/quickstart.md diff --git a/docs/_static/example_config.toml b/docs/_static/example_config.toml new file mode 100644 index 00000000..d140e4a9 --- /dev/null +++ b/docs/_static/example_config.toml @@ -0,0 +1,27 @@ +# Example tsinfer config for a small dataset. +# +# Run: +# vcf2zarr convert example_data.vcf.gz example_data.vcz +# tsinfer run example_config.toml -v + +[[source]] +name = "example" +path = "example_data.vcz" + +[ancestral_state] +path = "example_data.vcz" +field = "variant_AA" + +[[ancestors]] +name = "ancestors" +path = "example_ancestors.vcz" +sources = ["example"] + +[match] +output = "example_output.trees" + +[match.sources.ancestors] +node_flags = 0 +create_individuals = false + +[match.sources.example] diff --git a/docs/_static/example_data.vcf.gz b/docs/_static/example_data.vcf.gz new file mode 100644 index 0000000000000000000000000000000000000000..67ff20b8038e6768aa1b9277505818dea9b80cdb GIT binary patch literal 398 zcmb2|=3rp}f&Xj_PR>jWw;6i-H}W+Zh`2s~9xf2s$@n5nfmkD#qaheR!rx5KQ~j`T`&7v)R||K4N@ZY%v2-OG%bV9Ef&emO?jWP{NxlZ zQJyKz$<6yViRE9udLV^8tfJ}Ciw`;LGdzSk%llibl}kS>Uo?$M*Afl!TIjSs!FA{I zL*CmS@t+hcecpEAaOv{7IW{r&K0=yTdaEa>^tG4l)$X>H;pU6vdfLW#mUYfKllG5V zHP;W{;J0C$$1vaDF8}}BJ93N+4Dx7xo63;ezniZ?fyedvFT*yqN67^lGP`D-kl$fp z{rlrKCsv`hO$Al&-d5h_uyEq%alALTkGUg^jjy5qkFDQwYf+gtvllFvCbn`^sh?2G zVB3D_#hDj-6&P1AIJ2Ccae(u-Tlc@-$x&_h!jn$rwB3;@HVLci&Yq TO7Lr6^v2>UX$EF+1b_$tNtK-2 literal 0 HcmV?d00001 diff --git a/docs/_toc.yml b/docs/_toc.yml index 6ba9d46d..c9345995 100644 --- a/docs/_toc.yml +++ b/docs/_toc.yml @@ -10,7 +10,7 @@ parts: - file: installation - caption: Usage chapters: - - file: usage + - file: quickstart - caption: Inference chapters: - file: inference diff --git a/docs/quickstart.md b/docs/quickstart.md new file mode 100644 index 00000000..24434b3f --- /dev/null +++ b/docs/quickstart.md @@ -0,0 +1,337 @@ +(sec_quickstart)= + +# Quickstart + +_Tsinfer_ infers [tree sequences](https://tskit.dev/tutorials/what_is.html) +from phased genetic variation data. Input data is stored in +[VCF Zarr](https://github.com/sgkit-dev/vcf-zarr-spec/) (.vcz) format, and +the pipeline is controlled by a TOML configuration file. + +The typical workflow is: + +1. Convert a bgzipped VCF to VCZ format +2. Write a TOML config describing inputs, outputs, and parameters +3. Run the pipeline via the `tsinfer` CLI +4. Analyse the resulting tree sequence with [tskit](https://tskit.dev/) + + +(sec_quickstart_preparing_data)= + +## Preparing input data + +_Tsinfer_ reads phased genotype data from `.vcz` stores. If you have a +bgzipped, indexed VCF, convert it using +[vcf2zarr](https://sgkit-dev.github.io/bio2zarr/vcf2zarr/overview.html): + +```bash +vcf2zarr convert mydata.vcf.gz mydata.vcz +``` + +### Ancestral states + +Each site used for inference requires a known **ancestral allele**. There are +several ways to provide this: + +- **AA INFO field in the VCF.** If your VCF has an `AA` (ancestral allele) INFO + field, `vcf2zarr` will store it as `variant_AA` in the `.vcz` store. You can + then reference it directly in the config: + ```toml + [ancestral_state] + path = "mydata.vcz" + field = "variant_AA" + ``` + +- **Separate ancestral allele VCZ.** If ancestral alleles come from a different + source (e.g. an Ensembl ancestral allele VCF), convert that to `.vcz` too and + point to it: + ```toml + [ancestral_state] + path = "ancestral_alleles.vcz" + field = "variant_AA" + ``` + +Sites where the ancestral allele is unknown or does not match any allele in the +data are automatically treated as _non-inference_ sites (see +{ref}`sec_quickstart_inference_sites`). + + +(sec_quickstart_config)= + +## Writing the config + +The TOML config tells _tsinfer_ where to find inputs, where to write outputs, +and what parameters to use. Here is a minimal example: + +```toml +# -- Sources: one or more named VCZ stores --------------------------------- +[[source]] +name = "example" +path = "example_data.vcz" + +# -- Ancestral state -------------------------------------------------------- +[ancestral_state] +path = "example_data.vcz" +field = "variant_AA" + +# -- Ancestors: output store for inferred ancestors ------------------------- +[[ancestors]] +name = "ancestors" +path = "example_ancestors.vcz" +sources = ["example"] + +# -- Match: HMM matching and tree sequence output -------------------------- +[match] +output = "example_output.trees" + +[match.sources.ancestors] +node_flags = 0 # ancestors are not samples +create_individuals = false + +[match.sources.example] +# node_flags = 1 # default: mark as samples +# create_individuals = true # default: group into individuals +``` + +### Sources + +Each `[[source]]` block names a VCZ store. You can filter variants and samples +using bcftools-style expressions: + +```toml +[[source]] +name = "1kgp_chr20" +path = "data/1kgp_chr20.vcz" +include = "TYPE='snp' && N_ALT=1" # biallelic SNPs only +exclude = "FILTER != 'PASS'" +regions = "chr20:1000000-50000000" # restrict to a region +samples = "^NA12878" # exclude specific samples +``` + +### Ancestors + +The `[[ancestors]]` block configures ancestor generation. Key options: + +| Field | Default | Description | +|-------|---------|-------------| +| `name` | (required) | Unique name for this ancestor set | +| `path` | (required) | Output path for the ancestor `.vcz` store | +| `sources` | (required) | List of source names to build ancestors from | +| `max_gap_length` | 500,000 | Split intervals at gaps wider than this (bp) | +| `genotype_encoding` | `"eight_bit"` | `"one_bit"` uses ~8x less memory (biallelic only) | + +### Match + +The `[match]` section controls HMM matching. Each source that should appear +in the output tree sequence needs a `[match.sources.]` sub-table. +Ancestors should have `node_flags = 0` and `create_individuals = false`. + +| Field | Default | Description | +|-------|---------|-------------| +| `output` | (required) | Output `.trees` file path | +| `path_compression` | `true` | Viterbi path compression | +| `workdir` | — | Directory for checkpoints (enables resume) | +| `keep_intermediates` | `false` | Keep per-group checkpoint files | + +### Post-processing (optional) + +```toml +[post_process] +split_ultimate = true # split virtual root into per-tree roots +erase_flanks = true # erase ancestry outside informative sites +``` + + +(sec_quickstart_running)= + +## Running the pipeline + +### Full pipeline + +Run all steps (infer ancestors, match, post-process) in one command: + +```bash +tsinfer run config.toml --threads 4 -v +``` + +### Individual steps + +For large datasets, you may want to run steps separately: + +```bash +# Step 1: Build ancestors +tsinfer infer-ancestors config.toml --threads 4 -v + +# Step 2: Match ancestors and samples +tsinfer match config.toml --threads 4 -v +``` + +Post-processing and site augmentation can also be run separately: + +```bash +tsinfer post-process config.toml --input raw.trees -v +tsinfer augment-sites config.toml --input output.trees --output final.trees +``` + +### Validating the config + +Before running, check that your config is valid and all paths resolve: + +```bash +tsinfer config check config.toml +``` + +To see the fully resolved config with all defaults filled in: + +```bash +tsinfer config show config.toml +``` + +### Common CLI options + +| Flag | Description | +|------|-------------| +| `-t, --threads N` | Number of worker threads (default: 1) | +| `-f, --force` | Overwrite existing output files | +| `-p, --progress` | Show progress bars | +| `-v` | Verbose logging (repeat for more: `-vv` for debug) | +| `-l, --log-file FILE` | Write logs to a file | + + +(sec_quickstart_inspecting)= + +## Inspecting the result + +The output is a standard [tskit](https://tskit.dev/) tree sequence. Load it in +Python to explore: + +```python +import tskit + +ts = tskit.load("example_output.trees") +print(f"{ts.num_trees} trees, {ts.num_samples} samples, {ts.num_sites} sites") + +# Draw the trees +ts.draw_svg(size=(600, 300), y_axis=True) +``` + +Each sample in the original VCZ file corresponds to an _individual_ in the tree +sequence. Since diploid individuals have two haploid genomes, +`ts.num_samples` will be twice the number of diploid individuals. + +:::{note} +By default, internal node times in the inferred tree sequence are +_not_ in years or generations — they reflect allele frequencies. To get +meaningful dates, use [tsdate](https://tskit.dev/software/tsdate.html). +Calculating branch-length statistics on uncalibrated trees will raise an error. +::: + + +(sec_quickstart_inference_sites)= + +## Inference sites + +Not all sites are used by _tsinfer_ for inferring the genealogy. These +_non-inference_ sites are still included in the final tree sequence, but their +mutations are placed by +[parsimony](https://tskit.dev/tskit/docs/stable/python-api.html#tskit.Tree.map_mutations). +Non-inference sites include: + +- **Fixed sites** — no variation between samples +- **Singletons** — only one genome carries the derived allele +- **Unknown ancestral state** — the ancestral allele does not match any allele + at the site +- **Multiallelic sites** — more than two alleles + +Additional sites can be excluded from inference using the `exclude` filter +in the `[[source]]` config, or by adding an `[augment_sites]` section to +place them separately via parsimony. + + +(sec_quickstart_config_reference)= + +## Config reference + +### `[[source]]` + +| Field | Type | Default | Description | +|-------|------|---------|-------------| +| `name` | string | (required) | Unique name for this source | +| `path` | string | (required) | Path to VCZ store | +| `include` | string | — | bcftools include expression | +| `exclude` | string | — | bcftools exclude expression | +| `samples` | string | — | Sample filter (comma-separated, `^` to exclude) | +| `regions` | string | — | Genomic region (half-open) | +| `targets` | string | — | Exact target positions | +| `sample_time` | various | — | Per-sample times: constant, field name, or `{path, field}` | + +### `[ancestral_state]` + +| Field | Type | Default | Description | +|-------|------|---------|-------------| +| `path` | string | (required) | Path to VCZ containing ancestral alleles | +| `field` | string | (required) | Array name (e.g. `"variant_AA"`) | + +### `[[ancestors]]` + +| Field | Type | Default | Description | +|-------|------|---------|-------------| +| `name` | string | (required) | Unique ancestor set name | +| `path` | string | (required) | Output VCZ path | +| `sources` | list | (required) | Source names to build from | +| `max_gap_length` | int | 500,000 | Split at gaps wider than this (bp) | +| `samples_chunk_size` | int | 100 | Zarr chunk size (ancestor dim) | +| `variants_chunk_size` | int | 50,000 | Zarr chunk size (site dim) | +| `compressor` | string | `"zstd"` | Blosc compressor name | +| `compression_level` | int | 7 | Compression level (0–9) | + +### `[match]` + +| Field | Type | Default | Description | +|-------|------|---------|-------------| +| `output` | string | (required) | Output `.trees` path | +| `path_compression` | bool | `true` | Enable Viterbi path compression | +| `reference_ts` | string | — | Reference tree sequence path | +| `workdir` | string | — | Checkpoint directory (enables resume) | +| `keep_intermediates` | bool | `false` | Keep per-group checkpoints | + +### `[match.sources.]` + +| Field | Type | Default | Description | +|-------|------|---------|-------------| +| `node_flags` | int | 1 | tskit node flags (1 = `NODE_IS_SAMPLE`) | +| `create_individuals` | bool | `true` | Group sample nodes into individuals | + +### `[post_process]` + +| Field | Type | Default | Description | +|-------|------|---------|-------------| +| `split_ultimate` | bool | `true` | Split virtual root into per-tree roots | +| `erase_flanks` | bool | `true` | Erase ancestry outside informative sites | + +### `[augment_sites]` + +| Field | Type | Default | Description | +|-------|------|---------|-------------| +| `sources` | list | (required) | Source names for parsimony placement | + +### `[individual_metadata]` + +| Field | Type | Default | Description | +|-------|------|---------|-------------| +| `population` | string | — | VCZ array whose unique values become populations | +| `fields.` | string | — | Map tskit metadata keys to VCZ arrays | + + +(sec_quickstart_cli_reference)= + +## CLI reference + +| Command | Description | +|---------|-------------| +| `tsinfer run CONFIG` | Run the full pipeline | +| `tsinfer infer-ancestors CONFIG` | Build ancestor VCZ from sample data | +| `tsinfer match CONFIG` | Match ancestors and samples against the tree | +| `tsinfer post-process CONFIG --input FILE` | Post-process the matched tree sequence | +| `tsinfer augment-sites CONFIG --input FILE --output FILE` | Place non-inference sites by parsimony | +| `tsinfer config show CONFIG` | Print resolved config with defaults | +| `tsinfer config check CONFIG` | Validate config and verify paths | From 796b8d045a61cf064ed2cd087327b9f75650392e Mon Sep 17 00:00:00 2001 From: Jerome Kelleher Date: Mon, 30 Mar 2026 16:58:19 +0100 Subject: [PATCH 02/11] Get docs building cleanly with new CLI+config workflow - Add config.md: standalone config reference (extracted from quickstart) - Rewrite cli.rst: use sphinx-click to auto-document the Click CLI - Trim quickstart.md: minimal walkthrough with cross-refs to config/CLI - Add usage.md back to build: strip jupytext headers, convert code-cell to static code blocks, add legacy warning banner - Add development warning banner to index.md - Remove stale pages: api.rst, large_scale.md, file_formats.rst, simulation-example.py, example_ancestral_state.fa - Fix broken cross-references in inference.md and usage.md - Switch sphinx-argparse to sphinx-click in deps and config --- docs/_config.yml | 2 +- docs/_static/example_ancestral_state.fa | 2 - docs/_toc.yml | 7 +- docs/api.rst | 89 -------- docs/cli.rst | 49 +---- docs/config.md | 116 ++++++++++ docs/file_formats.rst | 43 ---- docs/index.md | 5 + docs/inference.md | 6 +- docs/large_scale.md | 205 ------------------ docs/quickstart.md | 276 +++--------------------- docs/simulation-example.py | 38 ---- docs/usage.md | 89 ++++---- pyproject.toml | 2 +- uv.lock | 29 +-- 15 files changed, 225 insertions(+), 733 deletions(-) delete mode 100644 docs/_static/example_ancestral_state.fa delete mode 100644 docs/api.rst create mode 100644 docs/config.md delete mode 100644 docs/file_formats.rst delete mode 100644 docs/large_scale.md delete mode 100644 docs/simulation-example.py diff --git a/docs/_config.yml b/docs/_config.yml index a99feb3c..f4775224 100644 --- a/docs/_config.yml +++ b/docs/_config.yml @@ -32,7 +32,7 @@ sphinx: - sphinx.ext.viewcode - sphinx.ext.intersphinx - sphinx_issues - - sphinxarg.ext + - sphinx_click - IPython.sphinxext.ipython_console_highlighting config: diff --git a/docs/_static/example_ancestral_state.fa b/docs/_static/example_ancestral_state.fa deleted file mode 100644 index 7c403500..00000000 --- a/docs/_static/example_ancestral_state.fa +++ /dev/null @@ -1,2 +0,0 @@ ->chr1 -nnnnnnnnnnnnnnnGnnnnnnnnnnnnnnnnnnnnnnnnnnnGnnnnnCnnnnTnnnnnnnnnnnnnnnCnnnAnnnnnnnnnTnnnnnnnnnAnnnn \ No newline at end of file diff --git a/docs/_toc.yml b/docs/_toc.yml index c9345995..24aa5f1a 100644 --- a/docs/_toc.yml +++ b/docs/_toc.yml @@ -11,17 +11,14 @@ parts: - caption: Usage chapters: - file: quickstart + - file: config + - file: usage - caption: Inference chapters: - file: inference - - file: large_scale - caption: Interfaces chapters: - - file: api - file: cli -- caption: File Formats - chapters: - - file: file_formats - caption: Miscellaneous chapters: - file: development diff --git a/docs/api.rst b/docs/api.rst deleted file mode 100644 index e8a6c3bf..00000000 --- a/docs/api.rst +++ /dev/null @@ -1,89 +0,0 @@ -.. _sec_api: - -================= -API Documentation -================= - -.. _sec_api_file_formats: - - -++++++++++++ -Variant data -++++++++++++ - -.. autoclass:: tsinfer.VariantData - - -.. autofunction:: tsinfer.add_ancestral_state_array - -+++++++++++++ -Ancestor data -+++++++++++++ - -.. autofunction:: tsinfer.load - -.. autoclass:: tsinfer.AncestorData - :inherited-members: - -.. todo:: - - 1. Add documentation for the data attributes in read-mode. - - -.. _sec_api_file_inference: - -***************** -Running inference -***************** - -.. autofunction:: tsinfer.infer - -.. autofunction:: tsinfer.generate_ancestors - -.. autoclass:: tsinfer.GenotypeEncoding - :members: - -.. autofunction:: tsinfer.match_ancestors - -.. autofunction:: tsinfer.match_samples - -.. autofunction:: tsinfer.augment_ancestors - -.. autofunction:: tsinfer.post_process - -***************** -Batched inference -***************** - -.. autofunction:: tsinfer.match_ancestors_batch_init - -.. autofunction:: tsinfer.match_ancestors_batch_groups - -.. autofunction:: tsinfer.match_ancestors_batch_group_partition - -.. autofunction:: tsinfer.match_ancestors_batch_group_finalise - -.. autofunction:: tsinfer.match_ancestors_batch_finalise - -.. autofunction:: tsinfer.match_samples_batch_init - -.. autofunction:: tsinfer.match_samples_batch_partition - -.. autofunction:: tsinfer.match_samples_batch_finalise - - -***************** -Container classes -***************** - -.. autoclass:: tsinfer.Variant - -.. autoclass:: tsinfer.Site - - -********** -Exceptions -********** - -.. autoexception:: tsinfer.FileFormatError - diff --git a/docs/cli.rst b/docs/cli.rst index b7127f8a..c1fc31d8 100644 --- a/docs/cli.rst +++ b/docs/cli.rst @@ -1,49 +1,18 @@ -.. _sec_cli: +.. _sec_cli_reference: ====================== Command line interface ====================== -.. warning:: - - The command line interface only supports the deprecated SampleData format - used in tsinfer<0.4.0. - -The command line interface in ``tsinfer`` is intended to provide a convenient -interface to the high-level :ref:`API functionality `. There are two -equivalent ways to invoke this program: - -.. code-block:: bash - - $ tsinfer - -or +The ``tsinfer`` command line interface runs the inference pipeline using a +TOML configuration file. See the :ref:`quickstart ` for an +introduction and the :ref:`config reference ` for all +available options. .. code-block:: bash - $ python3 -m tsinfer - -The first form is more intuitive and works well most of the time. The second -form is useful when multiple versions of Python are installed or if the -:command:`tsinfer` executable is not installed on your path. - -The :command:`tsinfer` program has five subcommands: :command:`list` prints a -summary of the data held in one of tsinfer's :ref:`file formats `; -:command:`infer` runs the complete :ref:`inference process ` for a given -input SampleData file; and -:command:`generate-ancestors`, :command:`match-ancestors` and -:command:`match-samples` run the three parts of this inference -process as separate steps. Running the inference as separate steps like this -is recommended for large inferences as it allows for greater control over -the inference process. - -++++++++++++++++ -Argument details -++++++++++++++++ - -.. argparse:: - :module: tsinfer - :func: get_cli_parser - :prog: tsinfer - :nodefault: + $ tsinfer run config.toml --threads 4 -v +.. click:: tsinfer.cli:main + :prog: tsinfer + :nested: full diff --git a/docs/config.md b/docs/config.md new file mode 100644 index 00000000..a272798c --- /dev/null +++ b/docs/config.md @@ -0,0 +1,116 @@ +(sec_config_reference)= + +# Configuration reference + +_Tsinfer_ is configured via a TOML file passed to the CLI. Paths in the config +are resolved relative to the config file's directory. + +A complete annotated example is in +[example_config.toml](https://github.com/tskit-dev/tsinfer/blob/main/example_config.toml). + + +## `[[source]]` + +Each `[[source]]` block defines a named view over a VCZ store. The same store +can appear multiple times with different filters. + +| Field | Type | Default | Description | +|-------|------|---------|-------------| +| `name` | string | (required) | Unique name for this source | +| `path` | string | (required) | Path to VCZ store | +| `include` | string | — | bcftools include expression (e.g. `"TYPE='snp'"`) | +| `exclude` | string | — | bcftools exclude expression | +| `samples` | string | — | Sample filter (comma-separated; prefix `^` to exclude) | +| `regions` | string | — | Genomic region, half-open (e.g. `"chr20:1000-50000"`) | +| `targets` | string | — | Exact target positions | +| `sample_time` | various | — | Per-sample times: constant, field name, or `{path, field}` dict | + + +## `[ancestral_state]` + +Specifies where to read the ancestral allele for each variant position. + +| Field | Type | Default | Description | +|-------|------|---------|-------------| +| `path` | string | (required) | Path to VCZ containing ancestral alleles | +| `field` | string | (required) | Array name in the store (e.g. `"variant_AA"`) | + + +## `[[ancestors]]` + +Controls the ancestor-generation step (`infer-ancestors`). At least one +`[[ancestors]]` block is required unless `[match]` specifies a `reference_ts`. + +| Field | Type | Default | Description | +|-------|------|---------|-------------| +| `name` | string | (required) | Unique ancestor set name | +| `path` | string | (required) | Output VCZ path | +| `sources` | list[str] | (required) | Source names to build ancestors from | +| `max_gap_length` | int | 500,000 | Split intervals at gaps wider than this (bp) | +| `samples_chunk_size` | int | 100 | Zarr chunk size (ancestor dimension) | +| `variants_chunk_size` | int | 50,000 | Zarr chunk size (site dimension) | +| `compressor` | string | `"zstd"` | Blosc compressor name | +| `compression_level` | int | 7 | Compression level (0–9) | +| `genotype_encoding` | string | `"eight_bit"` | `"one_bit"` uses ~8x less memory (biallelic only) | + + +## `[match]` + +Controls the HMM matching step. + +| Field | Type | Default | Description | +|-------|------|---------|-------------| +| `output` | string | (required) | Output `.trees` file path | +| `path_compression` | bool | `true` | Enable Viterbi path compression | +| `reference_ts` | string | — | Reference tree sequence (skip ancestor generation) | +| `workdir` | string | — | Checkpoint directory (enables resume) | +| `keep_intermediates` | bool | `false` | Keep per-group checkpoint files | + + +## `[match.sources.]` + +Per-source parameters. Every source that should appear in the output tree +sequence needs an entry here. + +| Field | Type | Default | Description | +|-------|------|---------|-------------| +| `node_flags` | int | 1 | tskit node flags (`1` = `NODE_IS_SAMPLE`, `0` for ancestors) | +| `create_individuals` | bool | `true` | Group sample nodes into tskit individuals | + + +## `[post_process]` + +Optional cleanup applied after matching. + +| Field | Type | Default | Description | +|-------|------|---------|-------------| +| `split_ultimate` | bool | `true` | Split virtual root into per-tree roots | +| `erase_flanks` | bool | `true` | Erase ancestry outside informative sites | + + +## `[augment_sites]` + +Place non-inference sites via parsimony. + +| Field | Type | Default | Description | +|-------|------|---------|-------------| +| `sources` | list[str] | (required) | Source names for parsimony placement | + + +## `[individual_metadata]` + +Map VCZ sample-dimensioned arrays into tskit individual metadata. + +| Field | Type | Default | Description | +|-------|------|---------|-------------| +| `population` | string | — | VCZ array whose unique values become tskit populations | + +### `[individual_metadata.fields]` + +Each key becomes a tskit metadata field; the value names the VCZ array. + +```toml +[individual_metadata.fields] +name = "sample_id" +sex = "sample_sex" +``` diff --git a/docs/file_formats.rst b/docs/file_formats.rst deleted file mode 100644 index faaae08b..00000000 --- a/docs/file_formats.rst +++ /dev/null @@ -1,43 +0,0 @@ -.. _sec_file_formats: - -============ -File formats -============ - -``tsinfer`` uses the excellent `zarr library `_ -to encode data in a form that is both compact and efficient to process. -See the :ref:`API documentation ` for details on -how to construct and manipulate these files using Python. The -:ref:`tsinfer list ` command provides a way to print out a -summary of these files. - - -.. _sec_file_formats_ancestors: - -************** -Ancestors File -************** - -The ancestors file contains the ancestral haplotype data inferred from the -sample data in the :ref:`sec_inference_generate_ancestors` step. - -.. todo:: Document the structure of the ancestors file. - - -.. _sec_file_formats_tree_sequences: - -************** -Tree sequences -************** - -The goal of ``tsinfer`` is to infer correlated genealogies from variation -data, and it uses the very efficient `succinct tree sequence -`_ data structure -to encode this output. Please see the `tskit documentation -`_ for details on how to -process and manipulate such tree sequences. - -The intermediate ``.ancestors.trees`` file produced by the -:ref:`sec_inference_match_ancestors` step is also a -tree sequence and can be loaded and analysed using the -`tskit API `_. diff --git a/docs/index.md b/docs/index.md index e44db54a..dacb3052 100644 --- a/docs/index.md +++ b/docs/index.md @@ -1,3 +1,8 @@ +```{warning} +This documentation is under active development and may be incomplete or +inaccurate. The software API is not yet stable. +``` + # Welcome to tsinfer's documentation! This is the documentation for {program}`tsinfer`, a method for inferring correlated diff --git a/docs/inference.md b/docs/inference.md index 8648ac71..81fbb715 100644 --- a/docs/inference.md +++ b/docs/inference.md @@ -129,7 +129,7 @@ as input. This provides very fast access to genotype data, compressed using cutt [compression methods](http://numcodecs.readthedocs.io). The input sample haplotypes and related metadata are a fraction of the size of a compressed VCF and can be processed efficiently. VCF can be converted to VCF Zarr by the (bio2zarr)[https://sgkit-dev.github.io/bio2zarr] -package. See {ref}`sec_usage` for examples. +package. See the {ref}`quickstart ` for examples. (sec_inference_generate_ancestors)= @@ -138,7 +138,7 @@ package. See {ref}`sec_usage` for examples. The first step in a `tsinfer` inference process is to generate a large number of potential ancestors and to store these in an -{ref}`ancestors file `. The ancestors +ancestors file. The ancestors file conventionally ends with `.ancestors`. The ancestor generation algorithm is described in the Methods section of @@ -153,7 +153,7 @@ Describe the ancestor generation algorithm in more detail here ## Matching ancestors & samples After we have generated a set of potential ancestors and stored them in -an {ref}`ancestors file `, we then +an ancestors file, we then run two matching steps. First we match the ancestors against each other to generate an "ancestors tree sequence", then we match the samples against this ancestors tree sequence to generate the final result. diff --git a/docs/large_scale.md b/docs/large_scale.md deleted file mode 100644 index 17d3ee69..00000000 --- a/docs/large_scale.md +++ /dev/null @@ -1,205 +0,0 @@ ---- -jupytext: - text_representation: - extension: .md - format_name: myst - format_version: 0.12 - jupytext_version: 1.9.1 -kernelspec: - display_name: Python 3 - language: python - name: python3 ---- - -:::{currentmodule} tsinfer -::: - -(sec_large_scale)= - -# Large scale inference - -Generally, for up to a few thousand samples a single multi-core machine -can infer a tree sequence in a few days, hours, or even minutes. -However, _tsinfer_ has been successfully used with datasets up to half a million -samples, where ancestor and sample matching can take several CPU-years. -At this scale inference must be scaled across many machines. -_Tsinfer_ provides specific APIs to enable this. -Here we detail considerations and tips for each step of the -inference process to help you scale up your analysis. A snakemake pipeline -which implements this parallelisation scheme is available as -[tsinfer-snakemake](https://github.com/benjeffery/tsinfer-snakemake). - -(sec_large_scale_ancestor_generation)= - -## Data preparation - -For large scale inference the data must be in [VCF Zarr](https://github.com/sgkit-dev/vcf-zarr-spec) -format, read by the {class}`VariantData` class. [Bio2zarr](https://github.com/sgkit-dev/bio2zarr) -is recommended for conversion from VCF, and [sgkit](https://github.com/sgkit-dev/sgkit) can then -be used to perform initial filtering. - -:::{todo} -An upcoming tutorial will detail conversion from VCF to a VCF Zarr suitable for tsinfer. -::: - - -## Ancestor generation - -Ancestor generation is generally the fastest step in inference. It is not yet -parallelised out-of-core in tsinfer and must be performed on a single machine. -However it scales well on machines with -many cores and hyperthreading via the `num_threads` argument to -{meth}`generate_ancestors`. The limiting factor is often that the -entire genotype array for the contig being inferred needs to fit in RAM. -This is the high-water mark for memory usage in tsinfer. - -If your data consists of only biallelic sites, with no missingness, -the `genotype_encoding` argument can be set to -{class}`tsinfer.GenotypeEncoding.ONE_BIT` which reduces the memory footprint of -the genotype array by a factor of 8, such that the RAM needed is roughly -`num_sites * num_samples * ploidy / 8 bytes`. This memory optimisation -results in a surprisingly small increase in runtime. - -## Ancestor matching - -Ancestor matching is one of the more time consuming steps of inference. It -proceeds in groups, progressively growing the tree sequence with younger -ancestors. At each stage the parallelism is limited to the number of ancestors -whose possible inheritors are already matched, as all possible inheritors -of a sample must be matched in an earlier group. For a typical human data set -the number of samples per group varies from single digits up to approximately -the number of samples. -The plot below shows the number of ancestors matched in each group for a typical -human data set, earlier groups are older ancestors: - -```{figure} _static/ancestor_grouping.png -:width: 80% -``` - -There are five tsinfer API methods that can be used to parallelise ancestor -matching. - -The five methods are: - -1. {meth}`match_ancestors_batch_init` -2. {meth}`match_ancestors_batch_groups` -3. {meth}`match_ancestors_batch_group_partition` -4. {meth}`match_ancestors_batch_group_finalise` -5. {meth}`match_ancestors_batch_finalise` - -Initially {meth}`match_ancestors_batch_init` should be called to -set up the batch matching and to determine the groupings of ancestors. -This method writes a file `metadata.json` to the `work_dir` that contains -a JSON encoded dictionary with configuration for later steps, and a key -`ancestor_grouping` which is a list of dictionaries, each containing the -list of ancestors in that group (key:`ancestors`) and a proposed partioning of -those ancestors into sets that can be matched in parallel (key:`partitions`). -The dictionary is also returned by the method. -The partitioning is controlled by the `min_work_per_job` and `max_num_partitions` -arguments. For each group, ancestors are placed in a partition until the sum of their -lengths exceeds `min_work_per_job`, when a new partition is started. However, the -number of partitions is not allowed to exceed `max_num_partitions`. It is suggested -to set `max_num_partitions` to around 3-4x the number of worker nodes available, -and `min_work_per_job` to around 2,000,000 for a typical human data set. - -Groups vs partitions is a point of common confusion. Note that groups of ancestors -are matched serially, and each group is split into partitions that can be -matched in parallel. - -Each group is matched in turn, either by calling {meth}`match_ancestors_batch_groups` -to match without partitioning, or by calling {meth}`match_ancestors_batch_group_partition` -many times in parallel followed by a single call to {meth}`match_ancestors_batch_group_finalise`. -Each call to {meth}`match_ancestors_batch_groups` or {meth}`match_ancestors_batch_group_finalise` -outputs the tree sequence to `work_dir`, which is then used by the next group. The length of -the `ancestor_grouping` in the metadata dictionary determines the group numbers that these methods -will need to be called for, and the length of the `partitions` list in each group determines -the number of calls to {meth}`match_ancestors_batch_group_partition` that are needed (if any). - -{meth}`match_ancestors_batch_groups` matches groups, without partitioning, from -`group_index_start` (inclusively) to `group_index_end` (exclusively). Combining -many groups into one call reduces the overhead from job submission and start -up times, but note on job failure the process can only be resumed from the -last `group_index_end`. - -To match a single group in parallel, call {meth}`match_ancestors_batch_group_partition` -once for each partition listed in the `ancestor_grouping[group_index]['partitions']` list, -incrementing `partition_index`. This will match the ancestors, placing the match data in -the `working_dir`. Once all are complete a single call to -{meth}`match_ancestors_batch_group_finalise` will then insert the matches and -output the tree sequence to `work_dir`. - -Each call to {meth}`match_ancestors_batch_groups` and {meth}`match_ancestors_batch_group_finalise` results in a tree sequence being written to `work_dir`. -These tree sequences are essentially checkpoints from with the batch matching workflow -can be resumed on job failure. - -Finally after the final group, call {meth}`match_ancestors_batch_finalise` to -combine the groups into a single tree sequence. - -The partitioning in `metadata.json` does not have to be used for every group. As early groups are -not matching to a large tree sequence it is often faster to not partition the first half of the -groups, depending on job set up and queueing time on your cluster. - -Calls to {meth}`match_ancestors_batch_group_partition` will only use a single core, but -{meth}`match_ancestors_batch_groups` will use as many cores as `num_threads` is set to. -Therefore this value and cluster resources requested should be scaled with the number of ancestors, -which can be read from the metadata dictionary. - -As an example of how the API methods can be used together, suppose the metadata dictionary -created by {meth}`match_ancestors_batch_init` contains the following: - -```python -{ - "ancestor_grouping": [ - {"ancestors": [0, ... 9], "partitions": None}, - { - "ancestors": [10, ... 15], - "partitions": [[10, 11, 12], [13, 14, 15]] - }, - {"ancestors": [16, ... 19], "partitions": None}, - {"ancestors": [20, ... 25], "partitions": None}, - {"ancestors": [26, ... 30], "partitions": None}, - { - "ancestors": [31, ... 41], - "partitions": [[31, 32, 33, 34, 35, 36], [37, 38, 39, 40, 41]] - }, - {"ancestors": [42, ... 45], "partitions": None}, - {"ancestors": [46, ... 50], "partitions": None}, - { - "ancestors": [51, ... 65], - "partitions": [ - [51, 52, 53, 54], - [55, 56, 57, 58], - [59, 60, 61, 62, 63, 64, 65] - ] - }, - ] -} -``` -Then the flow could look like the following diagram: (calls on the same horizontal line can be -done in parallel, note that method names are shortened): - -```{figure} _static/example_flow.svg -:width: 80% -``` - -Note that groups 1, 5 and 8 can be partitioned, but only groups 5 and 8 are actually partitioned in this example, as stated above partitioning for groups is optional. Groups 0-4 are matched in one call, groups 6 and 7 are matched in two calls, but -could have been matched in one. By splitting 6 and 7 the flow makes an additional resume point in the case of job failure at the cost of job start up and queueing time. - - -## Sample matching - -Sample matching is far simpler than ancestor matching as it is essentially the same as a single group -of ancestors. There are three API methods that work together to enable distributed sample matching. - -1. {meth}`match_samples_batch_init` -2. {meth}`match_samples_batch_partition` -3. {meth}`match_samples_batch_finalise` - -{meth}`match_samples_batch_init` should be called to set up the batch matching and to determine the -groupings of samples. Similar to {meth}`match_ancestors_batch_init` it has a `min_work_per_job` argument to control the level of parallelism. The method writes a file -`metadata.json` to the directory `work_dir` that contains a JSON encoded dictionary with -configuration for later steps. This is also returned by the call. The `num_partitions` key in -this dictionary is the number of times {meth}`match_samples_batch_partition` will need -to be called, with each partition index as the `partition_index` argument. These calls can happen -in parallel and write match data to the `work_dir` which is then used by -{meth}`match_samples_batch_finalise` to output the tree sequence. \ No newline at end of file diff --git a/docs/quickstart.md b/docs/quickstart.md index 24434b3f..4cddd38e 100644 --- a/docs/quickstart.md +++ b/docs/quickstart.md @@ -15,8 +15,6 @@ The typical workflow is: 4. Analyse the resulting tree sequence with [tskit](https://tskit.dev/) -(sec_quickstart_preparing_data)= - ## Preparing input data _Tsinfer_ reads phased genotype data from `.vcz` stores. If you have a @@ -27,35 +25,12 @@ bgzipped, indexed VCF, convert it using vcf2zarr convert mydata.vcf.gz mydata.vcz ``` -### Ancestral states - -Each site used for inference requires a known **ancestral allele**. There are -several ways to provide this: - -- **AA INFO field in the VCF.** If your VCF has an `AA` (ancestral allele) INFO - field, `vcf2zarr` will store it as `variant_AA` in the `.vcz` store. You can - then reference it directly in the config: - ```toml - [ancestral_state] - path = "mydata.vcz" - field = "variant_AA" - ``` - -- **Separate ancestral allele VCZ.** If ancestral alleles come from a different - source (e.g. an Ensembl ancestral allele VCF), convert that to `.vcz` too and - point to it: - ```toml - [ancestral_state] - path = "ancestral_alleles.vcz" - field = "variant_AA" - ``` +Each site used for inference requires a known **ancestral allele**. If your VCF +has an `AA` INFO field, `vcf2zarr` stores it as `variant_AA` in the `.vcz` +store and you can reference it directly in the config. Alternatively, ancestral +alleles can come from a separate VCZ store. See the +{ref}`config reference ` for details. -Sites where the ancestral allele is unknown or does not match any allele in the -data are automatically treated as _non-inference_ sites (see -{ref}`sec_quickstart_inference_sites`). - - -(sec_quickstart_config)= ## Writing the config @@ -63,166 +38,81 @@ The TOML config tells _tsinfer_ where to find inputs, where to write outputs, and what parameters to use. Here is a minimal example: ```toml -# -- Sources: one or more named VCZ stores --------------------------------- [[source]] -name = "example" -path = "example_data.vcz" +name = "mydata" +path = "mydata.vcz" -# -- Ancestral state -------------------------------------------------------- [ancestral_state] -path = "example_data.vcz" +path = "mydata.vcz" field = "variant_AA" -# -- Ancestors: output store for inferred ancestors ------------------------- [[ancestors]] name = "ancestors" -path = "example_ancestors.vcz" -sources = ["example"] +path = "ancestors.vcz" +sources = ["mydata"] -# -- Match: HMM matching and tree sequence output -------------------------- [match] -output = "example_output.trees" +output = "output.trees" [match.sources.ancestors] -node_flags = 0 # ancestors are not samples +node_flags = 0 create_individuals = false -[match.sources.example] -# node_flags = 1 # default: mark as samples -# create_individuals = true # default: group into individuals -``` - -### Sources - -Each `[[source]]` block names a VCZ store. You can filter variants and samples -using bcftools-style expressions: - -```toml -[[source]] -name = "1kgp_chr20" -path = "data/1kgp_chr20.vcz" -include = "TYPE='snp' && N_ALT=1" # biallelic SNPs only -exclude = "FILTER != 'PASS'" -regions = "chr20:1000000-50000000" # restrict to a region -samples = "^NA12878" # exclude specific samples +[match.sources.mydata] ``` -### Ancestors - -The `[[ancestors]]` block configures ancestor generation. Key options: - -| Field | Default | Description | -|-------|---------|-------------| -| `name` | (required) | Unique name for this ancestor set | -| `path` | (required) | Output path for the ancestor `.vcz` store | -| `sources` | (required) | List of source names to build ancestors from | -| `max_gap_length` | 500,000 | Split intervals at gaps wider than this (bp) | -| `genotype_encoding` | `"eight_bit"` | `"one_bit"` uses ~8x less memory (biallelic only) | - -### Match - -The `[match]` section controls HMM matching. Each source that should appear -in the output tree sequence needs a `[match.sources.]` sub-table. -Ancestors should have `node_flags = 0` and `create_individuals = false`. +The `[[source]]` block names a VCZ store. `[ancestral_state]` says where to +find ancestral alleles. `[[ancestors]]` configures ancestor generation. +`[match]` controls HMM matching and output. Each source that should appear in +the output needs a `[match.sources.]` entry — ancestors use +`node_flags = 0` (not samples). For the full set of options see the +{ref}`config reference `. -| Field | Default | Description | -|-------|---------|-------------| -| `output` | (required) | Output `.trees` file path | -| `path_compression` | `true` | Viterbi path compression | -| `workdir` | — | Directory for checkpoints (enables resume) | -| `keep_intermediates` | `false` | Keep per-group checkpoint files | - -### Post-processing (optional) - -```toml -[post_process] -split_ultimate = true # split virtual root into per-tree roots -erase_flanks = true # erase ancestry outside informative sites -``` - - -(sec_quickstart_running)= ## Running the pipeline -### Full pipeline - -Run all steps (infer ancestors, match, post-process) in one command: +Run all steps in one command: ```bash tsinfer run config.toml --threads 4 -v ``` -### Individual steps - -For large datasets, you may want to run steps separately: +Or run steps individually: ```bash -# Step 1: Build ancestors tsinfer infer-ancestors config.toml --threads 4 -v - -# Step 2: Match ancestors and samples tsinfer match config.toml --threads 4 -v ``` -Post-processing and site augmentation can also be run separately: - -```bash -tsinfer post-process config.toml --input raw.trees -v -tsinfer augment-sites config.toml --input output.trees --output final.trees -``` - -### Validating the config - -Before running, check that your config is valid and all paths resolve: +Validate a config before running: ```bash tsinfer config check config.toml ``` -To see the fully resolved config with all defaults filled in: - -```bash -tsinfer config show config.toml -``` - -### Common CLI options - -| Flag | Description | -|------|-------------| -| `-t, --threads N` | Number of worker threads (default: 1) | -| `-f, --force` | Overwrite existing output files | -| `-p, --progress` | Show progress bars | -| `-v` | Verbose logging (repeat for more: `-vv` for debug) | -| `-l, --log-file FILE` | Write logs to a file | +See the {ref}`CLI reference ` for all commands and options. -(sec_quickstart_inspecting)= - ## Inspecting the result -The output is a standard [tskit](https://tskit.dev/) tree sequence. Load it in -Python to explore: +The output is a standard [tskit](https://tskit.dev/) tree sequence: ```python import tskit -ts = tskit.load("example_output.trees") +ts = tskit.load("output.trees") print(f"{ts.num_trees} trees, {ts.num_samples} samples, {ts.num_sites} sites") - -# Draw the trees ts.draw_svg(size=(600, 300), y_axis=True) ``` -Each sample in the original VCZ file corresponds to an _individual_ in the tree -sequence. Since diploid individuals have two haploid genomes, -`ts.num_samples` will be twice the number of diploid individuals. +Each diploid individual in the VCZ file corresponds to an _individual_ in the +tree sequence with two haploid sample nodes, so `ts.num_samples` is twice the +number of diploid individuals. :::{note} -By default, internal node times in the inferred tree sequence are -_not_ in years or generations — they reflect allele frequencies. To get -meaningful dates, use [tsdate](https://tskit.dev/software/tsdate.html). -Calculating branch-length statistics on uncalibrated trees will raise an error. +Internal node times are allele frequencies, not years or generations. Use +[tsdate](https://tskit.dev/software/tsdate.html) to add meaningful dates. +Branch-length statistics on uncalibrated trees will raise an error. ::: @@ -230,108 +120,12 @@ Calculating branch-length statistics on uncalibrated trees will raise an error. ## Inference sites -Not all sites are used by _tsinfer_ for inferring the genealogy. These -_non-inference_ sites are still included in the final tree sequence, but their -mutations are placed by +Not all sites are used for inferring the genealogy. _Non-inference_ sites are +included in the final tree sequence with mutations placed by [parsimony](https://tskit.dev/tskit/docs/stable/python-api.html#tskit.Tree.map_mutations). -Non-inference sites include: +These include: - **Fixed sites** — no variation between samples - **Singletons** — only one genome carries the derived allele -- **Unknown ancestral state** — the ancestral allele does not match any allele - at the site +- **Unknown ancestral state** — ancestral allele does not match any allele - **Multiallelic sites** — more than two alleles - -Additional sites can be excluded from inference using the `exclude` filter -in the `[[source]]` config, or by adding an `[augment_sites]` section to -place them separately via parsimony. - - -(sec_quickstart_config_reference)= - -## Config reference - -### `[[source]]` - -| Field | Type | Default | Description | -|-------|------|---------|-------------| -| `name` | string | (required) | Unique name for this source | -| `path` | string | (required) | Path to VCZ store | -| `include` | string | — | bcftools include expression | -| `exclude` | string | — | bcftools exclude expression | -| `samples` | string | — | Sample filter (comma-separated, `^` to exclude) | -| `regions` | string | — | Genomic region (half-open) | -| `targets` | string | — | Exact target positions | -| `sample_time` | various | — | Per-sample times: constant, field name, or `{path, field}` | - -### `[ancestral_state]` - -| Field | Type | Default | Description | -|-------|------|---------|-------------| -| `path` | string | (required) | Path to VCZ containing ancestral alleles | -| `field` | string | (required) | Array name (e.g. `"variant_AA"`) | - -### `[[ancestors]]` - -| Field | Type | Default | Description | -|-------|------|---------|-------------| -| `name` | string | (required) | Unique ancestor set name | -| `path` | string | (required) | Output VCZ path | -| `sources` | list | (required) | Source names to build from | -| `max_gap_length` | int | 500,000 | Split at gaps wider than this (bp) | -| `samples_chunk_size` | int | 100 | Zarr chunk size (ancestor dim) | -| `variants_chunk_size` | int | 50,000 | Zarr chunk size (site dim) | -| `compressor` | string | `"zstd"` | Blosc compressor name | -| `compression_level` | int | 7 | Compression level (0–9) | - -### `[match]` - -| Field | Type | Default | Description | -|-------|------|---------|-------------| -| `output` | string | (required) | Output `.trees` path | -| `path_compression` | bool | `true` | Enable Viterbi path compression | -| `reference_ts` | string | — | Reference tree sequence path | -| `workdir` | string | — | Checkpoint directory (enables resume) | -| `keep_intermediates` | bool | `false` | Keep per-group checkpoints | - -### `[match.sources.]` - -| Field | Type | Default | Description | -|-------|------|---------|-------------| -| `node_flags` | int | 1 | tskit node flags (1 = `NODE_IS_SAMPLE`) | -| `create_individuals` | bool | `true` | Group sample nodes into individuals | - -### `[post_process]` - -| Field | Type | Default | Description | -|-------|------|---------|-------------| -| `split_ultimate` | bool | `true` | Split virtual root into per-tree roots | -| `erase_flanks` | bool | `true` | Erase ancestry outside informative sites | - -### `[augment_sites]` - -| Field | Type | Default | Description | -|-------|------|---------|-------------| -| `sources` | list | (required) | Source names for parsimony placement | - -### `[individual_metadata]` - -| Field | Type | Default | Description | -|-------|------|---------|-------------| -| `population` | string | — | VCZ array whose unique values become populations | -| `fields.` | string | — | Map tskit metadata keys to VCZ arrays | - - -(sec_quickstart_cli_reference)= - -## CLI reference - -| Command | Description | -|---------|-------------| -| `tsinfer run CONFIG` | Run the full pipeline | -| `tsinfer infer-ancestors CONFIG` | Build ancestor VCZ from sample data | -| `tsinfer match CONFIG` | Match ancestors and samples against the tree | -| `tsinfer post-process CONFIG --input FILE` | Post-process the matched tree sequence | -| `tsinfer augment-sites CONFIG --input FILE --output FILE` | Place non-inference sites by parsimony | -| `tsinfer config show CONFIG` | Print resolved config with defaults | -| `tsinfer config check CONFIG` | Validate config and verify paths | diff --git a/docs/simulation-example.py b/docs/simulation-example.py deleted file mode 100644 index 29afa099..00000000 --- a/docs/simulation-example.py +++ /dev/null @@ -1,38 +0,0 @@ -import builtins -import subprocess -import sys - -import msprime -import numpy as np - -if getattr(builtins, "__IPYTHON__", False): # if running IPython: e.g. in a notebook - num_diploids, seq_len = 100, 10_000 - name = "notebook-simulation" -else: # Take parameters from the command-line - num_diploids, seq_len = int(sys.argv[1]), float(sys.argv[2]) - name = "cli-simulation" - -ts = msprime.sim_ancestry( - num_diploids, - population_size=10**4, - recombination_rate=1e-8, - sequence_length=seq_len, - random_seed=6, -) -ts = msprime.sim_mutations(ts, rate=1e-8, random_seed=7) -ts_name = name + "-source.trees" -ts.dump(ts_name) -print( - f"Simulated {ts.num_samples} samples over {seq_len / 1e6} Mb:", - f"{ts.num_trees} trees and {ts.num_sites} sites", -) - -# Convert to a zarr file: this should be easier once a tskit2zarr utility is made, see -# https://github.com/sgkit-dev/bio2zarr/issues/232 -np.save(f"{name}-AA.npy", [s.ancestral_state for s in ts.sites()]) -ret = subprocess.run( - "python -m bio2zarr tskit2zarr convert --force".split() + [ts_name, name + ".vcz"], - stderr=subprocess.DEVNULL if name == "notebook-simulation" else None, -) -if ret.returncode == 0: - print(f"Converted to {name}.vcz") diff --git a/docs/usage.md b/docs/usage.md index c4b7a20f..8da51824 100644 --- a/docs/usage.md +++ b/docs/usage.md @@ -1,23 +1,16 @@ ---- -jupytext: - text_representation: - extension: .md - format_name: myst - format_version: 0.12 - jupytext_version: 1.9.1 -kernelspec: - display_name: Python 3 - language: python - name: python3 ---- - -:::{currentmodule} tsinfer -::: +(sec_usage)= +# Usage (legacy) -(sec_usage)= +```{warning} +This page documents the **old Python API** (`VariantData`, `tsinfer.infer`, +etc.) which is no longer the recommended workflow. It is retained as a +reference while the content is ported to the new TOML config + CLI +interface described in the {ref}`quickstart `. -# Usage +Code examples on this page are **not executed** and may not work with the +current version of tsinfer. +``` (sec_usage_toy_example)= @@ -30,15 +23,14 @@ For simplicity, we first demonstrate using a pre-generated .vcz file (later, in a VCF using [vcf2zarr](https://sgkit-dev.github.io/bio2zarr/vcf2zarr/overview.html)). -```{code-cell} ipython3 +```python import zarr vcf_zarr = zarr.open("_static/example_data.vcz") ``` Here's what the genotypes stored in that datafile look like: -```{code-cell} -:"tags": ["remove-input"] +```python import numpy as np G = vcf_zarr['call_genotype'][:] # read full genotype matrix into memory positions = vcf_zarr['variant_position'][:] @@ -75,7 +67,7 @@ alleles exist, and the ancestral state is known. _Tsinfer_ produces a genealogy that could have given rise to the data set above, based on the sites that vary between the samples. To provide extra information to the algorithm, you must wrap the .vcz file in a lightweight -{class}`tsinfer.VariantData` object, using syntax like: +`tsinfer.VariantData` object, using syntax like: ```{code} python vdata = tsinfer.VariantData("file.vcz", ancestral_state=***, ...) @@ -83,7 +75,7 @@ vdata = tsinfer.VariantData("file.vcz", ancestral_state=***, ...) #### Ancestral states -Importantly, the {class}`~tsinfer.VariantData` object requires an +Importantly, the `VariantData` object requires an *ancestral state* to be provided for each site used in inference. There are many methods for determing ancestral states: details are outside the scope of this manual, but we have started a @@ -113,7 +105,7 @@ Ancestral states can be specified in several ways: * Using a single string of ancestral states, e.g. from a FASTA file. The string should cover the entire genetic sequence, such that the `i`th character in the string is taken as the ancestral state for an inference site at position `i`. In this - case, the {meth}`add_ancestral_state_array` method can be used to extract the states + case, the `add_ancestral_state_array` method can be used to extract the states and save them to the VCF Zarr dataset, under the name `ancestral_state`. Note that if, as is common, variant positions in the .vcz file are one-based (starting at 1), rather than zero-based, you should add a padding character at the start of the string. @@ -122,7 +114,7 @@ Below we illustrate the single string method, using a stored FASTA file. In this file, the 16th, 44th, 50th, 55th, 71st, 75th, 85th, and 95th characters are `G`, `G`, `C`, `T`, `C`, `A`, `T`, and `A` (note that the 85th character, `T`, does not match any of the alleles in the .vcz genotypes for position 85). -```{code-cell} +```python import tsinfer import zarr import pyfaidx @@ -165,15 +157,15 @@ step of _tsinfer_. ### Topology inference -Once our data is wrapped in a {class}`~tsinfer.VariantData` object, we can infer +Once our data is wrapped in a `VariantData` object, we can infer a {ref}`tree sequence` e.g. using -_tsinfer_'s {ref}`Python API`. Note that each sample in the original +_tsinfer_'s Python API. Note that each sample in the original .vcz file will correspond to an *individual* in the resulting tree sequence. Since these three individuals are diploid, the resulting tree sequence will have `ts.num_samples == 6` (unlike in a .vcz file, a "sample" in tskit refers to a haploid genome, not a diploid individual). -```{code-cell} ipython3 +```python inferred_ts = tsinfer.infer(vdata) print(f"Inferred a genetic genealogy for {inferred_ts.num_samples} (haploid) genomes") ``` @@ -182,7 +174,7 @@ And that's it: we now have a fully functional {class}`tskit.TreeSequence` object that we can interrogate in the usual ways. For example, we can look at the {meth}`variants` in the tree sequence: -```{code-cell} ipython3 +```python print("TS sample", "\t".join(str(s) for s in inferred_ts.samples()), sep="\t") print("TS individual", "\t".join(str(inferred_ts.node(s).individual) for s in inferred_ts.samples()), sep="\t") print("-" * 60) @@ -201,7 +193,7 @@ are identical. Apart from the imputation of losslessly encode any genetic variation data, regardless of the inferred topology. You can check this programatically if you want: -```{code-cell} ipython3 +```python import numpy as np for v_orig, v_inferred in zip(vdata.variants(), inferred_ts.variants()): if any( @@ -216,7 +208,7 @@ We can examine the inferred genetic genealogy, in the form of {ref}`local trees`. _Tsinfer_ has placed mutations on the genealogy to explain the observed genetic variation: -```{code-cell} ipython3 +```python mut_labels = { m.id: "{:g}: {}→{}".format( s.position, @@ -251,8 +243,7 @@ is simply the frequency of the shared derived allele(s) on which the ancestral s is based. For this reason, the time units are described as "uncalibrated" in the plot, and trying to calculate statistics based on branch lengths will raise an error: -```{code-cell} ipython3 -:tags: ["raises-exception"] +```python inferred_ts.diversity(mode="branch") ``` @@ -278,7 +269,7 @@ the ancestral states array should only specify alleles for the unmasked sites. Below, for instance, is an example of including only sites up to position six in the contig labelled "chr1" in the `example_data.vcz` file: -```{code-cell} +```python import numpy as np import zarr @@ -313,7 +304,7 @@ For simplicity here we'll use Python to simulate some data under the coalescent with recombination, using [msprime](https://msprime.readthedocs.io/en/stable/api.html#msprime.simulate): -```{code-cell} ipython3 +```python import builtins import sys @@ -421,7 +412,7 @@ Once we have our `.vcz` file created, running the inference is straightforward. (sec_usage_simulation_example_inference)= -```{code-cell} ipython3 +```python # Infer & save a ts from the notebook simulation. vcf_zarr = zarr.load(f"{name}.vcz") # currently must load the zarr to get ancestral states vdata = tsinfer.VariantData(f"{name}.vcz", ancestral_state=vcf_zarr["variant_allele"][:, 0]) @@ -471,7 +462,7 @@ inferred within the the notebook, but exactly the same code will work with the larger version run from the cli simulation. -```{code-cell} ipython3 +```python import tskit subset = range(0, 6) # show first 6 samples @@ -484,7 +475,7 @@ print(f"True tree seq, simplified to {len(subset)} sampled genomes") source_subset.draw_svg(size=(800, 200), x_lim=limit, time_scale="rank") ``` -```{code-cell} ipython3 +```python inferred = tskit.load(prefix + "-simulation.trees") inferred_subset = inferred.simplify(subset, filter_sites=False) print(f"Inferred tree seq, simplified to {len(subset)} sampled genomes") @@ -541,20 +532,16 @@ For example data, we use a publicly available VCF file of the genetic variants from chromosome 24 of ten Norwegian and French house sparrows, *Passer domesticus* (thanks to Mark Ravinet for the data file): -```{code-cell} ipython3 -:tags: ["remove-output"] -import zarr - -vcf_location = "_static/P_dom_chr24_phased.vcf.gz" -!python -m bio2zarr vcf2zarr convert --force {vcf_location} sparrows.vcz +```bash +vcf2zarr convert P_dom_chr24_phased.vcf.gz sparrows.vcz ``` This creates the `sparrows.vcz` datastore, which we open using -{class}`tsinfer.VariantData`. The original VCF had the ancestral allelic +`tsinfer.VariantData`. The original VCF had the ancestral allelic state specified in the `AA` INFO field, so we can simply provide the string `"variant_AA"` as the ancestral_state parameter. -```{code-cell} ipython3 +```python # Do the inference: this VCF has ancestral states in the AA field vdata = tsinfer.VariantData("sparrows.vcz", ancestral_state="variant_AA") ts = tsinfer.infer(vdata) @@ -575,7 +562,7 @@ This can be done by adding some descriptive metadata for each population, and th each sample to one of those populations. In our case, the sample sparrow IDs beginning with "FR" are from France: -```{code-cell} ipython3 +```python import json import numpy as np import tskit @@ -610,7 +597,7 @@ zarr.save("sparrows.vcz/individuals_population", individuals_population) Note that the steps above to generate a .vcz file are not strictly part of _tsinfer_. We only invoke _tsinfer_ subsequently, when creating a -{class}`~tsinfer.VariantData` object. Moreover, _tsinfer_ treats the +`VariantData` object. Moreover, _tsinfer_ treats the .vcz information as read-only, and does not make a copy of it. This means _tsinfer_ is well-suited to using publicly provided, read-only .vcz datafiles. Furthermore, {ref}`sec_usage_toy_example_masks` make it easy @@ -620,7 +607,7 @@ or more samples than are required for your analysis. As the .vcz file we are now using contains population metadata, `tsinfer` will create a tree sequence whose sample nodes are correctly assigned to named populations: -```{code-cell} ipython3 +```python vdata = tsinfer.VariantData("sparrows.vcz", ancestral_state="variant_AA", individuals_population="individuals_population") sparrow_ts = tsinfer.infer(vdata) @@ -651,7 +638,7 @@ code demonstrates how to use the {meth}`tskit.TreeSequence.at` method to obtain 1Mb from the start of the sequence, and plot it, colouring the tips according to population: -```{code-cell} ipython3 +```python colours = {"Norway": "red", "France": "blue"} colours_for_node = {} for n in sparrow_ts.samples(): @@ -692,7 +679,7 @@ calculated locally (e.g. {ref}`per tree or in genomic windows Date: Mon, 30 Mar 2026 16:58:34 +0100 Subject: [PATCH 03/11] Reenable docs workflow --- .github/workflows/docs.yml | 9 +-------- 1 file changed, 1 insertion(+), 8 deletions(-) diff --git a/.github/workflows/docs.yml b/.github/workflows/docs.yml index 6190a6aa..61200366 100644 --- a/.github/workflows/docs.yml +++ b/.github/workflows/docs.yml @@ -9,11 +9,4 @@ on: jobs: Docs: - #Disabled until major code changes completed. - # uses: tskit-dev/.github/.github/workflows/docs.yml@v15 - - # Placeholder - runs-on: ubuntu-latest - steps: - - name: Placeholder - run: echo "Docs broken, fix when code is ready!" + uses: tskit-dev/.github/.github/workflows/docs.yml@v15 From 6cbc89840c5b6b5e4e1d376cf609399cbff5d26c Mon Sep 17 00:00:00 2001 From: Jerome Kelleher Date: Mon, 30 Mar 2026 16:03:44 +0100 Subject: [PATCH 04/11] Validate topology in matcher_indexes_alloc to prevent segfaults Add matcher_indexes_validate_tables() that checks edge/mutation node and site IDs are in bounds, and that node 0 has at most one child per tree interval. Without these checks, out-of-bounds IDs are used as array indices causing heap buffer overflows (confirmed via valgrind). Add four C test cases covering: edge child/parent out of bounds, mutation node out of bounds, mutation site out of bounds, and multiple overlapping root edges. --- lib/ancestor_matcher.c | 60 +++++++++++++ lib/err.c | 6 ++ lib/err.h | 2 + lib/tests/tests.c | 187 +++++++++++++++++++++++++++++++++++++++++ 4 files changed, 255 insertions(+) diff --git a/lib/ancestor_matcher.c b/lib/ancestor_matcher.c index 7bcfa166..66467dd6 100644 --- a/lib/ancestor_matcher.c +++ b/lib/ancestor_matcher.c @@ -211,6 +211,61 @@ matcher_indexes_copy_mutation_data(matcher_indexes_t *self, return ret; } +static int +matcher_indexes_validate_tables( + size_t num_nodes, size_t num_sites, const tsk_table_collection_t *tables) +{ + tsk_size_t j; + const tsk_size_t num_edges = tables->edges.num_rows; + const tsk_size_t num_mutations = tables->mutations.num_rows; + const tsk_id_t *restrict edges_parent = tables->edges.parent; + const tsk_id_t *restrict edges_child = tables->edges.child; + const tsk_id_t *restrict mutations_node = tables->mutations.node; + const tsk_id_t *restrict mutations_site = tables->mutations.site; + const double *restrict edges_left = tables->edges.left; + const double *restrict edges_right = tables->edges.right; + /* Track whether node 0 already has an edge as parent. We check for + * overlapping intervals among edges with parent==0, which would mean + * multiple roots in the same tree. We only need to track one interval + * at a time because a single overlapping pair is enough to flag. */ + double root_edge_left = -1; + double root_edge_right = -1; + + for (j = 0; j < num_edges; j++) { + if (edges_parent[j] < 0 || edges_parent[j] >= (tsk_id_t) num_nodes + || edges_child[j] < 0 || edges_child[j] >= (tsk_id_t) num_nodes) { + return TSI_ERR_BAD_EDGE_NODE; + } + if (edges_parent[j] == 0) { + if (root_edge_left < 0) { + root_edge_left = edges_left[j]; + root_edge_right = edges_right[j]; + } else { + /* Check if this edge overlaps with a previous root edge */ + if (edges_left[j] < root_edge_right && edges_right[j] > root_edge_left) { + return TSI_ERR_MULTIPLE_ROOTS; + } + /* Extend the tracked interval */ + if (edges_left[j] < root_edge_left) { + root_edge_left = edges_left[j]; + } + if (edges_right[j] > root_edge_right) { + root_edge_right = edges_right[j]; + } + } + } + } + for (j = 0; j < num_mutations; j++) { + if (mutations_node[j] < 0 || mutations_node[j] >= (tsk_id_t) num_nodes) { + return TSI_ERR_BAD_MUTATION_NODE; + } + if (mutations_site[j] < 0 || mutations_site[j] >= (tsk_id_t) num_sites) { + return TSI_ERR_BAD_MUTATION_SITE; + } + } + return 0; +} + int matcher_indexes_alloc(matcher_indexes_t *self, const tsk_table_collection_t *tables, const tsk_size_t *num_alleles, tsk_flags_t flags) @@ -225,6 +280,11 @@ matcher_indexes_alloc(matcher_indexes_t *self, const tsk_table_collection_t *tab * list, so *don't* set from the tables */ self->num_mutations = 0; + ret = matcher_indexes_validate_tables(self->num_nodes, self->num_sites, tables); + if (ret != 0) { + goto out; + } + self->left_index_edges = malloc(self->num_edges * sizeof(*self->left_index_edges)); self->right_index_edges = malloc(self->num_edges * sizeof(*self->right_index_edges)); self->sites.mutations = malloc(self->num_sites * sizeof(*self->sites.mutations)); diff --git a/lib/err.c b/lib/err.c index 4a6f3b82..b9d5dc2b 100644 --- a/lib/err.c +++ b/lib/err.c @@ -116,6 +116,12 @@ tsi_strerror(int err) ret = "Sites with multiple mutations are not supported " "in matcher_indexes"; break; + case TSI_ERR_BAD_EDGE_NODE: + ret = "Bad edge: parent or child node out of bounds"; + break; + case TSI_ERR_MULTIPLE_ROOTS: + ret = "Node 0 must have at most one child per tree interval"; + break; } return ret; } diff --git a/lib/err.h b/lib/err.h index 08eeca41..6ba90f51 100644 --- a/lib/err.h +++ b/lib/err.h @@ -29,6 +29,8 @@ #define TSI_ERR_IO -25 #define TSI_ERR_BAD_ANCESTRAL_STATE -26 #define TSI_ERR_MULTIPLE_MUTATIONS_AT_SITE -27 +#define TSI_ERR_BAD_EDGE_NODE -28 +#define TSI_ERR_MULTIPLE_ROOTS -29 // clang-format on #ifdef __GNUC__ diff --git a/lib/tests/tests.c b/lib/tests/tests.c index 48340998..40e15b9f 100644 --- a/lib/tests/tests.c +++ b/lib/tests/tests.c @@ -441,6 +441,186 @@ test_matcher_indexes_errors(void) tsk_table_collection_free(&tables); } +static void +test_matcher_indexes_edge_node_out_of_bounds(void) +{ + int ret; + matcher_indexes_t mi; + tsk_table_collection_t tables; + + memset(&mi, 0, sizeof(mi)); + + /* Test 1: child out of bounds. + * Build valid tables, sort/index, then overwrite child to bypass + * tskit's own validation in sort/build_index. */ + ret = tsk_table_collection_init(&tables, 0); + CU_ASSERT_EQUAL_FATAL(ret, 0); + tables.sequence_length = 100; + + tsk_node_table_add_row(&tables.nodes, 0, 2.0, TSK_NULL, TSK_NULL, NULL, 0); + tsk_node_table_add_row(&tables.nodes, 0, 1.0, TSK_NULL, TSK_NULL, NULL, 0); + tsk_node_table_add_row(&tables.nodes, 0, 0.5, TSK_NULL, TSK_NULL, NULL, 0); + + tsk_edge_table_add_row(&tables.edges, 0, 100, 0, 1, NULL, 0); + tsk_edge_table_add_row(&tables.edges, 0, 100, 1, 2, NULL, 0); + + tsk_site_table_add_row(&tables.sites, 10, "A", 1, NULL, 0); + tsk_mutation_table_add_row( + &tables.mutations, 0, 1, TSK_NULL, TSK_UNKNOWN_TIME, "T", 1, NULL, 0); + + ret = tsk_table_collection_sort(&tables, NULL, 0); + CU_ASSERT_EQUAL_FATAL(ret, 0); + ret = tsk_table_collection_build_index(&tables, 0); + CU_ASSERT_EQUAL_FATAL(ret, 0); + + /* Overwrite child of second edge to out-of-bounds value */ + tables.edges.child[1] = 5; + + ret = matcher_indexes_alloc(&mi, &tables, NULL, 0); + CU_ASSERT_EQUAL_FATAL(ret, TSI_ERR_BAD_EDGE_NODE); + matcher_indexes_free(&mi); + tsk_table_collection_free(&tables); + + /* Test 2: parent out of bounds */ + memset(&mi, 0, sizeof(mi)); + ret = tsk_table_collection_init(&tables, 0); + CU_ASSERT_EQUAL_FATAL(ret, 0); + tables.sequence_length = 100; + + tsk_node_table_add_row(&tables.nodes, 0, 2.0, TSK_NULL, TSK_NULL, NULL, 0); + tsk_node_table_add_row(&tables.nodes, 0, 1.0, TSK_NULL, TSK_NULL, NULL, 0); + tsk_node_table_add_row(&tables.nodes, 0, 0.5, TSK_NULL, TSK_NULL, NULL, 0); + + tsk_edge_table_add_row(&tables.edges, 0, 100, 0, 1, NULL, 0); + tsk_edge_table_add_row(&tables.edges, 0, 100, 1, 2, NULL, 0); + + tsk_site_table_add_row(&tables.sites, 10, "A", 1, NULL, 0); + tsk_mutation_table_add_row( + &tables.mutations, 0, 1, TSK_NULL, TSK_UNKNOWN_TIME, "T", 1, NULL, 0); + + ret = tsk_table_collection_sort(&tables, NULL, 0); + CU_ASSERT_EQUAL_FATAL(ret, 0); + ret = tsk_table_collection_build_index(&tables, 0); + CU_ASSERT_EQUAL_FATAL(ret, 0); + + /* Overwrite parent of second edge to out-of-bounds value */ + tables.edges.parent[1] = 10; + + ret = matcher_indexes_alloc(&mi, &tables, NULL, 0); + CU_ASSERT_EQUAL_FATAL(ret, TSI_ERR_BAD_EDGE_NODE); + matcher_indexes_free(&mi); + tsk_table_collection_free(&tables); +} + +static void +test_matcher_indexes_mutation_node_out_of_bounds(void) +{ + int ret; + matcher_indexes_t mi; + tsk_table_collection_t tables; + + memset(&mi, 0, sizeof(mi)); + ret = tsk_table_collection_init(&tables, 0); + CU_ASSERT_EQUAL_FATAL(ret, 0); + tables.sequence_length = 100; + + tsk_node_table_add_row(&tables.nodes, 0, 2.0, TSK_NULL, TSK_NULL, NULL, 0); + tsk_node_table_add_row(&tables.nodes, 0, 1.0, TSK_NULL, TSK_NULL, NULL, 0); + tsk_node_table_add_row(&tables.nodes, 0, 0.5, TSK_NULL, TSK_NULL, NULL, 0); + + tsk_edge_table_add_row(&tables.edges, 0, 100, 0, 1, NULL, 0); + + tsk_site_table_add_row(&tables.sites, 10, "A", 1, NULL, 0); + /* Add valid mutation, then overwrite node to out-of-bounds value */ + tsk_mutation_table_add_row( + &tables.mutations, 0, 1, TSK_NULL, TSK_UNKNOWN_TIME, "T", 1, NULL, 0); + + ret = tsk_table_collection_sort(&tables, NULL, 0); + CU_ASSERT_EQUAL_FATAL(ret, 0); + ret = tsk_table_collection_build_index(&tables, 0); + CU_ASSERT_EQUAL_FATAL(ret, 0); + + /* Overwrite mutation node to out-of-bounds value */ + tables.mutations.node[0] = 99; + + ret = matcher_indexes_alloc(&mi, &tables, NULL, 0); + CU_ASSERT_EQUAL_FATAL(ret, TSI_ERR_BAD_MUTATION_NODE); + matcher_indexes_free(&mi); + tsk_table_collection_free(&tables); +} + +static void +test_matcher_indexes_mutation_site_out_of_bounds(void) +{ + int ret; + matcher_indexes_t mi; + tsk_table_collection_t tables; + + memset(&mi, 0, sizeof(mi)); + ret = tsk_table_collection_init(&tables, 0); + CU_ASSERT_EQUAL_FATAL(ret, 0); + tables.sequence_length = 100; + + tsk_node_table_add_row(&tables.nodes, 0, 2.0, TSK_NULL, TSK_NULL, NULL, 0); + tsk_node_table_add_row(&tables.nodes, 0, 1.0, TSK_NULL, TSK_NULL, NULL, 0); + tsk_node_table_add_row(&tables.nodes, 0, 0.5, TSK_NULL, TSK_NULL, NULL, 0); + + tsk_edge_table_add_row(&tables.edges, 0, 100, 0, 1, NULL, 0); + + tsk_site_table_add_row(&tables.sites, 10, "A", 1, NULL, 0); + /* Add a valid mutation, then overwrite site to out-of-bounds value */ + tsk_mutation_table_add_row( + &tables.mutations, 0, 1, TSK_NULL, TSK_UNKNOWN_TIME, "T", 1, NULL, 0); + + ret = tsk_table_collection_sort(&tables, NULL, 0); + CU_ASSERT_EQUAL_FATAL(ret, 0); + ret = tsk_table_collection_build_index(&tables, 0); + CU_ASSERT_EQUAL_FATAL(ret, 0); + + /* Overwrite site ID after sort/index to bypass tskit's checks */ + tables.mutations.site[0] = 5; + + ret = matcher_indexes_alloc(&mi, &tables, NULL, 0); + CU_ASSERT_EQUAL_FATAL(ret, TSI_ERR_BAD_MUTATION_SITE); + matcher_indexes_free(&mi); + tsk_table_collection_free(&tables); +} + +static void +test_matcher_indexes_multiple_roots(void) +{ + int ret; + matcher_indexes_t mi; + tsk_table_collection_t tables; + + memset(&mi, 0, sizeof(mi)); + ret = tsk_table_collection_init(&tables, 0); + CU_ASSERT_EQUAL_FATAL(ret, 0); + tables.sequence_length = 100; + + tsk_node_table_add_row(&tables.nodes, 0, 2.0, TSK_NULL, TSK_NULL, NULL, 0); + tsk_node_table_add_row(&tables.nodes, 0, 1.0, TSK_NULL, TSK_NULL, NULL, 0); + tsk_node_table_add_row(&tables.nodes, 0, 0.5, TSK_NULL, TSK_NULL, NULL, 0); + + /* Two edges with parent=0 overlapping over [0,100): multiple roots */ + tsk_edge_table_add_row(&tables.edges, 0, 100, 0, 1, NULL, 0); + tsk_edge_table_add_row(&tables.edges, 0, 100, 0, 2, NULL, 0); + + tsk_site_table_add_row(&tables.sites, 10, "A", 1, NULL, 0); + tsk_mutation_table_add_row( + &tables.mutations, 0, 1, TSK_NULL, TSK_UNKNOWN_TIME, "T", 1, NULL, 0); + + ret = tsk_table_collection_sort(&tables, NULL, 0); + CU_ASSERT_EQUAL_FATAL(ret, 0); + ret = tsk_table_collection_build_index(&tables, 0); + CU_ASSERT_EQUAL_FATAL(ret, 0); + + ret = matcher_indexes_alloc(&mi, &tables, NULL, 0); + CU_ASSERT_EQUAL_FATAL(ret, TSI_ERR_MULTIPLE_ROOTS); + matcher_indexes_free(&mi); + tsk_table_collection_free(&tables); +} + static void test_packbits_1(void) { @@ -1370,6 +1550,13 @@ main(int argc, char **argv) test_ancestor_builder_one_bit_encoding }, { "test_ancestor_builder_mmap", test_ancestor_builder_mmap }, { "test_matcher_indexes_errors", test_matcher_indexes_errors }, + { "test_matcher_indexes_edge_node_out_of_bounds", + test_matcher_indexes_edge_node_out_of_bounds }, + { "test_matcher_indexes_mutation_node_out_of_bounds", + test_matcher_indexes_mutation_node_out_of_bounds }, + { "test_matcher_indexes_mutation_site_out_of_bounds", + test_matcher_indexes_mutation_site_out_of_bounds }, + { "test_matcher_indexes_multiple_roots", test_matcher_indexes_multiple_roots }, { "test_packbits_1", test_packbits_1 }, { "test_packbits_2", test_packbits_2 }, From bd982d3a17bbd125b93a2c9c0a02084ecd373c07 Mon Sep 17 00:00:00 2001 From: Jerome Kelleher Date: Mon, 30 Mar 2026 16:23:21 +0100 Subject: [PATCH 05/11] Remove false-positive multiple-roots check from topology validation The TSI_ERR_MULTIPLE_ROOTS check incorrectly rejected valid tree sequences where node 0 has multiple non-overlapping edges to different children (different roots in different trees). The union-based interval tracking created false positives. This is a logic invariant caught by debug asserts, not a memory-safety issue, so remove rather than fix. --- lib/ancestor_matcher.c | 26 -------------------------- lib/err.c | 3 --- lib/err.h | 1 - lib/tests/tests.c | 36 ------------------------------------ 4 files changed, 66 deletions(-) diff --git a/lib/ancestor_matcher.c b/lib/ancestor_matcher.c index 66467dd6..be4a2949 100644 --- a/lib/ancestor_matcher.c +++ b/lib/ancestor_matcher.c @@ -222,38 +222,12 @@ matcher_indexes_validate_tables( const tsk_id_t *restrict edges_child = tables->edges.child; const tsk_id_t *restrict mutations_node = tables->mutations.node; const tsk_id_t *restrict mutations_site = tables->mutations.site; - const double *restrict edges_left = tables->edges.left; - const double *restrict edges_right = tables->edges.right; - /* Track whether node 0 already has an edge as parent. We check for - * overlapping intervals among edges with parent==0, which would mean - * multiple roots in the same tree. We only need to track one interval - * at a time because a single overlapping pair is enough to flag. */ - double root_edge_left = -1; - double root_edge_right = -1; for (j = 0; j < num_edges; j++) { if (edges_parent[j] < 0 || edges_parent[j] >= (tsk_id_t) num_nodes || edges_child[j] < 0 || edges_child[j] >= (tsk_id_t) num_nodes) { return TSI_ERR_BAD_EDGE_NODE; } - if (edges_parent[j] == 0) { - if (root_edge_left < 0) { - root_edge_left = edges_left[j]; - root_edge_right = edges_right[j]; - } else { - /* Check if this edge overlaps with a previous root edge */ - if (edges_left[j] < root_edge_right && edges_right[j] > root_edge_left) { - return TSI_ERR_MULTIPLE_ROOTS; - } - /* Extend the tracked interval */ - if (edges_left[j] < root_edge_left) { - root_edge_left = edges_left[j]; - } - if (edges_right[j] > root_edge_right) { - root_edge_right = edges_right[j]; - } - } - } } for (j = 0; j < num_mutations; j++) { if (mutations_node[j] < 0 || mutations_node[j] >= (tsk_id_t) num_nodes) { diff --git a/lib/err.c b/lib/err.c index b9d5dc2b..7cb47118 100644 --- a/lib/err.c +++ b/lib/err.c @@ -119,9 +119,6 @@ tsi_strerror(int err) case TSI_ERR_BAD_EDGE_NODE: ret = "Bad edge: parent or child node out of bounds"; break; - case TSI_ERR_MULTIPLE_ROOTS: - ret = "Node 0 must have at most one child per tree interval"; - break; } return ret; } diff --git a/lib/err.h b/lib/err.h index 6ba90f51..5f396c2a 100644 --- a/lib/err.h +++ b/lib/err.h @@ -30,7 +30,6 @@ #define TSI_ERR_BAD_ANCESTRAL_STATE -26 #define TSI_ERR_MULTIPLE_MUTATIONS_AT_SITE -27 #define TSI_ERR_BAD_EDGE_NODE -28 -#define TSI_ERR_MULTIPLE_ROOTS -29 // clang-format on #ifdef __GNUC__ diff --git a/lib/tests/tests.c b/lib/tests/tests.c index 40e15b9f..4443b117 100644 --- a/lib/tests/tests.c +++ b/lib/tests/tests.c @@ -586,41 +586,6 @@ test_matcher_indexes_mutation_site_out_of_bounds(void) tsk_table_collection_free(&tables); } -static void -test_matcher_indexes_multiple_roots(void) -{ - int ret; - matcher_indexes_t mi; - tsk_table_collection_t tables; - - memset(&mi, 0, sizeof(mi)); - ret = tsk_table_collection_init(&tables, 0); - CU_ASSERT_EQUAL_FATAL(ret, 0); - tables.sequence_length = 100; - - tsk_node_table_add_row(&tables.nodes, 0, 2.0, TSK_NULL, TSK_NULL, NULL, 0); - tsk_node_table_add_row(&tables.nodes, 0, 1.0, TSK_NULL, TSK_NULL, NULL, 0); - tsk_node_table_add_row(&tables.nodes, 0, 0.5, TSK_NULL, TSK_NULL, NULL, 0); - - /* Two edges with parent=0 overlapping over [0,100): multiple roots */ - tsk_edge_table_add_row(&tables.edges, 0, 100, 0, 1, NULL, 0); - tsk_edge_table_add_row(&tables.edges, 0, 100, 0, 2, NULL, 0); - - tsk_site_table_add_row(&tables.sites, 10, "A", 1, NULL, 0); - tsk_mutation_table_add_row( - &tables.mutations, 0, 1, TSK_NULL, TSK_UNKNOWN_TIME, "T", 1, NULL, 0); - - ret = tsk_table_collection_sort(&tables, NULL, 0); - CU_ASSERT_EQUAL_FATAL(ret, 0); - ret = tsk_table_collection_build_index(&tables, 0); - CU_ASSERT_EQUAL_FATAL(ret, 0); - - ret = matcher_indexes_alloc(&mi, &tables, NULL, 0); - CU_ASSERT_EQUAL_FATAL(ret, TSI_ERR_MULTIPLE_ROOTS); - matcher_indexes_free(&mi); - tsk_table_collection_free(&tables); -} - static void test_packbits_1(void) { @@ -1556,7 +1521,6 @@ main(int argc, char **argv) test_matcher_indexes_mutation_node_out_of_bounds }, { "test_matcher_indexes_mutation_site_out_of_bounds", test_matcher_indexes_mutation_site_out_of_bounds }, - { "test_matcher_indexes_multiple_roots", test_matcher_indexes_multiple_roots }, { "test_packbits_1", test_packbits_1 }, { "test_packbits_2", test_packbits_2 }, From 093452390968bfdf367062aefb0828952beb9eb5 Mon Sep 17 00:00:00 2001 From: Jerome Kelleher Date: Mon, 30 Mar 2026 16:28:03 +0100 Subject: [PATCH 06/11] Remind claude to rebuild the C module. --- CLAUDE.md | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/CLAUDE.md b/CLAUDE.md index fe1d9a59..b050da5d 100644 --- a/CLAUDE.md +++ b/CLAUDE.md @@ -77,7 +77,9 @@ The public API is in `tsinfer/__init__.py`, exposing three main functions from ` Source in `lib/`. Three main classes exposed to Python: - `AncestorBuilder` — builds inferred ancestors from genotype data - `AncestorMatcher` — Li & Stephens HMM matching algorithm -- `TreeSequenceBuilder` — constructs tree sequences incrementally + +When changes are made to the C library, ensure that the ``_tskit`` module is rebuilt +before running Python tests. Vendored dependencies in `lib/subprojects/`: tskit C library and kastore. From 2f21031f0b9ee6d46223dff116bc4141c52fdd6b Mon Sep 17 00:00:00 2001 From: Jerome Kelleher Date: Mon, 30 Mar 2026 16:42:31 +0100 Subject: [PATCH 07/11] Validate that node 0 is root and add test cases for root topology Add TSI_ERR_NODE_0_NOT_ROOT check in matcher_indexes_validate_tables: verifies at least one edge has parent==0. Without this, parent-chain traversals reach NULL_NODE and segfault on array index -1. Add C tests: node-0-not-root error case and a valid multi-root tree sequence where node 0 has non-overlapping edges to different children. Add Python test for the error case. Fix existing Python tests that created MatcherIndexes without a vestigial root. --- lib/ancestor_matcher.c | 7 ++++ lib/err.c | 3 ++ lib/err.h | 1 + lib/tests/tests.c | 91 ++++++++++++++++++++++++++++++++++++++++++ tests/test_python_c.py | 11 +++++ 5 files changed, 113 insertions(+) diff --git a/lib/ancestor_matcher.c b/lib/ancestor_matcher.c index be4a2949..7c695873 100644 --- a/lib/ancestor_matcher.c +++ b/lib/ancestor_matcher.c @@ -222,12 +222,19 @@ matcher_indexes_validate_tables( const tsk_id_t *restrict edges_child = tables->edges.child; const tsk_id_t *restrict mutations_node = tables->mutations.node; const tsk_id_t *restrict mutations_site = tables->mutations.site; + bool node_0_is_parent = false; for (j = 0; j < num_edges; j++) { if (edges_parent[j] < 0 || edges_parent[j] >= (tsk_id_t) num_nodes || edges_child[j] < 0 || edges_child[j] >= (tsk_id_t) num_nodes) { return TSI_ERR_BAD_EDGE_NODE; } + if (edges_parent[j] == 0) { + node_0_is_parent = true; + } + } + if (num_edges > 0 && !node_0_is_parent) { + return TSI_ERR_NODE_0_NOT_ROOT; } for (j = 0; j < num_mutations; j++) { if (mutations_node[j] < 0 || mutations_node[j] >= (tsk_id_t) num_nodes) { diff --git a/lib/err.c b/lib/err.c index 7cb47118..334a23ad 100644 --- a/lib/err.c +++ b/lib/err.c @@ -119,6 +119,9 @@ tsi_strerror(int err) case TSI_ERR_BAD_EDGE_NODE: ret = "Bad edge: parent or child node out of bounds"; break; + case TSI_ERR_NODE_0_NOT_ROOT: + ret = "Node 0 must be the root: no edges have parent 0"; + break; } return ret; } diff --git a/lib/err.h b/lib/err.h index 5f396c2a..61196485 100644 --- a/lib/err.h +++ b/lib/err.h @@ -30,6 +30,7 @@ #define TSI_ERR_BAD_ANCESTRAL_STATE -26 #define TSI_ERR_MULTIPLE_MUTATIONS_AT_SITE -27 #define TSI_ERR_BAD_EDGE_NODE -28 +#define TSI_ERR_NODE_0_NOT_ROOT -29 // clang-format on #ifdef __GNUC__ diff --git a/lib/tests/tests.c b/lib/tests/tests.c index 4443b117..c0c26142 100644 --- a/lib/tests/tests.c +++ b/lib/tests/tests.c @@ -586,6 +586,42 @@ test_matcher_indexes_mutation_site_out_of_bounds(void) tsk_table_collection_free(&tables); } +static void +test_matcher_indexes_node_0_not_root(void) +{ + int ret; + matcher_indexes_t mi; + tsk_table_collection_t tables; + + memset(&mi, 0, sizeof(mi)); + ret = tsk_table_collection_init(&tables, 0); + CU_ASSERT_EQUAL_FATAL(ret, 0); + tables.sequence_length = 100; + + /* Node 0 is a leaf, node 2 is the root */ + tsk_node_table_add_row(&tables.nodes, 0, 0.5, TSK_NULL, TSK_NULL, NULL, 0); + tsk_node_table_add_row(&tables.nodes, 0, 0.3, TSK_NULL, TSK_NULL, NULL, 0); + tsk_node_table_add_row(&tables.nodes, 0, 2.0, TSK_NULL, TSK_NULL, NULL, 0); + + /* Edges: node 2 is parent of 0 and 1 (node 0 is NOT a parent) */ + tsk_edge_table_add_row(&tables.edges, 0, 100, 2, 0, NULL, 0); + tsk_edge_table_add_row(&tables.edges, 0, 100, 2, 1, NULL, 0); + + tsk_site_table_add_row(&tables.sites, 10, "A", 1, NULL, 0); + tsk_mutation_table_add_row( + &tables.mutations, 0, 0, TSK_NULL, TSK_UNKNOWN_TIME, "T", 1, NULL, 0); + + ret = tsk_table_collection_sort(&tables, NULL, 0); + CU_ASSERT_EQUAL_FATAL(ret, 0); + ret = tsk_table_collection_build_index(&tables, 0); + CU_ASSERT_EQUAL_FATAL(ret, 0); + + ret = matcher_indexes_alloc(&mi, &tables, NULL, 0); + CU_ASSERT_EQUAL_FATAL(ret, TSI_ERR_NODE_0_NOT_ROOT); + matcher_indexes_free(&mi); + tsk_table_collection_free(&tables); +} + static void test_packbits_1(void) { @@ -1417,6 +1453,59 @@ test_matching_root_switch(void) tsk_treeseq_free(&ts); } +/* Test that a valid tree sequence with non-overlapping edges from node 0 + * to different children (different roots in different trees) works correctly. + * This is a regression test for a previous false positive. */ +static void +test_matching_nonzero_root_valid(void) +{ + int ret; + tsk_treeseq_t ts; + tsk_table_collection_t tables; + allele_t h[] = { 0, 0 }; + allele_t match[2]; + tsk_id_t left[4], right[4], parent[4]; + tsk_size_t path_length; + + ret = tsk_table_collection_init(&tables, 0); + CU_ASSERT_EQUAL_FATAL(ret, 0); + tables.sequence_length = 100; + + /* Node 0 = virtual root, node 1 = root of tree 1, node 2 = root of tree 2, + * node 3 = leaf */ + tsk_node_table_add_row(&tables.nodes, 0, 3.0, TSK_NULL, TSK_NULL, NULL, 0); + tsk_node_table_add_row(&tables.nodes, 0, 2.0, TSK_NULL, TSK_NULL, NULL, 0); + tsk_node_table_add_row(&tables.nodes, 0, 2.0, TSK_NULL, TSK_NULL, NULL, 0); + tsk_node_table_add_row(&tables.nodes, 0, 1.0, TSK_NULL, TSK_NULL, NULL, 0); + + /* Node 0 is parent in both halves, but with different children */ + tsk_edge_table_add_row(&tables.edges, 0, 50, 0, 1, NULL, 0); + tsk_edge_table_add_row(&tables.edges, 50, 100, 0, 2, NULL, 0); + /* Node 3 is child of 1 in first half, child of 2 in second half */ + tsk_edge_table_add_row(&tables.edges, 0, 50, 1, 3, NULL, 0); + tsk_edge_table_add_row(&tables.edges, 50, 100, 2, 3, NULL, 0); + + /* One site per tree */ + tsk_site_table_add_row(&tables.sites, 25, "A", 1, NULL, 0); + tsk_site_table_add_row(&tables.sites, 75, "A", 1, NULL, 0); + tsk_mutation_table_add_row( + &tables.mutations, 0, 3, TSK_NULL, TSK_UNKNOWN_TIME, "T", 1, NULL, 0); + tsk_mutation_table_add_row( + &tables.mutations, 1, 3, TSK_NULL, TSK_UNKNOWN_TIME, "T", 1, NULL, 0); + + ret = tsk_table_collection_sort(&tables, NULL, 0); + CU_ASSERT_EQUAL_FATAL(ret, 0); + ret = tsk_treeseq_init(&ts, &tables, TSK_TS_INIT_BUILD_INDEXES); + CU_ASSERT_EQUAL_FATAL(ret, 0); + tsk_table_collection_free(&tables); + + run_match(&ts, 1e-8, 1e-20, h, match, &path_length, left, right, parent); + /* With different roots in each tree, the path may have 1 or 2 segments */ + CU_ASSERT_TRUE_FATAL(path_length >= 1 && path_length <= 2); + + tsk_treeseq_free(&ts); +} + static void test_strerror(void) { @@ -1521,6 +1610,7 @@ main(int argc, char **argv) test_matcher_indexes_mutation_node_out_of_bounds }, { "test_matcher_indexes_mutation_site_out_of_bounds", test_matcher_indexes_mutation_site_out_of_bounds }, + { "test_matcher_indexes_node_0_not_root", test_matcher_indexes_node_0_not_root }, { "test_packbits_1", test_packbits_1 }, { "test_packbits_2", test_packbits_2 }, @@ -1550,6 +1640,7 @@ main(int argc, char **argv) { "test_matching_impossible_extreme_mu", test_matching_impossible_extreme_mu }, { "test_matching_impossible_zero_recomb", test_matching_impossible_zero_recomb }, { "test_matching_root_switch", test_matching_root_switch }, + { "test_matching_nonzero_root_valid", test_matching_nonzero_root_valid }, { "test_strerror", test_strerror }, diff --git a/tests/test_python_c.py b/tests/test_python_c.py index a04e3409..b4237334 100644 --- a/tests/test_python_c.py +++ b/tests/test_python_c.py @@ -299,11 +299,20 @@ def test_get_traceback_bad_site(self): class TestMatcherIndexes: def test_single_tree(self): ts = tskit.Tree.generate_balanced(4).tree_sequence + ts = matching.add_vestigial_root(ts) tables = ts.dump_tables() ll_tables = _tsinfer.LightweightTableCollection(tables.sequence_length) ll_tables.fromdict(tables.asdict()) _ = _tsinfer.MatcherIndexes(ll_tables) + def test_node_0_not_root(self): + ts = tskit.Tree.generate_balanced(4).tree_sequence + tables = ts.dump_tables() + ll_tables = _tsinfer.LightweightTableCollection(tables.sequence_length) + ll_tables.fromdict(tables.asdict()) + with pytest.raises(_tsinfer.LibraryError, match="Node 0 must be the root"): + _tsinfer.MatcherIndexes(ll_tables) + def test_num_alleles(self): ts, mi, _ = make_matcher_indexes_and_matcher() num_alleles = np.array([2] * ts.num_sites, dtype=np.uint64) @@ -314,6 +323,7 @@ def test_num_alleles(self): def test_print_state(self, tmpdir): ts = tskit.Tree.generate_balanced(4).tree_sequence + ts = matching.add_vestigial_root(ts) tables = ts.dump_tables() ll_tables = _tsinfer.LightweightTableCollection(tables.sequence_length) ll_tables.fromdict(tables.asdict()) @@ -331,6 +341,7 @@ def test_print_state(self, tmpdir): def test_print_state_bad_file(self): ts = tskit.Tree.generate_balanced(4).tree_sequence + ts = matching.add_vestigial_root(ts) tables = ts.dump_tables() ll_tables = _tsinfer.LightweightTableCollection(tables.sequence_length) ll_tables.fromdict(tables.asdict()) From 27122b4f1fd4dd558eb91be9edd131bd0761955a Mon Sep 17 00:00:00 2001 From: Jerome Kelleher Date: Mon, 30 Mar 2026 16:59:31 +0100 Subject: [PATCH 08/11] Add validation that node 0 is never a child edge Check that no edge has child==0, which would mean node 0 is not always a root. Parent-chain traversals assume node 0 is always parentless; violating this causes incorrect results or segfaults. The full invariant (node 0 has exactly one child per tree) cannot be verified without building trees, so it remains guarded by the existing debug assert in the forward match. --- lib/ancestor_matcher.c | 3 +++ lib/tests/tests.c | 42 ++++++++++++++++++++++++++++++++++++++++++ 2 files changed, 45 insertions(+) diff --git a/lib/ancestor_matcher.c b/lib/ancestor_matcher.c index 7c695873..297aaa26 100644 --- a/lib/ancestor_matcher.c +++ b/lib/ancestor_matcher.c @@ -229,6 +229,9 @@ matcher_indexes_validate_tables( || edges_child[j] < 0 || edges_child[j] >= (tsk_id_t) num_nodes) { return TSI_ERR_BAD_EDGE_NODE; } + if (edges_child[j] == 0) { + return TSI_ERR_NODE_0_NOT_ROOT; + } if (edges_parent[j] == 0) { node_0_is_parent = true; } diff --git a/lib/tests/tests.c b/lib/tests/tests.c index c0c26142..15ce1821 100644 --- a/lib/tests/tests.c +++ b/lib/tests/tests.c @@ -622,6 +622,47 @@ test_matcher_indexes_node_0_not_root(void) tsk_table_collection_free(&tables); } +/* Node 0 appears as a child in an edge, even though it is also a parent. + * This means node 0 is not always a root, which violates the invariant. */ +static void +test_matcher_indexes_node_0_is_child(void) +{ + int ret; + matcher_indexes_t mi; + tsk_table_collection_t tables; + + memset(&mi, 0, sizeof(mi)); + ret = tsk_table_collection_init(&tables, 0); + CU_ASSERT_EQUAL_FATAL(ret, 0); + tables.sequence_length = 100; + + /* Node 0 has highest time but we'll make it a child via overwrite */ + tsk_node_table_add_row(&tables.nodes, 0, 3.0, TSK_NULL, TSK_NULL, NULL, 0); + tsk_node_table_add_row(&tables.nodes, 0, 2.0, TSK_NULL, TSK_NULL, NULL, 0); + tsk_node_table_add_row(&tables.nodes, 0, 1.0, TSK_NULL, TSK_NULL, NULL, 0); + + /* Valid edges: 0->1, 1->2 */ + tsk_edge_table_add_row(&tables.edges, 0, 100, 0, 1, NULL, 0); + tsk_edge_table_add_row(&tables.edges, 0, 100, 1, 2, NULL, 0); + + tsk_site_table_add_row(&tables.sites, 10, "A", 1, NULL, 0); + tsk_mutation_table_add_row( + &tables.mutations, 0, 2, TSK_NULL, TSK_UNKNOWN_TIME, "T", 1, NULL, 0); + + ret = tsk_table_collection_sort(&tables, NULL, 0); + CU_ASSERT_EQUAL_FATAL(ret, 0); + ret = tsk_table_collection_build_index(&tables, 0); + CU_ASSERT_EQUAL_FATAL(ret, 0); + + /* Overwrite child of first edge to make node 0 a child */ + tables.edges.child[0] = 0; + + ret = matcher_indexes_alloc(&mi, &tables, NULL, 0); + CU_ASSERT_EQUAL_FATAL(ret, TSI_ERR_NODE_0_NOT_ROOT); + matcher_indexes_free(&mi); + tsk_table_collection_free(&tables); +} + static void test_packbits_1(void) { @@ -1611,6 +1652,7 @@ main(int argc, char **argv) { "test_matcher_indexes_mutation_site_out_of_bounds", test_matcher_indexes_mutation_site_out_of_bounds }, { "test_matcher_indexes_node_0_not_root", test_matcher_indexes_node_0_not_root }, + { "test_matcher_indexes_node_0_is_child", test_matcher_indexes_node_0_is_child }, { "test_packbits_1", test_packbits_1 }, { "test_packbits_2", test_packbits_2 }, From 4d5c1e439754b6a7128da879df82c6821cb0daba Mon Sep 17 00:00:00 2001 From: Jerome Kelleher Date: Tue, 31 Mar 2026 10:03:27 +0100 Subject: [PATCH 09/11] Add banner. --- docs/_config.yml | 1 + docs/index.md | 5 ----- 2 files changed, 1 insertion(+), 5 deletions(-) diff --git a/docs/_config.yml b/docs/_config.yml index f4775224..fa3c1ac3 100644 --- a/docs/_config.yml +++ b/docs/_config.yml @@ -38,6 +38,7 @@ sphinx: config: html_theme: sphinx_book_theme html_theme_options: + announcement: "⚠ This documentation is under active development. The API is not yet stable and many elements are out of date." pygments_dark_style: monokai navigation_with_keys: false logo: diff --git a/docs/index.md b/docs/index.md index dacb3052..e44db54a 100644 --- a/docs/index.md +++ b/docs/index.md @@ -1,8 +1,3 @@ -```{warning} -This documentation is under active development and may be incomplete or -inaccurate. The software API is not yet stable. -``` - # Welcome to tsinfer's documentation! This is the documentation for {program}`tsinfer`, a method for inferring correlated From b9b1b9602eb8008f9d69dcc0a60704d9971388ba Mon Sep 17 00:00:00 2001 From: Jerome Kelleher Date: Tue, 31 Mar 2026 10:09:19 +0100 Subject: [PATCH 10/11] Reinstate large_scale.md as static legacy reference Restore the large scale inference page with jupytext headers stripped, {meth}/{class} references converted to plain text, and a legacy warning banner added. Builds cleanly without execution. --- docs/_toc.yml | 1 + docs/large_scale.md | 196 ++++++++++++++++++++++++++++++++++++++++++++ 2 files changed, 197 insertions(+) create mode 100644 docs/large_scale.md diff --git a/docs/_toc.yml b/docs/_toc.yml index 24aa5f1a..d31ae890 100644 --- a/docs/_toc.yml +++ b/docs/_toc.yml @@ -16,6 +16,7 @@ parts: - caption: Inference chapters: - file: inference + - file: large_scale - caption: Interfaces chapters: - file: cli diff --git a/docs/large_scale.md b/docs/large_scale.md new file mode 100644 index 00000000..ee58781d --- /dev/null +++ b/docs/large_scale.md @@ -0,0 +1,196 @@ +(sec_large_scale)= + +# Large scale inference (legacy) + +```{warning} +This page documents the **old Python batch matching API** which is no +longer the recommended workflow. It is retained as a reference while the +content is ported to the new TOML config + CLI interface described in the +{ref}`quickstart `. +``` + +Generally, for up to a few thousand samples a single multi-core machine +can infer a tree sequence in a few days, hours, or even minutes. +However, _tsinfer_ has been successfully used with datasets up to half a million +samples, where ancestor and sample matching can take several CPU-years. +At this scale inference must be scaled across many machines. +_Tsinfer_ provides specific APIs to enable this. +Here we detail considerations and tips for each step of the +inference process to help you scale up your analysis. A snakemake pipeline +which implements this parallelisation scheme is available as +[tsinfer-snakemake](https://github.com/benjeffery/tsinfer-snakemake). + +(sec_large_scale_ancestor_generation)= + +## Data preparation + +For large scale inference the data must be in [VCF Zarr](https://github.com/sgkit-dev/vcf-zarr-spec) +format, read by the `VariantData` class. [Bio2zarr](https://github.com/sgkit-dev/bio2zarr) +is recommended for conversion from VCF, and [sgkit](https://github.com/sgkit-dev/sgkit) can then +be used to perform initial filtering. + +:::{todo} +An upcoming tutorial will detail conversion from VCF to a VCF Zarr suitable for tsinfer. +::: + + +## Ancestor generation + +Ancestor generation is generally the fastest step in inference. It is not yet +parallelised out-of-core in tsinfer and must be performed on a single machine. +However it scales well on machines with +many cores and hyperthreading via the `num_threads` argument to +`generate_ancestors`. The limiting factor is often that the +entire genotype array for the contig being inferred needs to fit in RAM. +This is the high-water mark for memory usage in tsinfer. + +If your data consists of only biallelic sites, with no missingness, +the `genotype_encoding` argument can be set to +`GenotypeEncoding.ONE_BIT` which reduces the memory footprint of +the genotype array by a factor of 8, such that the RAM needed is roughly +`num_sites * num_samples * ploidy / 8 bytes`. This memory optimisation +results in a surprisingly small increase in runtime. + +## Ancestor matching + +Ancestor matching is one of the more time consuming steps of inference. It +proceeds in groups, progressively growing the tree sequence with younger +ancestors. At each stage the parallelism is limited to the number of ancestors +whose possible inheritors are already matched, as all possible inheritors +of a sample must be matched in an earlier group. For a typical human data set +the number of samples per group varies from single digits up to approximately +the number of samples. +The plot below shows the number of ancestors matched in each group for a typical +human data set, earlier groups are older ancestors: + +```{figure} _static/ancestor_grouping.png +:width: 80% +``` + +There are five tsinfer API methods that can be used to parallelise ancestor +matching. + +The five methods are: + +1. `match_ancestors_batch_init` +2. `match_ancestors_batch_groups` +3. `match_ancestors_batch_group_partition` +4. `match_ancestors_batch_group_finalise` +5. `match_ancestors_batch_finalise` + +Initially `match_ancestors_batch_init` should be called to +set up the batch matching and to determine the groupings of ancestors. +This method writes a file `metadata.json` to the `work_dir` that contains +a JSON encoded dictionary with configuration for later steps, and a key +`ancestor_grouping` which is a list of dictionaries, each containing the +list of ancestors in that group (key:`ancestors`) and a proposed partioning of +those ancestors into sets that can be matched in parallel (key:`partitions`). +The dictionary is also returned by the method. +The partitioning is controlled by the `min_work_per_job` and `max_num_partitions` +arguments. For each group, ancestors are placed in a partition until the sum of their +lengths exceeds `min_work_per_job`, when a new partition is started. However, the +number of partitions is not allowed to exceed `max_num_partitions`. It is suggested +to set `max_num_partitions` to around 3-4x the number of worker nodes available, +and `min_work_per_job` to around 2,000,000 for a typical human data set. + +Groups vs partitions is a point of common confusion. Note that groups of ancestors +are matched serially, and each group is split into partitions that can be +matched in parallel. + +Each group is matched in turn, either by calling `match_ancestors_batch_groups` +to match without partitioning, or by calling `match_ancestors_batch_group_partition` +many times in parallel followed by a single call to `match_ancestors_batch_group_finalise`. +Each call to `match_ancestors_batch_groups` or `match_ancestors_batch_group_finalise` +outputs the tree sequence to `work_dir`, which is then used by the next group. The length of +the `ancestor_grouping` in the metadata dictionary determines the group numbers that these methods +will need to be called for, and the length of the `partitions` list in each group determines +the number of calls to `match_ancestors_batch_group_partition` that are needed (if any). + +`match_ancestors_batch_groups` matches groups, without partitioning, from +`group_index_start` (inclusively) to `group_index_end` (exclusively). Combining +many groups into one call reduces the overhead from job submission and start +up times, but note on job failure the process can only be resumed from the +last `group_index_end`. + +To match a single group in parallel, call `match_ancestors_batch_group_partition` +once for each partition listed in the `ancestor_grouping[group_index]['partitions']` list, +incrementing `partition_index`. This will match the ancestors, placing the match data in +the `working_dir`. Once all are complete a single call to +`match_ancestors_batch_group_finalise` will then insert the matches and +output the tree sequence to `work_dir`. + +Each call to `match_ancestors_batch_groups` and `match_ancestors_batch_group_finalise` results in a tree sequence being written to `work_dir`. +These tree sequences are essentially checkpoints from with the batch matching workflow +can be resumed on job failure. + +Finally after the final group, call `match_ancestors_batch_finalise` to +combine the groups into a single tree sequence. + +The partitioning in `metadata.json` does not have to be used for every group. As early groups are +not matching to a large tree sequence it is often faster to not partition the first half of the +groups, depending on job set up and queueing time on your cluster. + +Calls to `match_ancestors_batch_group_partition` will only use a single core, but +`match_ancestors_batch_groups` will use as many cores as `num_threads` is set to. +Therefore this value and cluster resources requested should be scaled with the number of ancestors, +which can be read from the metadata dictionary. + +As an example of how the API methods can be used together, suppose the metadata dictionary +created by `match_ancestors_batch_init` contains the following: + +```python +{ + "ancestor_grouping": [ + {"ancestors": [0, ... 9], "partitions": None}, + { + "ancestors": [10, ... 15], + "partitions": [[10, 11, 12], [13, 14, 15]] + }, + {"ancestors": [16, ... 19], "partitions": None}, + {"ancestors": [20, ... 25], "partitions": None}, + {"ancestors": [26, ... 30], "partitions": None}, + { + "ancestors": [31, ... 41], + "partitions": [[31, 32, 33, 34, 35, 36], [37, 38, 39, 40, 41]] + }, + {"ancestors": [42, ... 45], "partitions": None}, + {"ancestors": [46, ... 50], "partitions": None}, + { + "ancestors": [51, ... 65], + "partitions": [ + [51, 52, 53, 54], + [55, 56, 57, 58], + [59, 60, 61, 62, 63, 64, 65] + ] + }, + ] +} +``` +Then the flow could look like the following diagram: (calls on the same horizontal line can be +done in parallel, note that method names are shortened): + +```{figure} _static/example_flow.svg +:width: 80% +``` + +Note that groups 1, 5 and 8 can be partitioned, but only groups 5 and 8 are actually partitioned in this example, as stated above partitioning for groups is optional. Groups 0-4 are matched in one call, groups 6 and 7 are matched in two calls, but +could have been matched in one. By splitting 6 and 7 the flow makes an additional resume point in the case of job failure at the cost of job start up and queueing time. + + +## Sample matching + +Sample matching is far simpler than ancestor matching as it is essentially the same as a single group +of ancestors. There are three API methods that work together to enable distributed sample matching. + +1. `match_samples_batch_init` +2. `match_samples_batch_partition` +3. `match_samples_batch_finalise` + +`match_samples_batch_init` should be called to set up the batch matching and to determine the +groupings of samples. Similar to `match_ancestors_batch_init` it has a `min_work_per_job` argument to control the level of parallelism. The method writes a file +`metadata.json` to the directory `work_dir` that contains a JSON encoded dictionary with +configuration for later steps. This is also returned by the call. The `num_partitions` key in +this dictionary is the number of times `match_samples_batch_partition` will need +to be called, with each partition index as the `partition_index` argument. These calls can happen +in parallel and write match data to the `work_dir` which is then used by +`match_samples_batch_finalise` to output the tree sequence. \ No newline at end of file From e4919f16d289383f6ed0eb0399e731a8dcf0d2a4 Mon Sep 17 00:00:00 2001 From: Jerome Kelleher Date: Tue, 31 Mar 2026 10:13:49 +0100 Subject: [PATCH 11/11] Improve C coverage, remind claude to do this. --- CLAUDE.md | 3 +++ lib/tests/tests.c | 39 +++++++++++++++++++++++++++++++++++++++ 2 files changed, 42 insertions(+) diff --git a/CLAUDE.md b/CLAUDE.md index b050da5d..6075a631 100644 --- a/CLAUDE.md +++ b/CLAUDE.md @@ -47,6 +47,9 @@ valgrind --leak-check=full --error-exitcode=1 ./build/tests Tests are in `lib/tests/tests.c` using the CUnit framework. The build uses `-Wall -Wextra -Werror -Wpedantic` and other strict warnings. +Ensure that all new C code is covered by tests in the C test suite by running +tests with coverage. + ## Architecture **tsinfer** infers tree sequences from genetic variation data stored in VCZ (Variant Call Zarr) format. diff --git a/lib/tests/tests.c b/lib/tests/tests.c index 15ce1821..cda33456 100644 --- a/lib/tests/tests.c +++ b/lib/tests/tests.c @@ -663,6 +663,43 @@ test_matcher_indexes_node_0_is_child(void) tsk_table_collection_free(&tables); } +/* Node 0 is disconnected: not a child, not a parent. Edges exist but none + * involve node 0. This hits the "no edges have parent==0" path. */ +static void +test_matcher_indexes_node_0_disconnected(void) +{ + int ret; + matcher_indexes_t mi; + tsk_table_collection_t tables; + + memset(&mi, 0, sizeof(mi)); + ret = tsk_table_collection_init(&tables, 0); + CU_ASSERT_EQUAL_FATAL(ret, 0); + tables.sequence_length = 100; + + /* Node 0 exists but has no edges */ + tsk_node_table_add_row(&tables.nodes, 0, 3.0, TSK_NULL, TSK_NULL, NULL, 0); + tsk_node_table_add_row(&tables.nodes, 0, 2.0, TSK_NULL, TSK_NULL, NULL, 0); + tsk_node_table_add_row(&tables.nodes, 0, 1.0, TSK_NULL, TSK_NULL, NULL, 0); + + /* Only edge is between nodes 1 and 2, node 0 not involved */ + tsk_edge_table_add_row(&tables.edges, 0, 100, 1, 2, NULL, 0); + + tsk_site_table_add_row(&tables.sites, 10, "A", 1, NULL, 0); + tsk_mutation_table_add_row( + &tables.mutations, 0, 2, TSK_NULL, TSK_UNKNOWN_TIME, "T", 1, NULL, 0); + + ret = tsk_table_collection_sort(&tables, NULL, 0); + CU_ASSERT_EQUAL_FATAL(ret, 0); + ret = tsk_table_collection_build_index(&tables, 0); + CU_ASSERT_EQUAL_FATAL(ret, 0); + + ret = matcher_indexes_alloc(&mi, &tables, NULL, 0); + CU_ASSERT_EQUAL_FATAL(ret, TSI_ERR_NODE_0_NOT_ROOT); + matcher_indexes_free(&mi); + tsk_table_collection_free(&tables); +} + static void test_packbits_1(void) { @@ -1653,6 +1690,8 @@ main(int argc, char **argv) test_matcher_indexes_mutation_site_out_of_bounds }, { "test_matcher_indexes_node_0_not_root", test_matcher_indexes_node_0_not_root }, { "test_matcher_indexes_node_0_is_child", test_matcher_indexes_node_0_is_child }, + { "test_matcher_indexes_node_0_disconnected", + test_matcher_indexes_node_0_disconnected }, { "test_packbits_1", test_packbits_1 }, { "test_packbits_2", test_packbits_2 },