henhaobio
Sign in

bioRxivMicrobiologyPreprint

Accurate detection of metagenomic strain-level associations using average nucleotide identity with StrainSpy

Sudaraka MallawaarachchiiD, Kshitij TandoniD, Nikhil RajaniD, Vanessa R. MarcelinoiD, Shahneen Kaur SandhuiD, Sammy BedouiiD, Danielle J. IngleiD, Ashray GunjuriD, Gerry Q. Tonkin‐HilliD

PDFDOICodeData
In-depth full-text digest available · Jump to digest ↓

Abstract

Genetic variation among microbial strains of the same species can profoundly influence their phenotypes, ecological functions, and impacts on human health. Traditionally, the relative abundance of a species has been used to identify associations between the microbiome and disease. However, this approach overlooks intra-species genetic variation and is susceptible to spurious correlations arising from the compositional nature of abundance data and microbial load. Fast, k-mer-based algorithms can now accurately estimate strain-level Average Nucleotide Identity (ANI) in metagenomes. Despite its value as an orthogonal metric for strain-level analysis, existing ANI association methods are limited to case-control designs and rely on arbitrary identity thresholds. To address this, we developed StrainSpy, a statistical algorithm that jointly captures strain variations and replacement across longitudinal and multi-cohort studies without predefined thresholds. In simulations, StrainSpy achieved >98% sensitivity versus 79% for comparable approaches at similar specificity. Analysis of a longitudinal antibiotic study in 12 healthy adults revealed novel strain-level associations across multiple species, including reduced strain-level diversity despite persistence of the species, and showed that a previously reported association with E. coli was driven by reduced microbial load. Across 3,414 metagenomes from 17 colorectal cancer cohorts, StrainSpy identified novel strain-level associations in E. coli accompanied by enrichment of adhesin genes in addition to recovering well established associations. Importantly, machine learning models built on StrainSpy-selected features matched the predictive accuracy of abundance-based approaches, whilst avoiding the confounding effects of microbial load. StrainSpy is available as an R package at github.com/gtonkinhill/strainspy.

In-depth summary

Claude Opus 5.5 · xhigh effort · 8 Nov 2026 · from the full text of version 2 (3 Oct 2026) · CC BY 4.0

StrainSpy fits zero-inflated beta mixed models to containment ANI, reaching >98% simulated sensitivity versus ~79% for logistic regression and detecting colorectal-cancer E. coli strain shifts.

Abundance-based microbiome association tests ignore within-species variation and are confounded by compositionality and microbial load, while existing cANI methods suit only case-control designs and need arbitrary identity thresholds. StrainSpy treats cANI as a threshold-free outcome in mixed models, supporting longitudinal and multi-cohort designs without per-species phylogenies.

Key results

  1. ShownIn signal-implantation simulations on healthy gut metagenomes, the zero-inflated beta (ZiB) model in StrainSpy detected strain identity shifts and depletion more sensitively than threshold-based logistic regression while keeping specificity near perfect; ordered beta failed on depletion.

    200 Zeevi et al. metagenomes, 100 strains modified. At cANI changes >1%, ZiB and ordered beta >98% sensitivity vs ~79% for logistic regression (specificity drop <0.05%); all models detected 0.5% changes and kept specificity >99.99%; ordered beta sensitivity <25% at depletions up to 75%. · StrainSpy provides robust results in simulated…; Using signal implantation to simulate strain…; Figure legends

  2. SuggestedReanalysis of a longitudinal antibiotic cohort revealed within-species strain restructuring in species that persisted, and indicated that a previously reported E. coli abundance increase reflected lower microbial load rather than a strain-level change.

    12 healthy men sampled at days 0, 4, 8, 42 and 180: reduced prevalence of 1,140 species at days 4 and 8; cANI shifts in 282 species, 96 still present at day 180 (A. putredinis cANI s.d. 0.40 baseline vs 0.83 at day 180). E. coli: StrainSpy adjusted p = 0.532; abundance test adjusted p = 0.028, then 0.44 with predicted load as covariate. · StrainSpy reveals strain dynamics following…; Testing for associations in gut metagenomic…; Estimating the microbial load

  3. ShownA single mixed model across 17 colorectal cancer (CRC) cohorts recovered known taxa such as Fusobacterium animalis and found new strain-level shifts in E. coli and Streptococcus vestibularis; among E. coli-positive samples, several E. coli adhesin genes were more prevalent in CRC.

    3,414 metagenomes (1,322 CRC, 614 adenoma, 1,478 controls; 10 countries): 64 species with differential prevalence and 9 strain-level associations; adenoma did not differ from controls; E. coli hit to GCF_000458275.1 (adjusted p = 0.04, 0.07% ANI difference, ~3−4 kbp divergence); in E. coli-positive paired-end samples, fimbrial, curli and autotransporter adhesins were more prevalent in CRC. · StrainSpy reveals novel strain-level…; Testing for associations in gut metagenomic…; Searching for E. coli adhesin genes in CRC…

  4. ShownIn read spike-in tests, StrainSpy recovered an implanted strain at low coverage where a StrainPhlAn-AnPan phylogenetic approach could not, and empirical Bayes priors removed false positives in a small-sample design.

    Escherichia sp. GCF_009820115.1 spiked into 20 cases at 1×: adjusted p = 1.2 × 10^−144, plus 17 related strains (>97.5% ANI). With 16 cases and 17 controls, the weak MAP prior gave false positives while empirical Bayes recovered the spike (adjusted p = 1.51 × 10^−7). P. distasonis spike found at 1×, 3×, 5×; AnPan missed it at 1×. · StrainSpy provides robust results in simulated…; Spiking reads from Escherichia sp…; Figure legends

  5. SuggestedElastic-net classifiers using StrainSpy-selected cANI features predicted CRC status in held-out cohorts at levels similar to published abundance-based models; immunotherapy response prediction across cohorts was inconsistent and no pooled response associations were found.

    Leave-one-dataset-out AUC 0.70 to 0.93 vs 0.71 to 0.95 reported for abundance-based models by Piccinno et al. (different classifier and cohorts). Immunotherapy: GYN-trained model AUC 0.77 in UGB and 0.45 in NEN; UGB-trained AUC 0.6 in GYN and 0.58 in NEN; melanoma AUC 0.21–0.82. · Containment ANI features identified by StrainSpy…; Pooled colorectal cancer dataset; Response to immunotherapy in pancancer

