# Metagenomics pipeline tutorial The metagenomics pipeline (`viralunity meta`) classifies reads of unknown content, optionally assembles them, and — if you ask it to — automatically picks reference genomes from the classification hits and runs a per-virus consensus assembly on the fly. This page builds it up incrementally on the same SARS-CoV-2 data you used in [the consensus tutorial](consensus.md), starting with a minimal Kraken2 run and ending with the full pipeline plus reference assembly. ```{note} This page assumes you have completed [Setup](setup.md), including downloading the metagenomic databases under `databases/`. ``` ## When to use it Use `viralunity meta` when the contents of your sample are unknown, when you want to screen for unexpected viruses, when you need a community profile across many samples, or when you want both a classification and per-virus consensus genomes in one workflow. If you already know what is in the sample, the [consensus pipeline](consensus.md) is simpler and faster. ## What it does The full pipeline strings together up to five conceptual phases, most of them optional and toggled by CLI flags. Anything in `[brackets]` is an optional step that only runs when you ask for it. ```text raw FASTQ ──► [fastp QC, Illumina only] ──► [host depletion: Deacon OR minimap2 vs --host-reference] │ │ summary filter chain, applied to every classifier track: │ RPM/RPKM ──► [.nr, contig tracks] ──► bleed ──► [.neg] ──► [.ictv] │ ├─► Kraken2 on reads ──► filter chain ──► Krona (raw + filtered) ├─► [DIAMOND blastx on reads] ──► filter chain ──► Krona (raw + filtered) ├─► [MEGAHIT de novo assembly] │ ├─► [racon / medaka polish — Nanopore only] │ ├─► [Kraken2 on contigs] ──► filter chain ──► Krona pair │ └─► [DIAMOND on contigs] ──► idxstats support filter ──► [NR validation] ──► filter chain ──► Krona pair ├─► [reference-assembly checkpoint] ──► per-(family, accession) consensus FASTA └─► MultiQC report (Illumina only) ``` The output of every classifier — reads or contigs — flows through the same cumulative cross-sample summary stack: a raw count table, an RPM-normalised table (and optionally RPKM), optionally an NR-validated table (contig tracks), a bleed-filtered table, optionally a negative-control-enrichment table, and optionally an ICTV vertebrate-virus table. Each step appends one suffix to a single growing filename — `…_taxa_summary_{RPM|RPKM}[.nr].bleed[.neg][.ictv].tsv` — so the longest-named file is the fully filtered summary. Each classifier also gets a *raw* and a *filtered* Krona HTML per sample so you can compare the unfiltered hits against what survives the filters. ## Decision matrix — what to turn on | If you want… | Add… | |------------------------------------------------------------------------------------|---------------------------------------------------------------------------------------------------------------| | "What's in my sample?" — fast first look | nothing extra (Kraken2 on reads is the default) | | Protein-level sensitivity for divergent viruses | `--run-diamond-reads --diamond-database … --taxids …` | | Recover near-full genomes via de novo assembly | `--run-denovo-assembly` (auto-runs Kraken2 on contigs; opt out with `--no-kraken2-contigs`) | | Find viruses missed by Kraken2 but visible via protein homology on assembled contigs | also `--run-diamond-contigs` | | Automatic per-virus consensus genome from any hit family | `--run-reference-assembly --method … --source … --families … --viral-genomes … --viral-taxids …` | | Remove host reads first (recommended) | `--deacon-index databases/deacon_indexes/panhuman-1.idx` (or `--host-reference host.fasta`) | | Suppress lab/blank background | `--negative-controls blank1,blank2` | | Keep only vertebrate-infecting viruses (drop phages, plant/fungal/invertebrate-only) | `--run-ictv-host-filter --ictv-vertebrate-taxids-file databases/ictv/vertebrate_virus_taxids.txt` | | Confirm de novo viral contigs against the full NCBI nr database | `--run-nr-validation --nr-diamond-database … ` (needs `--run-denovo-assembly --run-diamond-contigs`) | ## Worked example — Illumina, built up incrementally We will use the same `samples_illumina.csv` from [Setup](setup.md) — the two SARS-CoV-2 amplicon samples `itps-0001` and `itps-0002`. Each step below adds one slice of the pipeline. You can stop at any step that gives you what you need. ### Step A — minimal: Kraken2 on reads ```bash viralunity meta illumina \ --sample-sheet samples_illumina.csv \ --config-file results/meta_illumina/config.yml \ --output results/meta_illumina/ \ --run-name sarscov2 \ --kraken2-database databases/kraken2 \ --krona-database databases/krona/taxonomy \ --taxdump databases/taxdump \ --threads 2 --threads-total 4 \ --no-kraken2-contigs ``` This is the smallest useful invocation: trim with fastp, run Kraken2 on the trimmed reads, build a Krona plot per sample, and write the cross-sample summary tables. No dehosting, no assembly, no DIAMOND, no reference assembly. Look at the outputs under `results/meta_illumina/sarscov2/`: ```text qc/reports/multiqc_report.html samples// ├── fastp.html ├── kraken2_reads.report.txt # Kraken2 native report ├── kraken2_reads.krona.html # interactive sunburst — all hits └── kraken2_reads.filtered.krona.html # same, but pruned to taxa that pass the bleed filter metagenomics/taxonomic_assignments/kraken2_reads/ ├── results/.report.txt ├── kraken2_reads_taxa_summary.tsv # per-rank counts (family/genus/species) ├── kraken2_reads_taxa_summary_RPM.tsv # adds total_reads + RPM columns └── kraken2_reads_taxa_summary_RPM.bleed.tsv # adds bleed_threshold + bleed_pass ``` How to read them: - **`samples//kraken2_reads.report.txt`** — Kraken2's own per-taxon report (5 columns: percent, count, count-direct, rank, taxid, name). Use it as a familiar starting point if you know Kraken2 already. - **`samples//kraken2_reads.krona.html`** — open it in a browser. The interactive sunburst lets you drill from "Viruses" down to specific species. The `.filtered.krona.html` next to it is the same picture pruned to taxa that survived the bleed filter (see [Understanding the filters](#understanding-the-filters) below). Use the raw plot to explore, the filtered plot to act. - **The three TSVs** progress: raw counts → +RPM → +bleed filter. Each row is one `(sample, tool, mode, rank, taxid)` tuple. Open the `.bleed.tsv` in your spreadsheet of choice and filter to `bleed_pass == True` to get the call list. ```{tip} For a reference output, compare against `my_results/test_meta_illumina/`. Its config (`my_results/config_meta_illumina.yml`) corresponds to Step D below — the full pipeline. ``` ### Step B — add host depletion Even for non-human samples it is worth dehosting: human DNA from lab personnel and kit components is a common low-level contaminant. We use the pre-built Deacon human index because it is much faster than minimap2 against the full human genome: ```bash viralunity meta illumina \ --sample-sheet samples_illumina.csv \ --config-file results/meta_illumina/config.yml \ --output results/meta_illumina/ \ --kraken2-database databases/kraken2 \ --krona-database databases/krona/taxonomy \ --taxdump databases/taxdump \ --deacon-index databases/deacon_indexes/panhuman-1.idx \ --threads 2 --threads-total 4 \ --no-kraken2-contigs ``` What changed: - `--deacon-index` — host depletion runs *before* any classification. New per-sample outputs at `host_filtered/.R{1,2}.filtered.fastq.gz` and a merged `host_filtered/.merged.fastq.gz` (the merged file is the denominator used to compute RPM). - Pass `--host-reference host.fasta` instead if you want minimap2-based depletion against a non-human host you already have as a FASTA. They are mutually exclusive. ### Step C — add de novo assembly + contig classifiers + DIAMOND ```bash viralunity meta illumina \ --sample-sheet samples_illumina.csv \ --config-file results/meta_illumina/config.yml \ --output results/meta_illumina/ \ --kraken2-database databases/kraken2 \ --krona-database databases/krona/taxonomy \ --taxdump databases/taxdump \ --deacon-index databases/deacon_indexes/panhuman-1.idx \ --run-denovo-assembly \ --run-diamond-reads --run-diamond-contigs \ --diamond-database databases/diamond/viral.dmnd \ --taxids databases/diamond/protein2taxid.tsv \ --threads 4 --threads-total 8 ``` What this adds: - **MEGAHIT** runs on host-filtered reads, producing `denovo_assembly/megahit//final.contigs.fa`. - **Kraken2 on contigs** is on by default whenever `--run-denovo-assembly` is set; the contig-level summary tables land in `metagenomics/taxonomic_assignments/kraken2_contigs/`. (Opt out with `--no-kraken2-contigs`.) - **DIAMOND blastx** runs on both reads (`diamond_reads/`) and contigs (`diamond_contigs/`). DIAMOND complements Kraken2 because it works at the protein level — useful for divergent viruses whose nucleotide sequence has drifted away from anything in the Kraken2 index. - The contig-level DIAMOND output is post-filtered to *supported* contigs: reads are mapped back to assembled contigs (minimap2 + samtools idxstats), and a contig is only kept if at least one read maps to it. The output is named `*.diamond.supported.tsv` for this reason. You now have **four classification streams** running in parallel — `kraken2_reads`, `kraken2_contigs`, `diamond_reads`, `diamond_contigs` — each with its own raw → RPM → bleed (→ neg) summary stack and its own raw / filtered Krona pair under `samples//`. If you want the protein-level signal but only on contigs (and not on reads), drop `--run-diamond-reads`. Conversely, `--no-run-kraken2-reads` turns off the read-level Kraken2 stream. ### Step D — turn on reference assembly This is the headline feature that makes the meta pipeline more than a classifier. When `--run-reference-assembly` is on, a Snakemake checkpoint reads the post-filter taxa tables, picks one or more reference genomes per sample, and runs reference-guided consensus assembly automatically — the same alignment/consensus rules used by `viralunity consensus`. The checkpoint reads the **most-filtered** table in the chain — the longest-named `*_taxa_summary_*.tsv` for that track, i.e. it picks up `.nr`, `.neg`, and `.ictv` automatically when those steps ran; taxa that explicitly fail `bleed_pass` or `neg_pass` are dropped before selection (NA is kept). So contaminants suppressed by the cross-sample filters do **not** trigger reference assembly. (Older outputs that lack these tables fall back to the raw counts table.) ```bash viralunity meta illumina \ --sample-sheet samples_illumina.csv \ --config-file results/meta_illumina/config.yml \ --output results/meta_illumina/ \ --kraken2-database databases/kraken2 \ --krona-database databases/krona/taxonomy \ --taxdump databases/taxdump \ --deacon-index databases/deacon_indexes/panhuman-1.idx \ --run-denovo-assembly \ --run-diamond-reads --run-diamond-contigs \ --diamond-database databases/diamond/viral.dmnd \ --taxids databases/diamond/protein2taxid.tsv \ --run-reference-assembly \ --method kraken2 \ --source reads \ --families Coronaviridae \ --reads-count 100 \ --viral-genomes databases/virus_genomes/viral.genomes.fasta \ --viral-taxids databases/virus_genomes/genome2taxid.tsv \ --threads 4 --threads-total 8 ``` The new flags: - `--run-reference-assembly` — turn the feature on. - `--method` — which classifier's output drives the selection (`kraken2`, `diamond`, or `both`). Required. - `--source` — which level (`reads`, `contigs`, or `both`). Required. - `--families` — comma-separated list of viral families to target. Hits to any other family are ignored. Defaults to a broad list of common pathogens; here we restrict to coronaviruses. - `--reads-count` (or `--contigs-count` for contig sources) — minimum hits to trigger assembly for that family. Below the threshold, no reference is picked. - `--viral-genomes` + `--viral-taxids` — the genome FASTA and `genome2taxid.tsv` produced by `viralunity get-databases virus-genome`. The pipeline pulls reference sequences from these. Two new top-level outputs appear in `results/meta_illumina/sarscov2/`: - **`reference_targets.tsv`** — one row per `(sample, ref_key, accession)` selected by the checkpoint. The `ref_key` is `{family}_{accession}` (e.g. `Coronaviridae_NC_045512`), so two coronaviruses detected in one sample produce two separate ref_keys and two consensus runs. - **`assembly//.../consensus/final_consensus/.consensus.fasta`** — the per-sample consensus against the auto-selected reference. The full output tree under `assembly//` mirrors what `viralunity consensus` produces (references, BAMs, VCFs, coverage stats). See [Reference assembly under the hood](#reference-assembly-under-the-hood) below for what the checkpoint actually does. ## Nanopore variant For Nanopore data the same incremental story applies. The differences from Illumina: - **No fastp** — long reads skip QC; the variant caller is expected to absorb noise. - **Polishing is opt-in** — `--run-polish-racon` and/or `--run-polish-medaka --medaka-model r941_min_high_g360` polish the MEGAHIT contigs before contig-level classification. Pick a Medaka model that matches your basecaller. - **Variant caller in the reference-assembly path** is Clair3. Use `--clair3-model` to pick a model (default `r1041_e82_400bps_sup_v500`). A full Nanopore run with everything turned on: ```bash viralunity meta nanopore \ --sample-sheet samples_nanopore.csv \ --config-file results/meta_nanopore/config.yml \ --output results/meta_nanopore/ \ --kraken2-database databases/kraken2 \ --krona-database databases/krona/taxonomy \ --taxdump databases/taxdump \ --deacon-index databases/deacon_indexes/panhuman-1.idx \ --run-denovo-assembly \ --run-polish-medaka --medaka-model r941_min_high_g360 \ --run-diamond-reads --run-diamond-contigs \ --diamond-database databases/diamond/viral.dmnd \ --taxids databases/diamond/protein2taxid.tsv \ --run-reference-assembly \ --method kraken2 --source reads \ --families Coronaviridae \ --viral-genomes databases/virus_genomes/viral.genomes.fasta \ --viral-taxids databases/virus_genomes/genome2taxid.tsv \ --clair3-model r1041_e82_400bps_sup_v500 \ --threads 4 --threads-total 4 ``` Outputs follow the same layout as the Illumina full run; the differences are the missing `qc/reports/multiqc_report.html` and the presence of `denovo_assembly/megahit//polished.fasta` (when Medaka polishing is on). ## Understanding the filters Real metagenomic data is noisy. Without cross-sample filtering, every sample looks like it contains thirty different viruses — most of them index-hopping, kit contamination, or extremely low-level cross-sample bleed during sequencing. ViralUnity's filtering chain is designed to suppress those signals while preserving real positives. The taxa summary is a single filename that grows one suffix per step — `_taxa_summary_{RPM|RPKM}[.nr].bleed[.neg][.ictv].tsv` — so the longest-named file for a track is its fully-filtered summary. Each step: | Suffix (cumulative) | What it adds | |-----------------------|---------------------------------------------------------------------------------------------------------------| | `_taxa_summary.tsv` | Raw counts, one row per `(sample, tool, mode, rank, taxid)`. | | `…_RPM` / `…_RPKM` | The normalisation base (one or the other). `_RPM` adds `total_reads` + `rpm`; `_RPKM` also adds `genome_length_bp`, `n_genomes`, `rpkm` (only when `--viral-genomes`/`--viral-taxids` are set). | | `….nr` | **Contig tracks only, `--run-nr-validation`.** Adds `nr_pass`; contigs the NR LCA calls non-viral are removed (see below). | | `….bleed` | Adds `bleed_metric`, `bleed_max`, `bleed_threshold`, `bleed_applied`, `bleed_pass`. (Always produced.) | | `….neg` | `--negative-controls` only. Adds enrichment stats (`neg_metric`, `fold_enrichment`, `log10_ratio`, `z_score`, `neg_pass`, …). | | `….ictv` | `--run-ictv-host-filter` only. Removes taxa outside the ICTV vertebrate-infecting-virus allowlist (see below). | Row-removing steps (`.nr`, `.ictv`) also write a `*.dropped.tsv` sidecar listing what they removed, for audit. ### RPM normalisation Different samples have different total read counts. To compare a taxon's signal across samples we normalise by sequencing depth: `rpm = count / total_reads × 1,000,000`, where `total_reads` is the number of reads in the merged host-filtered FASTQ for that sample. The pipeline computes this per sample and joins it onto the summary. `total_reads` always reflects the correct denominator: dehosted reads (when dehosting is on), post-QC reads (Illumina with dehosting off), or raw reads (Nanopore with dehosting off). ### RPKM normalisation (optional) When you supply `--viral-genomes` (a RefSeq viral FASTA) and `--viral-taxids` (a genome2taxid mapping), the pipeline builds a per-taxon median genome-length table at every rank (family, genus, species) and computes: ```text rpkm = rpm × 1000 / genome_length_bp ``` RPKM at genus and family level is approximate (based on the median genome length of all accessions under that node). Taxons with no matching genome length get `rpkm = NA`. > **Segmented viruses.** `genome_length_bp` is the median length across the individual RefSeq *records* under a taxon, and each segment of a segmented virus (influenza, bunyaviruses, etc.) is a separate record. So for a segmented virus this is the median **segment** length, not the summed genome length — its RPKM is normalised per-segment and is therefore inflated and not directly comparable to a monopartite virus's RPKM. This does *not* affect the bleed or negative-control pass/fail (both are within-taxon ratios in which the constant genome length cancels); it affects absolute RPKM values, cross-taxon comparison, and which segmented taxa clear the `bleed_rpkm_floor` gate. The same assumption carries into `largest_contig_ref_coverage_pct` (see the `.ctgstats` note below). ### Bleed filter For each taxon, the pipeline looks at its maximum value across all samples, and sets a threshold at a small fraction of that maximum: ```text bleed_threshold = bleed_max * bleed_fraction # bleed_fraction = 0.005 by default bleed_pass = value >= bleed_threshold ``` The comparison metric (`bleed_metric` column) is **RPKM when `--viral-genomes` is supplied and the taxon has a genome length, otherwise RPM** — chosen per taxon. Because the bleed test is a *within-taxon* ratio across samples, and a taxon's genome length is constant, switching RPM→RPKM rescales every value in the group by the same factor and leaves `bleed_pass` unchanged; the only thing it changes is the floor gate below (RPKM values sit on a different scale). For example, if *Coronaviridae* hits 1000 RPM in sample A and 0.2 RPM in sample B: - `bleed_max` = 1000 - `bleed_threshold` = 1000 × 0.005 = 5 RPM - Sample A's row passes (`value=1000 ≥ 5`). - Sample B's row fails (`value=0.2 < 5`) — flagged as likely cross-sample bleed. If `bleed_max` is itself very small (below the metric's floor — `bleed_rpm_floor`, default 1.0, or `bleed_rpkm_floor`, default 0.1), the filter is a no-op and `bleed_applied` is `False` — there is no reliable signal to filter against, so every row is preserved. Tune the strictness with `--bleed-fraction` (default `0.005`). Lower values (`0.001`) are stricter; higher values (`0.01`) are more permissive. The floors are YAML-only knobs (`bleed_rpm_floor`, `bleed_rpkm_floor`); edit the generated config and rerun Snakemake to change them. ### Negative-control enrichment filter If you sequence blanks (no-template controls, extraction blanks, etc.) alongside real samples, declare them as negative controls: ```bash viralunity meta illumina --negative-controls blank1,blank2 … ``` The sample IDs must match the sample sheet exactly (no `sample-` prefix). For each taxon, the filter compares every biological sample to the distribution of the chosen metric (RPKM when available, else RPM) observed across the negative controls: | Controls | Gate | Option | |---|---|---| | 0 | no filter (`neg_pass = NA`) | — | | 1 | `log10_ratio ≥ threshold` | `--log10-ratio-threshold` (default 1.0, i.e. 10-fold) | | ≥ 2 | `z_score ≥ threshold` | `--z-score-threshold` (default 3.0) | | ≥ 2, SD = 0 | falls back to log10-ratio | — | `log10_ratio = log10((sample + pc) / (control_mean + pc))` where `pc` is `--enrichment-pseudocount` (default 1.0). All metrics are recorded in `*.neg.tsv` for full traceability: `fold_enrichment`, `log10_ratio`, `z_score`, `control_mean`, `control_sd`, `neg_metric`, `neg_decision`, and the thresholds and pseudocount used. Control statistics use **zero-fill**: a taxon not detected in a given control counts as `0` there, so `control_mean`/`control_sd` are computed over *all* declared controls (not only the ones where the taxon happens to appear), and the z-score uses that same denominator. When the control SD is `0` — a taxon absent from every control, or seen at an identical level in all of them — the z-score is `NA` (no division by zero) and the gate falls back to the log10-ratio. The bleed and negative filters compose: a row appears as a *call* only if it has `bleed_pass == True` **and** (when negative controls were configured) `neg_pass == True`. Taxa absent from the control rows are given a zero background — they pass the enrichment gate easily (conservative choice). **Aggregate (pooled) control.** Alongside the per-control z-score, the step also reports a complementary view that treats *all* negative controls as one pooled library: `pooled_control_metric` = Σ(control metric × control library size) / Σ(control library size) — i.e. raw reads pooled across controls, so a deeply-sequenced control legitimately counts more (a taxon absent from a control contributes 0). From it come `agg_fold_enrichment`, `agg_log10_ratio`, and `agg_fold_enrichment_10x_pass`/`agg_fold_enrichment_100x_pass`. This is *diagnostic only* — it does not change `neg_pass` — but it catches widespread, high-variance contamination whose per-control variance inflates the z-score denominator and can mask the signal. (If `total_reads` is unavailable, controls are weighted equally, making the pool a plain zero-filled mean.) ### NR validation (optional, contig tracks) `--run-nr-validation` adds a confirmation step for de novo **viral contigs** (contig tracks only; requires `--run-denovo-assembly` and `--run-diamond-contigs`). All samples' candidate viral contigs are combined and searched once with `diamond blastx` against a full NCBI **nr** database (`--nr-diamond-database` — a BLAST+ nr db read via `diamond prepdb`, or a native `.dmnd`; set one up with `viralunity get-databases nr`). Each contig gets an LCA consensus across its hits; it is kept only if at least a fraction `--nr-consensus-threshold` (default `0.5`) of its hits agree it is viral. The verdict is written into the contig tracks as an `nr_pass` column at the `.nr` step, and failing contigs are removed (logged to a `*.dropped.tsv` sidecar). Tune with `--nr-evalue`, `--nr-max-target-seqs`, `--nr-sensitivity` (default `fast`). This catches the common false positive where a de novo contig gets a viral hit from the (smaller) viral DIAMOND database but is really host/bacterial/vector sequence that only the full nr database can disambiguate. ### Largest-contig statistics (diamond_contigs, with `--viral-genomes`) When `--viral-genomes` is supplied, **both contig tracks** (`diamond_contigs` and `kraken2_contigs`) add three cheap columns per taxon (the `.ctgstats` chain step): `largest_contig_bp` — the length of the largest de novo viral contig assigned to that taxon — `largest_contig_ref_coverage_pct` — that length as a percentage of the reference genome (`largest_contig_bp / genome_length_bp × 100`) — and `largest_contig_median_depth` — that contig's median per-position depth. Each track remaps the host-filtered reads to its own viral contigs (the contigs it classified as viral) and runs `samtools depth -a`, so the depth reflects that track's own assignment. Together the three columns are a quick proxy for how much of the genome assembled and how deeply it was covered: e.g. a 7,000 bp contig for a 10 kb virus gives `largest_contig_ref_coverage_pct ≈ 70` at, say, 1,000× median depth — suggesting most of the genome is well covered without any extra heavy compute, a handy signal for deciding whether a reference assembly is worth attempting. The percentage is reported raw and **uncapped**, so a value above 100% means the largest contig is longer than the median reference (over-assembly, a chimeric contig, or a longer strain), and it is `NA` where no reference length is available (so it is approximate at family/genus ranks). *Caveat:* contig length is a proxy for, not a direct measure of, the fraction of the reference genome covered (a long contig can still be partial or chimeric) — this is a preliminary estimate, not a substitute for read remapping to the reference and consensus calling. It also assumes a **monopartite (non-segmented) genome**: `largest_contig_ref_coverage_pct` divides by `genome_length_bp`, which for a segmented virus is a median *segment* length (see the RPKM note above), while a de novo contig spans at most one segment. So for segmented viruses (influenza, bunyaviruses, etc.) this column reflects the largest single segment relative to the median segment — routinely >100% and *not* a whole-genome completeness measure. The read-level tracks (`kraken2_reads`, `diamond_reads`) have no contigs and so do not carry these columns. Alongside `nr_pass`, the `.nr` step appends `nr_is_virus`, `nr_species_correct`, `nr_correct_species` (NR's corrected species name where it confidently disagreed with the original call), and `final_species`. `final_species` is the confirmed species call — NR's correction when NR disagreed, otherwise the original `name` — so a single column carries the best available taxonomy for every row. ### ICTV vertebrate-virus filter (optional) `--run-ictv-host-filter` drops taxa outside a curated allowlist of vertebrate-infecting virus lineages — phages, plant/fungal/algal/invertebrate-only viruses, and giant viruses. It applies to **all four tracks** and is the last step in the chain (`.ictv`). The allowlist is a taxid file passed with `--ictv-vertebrate-taxids-file`, built from the ICTV Virus Metadata Resource (VMR) by `viralunity get-databases ictv-vertebrate-taxids` (which wraps `build_ictv_vertebrate_taxids.py`). Matching is lineage-aware: a taxon is kept when any ancestor in its lineage is on the allowlist, so strain- and species-level hits under an allowed genus/family survive. Removed rows are logged to a `*.dropped.tsv` sidecar. Use it when your samples are vertebrate (clinical/animal) and you want to suppress the large background of environmental phages and non-vertebrate viruses that classifiers surface. ### Filtered Krona plots Every classifier produces two Krona HTMLs per sample — `*.krona.html` (built directly from the raw classifier output) and `*.filtered.krona.html` (pruned to reads/contigs whose lineage has a passing ancestor at family/genus/species). The filtered plot is the visual equivalent of the `*_RPM.bleed[.neg].tsv` call list. Strain-level hits whose species or genus passes are preserved; rows with `taxid==0` are dropped. ## Reference assembly under the hood When you turn on `--run-reference-assembly`, a Snakemake checkpoint inspects the post-filter TSV for the source you chose (`--method` × `--source` selects one of the four classifier/level combinations) and decides per sample which reference genomes to assemble against. It reads the most-filtered table in that track's chain (the longest-named `*_taxa_summary_*.tsv`, including `.nr`/`.neg`/`.ictv` when enabled), and only considers taxa that did not explicitly fail `bleed_pass`/`neg_pass`. Two strategies are available: - **`taxid` (default)** — for each passing taxon assigned to a target family, look up the taxid directly in `viral_taxids` (`genome2taxid.tsv`). If the exact taxid is not present (often the classifier picked a strain that is below species in the NCBI taxonomy), resolve to the species-level ancestor via `taxdump` and retry. Fast and deterministic; assumes your classifier database and your `viral_genomes` database came from compatible RefSeq releases. - **`similarity`** — BLASTN each de novo contig (from MEGAHIT) against `viral_genomes`. For each contig, keep the best-bitscore hit that passes `--blast-qcov` and `--blast-pident`, then validate that its taxid lives under a target family via `taxdump`. Requires `--run-denovo-assembly`, and uses the BLAST index that `viralunity get-databases virus-genome` builds alongside the FASTA. Slower than `taxid` but finds closer relatives when the classifier database and the genome database diverge. Whichever strategy fires, the selected `(sample, family, accession)` triples are written to `reference_targets.tsv`. From there, each `ref_key = {family}_{accession}` triggers a full reference-guided consensus run, with output keyed under `assembly//`. When one sample has hits to multiple families, you get multiple `ref_key` directories. See [Notes — Reference selection strategies](../notes.md#reference-selection-strategies) for the full comparison table and worked invocations. ## Tuning knobs Sensitivity / specificity of detection: - **`--bleed-fraction`** (`0.005`) — lower = stricter cross-sample bleed filter. - **`--z-score-threshold`** (`3.0`) — negative-control gate with ≥ 2 controls; higher = stricter (see [Negative-control enrichment filter](#negative-control-enrichment-filter)). - **`--log10-ratio-threshold`** (`1.0`) — negative-control gate with 1 control (or when control SD = 0); higher = stricter. - **`--enrichment-pseudocount`** (`1.0`) — pseudocount added to sample and control means before fold-enrichment / log10-ratio. - **`--minimum-hit-group`** (`4`) — Kraken2's hit-group threshold; raising it makes Kraken2 more conservative. DIAMOND tuning: - **`--evalue`** (`1e-10`) — tighter = fewer false homologies. - **`--diamond-sensitivity`** (`sensitive`) — `mid-sensitive`, `more-sensitive`, `ultra-sensitive` trade compute time for sensitivity to divergent hits. Taxonomic false-positive filters: - **`--run-ictv-host-filter`** + **`--ictv-vertebrate-taxids-file`** — keep only vertebrate-infecting viruses (drops phages, plant/fungal/invertebrate-only viruses). Applies to all four tracks; see [ICTV vertebrate-virus filter](#ictv-vertebrate-virus-filter-optional). - **`--run-nr-validation`** + **`--nr-diamond-database`** — confirm de novo viral contigs against NCBI nr (contig tracks; needs `--run-denovo-assembly --run-diamond-contigs`). Tune with **`--nr-consensus-threshold`** (`0.5`), **`--nr-evalue`**, **`--nr-max-target-seqs`**, **`--nr-sensitivity`** (`fast`). Reference-assembly triggering: - **`--reads-count`** (`100`) / **`--contigs-count`** (`1`) — minimum hits to a family before reference assembly fires for that family. - **`--families`** — restrict reference assembly to a comma-separated list of families. - **`--reference-selection-strategy`** (`taxid` / `similarity`), **`--blast-qcov`**, **`--blast-pident`** — see [Notes](../notes.md#reference-selection-strategies). ## Common pitfalls - **Krona conda env must be active.** The `krona` taxonomy must be set up via `viralunity get-databases krona` *while the viralunity conda env is active*, because Krona looks for its taxonomy inside `$CONDA_PREFIX`. If `ktImportTaxonomy` runs with an empty taxonomy you get visually empty HTML plots. - **`--diamond-database` accepts either `.dmnd` or `.faa`.** If you pass a FASTA, the pipeline runs `diamond makedb` for you on the fly. The output `viral.dmnd` lands next to the FASTA. - **The `similarity` strategy needs the BLAST index.** `viralunity get-databases virus-genome` builds it alongside `viral.genomes.fasta`. If you brought your own viral genome FASTA, run `makeblastdb -dbtype nucl -in viral.genomes.fasta` manually before running `viralunity meta --reference-selection-strategy similarity`. - **`--negative-controls` IDs must match the sample sheet exactly** — no `sample-` prefix (that prefix is only used internally in the generated YAML). - **`--method` and `--source` are required when `--run-reference-assembly` is on.** They do not default; the pipeline will error out if you turn on reference assembly without specifying them. ## Reference - [Commands reference — `viralunity meta`](../commands.md#viralunity-meta): the complete flag table, including every resource and tuning option. - [Output layout — Metagenomics pipeline](../output.md#metagenomics-pipeline-viralunity-meta): the full output tree. - [Notes — Dynamic reference assembly](../notes.md#dynamic-reference-assembly): the full `taxid` vs `similarity` comparison. - [Notes — Reading the Krona "raw vs filtered" pair](../notes.md): how the filtered Krona is computed from the summary tables.