Validation, accuracy and limitations¶

This document covers both halves of the tool. Part A is metabolite gene capacity, validated against synthetic communities with known gene content. Part B is species composition and disease-pattern similarity, validated against public cohorts that carry disease labels.

Jump to: Part A headline figures · Part B headline figures · the specificity matrix · limitations

Everything about how well this test performs, and where it falls short, is in this one file. Nothing here is estimated — the figures come from openbiota validate, and the raw output is in validation_results.json.

Reproduce with:

openbiota validate                    # ~5 min with genomes cached
openbiota cohort                      # ~10 min, rebuilds the reference ranges

Headline figures¶

At the operating point the tool actually reports (confirmed tier):


Sensitivity	97.0% — genes present that were found
Specificity	100.0% — genes absent that were correctly called absent
Precision	100.0% — detections that were real
False positives	0 of 62 true negatives
False negatives	1 of 33 true positives
Accuracy	98.9%
Detection limit	carriers found down to 0.2% of cells
Quantitative calibration	slope 1.13, R² 0.72

Every pathway individually achieved 100% specificity and 100% precision.

Part B — profile similarity¶

A different kind of validation, because the ground truth is different. For gene capacity, truth is "does this genome carry this gene", which a synthetic community establishes exactly. For disease similarity there is no such truth — only whether the score separates labelled cases from labelled controls.

Only colorectal cancer has both a cross-cohort-validated microbiome association and public cohorts with disease labels attached, so it is the one profile on which the engine can be tested end to end. That is why it ships as a validation control rather than as a health readout.

The colorectal cancer result¶

258 labelled cases against 250 within-study controls (YachidaS 2019, via curatedMetagenomicData), scored through the identical engine, with the reference distribution built from the control arm only:


AUC	0.62 (95% CI 0.57–0.67)
Leave-one-study-out AUC	0.62
Ceiling for these six features	0.67, measured directly
Published cross-cohort performance	~0.80, from learned models over hundreds of features

Three things follow, and all three are stated in the report:

The engine recovers the signal. All five oral-origin markers are enriched in the correct direction — F. nucleatum in 18.6% of cases versus 4.4% of controls, P. micra 37.2% versus 14.8%, and similarly for P. stomatis, G. morbillorum and S. moorei. The confidence interval excludes chance.

0.62 is near the ceiling for six prespecified markers, not a shortfall in the engine. Counting how many of the five markers a sample carries gives AUC 0.671; summing their log abundances gives 0.674. The published 0.80+ comes from LASSO models fitted over the full profile, which is a different method that this tool deliberately does not use — naively transferred learned coefficients lose both accuracy and disease specificity.

The profile is not condition-specific, which is the most important finding in this document. See the matrix below.

Cross-profile specificity: the most informative table here¶

Every profile scored against every labelled condition. A profile should score highest against its own condition and near 0.5 against the others.

Profile	vs CRC	vs IBD	vs T2D
`crc`	0.62	0.76	0.52
`longcovid`	0.53	0.80	0.59
`mecfs`	0.56	0.49	0.60

crc separates IBD better (0.76) than the condition it was built for (0.62). Its oral-origin markers are a general-disturbance signal. Oral-gut translocation happens in periodontal disease, with reduced gastric acid, and in inflammatory bowel disease, not only in malignancy.
longcovid scores 0.80 against IBD — higher than anything else in the table. R. gnavus enrichment and F. prausnitzii depletion are classic IBD findings. Its own discrimination is unmeasured, because the accessions that would allow measuring it publish reads without a sample-to-outcome key.
mecfs is flat and near chance everywhere (0.49–0.60). Its features are weak discriminators of anything, which is consistent with the modest species-level effects both source cohorts reported. Notably it does not spike on IBD, so it is at least not a relabelled inflammation signal.

This is why GMWI2 accompanies every profile score. A community that is generally disturbed will resemble several of these patterns at once, and the anchor is what tells you whether that is what is happening.

What is not measured¶

