Package 'microbiomeMarker'

Description

Operators acting on marker_table to extract parts.

Usage

## S4 method for signature 'marker_table,ANY,ANY,ANY'
x[i, j, ..., drop = TRUE]

Arguments

Value

a marker_table object.

See Also

base::Extract()

Description

Extract taxa abundances from phyloseq objects.

Usage

abundances(object, transform = c("identity", "log10", "log10p"), norm = FALSE)

## S4 method for signature 'otu_table'
abundances(object, transform = c("identity", "log10", "log10p"), norm = FALSE)

## S4 method for signature 'phyloseq'
abundances(object, transform = c("identity", "log10", "log10p"), norm = FALSE)

## S4 method for signature 'microbiomeMarker'
abundances(object, transform = c("identity", "log10", "log10p"))

Arguments

Value

abundance matrix with taxa in rows and samples in columns.

See Also

otu_table, phyloseq, microbiomeMarker,transform_abundances

Examples

data(caporaso)
abd <- abundances(caporaso)

Description

Summarize phyloseq data into a higher phylogenetic level.

Usage

aggregate_taxa(x, level, verbose = FALSE)

Arguments

Details

This provides a convenient way to aggregate phyloseq OTUs (or other taxa) when the phylogenetic tree is missing. Calculates the sum of OTU abundances over all OTUs that map to the same higher-level group. Removes ambiguous levels from the taxonomy table. Returns a phyloseq object with the summarized abundances.

Value

Summarized phyloseq object

Author(s)

Contact: Leo Lahti [email protected]

References

See citation('microbiome')

Examples

data(caporaso)
caporaso_phylum <- aggregate_taxa(caporaso, "Phylum")

Description

Assign a new OTU table in microbiomeMarker object

Usage

## S4 replacement method for signature 'microbiomeMarker,otu_table'
otu_table(x) <- value

## S4 replacement method for signature 'microbiomeMarker,phyloseq'
otu_table(x) <- value

## S4 replacement method for signature 'microbiomeMarker,microbiomeMarker'
otu_table(x) <- value

Arguments

Value

a microbiomeMarker object.

Description

Calculating power, false discovery rates, false positive rates and auc ( area under the receiver operating characteristic (ROC) curve) for various DA methods.

Usage

compare_DA(
  ps,
  group,
  taxa_rank = "none",
  methods,
  args = list(),
  n_rep = 20,
  effect_size = 5,
  k = NULL,
  relative = TRUE,
  BPPARAM = BiocParallel::SnowParam(progressbar = TRUE)
)

Arguments

Details

To make this function support for different arguments for a certain DA method args allows list of list of list e.g. args = list(lefse = list(list(norm = "CPM"), list(norm = "TSS"))), which specify to compare the different norm arguments for lefse analysis.

For taxa_rank, only taxa_rank = "none" is supported, if this argument is not "none", it will be forced to "none" internally.

Value

an compareDA object, which contains a two-length list of:

• metrics: data.frame, FPR, AUC and spike detection rate for each run.

• mm: differential analysis results.

Description

Confounding variables may mask the actual differential features. This function utilizes constrained correspondence analysis (CCA) to measure the confounding factors.

Usage

confounder(
  ps,
  target_var,
  norm = "none",
  confounders = NULL,
  permutations = 999,
  ...
)

Arguments

Value

a data.frame contains three variables: confounder, pseudo-F and p value.

Examples

data(caporaso)
confounder(caporaso, "SampleType", confounders = "ReportedAntibioticUsage")

Description

16S read counts and phylogenetic tree file of 34 Illumina samples derived from Moving Pictures of the Human Microbiome (Caporaso et al.) Group label: gut, left palm, right palm, and tongue - indicating different sampled body sites.

Format

a phyloseq::phyloseq object

Author(s)

Yang Cao

Source

Data was downloaded from https://www.microbiomeanalyst.ca

References

Caporaso, et al. Moving pictures of the human microbiome. Genome Biol 12, R50 (2011).

https://doi.org/10.1186/gb-2011-12-5-r50

Description

Data from a cohort of 94 Bone Marrow Transplant patients previously published on in CID

Format

a phyloseq::phyloseq object

Author(s)

Yang Cao

Source

https://github.com/ying14/yingtools2/tree/master/data

References

Ying, et al. Intestinal Domination and the Risk of Bacteremia in Patients Undergoing Allogeneic Hematopoietic Stem Cell Transplantation, Clinical Infectious Diseases, Volume 55, Issue 7, 1 October 2012, Pages 905–914,

https://academic.oup.com/cid/article/55/7/905/428203

Description

The data from a subset of the Early Childhood Antibiotics and the Microbiome (ECAM) study, which tracked the microbiome composition and development of 43 infants in the United States from birth to 2 years of age, identifying microbiome associations with antibiotic exposure, delivery mode, and diet.

Format

a phyloseq::phyloseq object

References

Bokulich, Nicholas A., et al. "Antibiotics, birth mode, and diet shape microbiome maturation during early life." Science translational medicine 8.343 (2016): 343ra82-343ra82.

https://github.com/FrederickHuangLin/ANCOM/tree/master/data

Description

The data contains 22 European metagenomes from Danish, French, Italian, and Spanish individuals, and 13 Japanese and 4 American.

Format

a phyloseq::phyloseq object

Author(s)

Yang Cao

References

Arumugam, Manimozhiyan, et al. Enterotypes of the human gut microbiome. nature 473.7346 (2011): 174-180.

Description

The data from a study on colorectal cancer. Samples that had no DIAGNOSIS attribute assigned and with less than 500 reads (counts) were removed, and 191 samples remains (91 healthy and 86 Tumors).

Format

a phyloseq::phyloseq object

Author(s)

Yang Cao

References

Kostic et al. Genomic analysis identifies association of Fusobacterium with colorectal carcinoma. Genome research, 2012, 22(2), 292-298.

Description

A small subset of the HMP 16S dataset for finding biomarkers characterizing different level of oxygen availability in different bodysites

Format

a phyloseq::phyloseq object

Author(s)

Yang Cao

Source

http://huttenhower.sph.harvard.edu/webfm_send/129

Description

43 pediatric IBD stool samples obtained from the Integrative Human Microbiome Project Consortium (iHMP). Group label: CD and Controls.

