A sample of known answers: rebuilding a classic mouse gut study from its raw reads

A 2012 study of the mouse gut after weaning, rebuilt on AlkhemYst from the sequencer's own files.

---

The mice

In 2012 Patrick Schloss's laboratory at the University of Michigan followed a group of mice from weaning through their first year and collected their stool every day. The question was simple to state and hard to answer: once a young animal stops drinking milk, how long does its gut community take to settle, and what does settling look like?

The answer, in Gut Microbes [1], came in four parts. There were two community states, one in the first ten days after weaning and one from about two weeks on. The two states held the same bacterial populations in different proportions. The later state was the more stable one. And that stability did not come from the community becoming richer or more diverse.

The tube

A year later the same laboratory published the sequencing method behind those samples [2]. Kozich and colleagues showed how to read the V4 region of the 16S rRNA gene on an Illumina MiSeq with two barcodes per sample, and they did something that makes their data unusually useful to anyone checking a pipeline: they sequenced a mock community alongside the mice. A mock community is a tube of DNA from bacteria chosen in advance, twenty-one strains in this case, mixed and sequenced as if it were a sample. It is the one sample in the run where the right answer is known.

Twenty of those samples have been public since 2013: nineteen days from one mouse, nine of the first ten after weaning and ten around day 145, and the mock. They are the teaching set of mothur, Schloss's own software, and of DADA2 [3], whose tutorial runs on exactly these files. A great many people have learned amplicon analysis on them, which is the reason we chose them. A familiar dataset is the right first test of something new: when the answer is known, a different answer is a finding about the software, not the biology.

What we did

AlkhemYst now has a set of components for amplicon sequencing, from raw reads to predicted function, and this was their first full run. Fourteen steps on a canvas, no code written, every step from the standard library: a read quality report, primer trimming, quality filtering, denoising into exact sequence variants with DADA2, the mock community check, taxonomy against SILVA [4], a phylogenetic tree, a BIOM file for other tools, composition, alpha and beta diversity, differential abundance, a functional prediction with PICRUSt2 [5], and the pathways that prediction implies.

The fourteen steps became 1,259 lines of R, run on our servers in the EU with the network switched off, as every run is.

The canvas noticed one thing before it did anything else. Primer trimming found the primer in none of the reads. That is correct, and a reason to read the method rather than the defaults: in Kozich's design the custom sequencing primers contain the 16S primers themselves, so every read begins just after them and there is nothing to remove. The step said so and passed the reads on untouched. A step that had done as it was told would have discarded every read in the study.

The known answer

The mock community is where to start, because it is the only sample whose contents are not a matter of opinion.

DADA2 found 19 sequence variants in the mock sample. All 19 are exact, base-for-base matches to the reference sequences of the strains in the tube. None of them is unexpected: no contaminant, no chimera, no error that survived denoising. Twenty of the twenty-one strains are accounted for. One variant matches two of them, Staphylococcus aureus and S. epidermidis, whose V4 regions are identical, so no method reading this region could separate them. The strain missing is Propionibacterium acnes.

The DADA2 tutorial, run on the same reads, reports 20 variants in the mock, all 20 exact, and 232 variants in the study overall. We found 19 and 229. The difference has a plain cause: the tutorial truncates every read at fixed lengths, 240 and 160 bases, where our filtering step chooses the lengths from the quality profile and chose 243 and 159. Close, not identical, and for a reason that can be named, which is what a reproduction should look like.

The four claims

Then the mouse. One caution first: the paper's claims rest on several animals, and the public set holds one of them, so what follows tests each claim within a single mouse's time course. We grouped its samples as the paper does, the first ten days after weaning against days 141 to 150, and asked each of the four claims in turn.

Two states. They are there, and strongly. The early and late samples separate on every distance we computed: PERMANOVA [6] puts the time point behind 43 percent of the variation in Bray-Curtis distance (p = 0.001), and behind 32 to 40 percent on the other three.