No AUC is printed for mecfs or longcovid, and none should be. Several Long COVID accessions publish reads without a sample-to-outcome key; PRJNA878603 (ME/CFS) has published code and a key resources table but no disease-labelled cohort in curatedMetagenomicData. Where labels cannot be obtained, the honest statement is that discrimination is unmeasured — not a proxy figure.

Reference cohort stability¶

The species reference cohort holds 3,027 healthy adult stool metagenomes from 22 studies, median depth 38.2 million reads, all processed through one pipeline. Bootstrapping sub-cohorts of increasing size shows how quickly a percentile estimate settles:

Cohort size	SD of the median estimate
25	0.81
50	0.54
100	0.33
200	0.24
400	0.16
800	0.10

This is what justifies the sample-count target empirically rather than by assertion, and it is where the confidence grade's threshold of 200 comes from.

Host filtering, verified with a positive control¶

Human reads are removed against GRCh38 plus PhiX before taxonomic profiling. On sample A02 this removed 0 of 8,029,909 pairs, and a 0.00% rate is exactly what a silently broken index also looks like — so it was checked rather than assumed.

Two thousand synthetic read pairs cut directly from the GRCh38 sequence align at 100%, confirming the index resolves. Aligning the real sample finds 2 of 400,000 mates, a 0.00% overall rate. The sample genuinely carries negligible human DNA, which is what a well-prepared stool library should look like.

Reproduce with:

bowtie2 -p 10 --very-fast -x refs/host/host -1 <human_reads_1.fq> \
        -2 <human_reads_2.fq> -S /dev/null

Because host filtering was a no-op for this sample, the species profile is identical with and without it — a useful consistency check in its own right, since the two runs went through different code paths.

What was not validated, and why¶

Synthetic taxonomic communities. The plan was to extend the synthetic harness — which builds communities from reference genomes with known gene content — to spike profile taxa at defined abundances and check that the profiler recovers them. This was not done, because the colorectal validation above supersedes it: it tests the same machine end to end against clinical ground truth rather than against a simulator, and a simulator that produced reads from the same reference genomes MetaPhlAn's markers were drawn from would mostly measure that circularity. What synthetic communities would add that the labelled cohorts do not is a detection limit for individual species, which remains unmeasured.

Per-feature reproduction of the ME/CFS cohort. PRJNA878603 publishes code and a key resources table, so the reported per-feature directions could be re-derived. This was not done. The consequence is that the ME/CFS feature weights rest on the papers' reported effects rather than on independently reproduced ones, and the profile's discrimination is unmeasured either way.

Comparison against the published assembly-and-HMM pipeline for gene capacity remains open; the protocol is at the end of this document.

Two bugs this validation caught¶

Both produced plausible-looking numbers, which is what makes them worth recording.

A depth-dependent zero floor. The CLR's zero-replacement value was originally derived per sample as half that sample's smallest detected abundance — the textbook rule. That makes the floor depth-dependent: a deeply sequenced sample detects rarer taxa, gets a lower floor, and maps its zeros further down than a shallow one. For a taxon absent in most samples, that difference is the reported signal. The floor is now a fixed constant (1e-4 percent), applied identically to sample and reference.

Absence treated as a number. For a taxon absent from both the sample and most of the reference, the CLR values still differ, because they carry the geometric-mean denominator, which varies with how many other taxa each sample happened to detect. Comparing those values compares sequencing incidentals. On sample A02 it placed five absent colorectal markers at the 74th percentile and produced a spurious 74th-percentile similarity score. Absence is now a category: a sample lacking the taxon sits at the mid-rank of the reference samples that also lack it, which puts those five markers at the 48th percentile and the profile score at the 23rd, where it belongs.

Both are pinned as regression tests in tests/test_scoring.py.

How this was measured¶

The hard part of validating a metagenomic screen is knowing the right answer. This validation constructs the sample so that the answer is known.

