IBS7048 Metagenomics — Week 09: bioBakery and k-mer taxonomy

1 Before you start

This practical uses only small public HMP-derived demonstration reads distributed with official bioBakery projects. They are appropriate for learning commands, but not for biological inference.

Do not copy clinical or personally collected shotgun reads into this course folder or publish them to RPubs. Human-associated reads may retain host sequence even after routine host filtering. Use only data whose consent, IRB/data-use terms, and de-identification procedure permit teaching and public release.

1.1 Learning goals

• Explain the different questions answered by MetaPhlAn, HUMAnN, and ShortBRED.
• Read terminal commands and the main R functions used in this practical.
• Audit one sample before scaling to a six-sample batch.
• Use bounded parallelism instead of giving every sample all server cores.
• Merge taxonomic, pathway, and targeted AMR-family profiles.
• Compare marker-based profiling with k-mer/minimizer classification on the same reads.
• Interpret special rows without treating them as automatic pipeline failures.
• Separate AMR marker evidence from phenotypic antimicrobial resistance.

1.2 170-minute teaching plan

Time	Class activity	Live computation
0–60 min	Concepts, terminology, and the bioBakery workflow	None; use slides and bundled results
60–70 min	Break
70–130 min	bioBakery practical	MetaPhlAn one sample, HUMAnN mini smoke test, ShortBRED toy/one-sample demonstration
130–140 min	Break
140–170 min	Kraken2 + Bracken k-mer practical	One plaque sample; inspect bundled six-sample results

The remaining ten minutes in a three-hour class are a buffer for login, command typing, and questions. Do not wait for a six-sample HUMAnN or Kraken2 batch to finish live; start one sample and move to the validated bundled tables. The final deck uses slides 1–47 for this timed path. Slides 48–114 preserve the instructor’s additional k-mer/reference/assembly material as an optional appendix; select from it only by replacing, not adding to, the timed theory.

2 Terms and code notation used in this practical

2.1 Analysis terms

Term	Plain-language meaning
Read	The called nucleotide sequence from one end of a library fragment. Paired-end sequencing produces two reads from the same original fragment.
Sample	The complete read set analyzed together for one specimen or time point.
Mapping	Aligning a read to a reference sequence under stated matching rules.
Classification	Assigning a read or marker signal to a label in a tool-specific taxonomy. It is not automatically the same as mapping.
Taxon	A named microbial group, such as a species or genus.
SGB	Species-level genome bin. MetaPhlAn uses SGBs to represent both named and unnamed species-level groups.
Marker	A short sequence chosen because it provides specific evidence for a taxon or protein family.
k-mer	Every sequence word of length k. For ACGTAC and k = 3, the k-mers are ACG, CGT, GTA, and TAC.
Minimizer	A selected shorter word representing a k-mer. Kraken2 stores and queries minimizers rather than every full k-mer.
LCA	Lowest common ancestor. If one sequence word occurs in several reference genomes, Kraken2 stores the lowest taxonomic node shared by those genomes.
Feature / profile	A reported unit, such as a taxon or gene family, and the set of values reported for it across samples.
Gene family	Related genes grouped because their protein sequences or functions are similar.
Pathway	A set of biochemical reactions or gene families used to perform a larger cellular process.
Abundance	An estimate of how much of a taxon or feature is present. Units depend on the tool.
Coverage	How completely reads support a reference sequence or pathway. It is different from abundance.
Stratified row	A HUMAnN row such as feature|organism, which assigns part of a function to one taxon.
Reference database/index	The sequences and annotations searched by a tool, often converted to a program-specific format for speed.

Mapping rate and classification rate are not universal quality scores. A MetaPhlAn marker-mapping percentage, a Kraken classified percentage, and a ShortBRED marker-hit count have different denominators and reference spaces.

2.2 How to read a terminal command

metaphlan INPUT --input_type fasta --nproc 2 -o PROFILE

Piece	Role
metaphlan	Command: the program that runs.
INPUT	Positional argument: the input file path.
--input_type	Long option that changes program behavior.
fasta	Value supplied to the preceding option.
-o	Short option for the output path.
PROFILE	Output file path supplied to -o.

Other shell notation used below:

• bash script.sh runs a shell script.
• A backslash (\) at the end of a line continues the same command.
• $THREADS reads a shell variable. Double quotes, as in "$THREADS", keep spaces and special characters from splitting one value into several values.
• > writes standard output to a file; 2>&1 sends error messages to the same place.
• | sends the output of the command on the left into the command on the right.
• A relative path starts from the current directory. An absolute path starts from /.
• set -Eeuo pipefail makes a teaching script stop on failed commands, unset variables, or failed commands inside a pipeline.

2.3 Make your own .sh file and run it

The shared course folder is a read-only master copy. Do not edit or change permissions on the shared wrappers. Copy only the small starter script into your own home directory:

COURSE=/mnt/4T_samsung2/courses/IBS7048/2026F/Week_09_bioBakery

mkdir -p "$HOME/IBS7048/week09"
cp "$COURSE/03_scripts/16_student_metaphlan_starter.sh" \
  "$HOME/IBS7048/week09/run_metaphlan.sh"
cd "$HOME/IBS7048/week09"
nano run_metaphlan.sh

In nano, edit only your private copy. Press Ctrl+O, then Enter to save, and Ctrl+X to exit. The first line, #!/usr/bin/env bash, is the shebang: it tells Linux which interpreter should run the file when you use ./run_metaphlan.sh.

Newly copied teaching starters are deliberately not executable. Inspect the permission string, add execute permission only for yourself, inspect it again, and run it:

ls -l run_metaphlan.sh
chmod u+x run_metaphlan.sh
ls -l run_metaphlan.sh
./run_metaphlan.sh "$USER" practice01

