Shotgun mode overview¶

Shotgun mode profiles whole-metagenome Illumina libraries. There is no PCR and no primer: every read pair that comes off the sequencer is matched against a k-mer database, so what limits the result is sequencing depth and the reference database rather than the resolving power of one marker gene.

Is shotgun mode the right tool for your samples?

This path targets environmental communities. For host-associated or food samples a marker-gene profiler is usually more precise — see what shotgun mode is for.

The pipeline is Kraken2 for read classification and Bracken for abundance re-estimation, with decontamination and quality trimming in front of them and a combined OTU table at the end. It is selected with one config key:

mode: shotgun          # amplicon | shotgun

The mode can also be set on the command line without editing the file: --config mode=shotgun. Which mode suits which data is covered in Choosing a mode.

mode must match the data

Both modes accept {sample}_R1 / {sample}_R2 filenames, so pointing a shotgun run at amplicon data (or the reverse) does not always produce an error; the run can finish with the wrong analysis. Check mode before every run.

Shotgun mode needs reference data that MetaFlux does not fetch: a Kraken2 database for the default profiler, sylph databases and their lineage tables for profiler: sylph, and both for profiler: both. See Reference databases.

The chain¶

raw R1/R2 ──▶ BBDuk PhiX* ──▶ BBMap host* ──▶ fastp ──▶ Kraken2 ──▶ Bracken ──▶ OTU table
                _dephix         _dehost        _trim      report      report     .tsv + .biom
                                                            │
                                                            └──▶ KrakenTools per-taxon reads*

  * optional, toggled in config

sylph instead of Kraken2, or next to it

With shotgun.profiler: sylph, sylph replaces everything from Kraken2 onwards: it estimates genome abundances from k-mer sketches of the same trimmed reads and writes 03.abundance/sylph_table.tsv. With both, the two chains run side by side on the same reads and their tables are never merged. The cleaning steps are the same in every case. See sylph.

Each stage writes files tagged with the step it came off (_dephix, _dehost, _trim). When an optional step is switched off, the next step reads the output of the last step that ran, and no placeholder files are written.

The rules carry descending Snakemake priorities (download_sra 10 down to finalize_otu_table 1), so when many jobs are ready at once the scheduler works through the early stages first rather than starting classification while reads are still being fetched.

1. Input — local FASTQs or a BioProject¶

Reads come from one of two places.

Local files in input.fastq_dir. Shotgun mode accepts either naming convention and detects which one is in use:

Layout	Example	Where it usually comes from
`{sample}_R1` / `{sample}_R2`	`soil01_R1.fastq.gz`	Illumina/bcl2fastq output, MetaFlux's default
`{sample}_1` / `{sample}_2`	`SRR11605259_1.fastq.gz`	raw `fasterq-dump` output

Reads must carry the fastq.gz extension, or plain fastq if uncompressed. fq and fq.gz are accepted by the start-up check but are never resolved afterwards — the rules look only for .fastq.gz and .fastq — and in shotgun mode the missing .fastq.gz path matches the download_sra rule under _R1/_R2 naming, so the run tries to fetch each local sample from SRA rather than reporting that the extension is unsupported; with _1/_2 naming it stops with a missing-input error instead. The check itself sees only the R1 files: a mixture of extensions among them stops the run at start-up, but the R2 files are never inspected, so uniformity across the directory is not enforced.

Direct from NCBI SRA. Setting input.bioproject to a BioProject accession makes MetaFlux resolve that project's runs itself:

input:
  fastq_dir: /path/to/fastq_dir      # must be writable — downloads land here
  bioproject: PRJNA603575
  accession_list: null               # e.g. [SRR11605259, SRR11605260]

The accession is resolved through NCBI e-utils (esearch then efetch for the runinfo table) and the run list is cached in fastq_dir as {BIOPROJECT}_runinfo.csv, so a re-run does not query NCBI again. Only runs whose LibraryLayout is PAIRED are kept — single-end runs do not fit the rest of the pipeline. accession_list narrows that set further, which is how a project that mixes WGS with amplicon or RNA-Seq runs is reduced to the metagenomes.

Each run is then fetched with prefetch and unpacked with fasterq-dump (--split-files --skip-technical), and the two mates are compressed to {sample}_R1.fastq.gz / {sample}_R2.fastq.gz in fastq_dir — so downstream rules see the same filenames whether the data was downloaded or already local.

