Package 'strainspy'

Description

This function fits a linear regression model using stats::lm() or lmerTest::lmer(). The latter is used for models with random effects. It takes a SummarizedExperiment object as input, along with a user-defined formula, and fits a linear model on the assay data.

Usage

abundanceFit(
  se,
  design,
  nthreads = 1,
  scale_continuous = TRUE,
  transform = NULL,
  BPPARAM = NULL,
  progress = TRUE
)

Arguments

Details

Generally, these linear models are useful to model compositional abundance data (e.g. Taxonomic_abundance in Sylph Profile) with appropriate transforms.

Value

A strainspy_fit object with the following components:

Examples

library(strainspy)

example_meta_path <- system.file("extdata", "example_metadata.csv.gz", 
package = "strainspy")
example_meta <- read.csv(example_meta_path)
example_path <- system.file("extdata", "example_sylph_profile.tsv.gz", 
package = "strainspy")

se <- read_sylph(example_path, example_meta, 'Taxonomic_abundance')
se <- filter_by_presence(se, 20)

design <- as.formula(" ~ Case_status + Age_at_collection")
fit <- abundanceFit(se[1:10,],  design, nthreads=2, transform='CLR')
top_hits(fit, alpha = 0.5)

Description

This function calculates the distance between samples in a SummarizedExperiment object representing strain level ANI data.

Usage

ani_distance(
  se,
  taxonomy = NULL,
  tax_level = "Species",
  collapse_method = "representative"
)

Arguments

Value

A dist object of distances between samples.

Examples

example_path <- system.file("extdata", "example_sylph_profile.tsv.gz", 
package = "strainspy")
tax_path <- system.file("extdata", "example_taxonomy.tsv.gz", 
package = "strainspy")
se <- read_sylph(example_path)
tax <- read_taxonomy(tax_path)
distance_matrix <- ani_distance(se, tax)

Description

This function fits a logistic regression model similar to that published in the original Sylph paper. It takes a SummarizedExperiment object as input, along with a user-defined formula.

Usage

caseControlFit(
  se,
  design,
  min_identity = 0.98,
  nthreads = 1,
  scale_continuous = TRUE,
  BPPARAM = NULL,
  progress = TRUE
)

Arguments

Value

A strainspy_fit object with the following components:

Examples

library(strainspy)

example_meta_path <- system.file("extdata", "example_metadata.csv.gz", 
package = "strainspy")
example_meta <- read.csv(example_meta_path)
example_meta$Case_status <- factor(example_meta$Case_status)
example_path <- system.file("extdata", "example_sylph_profile.tsv.gz", 
package = "strainspy")
se <- read_sylph(example_path, example_meta)
se <- filter_by_presence(se)

design <- as.formula("Case_status ~ Value + Age_at_collection")

fit <- caseControlFit(se[1:10,], design, nthreads=2)

Description

In the zero-inflated beta model, beta effects represent changes in strain identity and are modeled on the logit scale. While this scale is effective for distinguishing phenotype groups statistically, the resulting effect sizes are not directly interpretable as differences in mean ANI. The same applies to effect sizes from the ordinal beta model.

Usage

comp_ani_diff_and_posthoc_test(
  se,
  fit,
  th = NULL,
  simplified_formula = NULL,
  beta_min_nz = 0.1,
  beta_min_ani_diff = 0.015,
  nthreads = 1L,
  reorder_hits = FALSE,
  BPPARAM = NULL,
  progress = TRUE
)

Arguments

Details

This function converts logit-scale beta effects into response-scale mean ANI differences using marginal means and pairwise contrasts (requires emmeans). It therefore provides an interpretable estimate of the magnitude of strain-level differences between phenotype groups, along with post hoc quality checks for poorly supported beta signals.

Posthoc testing uses two approaches. First, it counts the number of non-zero observations supporting each beta effect. Second, it compares the estimated mean ANI difference against a preset threshold. A comment is added to each hit depending on the outcome.

By default, this function uses the formula from getDesign(fit). However, if the formula contains multiple covariates, this can cause emmeans to behave unexpectedly. For example, if getDesign(fit) is ~ disease + age + sex + country + (1|study), this can cause this function to throw errors or warnings. Using a simplified_formula = ~ disease to override can be a useful approach.

This function is designed to directly work on the tibble returned from top_hits().

Value

A data.frame similar to th, with added columns giving details on: Component (Beta, ZI or both), Mean ANI_Difference, Contrast, Non-Zero Counts of contrast components, and a general Comment on each hit. The comment will be NA if the hit passes all posthoc tests.

Examples

if (requireNamespace("emmeans", quietly = TRUE)) {
  meta <- read.csv(system.file("extdata", "example_metadata.csv.gz", 
  package = "strainspy"))
  meta$Case_status <- factor(meta$Case_status)
  se <- read_sylph(system.file("extdata", "example_sylph_profile.tsv.gz", 
  package = "strainspy"), meta_data = meta)
  se <- filter_by_presence(se, min_nonzero = 30)[1:10, ]
  fit <- glmZiBFit(se, ~ Case_status + Age_at_collection, nthreads = 1)
  th <- top_hits(fit, coef = 2, alpha = 1)
  posthoc <- comp_ani_diff_and_posthoc_test(se, fit, th, nthreads = 1)
  head(posthoc$ANI_Difference)
}

Description

Fits fast binomial and beta models to estimate MAP prior SDs for use in zero-inflated and ordinal beta models fitted using glmmTMB.

Usage

