[R] 각 taxa에 대한 그룹간 비모수통계검정의 p-value값 정리 함수

작성: 2023-10-13-금

이전글

[R] Abundance의 평균, 표준편차, N, 95% CI(upper, lower)값 계산해주는 함수

[R] alpha diveristy의 각 index의 Mean, Sd, N, Upper, Lower, 95% CI값과 통계 분석 결과의 p-value값 저장하기

 

 

예제 데이터 

- QIIME2 tutorial인 moving-picture의 데이터입니다. 사람의 양 손바닥, 혀, 장의 마이크로바이옴 데이터를 담고 있습니다. 현재 저장된 형식은 phyloseq object입니다. phylsoeq형식은 biom format을 차용하여 R에서 편리하게 다루도록 개발된 것으로,  metadata와 otu table, taxonomy 그리고 계통정보를 한 번에 저장할 수 있습니다. 

ps.rds

0.08MB

 

 

함수 목표

- 각 taxa에서 body.site 에 따른 Abundance 비교하기

 

 

 

각 tax의 abundance비교 함수

library(dplyr)    # 데이터 처리 
library(phyloseq) # phyloseq
library(rlang)    # 함수 구조
library(combinat) # 순열 생산


ps <- readRDS("./ps.rds")
ps.rel <- phyloseq::transform_sample_counts(ps, fun = function(x) x/sum(x))

먼저 데이터를 불러오고, 상대적 분포도로 변형시켜줍니다.

 

Calculate_P_Values <- function(ps, tax_level, group_column, output_path) {
  # Phyloseq 데이터를 melt
  melt <- psmelt(ps)
  
  # 각 taxa 추출 
  unique_taxa <- unique(melt[[tax_level]])
  total_results <- data.frame()
  
  
  # taxon별로 계산하기  
  for (taxon in unique_taxa) {
    
    # Subset 데이터: 특정 Taxonomy level에 대한 데이터만 추출
    melt.2 <- melt %>% filter(!!sym(tax_level) == taxon)
    
    
    # Kruskal-Wallis 테스트를 수행하여 p-value 계산
    kruskal_results <- melt.2 %>%
      summarise(p_value = kruskal.test(as.formula(glue("Abundance ~", group_column)))$p.value)
    
    # 각 그룹 간 Wilcoxon 테스트를 수행하여 p-value 계산
    group_names <- unique(melt.2[[group_column]])
    statistic_results <- combn(group_names, 2, function(pair) { # 함수를 이용해서 순열 계산 (body.site가 짝지어서 for문처럼 돌아감. 그러나 for문보다 속도가 빠르다)
      subset_data <- melt.2 %>% filter(!!sym(group_column) %in% pair)
      p_value <- wilcox.test(subset_data$Abundance ~ subset_data[[group_column]])$p.value
      data.frame(Group1 = pair[1], Group2 = pair[2], wilcox_p_value = p_value)
    }, simplify = FALSE) %>%
      bind_rows()
    
    # p-value값 보정
    statistic_results[, "p_adjust"] <- p.adjust(statistic_results[, "wilcox_p_value"], method = "fdr")
    statistic_results[, "Kruskal-wallus_p_value"]  <- kruskal_results 
    statistic_results[, "taxonomy"] <- taxon
    
    total_results <- bind_rows(total_results, statistic_results)
  } 
  total_results
  
  # 결과 저장
  write.csv(total_results, paste0(output_path, "Compar_P_Values_", tax_level, ".csv"), row.names = F)
}

각 값을 계산한 후, 저장까지 완료하는 함수를 만들어 봤습니다. 

Phylum level에서는 시간이 얼마 걸리지 않았지만, Species 아래 레벨은 1분 이상의 시간이 걸릴 것으로 예상됩니다. 

 

결과 보기

- Phylum에서 body.site 에 따른 Abundance 비교해 보자

ps.rel.ph <- tax_glom(ps.rel, "Phylum")
Calculate_P_Values(ps.rel.ph, "Phylum", "body.site", "./")

Compar_P_Values_Phylum.csv

0.01MB

결과 미리보기

각 그룹에 따라서 wilcoxon rank sum와 kruskal wallis test의 p-value값이 표시된 것을 볼 수 있다. 물론 kruskal wallis와 taxonomy데이터가 중복되어 공간 낭비가 발생했다. R안에서는 이를 축약해서 볼 수 있지만, 엑셀에서 볼 때는 이렇게 데이터를 녹이는 것이 데이터를 찾기에 편할 것이라고 생각했다. 

 

읽어주셔서 감사합니다.