Format

a phyloseq::phyloseq object

Author(s)

Yang Cao

Source

https://www.microbiomeanalyst.ca/MicrobiomeAnalyst/resources

Description

The dataset contains 30 abundance profiles (obtained processing the 16S reads with RDP) belonging to 10 rag2 (control) and 20 truc (case) mice.

Format

a phyloseq::phyloseq object

Author(s)

Yang Cao

Source

http://www.huttenhower.org/webfm_send/73

Description

This function extracts the results of posthoc test.

Usage

extract_posthoc_res(object, features = NULL)

Arguments

Value

a IRanges::SimpleDFrameList object.

Examples

require(IRanges)
pht <- postHocTest(
    result = DataFrameList(
        featureA = DataFrame(
            comparisons = c("group2-group1", 
                "group3-group1", 
                "group3-group2"),
            diff_mean = runif(3),
            pvalue = rep(0.01, 3),
            ci_lower = rep(0.01, 3),
            ci_upper = rep(0.011, 3)
        ),
        featureB = DataFrame(
            comparisons = c("group2-group1", 
                "group3-group1", 
                "group3-group2"),
            diff_mean = runif(3),
            pvalue = rep(0.01, 3),
            ci_lower = rep(0.01, 3),
            ci_upper = rep(0.011, 3)
        )
    ),
    abundance = data.frame(
        featureA = runif(3),
        featureB = runif(3),
        group = c("group1", "group2", "grou3")
    )
)
extract_posthoc_res(pht, "featureA")[[1]]

Description

Import the output of dada2 into phyloseq object

Usage

import_dada2(
  seq_tab,
  tax_tab = NULL,
  sam_tab = NULL,
  phy_tree = NULL,
  keep_taxa_rows = TRUE
)

Arguments

Details

The output of the dada2 pipeline is a feature table of amplicon sequence variants (an ASV table): A matrix with rows corresponding to samples and columns to ASVs, in which the value of each entry is the number of times that ASV was observed in that sample. This table is analogous to the traditional OTU table. Conveniently, taxa names are saved as ASV1, ASV2, ..., in the returned phyloseq object.

Value

phyloseq::phyloseq object hold the taxonomy info, sample metadata, number of reads per ASV.

Examples

seq_tab <- readRDS(system.file("extdata", "dada2_seqtab.rds",
    package = "microbiomeMarker"
))
tax_tab <- readRDS(system.file("extdata", "dada2_taxtab.rds",
    package = "microbiomeMarker"
))
sam_tab <- read.table(system.file("extdata", "dada2_samdata.txt",
    package = "microbiomeMarker"
), sep = "\t", header = TRUE, row.names = 1)
ps <- import_dada2(seq_tab = seq_tab, tax_tab = tax_tab, sam_tab = sam_tab)
ps

Description

Import the output of picrust2 into phyloseq object

Usage

import_picrust2(
  feature_tab,
  sam_tab = NULL,
  trait = c("PATHWAY", "COG", "EC", "KO", "PFAM", "TIGRFAM", "PHENO")
)

Arguments

Details

PICRUSt2 is a software for predicting abundances of functional profiles based on marker gene sequencing data. The functional profiles can be predicted from the taxonomic profiles using PICRUSt2. "Function" usually refers to gene families such as KEGG orthologs and Enzyme Classification numbers, but predictions can be made for any arbitrary trait.

In the ⁠phyloseq object⁠, the predicted function abundance profile is stored in otu_table slot. And the functional trait is saved in tax_table slot, if the descriptions of function features is not added to the predicted table, tax_table will have only one rank Picrust_trait to represent the function feature id, or if the desciptions are added one more rank Picrust_description will be added to represent the description of function feature.

Value

phyloseq::phyloseq object.

Examples

sam_tab <- system.file(
    "extdata", "picrust2_metadata.tsv",
    package = "microbiomeMarker")
feature_tab <- system.file(
    "extdata", "path_abun_unstrat_descrip.tsv.gz",
    package = "microbiomeMarker")
ps <- import_picrust2(feature_tab, sam_tab, trait = "PATHWAY")
ps

Description

Import the qiime2 artifacts, including feature table, taxonomic table, phylogenetic tree, representative sequence and sample metadata into phyloseq object.

Usage

import_qiime2(
  otu_qza,
  taxa_qza = NULL,
  sam_tab = NULL,
  refseq_qza = NULL,
  tree_qza = NULL
)

Arguments

Value

phyloseq::phyloseq object.

Examples

otuqza_file <- system.file(
    "extdata", "table.qza",
    package = "microbiomeMarker"
)
taxaqza_file <- system.file(
    "extdata", "taxonomy.qza",
    package = "microbiomeMarker"
)
sample_file <- system.file(
    "extdata", "sample-metadata.tsv",
    package = "microbiomeMarker"
)
treeqza_file <- system.file(
    "extdata", "tree.qza",
    package = "microbiomeMarker"
)
ps <- import_qiime2(
    otu_qza = otuqza_file, taxa_qza = taxaqza_file,
    sam_tab = sample_file, tree_qza = treeqza_file
)
ps

Description

This is the recommended function for both building and accessing microbiome marker table (marker_table).

Usage

marker_table(object)

## S4 method for signature 'data.frame'
marker_table(object)

## S4 method for signature 'microbiomeMarker'
marker_table(object)

Arguments

Value

a marker_table object.

Examples

data(enterotypes_arumugam)
mm <- run_limma_voom(
    enterotypes_arumugam,
    "Enterotype",
    contrast = c("Enterotype 3", "Enterotype 2"),
    pvalue_cutoff = 0.05,
    p_adjust = "fdr"
)
marker_table(mm)

Description

This Class is inherit from data.frame. Rows represent the microbiome markers and variables represents feature of the marker.

Fields

names,row.names

a character vector, inherited from the input data.frame

.data

a list, each element corresponding the each column of the input data.frame

.S3Class

character, the S3 class marker_table inherited from: "data.frame"

Author(s)

Yang Cao

Description

This function replace the marker_table slot of object with value.

Usage

marker_table(object) <- value

Arguments

Value

a microbiomeMarker object.

Examples

