HTSSIP (version 1.4.1)

BD_shift: Assessing the magnitude of BD shifts with 16S rRNA community data by calculating the beta diversity between unlabeled control and labeled treatment gradient fraction communities.

Description

This function is meant to compare 16S rRNA sequence communities of gradient fractions from 2 gradients: a labeled treatment (eg., 13C-labeled DNA) and its corresponding unlabeled control. First, the beta-diversity (e.g, weighted-Unifrac) is calculated pairwise between fraction communities.

Usage

BD_shift(physeq, method = "unifrac", weighted = TRUE, fast = TRUE,
  normalized = FALSE, ex = "Substrate=='12C-Con'",
  perm_method = c("control", "treatment", "overlap", "adjacent"),
  nperm = 100, a = 0.1, parallel_perm = FALSE,
  parallel_dist = FALSE)

Arguments

physeq

phyloseq object

method

See phyloseq::distance

weighted

Weighted Unifrac (if calculating Unifrac)

fast

Fast calculation method

normalized

Normalized abundances

ex

Expression for selecting controls based on metadata

perm_method

"BD shift window" permutation method. See description.

nperm

Number of bootstrap permutations

a

The alpha for calculating confidence intervals

parallel_perm

Calculate bootstrap permutations in parallel

parallel_dist

Calculate beta-diveristy distances in parallel

Value

a data.frame object of weighted mean distances

Details

The sample_data table of the user-provided phyloseq object MUST contain the buoyant density (BD) of each sample (a "Buoyant_density" column in the sample_data table). The BD information is used to identify overlapping gradient fractions (gradient fractions usually only partially overlap in BD between gradients) between the labeled treatment gradient and the control gradient. Beta diversity between overlapping fractions is calculated. Then, to standardize the values relative to the unlabeled control (1 beta-diversity value for each control gradient fraction), the mean beta diversity of overlapping labeled treatment gradients is calculated for each unlabeled control, and the percent overlap of each labeled treatment fraction is used to weight the mean.

A permutation test is used to determine "BD shift windows". OTU abundances are permuted, and beta-diversity is calculated. The permutations are used to calculate confidence intervals. The possible permutation methods are:

• "control" = OTU abundances are permuted among all control samples, and these new samples are used as a null treatment. Thus, this provides a baseline beta-diversity distribution that would result from comparing the control fractions to a randomly shuffled version of themselves.

• "treatment" = OTU abundances are permuted among all treatment samples. Thus, a "homogenized" treatment gradient null model.

• "overlap" = OTU abundances are permuted among overlapping control & treatment fractions. Thus, is beta-diversity higher than if the overlapping treatment & control samples were homogenized. This method tends to be too permisive.

• "adjacent" = The null "treatment" communities are generated by permuting OTU abundances among adjacent control fractions. Thus, null model is local gradient region was homogenized.

Examples

# NOT RUN {
data(physeq_S2D2)
# }
# NOT RUN {
# Subsetting phyloseq by Substrate and Day
params = get_treatment_params(physeq_S2D2, c('Substrate', 'Day'))
params = dplyr::filter(params, Substrate!='12C-Con')
ex = "(Substrate=='12C-Con' & Day=='${Day}') | (Substrate=='${Substrate}' & Day == '${Day}')"
physeq_S2D2_l = phyloseq_subset(physeq_S2D2, params, ex)

# Calculating BD_shift on 1 subset (use lapply function to process full list)
wmean1 = BD_shift(physeq_S2D2_l[[1]], nperm=5)

ggplot(wmean1, aes(BD_min.x, wmean_dist)) +
   geom_point()

# Calculating BD_shift on all subsets; using just 5 permutations to speed up analysis
lapply(physeq_S2D2_l, BD_shift, nperm=5)
# }
# NOT RUN {
# }