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 distance matrix.

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)

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

Internal helper used by the fitting functions to check that progress is a single, non-missing logical value.

Usage

check_progress(progress)

Arguments

Value

progress, invisibly.

Description

read_sylph(), read_metaphlan() and read_sourmash() all reduce Sample_file to a bare basename when clean_names = TRUE. Two input paths that differ only by directory (runA/S1.fastq.gz and runB/S1.fastq.gz) therefore collapse to the same name. Duplicated names are not rejected by SummarizedExperiment, and every downstream match() on sample or feature names would then silently attach the wrong metadata to the wrong column, so fail here instead.

Usage

check_unique_names(x, what, clean_names = FALSE)

Arguments

Value

x, invisibly, if every name is unique.

Description

Automatically clean contig names for plots

Usage

clean_contig_names(contig_names, max_length = 50, return_map = FALSE)

Arguments

Value

A data.frame or original and shortened names, or a vector of shortened name in the same order

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. Same applies for 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 fit@design. However, if the formula contains multiple covariates, this can cause emmeans to behave unexpectedly. For example, if fit@design = ~ 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 tibble 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 commetn 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_disperion_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

Estimate beta diversity from ANI

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. Asmore 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

The fitting backends signal failure in two different ways. glmObFit(), glmQBFit() and abundanceFit() return NULL for a feature that could not be fitted, so the failure shows up as a zero-length entry. Both glmZiBFit() backends instead return a fully-populated list whose coefficient matrices are NA and whose log_likelihood is NA, so the zero-length test never fires there and the NA log-likelihood is the only signal.

Usage

failed_fit_index(results)

Arguments

Details

Detecting only one of the two leaves failed fits in the results, and applying the log_likelihood test blindly is worse: glmQBFit() does not report a log-likelihood at all, so every successful fit would look like a failure. This helper handles both protocols by treating a log-likelihood that is absent from every result as "not reported" rather than "all failed".

It deliberately uses vapply() rather than unlist(). unlist() silently drops the NULL log-likelihood of a zero-length entry, which shortens the vector and shifts every later index, so the returned positions would no longer line up with results or with rowData().

Value

An increasing integer vector of positions in results that should be dropped, or integer(0) if every fit succeeded.

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 >20 non-zero entries 
filtered_se <- filter_by_presence(se, min_nonzero = 20)

Description

This function is copied from filter_by_presence(). Additionally, it rescales abundances such that each sample sums back to 100. Since this behaviour is not useful for association testing, this function is not exported.

Usage

filter_by_presence_and_rescale(se, min_nonzero = 10, rescale_abundance = FALSE)

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 10 non-zero entries in the assay
filtered_se <- filter_by_presence(se)

# View the filtered object
filtered_se

Description

FThis function filters a SummarizedExperiment object to keep only rows (features/strains) present in at least min_subjects individuals and at least min_timepoints timepoints per subject.

Usage

filter_by_presence_longitudinal(
  se,
  min_subjects = 3,
  min_timepoints = 2,
  subject_col = "subject",
  time_col = "time"
)

Arguments

Value

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

Description

Fits a simple linear model for functions without random effects using stats::lm()

Usage

fit_lm(se_subset, col_data, combined_formula)

Arguments

Value

A list with model coefficients, zero-inflation coefficients, residuals, log-likelihood, and convergence status.

Description

Fits a linear model for a single feature in an assay using lmerTest::lmer()

Usage

fit_lmer(se_subset, col_data, combined_formula)

Arguments

Value

A list with model coefficients, zero-inflation coefficients, residuals, log-likelihood, and convergence status.

Description

Fits a logistic regression for a single feature in an assay using glmmTMB.

Usage

fit_logit_model(se_subset, col_data, design, fixed_priors, min_identity)

Arguments

Value

A list with model coefficients, zero-inflation coefficients, residuals, log-likelihood, and convergence status.

Description

Fits a zero-inflated beta regression for a single feature in an assay using glmmTMB.

Usage

fit_ob_model(se_subset, col_data, combined_formula, fixed_priors, family)

Arguments

Value

A list where each element contains per-feature fitted coefficients, residuals, and convergence status.

Description

Fits a quasibinomial regression for a single feature in an assay using stats::glm().

Usage

fit_qb_model(se_subset, col_data, combined_formula)

Arguments

Value

A list where each element contains per-feature fitted coefficients, residuals, and convergence status.

Description

Fits a zero-inflated beta regression for a single feature in an assay using glmmTMB.

Usage

fit_zero_inflated_beta(
  se_subset,
  col_data,
  combined_formula,
  design,
  fixed_priors
)

Arguments

Value

A list with model coefficients, zero-inflation coefficients, residuals, log-likelihood, and convergence status.

Description

Fits a zero-inflated beta regression for a single feature in an assay using gamlss.

