多组样本两两差异分析

使用 Wilcox test 对多组进行两两比较

phy为phyloseq对象，group为分组vector

DA_wilcox_paired <- function(phy, group){
  metadata <- sample_data(phy) %>% 
    as.matrix() %>%
    as.data.frame() %>%
    rownames_to_column('SampleID')
  comb_matrix <- combn(group,2)
  res_list <- list()
  for (i in 1:ncol(comb_matrix)){
    group1 <- comb_matrix[,i][1]
    group2 <- comb_matrix[,i][2]
    group1_samples <- metadata %>% filter(ExpGroup == group1) %>% .$SampleID
    group2_samples <- metadata %>% filter(ExpGroup == group2) %>% .$SampleID
    group_name <- paste0(group1, "_", group2)
    phy_tmp <- prune_samples(c(group1_samples, group2_samples), phy)
    phy_tmp <- prune_taxa(taxa_sums(phy_tmp) > 0, phy_tmp)
    metadata_tmp <- sample_data(phy_tmp) %>% 
      as.matrix() %>%
      as.data.frame() %>%
      rownames_to_column('SampleID')
    abundance_tab <- otu_table(phy_tmp) %>%
      as.data.frame() %>%
      rownames_to_column('OTUID') %>%
      pivot_longer(!OTUID, names_to = 'SampleID', values_to = 'value') %>%
      left_join(., metadata_tmp, by = 'SampleID')
    wilcox_res <- NULL
    for (j in unique(abundance_tab$OTUID)){
      input_data <- abundance_tab %>% filter(OTUID == j)
      input_group1 <- input_data %>% filter(ExpGroup == group1) %>% .$value
      input_group2 <- input_data %>% filter(ExpGroup == group2) %>% .$value
      input_group1_mean <- mean(input_group1)
      input_group2_mean <- mean(input_group2)
      enrihment_value <- ifelse(input_group1_mean > input_group2_mean, group1, group2)
      res <- wilcox.test(input_group1, input_group2)
      wilcox_res <- rbind(wilcox_res, c(group_name, j, enrihment_value, res$p.value))
    }
    wilcox_res <- as.data.frame(wilcox_res)
    colnames(wilcox_res) <- c('compare', 'OTUID', 'enrichment', 'pvalue')
    wilcox_res$p.adjust <- p.adjust(wilcox_res$pvalue, method = 'fdr')
    res_list[[group_name]] <- wilcox_res
  }
  return(res_list)
}

DA_paired_res <- DA_wilcox_paired(phy = phy_relab, group = c('A', 'B', 'C'))