compute_eb_priors(
  se,
  design,
  nthreads = 1L,
  BPPARAM = NULL,
  low_cutoff = 0.1,
  high_cutoff = 10,
  est_dispersion_prior = FALSE,
  progress = TRUE
)

Arguments

Value

An object of class strainspy_priors.

Examples

if (requireNamespace("fastglm", quietly = TRUE) && 
requireNamespace("betareg", quietly = TRUE)) {
  meta <- read.csv(system.file("extdata", "example_metadata.csv.gz", 
  package = "strainspy"))
  meta$Case_status <- factor(meta$Case_status)
  se <- read_sylph(system.file("extdata", "example_sylph_profile.tsv.gz", 
  package = "strainspy"), meta_data = meta)
  se <- filter_by_presence(se, min_nonzero = 30)
  pri <- compute_eb_priors(se[1:10, ], 
  ~ Case_status + Age_at_collection, nthreads = 1)
  class(pri)
}

Description

Constructor for strainspy_priors using preset or user defined priors

Usage

define_priors(
  se,
  design,
  method = "preset_weak",
  add_dispersion_prior = FALSE,
  priors_df = NULL
)

Arguments

Value

An object of class strainspy_priors.

Examples

meta <- read.csv(system.file("extdata", "example_metadata.csv.gz",
 package = "strainspy"))
meta$Case_status <- factor(meta$Case_status)
se <- read_sylph(system.file("extdata", "example_sylph_profile.tsv.gz", 
package = "strainspy"), meta_data = meta)
pri <- suppressWarnings(define_priors(se[1:10, ], 
~ Case_status + Age_at_collection, method = "preset_weak"))
class(pri)

Description

Computes between-sample distances from the ANI matrix and, unless distance_only = TRUE, runs NMDS, PERMANOVA and a beta-dispersion test grouped by phenotype. The vegan package is required for method = "jaccard" and for everything beyond the distance matrix.

Usage

estimate_beta_diversity(
  se,
  method = "manhattan",
  distance_only = FALSE,
  ani_threshold = 95,
  phenotype = NULL,
  show_plots = FALSE,
  return_plots = FALSE
)

Arguments

Value

If distance_only = TRUE, returns distance object. Otherwise a list with:

• distance — distance matrix

• nmds — NMDS object

• ordination — NMDS coordinates with group

• permanova — PERMANOVA results

• betadisper — betadisper object

• plots — list of ggplot2 objects (if return_plots = TRUE)

Examples

meta <- read.csv(system.file("extdata", "example_metadata.csv.gz", 
package = "strainspy"))
meta$Case_status <- factor(meta$Case_status)
se <- read_sylph(system.file("extdata", "example_sylph_profile.tsv.gz", 
package = "strainspy"), meta_data = meta)
se <- filter_by_presence(se, min_nonzero = 30)
d <- estimate_beta_diversity(se[1:15, ], distance_only = TRUE)

Description

This function estimates sample-level species richness and strain coverage using ANI-based presence/absence. A strain is considered present if ANI >= ani_threshold.

Usage

estimate_sample_richness(
  se,
  taxonomy,
  ani_threshold = 95,
  report_strain_richness = FALSE
)

Arguments

Details

Species richness is defined as the number of species with at least one strain exceeding the ANI threshold. This requires taxonomy data, which can be read in using read_taxonomy().

Experimental Strain coverage is a database-dependent, normalised score that increases with the number of reference strains detected within a species. When only one strain is present for a species, it contributes 1 towards strain coverage. As more strains are detected, the contribution increases toward 2, reaching 2 when all reference strains are detected. Although this is designed to prevent inflation for species with many reference strains, this measures the coverage of reference strain space in the database, not true biological strain richness. Therefore, this value is not returned by default.

The strain coverage metric is computed as:

$1 + (k - 1) / N$

where $k$ is the number of reference strains detected for a species and $N$ is the total number of reference strains for that species in the database.

This approach does not resolve true strain multiplicity and may overestimate coverage when a single biological strain matches many highly similar reference genomes. Alternative approaches are being investigated.

Value

A named list with:

species_richness

Integer vector giving the number of species detected per sample.

strain_coverage

NA or Numeric vector giving the database-normalised strain coverage per sample.

Examples

meta <- read.csv(system.file("extdata", "example_metadata.csv.gz", 
package = "strainspy"))
se <- read_sylph(system.file("extdata", "example_sylph_profile.tsv.gz", 
package = "strainspy"), meta_data = meta)
tax <- read_taxonomy(system.file("extdata", "example_taxonomy.tsv.gz", 
package = "strainspy"))
rich <- estimate_sample_richness(se[1:20, ], tax)

Description

Access fixef priors needed for ordinal beta regression

Usage

extract_fixef_priors(object)

## S4 method for signature 'strainspy_priors'
extract_fixef_priors(object)

Arguments

Value

A data frame of prior parameters for fixed effects.

Examples

meta <- read.csv(system.file("extdata", "example_metadata.csv.gz", 
package = "strainspy"))
meta$Case_status <- factor(meta$Case_status)
se <- read_sylph(system.file("extdata", "example_sylph_profile.tsv.gz", 
package = "strainspy"), meta_data = meta)
pri <- suppressWarnings(define_priors(se[1:10, ], 
~ Case_status + Age_at_collection, method = "preset_weak"))
head(extract_fixef_priors(pri))

Description

Extract data from a StrainSpy SummarizedExperiment object in long format

Usage

extract_strainspy(
  se,
  variables = NULL,
  contigs = NULL,
  taxonomy = NULL,
  drop_zeros = FALSE
)

Arguments

Value

