BBMap and the Human Reference Problem: Why Your Database Choice Matters as Much as Your Tool — Day 3 of 5

🧬 Day 3 of 5 — Benchmarking Host-Read Removal for High-Host-Content Metagenomes Day 1 framed the benchmark. Day 2 showed that two well-known alignment tools — BMTagger and Bowtie2 — produce noticeably different removal rates even when both use GRCh38 as the reference. The natural next question was: how much of that difference is the algorithm, and how much is the reference?

Today’s post answers that directly. I ran BBMap twice — same samples, same mapping parameters, same tool — against two different human references: the standard GRCh38.p12 and a repeat-masked hg19. The only thing that changed between the two runs was the database.

The result was a 3–5 percentage point difference in host removal across the high-host samples.

This is the most important lesson in the entire benchmark series: the reference is part of the method. Writing “human reads were removed with BBMap” in a Methods section is not a complete description. Stating the reference genome, the masking strategy, and the key alignment parameters is what makes a host-removal workflow reproducible.

What BBMap is and why it is used here

BBMap is a splice-aware short- and long-read aligner developed at the Joint Genome Institute. For host removal it offers a useful combination: a k-mer-seeded hybrid mapping strategy that is highly sensitive, configurable alignment thresholds that let you tune the sensitivity/specificity balance, and output flags that directly separate unmapped (retained microbial) reads from mapped (human) reads in a single command.

The key output flags are:

• outu1 / outu2 — reads that did not map to the human reference → your retained microbial reads
• outm1 / outm2 — reads that did map → your discarded human reads

This direct separation in a single pass is operationally cleaner than tools that require a post-processing extraction step.

The benchmark design for Day 3

The same fastp-preprocessed paired-end reads used on Days 1 and 2 were passed to BBMap twice, with only the path= argument changed between runs:

This is the controlled comparison that isolates the reference effect from the algorithm effect. When the only variable between two runs is the database, any difference in output is a reference effect.

BBMap: full command

Both runs used these parameters:

/usr/bin/time -v -o ${SAMPLE}_bbmap_${REF_NAME}.time.txt \
  bbmap.sh \
    minid=0.95 \
    maxindel=3 \
    bwr=0.16 \
    bw=12 \
    quickmatch \
    fast \
    minhits=2 \
    qtrim=rl \
    trimq=10 \
    untrim \
    path=${REF_DIR} \
    in1=${SAMPLE}_R1.fastp.fastq.gz \
    in2=${SAMPLE}_R2.fastp.fastq.gz \
    outu1=${OUTDIR}/${SAMPLE}_nonhuman_R1.fastq.gz \
    outu2=${OUTDIR}/${SAMPLE}_nonhuman_R2.fastq.gz \
    outm1=${OUTDIR}/${SAMPLE}_human_R1.fastq.gz \
    outm2=${OUTDIR}/${SAMPLE}_human_R2.fastq.gz \
    threads=16 \
    -Xmx48g

Key parameters explained:

The -Xmx48g flag deserves attention. BBMap stores the reference index in Java heap memory. For a full human genome reference, this requires tens of gigabytes. If you run BBMap without setting this appropriately for your reference size, it will fail mid-run with an out-of-memory error — sometimes after hours of processing. Always set this flag explicitly and allocate enough memory in your SLURM job request.

Building the BBMap reference database

For GRCh38:

bbmap.sh ref=GRCh38.p12.fa path=GRCh38_BBMap_DB -Xmx48g

For masked hg19:

bbmap.sh ref=hg19_masked.fa path=masked_hg19_DB -Xmx48g

The path= directory stores the pre-built index. Building the index from a full human genome takes 30–60 minutes and substantial memory, but it only needs to be done once. Store the index in a shared project location so all samples use the same database.

What masking actually does

Before looking at the results, it is worth understanding what a masked reference is and why someone would use one.

The human genome contains large amounts of:

• Repetitive sequence — transposable elements, satellite repeats, tandem repeats
• Low-complexity regions — homopolymers, simple sequence repeats
• Conserved sequence — regions similar to genes found in bacteria, archaea, or other eukaryotes

When a short read maps to one of these regions, the alignment may be ambiguous: the read matches the human reference, but it might equally match a microbial sequence. An aggressive mapper will classify it as human and remove it. A more conservative approach would leave it unclassified.

A repeat-masked reference replaces these problematic regions with N characters before indexing. A read that would have mapped to a repetitive or low-complexity human region can no longer map — because that region has been masked out. The mapper therefore retains reads it would otherwise have classified as human.

