11. Metagenomics & Microbiome Analysis

QIIME 2 is an extensible, decentralized microbiome analysis package that tracks analytical provenance through `.qza` and `.qzv` artifacts. It integrates raw demultiplexing, DADA2 ASV inference, taxonomic classification, and phylogenetic diversity metrics.

# 1. Import demultiplexed paired-end FASTQs into QIIME 2 artifact
qiime tools import \
  --type 'SampleData[PairedEndSequencesWithQuality]' \
  --input-path manifest.tsv \
  --output-path demux.qza \
  --input-format PairedEndFastqManifestPhred33V2

# 2. Denoise with DADA2 to resolve exact ASVs
qiime dada2 denoise-paired \
  --i-demultiplexed-seqs demux.qza \
  --p-trim-left-f 13 --p-trim-left-r 13 \
  --p-trunc-len-f 240 --p-trunc-len-r 200 \
  --o-table table.qza \
  --o-representative-sequences rep-seqs.qza \
  --o-denoising-stats stats.qza

# 3. Classify taxonomy against SILVA 138 Naive Bayes classifier
qiime feature-classifier classify-sklearn \
  --i-classifier silva-138-99-nb-classifier.qza \
  --i-reads rep-seqs.qza \
  --o-classification taxonomy.qza

# 4. Generate interactive stacked barplot
qiime taxa barplot --i-table table.qza --i-taxonomy taxonomy.qza --m-metadata-file sample-metadata.tsv --o-visualization taxa-bar-plots.qzv

DADA2 models Illumina amplicon sequencing errors to infer exact biological sample sequences (ASVs) without clustering by arbitrary 97% OTU similarity, outputting single-nucleotide resolved feature tables.

library(dada2)

# 1. Quality filtering and trimming
fnFs <- sort(list.files("raw_fastqs", pattern="_R1.fastq.gz", full.names=TRUE))
fnRs <- sort(list.files("raw_fastqs", pattern="_R2.fastq.gz", full.names=TRUE))
filtFs <- file.path("filtered", basename(fnFs))
filtRs <- file.path("filtered", basename(fnRs))

out <- filterAndTrim(fnFs, filtFs, fnRs, filtRs, truncLen=c(240, 160),
                     maxN=0, maxEE=c(2,2), truncQ=2, rm.phix=TRUE, compress=TRUE)

# 2. Learn error rates and run core sample inference
errF <- learnErrors(filtFs, multithread=TRUE)
errR <- learnErrors(filtRs, multithread=TRUE)
dadaFs <- dada(filtFs, err=errF, multithread=TRUE)
dadaRs <- dada(filtRs, err=errR, multithread=TRUE)

# 3. Merge pairs, construct sequence table, and remove chimeras
mergers <- mergePairs(dadaFs, filtFs, dadaRs, filtRs)
seqtab <- makeSequenceTable(mergers)
seqtab.nochim <- removeBimeraDenovo(seqtab, method="consensus", multithread=TRUE)

# 4. Assign taxonomy using SILVA reference
taxa <- assignTaxonomy(seqtab.nochim, "silva_nr99_v138.1_train_set.fa.gz", multithread=TRUE)
write.csv(as.data.frame(seqtab.nochim), "asv_abundance_table.csv")

Kraken 2 assigns taxonomic labels to short DNA reads by examining the lowest common ancestor (LCA) of exact 35-mers against a compact hash database. It classifies millions of reads in seconds with minimal memory overhead.

# 1. Run Kraken 2 on paired-end shotgun metagenome reads
kraken2 \
  --db /databases/kraken2_standard/ \
  --paired sample_R1.fq.gz sample_R2.fq.gz \
  --threads 16 \
  --confidence 0.1 \
  --report sample.kreport \
  --output sample.kraken

# 2. View top classified species in terminal
head -n 25 sample.kreport

Bayesian Re-estimation of Abundance with KrakEN (Bracken) redistributes reads assigned to higher taxonomic levels (genus, family) down to species or genus nodes using pre-computed genomic probabilities, producing accurate relative abundances.

# 1. Re-estimate species-level relative abundance from Kraken 2 report
bracken \
  -d /databases/kraken2_standard/ \
  -i sample.kreport \
  -o sample.bracken_species.tsv \
  -w sample.bracken_report.kreport \
  -l S \
  -t 10

# 2. Extract top 10 most abundant microbial species
sort -k7,7nr sample.bracken_species.tsv | head -n 11

MetaPhlAn 4 profiles the composition of microbial communities (Bacteria, Archaea, Eukaryotes, Viruses) from metagenomic shotgun sequencing data by mapping reads against millions of unique clade-specific marker genes.

# 1. Run MetaPhlAn 4 on paired-end metagenome reads
metaphlan \
  sample_R1.fq.gz,sample_R2.fq.gz \
  --input_type fastq \
  --nproc 16 \
  --bowtie2out sample.bowtie2.bz2 \
  -o sample_profiled_metagenome.txt \
  -t rel_ab_w_reads

# 2. Merge profiles across all cohort samples
merge_metaphlan_tables.py *_profiled_metagenome.txt > cohort_abundance_matrix.tsv

HUMAnN 3 (HMP Unified Metabolic Analysis Network) quantifies the presence and abundance of microbial pathways and gene families. It uses a tiered approach: nucleotide mapping against pangenomes followed by translated search against UniRef50/90.

# 1. Run HUMAnN 3 with MetaPhlAn taxonomic guidance
humann \
  --input sample_combined.fastq.gz \
  --output humann_out/ \
  --threads 16 \
  --taxonomic-profile sample_profiled_metagenome.txt \
  --protein-database /db/uniref90/

# 2. Normalize pathway abundances to copies per million (CPM)
humann_renorm_table \
  --input humann_out/sample_pathabundance.tsv \
  --output humann_out/sample_pathabundance_cpm.tsv \
  --units cpm

CheckM2 assesses the completeness and contamination of microbial isolate genomes and Metagenome-Assembled Genomes (MAGs) using machine learning (gradient boosted trees) trained on comprehensive genomic feature spaces.

# 1. Run CheckM2 predict model across bin directory
checkm2 predict \
  --threads 16 \
  --input /data/metagenome_bins/ \
  -x fa \
  --output-directory checkm2_results/ \
  --database_path /databases/CheckM2_database/uniref100.KO.1.dmnd

# 2. Filter for High-Quality MAGs (Completeness > 90%, Contamination < 5%)
awk -F'\t' '$2 >= 90 && $3 <= 5 {print $1, $2, $3}' \
  checkm2_results/quality_report.tsv > high_quality_mags.tsv