Kraken2 and Bracken Tutorial: Taxonomic Profiling of Shotgun Metagenomes

A shotgun metagenome is a soup: every organism in a sample, sheared into millions of short reads, sequenced together with no labels attached. Kraken2’s job is to assign each read to a taxon, and Bracken’s job is to convert those read counts into abundance estimates you can actually compare between samples.

Both tools are fast and widely used, which is exactly why so many published profiles are subtly wrong. The pitfalls are almost never computational — they are database choice, threshold settings, and comparing samples that were not comparable to begin with.

How Kraken2 classifies a read

Kraken2 uses exact k-mer matching against a database built from reference genomes. It does not do alignment, which is why it processes tens of millions of reads in minutes. Each k-mer in a read is looked up in a compact hash table that maps it to the lowest common ancestor (LCA) of all genomes containing that k-mer.

That design has one consequence you must hold in mind: when a k-mer is shared between several organisms, the read is assigned no lower than their common ancestor. Shared sequence pushes classification upward, so the tool is conservative by construction. It is not wrong, it is deliberately non-committal.

# What a classification output looks like
# C = classified, U = unclassified, then the taxon and the LCA mapping
C	SRR1234.1	Bacteria (taxid 2)	562:12 561:4
U	SRR1234.2	0	0

Choosing a database is the analytical decision

Everything downstream inherits the database’s biases. A database missing a genus cannot find it; a database over-representing one genus will over-report it. There is no neutral choice.

# Prebuilt standard database — verify the current date suffix
wget https://genome-idx.s3.amazonaws.com/kraken/k2_standard_20240605.tar.gz
mkdir -p db_standard
tar -xzf k2_standard_20240605.tar.gz -C db_standard

# Building your own, with an explicit library list
kraken2-build --download-taxonomy --db db_custom
kraken2-build --download-library bacteria --db db_custom
kraken2-build --download-library viral --db db_custom
kraken2-build --build --db db_custom --threads 16

Record the database build date in your methods section. k2_standard_20240605 and a database built six months later give different answers on the same reads, and a reader cannot reproduce your profile without knowing which you used.

Running the classification

kraken2 \
  --db db_standard \
  --threads 16 \
  --paired \
  --use-names \
  --confidence 0.10 \
  --report sample.report \
  --output sample.kraken2 \
  sample_R1.fastq.gz sample_R2.fastq.gz

The confidence threshold is the most consequential knob and the most frequently omitted from methods sections. With no threshold, a read is classified if even a single k-mer supports the assignment, which inflates species-level calls in low-biomass samples. Raising it to 0.10 or 0.20 trades sensitivity for specificity — and changes your headline result, so it must be stated.

# Compare thresholds rather than picking one blindly
for c in 0.0 0.05 0.10 0.20; do
  kraken2 --db db_standard --threads 8 --paired --confidence $c \
    --report "sample_conf${c}.report" --output /dev/null \
    sample_R1.fastq.gz sample_R2.fastq.gz
  echo -n "confidence $c: "
  grep -P '\tS\t' "sample_conf${c}.report" | wc -l
done

That loop prints the number of species detected at each threshold. If the count collapses between 0.05 and 0.10, your sample is dominated by weak, ambiguous assignments — worth knowing before you build a story on it.

Bracken: from reads to abundances

Kraken2 reports reads assigned at each taxonomic level. A read that matched a shared region sits at genus level, so summing species counts under-reports true abundance. Bracken redistributes those higher-level counts to species using the database’s k-mer distribution, then re-estimates abundance.

bracken \
  -d db_standard \
  -i sample.report \
  -o sample.bracken_species.tsv \
  -r 150 \
  -l S \
  -t 10

The read length mismatch is a quiet source of nonsense. If you trimmed to 140 bp but tell Bracken 150, the k-mer distribution it assumes does not match your data, and abundance estimates skew. It will not warn you.

The output has a column named new_est_reads and fraction_total_reads. Use the fraction, not the raw count, whenever you compare samples — and even then, see the next section.

The comparison problem nobody warns you about

Metagenomic relative abundance is compositional: every value is a proportion of a total that is itself arbitrary, driven by sequencing depth and by how much of the sample was host DNA. A sample that is 95 percent human reads and one that is 5 percent human reads can produce identical community profiles and wildly different absolute microbial loads.

Two practical consequences. First, remove host reads before classification rather than letting them consume depth:

# Map to the host and keep only unmapped reads
bowtie2 -x refs/human -1 sample_R1.fq.gz -2 sample_R2.fq.gz \
  --very-sensitive --threads 16 \
  | samtools view -bS -f 4 - \
  | samtools fastq -1 host_depleted_R1.fq.gz -2 host_depleted_R2.fq.gz -

Second, report unclassified and host fractions alongside your profiles. A profile that silently discards 60 percent of reads is not a profile of the community, it is a profile of the part you happened to recognise.

Do not compare raw percentages across studies

Relative abundances are closed data. If one taxon rises, others must fall, even with no biological change. When comparing groups, work at a level you can defend — genus rather than species for short reads — and use composition-aware methods rather than a t-test on percentages.

What a good result looks like

# Top genera by fraction, with the unclassified fraction kept visible
awk -F'\t' 'NR==1 || $2=="G"' sample.bracken_species.tsv | head
grep -P '\tU\t' sample.report    # what fraction is unclassified?
grep -P '\tS\t' sample.report | sort -k2,2nr | head -10

Sanity checks that cost seconds and catch most disasters: does the top taxon make biological sense for this sample type? Is the unclassified fraction under about 30 percent, or did something go wrong with adapters or host depletion? Do technical replicates agree within a few percent? A profile that fails these is a pipeline problem, not a discovery.

Reporting standards

Give the database build, the confidence threshold, the Bracken read length and level, whether host reads were removed and how, and the unclassified fraction per sample. With those five facts a reader can reproduce your table; without them your profile is an assertion.

If you want to take a metagenomic profile further — differential abundance, functional prediction, or linking taxa to phenotypes — the interpretation step is where most projects stall. Ampersand Academy covers that reasoning in depth, including the compositional statistics that make cross-sample comparison honest.

Practice

1. Classify one small paired-end sample against the standard-8 database and note the unclassified fraction.
2. Re-run at confidence 0.0 and 0.20 and list the species that disappear.
3. Run Bracken at species level and compare new_est_reads with the raw species counts from Kraken2.
4. Remove host reads with bowtie2 and quantify how much the community profile changes.

Related tutorials

• 16S rRNA Amplicon Sequencing Explained: How a Read Becomes a Taxonomy Label — marker-gene profiling compared with shotgun classification
• NGS Read Quality Control: FastQC and Trimming Explained — the QC and trimming that should happen before classification
• How to Download Public Sequencing Data: SRA, GEO, and ENA Step by Step — finding and pulling public metagenomes to practise on