작성: 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 그리고 계통정보를 한 번에 저장할 수 있습니다.
함수 목표
- 각 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", "./")
각 그룹에 따라서 wilcoxon rank sum와 kruskal wallis test의 p-value값이 표시된 것을 볼 수 있다. 물론 kruskal wallis와 taxonomy데이터가 중복되어 공간 낭비가 발생했다. R안에서는 이를 축약해서 볼 수 있지만, 엑셀에서 볼 때는 이렇게 데이터를 녹이는 것이 데이터를 찾기에 편할 것이라고 생각했다.
읽어주셔서 감사합니다.
작성: 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 그리고 계통정보를 한 번에 저장할 수 있습니다.
함수 목표
- 각 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", "./")각 그룹에 따라서 wilcoxon rank sum와 kruskal wallis test의 p-value값이 표시된 것을 볼 수 있다. 물론 kruskal wallis와 taxonomy데이터가 중복되어 공간 낭비가 발생했다. R안에서는 이를 축약해서 볼 수 있지만, 엑셀에서 볼 때는 이렇게 데이터를 녹이는 것이 데이터를 찾기에 편할 것이라고 생각했다.
읽어주셔서 감사합니다.