In a permission string such as -rwxr--r--, r means read, w means write, and x means execute. chmod u+x adds execute permission for the file’s owner (u) without opening the file to everyone. Do not use chmod 777.

• bash run_metaphlan.sh asks Bash to read the file; execute permission is not required.
• ./run_metaphlan.sh executes the file itself, so x permission and a valid shebang are required. ./ means “the file in this directory.”
• The first argument is your student/server ID. The second is a run tag. Use a new tag such as practice02 for a clean rerun instead of overwriting results.

The starter uses only 2 threads and 1 simultaneous job. It calls the validated shared wrappers and writes your output under /mnt/4T_samsung2/temp/IBS7048_week09/<student>/<run>/; it never writes to the shared course or database directories.

2.4 How to read R code

result <- metadata |>
  dplyr::select(sample_id, body_site) |>
  dplyr::filter(body_site == "stool")

• metadata and result are objects stored in memory.
• <- assigns the value on the right to the object on the left.
• select() and filter() are functions. A function accepts inputs, performs one job, and returns a value.
• Values inside parentheses are arguments. A named argument looks like na.rm = TRUE; an unnamed argument is matched by position.
• |> passes the result on its left to the first argument of the next function.
• dplyr::filter() means “use filter() from the dplyr package.” The :: notation avoids confusion when packages contain functions with the same name.

2.5 File formats

Suffix	Meaning in this practical
.fasta / .fna	Nucleotide sequences without per-base quality scores.
.fastq	Nucleotide reads plus per-base quality scores.
.faa	Protein sequences in FASTA format.
.gz / .bz2	Compressed files. The inner format still matters.
.tsv	A text table whose columns are separated by tabs.
.kreport	Kraken’s six-column hierarchical report. Bracken requires this default report format, not --use-mpa-style.
.Rmd / .html	Editable R Markdown source and its rendered webpage.
.pkl, .ffn.gz, .dmnd	MetaPhlAn metadata, HUMAnN nucleotide pangenomes, and DIAMOND protein indexes. Students use their paths but do not edit them.

3 Complementary questions from the same reads

Tool	Main question	Main output
MetaPhlAn 4	Who is there?	Relative abundance of taxa/SGBs from clade-specific markers
HUMAnN 3	What can they do?	Gene-family and pathway abundance, with organism stratification
ShortBRED	Are these selected protein families present?	Targeted family-level abundance from specific protein markers; here, CARD AMR homolog families
Kraken2 + Bracken	Which reference taxa receive whole-read k-mer/minimizer support?	Per-read taxonomic labels plus one-rank abundance estimates

MetaPhlAn does not try to classify every shotgun read. HUMAnN uses the taxonomic profile as a prescreen, searches detected-species pangenomes, and then searches remaining reads against a protein database. ShortBRED searches only a predefined set of protein-family markers. It complements HUMAnN rather than replacing it. Kraken2 takes a different route: it queries minimizers from all k-mers in each read against a genome-derived taxonomic database. Bracken then redistributes Kraken counts to one requested rank. Bracken does not classify reads again.

4 Files and fixed settings

The supplied metadata.tsv is the single list of samples used by the batch scripts.

find_week_root <- function(start = getwd(), max_up = 8) {
  p <- normalizePath(start, mustWork = TRUE)
  for (i in seq_len(max_up)) {
    if (file.exists(file.path(p, "00_README.md")) &&
        dir.exists(file.path(p, "02_practical"))) return(p)
    parent <- dirname(p)
    if (identical(parent, p)) break
    p <- parent
  }
  stop("Could not locate the Week 09 project root.")
}

week_root <- find_week_root()
metadata <- read.delim(
  file.path(week_root, "04_data", "metadata.tsv"),
  check.names = FALSE
)
knitr::kable(metadata, caption = "Six public MetaPhlAn tutorial samples")

Six public MetaPhlAn tutorial samples

sample_id	body_site	input_format	input_file	source_release
SRS014459-Stool	stool	fasta.gz	04_data/hmp_metaphlan_demo/SRS014459-Stool.fasta.gz	MetaPhlAn_4.0.2_demo
SRS014464-Anterior_nares	anterior_nares	fasta.gz	04_data/hmp_metaphlan_demo/SRS014464-Anterior_nares.fasta.gz	MetaPhlAn_4.0.2_demo
SRS014470-Tongue_dorsum	tongue_dorsum	fasta.gz	04_data/hmp_metaphlan_demo/SRS014470-Tongue_dorsum.fasta.gz	MetaPhlAn_4.0.2_demo
SRS014472-Buccal_mucosa	buccal_mucosa	fasta.gz	04_data/hmp_metaphlan_demo/SRS014472-Buccal_mucosa.fasta.gz	MetaPhlAn_4.0.2_demo
SRS014476-Supragingival_plaque	supragingival_plaque	fasta.gz	04_data/hmp_metaphlan_demo/SRS014476-Supragingival_plaque.fasta.gz	MetaPhlAn_4.0.2_demo
SRS014494-Posterior_fornix	posterior_fornix	fasta.gz	04_data/hmp_metaphlan_demo/SRS014494-Posterior_fornix.fasta.gz	MetaPhlAn_4.0.2_demo

4.1 Functions used in the setup chunk

find_week_root() is a custom function defined in this Rmd. Its inputs are start and max_up; its return value is the Week 09 folder path. It walks up from the current working directory until it finds both 00_README.md and 02_practical/. stop() produces an intentional, informative error if that folder cannot be found.

Function	Job
getwd()	Return the current working directory.
normalizePath()	Convert a path to a consistent absolute path.
file.path()	Join path components using the correct separator.
file.exists() / dir.exists()	Test whether a file or directory exists.
dirname()	Return the parent directory of a path.
read.delim()	Read a tab-separated table into an R data frame.
knitr::kable()	Display a data frame as a readable table in the HTML.