Methods, data and reproducibility

Fits a generalized linear mixed model per strain with cANI (from Sylph against a 99%-dereplicated GTDB catalogue) as the outcome, using a zero-inflated beta likelihood for depletion plus replacement, MAP estimation in glmmTMB with weak or empirical Bayes priors, and Benjamini-Hochberg or harmonic mean p-value correction.

  • Healthy gut metagenomes from Zeevi et al., split into simulated cases and controls for signal-implantation and read spike-in benchmarks · n = 200 metagenomes
  • Longitudinal faecal metagenomes from the Palleja et al. study of broad-spectrum antibiotic exposure (vancomycin, gentamicin, meropenem), five time points · human (healthy adult males) · n = 12 subjects (9 sampled at day 4)
  • Pooled public fecal metagenomes from a CRC meta-analysis: 17 datasets in 10 countries covering CRC, adenoma and controls · human · n = 3,414 (1,322 CRC, 614 adenoma, 1,478 controls)
  • Faecal shotgun metagenomes from immune checkpoint inhibitor cohorts: three rare cancer datasets (Gunjur et al.) plus six melanoma datasets · human · n = 640 metagenomes
  • GTDB-derived reference genome catalogue dereplicated at 99% identity, built as a Sylph database with GTDB taxonomy · n = 18,027 strains retained after filtering in simulations · Zenodo DOI 10.5281/zenodo.21796878
  • Simulated paired-end reads from Escherichia sp. and Parabacteroides distasonis genomes spiked into baseline metagenomes · Escherichia sp.; Parabacteroides distasonis · n = 20 spiked metagenomes per experiment · GCF_009820115.1; GCF_024791025.1
Methods & tools
StrainSpy (R package), Sylph v.0.9.0, glmmTMB, GTDB, art_modern v.1.3.4 (art 2.5.8), StrainPhlAn, AnPan, ARIBA v.2.14.7, adhesiomeR, caret v.7.0.1 with glmnet v.5.0 elastic net, Microbial Load Predictor (Galaxy model), skani
Software & environment
R package under MIT licence; Sylph v.0.9.0, art_modern v.1.3.4 (art 2.5.8), ARIBA v.2.14.7, caret v.7.0.1, glmnet v.5.0; MAP model fitting via glmmTMB.
Reproducibility gaps
Sample accessions and cohort metadata are provided only in Supplementary Tables 8-12; clinical cohorts are reused public metagenomes from earlier studies.
Limitations (6)
  • Cannot capture functional microbiome changes that occur without any shift in a strain's or species' ANI, so the authors present it as complementary to abundance-based methods rather than a replacement. · stated by the authors
  • Identity-driven batch effects, such as geographically distinct circulating lineages, are not removed by cANI modelling and must be specified as covariates or random effects in the mixed model. · stated by the authors
  • Does not reconstruct strain haplotypes or model species population structure; phylogeny-based tools such as AnPan are needed for those questions. · stated by the authors
  • With broad catalogues such as GTDB, a significant reference genome must be interpreted against the available reference diversity, since the precise strain carried by a sample is rarely represented in the database. · stated by the authors
  • CRC prediction AUCs are compared only indirectly with prior abundance-based results, since training cohorts and classifiers (elastic net vs Random Forest) differ. · stated by the authors
  • Microbial load in the antibiotic cohort was predicted computationally rather than measured, so the conclusion that the E. coli abundance signal was load-driven rests on a model-derived covariate in 12 subjects. · noted in this summary

Read this if you run strain-resolved microbiome association studies with longitudinal or multi-cohort designs and want a threshold-free alternative to relative abundance testing that is robust to microbial load.

Citations
0

Reproducibility ledger

Nobody has reported trying to reproduce this paper yet. Did you replicate an experiment or run the code? Your report saves the next lab a week.

Sign in to file a reproducibility report.

Related papers

Similarity Map
SelectedSame fieldCross-field
THISGhazi '25Galbraith '26Lin '19Luo '25Zhao '20Paz '26
66% semantic matchBioinformatics
Andrew R. Ghazi, Kelsey N. Thompson, Amrisha Bhosle+6 more·bioRxiv·Jan 6, 2025· 6 citations