GMSB文章七：微生物整合分析

library(readr)
library(openxlsx)
library(compositions)
library(tidyverse) 
library(mia)
library(ggpubr)

df_v1 <- read_csv("./data/GMSB-data/df_v1.csv", show_col_types = FALSE)
bias_corr_species <- read_csv("./data/GMSB-data/results/outputs/bias_corr_species.csv")
raw_species <- readRDS("./data/GMSB-data/rds/primary_ancombc2_species.rds")
sig_species <- read.xlsx("./data/GMSB-data/results/outputs/res_ancombc2.xlsx", sheet = 1) 

df_v1 <- df_v1 %>%
  dplyr::filter(group1 != "missing")
​
# Microbiome data
bias_corr_species <- bias_corr_species %>%
  dplyr::rowwise() %>%
  dplyr::filter(grepl("Species:", species)|grepl("Genus:", species)) %>%
  dplyr::mutate(species = ifelse(grepl("Genus:", species), 
                        paste(strsplit(species, ":")[[1]][2], "spp."),
                        strsplit(species, ":")[[1]][2])) %>%
  dplyr::ungroup() 
​
# Significant taxa from ANCOM-BC2
sig_species <- sig_species %>%
  dplyr::filter(p_val < 0.05) %>%
  .$taxon
​
# Subset significant taxa
df_da_species <- bias_corr_species %>%
  dplyr::filter(species %in% sig_species)
df_da_species <- t(df_da_species)
colnames(df_da_species) <- df_da_species[1, ]
df_da_species <- data.frame(df_da_species[-1, , drop = FALSE], check.names = FALSE) %>%
  rownames_to_column("sampleid") %>%
  dplyr::mutate(across(-1, as.numeric))
df_da_species[is.na(df_da_species)] <- 0
​
# 合并分组标签和物种丰度表
df_add_species <- df_v1 %>%
  dplyr::left_join(df_da_species, by = "sampleid")
df_add_species <- data.frame(df_add_species)
​
sig_species <- make.names(sig_species)
​
head(sig_species)

[1] "A.muciniphila"     "A.onderdonkii"     "Anaerovibrio.spp." "B.adolescentis"    "B.caccae"         
[6] "B.fragilis"

lm_eqn <- function(m){
  
  a <- unname(coef(m))[1]
  b <- unname(coef(m))[2]
  p_val <- summary(m)$coef[2, "Pr(>|t|)"]
  
  b <- ifelse(sign(b) >= 0, 
             paste0(" + ", format(b, digits = 2)), 
             paste0(" - ", format(-b, digits = 2)))
  
  eq <- substitute(paste(italic(y) == a, b, italic(x), ", ", italic(p) == p_val),
                  list(a = format(a, digits = 2), b = b,
                       p_val = format(p_val, digits = 2)))
  
  return(as.character(as.expression(eq)))              
}
​
plot_scatter <- function(df, var_name, tax_name, y_lab) {
  
  df_fig <- df %>%
    dplyr::select(all_of(c(var_name, tax_name))) %>%
    dplyr::rename(y = var_name) %>%
    tidyr::pivot_longer(-y, names_to = "tax", values_to = "x")
  
  df_lm <- df_fig %>% 
    dplyr::group_by(tax) %>%
    dplyr::do(fit = lm(formula = y ~ x, data = .)) %>%
    dplyr::summarise(eq = lm_eqn(m = fit))
  
  df_ann <- df_fig %>% 
    dplyr::group_by(tax) %>%
    dplyr::summarise(y = ifelse(mean(y, na.rm = TRUE) > 0, 
                         0.5 * max(y, na.rm = TRUE),
                         0.2 * abs(mean(y, na.rm = TRUE))),
              x = median(x, na.rm = TRUE)) %>%
      dplyr::mutate(eq = df_lm$eq,
             y_max = 1.05 * y)
  
  fig <- df_fig %>% 
    ggplot(aes(x = x, y = y)) + 
    geom_point(alpha = 0.8, color = "#BEAED4") +
    geom_smooth(method = "lm", se = TRUE, color = "skyblue", 
                formula = y ~ x) +
    facet_wrap(.~tax, scales = "free") +
    labs(x = "Micorbial abundances", 
         y = y_lab, 
         title = NULL) + 
    geom_blank(data = df_ann, aes(y = y_max)) +
    geom_text(data = df_ann, aes(x = x, y = y, label = eq), size = 3, 
              parse = TRUE, inherit.aes = FALSE) + 
    theme_bw() +
    theme(strip.background = element_rect(fill = "white"),
          legend.position = "bottom",
          plot.title = element_text(hjust = 0.5))
  
  return(fig)
}