A long-format data.frame with sample metadata, contig identity columns, ANI values, row-level annotations, and optional taxonomy columns.

Examples

library(strainspy)

example_meta_path <- system.file("extdata", "example_metadata.csv.gz", 
package = "strainspy")
example_meta <- read.csv(example_meta_path)
# Change variable from character to factor
example_meta$Case_status = factor(example_meta$Case_status)

example_path <- system.file("extdata", "example_sylph_profile.tsv.gz", 
package = "strainspy")
se <- read_sylph(example_path, example_meta)

tax_path <- system.file("extdata", "example_taxonomy.tsv.gz", 
package = "strainspy")
tax <- read_taxonomy(tax_path)

extract_strainspy(se, variables = c("Case_status", "Age_at_collection"), 
contigs = rownames(se)[1:5], taxonomy=tax)

Description

This function filters a SummarizedExperiment object to keep only rows (features/strains) that are present in a specified number of samples. The same rows are filtered from all assays in the SummarizedExperiment.

Usage

filter_by_presence(se, min_nonzero = 10)

Arguments

Value

A filtered SummarizedExperiment object with only the rows that meet the criteria.

Examples

example_path <- system.file("extdata", "example_sylph_profile.tsv.gz", 
package = "strainspy")
se <- read_sylph(example_path)

# Filter to keep only rows with at least 20 non-zero entries 
filtered_se <- filter_by_presence(se, min_nonzero = 20)

Description

Gets the ordered list of fitted contig names in a strainspy_fit

Usage

getContigNames(object)

## S4 method for signature 'strainspy_fit'
getContigNames(object)

Arguments

Value

A vector of contig names

Examples

meta <- read.csv(system.file("extdata", "example_metadata.csv.gz", 
package = "strainspy"))
meta$Case_status <- factor(meta$Case_status)
se <- read_sylph(system.file("extdata", "example_sylph_profile.tsv.gz", 
package = "strainspy"), meta_data = meta)
fit <- glmZiBFit(filter_by_presence(se, min_nonzero = 30)[1:10, ], 
~ Case_status + Age_at_collection, nthreads = 1, return_vcov = FALSE)
head(getContigNames(fit))

Description

Gets the ordered list of fitted genomes in a strainspy_fit

Usage

getGenomes(object)

## S4 method for signature 'strainspy_fit'
getGenomes(object)

Arguments

Value

A vector of fitted genomes

Examples

meta <- read.csv(system.file("extdata", "example_metadata.csv.gz", 
package = "strainspy"))
meta$Case_status <- factor(meta$Case_status)
se <- read_sylph(system.file("extdata", "example_sylph_profile.tsv.gz", 
package = "strainspy"), meta_data = meta)
fit <- glmZiBFit(filter_by_presence(se, min_nonzero = 30)[1:10, ], 
~ Case_status + Age_at_collection, nthreads = 1, return_vcov = FALSE)
head(getGenomes(fit))

Description

Returns the full prior data.frame stored in a strainspy_priors object, including both fixed-effect and zero-inflation rows. See extract_fixef_priors() for the fixed-effect rows only.

Usage

getPriorsDf(object)

## S4 method for signature 'strainspy_priors'
getPriorsDf(object)

Arguments

Value

A data.frame with columns prior, class and coef.

Examples

meta <- read.csv(system.file("extdata", "example_metadata.csv.gz",
package = "strainspy"))
meta$Case_status <- factor(meta$Case_status)
se <- read_sylph(system.file("extdata", "example_sylph_profile.tsv.gz",
package = "strainspy"), meta_data = meta)
pr <- define_priors(se, ~ Case_status + Age_at_collection, method = "preset_weak")
getPriorsDf(pr)

Description

Returns the residuals from a strainspy_fit object.

Usage

getResiduals(object)

## S4 method for signature 'strainspy_fit'
getResiduals(object)

Arguments

Value

A DataFrame containing residuals.

Examples

meta <- read.csv(system.file("extdata", "example_metadata.csv.gz", 
package = "strainspy"))
meta$Case_status <- factor(meta$Case_status)
se <- read_sylph(system.file("extdata", "example_sylph_profile.tsv.gz", 
package = "strainspy"), meta_data = meta)
fit <- glmZiBFit(filter_by_presence(se, min_nonzero = 30)[1:10, ], 
~ Case_status + Age_at_collection, nthreads = 1, return_vcov = FALSE)
class(getResiduals(fit))

Description

Gets the ordered list of samples used in strainspy_fit

Usage

getSamples(object)

## S4 method for signature 'strainspy_fit'
getSamples(object)

Arguments

Value

A vector of fitted sample names

Examples

meta <- read.csv(system.file("extdata", "example_metadata.csv.gz", 
package = "strainspy"))
meta$Case_status <- factor(meta$Case_status)
se <- read_sylph(system.file("extdata", "example_sylph_profile.tsv.gz", 
package = "strainspy"), meta_data = meta)
fit <- glmZiBFit(filter_by_presence(se, min_nonzero = 30)[1:10, ], 
~ Case_status + Age_at_collection, nthreads = 1, return_vcov = FALSE)
head(getSamples(fit))

Description

Returns the zero-inflation coefficients from a strainspy_fit object.

Usage

getZICoefficients(object)

## S4 method for signature 'strainspy_fit'
getZICoefficients(object)

Arguments

Value

A DataFrame containing zero-inflation coefficients.

Examples

