# Notes

## Negative-control enrichment filter

ViralUnity uses interpretable enrichment statistics to filter taxa that could represent contamination bleed-in from negative controls. The filter compares each sample's metric (RPKM when `--viral-genomes` is configured, otherwise RPM) to the distribution observed across negative-control samples.

**Decision tiers:**

| Controls (and taxon) | `neg_decision` | Gate |
|---|---|---|
| 0 configured | `none` | No filter; `neg_pass = NA` (treated as *keep*) |
| ≥ 1, taxon absent from every control (`control_max == 0`) | `absent_from_controls` | Ratios blanked to NA; all enrichment gates + `neg_pass` auto-pass True; z-score gates (`neg_pass_5`/`neg_pass_10`) stay NA |
| 1 | `log10_ratio` | `log10((sample + pc) / (control + pc)) ≥ --log10-ratio-threshold` |
| ≥ 2 | `z_score` | `(sample − control_mean) / control_sd ≥ --z-score-threshold` |
| ≥ 2, SD = 0 | `log10_ratio_fallback` | Same as single-control gate |

All metrics are recorded in the summary at the `.neg` step (the negative-control table in the cumulative chain, e.g. `*_RPKM.nr.bleed.neg.tsv`) for full traceability: `fold_enrichment`, `log10_ratio`, `z_score`, `control_mean`, `control_sd`, `control_median`, `control_max`, `neg_metric`, `neg_decision`, `neg_pass`, plus the thresholds and pseudocount used.

Taxa absent from *every* negative control (`control_max == 0`) are auto-passed with `neg_decision = "absent_from_controls"`: their enrichment ratios (`fold_enrichment`, `log10_ratio`, and the aggregate variants) would only reflect the sample's own abundance against the pseudocount, so those columns are blanked to NA and every enrichment gate — including `neg_pass` — is set True. Absence from the controls is itself evidence the taxon is not control-borne contamination. The z-score gates (`neg_pass_5`/`neg_pass_10`) stay NA, since a z-score is undefined without control signal. The lineage-aware Krona filter (`filter_krona_by_pass_taxids.py`) also treats NA `neg_pass` as *keep*, so increasing the number of controls can only make the filter stricter, never accidentally drop a taxon.

**RPKM and genome length:** when `--viral-genomes` points to a RefSeq FASTA (with a corresponding `--viral-taxids` genome2taxid mapping), ViralUnity builds a per-taxon genome-length table at all ranks (family, genus, species) using the median length across all accessions under each node. RPKM at higher ranks (genus, family) is an approximation based on this median.

## Segmented viruses

ViralUnity natively supports the assembly of segmented viral genomes. Instead of running the pipeline multiple times, pass multiple segment references with `--segmented-reference`:

```bash
viralunity consensus illumina \
    --sample-sheet samples.csv \
    --config-file config_segmented.yml \
    --output /path/to/output \
    --segmented-reference S=/path/to/s_segment.fasta \
    --segmented-reference L=/path/to/l_segment.fasta
```

Alternatively, pass a single multi-record FASTA to `--reference`: when it holds
more than one sequence ViralUnity treats it as segmented and splits it into one
reference per record automatically, deriving segment names from the FASTA headers
(first `|`/whitespace token, sanitised). A single combined `--gene-annotation` is
split by seqid to match. Use `--single-reference` to opt out and align all records
together as one reference. The split per-segment references are written under
`<output>/<run_name>/input_references/`.

When either form is detected, the pipeline dispatches a specialised modular workflow that processes each segment independently in parallel. Results are organised under `samples/{sample_name}/{segment_name}/`.

## Nanopore data

The metagenomics pipeline supports Nanopore data via `viralunity meta nanopore`. There is no fastp-based QC step; optional dehosting, MEGAHIT assembly, and Racon/Medaka polishing are available. Sample sheets have two columns (sample ID, path to one FASTQ/FASTA per sample).

## Reference header sanitization

