Metagenomics Shotgun

Metagenomics: Shotgun & 16S

When to Use

• Processing whole-metagenome shotgun sequencing (WMS) data for taxonomic and functional profiling
• Doing species/strain-resolution community profiling beyond what 16S can achieve
• Annotating community metabolic/functional pathways (UniRef, MetaCyc)
• Recovering and QC'ing metagenome-assembled genomes (MAGs)
• Running 16S amplicon analysis with QIIME2, or comparing alpha/beta diversity across groups

Version Compatibility

Kraken2 ≥2.1.3, Bracken ≥2.9, HUMAnN ≥3.9 (biobakery), MEGAHIT ≥1.2.9, MetaBAT2 ≥2.15/2.17, CheckM ≥1.2 (or CheckM2 ≥1.0), QIIME2 ≥2024.2, Bowtie2 ≥2.5, Python ≥3.10 (pandas, numpy, scipy), R ≥4.3 with vegan ≥2.6.

Prerequisites

• Tools: bowtie2, samtools, kraken2 + a reference DB (standard or PlusPF), bracken, humann (+ UniRef90/ChocoPhlAn DBs, ~25 GB), megahit, metabat2, checkm, qiime2
• Python: pandas, numpy, scipy; R: vegan
• Prior concepts: FASTQ QC and read alignment (see bio-read-qc, bio-read-alignment-bowtie2-alignment)

Quick Reference

Core Workflow