The validated server defaults are MetaPhlAn 4.2.2, HUMAnN 3.9, ShortBRED 0.9.5, Kraken2 2.1.3, Bracken 2.9, Jun23 CHOCOPhlAnSGB, CARD v4.0.1 protein-homolog markers, and the official June 2026 Kraken Standard-8 database. The classroom-safe defaults are two threads per sample and one sample at a time. The instructor can increase them only after checking total class load. The official MetaPhlAn release page lists 4.2.5 as of 2026-10-25; we keep the shared lab environment pinned to the version and database validated for this practical instead of silently changing a production environment. The lab’s validated CARD markers use v4.0.1; the current official CARD release is v4.0.2. The version difference must be reported rather than silently hidden. Standard-8 caps the Kraken database at 8 GB, so it trades some sensitivity and accuracy for a practical classroom footprint. The Bracken model is fixed at 100 bp, the closest bundled model to these 92.3–98.3 bp mean read lengths.

set -Eeuo pipefail
cd /mnt/4T_samsung2/courses/IBS7048/2026F/Week_09_bioBakery

# Individual accounts use $USER automatically. If a login is shared, set a
# unique student ID. Change the run tag when you want a clean rerun.
export IBS7048_USER_ID=20261234
export IBS7048_RUN_ID=practice01

bash 03_scripts/00_check_environment.sh

The scripts write to /mnt/4T_samsung2/temp/IBS7048_week09/${IBS7048_USER_ID}/${IBS7048_RUN_ID}/results by default. This keeps live output outside Dropbox and separates each student and each attempt. Do not use another student’s ID.

5 Part 1 — one MetaPhlAn sample

Always run one sample first. This checks the environment, database path, input format, log, mapping output, and profile before starting a batch.

bash 03_scripts/02_metaphlan_one.sh \
  SRS014476-Supragingival_plaque \
  04_data/hmp_metaphlan_demo/SRS014476-Supragingival_plaque.fasta.gz

5.1 MetaPhlAn options used by the wrapper

Option	Meaning
--input_type fasta	Interpret the input as nucleotide FASTA, even though the outer file is gzip-compressed.
--db_dir	Directory containing the MetaPhlAn marker database.
--index	Exact database release to use. Fixing this prevents silent version changes.
--nproc	Number of CPU threads used by this sample.
--mapout	Save the marker-alignment output for diagnostics and reuse.
-o	Write the taxonomic profile table to this path.

The profile contains estimated taxon abundances. The mapout contains marker-alignment evidence. The log records messages, and the versions file records software, database, and input provenance.

Inspect both the header and the first species rows:

profile=05_expected_results/live_run/metaphlan/SRS014476-Supragingival_plaque/\
SRS014476-Supragingival_plaque_profile.tsv

grep '^#' "$profile"
awk -F '\t' '$1 ~ /(^|[|])s__/ {print}' "$profile" | head

grep '^#' prints header lines beginning with #. awk -F '\t' tells awk that columns are separated by tabs. In the awk expression, $1 means the first column and ~ means “matches this regular expression.” head keeps only the first few matching rows so the terminal remains readable.

The profile header is provenance: it records the MetaPhlAn command, version, database release, and total reads. Keep it with the result.

6 Part 2 — six MetaPhlAn samples

The batch script reads metadata.tsv and uses xargs -P "$JOBS". On this server GNU parallel is not installed, so this portable pattern provides bounded parallelism without starting every sample at once.

bash 03_scripts/03_metaphlan_batch.sh
bash 03_scripts/04_merge_metaphlan.sh

The classroom defaults are THREADS=2 inside one sample and JOBS=1, meaning one sample at a time and a nominal maximum of two tool threads per student. In xargs -0 -n 2 -P "$JOBS", -0 reads null-separated values safely, -n 2 sends two arguments to each script call, and -P sets the number of concurrent calls.

6.1 Read the validated expected result

6.1.1 R function map for the next plot

Task	Functions used
Import and text matching	readLines, read.delim, grepl, sub, gsub
Select and reshape	select, pivot_longer
Group and summarize	group_by, summarise, slice_max, pull
Combine metadata	left_join, bind_rows
Plot	ggplot, aes, geom_col, scale_y_continuous, labs, theme
Control and type conversion	if, else, is.na, as.numeric, factor, %in%

The custom read_metaphlan_table(path) function accepts one profile path and returns a data frame. It first locates the real table header because MetaPhlAn places provenance lines above it.

library(dplyr)
library(tidyr)
library(ggplot2)
library(scales)

read_metaphlan_table <- function(path) {
  x <- readLines(path, warn = FALSE)
  h <- which(grepl("^#?clade_name\\t", x))[1]
  if (is.na(h)) stop("MetaPhlAn header not found: ", path)
  x[h] <- sub("^#", "", x[h])
  read.delim(
    text = paste(x[h:length(x)], collapse = "\n"),
    check.names = FALSE,
    quote = "",
    comment.char = ""
  )
}

expected_root <- file.path(week_root, "05_expected_results")
metaphlan_merged_path <- file.path(
  expected_root,
  "merged",
  "metaphlan_merged.tsv"
)