1. Ground truth from complete genomes. 24 reference genomes were downloaded from NCBI. For each, its own annotated proteome was searched at full length against the pathway reference sets, requiring a genuine ortholog: ≥70% identity across ≥70% of the reference and ≥60% of the query, plus the active-site residue criterion where a pathway defines one. That is the published method's decision procedure, applied where it works properly — on whole proteins rather than 50-residue fragments.

2. Synthetic samples. Reads were simulated from those same genomes at known cell abundances, with matched read lengths (100–151 bp), NovaSeq-style binned quality values, a 0.2% substitution error rate, and mate ids identical between R1 and R2 exactly as real NovaSeq output produces. Because the community composition is set, the true number of gene copies per 100 genomes is known exactly rather than estimated.

3. Measurement. The pipeline was run on the simulated reads with no knowledge of any of the above.

Five experiments, 3,000,000 read pairs each:

Experiment	Community	Measures
`balanced`	all 24 genomes, equal abundance	sensitivity, specificity
`negative`	6 organisms verified to lack the target genes	the false-positive rate, directly
`spike_0p05`	carrier at 5% of cells	detection at moderate abundance
`spike_0p01`	carrier at 1%	detection at low abundance
`spike_0p002`	carrier at 0.2%	the detection limit

Genomes used span the known carriers for every pathway — C. scindens, C. hylemonae and P. hiranonis for bile acids, K. pneumoniae and P. mirabilis for TMA, F. prausnitzii and R. intestinalis for butyrate, E. lenta, S. pasteurianus and L. plantarum for imidazole propionate, C. difficile for p-cresol, C. sporogenes for indole-3-propionate — plus six abundant gut organisms that carry none of them. The full list with accessions is in validation_results.json.

Operating points, and why the thresholds are what they are¶

Criterion	Sensitivity	Specificity	Precision	False positives
Any matching fragment	100.0%	93.5%	89.2%	4
≥2 fragments	100.0%	96.8%	94.3%	2
≥2 fragments and ≥75% identity	97.0%	100.0%	100.0%	0

The tool reports the third row as confirmed and everything else detected as provisional. The thresholds were not chosen in advance — they are where the measured false-positive count reaches zero.

Two things made this possible:

A single read is not evidence. Both single-fragment false positives disappear at two fragments, and no true call is affected: the weakest true positive had 64 fragments.
Identity separates the classes. Every false positive fell below 73% mean translated identity, while true positives clustered at 81–99% with a median of 98%. A read from an organism that genuinely carries the gene matches its reference closely; matches in the 60s are landing on distant relatives within the same protein fold.

The one call given up at the confirmed tier is the urdA carrier at 0.2% of cells — the hardest case in the set. That trade is the right way round: a reported detection is one you can act on, and borderline signals are still shown, labelled provisional.

Quantitative accuracy¶

Detection is one question; whether the number is right is another.

Overall: slope 1.13, intercept 0.50, R² 0.72 across 37 comparisons of reported against known copies per 100 genomes.

Per gene, as median reported ÷ true across confirmed calls:

Gene	Bias	Gene	Bias
`bai:baiF`	1.00×	`cutc:cutC`	0.84×
`bai:baiB`	1.02×	`indole:tnaA`	0.69×
`bai:baiE`	1.15×	`ipa:fldBC`	0.41×
`bai:baiCD`	1.16×	`pcresol:hpdB`	0.42×
`butyrate:buk`	1.17×	`propionate:pduP`	2.32×
`cutc:cntA`	1.17×	`propionate:pct`	2.42×
`bai:baiA`, `bai:baiH`	1.21×	`urda:urdA`	2.43×
`ipa:fldH`	1.29×
`bsh:bsh`	1.31×
`butyrate:but`	1.66×

Most genes land within about 1.7× of truth. Four are off by roughly 2.4× — two over (pct, pduP, urdA) and two under (fldBC, hpdB, the latter pair having the smallest reference sets, so reads from carriers not represented in the set are lost).