The same populations, in different proportions. This one we cannot confirm cleanly. Membership alone, which bacteria are present regardless of how many, also separates the two states (Jaccard, 32 percent of the variation), less than proportions do but clearly. The likely reason is resolution. The paper counted operational taxonomic units, sequences grouped at 97 percent identity; we counted exact variants, which split some of those units in two. A finer ruler finds differences in membership that a coarser one folds together. We are reporting it as a difference rather than explaining it away.

The later state is more stable. Yes. Measured as the average distance of each sample from its group's centre, the early samples are spread more widely than the late ones: 0.224 against 0.151 in Bray-Curtis (p = 0.005), 0.104 against 0.063 in weighted UniFrac (p = 0.029). On membership alone the two spreads do not differ. What steadies is the balance of the community, not its roster, which fits the paper's second claim even where our first test of it did not.

Stability without more diversity. Also yes. Richness, Shannon diversity, Simpson diversity and evenness do not differ between the two states (p between 0.12 and 0.68). Phylogenetic diversity sits on the boundary, p = 0.05, lower in the late samples if anything. The late community is steadier without being richer.

Six methods, one question

Which bacteria actually changed between the two states is the question a reader most wants answered, and the one where the method matters most. We ran it six ways on the same table, at genus level, on the same 52 genera.

The answers range from five genera to twenty-four. Nothing here is broken. Nearing and colleagues compared fourteen such methods on thirty-eight datasets in 2022 [7] and found the same spread, with ALDEx2 and ANCOM among the most consistent, and recommended reading several methods together rather than trusting any one.

Read together, these six say something clear. All of them agree on three genera, and agree on the direction: Lactobacillus and Turicibacter are more abundant in the late state, and a group of the Lachnospiraceae (UCG-001) less. The two methods Nearing found most consistent, ALDEx2 and ANCOM-BC2, name five genera each and share three of them, the same three. Everything beyond those three depends on the method, which is worth knowing before writing any of it into a results section. A platform that offers one method by default and hides the others makes that check hard to do. Having all six one click apart makes it the obvious thing to do.

The part that was actually hard

Almost none of it was the analysis.

The first run through the platform stopped at the taxonomy step. The classifier loads SILVA's reference into memory, and the sandbox every run executes in capped memory at four gigabytes, a limit set for tabular work long before the platform read sequencing data. Every test of the new steps had run outside that sandbox, which is why none of them saw it. Sequencing canvases now run with the memory and time they need, and every other canvas keeps the limits it had.

The functional prediction needed its own change. PICRUSt2 places each variant into a tree of twenty thousand reference genomes, and its default placement tool wanted more memory than the server has. It now uses the other placement method PICRUSt2 offers, which divides the tree into parts.

And the BIOM file, the format other microbiome tools read, came out of the R library that writes it in a shape QIIME 2's reader refuses. We now write it ourselves and check it with the reader that matters.

And one we found while writing this article. Run on the table the first canvas exported, ANCOM-BC2 named seven genera; run inside the first canvas, five. The counts were identical. What differed was the order of the bacteria in the table, which a phylogenetic tree imposes when it is loaded with them, and the step that groups variants into genera named and merged the unclassified ones in that order. The differential abundance step now sorts the table before it tests, and both routes give the same answer. The same data has to give the same result however it arrives, or a reproduction proves nothing.

None of these is interesting on its own. They are here because an account in which everything worked the first time would be a worse guide to what bringing a new kind of data to a platform actually costs.

Why this story

There is nothing to correct here, and we did not go looking for anything. The paper's claims mostly hold on a different method six years later, the one that does not hold cleanly has a reason we can state, and the mock community came back exact.

The story is the tube. Somebody decided, in 2013, to spend sequencing capacity on a sample that tells you nothing about mice and everything about the measurement, and then to publish it beside the mice. That decision is why this reproduction can say "correct" rather than "plausible". A mock community costs one sample per run. If your laboratory sequences amplicons and does not include one, it is the cheapest check there is.

---

Evidence

Papers. Schloss PD, Schubert AM, Zackular JP, Iverson KD, Young VB, Petrosino JF (2012). Stabilization of the murine gut microbiome following weaning. Gut Microbes 3(4):383-393. doi 10.4161/gmic.21008.

