Vignette: ZymoBIOMICS Community Standard¶

1. About¶

ZymoBIOMICS® Microbial Community Standard (catalog D6300) is a commercial mock community of eight bacteria and two yeasts, mixed by the vendor to a known composition and sold specifically to expose bias in microbiome/metagenomics workflows — extraction, library prep, sequencing, and bioinformatics alike. Because the "right answer" is known in advance, it's a useful stand-in for a real sample when you want to sanity-check a pipeline rather than a biological hypothesis.

This vignette runs metamorph's combined workflow — read-based classification and assembly/binning/genome-level classification together, in one invocation — against public sequencing data generated from this standard, and asks a single question twice, once per method: taking everything a run produces — not any one tool's output in isolation — how close does metamorph's picture of this community land next to the vendor's defined ground truth, and do the two independent methods agree with each other? It's a notebook-style walk from raw data to a real answer, in the same spirit as the walkthrough but focused on interpreting output rather than touring --workflow modes.

2. Vendor ground truth¶

The standard's defined composition (Table 1 of the D6300 instruction manual, Ver. 1.3.0 — linked from the product page) is given in two units that matter for different kinds of sequencing:

¹ Reclassified as Bacillus spizizenii in 2020; both names refer to the same strain (NRRL B-354) in this standard. Lactobacillus fermentum underwent an analogous 2020 genus-level split (into Limosilactobacillus and several others); metamorph's own classifiers report both of this standard's organisms under their current GTDB names, not the vendor's. Section 5 covers exactly how that's handled, and treats it as one instance of a single general rule, not a Zymo-specific fix. ² The vendor's current lot ships Cryptococcus deneoformans (NRRL Y-2534) under the same "neoformans" species-complex label; treated as equivalent here.

The manual also ships its own real benchmark of this exact bias question (Figure 1): four DNA extraction methods run against the same standard and sequenced identically, compared to the theoretical composition. Only the vendor's own kit reproduces the target composition; the other three (including the Human Microbiome Project's fecal protocol) are visibly skewed toward the Gram-negative, easy-to-lyse organisms. Section 5 recreates that same chart with both of metamorph's results added to it.

3. Data used in this run¶

Public data don't come pre-packaged as ENA's ERR16012159–ERR16012165, seven replicates deposited under the "Automated Enzymatic Prep" library preparation in ENA study PRJEB105240, all sequenced from the same defined ZymoBIOMICS standard. This run uses the full-depth FASTQs (not a subsample) — roughly 2–2.3 GB per mate per replicate, ~27 GB total across all 7.

This run also passes the standard's own two yeast reference genomes (Saccharomyces cerevisiae, Cryptococcus neoformans/deneoformans, from the vendor's reference genome bundle) to metamorph's --host-genome flag. metamorph has no yeast-calling capability, so without this flag any read belonging to either yeast that happens to survive default (human-only) dehosting is left for the bacterial classifiers to deal with. --host-genome is a generic pipeline feature, not specific to this standard or to yeasts — it screens reads against any additional FASTA genome(s) given to it, strictly on top of the default human screen, and here it's pointed at this sample's own known non-bacterial content so that signal is removed before classification rather than left for a bacterial classifier to (mis)handle.

Sample sheet (DNA-only, 7 samples):

DNA
/path/to/ERR16012159_R1.fastq.gz
/path/to/ERR16012159_R2.fastq.gz
/path/to/ERR16012160_R1.fastq.gz
/path/to/ERR16012160_R2.fastq.gz
...
/path/to/ERR16012165_R1.fastq.gz
/path/to/ERR16012165_R2.fastq.gz

./metamorph run \
    --samplesheet samplesheet.tsv \
    --output output/ \
    --workflow combined \
    --host-genome /path/to/ZymoBIOMICS.STD.refseq.v3/Genomes/yeast.fasta.gz \
    --sif-cache /data/OpenOmics/SIFs

One command, one sample sheet, no per-sample babysitting, and no second invocation to get assembly-based results — --workflow combined runs both the read-based stack and the assembly/binning/genome-classification stack in the same run. The rest of this page is everything that one invocation produced.

4. What a single run produces¶