This is why results are reported as percentiles. The bias is systematic per gene — the same gene measured the same way in two samples is off by the same factor in both, so it cancels when comparing samples to each other or to the reference cohort. It does not cancel in an absolute number. Read the percentile; treat the absolute value as accurate to within roughly a factor of two.

The active-site residue filter, measured¶

urdA is the pathway whose specificity is hardest, because fumarate reductase flavoprotein is a close relative. The published method separates them by requiring tyrosine or methionine at FAD-site residue 373 of the full-length protein, which a 151 bp read mostly cannot reach.

The pipeline's reference-side substitute — checking whether the matched reference itself carries an acceptable residue — was tested against residue-filtered truth:

Figure	Slope vs truth	R²
Headline count	4.31×	0.991
Residue-consistent subset	1.65×	1.000

The filter cuts overcounting by a factor of 2.6, with essentially perfect linearity. This is why both figures are reported: the headline is an upper bound, the residue-consistent subset is the tighter estimate, and the truth sits between them.

Two supporting checks:

Check	Result
Anchor self-test — does position 373 map onto itself in `Q8CVD0`?	Passes on every build; failure aborts the run
Curated Swiss-Prot urocanate reductases carrying Y/M at their own mapped position	6 of 6 (100%), at correctly distinct positions 373, 594, 394, 598, 594, 408
Automatically annotated (TrEMBL) entries carrying Y/M	~47% of 1,537 mapped
Read-side vs reference-side estimate on the real sample	~477 vs 514 fragments — 7% apart, from two checks that share no evidence

Was the sample sequenced deeply enough?¶

The pipeline reads every fragment in the input — nothing is subsampled. The separate question is whether the sequencer sampled the stool deeply enough that a rarer organism could not have been missed. openbiota depth-check answers it by screening the same sample at increasing depths and watching for a plateau.

A pathway counts as converged when its value changes between half and full depth by less than the Poisson counting error implied by its own fragment count — i.e. by less than its own noise floor, leaving no recoverable signal.

Measured on the 8,029,909-read-pair sample in this repository:

Pathway	502k	1.0M	2.0M	4.0M	8.0M	Change	Noise	Settled
`bai`	7.60	8.61	8.58	8.42	8.20	2.7%	±3.0%	yes
`bsh`	67.49	59.08	59.26	58.01	55.62	4.3%	±3.0%	yes
`butyrate`	87.80	87.69	85.98	83.27	82.77	0.6%	±2.3%	yes
`cutc`	3.39	3.57	3.63	3.87	4.05	4.3%	±7.0%	yes
`indole`	9.09	10.43	11.50	11.40	11.55	1.3%	±5.5%	yes
`ipa`	12.75	11.17	10.68	9.87	9.53	3.5%	±5.6%	yes
`pcresol`	9.61	9.98	10.23	8.44	8.61	2.0%	±4.7%	yes
`propionate`	32.07	36.19	35.35	32.30	32.99	2.1%	±2.7%	yes
`urda`	35.69	36.86	35.26	35.30	35.07	0.6%	±2.7%	yes

All nine pathways had plateaued well before full depth — most are within a few percent of their final value by 500,000 read pairs, which is 6% of the data. Deeper sequencing would not change these results.

This also demonstrates the depth-independence that the rpoB normalisation is supposed to provide, on real data rather than by argument: a 16-fold change in depth moves the reported values by single-digit percentages, while the raw fragment counts scale linearly with depth as they must.

Two consequences worth drawing out:

Deeper sequencing buys precision, not new detections. Fragment counts scale with depth, so counting error shrinks as 1/√n. At 8M read pairs the headline pathways carry 2-3% counting error, which is already well below the per-gene systematic bias measured above.
A shallower sample is still usable, provided it is sampled uniformly. This is why the reference cohort can be built from 2M-pair subsamples rather than whole runs: with uniform sampling the expected values match and only the per-sample noise is wider. A prefix is not a uniform sample — see the reference-cohort section below for the measurement.

The reference cohort¶