meta <- read.csv(system.file("extdata", "example_metadata.csv.gz", 
package = "strainspy"))
meta$Case_status <- factor(meta$Case_status)
se <- read_sylph(system.file("extdata", "example_sylph_profile.tsv.gz", 
package = "strainspy"), meta_data = meta)
fit <- glmZiBFit(filter_by_presence(se, min_nonzero = 30)[1:10, ], 
~ Case_status + Age_at_collection, nthreads = 1, return_vcov = FALSE)
class(getZICoefficients(fit))

Description

This function fits an ordinal beta regression model using the glmmTMB package. It takes a SummarizedExperiment object as input, along with a user-defined formula, and fits an ordinal beta regression model on the assay data.

Usage

glmObFit(
  se,
  design,
  nthreads = 1L,
  scale_continuous = TRUE,
  MAP_prior = "preset_weak",
  family = glmmTMB::ordbeta(),
  BPPARAM = NULL,
  progress = TRUE
)

Arguments

Details

For robust model fitting, it is recommended to use maximum a posteriori (MAP) estimation. Empirical Bayes MAP priors can be computed using strainspy::compute_eb_priors(). They can also be manually defined using strainspy::define_priors(). The resulting strainspy_priors object can be passed to the MAP_prior argument.

Alternatively, a character string ("preset_weak" or "preset_strong") or a valid prior data.frame (see ?glmmTMB::priors) can be provided instead. This will result in a call to strainspy::define_priors() to generate priors. To disable MAP estimation and fit by maximum likelihood without priors, set MAP_prior = NULL.

Variance-covariance matrices are not returned by this fitter: the vcov slot of the result is NULL, so the fit cannot be used with comp_ani_diff_and_posthoc_test().

Value

A strainspy_fit object with the following components:

Examples

library(strainspy)

example_meta_path <- system.file("extdata", "example_metadata.csv.gz", 
package = "strainspy")
example_meta <- read.csv(example_meta_path)
example_path <- system.file("extdata", "example_sylph_profile.tsv.gz", 
package = "strainspy")
se <- read_sylph(example_path, example_meta)
se <- filter_by_presence(se)

design <- as.formula(" ~ Case_status + Age_at_collection")

fit <- glmObFit(se[1:10,],  design, family=glmmTMB::ordbeta())

Description

This function fits a quasibinomial regression model using stats::glm().

Usage

glmQBFit(
  se,
  design,
  nthreads = 1L,
  scale_continuous = TRUE,
  BPPARAM = NULL,
  progress = TRUE
)

Arguments

Details

It takes a SummarizedExperiment object as input, along with a user-defined formula containing only fixed effect parameters, and fits a quasibinomial regression model on the assay data.

Variance-covariance matrices are not returned by this fitter: the vcov slot of the result is NULL, so the fit cannot be used with comp_ani_diff_and_posthoc_test().

Value

A strainspy_fit object with the following components:

Examples

library(strainspy)

example_meta_path <- system.file("extdata", "example_metadata.csv.gz", 
package = "strainspy")
example_meta <- read.csv(example_meta_path)
example_path <- system.file("extdata", "example_sylph_profile.tsv.gz",
package = "strainspy")
se <- read_sylph(example_path, example_meta)
se <- filter_by_presence(se)

design <- as.formula(" ~ Case_status")

fit <- glmQBFit(se[1:10,], design, nthreads=2)

Description

This function fits a zero-inflated beta regression model using either glmmTMB or gamlss packages. It takes a SummarizedExperiment object as input, along with a user-defined formula, and fits a zero-inflated beta regression model on the assay data.

Usage

glmZiBFit(
  se,
  design,
  nthreads = 1,
  scale_continuous = TRUE,
  MAP_prior = "preset_weak",
  BPPARAM = NULL,
  method = "glmmTMB",
  return_vcov = TRUE,
  progress = TRUE
)

Arguments

Value

A strainspy_fit object with the following components:

Examples

library(strainspy)

example_meta_path <- system.file("extdata", "example_metadata.csv.gz", 
package = "strainspy")
example_meta <- read.csv(example_meta_path)
example_path <- system.file("extdata", "example_sylph_profile.tsv.gz", 
package = "strainspy")
se <- read_sylph(example_path, example_meta)
se <- filter_by_presence(se)

design <- as.formula(" ~ Case_status")

fit <- glmZiBFit(se[1:10,], design, nthreads=2)

Description

This function performs multiple testing adjustment on the results of a strainspy_fit object at each taxonomic level, using either the Harmonic Mean P-value (HMP) or Bonferroni correction.

Usage

hadjust(object, coef = 2, method = "HMP", taxonomy, index_range = FALSE)

Arguments

Value

A tibble with one row per taxon and level, with columns Level, Model, Name, the number of strains (strain_count for HMP, sequence_count for Bonferroni), mean_coefficient and the adjusted p-value p_adjust.

Examples

meta <- read.csv(system.file("extdata", "example_metadata.csv.gz", 
package = "strainspy"))
meta$Case_status <- factor(meta$Case_status)
se <- read_sylph(system.file("extdata", "example_sylph_profile.tsv.gz", 
package = "strainspy"), meta_data = meta)
se <- filter_by_presence(se, min_nonzero = 30)
tax <- read_taxonomy(system.file("extdata", "example_taxonomy.tsv.gz",
package = "strainspy"))
fit <- glmZiBFit(se[1:10, ], ~ Case_status + Age_at_collection, 
nthreads = 1, return_vcov = FALSE)
out <- hadjust(fit, coef = 2, taxonomy = tax)
head(out)

Description

Merge MetaPhlAn >=4.0 profiles into a single file to use with read_metaphlan() and generate the taxonomy file required for downstream analysis.