Goal: Classify reads taxonomically without host contamination inflating the counts. Approach: align to the host genome and keep only unmapped read pairs, then classify with Kraken2 and re-estimate species-level abundance with Bracken (Kraken2's LCA algorithm over-assigns reads to higher ranks).

# Remove host reads
bowtie2 -x hg38 -1 R1.fq.gz -2 R2.fq.gz \
    --un-conc-gz decontam_%.fq.gz > /dev/null

# Kraken2 classification
kraken2 --db standard/ --paired --gzip-compressed \
    decontam_1.fq.gz decontam_2.fq.gz \
    --report kraken2_report.txt --output kraken2_out.txt

# Bracken species-level re-estimation (requires >= 10 reads assigned, -t 10)
bracken -d standard/ -i kraken2_report.txt \
    -o bracken_species.txt -r 150 -l S -t 10

Goal: Determine what the community is doing metabolically, not just who is there. Approach: run HUMAnN3 on decontaminated, merged reads to get gene-family and pathway abundances, then normalize and join across samples for comparison.

humann --input decontam_merged.fq.gz \
    --output humann3_out/ --threads 8

# Normalize pathways to copies-per-million
humann_renorm_table --input humann3_out/sample_pathabundance.tsv \
    --output pathways_cpm.tsv --units cpm

# Join multiple samples into one table
humann_join_tables --input humann3_outputs/ \
    --output all_pathways.tsv --file_name pathabundance

Goal: Recover genome-resolved metagenome-assembled genomes (MAGs). Approach: assemble contigs de novo, compute per-contig coverage depth from mapped reads, bin contigs by tetranucleotide composition + coverage, then QC each bin's completeness/contamination.

# Assembly
megahit -1 decontam_1.fq.gz -2 decontam_2.fq.gz \
    -o megahit/ --min-contig-len 500 -t 16

# Coverage depth for binning
bowtie2-build megahit/final.contigs.fa contigs_index
bowtie2 -x contigs_index -1 decontam_1.fq.gz -2 decontam_2.fq.gz | \
    samtools sort -o contigs.bam && samtools index contigs.bam
jgi_summarize_bam_contig_depths --outputDepth depths.txt contigs.bam

# Binning + quality control (MIMAG standard: HQ >=90% complete, <5% contam)
metabat2 -i megahit/final.contigs.fa -a depths.txt -o bins/bin --minContig 1500
checkm lineage_wf bins/ checkm_out/ -t 8 -x fa

16S QIIME2 Workflow

qiime tools import \
    --type 'SampleData[PairedEndSequencesWithQuality]' \
    --input-path manifest.csv --output-path reads.qza \
    --input-format PairedEndFastqManifestPhred33V2

qiime dada2 denoise-paired \
    --i-demultiplexed-seqs reads.qza \
    --p-trunc-len-f 250 --p-trunc-len-r 200 \
    --o-table table.qza --o-representative-sequences rep_seqs.qza

qiime feature-classifier classify-sklearn \
    --i-classifier silva138_classifier.qza \
    --i-reads rep_seqs.qza --o-classification taxonomy.qza

qiime diversity core-metrics-phylogenetic \
    --i-table table.qza --i-phylogeny rooted_tree.qza \
    --p-sampling-depth 5000 --m-metadata-file metadata.tsv \
    --output-dir diversity/

MAG Quality Standards (MIMAG)

Diversity Analysis (Python)

Goal: Quantify within-sample (alpha) and between-sample (beta) diversity from a taxa-by-sample abundance table. Approach: compute Shannon/Simpson/richness per sample, then build a Bray-Curtis distance matrix, ordinate it with classical PCoA, and test group separation with a permutation-based PERMANOVA.

import numpy as np
import pandas as pd
from scipy.spatial.distance import braycurtis


def shannon_diversity(counts) -> float:
    """Shannon H' — richness + evenness; 0 = no diversity."""
    counts = np.array(counts, dtype=float)
    p = counts[counts > 0] / counts[counts > 0].sum()
    return -np.sum(p * np.log(p))


def simpson_diversity(counts) -> float:
    """Simpson 1-D — probability two random reads differ; robust to rare taxa."""
    counts = np.array(counts, dtype=float)
    p = counts[counts > 0] / counts[counts > 0].sum()
    return 1 - np.sum(p ** 2)


def observed_richness(counts) -> int:
    return int(np.sum(np.array(counts) > 0))


def pcoa(dm: np.ndarray) -> tuple[np.ndarray, np.ndarray]:
    """Classical PCoA on a distance matrix. Returns (coords n x 2, variance_explained[:2])."""
    n = len(dm)
    H = np.eye(n) - np.ones((n, n)) / n
    B = -0.5 * H @ (dm ** 2) @ H
    eigvals, eigvecs = np.linalg.eigh(B)
    order = np.argsort(eigvals)[::-1]
    eigvals, eigvecs = eigvals[order], eigvecs[:, order]
    pos = eigvals > 0
    coords = eigvecs[:, pos] * np.sqrt(eigvals[pos])
    prop = eigvals[pos] / eigvals[pos].sum()
    return coords[:, :2], prop[:2]


def permanova(dm: np.ndarray, grouping: np.ndarray, n_perm: int = 999) -> tuple[float, float]:
    """PERMANOVA on a distance matrix. Returns (F-statistic, p-value)."""
    unique = np.unique(grouping)
    n = len(grouping)

    def f_stat(g):
        ss_tot = np.sum(dm ** 2) / n
        ss_w = sum(
            np.sum(dm[np.ix_(g == grp, g == grp)] ** 2) / (2 * (g == grp).sum())
            for grp in unique if (g == grp).sum() > 1
        )
        df_b, df_w = len(unique) - 1, n - len(unique)
        return ((ss_tot - ss_w) / df_b) / (ss_w / df_w) if df_w and ss_w else 0.0

    obs = f_stat(grouping)
    perm_count = sum(f_stat(np.random.permutation(grouping)) >= obs for _ in range(n_perm))
    return obs, (perm_count + 1) / (n_perm + 1)


otu = pd.read_csv("otu_table.tsv", sep="\t", index_col=0)  # samples x taxa
alpha = pd.DataFrame({
    "shannon": otu.apply(shannon_diversity, axis=1),
    "simpson": otu.apply(simpson_diversity, axis=1),
    "richness": otu.apply(observed_richness, axis=1),
})

otu_rel = otu.div(otu.sum(axis=1), axis=0)
n = len(otu_rel)
dist = np.zeros((n, n))
for i in range(n):
    for j in range(i + 1, n):
        dist[i, j] = dist[j, i] = braycurtis(otu_rel.iloc[i], otu_rel.iloc[j])

coords, prop_explained = pcoa(dist)
f_stat, p_value = permanova(dist, grouping=np.array(["A", "B"] * (n // 2)))

Diversity Analysis (R / vegan)

Goal: Same diversity/PERMANOVA workflow using the standard R ecology stack. Approach: vegan::diversity/specnumber for alpha metrics, vegdist + adonis2 for beta diversity and group testing, betadisper to check the equal-dispersion assumption PERMANOVA relies on.

library(vegan)

otu <- read.delim("otu_table.tsv", row.names = 1, check.names = FALSE)  # samples x taxa
meta <- read.delim("metadata.tsv", row.names = 1)

alpha <- data.frame(
  shannon  = diversity(otu, index = "shannon"),
  simpson  = diversity(otu, index = "simpson"),
  richness = specnumber(otu)
)

bc_dist <- vegdist(otu, method = "bray")
perm <- adonis2(bc_dist ~ group, data = meta, permutations = 999)
print(perm)

# PERMANOVA is sensitive to within-group dispersion, not just centroids
disp <- betadisper(bc_dist, meta$group)
anova(disp)

Parsing Outputs (Python)

import pandas as pd


def read_kraken2_report(path: str) -> pd.DataFrame:
    """Parse a Kraken2 report (6-column TSV: pct, clade_reads, direct_reads, rank, taxid, name)."""
    cols = ["pct", "clade_reads", "direct_reads", "rank", "taxid", "name"]
    df = pd.read_csv(path, sep="\t", header=None, names=cols)
    df["name"] = df["name"].str.strip()
    return df


def read_checkm(path: str) -> pd.DataFrame:
    """Parse CheckM lineage_wf qa output and flag MIMAG quality tiers."""
    df = pd.read_csv(path, sep="\t")
    df.columns = [c.strip().lower().replace(" ", "_") for c in df.columns]
    high_q = df[(df["completeness"] >= 90) & (df["contamination"] < 5)]
    med_q = df[(df["completeness"] >= 50) & (df["contamination"] < 10)]
    print(f"High-quality MAGs: {len(high_q)}")
    print(f"Medium-quality MAGs: {len(med_q)}")
    return df


report = read_kraken2_report("kraken2_report.txt")
species = report[report["rank"] == "S"].sort_values("pct", ascending=False)
print(species[["name", "pct"]].head(10).to_string(index=False))

Pitfalls

• Host decontamination is critical — failure to remove host reads inflates classification rates and skews abundances
• Kraken2 database choice — standard (archaea + bacteria + viral) vs PlusPF (adds protozoa/fungi) vs custom database changes what can be detected
• Bracken threshold — -t 10 requires ≥10 reads assigned to a taxon; lower it for low-coverage samples, but expect noisier low-abundance calls
• HUMAnN3 databases — requires UniRef90 and ChocoPhlAn reference databases (~25 GB) downloaded separately before running
• Binning quality — MetaBAT2 needs ≥2x coverage and ≥500 bp contigs; combine with CONCOCT/MaxBin2 and DAS Tool for bin refinement
• Rarefying without checking curves — rarefy only when rarefaction curves plateau, otherwise you discard real signal
• Comparing raw alpha diversity without rarefaction — deeper-sequenced samples will always appear more diverse
• PERMANOVA sensitivity to dispersion — significant adonis2 results can reflect unequal within-group variance, not just centroid differences; always pair with betadisper

See Also

• bio-read-qc — FASTQ quality control and adapter trimming before decontamination
• bio-read-alignment-bowtie2-alignment — host-genome alignment mechanics used for decontamination
• bio-microbiome-diversity-analysis — deeper alpha/beta diversity statistics
• bio-genome-assembly-metagenome-assembly — MEGAHIT/binning assembly details