Whatever the source, the sample list is written to samples.tsv at the top of the output directory, with one row per sample: sample, source (local or bioproject), bioproject, r1, r2. It exists so downstream R or Python scripts can read the manifest instead of re-parsing the config or globbing the FASTQ directory again.

2. PhiX removal (optional)¶

Rule decontam_phix — off with shotgun.decontamination.remove_phix: false.

PhiX is the control genome spiked into most Illumina runs; those reads are real sequence but they are not part of the sample. BBDuk matches the reads against the PhiX reference bundled inside the BBTools package, so nothing is downloaded or indexed.

Input: the raw FASTQs. Output: {sample}_dephix_R1/_R2.fastq.gz (temporary — see below) plus {sample}_dephix_stats.txt, which feeds MultiQC; the read tracking table takes its PhiX counts from BBDuk's log. Consumed by the host step, or by fastp when host removal is off.

BBDuk needs R1 and R2 to be properly paired

BBDuk reads the two files in lockstep and requires them to hold the same number of records in the same order. Read length does not matter — ragged reads, and mates of a pair with different lengths, are processed normally. A crash with There appear to be different numbers of reads in the paired input files means the pairing is broken, usually because reads were filtered before deposition and orphaned mates were left in one file — common for SRA data. repair.sh, from the same BBTools package, puts the files back in step. A crash with AssertionError: List size mismatch: N vs M is different: it comes from newer bbmap releases (39.76 to 40.02), which print it on correctly paired reads too, and the cure is MetaFlux's pinned environment, not repair.sh. Details are in Decontamination.

3. Host and contaminant removal (optional)¶

Rules fetch_host_refs, build_host_index, decontam_host — off when shotgun.decontamination.host_genomes is an empty list.

Host DNA dominates many metagenomes (gut, skin, plant tissue) and, left in, it costs classification time and can be mis-assigned to microbial taxa. BBMap aligns the reads against one or more user-supplied genomes and keeps only the reads that do not map.

Input: the PhiX-filtered reads, or the raw reads when PhiX removal is off. Output: {sample}_dehost_R1/_R2.fastq.gz plus {sample}_dehost_stats.txt with the surviving pair count. Consumed by fastp.

Which reference to use — and why a masked human genome matters — is covered in Decontamination.

4. Adapter and quality trimming¶

Rule trim_adapters, always runs.

fastp removes adapter read-through, trims poor quality from both ends with a sliding window (--cut_front --cut_right), and drops reads that end up shorter than shotgun.fastp.min_read_length (default 100; around 70 suits HiSeq 2×100 data). Adapters are detected from the data itself (--detect_adapter_for_pe), so no adapter file is needed.

Input: whichever decontaminated stage ran last, or the raw reads if both are off. Output: {sample}_trim_R1/_R2.fastq.gz (temporary) and the {sample}_fastp.html / {sample}_fastp.json reports. The JSON is used three times over: by MultiQC, by the read tracking table, and by Bracken, which reads the post-trim mean read length out of it.

5. Kraken2 classification¶

Rule kraken2.

Kraken2 walks each read pair k-mer by k-mer, looking up the minimizer of each k-mer in the database — every k-mer sharing a minimizer inherits that minimizer's stored lowest common ancestor, which is why one database entry can account for a whole run of apparent k-mer hits. The taxa hit form a small tree, each node weighted by how many hits it received, and the read is labelled with the leaf of the highest-scoring root-to-leaf path (the LCA of the leaves if two paths tie). A handful of stray hits to a distant taxon therefore does not drag the whole read up the taxonomy.

shotgun:
  kraken:
    confidence: 0.15      # higher = stricter
    hit_groups: 3         # minimum distinct hit groups per read
    memory_mapping: true
    threads: null         # null = use resources.threads.kraken2

confidence is the fraction of a read's k-mers that must fall inside the clade of the taxon the read is assigned to. A read that falls short is not simply rejected: Kraken2 moves the label up the taxonomy until the score clears the threshold, and calls the read unclassified only when not even root qualifies. Raising it therefore mostly pushes assignments to coarser ranks — which leaves more work for Bracken in step 6 — and only secondarily increases the unclassified fraction. A sweep from 0.00 to 1.00 on simulated samples (Wright et al. 2023) found the best setting depends on what is being optimised: 0.60 gave the highest mean F1, while 0.15 gave the lowest L1 distance and the best recall and precision over the reads actually classified. MetaFlux defaults to 0.15, which favours accurate relative abundances; 0.15 is also the threshold used for CAMI2 benchmarking by Nyström-Persson et al. (2025).