Percentiles come from 91 adult stool shotgun metagenomes from ENA study PRJNA1271016. ENA describes the study as "MARS shotgun metagenomics — Gut microbiome and AD study in MARS cohort" (University of Wisconsin–Madison, USA, 2023): an older-adult, Alzheimer's-disease study population. It is the closest available technical match to this pipeline's inputs — same body site, same library strategy, same NovaSeq X platform family — and the report names it on every page that quotes a percentile, because a reader comparing a child against it needs to know that.


Samples	91 (one of 92 excluded: 580 `rpoB` fragments, a low-bacterial specimen)
Reads per sample	2,000,000 pairs, sampled uniformly across the whole run
Sampling	30 evenly spaced spot windows per run, seeded (`20260909`), reproducible
`rpoB` depth	median 2,219, range 1,820–2,674
Provenance	`refs/reference_ranges.json` records the 91 run accessions, the seed and the study title; `openbiota cohort --runs` rebuilds it exactly

Observed ranges, copies per 100 bacterial genomes:

Pathway	5th	25th	Median	75th	95th
`b12`	61.06	74.78	87.76	101.32	119.17
`bai`	1.08	1.95	2.66	4.06	6.57
`bcaa`	136.47	165.37	184.13	197.45	215.17
`bsh`	51.02	70.43	87.20	99.93	122.72
`butyrate`	17.02	22.91	27.20	32.00	39.30
`cutc`	1.09	1.75	2.47	3.87	5.57
`gaba`	5.83	11.35	18.85	24.74	36.61
`indole`	5.27	9.51	14.46	20.45	30.28
`ipa`	1.33	3.09	4.65	7.15	11.17
`pcresol`	2.72	4.92	6.52	8.87	15.29
`propionate`	12.10	15.71	19.72	25.62	38.48
`urda`	18.00	29.05	37.73	47.39	62.40

Why the reads are sampled across the run, and what that changed¶

Until this build the cohort was made from the first 600,000 read pairs of each of 34 runs — a file prefix — on the argument that every figure is a ratio to rpoB, so depth divides out. That was asserted, not measured. Measured, on one sample screened three ways at the same 600k budget against its full-depth (8.03M) answer:

600k arm	median error vs full depth	panels >10% off
first 600k read pairs (the old cohort's method)	19.6%	16 of 24
a uniform random 600k	11.6%	13 of 24

and the random-arm deviations are pure counting noise: median |z| = 0.87 against Poisson expectation 1.0, no panel beyond 3σ, sign split 16 high / 7 low. Two conclusions follow, and they decided the design:

The head of a FASTQ is a biased sample of it. Counting by flowcell tile, the first 600k pairs of SAMPLE2_A02 carry 28.98 panel fragments per rpoB fragment against 27.26 for the rest of the file. The first tiles have a different quality profile, and at a fixed identity threshold that changes how many fragments each panel accepts. Uniform sampling removes it.
Once sampling is uniform, depth is only precision. There is no directional drift with depth, so the cohort does not have to be depth-matched to the query. It has to be deep enough that its own counting noise is small against the biological spread it describes — at 2M pairs that noise is ~5% per panel against relative interquartile ranges of 0.3–0.9.

The rebuild did what the diagnosis predicted, and one thing it did not:

	old (34 × first-600k)	new (91 × random-2M)
Panel medians	—	0.929× the old, every panel lower: the prefix cohort was ~7% inflated
Relative IQR (spread)	0.54	0.52 — only 4% narrower. The old spread was already dominated by biology, not noise; the precision gain I projected from noise arithmetic was overstated
Our five samples, median percentile	64.1	73.0 — higher, because the inflated scale came down. Correct direction; it unmasks the population mismatch below rather than reducing it

Sampling itself is fastq-dump spot-range access into the SRA archive, which seeks rather than transfers: a window at spot 18,000,000 costs the same as one at spot 1, so an unbiased 2M-pair sample of a 32M-pair run takes about four minutes rather than the fifteen it takes to stream and discard the whole file. openbiota cohort --sampling prefix is retained for cheap builds and is labelled as biased in the output.

Depth still matters per panel¶

Low-count panels have not converged at 2M and the reader should know which: openbiota depth-check screens a sample at increasing depths and names them. Comparing SAMPLE2_A02 at full depth against a 600k cohort moved its median panel percentile by only −1.3 points — but 7 of 25 panels moved by more than 10, and histamine, ipa, butyrate, oxalate, pcresol are the usual suspects.

The population, honestly¶

Four of the five samples in this repository are children aged 8–18 and the fifth is 48; the reference population is older adults in an Alzheimer's study. Their percentiles run high (median 73, not 50) and no amount of sampling or depth changes that, because it is a real difference in who is being compared. A paediatric functional cohort would; it is the next piece of work.

Cross-platform check: Tiny Health on the same FASTQ files¶

Every validation above is against simulated reads, reference genomes or this pipeline's own cohort. The one external check on real samples is that Tiny Health, a CLIA-certified shotgun metagenomics service, ran the same five stool specimens. Two independent pipelines, one set of DNA.

Tiny Health is not treated as ground truth. Where the two disagree, the question is which is right, and the answer comes from the sequence evidence.

Reproduce with python -m tools.compare_tinyhealth against your own Tiny Health report and run. The comparison is a tool, not a fixture: nothing from these samples is pinned in the test suite.

What is comparable, and what it shows¶

Absolute values are not comparable — Tiny Health reports RPKM against their reference population, this pipeline reports copies per 100 bacterial genomes against PRJNA1271016. What is comparable is where a reading sits on each platform's own scale, and relative species abundance, which is unitless.

Check	Result
Like-for-like functional readings, 5 samples	89 pairs
Same band (low / typical / high) on both	45 (51%)
Within one band	83 (93%)
Opposite at a band edge (see below)	3
Opposite	3 (3%), all three explained below
Their top-20 species also detected here (`SAMPLE2_A02`)	17 / 17 (100%)
Abundance rank correlation on those species	Spearman ρ = 0.775
Median abundance ratio between platforms	1.41x

The 1.41x median abundance ratio is the expected size of a MetaPhlAn-versus- GTDB disagreement on the same reads: the two use different marker sets and different reference catalogues, and Tiny Health's figures run consistently higher because their denominator excludes more unclassified sequence.

Band-edge straddles¶

A three-band comparison is fragile exactly at the lines. Three pairs read as "opposite" only because Tiny Health's value sits 0.5–1.5% below their own cutoff while ours sits 1–3 percentile points above our 75th line: SAMPLE1_A01 propionate (theirs 546 vs cutoff 549; ours 76.3rd), SAMPLE1_A01 vitamin B9 (268.7 vs 271; 77.4th) and SAMPLE6_A06 hydrogen sulfide index (9.01 vs 9.15; 77.5th). For propionate and B9 the rank order of the three same-template samples is identical on both platforms — the two agree about the values and differ about where a band starts. The comparison tool reports these separately and applies the rank-order check to every straddle.

The three opposite calls, resolved¶

Branched-chain amino acids, SAMPLE1_A01 and SAMPLE2_A02 — reference population, not measurement. On the three samples that share a Tiny Health report template, BCAA capacity ranks identically on both platforms (theirs 2507 / 3064 / 3025 RPKM; ours 209 / 235 / 210 copies per 100 genomes). The two pipelines measure the same thing and agree on the ordering; they disagree only about whether that value is high, because they compare against different populations. Neither is wrong. Worth knowing: this panel is close to saturation — ilvC and ilvD are near-universal in gut bacteria, so the whole cohort spans 123–234 and a percentile inside that range carries little information.

GABA production, SAMPLE3_A03 — a Tiny Health false negative. They report 0.0 rpkm. This pipeline finds 1,940 gadB fragments at 100% median identity, best-matching Bacteroides clarus (592), Phocaeicola dorei (182), Bacteroides stercoris (178) and Phocaeicola vulgatus (146), in a sample that is 34% Bacteroidaceae. Bacteroides is the dominant genus of gut GABA producers (Strandwitz et al., Nature Microbiology 4:396–403, 2019, doi:10.1038/s41564-018-0307-3). A third of the community cannot carry the gene and produce a zero; this pipeline counts the gene whatever organism carries it, so it cannot report that combination.

One pair excluded as not like-for-like¶

Histamine. Tiny Health's v4.8.3 metric is "Histamine-producing species": the summed abundance of a curated species list. This pipeline's is hdcA gene capacity, whatever organism carries it. For SAMPLE1_A01 and SAMPLE3_A03 they report 0.0% while this pipeline finds hdcA at 100% identity in Eggerthella lenta (14 and 66 fragments) and Bacteroides fragilis — neither on their list. This is the structural argument for gene-based over taxon-based functional calls: a curated list can only find carriers someone already thought of. The two are counted separately rather than scored as a disagreement, because they are measurements of different quantities.

A reproducibility problem on their side¶

Tiny Health's RPKM values are not comparable across their own report versions. Between template v4.8.3 (samples SAMPLE1, SAMPLE3, SAMPLE6) and v5.x (SAMPLE2, SAMPLE4), their median value shifted by 2.41x — and in the same direction for all 12 comparable panels (range 1.74–3.02x). Twelve of twelve in one direction is p ≈ 0.0002 under a sign test.

That could be biology if the two groups of people simply differed. They do not: this pipeline measured the same specimens and puts the two groups 1.12x apart, with the ratios scattered on both sides of 1.0 (0.79–1.65). A consistent multiplier on one platform and scatter on the other is what a silent renormalisation looks like, not a cohort difference. The practical consequence is that a Tiny Health customer comparing a 2026 report to an earlier one would see every functional number fall by roughly half with no biological change. test_our_units_are_stable_across_pipeline_revisions asserts this pipeline does not do the same.

What this check does not cover¶

The comparison touches 89 of the ~1,390 quantities this pipeline reports per sample. Disease-pattern resemblance, biological age, resistance and virulence determinants, skin research panels and the species-level percentiles have no Tiny Health counterpart, so nothing here validates them. They are covered, to the extent they are covered at all, by the sections above and by docs/VALIDATION.md's per-module limitations.

Limitations¶

Grouped by what they affect.

What the measurement is¶

This measures gene capacity, not metabolite concentration. Carrying the gene means the bacteria can make the metabolite; how much they actually make depends on substrate availability, transit time and competition. Measuring the metabolites requires mass spectrometry on blood or stool.
Gene carriage predicts metabolite level weakly. In the source study of 294 samples, of all urdA-positive genomes only Streptococcus pasteurianus abundance was significantly associated with plasma imidazole propionate. This is a property of the biology, not of this implementation.
No clinical thresholds exist. No value in any output corresponds to a diagnosis, because no such threshold has been published for any of these pathways. The percentile says where you sit in a population; it does not say what that means for your health.

Accuracy boundaries¶

Absolute values are accurate to roughly a factor of two, with four genes at ~2.4×. Percentiles and longitudinal changes are much more reliable than absolute numbers.
The detection limit is about 0.2% of cells for a well-referenced gene. Rarer carriers may be missed.
urdA loses sensitivity at the confirmed tier (75% versus 100% permissive), because its hardest case sits right at the identity threshold.
Small reference sets bound sensitivity. bai:baiCD and bai:baiH have four references each, propionate:lcdA five, pcresol:hpdB 77. A low value for these is weak evidence of absence, and the few-refs flag marks them.

Profile similarity boundaries¶

A percentile is not a probability. It says the sample is more concordant with a published pattern than that fraction of a matched reference group. Nothing more.
The associations are mostly not condition-specific. The specificity matrix above shows the colorectal pattern separating IBD better than colorectal cancer, and the Long COVID pattern scoring 0.80 against IBD. Read every profile score beside the dysbiosis anchor.
Two of the three profiles have unmeasured discrimination. No labelled cohort exists for ME/CFS or Long COVID that this tool can reach, so no AUC is printed for them.
No exposure-matched contrast exists for ME/CFS or Long COVID. The right comparison is people with the same exposure who did not develop the syndrome; the published cohorts use healthy or uninfected controls instead, so syndrome-specific signal cannot be separated from post-exposure signal. Both profiles declare this and take a confidence penalty for it.
The combined 50/35/15 module split is an engineering starting point, not an evidence-derived weighting. Report the module scores separately, which the output does.
Profiler versions must match. Species abundances from MetaPhlAn 3 and MetaPhlAn 4 are not interchangeable: version 4 reorganised the taxonomy into SGBs and renamed genera. The sample is profiled with the same version as the reference cohort, both are recorded in results.json, and a mismatch downgrades confidence.
Long COVID scores chronic state only. Acute future-risk prediction from stool at 0–2 months post-infection is a different task, and features from it are excluded rather than reused.
Effect sizes in this literature are small. The largest independent prospective Long COVID study (799 outpatients) reached species-level AUROC 0.62 against 0.72 for clinical variables alone; a 2,561-participant community study found no association with illness duration. The frequently quoted AUC of 0.96 classifies a study's own internally defined enterotypes, not an externally validated diagnosis.

Method boundaries¶

No per-genome attribution for gene capacity. Read-level search cannot say which organism carries a gene. Named organisms in the gene-capacity output are the nearest reference sequence, not identifications. Species in Part B are measured identifications, from marker genes.
Synteny is unavailable. cutC/cutD gene neighbourhood is a strong specificity signal in assembly-based work and cannot be used here.
Taxonomy is coarse. The community profile is reliable at phylum level only, and is bounded by the reviewed rpoB reference set.
The succinate route to propionate is not screened, so propionate capacity is systematically understated. See PANELS.md.

Validation boundaries¶

Depth is not a limiting factor for this sample, but it would be for a much shallower one. Below roughly 500,000 read pairs the rarer pathways (bai, cutc, ipa) fall under the fragment counts needed for a confirmed call.
Validation used simulated reads. Simulation reproduces read length, quality binning, error rate and mate structure, but not every artefact of a real library — chimeras, uneven coverage, PCR duplicates, strain-level variation within a species.
24 genomes is not the gut. A real sample contains hundreds of species, many unsequenced. Some will carry homologs the validation set does not represent.
The reference cohort is 33 samples of mixed clinical status. ENA metadata for PRJNA1271016 does not distinguish case from control, so the percentiles describe adult stool generally rather than a screened-healthy population. Widen it with openbiota cohort --samples 100.
Not compared against the published pipeline on shared samples. The strongest possible validation would run both methods on the same real metagenomes and compare per-genome calls, using the MAGs in PRJNA1469852. That has not been done. The protocol is below.

Remaining work: comparison against the published method¶

The one validation not yet done, and how to do it.

1. Fetch the MAGs. PRJNA1469852 holds the metagenome-assembled genomes the source study called urdA-positive and negative, per sample.

2. Reproduce their ground truth. Call genes on each MAG with Prodigal, score against their ImP-producer HMM, keep hits above 60% of the lowest-scoring model sequence, align with Clustal Omega and require Y or M at residue 373.

3. Run this pipeline over the matching raw metagenomes from PRJNA1271016.

for r1 in fastq/*_1.fastq.gz; do
  openbiota run --r1 "$r1" --r2 "${r1/_1/_2}" --out results/comparison
done

4. Build the confusion matrix against their per-sample calls, sweeping a threshold over both copies_per_100_genomes and copies_per_100_genomes_residue_consistent. The difference between those two ROC curves quantifies what the reference-side residue filter buys on real data rather than synthetic.

5. If plasma imidazole propionate is available, correlate it against the pipeline's output. Expect a weak relationship — that is the correct result, not a failure, and it is the point about gene carriage made above.

Until this is done, accuracy is characterised against synthetic ground truth and against curated reference proteins, which is what the figures at the top of this file represent.