strainspy Workflow
28 October 2026
Contents
- 1 Why strainspy?
- 2 Installation
- 3 Load package and example data
- 4 Choosing key wrapper inputs
- 5 Fit the main model (
glmZiBFit) - 6 Find top hits
- 7 Add taxonomy context
- 8 Generate plot-ready outputs
- 9 Diversity and prediction preparation
- 10 Working with other profile inputs
- 11 Notes
- 12 Next steps
- 13 Session information
1 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.
2 Installation
Install the release version from Bioconductor:
if (!requireNamespace("BiocManager", quietly = TRUE)) {
install.packages("BiocManager")
}
BiocManager::install("strainspy")Conda, container and GitHub installation options are described in the package README.
3 Load package and example data
The bundled files are described in inst/extdata/README.md
(system.file("extdata", "README.md", package = "strainspy")).
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: 200cat("Retained strains:", nrow(se), "\n")## Retained strains: 634cat("min_nonzero used:", min_nonzero, "\n")## min_nonzero used: 204 Choosing 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)).
5 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%fit## 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
## ------------------------------6 Find top hits
top_hits() is usually the first table you inspect after fitting.
It highlights contigs with the 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 17 Add 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.03918 Generate plot-ready outputs
The README demonstrates visualization early because it helps with interpretation.
The two main plots are rendered below. Set plot = FALSE in either call to get the underlying table instead of the figure.
plot_volcano(fit, coef = 2)
plot_manhattan(fit, coef = 2, taxonomy = taxonomy)
9 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] 200class(beta_dist)## [1] "dist"dim(pred_df)## [1] 200 510 Working 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"class(mp)## [1] "SummarizedExperiment"
## attr(,"package")
## [1] "SummarizedExperiment"11 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.
12 Next steps
After this walkthrough, common next analyses are:
- Tighten
alphaintop_hits()to focus on the strongest hits. - Use taxonomy-informed Manhattan summaries to prioritize clades.
- Export
prep_for_prediction()outputs to train external classifiers.
13 Session information
sessionInfo()## R version 4.6.1 Patched (2026-06-24 r90190)
## Platform: x86_64-apple-darwin20
## Running under: macOS Ventura 13.7.8
##
## Matrix products: default
## BLAS: /Library/Frameworks/R.framework/Versions/4.6-x86_64/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/4.6-x86_64/Resources/lib/libRlapack.dylib; LAPACK version 3.12.1
##
## locale:
## [1] C/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
##
## time zone: America/New_York
## tzcode source: internal
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] strainspy_0.99.3 BiocStyle_2.41.0
##
## loaded via a namespace (and not attached):
## [1] RColorBrewer_1.1-3 jsonlite_2.0.0
## [3] magrittr_2.0.5 magick_2.9.1
## [5] TH.data_1.1-5 estimability_2.0.0
## [7] farver_2.1.2 nloptr_2.2.1
## [9] rmarkdown_2.32 fs_2.1.0
## [11] vctrs_0.7.3 minqa_1.2.8
## [13] ggtree_4.3.0 tinytex_0.61
## [15] htmltools_0.5.9 S4Arrays_1.13.2
## [17] SparseArray_1.13.4 gridGraphics_0.5-1
## [19] sass_0.4.10 bslib_0.12.0
## [21] htmlwidgets_1.6.4 sandwich_3.1-3
## [23] emmeans_2.0.4 zoo_1.9-1
## [25] cachem_1.1.0 TMB_1.9.25
## [27] lifecycle_1.0.5 pkgconfig_2.0.3
## [29] Matrix_1.7-6 R6_2.6.1
## [31] fastmap_1.2.0 rbibutils_2.4.1
## [33] MatrixGenerics_1.25.0 digest_0.6.39
## [35] numDeriv_2016.8-1.1 aplot_0.3.2
## [37] patchwork_1.3.2 S4Vectors_0.51.10
## [39] GenomicRanges_1.65.4 labeling_0.4.3
## [41] polyclip_1.10-7 abind_1.4-8
## [43] mgcv_1.9-4 compiler_4.6.1
## [45] fontquiver_0.2.1 withr_3.0.3
## [47] S7_0.2.2 BiocParallel_1.47.0
## [49] viridis_0.6.5 ggforce_0.5.0
## [51] R.utils_2.13.0 MASS_7.3-66
## [53] rappdirs_0.3.4 DelayedArray_0.39.8
## [55] harmonicmeanp_3.0.1 tools_4.6.1
## [57] otel_0.2.0 ape_5.8-1
## [59] R.oo_1.27.1 glue_1.8.1
## [61] nlme_3.1-171 grid_4.6.1
## [63] generics_0.1.4 snow_0.4-4
## [65] glmmTMB_1.1.15.2 gtable_0.3.6
## [67] R.methodsS3_1.8.2 tidyr_1.3.2
## [69] data.table_1.18.6.1 utf8_1.2.6
## [71] XVector_0.53.0 BiocGenerics_0.59.12
## [73] ggrepel_0.9.8 pillar_1.11.1
## [75] stringr_1.6.0 yulab.utils_0.2.5
## [77] splines_4.6.1 dplyr_1.2.1
## [79] tweenr_2.0.3 treeio_1.37.1
## [81] lattice_0.23-1 survival_3.8-12
## [83] tidyselect_1.2.1 fontLiberation_0.1.0
## [85] knitr_1.52 RTMB_2.0
## [87] fontBitstreamVera_0.1.1 reformulas_0.4.4
## [89] gridExtra_2.3.1 bookdown_0.48
## [91] IRanges_2.47.5 Seqinfo_1.3.2
## [93] SummarizedExperiment_1.43.0 stats4_4.6.1
## [95] xfun_0.62 Biobase_2.73.2
## [97] matrixStats_1.5.0 stringi_1.8.9
## [99] lazyeval_0.2.3 ggfun_0.2.2
## [101] yaml_2.3.12 boot_1.3-32
## [103] evaluate_1.0.5 codetools_0.2-20
## [105] gdtools_0.5.2 tibble_3.3.1
## [107] BiocManager_1.30.27 ggplotify_0.1.4
## [109] cli_3.6.6 FMStable_0.1-4
## [111] xtable_1.8-8 systemfonts_1.3.2
## [113] Rdpack_2.6.6 jquerylib_0.1.4
## [115] dichromat_2.0-1 Rcpp_1.1.2
## [117] coda_0.19-4.1 parallel_4.6.1
## [119] ggplot2_4.0.3 lme4_2.1-0
## [121] ggthemes_7.0.0 viridisLite_0.4.3
## [123] mvtnorm_1.4-2 tidytree_0.4.8
## [125] ggiraph_0.9.6 lmerTest_3.2-1
## [127] scales_1.4.0 purrr_1.2.2
## [129] rlang_1.3.0 multcomp_1.4-32