The nanopore consensus workflow automatically sanitizes reference FASTA headers before use. Special characters (`/`, `\`, `|`, `,`, `~`, and spaces) in sequence identifiers are replaced with underscores (`_`). This prevents issues with downstream tools such as clair3 that use the sequence ID to create output directories.

## Dynamic reference assembly

When `--run-reference-assembly` is enabled, the metagenomics pipeline adds a post-classification step that automatically selects reference genomes and runs consensus assembly for each sample that meets the hit threshold.

### How it works

1. **Selection checkpoint** — after taxonomic classification, a Snakemake checkpoint reads the taxa summary TSVs and selects candidate references. The `--method` and `--source` flags control which classification results are used (e.g., `kraken2` on `reads`). Only samples with at least `--reads-count` reads (or `--contigs-count` contigs) assigned to one of the `--families` are selected.

2. **Reference extraction** — for each selected sample, a unique assembly key (`ref_key`) is constructed from the matched family and the reference accession (`{family}_{accession}`). The corresponding genome is extracted from `--viral-genomes`. The selection strategy controls how the reference accession is chosen (see below).

3. **Consensus assembly** — uses the same alignment + consensus rules as `viralunity consensus`, with the dynamically-selected reference. Output goes to `assembly/{ref_key}/`.

### Reference selection strategies

#### `--reference-selection-strategy taxid` (default)

For each sample that has classification hits to a target family, every taxid from the summary (Kraken2 or Diamond) is looked up in `--viral-taxids` (`genome2taxid.tsv`). The lookup is performed in two steps:

1. **Exact match** — the taxid value from the summary row is searched directly in `genome2taxid.tsv`. If one or more accessions are found, all of them are used as reference targets.
2. **Species-level fallback** — if the exact taxid is absent (e.g. the classifier assigned a strain or subspecies taxid that is below species in the NCBI taxonomy), the species-level ancestor is resolved via `--taxdump` and retried against `genome2taxid.tsv`.

`--taxdump` is used only to (a) validate that the matched taxid belongs to a target family and (b) build the reference key label — it is not used to normalise the rank of the lookup itself. A warning is emitted once per family if no accession could be found for that family in the sample.

**When to use:** when your classifier database and your `--viral-genomes` database were both built from the same RefSeq release. The taxid linkage is direct and requires no sequence comparison. Fast and deterministic.

**Limitation:** a taxid match does not guarantee the reference is the closest relative in the database — it only guarantees taxonomic identity. If you need the closest possible reference rather than any representative of the species, use `similarity` instead.

**Required databases:** `--viral-genomes`, `--viral-taxids`, and `--taxdump` (all produced by `viralunity get-databases virus-genome` / NCBI taxdump).

```bash
viralunity meta illumina \
    --sample-sheet samples.csv \
    --config-file config.yaml \
    --output /path/to/output \
    --kraken2-database /path/to/kraken2_db \
    --krona-database /path/to/krona_taxonomy \
    --taxdump /path/to/taxdump \
    --run-reference-assembly \
    --reference-selection-strategy taxid \
    --method kraken2 \
    --source reads \
    --families Coronaviridae,Orthomyxoviridae \
    --reads-count 100 \
    --viral-genomes databases/virus_genomes/viral.genomes.fasta \
    --viral-taxids databases/virus_genomes/genome2taxid.tsv \
    --threads 4 --threads-total 8
```

#### `--reference-selection-strategy similarity`

De novo assembled contigs from each sample are BLASTed against `--viral-genomes` using `blastn`. For each contig, BLAST returns hits sorted by bitscore (best first). The first hit per contig that passes both `--blast-qcov` and `--blast-pident` thresholds is selected; subsequent hits for the same contig are ignored. The selected accession is then validated against `--viral-taxids` and `--taxdump`: only hits whose taxid (at any rank — strain, species, genus, or family) traces back to a target family are kept as reference targets.

**When to use:** when sequence identity to a specific genome matters more than taxonomic label — for example, when working with divergent or novel strains where taxid-based lookup may return a reference that is genetically distant. Requires de novo assembly (`--run-denovo-assembly`) to produce the contigs used as BLAST queries.

**Limitation:** requires a pre-built BLAST index alongside `--viral-genomes` (created automatically by `viralunity get-databases virus-genome`). If no contigs pass the identity/coverage thresholds, no reference assembly is triggered for that sample.

**Required databases:** `--viral-genomes` with its BLAST index, `--viral-taxids`, and `--taxdump` (produced by `viralunity get-databases virus-genome` / NCBI taxdump).

```bash
viralunity meta illumina \
    --sample-sheet samples.csv \
    --config-file config.yaml \
    --output /path/to/output \
    --kraken2-database /path/to/kraken2_db \
    --krona-database /path/to/krona_taxonomy \
    --taxdump /path/to/taxdump \
    --run-denovo-assembly \
    --run-reference-assembly \
    --reference-selection-strategy similarity \
    --method kraken2 \
    --source reads \
    --families Coronaviridae,Orthomyxoviridae \
    --reads-count 100 \
    --blast-qcov 80 \
    --blast-pident 80 \
    --viral-genomes databases/virus_genomes/viral.genomes.fasta \
    --threads 4 --threads-total 8
```

### Strategy comparison

| | `taxid` | `similarity` |
|---|---|---|
| Requires de novo assembly | No | Yes |
| Requires `--viral-taxids` | Yes | Yes (family validation) |
| Requires `--taxdump` | Yes (family validation + species fallback) | Yes (family validation) |
| Requires BLAST index | No | Yes |
| Selection basis | Taxonomic ID match (exact → species fallback) | Best BLAST hit per contig |
| Speed | Fast | Slower (BLAST per sample) |
| Best for | Known strains with good RefSeq coverage | Divergent or novel strains |

### Required databases

| Option | File produced by | Used by |
|--------|-----------------|---------|
| `--viral-genomes` | `viralunity get-databases virus-genome` | Both strategies |
| `--viral-taxids` | `viralunity get-databases virus-genome` | Both (`taxid`: primary lookup; `similarity`: family validation) |
| `--taxdump` | NCBI taxdump | Both (family validation; species-level fallback for `taxid`) |
| BLAST index (`.nhr`/`.nin`/`.nsq`) | `viralunity get-databases virus-genome` | `similarity` only |

## Running tests

```bash
pytest -v test
```