t_formula <- as.formula(paste0("cbind(crp, cd14, cd163) ~ ",
                              paste0(sig_species, collapse = " + ")))
fit <- manova(t_formula, data = df_add_species)
​
summ <- summary(fit, test = "Pillai")
​
df_summ <- data.frame(summ$stats) %>%
  rownames_to_column("Taxon") %>%
  dplyr::transmute(Taxon,
                   approx.F = signif(approx.F, 2),
                   num.Df, den.Df, 
                   P = round(`Pr..F.`, 2)) %>%
  dplyr::filter(Taxon != "Residuals") %>%
  dplyr::arrange(P)
​
head(df_summ)
​
cowplot::plot_grid(
  plot_scatter(
    df = df_add_species,
    var_name = "crp",
    tax_name = "Lachnospira.spp.",
    y_lab = "crp"),
  plot_scatter(
    df = df_add_species,
    var_name = "cd14",
    tax_name = "Lachnospira.spp.",
    y_lab = "cd14"),
  plot_scatter(
    df = df_add_species,
    var_name = "cd163",
    tax_name = "Lachnospira.spp.",
    y_lab = "cd163"),
  ncol = 2, align = "hv")

t_formula <- as.formula(paste0("cbind(acetate, valerate) ~ ",
                              paste0(sig_species, collapse = " + ")))
fit <- manova(t_formula, data = df_add_species)
​
summ <- summary(fit, test = "Pillai")
​
df_summ <- data.frame(summ$stats) %>%
  rownames_to_column("Taxon") %>%
  dplyr::transmute(Taxon,
                   approx.F = signif(approx.F, 2),
                   num.Df, den.Df, 
                   P = round(`Pr..F.`, 2)) %>%
  dplyr::filter(Taxon != "Residuals") %>%
  dplyr::arrange(P)
​
head(df_summ)
​
cowplot::plot_grid(
  plot_scatter(
    df = df_add_species,
    var_name = "acetate",
    tax_name = "B.uniformis",
    y_lab = "acetate"),
  plot_scatter(
    df = df_add_species,
    var_name = "valerate",
    tax_name = "B.uniformis",
    y_lab = "valerate"),
  ncol = 2, align = "hv")

fit <- manova(cbind(crp, cd14, cd163) ~ 
               acetate + valerate, 
             data = df_v1)
​
summ <- summary(fit, test = "Pillai")
​
df_summ <- data.frame(summ$stats) %>%
  rownames_to_column("Taxon") %>%
  dplyr::transmute(Taxon,
                   approx.F = signif(approx.F, 2),
                   num.Df, den.Df, 
                   P = round(`Pr..F.`, 2)) %>%
  dplyr::filter(Taxon != "Residuals") %>%
  dplyr::arrange(P)
​
head(df_summ)
​
cowplot::plot_grid(
  plot_scatter(
    df = df_add_species,
    var_name = "crp",
    tax_name = "acetate",
    y_lab = "crp"),
  plot_scatter(
    df = df_add_species,
    var_name = "cd14",
    tax_name = "acetate",
    y_lab = "cd14"),
  ncol = 2, align = "hv")