--workflow combined runs the full pre-screening stack (QC, host removal, Centrifuger) on every sample, layers HUMAnN3/MetaPhlAn4 community profiling and cohort-level diversity analysis on top (same as read-based), and assembles, bins, dereplicates, and taxonomically classifies recovered genomes across the whole cohort. For these 7 samples, that's a real, immediately usable set of outputs sitting in metagenome_results/:

Sections 5–7 build the rest of this page entirely out of these real files — no additional tools, no manual re-analysis, just combining outputs metamorph already wrote to disk.

4.1 Where every number on this page comes from¶

Every statistic in Sections 5–7 traces back to exactly one or two files metamorph wrote under metagenome_results/, run through one specific program, with one documented transformation applied before it reaches a table on this page. Nothing below is tuned per-organism — the same code path runs for all 8 target bacteria in both methods.

4.2 Dehosting: how much was screened out, per replicate¶

Every read pair in this run passes through two sequential dehosting filters before it ever reaches a classifier: the default human screen (hg38), then the injected yeast screen (Section 3). Each filter now writes its own plain-text summary alongside the dehosted FASTQs — trimmed_reads/<replicate>/<replicate>_hg38_dehost_summary.txt for the default screen (always written), and trimmed_reads/<replicate>/<replicate>_injected_host_dehost_summary.txt for the injected-genome screen (written whenever --host-genome is used). Each lists the reference and program used, the exact retention criterion, and that replicate's input/removed/retained read-pair counts and percentages — this is the file to open for a direct answer to "how much did dehosting remove for this sample." (These two files didn't exist yet when this run itself executed; both dehosting rules were extended to write them afterward, and the ones in this output directory were backfilled from the same read-pair counts reported below, which were already independently verified against the surviving screen_injected_host_dna Slurm logs. Every run from here on writes them natively, no backfill needed.)

Read pairs removed at each stage, per replicate (removed = present going in, absent coming out of that stage; the remainder is retained and carried forward to the next stage):

The hg38 screen removes almost nothing (0.11–0.17% of reads, mean 0.14%) — consistent with this being a defined mock community with no human-derived material, not a clinical or environmental sample where a much larger human fraction would be expected. The injected yeast screen removes an order of magnitude more (1.56–1.93%, mean 1.79%) — a real, consistently-sized fraction of reads in every replicate that would otherwise have been left for a bacterial classifier to (mis)handle, had --host-genome not been pointed at this standard's own yeast genomes. Section 9 covers the direct consequence of removing that fraction: read-based classification can no longer detect either yeast once their reads are screened out this early.

5. Two independent composition estimates vs. vendor target¶

metamorph's combined workflow produces two genuinely independent pictures of this community's composition, built from different evidence:

• Read-based: Centrifuger's per-replicate species-level abundance (GTDB+RefSeq k-mer classification of every read).
• Assembly-based: each replicate's reads mapped back against the cohort-wide dereplicated genome set, converted to a genome-copy-equivalent abundance (recruited bases ÷ that genome's own size — the same size-correction logic the vendor's own manual uses to turn mass into genome copies, see Section 2's note).

Read-based (Centrifuger, all 7 replicates, both rules applied; renormalized across the 8 bacteria to 100% per replicate):

Assembly-based (GTDB-Tk + genome-copy-normalized MAG mapping, all 7 replicates, both rules applied; renormalized across the 8 bacteria to 100% per replicate):

No yeasts appear in this estimate: neither yeast assembled into a recoverable bin at this coverage (expected — see Section 9), and screening their reads out at the dehosting stage (Section 3) means there's no leftover yeast signal for assembly to pick up either.

Both estimates plug into the manual's own Figure 1 — the vendor's benchmark of four DNA extraction methods against the theoretical composition — putting both of metamorph's results in the same frame as the vendor's own published comparison:

Both metamorph bars are recognizably the same community as the vendor's five, with one deviation both estimates agree on and share with each other more than with any vendor bar: Lactobacillus fermentum visibly larger than any vendor bar, Staphylococcus aureus visibly smaller. Enterococcus faecalis is also smaller than target in both, though less dramatically in the assembly-based bar than the read-based one. Section 6 checks whether these gaps are real and whether the two methods agree with each other.

