一條函數使用三種方法檢測群落差異(MRPP,anosim,adonis)

微生物生態學領域經常要使用三種差異分析方法，分別似乎MRPP，anosim，adonis。很多人也問，到底使用哪種方法比較好，當然也有一些博文對他們有過簡單的區別，但是目前沒什麼大的不同。這三種方法是可以一起使用的。

不僅僅要一起使用，還要再有多個分組的情況下，進行兩兩比較。那這便是咱們這次的故事了：

pairMicroTest(ps = ps,dist = "bray",Micromet = "MRPP")

運行結果
        ID             stat        p
1 KO_VS_OE   MRPP.delta 0.3 p: 0.002
2 KO_VS_WT MRPP.delta 0.295 p: 0.002
3 OE_VS_WT MRPP.delta 0.275 p: 0.008
pairMicroTest = function(ps = ps,dist = "bray",Micromet = "anosim"){
  
  ps1_rela  = transform_sample_counts(ps, function(x) x / sum(x) );ps1_rela
  library(vegan)
  #-準備矩陣和分組文件
  map = as.data.frame(sample_data(ps1_rela))
  
  aa = levels(map$Group)
  aa
  aaa = combn(aa,2)
  aaa
  dim(aaa)[2]
  
  # 構建三個空列
  ID = rep("a",dim(aaa)[2])
  R = rep("a",dim(aaa)[2])
  P = rep("a",dim(aaa)[2])
  # i = 2
  for (i in 1:dim(aaa)[2]) {
    print(i)
    Desep_group = aaa[,i]
    map = as.data.frame(sample_data(ps1_rela))
    head(map)
    map$ID = row.names(map)
    maps<- dplyr::filter(map,Group %in% Desep_group)
    row.names(maps) = maps$ID
    ps_sub = ps1_rela
    sample_data( ps_sub ) = maps
    # ps_sub <- phyloseq::subset_samples(ps1_rela,Group %in% Desep_group);ps_sub
    
    
    ps_sub = phyloseq::filter_taxa(ps_sub, function(x) sum(x ) > 0 , TRUE);ps_sub
    map = as.data.frame(sample_data(ps_sub))
    unif <- phyloseq::distance(ps_sub  , method=dist, type="samples")
    print(ps_sub)
    
    if (Micromet == "MRPP") {
      mrpp = vegan::mrpp(unif, map$Group)
      as = round(mrpp$delta,3)
      R2 <- paste("MRPP.delta ",as, sep = "")
      R2
      p_v = paste("p: ",round(mrpp$Pvalue,3), sep = "")
      p_v

    }
    
    if (Micromet == "anosim") {
      dat.ano = anosim(unif, map$Group)
      a = round(dat.ano$statistic,3)
      R2 <- paste("ANOSIM.r ",a, sep = "")
      p_v = paste("p: ",round(dat.ano$signif,3), sep = "")

    }
    if (Micromet == "adonis") {
      ado =  vegan::adonis(unif  ~ map$Group,permutations = 999)
      a = round(as.data.frame(ado$aov.tab[5])[1,1],3)
      R2 <- paste("adonis:R ",a, sep = "")
      b = as.data.frame(ado$aov.tab[6])[1,1]
      p_v = paste("p: ",b, sep = "")
    }
    ID[i] = paste(Desep_group[1],Desep_group[2],sep = "_VS_")
    P[i] = p_v
    R[i] = R2
  }

  result = data.frame(ID = ID,stat = R,p = P)
  head(result)
  
  return(result)
}

我雖然使用了替代的方法可以繼續按想法實現這個功能，但是%in%的問題並沒有解決，不僅是這個函數，還有我曾今推送的maptre中的函數：
# 提取otu表格
  otu_table = as.data.frame(t(vegan_otu(ps1_rela)))
  otu_table$mean = rowMeans(otu_table)
  otu_table$ID = row.names(otu_table)
  head(otu_table)
  #按照從大到小排序
  otu_table<- arrange(otu_table, desc(mean))
  subtab = head(otu_table,N)
  head(subtab)
  row.names(subtab) =subtab$ID
  subtab = subtab[,1:(dim(subtab)[2]-2)]
  subtab = as.matrix(subtab)

  #對phyloseq取子集
  ps_sub <- phyloseq(otu_table(subtab, taxa_are_rows=TRUE),
                     tax_table(tax))
  ps_sub
本來極提取子集使用一下函數是最簡單的了，但是還是一樣的問題，這個函數再function中不起作用。兩個命令相似的特徵都是使用了%in%符號。
ps <- ps %>%
  subset_taxa(
    Kingdom %in% "Bacteria" 
  )
所以%in%符號和function之間有什麼隱情呢? 希望可以得到大家的幫助。