Kozich JJ, Westcott SL, Baxter NT, Highlander SK, Schloss PD (2013). Development of a dual-index sequencing strategy and curation pipeline for analyzing amplicon sequence data on the MiSeq Illumina sequencing platform. Applied and Environmental Microbiology 79(17):5112-5120. doi 10.1128/AEM.01043-13.

Data. The MiSeq SOP set distributed by the mothur project (mothur.org/wiki/MiSeqSOP): 20 samples (19 days from one mouse, and the mock community), 40 FASTQ files, 2 x 250 bases, V4, with the reference sequences of the HMP mock community (HMPMOCK.v35.fasta). Time points from the file names: days 0 to 9 after weaning as early, days 141 to 150 as late.

Pipeline. pipe_LyAxasBFT0kf, fourteen library components, language R. Read Quality Report, Primer Trimming (Cutadapt), Filter & Trim Reads, ASV Inference (DADA2), Mock Community Check, Taxonomy Assignment (SILVA 138.1), Phylogenetic Tree, BIOM Export, Taxonomy Composition, Alpha Diversity, Beta Diversity, Differential Abundance (ANCOM-BC2), Functional Prediction (PICRUSt2), Pathway Abundance.

Runs. run_D1PqKm0og5iL, 14 of 14 steps, 1,259 lines of generated R, every number above except the table of six methods. run_QOq7b0idhn4P (pipeline pipe_BTvrsxGLbjZR), the six methods on the BIOM file, sequences and tree the first run exported, imported with Import Feature Table.

Measured. 152,360 read pairs; 91.5 percent passed filtering; 229 exact sequence variants; 80.7 percent of reads reached the final chimera-free table; 100 percent classified to phylum, 55 percent to genus. Mock community: 19 variants, 19 exact matches, 20 of 21 strains, none unexpected.

Credit. The study, the method and the data are their authors', in work led by Patrick Schloss's laboratory at the University of Michigan. The mock community reference is the Human Microbiome Project's. This article describes a second analysis of their data and claims nothing further.

References

[1] Schloss PD et al. Stabilization of the murine gut microbiome following weaning. Gut Microbes 2012;3(4):383-393. doi:10.4161/gmic.21008

[2] Kozich JJ et al. Development of a dual-index sequencing strategy and curation pipeline for analyzing amplicon sequence data on the MiSeq Illumina sequencing platform. Applied and Environmental Microbiology 2013;79(17):5112-5120. doi:10.1128/AEM.01043-13

[3] Callahan BJ, McMurdie PJ, Rosen MJ, Han AW, Johnson AJA, Holmes SP. DADA2: High-resolution sample inference from Illumina amplicon data. Nature Methods 2016;13(7):581-583. doi:10.1038/nmeth.3869

[4] Quast C, Pruesse E, Yilmaz P, Gerken J, Schweer T, Yarza P, Peplies J, Glöckner FO. The SILVA ribosomal RNA gene database project: improved data processing and web-based tools. Nucleic Acids Research 2013;41(D1):D590-D596. doi:10.1093/nar/gks1219

[5] Douglas GM, Maffei VJ, Zaneveld JR, Yurgel SN, Brown JR, Taylor CM, Huttenhower C, Langille MGI. PICRUSt2 for prediction of metagenome functions. Nature Biotechnology 2020;38(6):685-688. doi:10.1038/s41587-020-0548-6

[6] Anderson MJ. A new method for non-parametric multivariate analysis of variance. Austral Ecology 2001;26(1):32-46. doi:10.1111/j.1442-9993.2001.01070.pp.x

[7] Nearing JT, Douglas GM, Hayes MG, MacDonald J, Desai DK, Allward N, Jones CMA, Wright RJ, Dhanani AS, Comeau AM, Langille MGI. Microbiome differential abundance methods produce different results across 38 datasets. Nature Communications 2022;13:342. doi:10.1038/s41467-022-28034-z