LorMe (version 1.1.0)

tbRDA_analysis: tbRDA analysis

Description

RDA analysis including co-linearity diagnostics and necessary statistics.

Usage

tbRDA_analysis(otudata, envdata, collinearity, perm.test = TRUE)

Value

Three permutation test result print ,one preview plot ,a RDA object(default name:otu.tab.1) and a summary of RDA object

Arguments

otudata

Feature table of all numeric variable, with annotation in row names

envdata

Environmental factor of all numeric variable,with sample-ID in row names and environmental factor in column names

collinearity

If done collinearity diagnostics. Default,TRUE.

perm.test

Logical. If conduct permutation test. Default:TRUE.

Author

Wang Ningqi

Examples

  ### Data preparation ###
  library(vegan)
  data(varechem)
  head(varechem)
  data(testotu)
  require(tidyr); require(magrittr)  ## Or use pipe command in "dplyr"

  sep_testotu <- Filter_function(
    input = testotu,
    threshold = 0.0001,
    format = 1
  ) %>%
  separate(
    ., col = taxonomy,
    into = c("Domain", "Phylum", "Order", "Family", "Class", "Genus", "Species"),
    sep = ";"
  )

  top10phylum <- aggregate(
    sep_testotu[, 2:21],
    by = list(sep_testotu$Phylum),
    FUN = sum
  ) %>%
  Top_taxa(
    input = .,
    n = 10,
    inputformat = 2,
    outformat = 1
  )
  rownames(top10phylum) <- top10phylum[, 1]
  top10phylum <- top10phylum[, -1]

  group <- data.frame(
    group = c(rep("a", 10), rep("b", 10)),
    factor1 = rnorm(10),
    factor2 = rnorm(mean = 100, 10)
  )

  ### RDA analysis ###
  set.seed(999)
  RDAresult <- tbRDA_analysis(
    top10phylum,
    varechem[1:20, ],
    TRUE
  )

  # Environmental statistics
  print(RDAresult$factor_statistics)

  # Visualization using ggplot
  rda_object <- RDAresult$rda_object
  rda_summary <- RDAresult$rdasummary
  rda_env <- as.data.frame(rda_summary$biplot)
  rda_sample <- as.data.frame(rda_summary$sites)
  rda_otu <- as.data.frame(rda_summary$species)

  xlab <- paste0("RDA1:", round(RDAresult$rdasummary$concont$importance[2, 1], 4) * 100, "%")
  ylab <- paste0("RDA2:", round(RDAresult$rdasummary$concont$importance[2, 2], 4) * 100, "%")

  library(ggrepel)
  # Create a sample RDA plot
  RDAplot <- ggplot(data = rda_sample, aes(RDA1, RDA2)) +
    geom_point(aes(color = group$group), size = 2) +
    geom_point(data = rda_otu, pch = "+", color = "orange", size = 4) +
    geom_hline(yintercept = 0) +
    geom_vline(xintercept = 0) +
    geom_segment(data = rda_env, aes(x = 0, y = 0, xend = RDA1 * 0.8, yend = RDA2 * 0.8),
                 arrow = arrow(angle = 22.5, length = unit(0.35, "cm")),
                 linetype = 1, size = 0.6, colour = "red") +
    geom_text_repel(color = "red", data = rda_env,
                    aes(RDA1, RDA2, label = row.names(rda_env))) +
    labs(x = xlab, y = ylab, color = "Treatment",
         title = paste0("p = ", anova.cca(rda_object)["Model", "Pr(>F)"])) +
    stat_ellipse(aes(color = group$group), level = 0.95) +
    geom_text_repel(size = 3, color = "orange",
                    data = subset(rda_otu, RDA1 > 0.1 | RDA1 < (-0.1)),
                    aes(RDA1, RDA2, label = rownames(subset(rda_otu, RDA1>0.1|RDA1<(-0.1))))) +
    theme_zg()

  # Print the RDA plot
  print(RDAplot)