if (!file.exists(metaphlan_merged_path)) {
  cat("Validated merged MetaPhlAn table has not been bundled yet.\n")
} else {
  mpa <- read_metaphlan_table(metaphlan_merged_path)
  id_cols <- intersect(
    c("clade_name", "NCBI_tax_id", "additional_species"),
    names(mpa)
  )
  abundance_cols <- setdiff(names(mpa), id_cols)

  species_long <- mpa |>
    filter(grepl("(^|\\|)s__[^|]+$", clade_name)) |>
    mutate(
      species = sub(".*\\|s__", "", clade_name),
      species = sub("\\|.*$", "", species),
      species = gsub("_", " ", species)
    ) |>
    select(species, all_of(abundance_cols)) |>
    pivot_longer(-species, names_to = "sample", values_to = "relative_abundance") |>
    mutate(relative_abundance = as.numeric(relative_abundance))

  unclassified_long <- mpa |>
    filter(clade_name == "UNCLASSIFIED") |>
    select(all_of(abundance_cols)) |>
    pivot_longer(everything(), names_to = "sample", values_to = "relative_abundance") |>
    mutate(species = "Unclassified",
           relative_abundance = as.numeric(relative_abundance))

  top_species <- species_long |>
    group_by(species) |>
    summarise(mean_abundance = mean(relative_abundance), .groups = "drop") |>
    slice_max(mean_abundance, n = 12, with_ties = FALSE) |>
    pull(species)

  plot_data <- bind_rows(species_long, unclassified_long) |>
    mutate(species_plot = if_else(species %in% top_species, species, "Other species")) |>
    mutate(species_plot = if_else(species == "Unclassified", "Unclassified", species_plot)) |>
    group_by(sample, species_plot) |>
    summarise(relative_abundance = sum(relative_abundance), .groups = "drop") |>
    left_join(metadata |> select(sample_id, body_site), by = c("sample" = "sample_id")) |>
    mutate(body_site = factor(body_site, levels = metadata$body_site))

  ggplot(plot_data, aes(body_site, relative_abundance, fill = species_plot)) +
    geom_col(width = 0.78) +
    scale_y_continuous(labels = label_percent(scale = 1),
                       expand = expansion(mult = c(0, 0.03))) +
    labs(x = NULL, y = "Relative abundance", fill = "Species") +
    theme_minimal(base_size = 12) +
    theme(axis.text.x = element_text(angle = 35, hjust = 1),
          panel.grid.major.x = element_blank())
}

Species-level profiles from six 20k-read tutorial subsets

Read the pipeline left to right: keep species rows, move sample columns into rows with pivot_longer(), select the most abundant taxa, join the body-site metadata, and draw one stacked bar per sample.

Questions:

1. Which body sites are dominated by one or two taxa?
2. Does a 20k-read subset support a strong biological conclusion? Why not?
3. Why is MetaPhlAn marker mapping not directly comparable to Kraken classified %?

7 Part 3 — HUMAnN smoke test

The official HUMAnN demo uses tiny bundled databases and should finish quickly. It confirms that the executable, Python environment, and basic output writers work before using the full databases.

bash 03_scripts/05_humann_demo.sh

Expected files include *_genefamilies.tsv, *_pathabundance.tsv, and *_pathcoverage.tsv.

7.1 HUMAnN options and utilities

Option or command	Meaning
--input / --output	Input reads and output directory.
--output-basename	Prefix shared by the sample’s HUMAnN result files.
--threads	CPU threads used for this sample.
--taxonomic-profile	Reuse the validated MetaPhlAn profile during the full run.
--nucleotide-database	Species pangenomes searched before translated search.
--protein-database	UniRef protein database searched for remaining reads.
--remove-temp-output	Remove large intermediate files after success.
humann_join_tables	Join the same output type from several samples.
humann_renorm_table --units relab	Convert a table to sample-wise relative abundance.

The tiny official smoke test may bypass steps to test installation quickly. Do not copy smoke-test shortcuts into a full biological run without understanding which search stage they skip.

8 Part 4 — full HUMAnN analysis

Run the plaque sample first. The wrapper reuses its MetaPhlAn profile, writes intermediate files to private scratch outside Dropbox, normalizes gene-family and pathway tables to relative abundance, and moves only completed results into the results folder.

bash 03_scripts/06_humann_one_full.sh \
  SRS014476-Supragingival_plaque \
  04_data/hmp_metaphlan_demo/SRS014476-Supragingival_plaque.fasta.gz

Then run the six samples with two concurrent HUMAnN jobs and join the tables.

bash 03_scripts/07_humann_batch_full.sh
bash 03_scripts/08_join_humann_tables.sh

8.1 Functional results

feature|organism rows are stratified contributions. Rows without | are the community totals. UNMAPPED and UNINTEGRATED are special summaries rather than ordinary pathways.

In the validated full-database run, the buccal subset reconstructed five unstratified pathways; the other five tiny subsets produced gene-family signal but no reconstructed pathway. This sparse result is itself a useful depth and reference-representation check.

read_humann_table <- function(path) {
  x <- read.delim(path, check.names = FALSE, quote = "", comment.char = "")
  names(x)[1] <- "feature"
  x
}

path_candidates <- file.path(
  expected_root,
  "merged",
  "humann_pathabundance_relab.tsv"
)
path_candidates <- path_candidates[file.exists(path_candidates)]

pathway_long <- NULL
if (length(path_candidates)) {
  pathways <- read_humann_table(path_candidates[1])
  pathway_long <- pathways |>
    filter(!grepl("\\|", feature),
           !feature %in% c("UNMAPPED", "UNINTEGRATED")) |>
    pivot_longer(-feature, names_to = "sample", values_to = "abundance") |>
    mutate(
      sample = sub("_Abundance$", "", sample),
      abundance = as.numeric(abundance)
    ) |>
    left_join(metadata |> select(sample_id, body_site),
              by = c("sample" = "sample_id")) |>
    mutate(body_site = factor(body_site, levels = metadata$body_site))

}