6. Is the gap real, and do the two methods agree?¶

Is each deviation bigger than replicate-to-replicate noise? A one-sample t-test per species (the 7 replicate values against the target, df=6) says yes for every species in the read-based estimate, and for 7 of 8 species in the assembly-based estimate:

Read-based classification puts every one of the 8 bacteria statistically significantly off the genome-copy target, reproducibly across 7 independent replicates. The assembly-based estimate agrees for 7 of the 8 — but not for Bacillus subtilis, whose assembly-based mean (10.67%) lands within replicate-to-replicate noise of its 10.41% target (p = 0.21). That single exception is itself informative: it shows the assembly-based method isn't just uniformly biased in the same direction as read-based classification — for at least this one organism, it recovers the target closely enough that the two methods' disagreement with each other (not just with the vendor) is the more interesting fact.

Genome-level classification resolves the B. subtilis/spizizenii naming question directly: GTDB-Tk classifies the recovered genome as Bacillus spizizenii with no cross-tool translation needed, the same current name read-based classification already reports (Section 2, footnote 1). The GTDB placeholder-suffix split does occur at the genome level, though — the recovered Listeria genome is classified s__Listeria monocytogenes_B by GTDB-Tk itself, the same split read classifiers show — so the generic collapse rule from Section 5 is applied identically to both read- and genome-level calls, not just read-level ones.

Do the two independent methods agree with each other, even where they both disagree with the vendor? Bray–Curtis dissimilarity (8-bacteria basis) between the two mean compositions is 0.155 — smaller than either method's distance to the vendor target (0.301 read-based, 0.181 assembly-based). The assembly-based estimate's remaining distance from the vendor target is concentrated in just three species: L. fermentum (over-called, 1.55×), S. aureus, and E. faecalis (both under-called, 0.23× and 0.58×) together account for roughly 83% of that total dissimilarity — the other five species are each within about 1–2 percentage points of their target. The same two species dominate in both methods: L. fermentum over-called by both (2.11× read-based, 1.55× assembly-based), S. aureus under-called by both (0.19× and 0.23×). Two methods built from different evidence (reads vs. recovered genomes) landing closer to each other than either lands to the vendor's defined composition is a real signal: whatever is causing these two shared deviations is upstream of classification — it survives a complete change of classification strategy — which points toward extraction or library prep rather than a database or algorithm artifact in either tool.

Would averaging in more replicates close the gap, or just narrow the uncertainty around it? Bray–Curtis dissimilarity between the cumulative mean composition (replicates 1, then 1–2, then 1–3, … 1–7) and the full 10-taxon vendor Genome Copy % target (both yeasts counted), for both methods:

Neither line trends toward 0 as replicates accumulate — both are essentially flat from replicate 1 onward (read-based settles at ~0.30–0.31, assembly-based at ~0.17–0.19). Both estimates converge — replicate-to-replicate variance is genuinely small (Section 8) — but both converge to a point well short of the vendor's defined composition, not to it. More replicates sharpen confidence in where that point is for each method; on this evidence, they don't move either point closer to the vendor's target, and the assembly-based point sits substantially closer to the vendor target than the read-based one across the full replicate range.

Does Escherichia coli's shortfall point in a consistent direction between the two methods? No — the two methods disagree about which direction the bias runs for this one organism: read-based puts it at 0.28× target (under-called), assembly-based puts it at 1.17× target (over-called). Section 7 adds a third, fully independent method to this specific question — mapping directly against the vendor's own reference genomes, bypassing both Centrifuger's database and metamorph's own recovered genomes — and it lands on the same side as assembly-based (1.28× target, over-called), not read-based. That's evidence, not yet proof, that this particular disagreement is a read-based/Centrifuger-specific artifact rather than a real DNA-abundance effect.

7. A third, independent check: mapping directly against the vendor's reference genomes¶

Sections 5–6 build two estimates from what --workflow combined itself produces: Centrifuger's database-driven read classification, and metamorph's own assembled/dereplicated genomes. Both still share one thing in common — they classify reads against something built from the data itself (a general-purpose database, or genomes assembled from these exact samples). The D6300 manual's own Appendix A describes a third option that depends on neither:

The manual's own suggested mitigation — "a mapping tool that can choose to ignore reads that map to more than one genome" — is exactly the unambiguous-reads criterion metamorph's own assembly-based estimate already uses (Section 4.1). This cross-check reuses that same criterion and the same coverage-normalization logic, with one variable changed: the reference set.

All three estimates, per species (renormalized across the 8 bacteria; error bars are SD across the 7 replicates):

The two genome-level methods — metamorph's own recovered MAGs and the vendor's reference genomes, built from completely different evidence — agree with each other far more than either agrees with read-based classification. Bray–Curtis dissimilarity (8-bacteria basis): assembly-based vs. direct reference-mapping is 0.039 — roughly a quarter of assembly-based's distance to read-based (0.155) and about a fifth of read-based's distance to the vendor target (0.301). Direct reference-mapping's own distance to the vendor target (0.175) is close to assembly-based's (0.181) — both genome-level methods land in essentially the same place relative to the vendor's defined composition.

That tight two-way agreement makes the handful of real disagreements easy to spot:

• Escherichia coli* and Salmonella enterica are the two cases where read-based disagrees in direction with both genome-level methods. *E. coli: read-based 0.28× target (under-called) vs. 1.17× (assembly-based) and 1.28× (direct reference-mapping) — both genome-level methods over-call it. S. enterica: read-based 0.84× (under-called) vs. 1.23× and 1.34× — the same pattern, not previously flagged on this page because the assembly-based/read-based comparison alone couldn't distinguish "both non-read-based methods agree" from "coincidence." With a third, fully independent method landing on the same side as assembly-based for both organisms, the more likely explanation is a read-based/Centrifuger-specific classification effect for these two Enterobacterales — not a real DNA-abundance effect, and not an assembly-specific artifact either, since direct reference-mapping never touches metamorph's own assembled genomes at all.
• Bacillus subtilis* is a borderline third case. Read-based over-calls it clearly (1.31×, *p = 0.00003 in Section 6). Both genome-level methods put it much closer to target — assembly-based not significantly different (Section 6), and direct reference-mapping only marginally so (1.04×, p = 0.056, just short of the conventional 0.05 cutoff). Weaker than the E. coli/S. enterica pattern, but the same shape: read-based diverges, the two genome-level methods agree with each other and with the vendor target.
• The other five species — P. aeruginosa, L. fermentum, Listeria, E. faecalis, S. aureus — show all three methods agreeing on direction, even where the magnitude differs (most visibly L. fermentum: 2.11× read-based vs. 1.55–1.37× for the two genome-level methods — still a real over-call by every method, just a smaller one once classification is bypassed entirely). This is the strongest version yet of Section 6's "upstream of classification" argument: a bias three independently-built methods agree on is harder to attribute to any one tool's database or algorithm.
• Both yeasts remain undetectable here too (S. cerevisiae ~0.0001%, C. neoformans ~0.00003% of total mapped coverage) — consistent with Section 9's point that --host-genome removes their reads before any classifier or mapper, including this one, ever sees them. This isn't a new finding; it's a confirmation that the near-zero yeast signal elsewhere on this page reflects reads actually being gone, not a database that happens not to contain these yeasts (it does — they're in this exact reference set — and still almost nothing maps).

7.1 How closely related are these genomes? (Mash distance)¶

Section 7's scope note above raises a specific question: is E. coli/S. enterica's unusual behavior plausibly related to genomic similarity between the two of them, as Appendix A's own conserved-region warning describes? That's a question about the reference genomes themselves, independent of any sequencing data — answerable directly with Mash, which estimates genome-to-genome distance from k-mer sketches rather than full alignment.

Mash's distance estimate depends on sketch size (-s): too few hashes and a real-but-small similarity between distant genomes can be missed entirely (zero shared hashes just saturates the estimate at a distance of 1, which reads as "maximally dissimilar" even when it really just means "undetected at this resolution"). A sketch-size grid (s = 10,000 to 500,000, all 10×10 genome pairs at each) was run to check for exactly that — the figure below uses the largest, s = 500,000, since more hashes give the most stable distance estimate and the least of that saturation artifact; the full grid is saved at enzymatic_community_data/ZymoBIOMICS.STD.refseq.v3/Genomes/mash_grid_results.tsv if the other sketch sizes are of interest.

