Bio Applied Taxonomic Profiling

Taxonomic Profiling of Shotgun Metagenomes

When to Use

• Classifying shotgun metagenomic shotgun reads to species/strain level (vs 16S amplicon genus-level resolution)
• Removing host (human/animal) contamination from gut, skin, BAL, or blood metagenome samples before classification
• Running Kraken2 → Bracken to correct LCA-inflated abundance estimates
• Deciding between Kraken2+Bracken (sensitivity, speed) and MetaPhlAn4 (low RAM, marker-gene precision)
• Computing alpha diversity (Shannon, Simpson, richness) from a Bracken/MetaPhlAn abundance table

Version Compatibility

Kraken2 ≥2.1.3, Bracken ≥2.9, Bowtie2 ≥2.5, MetaPhlAn ≥4.1, Python ≥3.10 with pandas ≥2.0, numpy ≥1.26, scipy ≥1.11.

Prerequisites

• Kraken2 reference database downloaded locally (standard ~55 GB RAM footprint; PlusPF ~100 GB adds protozoa+fungi)
• Bowtie2 host genome index (e.g. hg38) for decontamination
• QC'd, adapter-trimmed paired-end FASTQ (see bio-read-qc-fastp-workflow)
• Concepts: k-mer LCA classification, relative abundance vs read count

16S vs Shotgun

Goal: Remove host DNA before classification so it doesn't inflate counts or cause false bacterial assignments. Approach: Align reads to a host genome index with Bowtie2 and keep only the unmapped pairs (--un-conc-gz).

# Build once: bowtie2-build hg38.fa hg38_index/hg38
bowtie2 -x hg38_index/hg38 \
    -1 sample_R1.fastq.gz -2 sample_R2.fastq.gz \
    -p 8 --very-sensitive \
    --un-conc-gz decontam_%.fastq.gz \
    > /dev/null
# decontam_1.fastq.gz / decontam_2.fastq.gz = non-host read pairs

Typical host contamination: gut 1-10%, skin 10-30%, BAL 40-80%, blood/PBMC 60-95%.

import subprocess


def parse_bowtie2_contamination(bowtie2_stderr: str) -> float:
    """Extract the host-mapped read fraction from Bowtie2's stderr summary.

    Bowtie2 writes its alignment summary to stderr, e.g. a line like
    "93.42% overall alignment rate" — that rate IS the host contamination
    fraction when the index is a host genome.
    """
    for line in bowtie2_stderr.splitlines():
        if "overall alignment rate" in line:
            return float(line.strip().split("%")[0]) / 100.0
    raise ValueError("Could not find alignment rate in Bowtie2 output")

Kraken2 Classification

Goal: Assign each decontaminated read to the lowest common ancestor (LCA) of all reference genomes sharing its k-mers. Approach: Run Kraken2 against a prebuilt database with a confidence filter, then load the tab-delimited report into a DataFrame for downstream filtering.

kraken2 --db kraken2_db/ --paired --gzip-compressed --threads 16 \
    --confidence 0.1 --minimum-hit-groups 3 \
    --report kraken2_report.txt --output kraken2_output.txt \
    decontam_1.fastq.gz decontam_2.fastq.gz

• --confidence 0.1: minimum fraction of k-mers classified to a node (reduces false assignments; default 0)
• --minimum-hit-groups 3: minimum minimizer groups required before assignment
• Database: standard (55 GB RAM) vs PlusPF (100 GB, adds protozoa+fungi)

import pandas as pd

KRAKEN2_COLUMNS = [
    "pct_clade", "reads_clade", "reads_direct", "rank_code", "taxid", "name",
]


def load_kraken2_report(path: str, min_pct: float = 1.0) -> pd.DataFrame:
    """Load a Kraken2 report.txt and keep clades above `min_pct` of reads.

    Report columns: % reads in clade, reads at exact taxon, reads in clade,
    rank code (S/G/F/O/C/P/K), NCBI taxid, indented scientific name.
    """
    df = pd.read_csv(path, sep="\t", header=None, names=KRAKEN2_COLUMNS)
    df["name"] = df["name"].str.strip()
    return df[df["pct_clade"] >= min_pct].reset_index(drop=True)

Bracken Abundance Re-estimation

Goal: Correct LCA inflation — reads pushed up to genus/family because they share k-mers across species — by redistributing them back to species level. Approach: Run Bracken on the Kraken2 report with the actual read length, then use fraction_total_reads for downstream analysis.

bracken -d kraken2_db/ -i kraken2_report.txt \
    -o bracken_species.txt -w bracken_species_report.txt \
    -r 150 -l S -t 10

Alternative: MetaPhlAn4 aligns reads to clade-specific marker genes instead of k-mer LCA — ~3 GB RAM vs 55-100 GB, fewer false positives, direct compatibility with human-microbiome cohorts (HMP, curatedMetagenomicData), but lower sensitivity for organisms outside its marker set.

import numpy as np
import pandas as pd
from scipy.stats import entropy


def load_bracken_abundance(path: str) -> pd.DataFrame:
    """Load a bracken_species.txt file and return it indexed by species name."""
    df = pd.read_csv(path, sep="\t")
    return df.set_index("name")


def shannon_diversity(fractions: pd.Series) -> float:
    """Shannon index H = -sum(p_i * ln(p_i)) from a relative-abundance column
    (e.g. Bracken's fraction_total_reads or MetaPhlAn's relative_abundance/100).
    """
    p = fractions[fractions > 0].to_numpy()
    p = p / p.sum()
    return float(entropy(p))


def simpson_diversity(fractions: pd.Series) -> float:
    """Simpson's diversity 1 - D, where D = sum(p_i**2)."""
    p = fractions[fractions > 0].to_numpy()
    p = p / p.sum()
    return float(1.0 - np.sum(p ** 2))

Pitfalls

• Database choice matters: organisms absent from the reference DB are missed or misclassified to a related taxon
• LCA inflation without Bracken: reads sharing k-mers across species get pushed up to genus/family — always run Bracken before comparing species-level abundances
• Host contamination: if not removed, host reads inflate library size and cause false-positive bacterial hits from human sequence similar to microbial regions
• Read-length mismatch: Bracken's -r must match your actual sequencing read length, or redistribution weights will be wrong
• Depth-sensitive diversity: Shannon/Simpson/richness all depend on sequencing depth — rarefy to equal depth before comparing samples

See Also

• bio-metagenomics-kraken-classification
• bio-metagenomics-abundance-estimation
• bio-metagenomics-metaphlan-profiling
• bio-microbiome-diversity-analysis