if (!is.null(pathway_long) && nrow(pathway_long)) {
  top_pathways <- pathway_long |>
    group_by(feature) |>
    summarise(mean_abundance = mean(abundance), .groups = "drop") |>
    slice_max(mean_abundance, n = 15, with_ties = FALSE) |>
    pull(feature)

  print(pathway_long |>
    filter(feature %in% top_pathways) |>
    ggplot(aes(body_site, reorder(feature, abundance), fill = abundance)) +
    geom_tile() +
    scale_fill_viridis_c(trans = "sqrt") +
    labs(x = NULL, y = NULL, fill = "Relative\nabundance") +
    theme_minimal(base_size = 11) +
    theme(axis.text.x = element_text(angle = 35, hjust = 1)))
} else {
  demo_path <- file.path(
    expected_root,
    "humann_demo",
    "humann_demo",
    "demo_pathabundance.tsv"
  )
  demo_raw <- read_humann_table(demo_path)
  demo_value <- setdiff(names(demo_raw), "feature")[1]
  demo <- demo_raw |>
    filter(!grepl("\\|", feature),
           !feature %in% c("UNMAPPED", "UNINTEGRATED")) |>
    mutate(abundance = as.numeric(.data[[demo_value]])) |>
    slice_max(abundance, n = 15, with_ties = FALSE) |>
    mutate(relative = abundance / sum(abundance),
           feature = reorder(sub(":.*$", "", feature), relative))

  cat("No non-special pathway rows were found in the merged full run; showing the official HUMAnN demo instead.\n\n")
  print(ggplot(demo, aes(relative, feature)) +
    geom_col(fill = "#4E79A7", width = 0.75) +
    scale_x_continuous(labels = label_percent()) +
    labs(x = "Share among displayed pathways", y = NULL) +
    theme_minimal(base_size = 12))
}

Full-database HUMAnN result from six 20k-read tutorial subsets; pathway recovery was limited to the buccal subset

9 Part 5 — targeted AMR profiling with ShortBRED

9.1 Where ShortBRED fits

HUMAnN profiles broad community functional potential. ShortBRED asks a narrower question: do the reads contain markers for a predefined set of protein families? Here the targets are CARD antimicrobial-resistance protein homolog families.

Stage	Input	Output	How often it runs
ShortBRED-Identify	Proteins of interest plus a broad protein background	Family-specific peptide markers	Once when building or revising a marker database
ShortBRED-Quantify	Sample nucleotide reads plus prepared markers	Family, Count, Hits, and TotMarkerLength	Once per sample

The lab’s CARD marker set is already built against UniRef90 2026_01. Students should not rebuild it during class. The full marker-building reference area is about 225 GB; per-sample Quantify needs only the 688 KB marker FASTA.

ShortBRED term	Meaning
CARD	Comprehensive Antibiotic Resistance Database. Here it supplies curated protein homolog models used to build the lab marker set.
ARO	Antibiotic Resistance Ontology identifier used by CARD to name and relate resistance concepts.
GOI	“Genes/proteins of interest”: the protein FASTA from which markers will be built.
UniRef90	UniProt Reference Clusters grouped at 90% sequence identity. Here it is the broad protein background used to remove nonspecific regions.
WGS reads	Whole-genome shotgun nucleotide reads supplied to ShortBRED-Quantify.
TM	True marker: a peptide region specific to the target family under the chosen references and thresholds.
JM	Junction marker: a contiguous peptide region partially, but not completely, overlapped by another reference/consensus sequence; the discriminating evidence lies across that boundary.
QM	Quasi marker: a less-specific fallback when no fully unique marker is available; interpret cautiously.
RPKM-like Count	A normalized estimate related to reads per kilobase per million input reads. ShortBRED adjusts for marker length, input depth, and average read length, then reports the family median of marker estimates. It is not a raw hit count or ordinary full-gene RPKM.
--clustid	Protein identity threshold used by Identify to group GOI sequences into families.
--id / --pctlength	Minimum read-marker identity and aligned fraction used by Quantify.
--avgreadBP / --minreadBP	Average read length used in normalization and the minimum read length processed.

The six public subsets have mean read lengths of 92.3–98.3 bp, so this teaching workflow uses a common average of 100 bp and a minimum of 90 bp. Depending on the subset, 6.3–26.5% of reads are shorter than 90 bp and are therefore excluded from accepted ShortBRED hits. The settings are retained here only for a tiny workflow exercise, not a biological comparison. For NovaSeq data, calculate both SHORTBRED_AVG_READ_BP and SHORTBRED_MIN_READ_BP from the post-QC read-length distribution rather than assuming every retained read is still 150 bp.

9.2 Official mini demo

This first command uses the official ShortBRED toy proteins, reference proteins, and 100 simulated reads. It runs both Identify and Quantify and must produce at least one positive family.

bash 03_scripts/09_shortbred_demo.sh

9.3 CARD markers: one sample first

The class wrapper decompresses the input to a private regular FASTA under /mnt/4T_samsung2/temp; it does not create a large decompressed copy inside Dropbox. A regular file is required because USEARCH seeks within FASTA input. The wrapper removes the temporary FASTA when the command ends. The buccal subset is used first because the validated tiny-data run contains a small amount of CARD marker evidence.

bash 03_scripts/10_shortbred_card_one.sh \
  SRS014472-Buccal_mucosa \
  04_data/hmp_metaphlan_demo/SRS014472-Buccal_mucosa.fasta.gz

Then run six samples with bounded parallelism and merge both the normalized Count values and accepted Hits values. The merge script writes shortbred_CARD_AMR_counts.tsv plus the companion shortbred_CARD_AMR_hits.tsv.

THREADS=2 JOBS=1 bash 03_scripts/11_shortbred_card_batch.sh

source 03_scripts/00_config.sh
python3 03_scripts/12_merge_shortbred_card.py \
  "$RESULTS_DIR/shortbred_card" \
  "$RESULTS_DIR/merged/shortbred_CARD_AMR_counts.tsv"