One pair stands out by a wide margin: E. coli/S. enterica at 0.16 — the closest relationship anywhere in this heatmap, well below the next-closest bacterial pair (E. faecalis/L. monocytogenes, 0.20) and roughly a third of the typical cross-genus bacterial distance here (~0.33–0.45). That's exactly consistent with their shared Enterobacteriaceae ancestry, and it means that of all 28 bacterial pairs in this standard, E. coli and S. enterica are each other's single most plausible source of exactly the kind of shared, highly-conserved sequence Appendix A's note warns about (Section 7). It doesn't by itself explain the specific pattern found in Section 7 — Appendix A's mechanism predicts conserved regions pulling reads from a high-abundance organism toward a low-abundance one (an asymmetric over/under pair), while Section 7 found both organisms over-called by the genome-level methods relative to read-based, not one over- and one under-called — but it's real, independent evidence that these two specific organisms are genomically closer to each other than any other pair in the standard, which is a prerequisite for that kind of cross-mapping effect even if it isn't full proof of it here.

The two yeasts' distances are the least trustworthy numbers in this heatmap: several yeast-vs-bacterium cells are pinned at exactly 1.0 even at this largest sketch size (e.g. L. fermentum vs. either yeast), which given how distantly related eukaryotes and bacteria are is at least plausible as a real result — but it's indistinguishable from the sketch still being too small to detect a genuine, very weak signal. Nothing in this vignette's composition estimates depends on the yeast distances, so this is noted for completeness, not because it affects any other finding on this page.

8. Inter-replicate consistency¶

dna_humann_diversity_calculation (only runs when a sample sheet has more than one sample) gives a real look at replicate-to-replicate variance, not just a point estimate:

9. Caveats and what's still open¶

• The E. coli/S. enterica read-based divergence (Section 6, Section 7). Read-based puts both organisms on the under-called side (0.28× and 0.84× target); the two genome-level methods (assembly-based and direct reference-mapping, Section 7) both put them over-called instead (1.17–1.28× and 1.23–1.34×). The agreement between two independently-built genome-level methods against a disagreeing read-based classifier is evidence this is a read-based/Centrifuger-specific effect rather than a real DNA-abundance effect — but it isn't a diagnosed root cause. Still open: why Centrifuger's classification specifically diverges for these two Enterobacterales and not the other six organisms.
• Yeast reads are screened out before classification by design — this run evaluates bacterial detection, not yeast detection. Removing reads that map to this standard's own yeast genomes before they reach a classifier is the explicit purpose of pointing --host-genome at them here (Section 3, Section 4.2): a clean read pool for the 8 target bacteria, with no risk of fungal reads leaking into a bacterial call. Across the 7 replicates, 1.56–1.93% of reads that had already survived default (hg38) dehosting mapped to the injected yeast genomes and were removed at this stage (Section 4.2). The expected, direct consequence is that Centrifuger reports essentially no Saccharomyces cerevisiae signal in this run (~0.0003%) — those reads were removed before they could ever reach it. Section 5's composition estimates are built and renormalized over the 8 bacteria only; the yeast rows shown there are informational, not a measure of this run's yeast-detection accuracy, which was never the question this run is set up to answer.

10. Reproducing this vignette¶

The 7 "Automated Enzymatic Prep" runs are ERR16012159 through ERR16012165 in ENA study PRJEB105240. Fetch metadata (including the exact FASTQ paths) for any of them with:

curl -sL "https://www.ebi.ac.uk/ena/portal/api/filereport?accession=ERR16012159&result=read_run&fields=run_accession,fastq_ftp,library_name&format=tsv"

Download each FASTQ in full (no byte-range truncation). For the --host-genome yeast references, download and extract the vendor's reference genome bundle and point --host-genome at its combined yeast FASTA (Genomes/yeast.fasta.gz, or any FASTA containing both yeast genomes). Run with --workflow combined as shown in Section 3 to reproduce every number and figure on this page.