In HuaZou/MicrobiomeAnalysis: Data analysis toolkits in metagenomics

library(knitr)
knit_hooks$set(crop = hook_pdfcrop)
knitr::opts_chunk$set(crop=TRUE, tidy=FALSE, warning=FALSE, message=FALSE, fig.align="center")
Biocpkg <- function (pkg){
    sprintf("[%s](http://bioconductor.org/packages/%s)", pkg, pkg)
}
CRANpkg <- function(pkg){
    cran <- "https://CRAN.R-project.org/package" 
    fmt <- "[%s](%s=%s)"
    sprintf(fmt, pkg, cran, pkg) 
}

library(ggplot2)
library(phyloseq)
library(dplyr)
library(tibble)
library(MicrobiomeAnalysis)
library(conflicted)
conflict_prefer("filter", "dplyr")
conflict_prefer("select", "dplyr")

Overview of MicrobiomeAnalysis

knitr::include_graphics("../man/figures/Schematic.png")

Convert inputs into phyloseq-class object

inputs from dada2

seq_tab <- readRDS(
  system.file("extdata", "dada2_seqtab.rds",
     package = "MicrobiomeAnalysis"))

tax_tab <- readRDS(
  system.file("extdata", "dada2_taxtab.rds",
    package = "MicrobiomeAnalysis"))

sam_tab <- read.table(
  system.file("extdata", "dada2_samdata.txt",
     package = "MicrobiomeAnalysis"), 
  sep = "\t", header = TRUE, row.names = 1)

ps <- import_dada2(seq_tab = seq_tab, 
                   tax_tab = tax_tab, 
                   sam_tab = sam_tab)

ps

inputs from qiime2

otuqza_file <- system.file(
     "extdata", "table.qza",
     package = "MicrobiomeAnalysis")

taxaqza_file <- system.file(
     "extdata", "taxonomy.qza",
     package = "MicrobiomeAnalysis")

sample_file <- system.file(
     "extdata", "sample-metadata.tsv",
     package = "MicrobiomeAnalysis")

treeqza_file <- system.file(
     "extdata", "tree.qza",
     package = "MicrobiomeAnalysis")

ps_object <- import_qiime2(
     otu_qza = otuqza_file, 
     taxa_qza = taxaqza_file,
     sam_tab = sample_file, 
     tree_qza = treeqza_file)
ps_object

Data Processing

Extract specific levels phyloseq-class object

• Extracting "Genus" levels phyloseq object

ps_genus <- aggregate_taxa(ps_object, 
                           level = "Genus")
ps_genus

head(ps_genus@otu_table@.Data, 3)

Summarize the specific taxonomic levels

• Phyloseq object contains from Kingdom to the the specific taxonomic levels (Phylum)

ps_summarize_genus <- summarize_taxa(
    ps_object, 
    level = "Genus")
ps_summarize_genus

Data Transformation

• Transform the taxa abundances in otu_table, which applying a mathematical transformation on individual values themselves.

ps_genus_transform <- transform_abundances(
    object = ps_genus,
    transform = "log10p")

head(ps_genus_transform@otu_table@.Data, 3)

Data Imputation

• Imputing missing value of OTU_table in phyloseq-class object

ps_genus_impute <- impute_abundance(
    object = ps_genus,
    group = "body.site",
    ZerosAsNA = TRUE,
    RemoveNA = TRUE,
    cutoff = 20,    
    method = "half_min")

head(ps_genus_impute@otu_table@.Data, 3)

Data Normalization

• Normalizing the OTU_table in phyloseq-class object sample by sample to reduce the effects of systematic differences such as library size.

ps_genus_norm <- normalize(
    object = ps_genus,
    method = "TSS")

head(ps_genus_norm@otu_table@.Data, 3)

Data Scaling

• Data scaling adjusts each variable/feature by a scaling factor computed based on the dispersion of the variable.

ps_genus_scale <- scale_variables(
    object = ps_genus,
    method = "zscore")

head(ps_genus_scale@otu_table@.Data, 3)

Data Trimming

• Trimming samples or taxa whose prevalence is less than threshold

ps_genus_trim <- trim_prevalence(
    object = ps_genus,
    level = NULL,
    cutoff = 0.1,
    trim = "feature")

ps_genus_trim

head(ps_genus_trim@otu_table@.Data, 3)

Data Filtering

• Filtering taxa who is low relative abundance or unclassified (Ref: [@thingholm2019obese])

• Taxa more than Mean relative abundance across all samples: 100;

• Taxa more than Minimum relative abundance at least one sample: 1000.

ps_genus_filter <- filter_abundance(
    object = ps_genus,
    level = NULL,
    cutoff_mean = 100,
    cutoff_one = 1000,
    unclass = TRUE)

ps_genus_filter

head(ps_genus_filter@otu_table@.Data, 3)

Multivariate Community Analysis

Know more details of the aftermentioned statistical methods to see [@xia2018statistical].

Multivariate homogeneity of groups dispersions (variances)

[@anderson2006multivariate]