Usage

fit_zero_inflated_beta_gamlss(
  se_subset,
  col_data,
  combined_formula,
  design,
  fixed_priors,
  progress = TRUE
)

Arguments

Value

A list with model coefficients, zero-inflation coefficients, residuals, log-likelihood, and convergence status.

Description

This function generates a vector of colours for plotting, prioritising a fixed palette of 21 distinct colours. If more than 21 colours are required, additional colours are generated using grDevices::rainbow().

Usage

get_colors(n)

Arguments

Value

A character vector of n colour hex codes.

Description

Generate a confusion matrix for a fitted model

Usage

get_confusion_mx(top_hits, gt_contigs, all_contigs, print_cm = TRUE)

Arguments

Value

Confusion Matrix.

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

TODO: Implement return of vcov info. Patched for compatibility (returns NULL)

Usage

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

Arguments

Details

This function fits a ordinalbeta 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.

For robust model fitting, it is recommended to use maximum a posteriori (MAP) estimation. Emprical 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 perform standard REML fitting, set MAP_prior = NULL.

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_continous = 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.

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 betaGLM object integrating taxonomic information. Various correction methods are available, including Harmonic Mean P-value (HMP), Bonferroni and Benjamini-Hochberg (BH).

Usage

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

Arguments

Value

A tibble with original and adjusted p-values.

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 NULL. The merged sylph table is written to output.

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 less than 5 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

nobars_ copied from https://github.com/bbolker/reformulas This function can work with formulas like ~cohort + timepoint + (1 + timepoint | patientID), above fails We can either use it like this, or add reformulas (or lme4) as a dependency

Usage

nobars_(term)

Arguments

Value

A formula without mixed components

Description

This function reduces all non-zero ANI values by eps

Usage

offset_ANI(x, eps = 0.01)

Arguments

Value

Adjusted ANI vector

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 = "holm",
  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 data.frame with 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

This function is copied from plot_manhattan()

Usage

plot_manhattan_gt(
  object,
  coef = 2,
  taxonomy = NULL,
  method = "HMP",
  alpha = 0.05,
  levels = c("Phylum", "Genus", "Species", "Strain"),
  colour_level = "Phylum",
  plot = TRUE,
  ground_truth = NULL
)

Arguments

Value

A ggplot object showing the Manhattan plot.

Description

Generate a conventional Manhattan plot for the strain all strains.

Usage

plot_manhattan_tree(object, taxonomy, coef = 2)

Arguments

Value

A ggplot object showing the Manhattan plot combined with the taxonomic tree

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 method for strainspy_priors objects

Usage

plot_prior_bootstrap(object, prior)

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

Arguments

Value

Invisibly returns the object

Examples

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

Description

This function creates a volcano plot from the results of a betaGLM analysis, 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 SummarizedExperiement::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 create 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 create 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 of 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 taxonomy file in Sylph format.

Usage

read_taxonomy(file_path)

Arguments

Value

A data.table containing parsed taxonomy information.

Examples

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

Description

This function applies a simple linear transformation to simulate a treatment effect.

Usage

rescale_beta(x, beta = 0.98, zi = 0.1)

Arguments

Value

a vector of transformed values

Description

Ensures MAP_prior is interpreted correctly and returns a valid strainspy_priors object or NULL (for unpenalised ML fits).

Usage

resolve_priors(MAP_prior, se, design)

Arguments

Value

A strainspy_priors object or NULL.

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

The strainspy_fit class is designed to store the results of strainspy

Slots

row_data

A DataFrame containing row data.

priors

A strainspy_priors object or NULL for REML.

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 DataFrame containing the variance-covariance matrix.

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 DataFrame containing the zero-inflated variance-covariance matrix.

residuals

A DataFrame containing residuals from the model.

convergence

A logical indicating whether model fitting converged.

design

A matrix representing the design matrix of the model.

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 with effect terms used

call

The matched call of the model.

Description

Provides a high-level summary of a strainspy_fit object, including its call, coefficients, zero-inflation coefficients, dispersion, and the availability of FDR and residuals.

Usage

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

Arguments

Value

A list summarizing the strainspy_fit object.

Description

This function performs multiple testing adjustment on the results of a strainspy_fit object integrating taxonomic information. Various correction methods are available, including Harmonic Mean P-value (HMP), Bonferroni and Benjamini-Hochberg (BH).

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)

Description

Update the model fit by calling the specified fit function only a subset of features. Same design will be fitted.

Usage

update_fit(
  fit,
  update_idx,
  se,
  scale_continuous = TRUE,
  min_identity = 0.98,
  method = "glmmTMB",
  family = glmmTMB::ordbeta(),
  nthreads = 1,
  BPPARAM = NULL,
  progress = TRUE
)

Arguments

Value

Updated version of the fit object