9.4 Reading the ShortBRED output

Column	Meaning
Family	CARD/ARO-tagged protein family represented by one or more markers. The organism in this header is the reference sequence source, not proof of the organism carrying the gene in this sample.
Count	Family-level normalized abundance in RPKM-like units. ShortBRED derives marker-level estimates and reports the family median.
Hits	Number of accepted read-to-marker hits. Use it as supporting evidence, not as the main abundance value.
TotMarkerLength	Total amino-acid length of the markers associated with the family.

Hits > 0 can coexist with Count = 0 because family-level aggregation uses marker-level evidence. Conversely, one hit in a 20k-read subset can produce a large normalized Count. Always examine both columns and the original depth.

read_shortbred_merged <- function(path) {
  x <- read.delim(path, check.names = FALSE, quote = "", comment.char = "")
  stopifnot(names(x)[1] == "Family")
  x
}

shortbred_count_path <- file.path(
  expected_root,
  "merged",
  "shortbred_CARD_AMR_counts.tsv"
)

shortbred_hit_path <- file.path(
  expected_root,
  "merged",
  "shortbred_CARD_AMR_hits.tsv"
)

shortbred_counts <- read_shortbred_merged(shortbred_count_path) |>
  pivot_longer(-Family, names_to = "sample", values_to = "count") |>
  mutate(count = as.numeric(count))

shortbred_hits <- read_shortbred_merged(shortbred_hit_path) |>
  pivot_longer(-Family, names_to = "sample", values_to = "hits") |>
  mutate(hits = as.integer(hits))

shortbred_evidence <- shortbred_counts |>
  left_join(shortbred_hits, by = c("Family", "sample")) |>
  mutate(
    ARO = sub(".*[|](ARO_[0-9]+)[|].*", "\\1", Family),
    gene = sub(".*[|]ARO_[0-9]+[|]", "", Family),
    gene = sub("__.*$", "", gene),
    label = paste0(gene, " (", ARO, ")")
  ) |>
  filter(count > 0 | hits > 0) |>
  left_join(metadata |> select(sample_id, body_site),
            by = c("sample" = "sample_id"))

knitr::kable(
  shortbred_evidence |>
    select(body_site, label, count, hits) |>
    arrange(body_site, desc(count), label),
  digits = 2,
  caption = "Families with a nonzero Count or accepted Hits in the validated tiny subsets"
)

Families with a nonzero Count or accepted Hits in the validated tiny subsets

body_site	label	count	hits
anterior_nares	BRO-2 (ARO_3004762)	76.06	2
buccal_mucosa	patA (ARO_3000024)	641.03	1
buccal_mucosa	patB (ARO_3000025)	0.00	2

shortbred_positive <- shortbred_evidence |>
  filter(count > 0)

ggplot(shortbred_positive, aes(body_site, count, fill = label)) +
  geom_col(width = 0.68) +
  scale_y_continuous(expand = expansion(mult = c(0, 0.05))) +
  labs(x = NULL, y = "ShortBRED Count (RPKM-like)", fill = "CARD family") +
  theme_minimal(base_size = 12) +
  theme(axis.text.x = element_text(angle = 35, hjust = 1),
        panel.grid.major.x = element_blank())

Validated CARD ShortBRED result from six 20k-read public tutorial subsets

In this validation, only the anterior-nares and buccal subsets had nonzero family estimates: BRO-2 had two accepted hits, and patA had one. The buccal subset also had two patB hits but a family-level Count of zero, illustrating why both columns must be inspected. These sparse observations verify the workflow. They do not establish sample-level resistance, clinical susceptibility, gene host, or between-body-site differences.

9.5 CARD and MEGARes are different reference workflows

Reference	Data type	Appropriate use here
CARD v4.0.1 protein homolog models	Protein sequences	The lab generated ShortBRED protein-family markers and uses them in this practical.
MEGARes	Nucleotide sequences plus hierarchical AMR annotation	Use with AMR++ or a documented nucleotide-mapping workflow, including SNP confirmation where required. Do not substitute its FASTA for the ShortBRED protein marker file.

The current lab server keeps MEGARes v3.00 for its separate nucleotide workflow; the official MEGARes site now lists v4.0. Translating every MEGARes nucleotide record in six frames would mix coding sequences, rRNA targets, partial genes, and SNP-dependent models. It would not create a valid ShortBRED database.

ShortBRED with the present CARD homolog-marker set also does not test SNP-dependent variant models, rRNA mutations, gene knockouts, or resistance by overexpression. A marker hit describes sequence evidence, not phenotype.

10 Part 6 — k-mer taxonomy with Kraken2 and Bracken

10.1 Marker mapping and k-mer classification

MetaPhlAn and Kraken2 can analyze the same nucleotide reads, but their numbers answer different questions.

Method	Evidence searched	Main reported number	Denominator / reference space
MetaPhlAn	Selected clade-specific marker sequences	Taxon relative abundance and marker mapping evidence	Marker database and MetaPhlAn model
Kraken2	Minimizers selected from every read k-mer	Classified and unclassified reads plus taxonomic labels	All input reads searched against the chosen Kraken genome DB
Bracken	Kraken assignments plus a DB/read-length distribution	Abundance estimate at one chosen rank	Kraken report and the matching Bracken model

A higher Kraken classified percentage does not prove that Kraken is more accurate than MetaPhlAn. Database content, classification rules, and the denominator all differ. Compare biological patterns and provenance, not the percentages as if they were the same measurement.

10.2 One plaque sample live

