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.
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.
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: 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 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 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. The two main plots are rendered below. Set
plot = FALSE in either call to get the underlying table
instead of the figure.
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:
- 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.
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.3 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.61 cachem_1.1.0
## [23] aplot_0.3.2 jsonlite_2.0.0
## [25] FMStable_0.1-4 DelayedArray_0.39.8
## [27] BiocParallel_1.47.0 tweenr_2.0.3
## [29] RTMB_2.0 parallel_4.6.1
## [31] R6_2.6.1 stringi_1.8.9
## [33] bslib_0.12.0 RColorBrewer_1.1-3
## [35] boot_1.3-32 estimability_2.0.0
## [37] numDeriv_2016.8-1.1 GenomicRanges_1.65.4
## [39] jquerylib_0.1.4 Rcpp_1.1.2
## [41] Seqinfo_1.3.2 SummarizedExperiment_1.43.0
## [43] knitr_1.52 snow_0.4-4
## [45] zoo_1.9-1 R.utils_2.13.0
## [47] IRanges_2.47.5 Matrix_1.7-6
## [49] splines_4.6.1 glmmTMB_1.1.15.2
## [51] tidyselect_1.2.1 viridis_0.6.5
## [53] abind_1.4-8 yaml_2.3.12
## [55] TMB_1.9.25 codetools_0.2-20
## [57] lmerTest_3.2-1 lattice_0.23-1
## [59] tibble_3.3.1 Biobase_2.73.2
## [61] treeio_1.37.1 withr_3.0.3
## [63] S7_0.2.2 evaluate_1.0.5
## [65] gridGraphics_0.5-1 polyclip_1.10-7
## [67] pillar_1.11.1 BiocManager_1.30.27
## [69] ggtree_4.3.0 MatrixGenerics_1.25.0
## [71] stats4_4.6.1 reformulas_0.4.4
## [73] ggfun_0.2.2 generics_0.1.4
## [75] S4Vectors_0.51.10 ggplot2_4.0.3
## [77] scales_1.4.0 tidytree_0.4.8
## [79] minqa_1.2.8 harmonicmeanp_3.0.1
## [81] glue_1.8.1 gdtools_0.5.1
## [83] emmeans_2.0.4 lazyeval_0.2.3
## [85] maketools_1.3.2 tools_4.6.1
## [87] sys_3.4.3 data.table_1.18.6.1
## [89] lme4_2.0-6 ggiraph_0.9.6
## [91] mvtnorm_1.4-2 buildtools_1.0.0
## [93] fs_2.1.0 grid_4.6.1
## [95] tidyr_1.3.2 ape_5.8-1
## [97] rbibutils_2.4.1 nlme_3.1-171
## [99] patchwork_1.3.2 ggforce_0.5.0
## [101] cli_3.6.6 rappdirs_0.3.4
## [103] fontBitstreamVera_0.1.1 viridisLite_0.4.3
## [105] ggthemes_7.0.0 S4Arrays_1.13.2
## [107] dplyr_1.2.1 gtable_0.3.6
## [109] R.methodsS3_1.8.2 yulab.utils_0.2.5
## [111] sass_0.4.10 digest_0.6.39
## [113] fontquiver_0.2.1 BiocGenerics_0.59.12
## [115] SparseArray_1.13.4 ggrepel_0.9.8
## [117] ggplotify_0.1.4 htmlwidgets_1.6.4
## [119] farver_2.1.2 htmltools_0.5.9
## [121] R.oo_1.27.1 lifecycle_1.0.5
## [123] fontLiberation_0.1.0 MASS_7.3-66