Usage

merge_metaphlan_files(
  metaphlan_files,
  output_folder,
  output_filename = "metaphlan",
  sgb_database = TRUE
)

Arguments

Value

Invisibly returns NULL. The function writes two files to output_folder: a merged profile table (⁠*_merged.tsv.gz⁠) and a taxonomy table (⁠*_taxonomy.tsv.gz⁠).

Examples

files <- character(0)
is.character(files)

Description

Merge the outputs of ⁠sourmash search⁠ or ⁠sourmash gather⁠ into a single file. The output of this function is suitable to use in read_sourmash()

Usage

merge_sourmash_files(
  sourmash_files,
  manifest_file = NULL,
  output,
  strip_unusable = TRUE
)

Arguments

Details

Refer to sourmash documentation for more information.

Value

Invisibly returns NULL. The merged sourmash table is written to output.

Examples

# Split the bundled sourmash output into two per-query files, then merge them.
hits <- data.table::fread(
  system.file("extdata", "example_sourmash.csv.gz", package = "strainspy")
)
queries <- head(unique(hits$query_name), 2)
parts <- vapply(queries, function(q) {
  f <- tempfile(fileext = ".csv")
  data.table::fwrite(hits[hits$query_name == q, ], f)
  f
}, character(1))

merged <- tempfile(fileext = ".csv.gz")
merge_sourmash_files(parts, output = merged)

sm <- read_sourmash(merged)
dim(sm)

Description

Merge the outputs of ⁠sylph query⁠ or ⁠sylph profile⁠ into a single file. The output of this function is suitable to use in read_sylph()

Usage

merge_sylph_files(sylph_files, output)

Arguments

Value

Invisibly returns output. The merged sylph table is written to that path.

Examples

# Split the bundled profile into two per-sample files, then merge them back.
profile <- data.table::fread(
  system.file("extdata", "example_sylph_profile.tsv.gz", package = "strainspy")
)
samples <- head(unique(profile$Sample_file), 2)
parts <- vapply(samples, function(s) {
  f <- tempfile(fileext = ".tsv")
  data.table::fwrite(profile[profile$Sample_file == s, ], f, sep = "\t")
  f
}, character(1))

merged <- tempfile(fileext = ".tsv.gz")
merge_sylph_files(parts, merged)

# The merged file carries one header and both samples
se <- read_sylph(merged)
dim(se)

Description

This function adds or modifies metadata in a SummarizedExperiment object.

Usage

modify_metadata(se, meta_data, replace = FALSE, convert = TRUE)

Arguments

Details

If convert = TRUE, character columns with 5 or fewer unique (non-NA) values will be automatically converted to factors, with levels ordered alphabetically. For meaningful association testing, however, factor levels should reflect the intended reference level. It is recommended to manually convert such variables to factors with explicitly set levels before calling modify_metadata().

Value

Updated SummarizedExperiment object.

Examples

# Add or modify metadata subsequently
  example_meta_path <- system.file("extdata", "example_metadata.csv.gz", 
  package = "strainspy")
  example_meta <- read.csv(example_meta_path)
  example_path <- system.file("extdata", "example_sylph_profile.tsv.gz", 
  package = "strainspy")
  se <- read_sylph(example_path)
  se <- modify_metadata(se, example_meta)

Description

Generate a violin or box plot for a subset of contigs for a categorical phenotype.

Usage

plot_ani_dist(
  se,
  phenotype,
  contigs,
  facet_phenotype = NULL,
  contig_names = NULL,
  show_zero_plot = TRUE,
  show_points = TRUE,
  plot = TRUE,
  plot_type = "box"
)

Arguments

Value

If plot = TRUE, returns a ggplot object (or a combined patchwork plot when show_zero_plot = TRUE). If plot = FALSE, returns a long-format table used to generate the ANI distribution plot.

Examples

library(strainspy)

example_meta_path <- system.file("extdata", "example_metadata.csv.gz", package = "strainspy")
example_meta <- read.csv(example_meta_path)
# Change variable from character to factor
example_meta$Case_status = factor(example_meta$Case_status)

example_path <- system.file("extdata", "example_sylph_profile.tsv.gz", package = "strainspy")
se <- read_sylph(example_path, example_meta)

plot_ani_dist(se, phenotype = "Case_status", contigs = rownames(se)[1:5])

Description

Generate a histogram for a single contig for a categorical phenotype.

Usage

plot_histogram(
  se,
  phenotype,
  contig,
  bins = 250,
  drop_zeros = FALSE,
  show_quantiles = TRUE,
  fit_spline = FALSE,
  fit_control = 100,
  separate_facets = TRUE
)

Arguments

Value

A ggplot histogram object. A summary table of per-group statistics is attached as attr(plot, "summary_table").

Examples

library(strainspy)

example_meta_path <- system.file("extdata", "example_metadata.csv.gz", 
package = "strainspy")
example_meta <- read.csv(example_meta_path)
# Change variable from character to factor
example_meta$Case_status = factor(example_meta$Case_status)

example_path <- system.file("extdata", "example_sylph_profile.tsv.gz", 
package = "strainspy")
se <- read_sylph(example_path, example_meta)

plot_histogram(se, phenotype = "Case_status", contig = rownames(se)[1])

Description

This function can generate multiple Manhattan plots from a strainspy_fit object.

Usage

plot_manhattan(
  object,
  coef = 2,
  taxonomy = NULL,
  aggregate_by_taxa = NULL,
  method = NULL,
  alpha = 0.05,
  tax_levels = c("Phylum", "Genus", "Species", "Strain"),
  plot = TRUE
)