data(enterotypes_arumugam)
mm <- run_limma_voom(
    enterotypes_arumugam,
    "Enterotype",
    contrast = c("Enterotype 3", "Enterotype 2"),
    pvalue_cutoff = 0.1,
    p_adjust = "fdr"
)
mm_marker <- marker_table(mm)
mm_marker
marker_table(mm) <- mm_marker[1:2, ]
marker_table(mm)

Description

This the constructor to build the microbiomeMarker object, don't use the new() constructor.

Usage

microbiomeMarker(
  marker_table = NULL,
  norm_method = NULL,
  diff_method = NULL,
  ...
)

Arguments

Value

a microbiomeMarker object.

See Also

phyloseq::phyloseq()

Examples

microbiomeMarker(
    marker_table = marker_table(data.frame(
        feature = c("speciesA", "speciesB"),
        enrich_group = c("groupA", "groupB"),
        ef_logFC = c(-2, 2),
        pvalue = c(0.01, 0.01),
        padj = c(0.01, 0.01),
        row.names = c("marker1", "marker2")
    )),
    norm_method = "TSS",
    diff_method = "DESeq2",
    otu_table = otu_table(matrix(
        c(4, 1, 1, 4),
        nrow = 2, byrow = TRUE,
        dimnames = list(c("speciesA", "speciesB"), c("sample1", "sample2"))
    ),
    taxa_are_rows = TRUE
    ),
    tax_table = tax_table(matrix(
        c("speciesA", "speciesB"),
        nrow = 2,
        dimnames = list(c("speciesA", "speciesB"), "Species")
    )),
    sam_data = sample_data(data.frame(
        group = c("groupA", "groupB"),
        row.names = c("sample1", "sample2")
    ))
)

Description

microbiomeMarker-class is inherited from the phyloseq::phyloseq by adding a custom slot microbiome_marker to save the differential analysis results. And it provides a seamless interface with phyloseq, which makes microbiomeMarker simple and easy to use. For more details on see the document of phyloseq::phyloseq.

Usage

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

Arguments

Value

a microbiomeMarker object.

Slots

marker_table

a data.frame, a marker_table object.

norm_method

character, method used to normalize the input phyloseq object.

diff_method

character, method used for marker identification.

See Also

phyloseq::phyloseq, marker_table, summarize_taxa()

Description

Get the number of microbiome markers

Usage

nmarker(object)

## S4 method for signature 'microbiomeMarker'
nmarker(object)

## S4 method for signature 'marker_table'
nmarker(object)

Arguments

Value

an integer, the number of microbiome markers

Examples

mt <- marker_table(data.frame(
    feature = c("speciesA", "speciesB"),
    enrich_group = c("groupA", "groupB"),
    ef_logFC = c(-2, 2),
    pvalue = c(0.01, 0.01),
    padj = c(0.01, 0.01),
    row.names = c("marker1", "marker2")
))
nmarker(mt)

Description

It is critical to normalize the feature table to eliminate any bias due to differences in the sampling sequencing depth.This function implements six widely-used normalization methods for microbial compositional data.

For rarefying, reads in the different samples are randomly removed until the same predefined number has been reached, to assure all samples have the same library size. Rarefying normalization method is the standard in microbial ecology. Please note that the authors of phyloseq do not advocate using this rarefying a normalization procedure, despite its recent popularity

TSS simply transforms the feature table into relative abundance by dividing the number of total reads of each sample.

CSS is based on the assumption that the count distributions in each sample are equivalent for low abundant genes up to a certain threshold. Only the segment of each sample’s count distribution that is relatively invariant across samples is scaled by CSS

RLE assumes most features are not differential and uses the relative abundances to calculate the normalization factor.

TMM calculates the normalization factor using a robust statistics based on the assumption that most features are not differential and should, in average, be equal between the samples. The TMM scaling factor is calculated as the weighted mean of log-ratios between each pair of samples, after excluding the highest count OTUs and OTUs with the largest log-fold change.

In CLR, the log-ratios are computed relative to the geometric mean of all features.

norm_cpm: This normalization method is from the original LEfSe algorithm, recommended when very low values are present (as shown in the LEfSe galaxy).

Usage

## S4 method for signature 'phyloseq'
normalize(object, method = "TSS", ...)

## S4 method for signature 'otu_table'
normalize(object, method = "TSS", ...)

## S4 method for signature 'data.frame'
normalize(object, method = "TSS", ...)

## S4 method for signature 'matrix'
normalize(object, method = "TSS", ...)

norm_rarefy(
  object,
  size = min(sample_sums(object)),
  rng_seed = FALSE,
  replace = TRUE,
  trim_otus = TRUE,
  verbose = TRUE
)

norm_tss(object)

norm_css(object, sl = 1000)

norm_rle(
  object,
  locfunc = stats::median,
  type = c("poscounts", "ratio"),
  geo_means = NULL,
  control_genes = NULL
)

norm_tmm(
  object,
  ref_column = NULL,
  logratio_trim = 0.3,
  sum_trim = 0.05,
  do_weighting = TRUE,
  Acutoff = -1e+10
)

norm_clr(object)

norm_cpm(object)

Arguments

Value

the same class with object.

See Also

edgeR::calcNormFactors(),DESeq2::estimateSizeFactorsForMatrix(), metagenomeSeq::cumNorm()

phyloseq::rarefy_even_depth()

metagenomeSeq::calcNormFactors()

DESeq2::estimateSizeFactorsForMatrix()

edgeR::calcNormFactors()

Examples

data(caporaso)
normalize(caporaso, "TSS")

Description

