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 diff --git a/CLAUDE.md b/CLAUDE.md index fe1d9a59..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. @@ -77,7 +80,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. diff --git a/docs/_config.yml b/docs/_config.yml index a99feb3c..fa3c1ac3 100644 --- a/docs/_config.yml +++ b/docs/_config.yml @@ -32,12 +32,13 @@ sphinx: - sphinx.ext.viewcode - sphinx.ext.intersphinx - sphinx_issues - - sphinxarg.ext + - sphinx_click - IPython.sphinxext.ipython_console_highlighting 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/_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/_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 00000000..67ff20b8 Binary files /dev/null and b/docs/_static/example_data.vcf.gz differ diff --git a/docs/_toc.yml b/docs/_toc.yml index 6ba9d46d..d31ae890 100644 --- a/docs/_toc.yml +++ b/docs/_toc.yml @@ -10,6 +10,8 @@ parts: - file: installation - caption: Usage chapters: + - file: quickstart + - file: config - file: usage - caption: Inference chapters: @@ -17,11 +19,7 @@ parts: - 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/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 index 17d3ee69..ee58781d 100644 --- a/docs/large_scale.md +++ b/docs/large_scale.md @@ -1,22 +1,13 @@ ---- -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 +# 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. @@ -34,7 +25,7 @@ which implements this parallelisation scheme is available as ## 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) +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. @@ -49,13 +40,13 @@ 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 +`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 +`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. @@ -81,13 +72,13 @@ 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` +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 {meth}`match_ancestors_batch_init` should be called to +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 @@ -106,46 +97,46 @@ Groups vs partitions is a point of common confusion. Note that groups of ancesto 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` +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 {meth}`match_ancestors_batch_group_partition` that are needed (if any). +the number of calls to `match_ancestors_batch_group_partition` that are needed (if any). -{meth}`match_ancestors_batch_groups` matches groups, without partitioning, from +`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` +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 -{meth}`match_ancestors_batch_group_finalise` will then insert the matches and +`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`. +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 {meth}`match_ancestors_batch_finalise` to +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 {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. +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 {meth}`match_ancestors_batch_init` contains the following: +created by `match_ancestors_batch_init` contains the following: ```python { @@ -191,15 +182,15 @@ could have been matched in one. By splitting 6 and 7 the flow makes an additiona 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` +1. `match_samples_batch_init` +2. `match_samples_batch_partition` +3. `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 +`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 {meth}`match_samples_batch_partition` will need +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 -{meth}`match_samples_batch_finalise` to output the tree sequence. \ No newline at end of file +`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 new file mode 100644 index 00000000..4cddd38e --- /dev/null +++ b/docs/quickstart.md @@ -0,0 +1,131 @@ +(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/) + + +## 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 +``` + +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. + + +## 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 +[[source]] +name = "mydata" +path = "mydata.vcz" + +[ancestral_state] +path = "mydata.vcz" +field = "variant_AA" + +[[ancestors]] +name = "ancestors" +path = "ancestors.vcz" +sources = ["mydata"] + +[match] +output = "output.trees" + +[match.sources.ancestors] +node_flags = 0 +create_individuals = false + +[match.sources.mydata] +``` + +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 `. + + +## Running the pipeline + +Run all steps in one command: + +```bash +tsinfer run config.toml --threads 4 -v +``` + +Or run steps individually: + +```bash +tsinfer infer-ancestors config.toml --threads 4 -v +tsinfer match config.toml --threads 4 -v +``` + +Validate a config before running: + +```bash +tsinfer config check config.toml +``` + +See the {ref}`CLI reference ` for all commands and options. + + +## Inspecting the result + +The output is a standard [tskit](https://tskit.dev/) tree sequence: + +```python +import tskit + +ts = tskit.load("output.trees") +print(f"{ts.num_trees} trees, {ts.num_samples} samples, {ts.num_sites} sites") +ts.draw_svg(size=(600, 300), y_axis=True) +``` + +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} +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. +::: + + +(sec_quickstart_inference_sites)= + +## Inference sites + +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). +These include: + +- **Fixed sites** — no variation between samples +- **Singletons** — only one genome carries the derived allele +- **Unknown ancestral state** — ancestral allele does not match any allele +- **Multiallelic sites** — more than two alleles 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 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; + 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_child[j] == 0) { + return TSI_ERR_NODE_0_NOT_ROOT; + } + 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) { + 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 +264,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..334a23ad 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_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 08eeca41..61196485 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_NODE_0_NOT_ROOT -29 // clang-format on #ifdef __GNUC__ diff --git a/lib/tests/tests.c b/lib/tests/tests.c index 48340998..cda33456 100644 --- a/lib/tests/tests.c +++ b/lib/tests/tests.c @@ -441,6 +441,265 @@ 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_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); +} + +/* 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); +} + +/* 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) { @@ -1272,6 +1531,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) { @@ -1370,6 +1682,16 @@ 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_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 }, @@ -1399,6 +1721,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/pyproject.toml b/pyproject.toml index cee8bf19..d99551ac 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -90,7 +90,7 @@ docs = [ "jupyter-book<2", "msprime", "pyfaidx", - "sphinx-argparse", + "sphinx-click", "sphinx-book-theme", "sphinx-issues", ] 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()) diff --git a/uv.lock b/uv.lock index 1f1c19a2..08ecdede 100644 --- a/uv.lock +++ b/uv.lock @@ -2449,29 +2449,30 @@ wheels = [ ] [[package]] -name = "sphinx-argparse" -version = "0.5.2" +name = "sphinx-book-theme" +version = "1.1.4" source = { registry = "https://pypi.org/simple" } dependencies = [ - { name = "docutils" }, + { name = "pydata-sphinx-theme" }, { name = "sphinx" }, ] -sdist = { url = "https://files.pythonhosted.org/packages/3b/21/a8c64e6633652111e6e4f89703182a53cbc3ed67233523e47472101358b6/sphinx_argparse-0.5.2.tar.gz", hash = "sha256:e5352f8fa894b6fb6fda0498ba28a9f8d435971ef4bbc1a6c9c6414e7644f032", size = 27838, upload-time = "2024-07-17T12:08:08.219Z" } +sdist = { url = "https://files.pythonhosted.org/packages/45/19/d002ed96bdc7738c15847c730e1e88282d738263deac705d5713b4d8fa94/sphinx_book_theme-1.1.4.tar.gz", hash = "sha256:73efe28af871d0a89bd05856d300e61edce0d5b2fbb7984e84454be0fedfe9ed", size = 439188, upload-time = "2025-02-20T16:32:32.581Z" } wheels = [ - { url = "https://files.pythonhosted.org/packages/e5/43/9f0e9bfb3ce02cbf7747aa2185c48a9d6e42ba95736a5e8f511a5054d976/sphinx_argparse-0.5.2-py3-none-any.whl", hash = "sha256:d771b906c36d26dee669dbdbb5605c558d9440247a5608b810f7fa6e26ab1fd3", size = 12547, upload-time = "2024-07-17T12:08:06.307Z" }, + { url = "https://files.pythonhosted.org/packages/51/9e/c41d68be04eef5b6202b468e0f90faf0c469f3a03353f2a218fd78279710/sphinx_book_theme-1.1.4-py3-none-any.whl", hash = "sha256:843b3f5c8684640f4a2d01abd298beb66452d1b2394cd9ef5be5ebd5640ea0e1", size = 433952, upload-time = "2025-02-20T16:32:31.009Z" }, ] [[package]] -name = "sphinx-book-theme" -version = "1.1.4" +name = "sphinx-click" +version = "6.2.0" source = { registry = "https://pypi.org/simple" } dependencies = [ - { name = "pydata-sphinx-theme" }, + { name = "click" }, + { name = "docutils" }, { name = "sphinx" }, ] -sdist = { url = "https://files.pythonhosted.org/packages/45/19/d002ed96bdc7738c15847c730e1e88282d738263deac705d5713b4d8fa94/sphinx_book_theme-1.1.4.tar.gz", hash = "sha256:73efe28af871d0a89bd05856d300e61edce0d5b2fbb7984e84454be0fedfe9ed", size = 439188, upload-time = "2025-02-20T16:32:32.581Z" } +sdist = { url = "https://files.pythonhosted.org/packages/9a/ed/a9767cd1b8b7fbdf260a89d5c8c86e20e3536b9878579e5ab7965a291e55/sphinx_click-6.2.0.tar.gz", hash = "sha256:fc78b4154a4e5159462e36de55b8643747da6cda86b3b52a8bb62289e603776c", size = 27035, upload-time = "2025-12-04T19:33:05.437Z" } wheels = [ - { url = "https://files.pythonhosted.org/packages/51/9e/c41d68be04eef5b6202b468e0f90faf0c469f3a03353f2a218fd78279710/sphinx_book_theme-1.1.4-py3-none-any.whl", hash = "sha256:843b3f5c8684640f4a2d01abd298beb66452d1b2394cd9ef5be5ebd5640ea0e1", size = 433952, upload-time = "2025-02-20T16:32:31.009Z" }, + { url = "https://files.pythonhosted.org/packages/44/bd/cb244695f67f77b0a36200ce1670fc42a6fe2770847e870daab99cc2b177/sphinx_click-6.2.0-py3-none-any.whl", hash = "sha256:1fb1851cb4f2c286d43cbcd57f55db6ef5a8d208bfc3370f19adde232e5803d7", size = 8939, upload-time = "2025-12-04T19:33:04.037Z" }, ] [[package]] @@ -2866,8 +2867,8 @@ dev = [ { name = "pytest-xdist" }, { name = "ruff" }, { name = "setuptools" }, - { name = "sphinx-argparse" }, { name = "sphinx-book-theme" }, + { name = "sphinx-click" }, { name = "sphinx-issues" }, { name = "twine" }, { name = "validate-pyproject", extra = ["all"] }, @@ -2880,8 +2881,8 @@ docs = [ { name = "jupyter-book" }, { name = "msprime" }, { name = "pyfaidx" }, - { name = "sphinx-argparse" }, { name = "sphinx-book-theme" }, + { name = "sphinx-click" }, { name = "sphinx-issues" }, ] lint = [ @@ -2939,8 +2940,8 @@ dev = [ { name = "pytest-xdist" }, { name = "ruff", specifier = "==0.15.1" }, { name = "setuptools" }, - { name = "sphinx-argparse" }, { name = "sphinx-book-theme" }, + { name = "sphinx-click" }, { name = "sphinx-issues" }, { name = "twine" }, { name = "validate-pyproject", extras = ["all"] }, @@ -2953,8 +2954,8 @@ docs = [ { name = "jupyter-book", specifier = "<2" }, { name = "msprime" }, { name = "pyfaidx" }, - { name = "sphinx-argparse" }, { name = "sphinx-book-theme" }, + { name = "sphinx-click" }, { name = "sphinx-issues" }, ] lint = [