Arguments

Details

The generated plot will be more informative if taxonomic information is provided. Two types of plots can be generated with taxonomic information.

1. Set aggregate_by_taxa = TRUE to aggregate p-values across taxonomic levels specified by tax_levels. This generates a stacked segment Manhattan plot. In this plot, -log10 p-values adjusted using method will be plotted in taxonomic ordering. Significance threshold will be alpha.

2. Set aggregate_by_taxa = FALSE to generate a traditional Manhattan plot for all strains coloured by phylum. Unadjusted p-values will be plotted and Bonferroni significance thresholds (0.05, 0.01) will be shown. Strains will be in taxonomic ordering and a tree showing the taxonomic relationship between strains will be plotted below the Manhattan plot.

If taxonomic information is not provided (taxonomy = NULL), a plot with -log10 p-values adjusted using method will be generated. Strains will be plotted in alphabetical order. Significance threshold will be alpha.

If Zero-Inflated Beta regression was used to model object, separate facets will be used for zero-inflated and beta components of the model.

Value

A ggplot object showing the Manhattan plot if plot = TRUE, else a tibble with the plot data.

Examples

meta <- read.csv(system.file("extdata", "example_metadata.csv.gz", 
package = "strainspy"))
meta$Case_status <- factor(meta$Case_status)
se <- read_sylph(system.file("extdata", "example_sylph_profile.tsv.gz", 
package = "strainspy"), meta_data = meta)
se <- filter_by_presence(se, min_nonzero = 30)
tax <- read_taxonomy(system.file("extdata", "example_taxonomy.tsv.gz", 
package = "strainspy"))
fit <- glmZiBFit(se[1:10, ], ~ Case_status + Age_at_collection, 
nthreads = 1, return_vcov = FALSE)
mh_data <- plot_manhattan(fit, taxonomy = tax, plot = FALSE)
head(mh_data)

Description

Performs classical multidimensional scaling (MDS) on a SummarizedExperiment object and plots the result, optionally colored by a phenotype.

Usage

plot_mds(
  se,
  coef = 2,
  taxonomy = NULL,
  palette = c("#DAA520", "#3CB371", "#483D8B"),
  plot = TRUE
)

Arguments

Value

A ggplot object showing the MDS plot, or a tibble of MDS coordinates if plot = FALSE.

Examples

example_meta_path <- system.file("extdata", "example_metadata.csv.gz", 
package = "strainspy")
example_meta <- read.csv(example_meta_path)
example_path <- system.file("extdata", "example_sylph_profile.tsv.gz", 
package = "strainspy")
se <- read_sylph(example_path, example_meta)
plot_mds(se, coef = 2)

Description

Perform PCA and Generate a Biplot

Usage

plot_pca_biplot(se, phenotype, plot = TRUE)

Arguments

Details

This function performs Principal Component Analysis (PCA) on a SummarizedExperiment object and generates a PCA biplot colored by a specified phenotype.

Value

If plot = TRUE, returns a ggplot2 PCA plot. If plot = FALSE, returns a tibble with PCA results.

Examples

example_meta_path <- system.file("extdata", "example_metadata.csv.gz", 
package = "strainspy")
example_meta <- read.csv(example_meta_path)
example_path <- system.file("extdata", "example_sylph_profile.tsv.gz", 
package = "strainspy")
se <- read_sylph(example_path, example_meta)
plot_pca_biplot(se, phenotype = 3)

Description

Plot the bootstrap distributions behind a strainspy_priors object

Usage

plot_prior_bootstrap(object, prior)

## S4 method for signature 'strainspy_priors'
plot_prior_bootstrap(object, prior)

Arguments

Value

A patchwork of ggplot2 histograms.

Examples

methods::hasMethod("plot_prior_bootstrap", "strainspy_priors")

Description

This function creates a volcano plot from the results of a strainspy model fit, highlighting the relationship between effect size (coefficients) and statistical significance (-log10 p-values).

Usage

plot_volcano(
  object,
  coef = 2,
  alpha = 0.05,
  method = "holm",
  palette = rev(viridis::viridis(3)),
  label = FALSE,
  contig_names = NULL,
  plot = TRUE
)

Arguments

Value

A ggplot2 object if plot = TRUE, otherwise a dataframe containing the relevant statistics.

Examples

example_meta_path <- system.file("extdata", "example_metadata.csv.gz", 
package = "strainspy")
example_meta <- read.csv(example_meta_path)
example_path <- system.file("extdata", "example_sylph_profile.tsv.gz", 
package = "strainspy")
se <- read_sylph(example_path, example_meta)
se <- filter_by_presence(se)
design <- as.formula(" ~ Case_status + Sex")
fit <- glmZiBFit(se[1:10,], design, nthreads=2)
plot_volcano(fit, coef = 2, label = TRUE)

Description

Formats strain-level data and optional covariates into a data.frame object ready for modeling with caret. The function extracts selected contigs from the sy assay matrix, combines them with the specified outcome and optional covariates from the metadata, and performs basic validation and cleaning. Metadata is pulled from SummarizedExperiment::colData(sy), and generally assumes strainspy::modify_metadata() was used upstream.

Usage

prep_for_prediction(
  sy,
  outcome,
  contigs,
  covariates = NULL,
  use_genome_names = TRUE
)

Arguments

Details

Character variables are automatically converted to factors. Numeric variables are returned unchanged. For pre processing, refer to caret::train().

Value

A data.frame containing the specified outcome variable, optional covariates, and selected contig features. The object is ready for direct use in caret functions such as caret::train().