This function convert [phyloseq::phyloseq-class⁠] to [⁠DESeq2::DESeqDataSet-class⁠], which can then be tested using [⁠DESeq2::DESeq()'].

Usage

phyloseq2DESeq2(ps, design, ...)

Arguments

Value

a DESeq2::DESeqDataSet object.

See Also

DESeq2::DESeqDataSetFromMatrix(),DESeq2::DESeq()

Examples

data(caporaso)
phyloseq2DESeq2(caporaso, ~SampleType)

Description

This function convert phyloseq::phyloseq object to edgeR::DGEList object, can then can be used to perform differential analysis using the methods in edgeR.

Usage

phyloseq2edgeR(ps, ...)

Arguments

Value

A edgeR::DGEList object.

Examples

data(caporaso)
dge <- phyloseq2edgeR(caporaso)

Description

The phyloseq data is converted to the relevant metagenomeSeq::MRexperiment object, which can then be tested in the zero-inflated mixture model framework in the metagenomeSeq package.

Usage

phyloseq2metagenomeSeq(ps, ...)

otu_table2metagenomeSeq(ps, ...)

Arguments

Value

A metagenomeSeq::MRexperiment object.

See Also

metagenomeSeq::fitTimeSeries(), metagenomeSeq::fitLogNormal(),metagenomeSeq::fitZig(), metagenomeSeq::MRtable(),metagenomeSeq::MRfulltable()

Examples

data(caporaso)
phyloseq2metagenomeSeq(caporaso)

Description

plot the abundances of markers

Usage

plot_abundance(mm, label_level = 1, max_label_len = 60, markers = NULL, group)

Arguments

Value

a ggplot2::ggplot object.

Examples

data(enterotypes_arumugam)
mm <- run_limma_voom(
    enterotypes_arumugam,
    "Enterotype",
    contrast = c("Enterotype 3", "Enterotype 2"),
    pvalue_cutoff = 0.01,
    p_adjust = "none"
)
plot_abundance(mm, group = "Enterotype")

Description

plot cladogram of micobiomeMaker results

Usage

plot_cladogram(
  mm,
  color,
  only_marker = FALSE,
  branch_size = 0.2,
  alpha = 0.2,
  node_size_scale = 1,
  node_size_offset = 1,
  clade_label_level = 4,
  clade_label_font_size = 4,
  annotation_shape = 22,
  annotation_shape_size = 5,
  group_legend_param = list(),
  marker_legend_param = list()
)

Arguments

Value

a ggtree object

Author(s)

Chenhao Li, Guangchuang Yu, Chenghao Zhu, Yang Cao

References

This function is modified from clada.anno from microbiomeViz.

See Also

ggtree::ggtree()

Examples

data(kostic_crc)
kostic_crc_small <- phyloseq::subset_taxa(
    kostic_crc,
    Phylum %in% c("Firmicutes")
)
mm_lefse <- run_lefse(
    kostic_crc_small,
    wilcoxon_cutoff = 0.01,
    group = "DIAGNOSIS",
    kw_cutoff = 0.01,
    multigrp_strat = TRUE,
    lda_cutoff = 4
)
plot_cladogram(mm_lefse, color = c("darkgreen", "red"))

Description

bar and dot plot of effect size microbiomeMarker data. This function returns a ggplot2 object that can be saved or further customized using ggplot2 package.

Usage

plot_ef_bar(mm, label_level = 1, max_label_len = 60, markers = NULL)

plot_ef_dot(mm, label_level = 1, max_label_len = 60, markers = NULL)

Arguments

Value

a ggplot project

Examples

data(enterotypes_arumugam)
mm <- run_limma_voom(
    enterotypes_arumugam,
    "Enterotype",
    contrast = c("Enterotype 3", "Enterotype 2"),
    pvalue_cutoff = 0.01,
    p_adjust = "none"
)
plot_ef_bar(mm)

Description

Display the microbiome marker using heatmap, in which rows represents the marker and columns represents the samples.

Usage

plot_heatmap(
  mm,
  transform = c("log10", "log10p", "identity"),
  cluster_marker = FALSE,
  cluster_sample = FALSE,
  markers = NULL,
  label_level = 1,
  max_label_len = 60,
  sample_label = FALSE,
  scale_by_row = FALSE,
  annotation_col = NULL,
  group,
  ...
)

Arguments

Value

a ComplexHeatmap::Heatmap object.

See Also

transform_abundances,ComplexHeatmap::Heatmap()

Examples

data(kostic_crc)
kostic_crc_small <- phyloseq::subset_taxa(
    kostic_crc,
    Phylum %in% c("Firmicutes")
)
mm_lefse <- run_lefse(
    kostic_crc_small,
    wilcoxon_cutoff = 0.01,
    group = "DIAGNOSIS",
    kw_cutoff = 0.01,
    multigrp_strat = TRUE,
    lda_cutoff = 4
)
plot_heatmap(mm_lefse, group = "DIAGNOSIS")

Description

Visualize the result of post-hoc test using ggplot2

Usage

plot_postHocTest(pht, feature, step_increase = 0.12)

Arguments

Value

a ggplot object

Examples

data(enterotypes_arumugam)
ps <- phyloseq::subset_samples(
    enterotypes_arumugam,
    Enterotype %in% c("Enterotype 3", "Enterotype 2", "Enterotype 1")
) %>%
    phyloseq::subset_taxa(Phylum == "Bacteroidetes")
pht <- run_posthoc_test(ps, group = "Enterotype")
plot_postHocTest(pht, feature = "p__Bacteroidetes|g__Alistipes")

Description

Show the ROC curve of the microbiome marker calculated by run_sl.

Usage

plot_sl_roc(mm, group, nfolds = 3, nrepeats = 3, tune_length = 5, ...)

Arguments

Value

a ggplot2::ggplot object.

See Also

run_sl()

Examples

data(enterotypes_arumugam)
# small example phyloseq object for test
ps_s <- phyloseq::subset_taxa(
    enterotypes_arumugam,
    Phylum %in% c("Firmicutes", "Bacteroidetes")
)

set.seed(2021)
mm <- run_sl(
    ps_s,
    group = "Gender",
    taxa_rank = "Genus",
    nfolds = 2,
    nrepeats = 1,
    top_n = 15,
    norm = "TSS",
    method = "LR",
)
plot_sl_roc(mm, group = "Gender")

Description

Plotting DA comparing result

Usage

## S3 method for class 'compareDA'
plot(x, sort = c("score", "auc", "fpr", "power"), ...)

Arguments

Value

a ggplot2::ggplot object containing 4 subplots: "auc", "fdr", "power"and "score" plot.

Description

This function is used for create postHocTest object, and is only used for developers.

Usage

postHocTest(
  result,
  abundance,
  conf_level = 0.95,
  method = "tukey",
  method_str = paste("Posthoc multiple comparisons of means: ", method)
)

Arguments

Value

a postHocTest object.

Examples

require(IRanges)
pht <- postHocTest(
    result = DataFrameList(
        featureA = DataFrame(
            comparisons = c("group2-group1", 
                "group3-group1", 
                "group3-group2"),
            diff_mean = runif(3),
            pvalue = rep(0.01, 3),
            ci_lower = rep(0.01, 3),
            ci_upper = rep(0.011, 3)
        ),
        featureB = DataFrame(
            comparisons = c("group2-group1", 
                "group3-group1", 
                "group3-group2"),
            diff_mean = runif(3),
            pvalue = rep(0.01, 3),
            ci_lower = rep(0.01, 3),
            ci_upper = rep(0.011, 3)
        )
    ),
    abundance = data.frame(
        featureA = runif(3),
        featureB = runif(3),
        group = c("group1", "group2", "grou3")
    )
)
pht

Description

The postHocTest Class, represents the result of post-hoc test result among multiple groups

Usage

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

Arguments

Value

a postHocTest object.

Slots

result

a IRanges::DataFrameList, each DataFrame consists of five variables:

• comparisons: character, specify which two groups to test (the group names are separated by "_)

• diff_mean: numeric, difference in mean abundances

• pvalue: numeric, p values

• ci_lower and ci_upper: numeric, lower and upper confidence interval of difference in mean abundances

abundance

abundance of each feature in each group

conf_level

confidence level

method

method used for post-hoc test

method_str

method illustration

Author(s)

Yang Cao

Description

Perform differential analysis using ALDEx2

Usage

run_aldex(
  ps,
  group,
  taxa_rank = "all",
  transform = c("identity", "log10", "log10p"),
  norm = "none",
  norm_para = list(),
  method = c("t.test", "wilcox.test", "kruskal", "glm_anova"),
  p_adjust = c("none", "fdr", "bonferroni", "holm", "hochberg", "hommel", "BH", "BY"),
  pvalue_cutoff = 0.05,
  mc_samples = 128,
  denom = c("all", "iqlr", "zero", "lvha"),
  paired = FALSE
)

Arguments

Value

a microbiomeMarker object.

References

Fernandes, A.D., Reid, J.N., Macklaim, J.M. et al. Unifying the analysis of high-throughput sequencing datasets: characterizing RNA-seq, 16S rRNA gene sequencing and selective growth experiments by compositional data analysis. Microbiome 2, 15 (2014).

See Also

ALDEx2::aldex()

Examples

data(enterotypes_arumugam)
ps <- phyloseq::subset_samples(
    enterotypes_arumugam,
    Enterotype %in% c("Enterotype 3", "Enterotype 2")
)
run_aldex(ps, group = "Enterotype")

Description

Perform significant test by comparing the pairwise log ratios between all features.

Usage

run_ancom(
  ps,
  group,
  confounders = character(0),
  taxa_rank = "all",
  transform = c("identity", "log10", "log10p"),
  norm = "TSS",
  norm_para = list(),
  p_adjust = c("none", "fdr", "bonferroni", "holm", "hochberg", "hommel", "BH", "BY"),
  pvalue_cutoff = 0.05,
  W_cutoff = 0.75
)

Arguments

Details

In an experiment with only two treatments, this tests the following hypothesis for feature $i$:

$H_{0i}: E(log(\mu_i^1)) = E(log(\mu_i^2))$

where $\mu_i^1$ and $\mu_i^2$ are the mean abundances for feature $i$ in the two groups.

The developers of this method recommend the following significance tests if there are 2 groups, use non-parametric Wilcoxon rank sum test stats::wilcox.test(). If there are more than 2 groups, use nonparametric stats::kruskal.test() or one-way ANOVA stats::aov().

Value

a microbiomeMarker object, in which the slot of marker_table contains four variables:

• feature, significantly different features.

• enrich_group, the class of the differential features enriched.

• effect_size, differential means for two groups, or F statistic for more than two groups.

• W, the W-statistic, number of features that a single feature is tested to be significantly different against.

Author(s)

Huang Lin, Yang Cao

References

Mandal et al. "Analysis of composition of microbiomes: a novel method for studying microbial composition", Microbial Ecology in Health & Disease, (2015), 26.

Examples

data(enterotypes_arumugam)
ps <- phyloseq::subset_samples(
    enterotypes_arumugam,
    Enterotype %in% c("Enterotype 3", "Enterotype 2")
)
run_ancom(ps, group = "Enterotype")

Description

Differential abundance analysis for microbial absolute abundance data. This function is a wrapper of ANCOMBC::ancombc().

Usage

run_ancombc(
  ps,
  group,
  confounders = character(0),
  contrast = NULL,
  taxa_rank = "all",
  transform = c("identity", "log10", "log10p"),
  norm = "none",
  norm_para = list(),
  p_adjust = c("none", "fdr", "bonferroni", "holm", "hochberg", "hommel", "BH", "BY"),
  prv_cut = 0.1,
  lib_cut = 0,
  struc_zero = FALSE,
  neg_lb = FALSE,
  tol = 1e-05,
  max_iter = 100,
  conserve = FALSE,
  pvalue_cutoff = 0.05
)

Arguments

Details

contrast must be a two length character or NULL (default). It is only required to set manually for two groups comparison when there are multiple groups. The order determines the direction of comparison, the first element is used to specify the reference group (control). This means that, the first element is the denominator for the fold change, and the second element is used as baseline (numerator for fold change). Otherwise, users do required to concern this parameter (set as default NULL), and if there are two groups, the first level of groups will set as the reference group; if there are multiple groups, it will perform an ANOVA-like testing to find markers which difference in any of the groups.

Value

a microbiomeMarker object.

References

Lin, Huang, and Shyamal Das Peddada. "Analysis of compositions of microbiomes with bias correction." Nature communications 11.1 (2020): 1-11.

See Also

ANCOMBC::ancombc

Examples

data(enterotypes_arumugam)
ps <- phyloseq::subset_samples(
    enterotypes_arumugam,
    Enterotype %in% c("Enterotype 3", "Enterotype 2")
)
run_ancombc(ps, group = "Enterotype")

Description

Differential expression analysis based on the Negative Binomial distribution using DESeq2.

Usage

run_deseq2(
  ps,
  group,
  confounders = character(0),
  contrast = NULL,
  taxa_rank = "all",
  norm = "RLE",
  norm_para = list(),
  transform = c("identity", "log10", "log10p"),
  fitType = c("parametric", "local", "mean", "glmGamPoi"),
  sfType = "poscounts",
  betaPrior = FALSE,
  modelMatrixType,
  useT = FALSE,
  minmu = ifelse(fitType == "glmGamPoi", 1e-06, 0.5),
  p_adjust = c("none", "fdr", "bonferroni", "holm", "hochberg", "hommel", "BH", "BY"),
  pvalue_cutoff = 0.05,
  ...
)

Arguments

Details

Note: DESeq2 requires the input is raw counts (un-normalized counts), as only the counts values allow assessing the measurement precision correctly. For more details see the vignette of DESeq2 (vignette("DESeq2")).

Thus, this function only supports "none", "rarefy", "RLE", "CSS", and "TMM" normalization methods. We strongly recommend using the "RLE" method (default normalization method in the DESeq2 package). The other normalization methods are used for expert users and comparisons among different normalization methods.

For two groups comparison, this function utilizes the Wald test (defined by DESeq2::nbinomWaldTest()) for hypothesis testing. A Wald test statistic is computed along with a probability (p-value) that a test statistic at least as extreme as the observed value were selected at random. contrasts are used to specify which two groups to compare. The order of the names determines the direction of fold change that is reported.

Likelihood ratio test (LRT) is used to identify the genes that significantly changed across all the different levels for multiple groups comparisons. The LRT identified the significant features by comparing the full model to the reduced model. It is testing whether a feature removed in the reduced model explains a significant variation in the data.

contrast must be a two length character or NULL (default). It is only required to set manually for two groups comparison when there are multiple groups. The order determines the direction of comparison, the first element is used to specify the reference group (control). This means that, the first element is the denominator for the fold change, and the second element is used as baseline (numerator for fold change). Otherwise, users do required to concern this parameter (set as default NULL), and if there are two groups, the first level of groups will set as the reference group; if there are multiple groups, it will perform an ANOVA-like testing to find markers which difference in any of the groups.

Value

a microbiomeMarker object.

References

Love, Michael I., Wolfgang Huber, and Simon Anders. "Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2." Genome biology 15.12 (2014): 1-21.

See Also

DESeq2::results(),DESeq2::DESeq()

Examples

data(enterotypes_arumugam)
ps <- phyloseq::subset_samples(
    enterotypes_arumugam,
    Enterotype %in% c("Enterotype 3", "Enterotype 2")) %>%
    phyloseq::subset_taxa(Phylum %in% c("Firmicutes"))
run_deseq2(ps, group = "Enterotype")

Description

Differential expression analysis based on the Negative Binomial distribution using edgeR.

Usage

run_edger(
  ps,
  group,
  confounders = character(0),
  contrast = NULL,
  taxa_rank = "all",
  method = c("LRT", "QLFT"),
  transform = c("identity", "log10", "log10p"),
  norm = "TMM",
  norm_para = list(),
  disp_para = list(),
  p_adjust = c("none", "fdr", "bonferroni", "holm", "hochberg", "hommel", "BH", "BY"),
  pvalue_cutoff = 0.05,
  ...
)

Arguments

Details

Note that edgeR is designed to work with actual counts. This means that transformation is not required in any way before inputting them to edgeR.

There are two test methods for differential analysis in edgeR, likelihood ratio test (LRT) and quasi-likelihood F-test (QLFT). The QLFT method is recommended as it allows stricter error rate control by accounting for the uncertainty in dispersion estimation.

contrast must be a two length character or NULL (default). It is only required to set manually for two groups comparison when there are multiple groups. The order determines the direction of comparison, the first element is used to specify the reference group (control). This means that, the first element is the denominator for the fold change, and the second element is used as baseline (numerator for fold change). Otherwise, users do required to concern this parameter (set as default NULL), and if there are two groups, the first level of groups will set as the reference group; if there are multiple groups, it will perform an ANOVA-like testing to find markers which difference in any of the groups.

Value

a microbiomeMarker object.

Author(s)

Yang Cao

References

Robinson, Mark D., and Alicia Oshlack. "A scaling normalization method for differential expression analysis of RNA-seq data." Genome biology 11.3 (2010): 1-9.

Robinson, Mark D., Davis J. McCarthy, and Gordon K. Smyth. "edgeR: a Bioconductor package for differential expression analysis of digital gene expression data." Bioinformatics 26.1 (2010): 139-140.

See Also

edgeR::glmFit(),edgeR::glmQLFit(),edgeR::estimateDisp() ,normalize()

Examples

data(enterotypes_arumugam)
ps <- phyloseq::subset_samples(
    enterotypes_arumugam,
    Enterotype %in% c("Enterotype 3", "Enterotype 2")
)
run_edger(ps, group = "Enterotype")

Description

Perform Metagenomic LEFSe analysis based on phyloseq object.

Usage

run_lefse(
  ps,
  group,
  subgroup = NULL,
  taxa_rank = "all",
  transform = c("identity", "log10", "log10p"),
  norm = "CPM",
  norm_para = list(),
  kw_cutoff = 0.05,
  lda_cutoff = 2,
  bootstrap_n = 30,
  bootstrap_fraction = 2/3,
  wilcoxon_cutoff = 0.05,
  multigrp_strat = FALSE,
  strict = c("0", "1", "2"),
  sample_min = 10,
  only_same_subgrp = FALSE,
  curv = FALSE
)

Arguments

Value

a microbiomeMarker object, in which the slot of marker_table contains four variables:

• feature, significantly different features.

• enrich_group, the class of the differential features enriched.

• lda, logarithmic LDA score (effect size)

• pvalue, p value of kw test.

Author(s)

Yang Cao

References

Segata, Nicola, et al. Metagenomic biomarker discovery and explanation. Genome biology 12.6 (2011): R60.

See Also

normalize

Examples

data(kostic_crc)
kostic_crc_small <- phyloseq::subset_taxa(
    kostic_crc,
    Phylum == "Firmicutes"
)
mm_lefse <- run_lefse(
    kostic_crc_small,
    wilcoxon_cutoff = 0.01,
    group = "DIAGNOSIS",
    kw_cutoff = 0.01,
    multigrp_strat = TRUE,
    lda_cutoff = 4
)

Description

Differential analysis using limma-voom

Usage

run_limma_voom(
  ps,
  group,
  confounders = character(0),
  contrast = NULL,
  taxa_rank = "all",
  transform = c("identity", "log10", "log10p"),
  norm = "none",
  norm_para = list(),
  voom_span = 0.5,
  p_adjust = c("none", "fdr", "bonferroni", "holm", "hochberg", "hommel", "BH", "BY"),
  pvalue_cutoff = 0.05,
  ...
)

Arguments

Details

contrast must be a two length character or NULL (default). It is only required to set manually for two groups comparison when there are multiple groups. The order determines the direction of comparison, the first element is used to specify the reference group (control). This means that, the first element is the denominator for the fold change, and the second element is used as baseline (numerator for fold change). Otherwise, users do required to concern this parameter (set as default NULL), and if there are two groups, the first level of groups will set as the reference group; if there are multiple groups, it will perform an ANOVA-like testing to find markers which difference in any of the groups.

Value

a microbiomeMarker object.

References

Law, C. W., Chen, Y., Shi, W., & Smyth, G. K. (2014). voom: Precision weights unlock linear model analysis tools for RNA-seq read counts. Genome biology, 15(2), 1-17.

Examples

data(enterotypes_arumugam)
mm <- run_limma_voom(
    enterotypes_arumugam,
    "Enterotype",
    contrast = c("Enterotype 3", "Enterotype 2"),
    pvalue_cutoff = 0.01,
    p_adjust = "none"
)
mm

Description

run_marker is a wrapper of all differential analysis functions.

Usage

run_marker(
  ps,
  group,
  da_method = c("lefse", "simple_t", "simple_welch", "simple_white", "simple_kruskal",
    "simple_anova", "edger", "deseq2", "metagenomeseq", "ancom", "ancombc", "aldex",
    "limma_voom", "sl_lr", "sl_rf", "sl_svm"),
  taxa_rank = "all",
  transform = c("identity", "log10", "log10p"),
  norm = "none",
  norm_para = list(),
  p_adjust = c("none", "fdr", "bonferroni", "holm", "hochberg", "hommel", "BH", "BY"),
  pvalue_cutoff = 0.05,
  ...
)

Arguments

Details

This function is only a wrapper of all differential analysis functions, We recommend to use the corresponding function, since it has a better default arguments setting.

Value

a microbiomeMarker object.

See Also

run_lefse(),run_simple_stat(),run_test_two_groups(), run_test_multiple_groups(),run_edger(),run_deseq2(), run_metagenomeseq,run_ancom(),run_ancombc(),run_aldex(), run_limma_voom(),run_sl()

Description

Differential expression analysis based on the Zero-inflated Log-Normal mixture model or Zero-inflated Gaussian mixture model using metagenomeSeq.

Usage

run_metagenomeseq(
  ps,
  group,
  confounders = character(0),
  contrast = NULL,
  taxa_rank = "all",
  transform = c("identity", "log10", "log10p"),
  norm = "CSS",
  norm_para = list(),
  method = c("ZILN", "ZIG"),
  p_adjust = c("none", "fdr", "bonferroni", "holm", "hochberg", "hommel", "BH", "BY"),
  pvalue_cutoff = 0.05,
  ...
)

Arguments

Details

metagnomeSeq provides two differential analysis methods, zero-inflated log-normal mixture model (implemented in metagenomeSeq::fitFeatureModel()) and zero-inflated Gaussian mixture model (implemented in metagenomeSeq::fitZig()). We recommend fitFeatureModel over fitZig due to high sensitivity and low FDR. Both metagenomeSeq::fitFeatureModel() and metagenomeSeq::fitZig() require the abundance profiles before normalization.

For metagenomeSeq::fitZig(), the output column is the coefficient of interest, and logFC column in the output of metagenomeSeq::fitFeatureModel() is analogous to coefficient. Thus, logFC is really just the estimate the coefficient of interest in metagenomeSeq::fitFeatureModel(). For more details see these question Difference between fitFeatureModel and fitZIG in metagenomeSeq.

contrast must be a two length character or NULL (default). It is only required to set manually for two groups comparison when there are multiple groups. The order determines the direction of comparison, the first element is used to specify the reference group (control). This means that, the first element is the denominator for the fold change, and the second element is used as baseline (numerator for fold change). Otherwise, users do required to concern this paramerter (set as default NULL), and if there are two groups, the first level of groups will set as the reference group; if there are multiple groups, it will perform an ANOVA-like testing to find markers which difference in any of the groups.

Of note, metagenomeSeq::fitFeatureModel() is not allows for multiple groups comparison.

Value

a microbiomeMarker object.

Author(s)

Yang Cao

References

Paulson, Joseph N., et al. "Differential abundance analysis for microbial marker-gene surveys." Nature methods 10.12 (2013): 1200-1202.

Examples

data(enterotypes_arumugam)
ps <- phyloseq::subset_samples(
    enterotypes_arumugam,
    Enterotype %in% c("Enterotype 3", "Enterotype 2")
)
run_metagenomeseq(ps, group = "Enterotype")

Description

Multiple group test, such as anova and Kruskal-Wallis rank sum test, can be used to uncover the significant feature among all groups. Post hoc tests are used to uncover specific mean differences between pair of groups.

Usage

run_posthoc_test(
  ps,
  group,
  transform = c("identity", "log10", "log10p"),
  norm = "TSS",
  norm_para = list(),
  conf_level = 0.95,
  method = c("tukey", "games_howell", "scheffe", "welch_uncorrected")
)

Arguments

Value

a postHocTest object

See Also

postHocTest, run_test_multiple_groups()

Examples

data(enterotypes_arumugam)
ps <- phyloseq::subset_samples(
    enterotypes_arumugam,
    Enterotype %in% c("Enterotype 3", "Enterotype 2", "Enterotype 1")
) %>%
    phyloseq::subset_taxa(Phylum == "Bacteroidetes")
pht <- run_posthoc_test(ps, group = "Enterotype")
pht

Description

Perform simple statistical analysis of metagenomic profiles. This function is a wrapper of run_test_two_groups and run_test_multiple_groups.

Usage

run_simple_stat(
  ps,
  group,
  taxa_rank = "all",
  transform = c("identity", "log10", "log10p"),
  norm = "TSS",
  norm_para = list(),
  method = c("welch.test", "t.test", "white.test", "anova", "kruskal"),
  p_adjust = c("none", "fdr", "bonferroni", "holm", "hochberg", "hommel", "BH", "BY"),
  pvalue_cutoff = 0.05,
  diff_mean_cutoff = NULL,
  ratio_cutoff = NULL,
  eta_squared_cutoff = NULL,
  conf_level = 0.95,
  nperm = 1000,
  ...
)

Arguments

Value

a microbiomeMarker object.

See Also

run_test_two_groups(),run_test_multiple_groups()

Examples

data(enterotypes_arumugam)
ps <- phyloseq::subset_samples(
    enterotypes_arumugam,
    Enterotype %in% c("Enterotype 3", "Enterotype 2")
)
run_simple_stat(ps, group = "Enterotype")

Description

Identify biomarkers using logistic regression, random forest, or support vector machine.

Usage

run_sl(
  ps,
  group,
  taxa_rank = "all",
  transform = c("identity", "log10", "log10p"),
  norm = "none",
  norm_para = list(),
  nfolds = 3,
  nrepeats = 3,
  sampling = NULL,
  tune_length = 5,
  top_n = 10,
  method = c("LR", "RF", "SVM"),
  ...
)

Arguments

Details

Only support two groups comparison in the current version. And the marker was selected based on its importance score. Moreover, The hyper-parameters are selected automatically by a grid-search based method in the N-time K-fold cross-validation. Thus, the identified biomarker based can be biased due to model overfitting for small datasets (e.g., with less than 100 samples).

The argument top_n is used to denote the number of markers based on the importance score. There is no rule or principle on how to select top_n, however, usually it is very useful to try a different top_n and compare the performance of the marker predictions for the testing data.

Value

a microbiomeMarker object.

Author(s)

Yang Cao

See Also

caret::train(),caret::trainControl()

Examples

data(enterotypes_arumugam)
# small example phyloseq object for test
ps_small <- phyloseq::subset_taxa(
    enterotypes_arumugam,
    Phylum %in% c("Firmicutes", "Bacteroidetes")
)

set.seed(2021)
mm <- run_sl(
    ps_small,
    group = "Gender",
    taxa_rank = "Genus",
    nfolds = 2,
    nrepeats = 1,
    top_n = 15,
    norm = "TSS",
    method = "LR",
)
mm

Description

Statistical test for multiple groups

Usage

run_test_multiple_groups(
  ps,
  group,
  taxa_rank = "all",
  transform = c("identity", "log10", "log10p"),
  norm = "TSS",
  norm_para = list(),
  method = c("anova", "kruskal"),
  p_adjust = c("none", "fdr", "bonferroni", "holm", "hochberg", "hommel", "BH", "BY"),
  pvalue_cutoff = 0.05,
  effect_size_cutoff = NULL
)

Arguments

Value

a microbiomeMarker object.

See Also

run_posthoc_test(),run_test_two_groups(),run_simple_stat()

Examples

data(enterotypes_arumugam)
ps <- phyloseq::subset_samples(
    enterotypes_arumugam,
    Enterotype %in% c("Enterotype 3", "Enterotype 2", "Enterotype 1")
)
mm_anova <- run_test_multiple_groups(
    ps,
    group = "Enterotype",
    method = "anova"
)

Description

Statistical test between two groups

Usage

run_test_two_groups(
  ps,
  group,
  taxa_rank = "all",
  transform = c("identity", "log10", "log10p"),
  norm = "TSS",
  norm_para = list(),
  method = c("welch.test", "t.test", "white.test"),
  p_adjust = c("none", "fdr", "bonferroni", "holm", "hochberg", "hommel", "BH", "BY"),
  pvalue_cutoff = 0.05,
  diff_mean_cutoff = NULL,
  ratio_cutoff = NULL,
  conf_level = 0.95,
  nperm = 1000,
  ...
)

Arguments

Value

a microbiomeMarker object.

Author(s)

Yang Cao

See Also

run_test_multiple_groups(),run_simple_stat

Examples

data(enterotypes_arumugam)
mm_welch <- run_test_two_groups(
    enterotypes_arumugam,
    group = "Gender",
    method = "welch.test"
)
mm_welch

Description

Subset markers based on an expression related to the columns and values within the marker_table slot of mm.

Usage

subset_marker(mm, ...)

Arguments

Value

a subset object in the same class with mm.

Examples

data(enterotypes_arumugam)
mm <- run_limma_voom(
    enterotypes_arumugam,
    "Enterotype",
    contrast = c("Enterotype 3", "Enterotype 2"),
    pvalue_cutoff = 0.01,
    p_adjust = "none"
)
subset_marker(mm, pvalue < 0.005)

Description

Provides summary information of the representation of a taxonomic levels within each sample.

Usage

summarize_taxa(ps, level = rank_names(ps)[1], absolute = TRUE, sep = "|")

Arguments

Value

a phyloseq::phyloseq object, where each row represents a taxa, and each col represents the taxa abundance of each sample.

Examples

data(enterotypes_arumugam)
summarize_taxa(enterotypes_arumugam)

Description

Summary differential analysis methods comparison results

Usage

## S3 method for class 'compareDA'
summary(
  object,
  sort = c("score", "auc", "fpr", "power"),
  boot = TRUE,
  boot_n = 1000L,
  prob = c(0.05, 0.95),
  ...
)

Arguments

Value

a data.frame containing measurements for differential analysis methods:

• call: differential analysis commands.

• auc: area under curve of ROC.

• fpr: false positive rate

• power: empirical power.

• fdr: false discover7y rate.

• score: score whch is calculated as $(auc - 0.5) * power - fdr$.

• ⁠score_*⁠: confidence limits of score.

Description

Transform the taxa abundances in otu_table sample by sample, which means the counts of each sample will be transformed individually.

Usage

transform_abundances(object, transform = c("identity", "log10", "log10p"))

Arguments

Value

A object matches the class of argument object with the transformed otu_table.

See Also

abundances()

Examples

data(oxygen)
x1 <- transform_abundances(oxygen)
head(otu_table(x1), 10)
x2 <- transform_abundances(oxygen, "log10")
head(otu_table(x2), 10)
x3 <- transform_abundances(oxygen, "log10p")
head(otu_table(x3), 10)