The wrapper below uses the same plaque FASTA used in Part 1. The file is a single-end teaching subset, so the command deliberately does not use --paired. The comments inside 13_kraken_bracken_one.sh follow a fixed order: purpose, input, output, validation, Kraken2, Bracken, and output checks.

bash 03_scripts/13_kraken_bracken_one.sh \
  SRS014476-Supragingival_plaque \
  04_data/hmp_metaphlan_demo/SRS014476-Supragingival_plaque.fasta.gz

10.2.1 Kraken2 options used by the wrapper

Option	Meaning and teaching cue
--db	Exact versioned Standard-8 database. The older shared kraken2_db is incomplete and must not be substituted.
--threads	CPU threads used by this one sample. It is different from the number of concurrent samples.
--memory-mapping	Reads database pages through the operating system so concurrent jobs can share them.
--gzip-compressed	Reads the supplied .fasta.gz directly.
--minimum-hit-groups 2	Requires support from at least two hit groups. The current default is written explicitly for provenance.
--confidence 0	Baseline confidence threshold. It is not a probability and does not mean the result has 0% confidence.
--use-names	Writes taxon names into the per-read output for easier inspection.
--report / --output	Writes the hierarchical sample report and the per-read classification table.

The per-read Kraken file begins with C for classified or U for unclassified. The six-column .kreport summarizes clade counts through the taxonomy. Keep the default report format because Bracken cannot read an --use-mpa-style report.

The teaching Standard-8 database contains a human library, but a Kraken human label is not a substitute for validated host-read removal or a privacy screen. Host removal must be performed and audited earlier in the production pipeline with an appropriate, versioned human reference.

10.2.2 Bracken options used by the wrapper

Option	Meaning and teaching cue
-d	The same Kraken database used to classify the sample.
-i	The default six-column Kraken report.
-r 100	The bundled 100 bp read-length distribution, closest to this tutorial data.
-l S	Estimate abundance at species rank.
-t 10	Require at least ten reads assigned at or below a taxon before reporting it.
-o / -w	Write a tabular species result and a Kraken-style Bracken report.

The Bioconda Bracken launcher calls python internally. The wrapper therefore adds the pinned Kraken/Bracken environment to PATH for the Bracken command; this makes the result reproducible without asking students to activate conda.

Bracken recalculates percentages among the reported classification levels and does not include unclassified reads in its rewritten report. Therefore fraction_total_reads must not be read as a fraction of every raw input read. Keep the Kraken classified/unclassified counts beside the Bracken table.

10.3 Six samples and merged tables

Run the batch only after the one-sample output passes. THREADS controls one Kraken2 process; JOBS controls how many samples run at once. During class, show this code and use the precomputed six-sample tables below.

THREADS=2 JOBS=1 bash 03_scripts/14_kraken_bracken_batch.sh

source 03_scripts/00_config.sh
python3 03_scripts/15_merge_kraken_bracken.py \
  "$RESULTS_DIR/kraken_bracken" \
  "$RESULTS_DIR/merged/kraken_bracken"

kraken_summary_path <- file.path(
  expected_root,
  "merged",
  "kraken_bracken_summary.tsv"
)
bracken_species_path <- file.path(
  expected_root,
  "merged",
  "kraken_bracken_species_fraction.tsv"
)

if (!file.exists(kraken_summary_path) || !file.exists(bracken_species_path)) {
  cat("Validated Kraken2/Bracken tables have not been bundled yet.\n")
} else {
  kraken_summary <- read.delim(
    kraken_summary_path,
    check.names = FALSE,
    quote = "",
    comment.char = ""
  ) |>
    left_join(metadata |> select(sample_id, body_site), by = "sample_id") |>
    mutate(body_site = factor(body_site, levels = metadata$body_site))

  print(knitr::kable(
    kraken_summary |>
      select(body_site, total_reads, classified_reads,
             classified_percent, bracken_species_reported),
    digits = 2,
    caption = "Kraken2/Bracken six-sample workflow check"
  ))

  print(ggplot(kraken_summary, aes(body_site, classified_percent)) +
    geom_col(fill = "#4E79A7", width = 0.72) +
    scale_y_continuous(labels = label_percent(scale = 1),
                       expand = expansion(mult = c(0, 0.04))) +
    labs(x = NULL, y = "Kraken2 classified reads") +
    theme_minimal(base_size = 12) +
    theme(axis.text.x = element_text(angle = 35, hjust = 1),
          panel.grid.major.x = element_blank()))

  bracken_matrix <- read.delim(
    bracken_species_path,
    check.names = FALSE,
    quote = "",
    comment.char = ""
  )
  stopifnot(names(bracken_matrix)[1] == "species")

  bracken_long <- bracken_matrix |>
    pivot_longer(-species, names_to = "sample_id", values_to = "fraction") |>
    mutate(fraction = as.numeric(fraction)) |>
    left_join(metadata |> select(sample_id, body_site), by = "sample_id")

  top_bracken <- bracken_long |>
    group_by(species) |>
    summarise(mean_fraction = mean(fraction), .groups = "drop") |>
    slice_max(mean_fraction, n = 10, with_ties = FALSE) |>
    pull(species)

  bracken_plot <- bracken_long |>
    mutate(species_plot = if_else(species %in% top_bracken,
                                  species, "Other reported species")) |>
    group_by(body_site, species_plot) |>
    summarise(fraction = sum(fraction), .groups = "drop") |>
    mutate(body_site = factor(body_site, levels = metadata$body_site))

  print(ggplot(bracken_plot, aes(body_site, fraction, fill = species_plot)) +
    geom_col(width = 0.72) +
    scale_y_continuous(labels = label_percent(),
                       expand = expansion(mult = c(0, 0.04))) +
    labs(x = NULL, y = "Bracken fraction_total_reads",
         fill = "Species") +
    theme_minimal(base_size = 12) +
    theme(axis.text.x = element_text(angle = 35, hjust = 1),
          panel.grid.major.x = element_blank()))
}