Examples

library(strainspy)

example_meta_path <- system.file("extdata", "example_metadata.csv.gz", 
package = "strainspy")
example_meta <- read.csv(example_meta_path)
example_path <- system.file("extdata", "example_sylph_profile.tsv.gz", 
package = "strainspy")
sy <- read_sylph(example_path, example_meta)
contigs <- rownames(sy)[1:10]
df <- prep_for_prediction(
  sy = sy,
  outcome = "Case_status",
  contigs = contigs,
  covariates = c("Sex", "BMI")
)

Description

This function reads a merged MetaPhlAn profile output and creates a SummarizedExperiment object. An optional metadata table can be provided and added directly to colData via an internal call to modify_metadata().

Usage

read_metaphlan(
  file_path,
  meta_data = NULL,
  variable = "relative_abundance",
  clean_names = TRUE
)

Arguments

Value

A SummarizedExperiment object with the following components:

Note

Since MetaPhlAn profiles are usually generated per sample, first use merge_metaphlan_files() to merge and generate the input for this function.

If metadata is provided, it must meet the input requirements described in modify_metadata(). If appending metadata fails due to pre processing requirements, this function will issue a detailed warning and still return the SummarizedExperiment object without metadata. In that case, apply necessary fixes and subsequently call modify_metadata().

Examples

# Read a merged metaphlan profile file into a SummarizedExperiment object
  example_path <- system.file("extdata", "metaphlan_merged.tsv.gz", 
  package = "strainspy")
  mp <- read_metaphlan(example_path)

Description

This function reads a ⁠sourmash gather⁠ or ⁠sourmash search⁠ -like output CSV file and creates a SummarizedExperiment object. An optional metadata table can be provided and added directly to colData via an internal call to modify_metadata().

Usage

read_sourmash(
  file_path,
  meta_data = NULL,
  variable = "match_containment_ani",
  clean_names = TRUE
)

Arguments

Value

A SummarizedExperiment object with the following components:

Note

Some sourmash subcommands, like gather, are designed to take exactly one query file. To merge and save a CSV file suitable to use with this function, use merge_sourmash_files().

If metadata is provided, it must meet the input requirements described in modify_metadata(). If appending metadata fails due to pre processing requirements, this function will issue a detailed warning and still return the SummarizedExperiment object without metadata. In that case, apply necessary fixes and subsequently call modify_metadata().

Examples

# Read a merged sourmash file into a SummarizedExperiment object
  example_path <- system.file("extdata", "example_sourmash.csv.gz", 
  package = "strainspy")
  sm <- read_sourmash(example_path)

Description

This function reads a Sylph output file (either query or profile), parses the data, and returns a SummarizedExperiment object. An optional metadata table can be provided and added directly to colData via an internal call to modify_metadata().

Usage

read_sylph(
  file_path,
  meta_data = NULL,
  variable = "Adjusted_ANI",
  min_identity = 95,
  clean_names = TRUE
)

Arguments

Details

If metadata is provided, it must meet the input requirements described in modify_metadata(). If appending metadata fails due to pre processing requirements, this function will issue a detailed warning and still return the SummarizedExperiment object without metadata. In that case, apply necessary fixes and subsequently call modify_metadata().

Value

A SummarizedExperiment object with the following components:

Examples

# Read a Sylph file 
  example_path <- system.file("extdata", "example_sylph_profile.tsv.gz", 
  package = "strainspy")
  sy <- read_sylph(example_path)

Description

This function reads a headerless, two-column, tab-separated taxonomy file: genome accession and GTDB-style lineage string (⁠d__...;p__...;...;s__...⁠), as distributed with the sylph GTDB databases.

Usage

read_taxonomy(file_path)

Arguments

Value

A data.table with columns Genome, Taxonomy, Domain, Phylum, Class, Order, Family, Genus, Species and index (the row order after sorting by lineage).

Examples

tax <- read_taxonomy(system.file("extdata", "example_taxonomy.tsv.gz", 
package = "strainspy"))
head(tax)

Description

Displays a summary of a strainspy_fit object, including model design, convergence, and coefficient matrix dimensions.

Usage

## S4 method for signature 'strainspy_fit'
show(object)

Arguments

Value

Invisibly returns NULL after printing a compact summary.

Description

Show method for strainspy_priors objects

Usage

## S4 method for signature 'strainspy_priors'
show(object)

Arguments

Value

Invisibly returns the object

Description

A wrapper to run a full strainspy workflow with explicit input validation.

Usage

strainspy(
  meta_path,
  sylph_path,
  design_formula,
  coef,
  min_nonzero = NULL,
  taxonomy_path = NULL,
  output_path = NULL,
  alpha = 0.05,
  p_adjust_method = "holm",
  model_method = "glmmTMB",
  scale_continuous = TRUE,
  nthreads = 2,
  return_vcov = FALSE,
  attach_taxonomy = TRUE,
  overwrite = TRUE,
  progress = TRUE
)

Arguments

Value

A data frame of top hits. If taxonomy_path is supplied and attach_taxonomy = TRUE, taxonomy columns are appended.

Examples

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")

# Decision 1: choose predictors for the design formula
# Decision 2: choose which coefficient index to report in top_hits()
# For ~ Case_status + Age_at_collection, coef = 2 corresponds to Case_status

out <- strainspy(
  meta_path = meta_path,
  sylph_path = sylph_path,
  design_formula = "~ Case_status + Age_at_collection",
  coef = 2,
  taxonomy_path = taxonomy_path,
  output_path = tempfile(fileext = ".csv"),
  nthreads = 2,
  alpha = 1
)
head(out)