The trade-off is explicit:

Standard reference (GRCh38):
  Higher sensitivity for human reads
  Risk: some ambiguous/microbial reads classified as human

Masked reference:
  Lower sensitivity for human reads in repetitive regions
  Benefit: reduced false-positive removal of microbial reads

Neither is universally correct. The right choice depends on what matters most for your biological question.

Benchmark results

Across the five samples:

The masked reference consistently removed fewer reads. The difference is largest in high-host samples (Samples 1, 2, and 5) and smallest in the low-host sample (Sample 4) — which makes intuitive sense. In high-host samples, a larger fraction of reads are genuinely human, so the effect of masking out repetitive regions is proportionally larger.

The tempting but wrong conclusion

Looking at this table, a natural reaction is:

“GRCh38 removes more reads, so it is better.”

This is the wrong conclusion — or at least, it is an incomplete one. The additional reads removed by GRCh38 could represent:

1. True human reads that the masked reference missed — in which case GRCh38 is more sensitive and the additional removals are correct
2. Ambiguous reads that map to repetitive regions shared between human and microbial genomes — in which case both references are removing something real, but the GRCh38 run is removing more of the microbial co-mapping reads
3. Microbial false positives — reads from vaginal community members that share short regions of sequence similarity with unmasked repetitive elements

The removal percentage alone cannot distinguish between these three possibilities. That requires the read-level discordance analysis coming on Day 5, where I map discordant reads back to taxonomic databases and ask what organisms they actually come from.

Computational performance

Approximate mean performance across all five samples:

Two things stand out here.

First, BBMap required substantially more RAM than any other tool in this benchmark. BMTagger used ~8 GB, Bowtie2 used ~4 GB, and BBMap needed ~45 GB. This is the most practically significant limitation of BBMap for host removal: on a shared HPC cluster, allocating 45–50 GB per sample limits how many simultaneous jobs you can run on a given node. At scale — 100 samples — this becomes a scheduling bottleneck.

Second, the masked reference was ~27% faster than GRCh38. Since the masked reference has fewer indexable positions (the N-masked regions are ignored), the mapper processes fewer potential alignment targets per read. Faster runtime with lower sensitivity in repetitive regions — a consistent package.

Why the reference is part of the method

This benchmark illustrates a principle that applies well beyond host removal:

The output of a bioinformatics tool is a joint product of the algorithm, the database, the parameters, and the preprocessing. Naming only the tool is an incomplete description.

For host removal specifically, a reproducible Methods section should include:

✓ Tool name and version
✓ Human reference genome (GRCh38, hg19, T2T, etc.)
✓ Whether the reference was masked, and how
✓ Key alignment parameters (identity threshold, minhits, etc.)
✓ How paired-end reads were handled
✓ How unmapped reads were retained

A Methods section that says only “Human reads were removed using BBMap” does not contain enough information for another researcher to reproduce the analysis — or to understand why their results differ if they attempt to.

A complete description would read more like:

Human reads were removed using BBMap (v39.x) mapped against GRCh38.p12, requiring ≥95% nucleotide identity (minid=0.95) and a minimum of 2 k-mer seed hits (minhits=2). Unmapped read pairs were retained for downstream analysis.

That is a reproducible methods statement.

A framework for thinking about host removal

The benchmark has now tested four configurations across three tools. A useful way to organize them:

Host-removal result =
    f(algorithm, reference, parameters, preprocessing)

To isolate algorithm effects:
    → keep reference and parameters constant, vary the tool

To isolate reference effects:
    → keep algorithm and parameters constant, vary the database
      (this is what Day 3 tested)

To evaluate end-to-end performance:
    → use each tool's recommended reference and parameters,
      compare realistic production configurations

All three comparisons answer different questions. A benchmark that only runs each tool with its default settings on a single reference is measuring something, but it is conflating algorithm, reference, and parameter effects in a way that makes it hard to interpret differences.

What comes tomorrow

Day 4 introduces Cleanifier — the newest method in this benchmark and the one that most changed my thinking about the tradeoffs involved.

In the benchmark, Cleanifier produced host-removal percentages close to BMTagger while running in approximately 1.5 minutes per sample. That is not a typo. Where BBMap takes 77 minutes and requires 45 GB of RAM, Cleanifier finished the same sample in under two minutes.

That result was interesting enough that I did not stop at the runtime table. I asked the more important question: what exactly are the reads that Cleanifier removes that the alignment-based tools miss — and vice versa? That read-level investigation became the most analytically interesting part of the benchmark, and it is the focus of Day 4.