## 
## 
## Table: Kraken2/Bracken six-sample workflow check
## 
## |body_site            | total_reads| classified_reads| classified_percent| bracken_species_reported|
## |:--------------------|-----------:|----------------:|------------------:|------------------------:|
## |stool                |       20000|             5587|              27.93|                       36|
## |anterior_nares       |       20257|             8737|              43.13|                       11|
## |tongue_dorsum        |       20000|             4024|              20.12|                       32|
## |buccal_mucosa        |       20000|             5962|              29.81|                       29|
## |supragingival_plaque |       20000|             3975|              19.88|                       40|
## |posterior_fornix     |       20000|             7284|              36.42|                       15|

Kraken2 classified-read percentage for the same six 20k-read tutorial subsets

Kraken2 classified-read percentage for the same six 20k-read tutorial subsets

Teaching questions:

1. Why can the Kraken classified percentage exceed a MetaPhlAn marker mapping percentage without indicating better accuracy?
2. What changes if the database release changes while the reads stay fixed?
3. Why must the Bracken read-length model match the post-QC reads?
4. Which step classifies reads, and which step estimates abundance at one rank?

11 Interpretation checklist

• UNCLASSIFIED: MetaPhlAn did not report a supported clade assignment under the current marker database/model.
• UNMAPPED: HUMAnN did not assign the sequence to a known gene family.
• UNINTEGRATED: gene-family signal was found but not reconstructed into a reported pathway.
• ShortBRED Count is a normalized family estimate; Hits is supporting raw read evidence. Zero in these tiny subsets does not prove biological absence.
• A CARD family hit is not an antimicrobial-susceptibility result and does not identify the organism carrying the gene.
• CARD protein markers and MEGARes nucleotide mapping answer related but non-interchangeable AMR questions.
• Kraken2 classified percentage and MetaPhlAn marker mapping are different metrics. Bracken is an abundance estimator based on Kraken output, not a new classifier.
• Standard-8 is a capped teaching database. Results must not be presented as a production full-database profile.
• Taxon abundance and pathway abundance are not one-to-one: genes can be shared, pathways can be incomplete, and functions can be contributed by several taxa.
• Tool version, database release, preprocessing, and input read set must all be fixed before comparing runs.

12 Exercises

1. Change JOBS=1 and JOBS=2. Which controls concurrency and which controls threads per sample?
2. Repeat the species plot with the top five species only.
3. Split a HUMAnN table into stratified and unstratified components.
4. Explain why an AMR marker hit is not the same as phenotypic resistance.
5. Explain why the MEGARes nucleotide FASTA cannot replace the CARD ShortBRED protein-marker FASTA.
6. Explain why these tiny subsets should not be used for publication-grade body-site comparisons.
7. Write a provenance block containing input checksum, tool version, database index, command, and completion status.
8. Explain k-mer, minimizer, and lowest common ancestor using one short DNA sequence.
9. Explain why the course script does not use --paired for these FASTA files.

13 Sources and version notes

• ShortBRED repository and official example files: https://github.com/biobakery/shortbred.
• ShortBRED method: Kaminski et al., PLoS Computational Biology (2015), https://doi.org/10.1371/journal.pcbi.1004557.
• CARD data and release history: https://card.mcmaster.ca/download/.
• MEGARes and AMR++: https://www.meglab.org/megares/ and https://github.com/Microbial-Ecology-Group/AMRplusplus.
• Kraken2 manual: https://github.com/DerrickWood/kraken2/wiki/Manual.
• Bracken documentation: https://github.com/jenniferlu717/Bracken.
• Official June 2026 prebuilt Kraken2/Bracken indexes: https://benlangmead.github.io/aws-indexes/k2.

ShortBRED is MIT-licensed. This course bundles only its official toy example. The CARD-derived production marker FASTA remains on the lab server; it is not republished with the public Rmd/HTML.

14 Session information

sessionInfo()

## R version 4.4.2 (2024-10-31)
## Platform: aarch64-apple-darwin20
## Running under: macOS 27.0.1
## 
## Matrix products: default
## BLAS:   /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/lib/libRblas.0.dylib 
## LAPACK: /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/lib/libRlapack.dylib;  LAPACK version 3.12.0
## 
## locale:
## [1] C.UTF-8/C.UTF-8/C.UTF-8/C/C.UTF-8/C.UTF-8
## 
## time zone: Asia/Seoul
## tzcode source: internal
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
## [1] scales_1.4.0  ggplot2_4.0.0 tidyr_1.3.1   dplyr_1.1.4  
## 
## loaded via a namespace (and not attached):
##  [1] vctrs_0.6.5        cli_3.6.5          knitr_1.50         rlang_1.1.6       
##  [5] xfun_0.53          purrr_1.1.0        generics_0.1.4     S7_0.2.0          
##  [9] jsonlite_2.0.0     labeling_0.4.3     glue_1.8.0         htmltools_0.5.8.1 
## [13] sass_0.4.10        rmarkdown_2.30     grid_4.4.2         evaluate_1.0.5    
## [17] jquerylib_0.1.4    tibble_3.3.0       fastmap_1.2.0      yaml_2.3.10       
## [21] lifecycle_1.0.4    compiler_4.4.2     RColorBrewer_1.1-3 pkgconfig_2.0.3   
## [25] farver_2.1.2       digest_0.6.37      viridisLite_0.4.2  R6_2.6.1          
## [29] tidyselect_1.2.1   pillar_1.11.1      magrittr_2.0.4     bslib_0.9.0       
## [33] withr_3.0.2        gtable_0.3.6       tools_4.4.2        cachem_1.1.0