Skip to content
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
45 changes: 45 additions & 0 deletions docs/classify.md
Original file line number Diff line number Diff line change
Expand Up @@ -20,6 +20,51 @@ pathotypr classify \

Exactly one input source is required: a single FASTA (`-i`), a TSV list of samples (`-l`), or several FASTA files (`--input-files`). With `-i` and `--input-files` each FASTA **record** is treated as its own genome, so a multi-contig assembly yields one row per contig; `-l` aggregates every contig of a sample under its sample name.

## Inputs

`classify` needs three things, plus an optional annotation.

| What | Flag | Format |
|---|---|---|
| Marker panel | `-m, --markers` | [Marker TSV](marker_format.md), tab-separated |
| Reference genome | `-r, --reference` | **Single-record** FASTA that the marker positions are numbered against |
| Genomes to type | one of `-i` / `-l` / `--input-files` | FASTA, plain or gzipped |
| Annotation *(optional)* | `--gff` / `--gff-files` / list column 3 | [GFF3](input-formats.md#gff3-annotation-classify); adds gene and amino-acid columns |

### Choosing the input source

The three sources are **not** interchangeable, and the difference decides how your results are grouped:

| Source | One row per | Use when |
|---|---|---|
| `-i, --input` | FASTA **record** | A single genome, or you want per-contig results |
| `--input-files` | FASTA **record**, named `[file] record` | Several assemblies, per-contig results |
| `-l, --input-list` | **sample** | Draft assemblies: every contig of a sample is aggregated under one name |

!!! warning "Draft assemblies belong in `--input-list`"
With `-i` or `--input-files`, a 200-contig draft produces **200 rows**, each
typed on the handful of markers that happen to fall on that contig, and none
of them representing the isolate. Use `-l` and give the sample a name, so its
contigs are counted together.

### Requirements

The run **fails** unless all of these hold:

- [x] The reference is a **single record**. A multi-FASTA is rejected with `contains multiple records; provide a single-record FASTA`.
- [x] Marker positions are **1-based coordinates on that exact reference**. Nothing cross-checks this, so a panel built against a different assembly produces silently wrong calls.
- [x] `--kmer-size` between **1 and 31**.
- [x] `--min-flank-bases` below **half** of `--kmer-size`, or no marker k-mer can be built.
- [x] Sample names in `-l` are **unique**; a duplicate aborts the run.
- [x] Every path in the list exists, checked before any work starts.

Rows in the marker TSV are skipped, with a log message, when they have fewer than 4 columns, a non-numeric position, an empty REF/ALT/lineage, or an allele too long for the k-mer window (`max(len(REF), len(ALT)) ≤ kmer_size − 2 × min_flank_bases`, so **11 bp** at defaults).

!!! tip "Both strands are scanned"
Contig orientation in a draft assembly is arbitrary, so each contig is
scanned forward and reverse-complemented. A marker found on several contigs,
or on both strands, still counts **once**.

## Options

| Option | Default | Required | Description |
Expand Down
7 changes: 7 additions & 0 deletions docs/input-formats.md
Original file line number Diff line number Diff line change
Expand Up @@ -2,6 +2,13 @@

pathotypr reads a small set of plain-text formats. Every file is **UTF-8**; every tabular file is **tab-separated** (`\t`) with **no spaces around the delimiters**. This page is the authoritative reference for each one — which command consumes it, the exact column semantics, and the parsing rules the tools apply.

!!! tip "Looking for one command's inputs?"
This page is organised by **format**. Each command page also carries an
**Inputs** section listing exactly what that command needs and what it
rejects: [train](train.md#inputs) · [predict](predict.md#inputs) ·
[classify](classify.md#inputs) · [split-fastq](split-fastq.md#inputs) ·
[match](match.md#inputs).

!!! abstract "Formats at a glance"

| Format | Consumed by | Flag |
Expand Down
28 changes: 28 additions & 0 deletions docs/match.md
Original file line number Diff line number Diff line change
Expand Up @@ -15,6 +15,34 @@ pathotypr match -r <references.fasta> (-i <reads.fq> ... | -l <list.tsv>) [-o <r

You must supply exactly one input source (`-i/--input` or `-l/--input-list`) and a reference database (`-r/--references`).

## Inputs

| What | Flag | Format |
|---|---|---|
| Reference database | `-r, --references` | **Multi-record** FASTA: one record per candidate genome |
| Reads | `-i, --input` or `-l, --input-list` | FASTQ, plain or gzipped |

!!! danger "`match` pools everything into a single query"
This is the one place where a sample list does **not** mean one result per
row. `match` reads every FASTQ from every row, merges them into **one**
k-mer set, and reports **one** best-matching reference. The sample-name
column is read but ignored for grouping.

List only the reads belonging to the single query you want matched. To match
several samples, run `match` once per sample.

### Requirements

- [x] The reference FASTA is **multi-record**, the opposite of what [`classify`](classify.md) and [`split-fastq`](split-fastq.md) want. Each record is one candidate genome, and its header is what gets reported.
- [x] `-k` between **1 and 31**.
- [x] Every FASTQ path exists, checked before counting starts.

!!! info "No markers and no reference coordinates"
`match` compares raw k-mer content, so it needs no marker panel, no GFF and
no shared coordinate system. That is what lets it answer "which of these
genomes is my sample closest to" for an unknown organism, before you have a
panel to type it with.

## How it works

1. **Count** k-mers across every input FASTQ file in parallel, tallying occurrence counts.
Expand Down
28 changes: 28 additions & 0 deletions docs/predict.md
Original file line number Diff line number Diff line change
Expand Up @@ -10,6 +10,34 @@
pathotypr predict -i <query.fasta> -m <model.pathotypr.zst> -o <predictions.tsv> [OPTIONS]
```

## Inputs

`predict` needs two files.

| What | Flag | Format |
|---|---|---|
| Genomes to classify | `-i, --input` | FASTA, plain or gzipped. Multi-record is fine, and **each record is classified separately** |
| Trained model | `-m, --model` | A `.pathotypr.zst` bundle written by [`train`](train.md) |

### Requirements

- [x] The model must be a **current-format bundle**. An incompatible format version is a hard error and the model has to be retrained. A bundle written by a different pathotypr version only warns.
- [x] Records **shorter than the model's k** are skipped with a warning: no k-mer can be extracted from them.
- [x] No labels are needed or read. The **full header** is carried through to the output.

!!! warning "The model dictates k, not you"
`predict` has no `-k`. The k-mer size is whatever `train` used and is stored
in the bundle, so a model trained at `k=21` always vectorises queries at
`k=21`. This is also why the minimum usable sequence length depends on the
model rather than on the query file.

!!! info "Every record gets a call"
There is no "unknown" class and no confidence threshold. A query from a
lineage absent from the training set still receives the closest label the
forest can find, distinguishable only by a **low confidence and margin**.
Read those two columns before acting on a call: see
[Output](#predictionstsv-always-written) for what they mean.

## Options

| Option | Default | Required | Description |
Expand Down
52 changes: 52 additions & 0 deletions docs/split-fastq.md
Original file line number Diff line number Diff line change
Expand Up @@ -15,6 +15,58 @@ pathotypr split-fastq -m <markers.tsv> -r <reference.fasta> (-i <reads.fq> ... |

You must supply markers (`-m`), a reference (`-r`), and exactly one input source: one or more files with `-i/--input`, **or** a sample list with `-l/--input-list`.

## Inputs

| What | Flag | Format |
|---|---|---|
| Marker panel | `-m, --markers` | [Marker TSV](marker_format.md), tab-separated |
| Reference genome | `-r, --reference` | **Single-record** FASTA the marker positions are numbered against |
| Reads | `-i, --input` or `-l, --input-list` | FASTQ, plain or gzipped |

`-i` takes one or more files that together form **one** sample. `-l` takes a
[sample list](input-formats.md#sample-list-tsv-split-fastq-and-match-input-list)
where each row is an **independent** sample with its own report: `sample_name`
followed by one or more FASTQ paths.

### Requirements

The run **fails** unless all of these hold:

- [x] The reference is a **single record**, same rule as [`classify`](classify.md).
- [x] Marker positions are **1-based on that exact reference**.
- [x] `-k` between **1 and 31**.
- [x] Sample names in `-l` are **unique**; a duplicate aborts the run.
- [x] Every FASTQ path exists, checked before scanning starts.

!!! warning "Indel markers are skipped here"
Any marker whose REF and ALT differ in length is dropped, and the count is
logged as `indels skipped`. Short reads across repetitive regions (PE/PPE
families, IS elements) give unreliable k-mer matches for insertions and
deletions. If your panel leans on indels, type assemblies with
[`classify`](classify.md) instead, which has full-length context.

!!! warning "Use plain ACGT alleles"
Marker k-mers are 2-bit encoded. A marker containing an ambiguity code or
any other non-ACGT character is **never indexed**, so it silently becomes
unmatchable rather than raising an error.

### Depth and allele fraction

Two thresholds decide whether a marker is called, and they are the ones worth thinking about before a run:

| Option | Default | What it means |
|---|---|---|
| `--min-depth` | `10` | REF + ALT k-mer hits needed at a position before it is called at all |
| `--min-alt-percent` | `95` | Share of those hits that must carry the ALT allele |

!!! info "The default is deliberately strict about mixtures"
At `95`, a position sitting at 60% ALT is **not called**, so intermediate
frequencies never reach the output. That is what you want for a clean
consensus call, and it is why mixed infections surface as **fixed markers on
incompatible lineage branches** rather than as intermediate fractions. Lower
it only if you intend to inspect within-sample frequencies, and expect
sequencing noise to come with them.

## How it works

The marker database is built once and cached in memory, then reused for every sample in the run.
Expand Down
35 changes: 35 additions & 0 deletions docs/train.md
Original file line number Diff line number Diff line change
Expand Up @@ -19,6 +19,41 @@ ACGTACGT... → label = "L2.2.1"
ACGTACGT... → label = "L4.3.4.2"
```

## Inputs

`train` needs exactly one file.

| What | Flag | Format |
|---|---|---|
| Labeled genomes | `-i, --input` | [Training FASTA](input-formats.md#training-fasta-train), plain or gzipped |

### Requirements

The run **fails** unless all of these hold:

- [x] At least **two distinct labels**. One class cannot be trained against nothing.
- [x] At least **two records**, so the train/test split leaves a non-empty training set. With `--test-split 0` a single record is accepted, but the accuracy estimate is then meaningless.
- [x] `-k` between **1 and 31**.

These are tolerated, with a warning in the log:

- Records with an **empty sequence** are skipped.
- Extra header text after the label is ignored.

!!! tip "How many genomes per class"
Nothing enforces a minimum per class, but a class represented by one or two
genomes contributes almost nothing the forest can generalise from. Aim for
**10 or more per class**, and check the per-class metrics the run prints
before trusting a model built on fewer.

!!! warning "Labels are taken verbatim"
`L4` and `L4 ` (trailing space aside) are one label, but `L4` and `L4.1` are
two unrelated classes: `train` sees flat strings and knows nothing about the
hierarchy in the name. If you want a model that separates sublineages, label
the records at the depth you want predicted. For hierarchical calls from
marker panels instead, use [`classify`](classify.md) or
[`split-fastq`](split-fastq.md).

## Options

| Option | Default | Required | Description |
Expand Down
3 changes: 2 additions & 1 deletion pathotypr-core/src/classify_split_fastq.rs
Original file line number Diff line number Diff line change
Expand Up @@ -73,7 +73,8 @@ pub struct SplitFastqArgs {
#[arg(long)]
pub nested_classification: bool,

/// Size of diagnostic k-mers (must be odd, between 11 and 31). Default: 31.
/// Size of diagnostic k-mers (1-31). An odd value centres the variant in the
/// k-mer, so 31 or 21 are the usual choices. Default: 31.
#[arg(short = 'k', long, default_value_t = split_kmer::DEFAULT_MARKER_KMER_LEN)]
pub kmer_size: usize,

Expand Down
Loading