Description

Getters for the per-feature statistics and fitting metadata stored in a strainspy_fit object. Use these rather than reaching into the object with @.

Usage

getRowData(object)

## S4 method for signature 'strainspy_fit'
getRowData(object)

getCoefficients(object)

## S4 method for signature 'strainspy_fit'
getCoefficients(object)

getStdErrors(object)

## S4 method for signature 'strainspy_fit'
getStdErrors(object)

getPValues(object)

## S4 method for signature 'strainspy_fit'
getPValues(object)

getZIStdErrors(object)

## S4 method for signature 'strainspy_fit'
getZIStdErrors(object)

getZIPValues(object)

## S4 method for signature 'strainspy_fit'
getZIPValues(object)

getVcov(object)

## S4 method for signature 'strainspy_fit'
getVcov(object)

getZIVcov(object)

## S4 method for signature 'strainspy_fit'
getZIVcov(object)

getDesign(object)

## S4 method for signature 'strainspy_fit'
getDesign(object)

getFamily(object)

## S4 method for signature 'strainspy_fit'
getFamily(object)

getFitCall(object)

## S4 method for signature 'strainspy_fit'
getFitCall(object)

isScaleContinuous(object)

## S4 method for signature 'strainspy_fit'
isScaleContinuous(object)

getConvergence(object)

## S4 method for signature 'strainspy_fit'
getConvergence(object)

Arguments

Value

• getRowData(): a DataFrame of feature annotations, one row per feature.

• getCoefficients(), getStdErrors(), getPValues(): DataFrames of per-feature estimates, one column per model term.

• getZIStdErrors(), getZIPValues(): the zero-inflation counterparts, or NULL when the model has no zero-inflation component.

• getVcov(), getZIVcov(): a List of per-feature variance-covariance matrices, or NULL when they were not stored.

• getDesign(): the model formula.

• getFamily(): the model family object.

• getFitCall(): the call that produced the fit.

• isScaleContinuous(): TRUE if continuous covariates were z-scored before fitting.

• getConvergence(): a logical vector, one entry per feature, TRUE where the model converged.

Examples

meta <- read.csv(system.file("extdata", "example_metadata.csv.gz",
package = "strainspy"))
meta$Case_status <- factor(meta$Case_status)
se <- read_sylph(system.file("extdata", "example_sylph_profile.tsv.gz",
package = "strainspy"), meta_data = meta)
fit <- glmZiBFit(filter_by_presence(se, min_nonzero = 30)[1:10, ],
~ Case_status + Age_at_collection, nthreads = 1, return_vcov = FALSE)
head(getCoefficients(fit))
getDesign(fit)
isScaleContinuous(fit)

Description

The strainspy_fit class stores the results of a strainspy model fit.

Slots

row_data

A DataFrame containing row data.

priors

A strainspy_priors object, or NULL for an unpenalised maximum-likelihood fit.

scale_continuous

A logical indicating whether all numeric columns were z-score standardized

coefficients

A DataFrame containing model coefficients.

std_errors

A DataFrame containing standard errors of the coefficients.

p_values

A DataFrame containing p-values for the coefficients.

vcov

A List of per-feature variance-covariance matrices, or NULL.

zi_coefficients

A DataFrame containing zero-inflation coefficients.

zi_std_errors

A DataFrame containing zero-inflation standard errors.

zi_p_values

A DataFrame containing zero-inflation p-values.

zi_vcov

A List of per-feature zero-inflation variance-covariance matrices, or NULL.

residuals

A DataFrame containing residuals from the model.

convergence

A logical indicating whether model fitting converged.

design

The model formula.

call

The matched call of the model.

family

The family of the fit.

Description

The strainspy_priors class is designed to store MAP priors used for ordinal and zero-inflated beta model fitting using glmmTMB.

Slots

method

Character. Prior estimation method used ("preset_weak", "preset_strong", "empirical", "manual").

priors_df

Data frame. Data frame containing prior parameters.

boot_fixef

Numeric matrix. Bootstrap SD values for fixed effect priors

boot_fixef_ZI

Numeric matrix. Bootstrap SD values for zero inflated fixed effect priors

boot_disp

Numeric matrix. Bootstrap mean and SD values for the dispersion term of beta component.

design

The formula with the effect terms used

call

The matched call of the model.

Description

Provides a high-level summary of a strainspy_fit object: its call, the first rows of the coefficients and zero-inflation coefficients, and whether residuals are available.

Usage

## S4 method for signature 'strainspy_fit'
summary(object)

Arguments

Value

A list summarizing the strainspy_fit object.

Description

Adjusts the p-values of one coefficient of a strainspy_fit for multiple testing with p.adjust() and returns the features that pass alpha, sorted by adjusted p-value. For taxonomy-aware adjustment see hadjust().

Usage

top_hits(fit, coef = 2, method = "holm", alpha = 0.05)

Arguments

Value

A tibble with the top hits, sorted by adjusted p-value.

Examples

meta <- read.csv(system.file("extdata", "example_metadata.csv.gz",
 package = "strainspy"))
meta$Case_status <- factor(meta$Case_status)
se <- read_sylph(system.file("extdata", "example_sylph_profile.tsv.gz", 
package = "strainspy"), meta_data = meta)
fit <- glmZiBFit(filter_by_presence(se, min_nonzero = 30)[1:10, ], 
~ Case_status + Age_at_collection, nthreads = 1, return_vcov = FALSE)
top_hits(fit, coef = 2, alpha = 1)