strainspy Workflow
Why strainspy?
strainspy fits strain-level association models from
ANI/abundance profiles and metadata. This guide follows the README
workflow and is written as a practical first analysis.
You will walk through the core path most users care about:
- Load metadata and ANI profiles.
- Filter to informative strains.
- Fit a zero-inflated beta model with
glmZiBFit(). - Rank and interpret hits with
top_hits(). - Explore taxonomy-aware summaries and plot-ready outputs.
- Prepare data for downstream prediction and diversity analyses.
All examples use bundled data so the workflow runs end-to-end.
Load package and example data
library(strainspy)
meta_path <- system.file("extdata", "example_metadata.csv.gz", package = "strainspy")
sylph_path <- system.file("extdata", "example_sylph_profile.tsv.gz", package = "strainspy")
taxonomy_path <- system.file("extdata", "example_taxonomy.tsv.gz", package = "strainspy")
meta <- read.csv(meta_path)
meta$Case_status <- factor(meta$Case_status)
taxonomy <- read_taxonomy(taxonomy_path)
se <- read_sylph(sylph_path, meta_data = meta)
# Recommended default: keep strains observed in at least 10% of samples
min_nonzero <- ceiling(0.10 * ncol(se))
se <- filter_by_presence(se, min_nonzero = min_nonzero)
cat("Samples:", ncol(se), "\n")## Samples: 200## Retained strains: 634## min_nonzero used: 20Choosing key wrapper inputs
If you use the high-level strainspy() wrapper, three
inputs should be chosen deliberately:
design_formula: the predictors you want to test.coef: which coefficient index from that design to report intop_hits().min_nonzero: prevalence threshold (recommended default:ceiling(0.10 * n_samples)).
For the design ~ Case_status + Age_at_collection,
coef = 2 targets Case_status. You can verify
term order with
colnames(model.matrix(~ Case_status + Age_at_collection, meta)).
Fit the main model (glmZiBFit)
design <- as.formula("~ Case_status + Age_at_collection")
fit <- glmZiBFit(
se,
design,
nthreads = 2L,
return_vcov = FALSE
)## | | | 0% | |========== | 14% | |==================== | 29% | |============================== | 43% | |======================================== | 57% | |================================================== | 71% | |============================================================ | 86% | |======================================================================| 100%## strainspy_fit object
## ------------------------------
## Design formula: ~Case_status + Age_at_collection
## ------------------------------
## Converged=634, Failed=0
## ------------------------------
## Coefficients: 634 features x 3 terms
## ------------------------------
## ZI Coefficients: 634 features x 3 terms
## ------------------------------
## Residuals available: Yes
## Call:
## glmZiBFit(se = se, design = design, nthreads = 2L, return_vcov = FALSE)
## ------------------------------
## Priors used:
## strainspy_priors object
## Method: preset_weak
##
## Priors:
## - fixef : 2 coefficients
## - fixef_zi : 2 coefficients
##
## Example priors:
## prior class coef
## 1 normal(0,5) fixef Case_statusPD
## 2 normal(0,5) fixef Age_at_collection
## 3 normal(0,5) fixef_zi Case_statusPD
## 4 normal(0,5) fixef_zi Age_at_collection
##
## Bootstrap summary for fixef priors:
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## NA NA NA NaN NA NA
##
## Bootstrap summary for fixef_zi priors:
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## NA NA NA NaN NA NA
## ------------------------------Find top hits
top_hits() is usually the first table you inspect after
fitting. It highlights contigs with strongest phenotype association and
provides adjusted p-values.
hits <- top_hits(fit, coef = 2, alpha = 1)
head(hits[, c("Contig_name", "coefficient", "p_value", "p_adjust")])## # A tibble: 6 × 4
## Contig_name coefficient p_value p_adjust
## <chr> <dbl> <dbl> <dbl>
## 1 NZ_WSNW01000023.1 Coprobacillus cateniformis str… 0.167 1.80e-1 1
## 2 UREB01000001.1 uncultured Clostridiales bacteriu… 0.248 1.01e-4 0.0638
## 3 DVOB01000010.1 MAG TPA_asm: Candidatus Copromorp… 0.231 1.08e-4 0.0681
## 4 JAHZEH010000001.1 MAG: Clostridia bacterium isol… 0.350 1.40e-4 0.0882
## 5 NZ_JXXK01000001.1 Ruthenibacterium lactatiforman… 0.121 1.41e-4 0.0887
## 6 NZ_VICC01000001.1 Actinomyces oris strain CCUG 3… 0.111 1.43e-1 1Add taxonomy context
When taxonomy is available, hadjust() summarizes
significance across taxonomic levels.
hier_hits <- hadjust(fit, coef = 2, taxonomy = taxonomy)
head(hier_hits[, c("Level", "Model", "Name", "p_adjust")])## # A tibble: 6 × 4
## Level Model Name p_adjust
## <chr> <chr> <chr> <dbl>
## 1 Phylum Beta Bacillota_A 0.0151
## 2 Class Beta Clostridia 0.0151
## 3 Order Beta Christensenellales 0.0320
## 4 Phylum Zero-Inflated Bacillota_I 0.0391
## 5 Class Zero-Inflated Bacilli_A 0.0391
## 6 Order Zero-Inflated Erysipelotrichales 0.0391Generate plot-ready outputs
The README demonstrates visualization early because it helps with
interpretation. To keep this vignette stable in automated builds, we
return data frames (plot = FALSE) rather than rendering all
figures. You can set plot = TRUE interactively.
Diversity and prediction preparation
Beyond hit calling, strainspy supports ecological
summaries and model-ready tables.
rich <- estimate_sample_richness(se, taxonomy)
beta_dist <- estimate_beta_diversity(se, distance_only = TRUE)
pred_df <- prep_for_prediction(
sy = se,
outcome = "Case_status",
contigs = rownames(se)[1:3],
covariates = "Age_at_collection",
use_genome_names = FALSE
)
length(rich$species_richness)## [1] 200## [1] "dist"## [1] 200 5Working with other profile inputs
The same downstream workflow can be used with other supported profile readers.
example_sourmash_path <- system.file("extdata", "example_sourmash.csv.gz", package = "strainspy")
example_metaphlan_path <- system.file("extdata", "metaphlan_merged.tsv.gz", package = "strainspy")
sm <- read_sourmash(example_sourmash_path, meta)
mp <- read_metaphlan(example_metaphlan_path, meta)
class(sm)## [1] "SummarizedExperiment"
## attr(,"package")
## [1] "SummarizedExperiment"## [1] "SummarizedExperiment"
## attr(,"package")
## [1] "SummarizedExperiment"Notes
- If taxonomy is not available, sequence-level workflows
(
top_hits,plot_volcano) still run. - For optional features requiring additional packages, see individual function help pages.
Next steps
After this walkthrough, common next analyses are:
- Re-fit on the full filtered dataset (
se) and tightenalphaintop_hits(). - Use taxonomy-informed Manhattan summaries to prioritize clades.
- Export
prep_for_prediction()outputs to train external classifiers.
Session information
## R version 4.6.1 (2026-06-24)
## Platform: x86_64-pc-linux-gnu
## Running under: Ubuntu 26.04.1 LTS
##
## Matrix products: default
## BLAS: /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3
## LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.32.so; LAPACK version 3.12.0
##
## locale:
## [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
## [3] LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8
## [5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8
## [7] LC_PAPER=en_US.UTF-8 LC_NAME=C
## [9] LC_ADDRESS=C LC_TELEPHONE=C
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
##
## time zone: Etc/UTC
## tzcode source: system (glibc)
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] strainspy_0.99.0 BiocStyle_2.41.0
##
## loaded via a namespace (and not attached):
## [1] Rdpack_2.6.6 gridExtra_2.3.1
## [3] sandwich_3.1-3 rlang_1.3.0
## [5] magrittr_2.0.5 otel_0.2.0
## [7] matrixStats_1.5.0 compiler_4.6.1
## [9] mgcv_1.9-4 systemfonts_1.3.2
## [11] vctrs_0.7.3 stringr_1.6.0
## [13] pkgconfig_2.0.3 fastmap_1.2.0
## [15] XVector_0.53.0 labeling_0.4.3
## [17] utf8_1.2.6 rmarkdown_2.32
## [19] nloptr_2.2.1 purrr_1.2.2
## [21] xfun_0.60 cachem_1.1.0
## [23] aplot_0.3.1 jsonlite_2.0.0
## [25] FMStable_0.1-4 DelayedArray_0.39.6
## [27] BiocParallel_1.47.0 tweenr_2.0.3
## [29] parallel_4.6.1 R6_2.6.1
## [31] bslib_0.12.0 stringi_1.8.9
## [33] RColorBrewer_1.1-3 boot_1.3-32
## [35] estimability_2.0.0 GenomicRanges_1.65.4
## [37] jquerylib_0.1.4 numDeriv_2016.8-1.1
## [39] Rcpp_1.1.2 Seqinfo_1.3.2
## [41] SummarizedExperiment_1.43.0 knitr_1.52
## [43] snow_0.4-4 zoo_1.9-0
## [45] R.utils_2.13.0 IRanges_2.47.5
## [47] Matrix_1.7-6 splines_4.6.1
## [49] glmmTMB_1.1.14 tidyselect_1.2.1
## [51] viridis_0.6.5 abind_1.4-8
## [53] yaml_2.3.12 TMB_1.9.25
## [55] codetools_0.2-20 lmerTest_3.2-1
## [57] lattice_0.23-1 tibble_3.3.1
## [59] Biobase_2.73.2 treeio_1.37.0
## [61] withr_3.0.3 S7_0.2.2
## [63] evaluate_1.0.5 gridGraphics_0.5-1
## [65] polyclip_1.10-7 pillar_1.11.1
## [67] BiocManager_1.30.27 ggtree_4.3.0
## [69] MatrixGenerics_1.25.0 stats4_4.6.1
## [71] reformulas_0.4.4 ggfun_0.2.1
## [73] generics_0.1.4 S4Vectors_0.51.10
## [75] ggplot2_4.0.3 scales_1.4.0
## [77] tidytree_0.4.8 minqa_1.2.8
## [79] harmonicmeanp_3.0.1 glue_1.8.1
## [81] gdtools_0.5.1 emmeans_2.0.4
## [83] lazyeval_0.2.3 maketools_1.3.2
## [85] tools_4.6.1 sys_3.4.3
## [87] data.table_1.18.6.1 lme4_2.0-6
## [89] ggiraph_0.9.6 mvtnorm_1.4-2
## [91] buildtools_1.0.0 fs_2.1.0
## [93] grid_4.6.1 tidyr_1.3.2
## [95] ape_5.8-1 rbibutils_2.4.1
## [97] nlme_3.1-171 patchwork_1.3.2
## [99] ggforce_0.5.0 cli_3.6.6
## [101] rappdirs_0.3.4 fontBitstreamVera_0.1.1
## [103] viridisLite_0.4.3 ggthemes_6.0.0
## [105] S4Arrays_1.13.0 dplyr_1.2.1
## [107] gtable_0.3.6 R.methodsS3_1.8.2
## [109] yulab.utils_0.2.5 sass_0.4.10
## [111] digest_0.6.39 fontquiver_0.2.1
## [113] BiocGenerics_0.59.12 SparseArray_1.13.2
## [115] ggplotify_0.1.3 htmlwidgets_1.6.4
## [117] farver_2.1.2 R.oo_1.27.1
## [119] htmltools_0.5.9 lifecycle_1.0.5
## [121] fontLiberation_0.1.0 MASS_7.3-66