hit_groups counts hit groups, not k-mer hits. A hit group is a run of overlapping k-mers that share one minimizer found in the database — one lookup, at one place in the read — so a read can show dozens of k-mer hits and still amount to a single group. hit_groups: 3 requires three such separate hits before a classification stands (Kraken2's own default is 2), leaving a read unclassified when all its evidence sits at one or two places. MetaFlux uses 3, following the Kraken protocol (Lu et al. 2022). On the D6331 mock, every value tested from 1 to 8 produced the same species list; see Kraken2 hit groups. The depth-scaled Bracken threshold and k-mer floor were both fitted at 3.

Input: the fastp-trimmed reads and the Kraken2 database. Output, in 02.classification/kraken2/:

File	Contents
`{sample}_report.txt`	the per-taxon summary tree, written with `--report-minimizer-data`
`{sample}_output.txt`	one line per read pair with its assignment
`{sample}_R1.fastq.gz` / `_R2.fastq.gz`	the reads Kraken2 could classify

The classified-read FASTQs are kept, not deleted, because per-taxon read extraction draws from them. The report is consumed by Bracken, by MultiQC and by read extraction, which finds the taxon and everything below it there; the per-read output is consumed by read extraction.

The report format is not the standard one

--report-minimizer-data adds two columns to the usual Kraken2 report, which shifts the rank and taxid columns. KrakenTools reads the report as it stands, but anything else parsing these reports needs to account for the shift.

Memory behaviour is the main practical concern here and is covered in Reference databases.

6. Bracken abundance re-estimation¶

Rule bracken.

Kraken2 leaves many reads assigned above species level, at whatever ancestor the k-mers agreed on. Bracken pushes those reads back down using the database's k-mer distributions, giving abundances at a chosen rank.

shotgun:
  bracken:
    tax_lev: S            # P | C | O | F | G | S (D and K are refused)
    threshold: auto       # auto = scaled to each sample's depth; or a fixed read count

Bracken's estimates depend on read length, and the database ships pre-computed distributions for several lengths (50, 75, 100, 150, 200, 250 and 300 bp). MetaFlux reads the mean post-trim read length from the sample's fastp JSON and picks the closest available distribution automatically — using a 150 bp distribution on 100 bp reads would bias the redistribution. With threshold: auto the minimum read count is computed from the sample's classified read pairs in the same step. All three small calculations (mean length, closest distribution, threshold) are done with jq from the Bracken environment, so they give the same answer on any machine.

MetaFlux pins Bracken 3.1p1, bioconda's build of 3.1 with a fix for the report: with tax_lev other than S, plain 3.1 writes 0 reads assigned directly to a taxon (third column of {sample}_report.txt) on almost every row. The abundances MetaFlux reports come from the clade counts and are the same with either build.

Input: the Kraken2 report plus the fastp JSON. Output, in 03.abundance/: {sample}_report.txt (Bracken's re-estimated report) and {sample}_output.txt (the per-taxon abundance table). The reports are consumed by the OTU-table step and by MultiQC.

Drop to genus on a shallow database

Species-level re-estimation is only as good as the database's species representation. With a small or size-capped index, tax_lev: G gives a more reliable table. The OTU table then lists genera, one row each, with an empty species slot.

The numbers in the OTU table are estimates, and they do not sum to the classified reads

Two things about Bracken's output are easy to misread, and both matter when the numbers reach a methods section or a normalisation step.

They are estimated, not observed. A value in otu_table.tsv is not a count of reads that Kraken2 assigned to that taxon. It is Bracken's estimate of how many reads originated from it, after redistributing reads that Kraken2 had left at coarser ranks. Reporting Bracken's estimates as observed assignments overstates what was measured.

They do not account for every classified read. Bracken reports only at the target rank. Reads classified above it that cannot be pushed down — because the evidence does not separate the descendants — are absent from the output altogether. So the column total of otu_table.tsv is smaller than the classified column in stats/read_tracking.txt, and the two are not meant to reconcile. When a relative abundance is computed, be explicit about which denominator was used: the OTU-table total, or the classified reads.

The four stringency settings at a glance¶

Four settings control how much evidence a species needs to be reported. All four are in the shotgun: section of the config.

Setting	Applies to	Counts	Default
`kraken.confidence`	each read pair, in Kraken2	share of the pair's k-mers inside the assigned clade	0.15
`kraken.hit_groups`	each read pair, in Kraken2	hit groups matching the database	3
`kmer_evidence`, the k-mer evidence gate	each taxon, between Kraken2 and Bracken	distinct minimizers specific to the taxon in the sample	off; when enabled, floor 300 by default, or `auto`
`bracken.threshold`, Bracken's `-t`	each taxon, in Bracken	reads assigned to the taxon's clade	`auto`: 0.00005 × classified read pairs, minimum 10

The settings act in this order. Kraken2 moves the label of a read pair up the taxonomy until its confidence reaches the threshold; a pair with fewer hit groups than the minimum is left unclassified. The gate removes taxa with too few distinct minimizers and assigns their reads to the parent taxon. Bracken removes taxa below its read threshold, together with their reads, and redistributes the reads assigned at higher ranks among the remaining taxa.

The gate targets false calls supported by many reads that match the same few conserved or shared sequences. Such calls can pass the read threshold and still fail the gate.

The k-mer length k and the minimizer length l are set when the Kraken2 database is built (35 and 31 in the standard indices; see the first line of inspect.txt). None of the four settings changes them.

Of the settings tested in the benchmark, the Bracken threshold had the largest effect on false species calls; see Choosing confidence and threshold. The gate is optional; see Species-level k-mer evidence. Both depth-scaled rules were fitted at hit groups 3 and confidence 0.15. The k-mer floor was also fitted on one set of runs at confidence 0.20. Confidence and hit groups should therefore stay at their defaults.

7. OTU table¶

Rule finalize_otu_table.

MetaFlux builds the OTU table itself, straight from the per-sample Bracken reports: one row per taxon at shotgun.bracken.tax_lev, one column per sample, and in each cell Bracken's re-estimated read count for that taxon (the clade count of its row in the sample's report). Taxa appear in the order they are first met, reading the samples in the order of the table's columns, which are sorted by sample name (as in samples.tsv).

Every row carries a lineage of eight fixed slots, one for each rank code Kraken2 prints in its reports, on any database:

d__ k__ p__ c__ o__ f__ g__ s__      domain, kingdom, phylum, class, order, family, genus, species

Where the names come from. A slot takes the name of the row, on the path from the report's root down to the taxon, whose code is exactly the slot's letter (D, K, P …). Rows with a numbered code (K1, P3, R2 …) are in-between nodes and get no slot, with one exception: a virus realm, which Kraken2 also prints with a numbered code (R2; D1 on a database built when NCBI still ranked Viruses a superkingdom), fills d__ (below).
Empty slots. A rank that the lineage does not have leaves its slot empty: k__ is empty in every row on a GTDB database, which has no kingdom rank, and on NCBI databases for the taxa NCBI places in no kingdom, most protists among them; at tax_lev: G the s__ slot is empty.
Viruses carry their realm in d__. A virus has no domain. Its d__ slot holds its realm, the top rank of virus taxonomy, and its viral kingdom sits in k__, on every NCBI database. Kraken2 has no rank code for a realm, so MetaFlux recognises realms by their rank in the database's own taxonomy file, taxo.k2d (see Reference databases). A virus without a realm keeps an empty d__, as do plasmids. Viruses itself is never written as a domain, also on a database built when NCBI still ranked it a superkingdom.
Species names are whole, exactly as the database writes them.

d__Bacteria;k__Pseudomonadati;p__Pseudomonadota;c__Gammaproteobacteria;o__Enterobacterales;f__Enterobacteriaceae;g__Escherichia;s__Escherichia coli
d__Eukaryota;k__Fungi;p__Ascomycota;c__Pichiomycetes;o__Serinales;f__Debaryomycetaceae;g__Candida;s__Candida albicans
d__Duplodnaviria;k__Heunggongvirae;p__Uroviricota;c__Caudoviricetes;o__;f__Peduoviridae;g__Peduovirus;s__Peduovirus P37

While it builds the table, the step also checks the reports and stops on two things it cannot trust: a D or K code that disagrees with the database (reports written by a Kraken2 older than 2.17, which codes NCBI's domains R1 or R2, or reports made with another database), and one taxon with two different lineages. See Troubleshooting.

Then it applies the optional taxon filter:

shotgun:
  taxonomy_filter:
    enabled: false
    keep:    []           # e.g. [d__Bacteria, d__Archaea]
    discard: []           # e.g. [g__Homo]

The filter uses the same rank-aware token matching as the amplicon side — each token must match a whole segment of the lineage, never a loose substring — and it runs on the final table only, after Bracken. It therefore changes which OTUs appear in the table, never how a read was classified. It is off by default because metagenomic sample types vary too much for a sensible universal setting. Common uses are keeping prokaryotes (keep: [d__Bacteria, d__Archaea], which drops viruses too, since their d__ slot holds a realm) and dropping residual host (discard: [g__Homo]). The tokens are checked against the eight slots and the Kraken2 database when the run starts, so a token that can never match (k__Bacteria, now d__Bacteria) stops the run before Kraken2 does any work. If the filter would remove every taxon, the step stops with an error that names the tokens and lists the d__ values in the table, rather than writing an empty table. The machinery is described in full under Keeping and discarding taxa.

Output: 03.abundance/otu_table.biom and 03.abundance/otu_table.tsv, written from the same rows and numbers, so they always agree. These are the main deliverable of a shotgun run. With profiler: both the sylph table next to them has seven slots and no k__, so compare the two by rank name (see taxonomy strings in the final tables).

8. Per-taxon read extraction (optional)¶

Rules extract_taxon_reads and compress_extracted — only run when shotgun.extract_taxa is non-empty.

KrakenTools pulls the reads assigned to a named taxon, and everything below it, into per-sample FASTQ files — the starting point for a targeted assembly or a closer look at one organism. See Per-taxon read extraction.

What lands on disk¶

out_dir/
├── 01.preprocessing/     # dephix / dehost stats, fastp JSON + HTML, host index
├── 02.classification/
│   ├── kraken2/          # Kraken2 reports, per-read output, classified reads, k-mer evidence
│   └── sylph/            # sylph sketches and sylph-tax files (profiler sylph | both)
├── 03.abundance/         # Bracken reports, otu_table.tsv / otu_table.biom
├── 04.extracted_reads/   # KrakenTools per-taxon reads (only if extract_taxa set)
├── stats/                # read_tracking.txt; kmer_evidence/ plots (profiler kraken2 | both);
│                         # sylph/ summary (profiler sylph | both)
├── multiqc/              # multiqc_report.html
├── samples.tsv           # sample → source → FASTQ manifest
└── logs/                 # per-rule logs

The cleaned reads themselves are marked temporary: Snakemake deletes the _dephix, _dehost and _trim FASTQs as soon as the rule that needed them has finished, so 01.preprocessing/ keeps QC and statistics rather than several copies of the library. What stays there is:

File	Present when
`{sample}_fastp.html` / `_fastp.json`	always
`{sample}_dephix_stats.txt`	`remove_phix: true`
`{sample}_dehost_stats.txt`	`host_genomes` is set
`refs/bbmap_index/` (BBMap host index)	`host_genomes` is set

The concatenated host reference, refs/host_combined.fna.gz, is temporary as well: it is deleted once the index has been built, which is why a URL in host_genomes is fetched again on every fresh run.

Read stats/read_tracking.txt first after a run. One row per sample, counts in read pairs throughout, columns in pipeline order:

sample → raw → [nophix] → [no_host] → trimmed → classified → unclassified

The bracketed columns appear only when the corresponding step ran. A large drop at no_host means the sample was host-dominated — but only when PhiX removal ran too. With remove_phix: false the raw column is taken from the fastp report, that is, from the reads entering fastp, which is the output of host removal; raw and no_host then count exactly the same reads and no host drop can appear however host-dominated the sample is. In that configuration read the mapping rate out of logs/decontam_host/{sample}.log instead.

A large unclassified fraction usually points at the database rather than the data — an environment poorly covered by the chosen index, not a failed run.

multiqc/multiqc_report.html collects the fastp, Kraken2, Bracken and BBTools statistics into one page.

Where to go next¶

Reference databases — picking a Kraken2 index and keeping its memory use under control.
Decontamination — PhiX and host removal in detail.
Per-taxon read extraction — pulling reads for one taxon.