run_betadisper(object = ps_object,
               level = "Genus",
               variable = "body.site", 
               method = "bray")

Permutational Multivariate Analysis of Variance (PERMANOVA)

[@anderson2014permutational]

run_PERMANOVA(object = ps_object,
              level = NULL,
              method = "bray")

Mantel Test (MANTEL)

[@mantel1967detection]

run_MANTEL(object = ps_object,
           y_variables = "body.site",
           z_variables = c("subject", "reported.antibiotic.usage"),
           norm = FALSE,
           method = "mantel.partial",
           method_dist = c("bray", "euclidean", "jaccard"))

Analysis of Similarity (ANOSIM)

[@clarke1993non]

run_ANOSIM(object = ps_object,
           level = "Genus",
           variable = "body.site", 
           method = "bray")

Multi-response Permutation Procedures (MRPP)

[@mielke1991application]

run_MRPP(object = ps_object,
         level = "Genus",
         variable = "body.site", 
         method = "bray")

Alpha diversity

• calculating alpha diversity index on microbiota data

alphaindex <- get_alphaindex(
  ps = ps_object,
  level = "Genus",
  indices = c("Shannon", "Chao1"))

head(alphaindex)

Beta dispersion

• Permutation test for homogeneity of multivariate dispersions on microbiota data

data("Zeybel_2022_gut")

run_betadisper(
  object = Zeybel_2022_gut,
  level = "Phylum",
  variable = "LiverFatClass",
  method = "bray")

Beta diversity

• Principal Coordinate Analysis

data("Zeybel_2022_gut")
ps_zeybel <- summarize_taxa(Zeybel_2022_gut, level = "Genus")

ord_result <- run_ord(
  object = ps_zeybel,
  variable = "LiverFatClass",
  method = "PCoA")

plot_ord(
  reslist = ord_result,
  variable = "LiverFatClass",
  ellipse_type = "ellipse_groups",
  sideboxplot = TRUE,
  sample_label = FALSE)

Differential analysis

Aldex

ps_genus_group <- phyloseq::subset_samples(
     ps_genus,
     body.site %in% c("gut", "tongue"))

run_aldex(ps = ps_genus_group, 
          group = "body.site", 
          taxa_rank = "Genus")

limma-voom

run_limma_voom(ps = ps_genus_group, 
               group = "body.site", 
               taxa_rank = "Phylum")

ANCOM

data("enterotypes_arumugam")
ps_enterotypes <- phyloseq::subset_samples(
     enterotypes_arumugam,
     Enterotype %in% c("Enterotype 3", "Enterotype 2"))

run_ancom(ps = ps_enterotypes, 
          group = "Enterotype")

ANCOMBC

run_ancombc(ps = ps_enterotypes, 
            group = "Enterotype", 
            confounders = "Gender")

DESeq2

run_deseq2(ps = ps_enterotypes, 
           group = "Enterotype")

edgeR

run_edger(ps_enterotypes, group = "Enterotype")

lefse

run_lefse(ps = ps_enterotypes, 
          group = "Enterotype")

metagenomeSeq

run_metagenomeseq(ps = ps_enterotypes, 
                  group = "Enterotype")

test_multiple_groups

• anova

• kruskal

ps <- phyloseq::subset_samples(
    enterotypes_arumugam,
    Enterotype %in% c("Enterotype 3", "Enterotype 2", "Enterotype 1"))

run_test_multiple_groups(ps = ps, 
                         group = "Enterotype", 
                         method = "anova")

supervised leaning (SL) methods

• logistic regression

ps_small <- phyloseq::subset_taxa(
    enterotypes_arumugam,
    Phylum %in% c("Firmicutes", "Bacteroidetes")
)

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

mm

Differential analysis in metabolomics

• Differential approaches:

• FoldChange analysis

• VIP (Variable important in projection) analysis

• T-test

data("Zeybel_2022_protein")
Zeybel_2022_protein_imput <- impute_abundance(
      Zeybel_2022_protein,
      group = "LiverFatClass",
      method = "knn")
Zeybel_2022_protein_norm <- scale_variables(
      Zeybel_2022_protein_imput,
      method == "zscore")
DA_results <- run_metabolomeDA(
  object_raw = Zeybel_2022_protein,
  object_norm = Zeybel_2022_protein_norm,
  variable = "LiverFatClass",
  variable_name = c("None", "Severe"))

head(DA_results[, 1:4], 4)

• Volcano to display the significant results

plot_volcano(
   da_res = DA_results,
   group_names = c("None", "Severe"),
   x_index = "Log2FoldChange",
   x_index_cutoff = 0.1,
   y_index = "AdjustedPvalue",
   y_index_cutoff = 0.5,
   group_color = c("red", "grey", "blue"),
   add_enrich_arrow = TRUE,
   topN = 10)

Session information {-}

This vignette was created under the following conditions:

sessionInfo()

References {-}

HuaZou/MicrobiomeAnalysis documentation built on May 13, 2024, 11:10 a.m.