작성 : 2023-06-07
수정 : 2024-01-26
PICRUSt2 결과를 어떻게 보여줘야 보는 이로 하여금 이해가 수월할지 고민해 보았다.
마이크로바이옴 논문에서는 Erro bar를 사용하여 각 비교 그룹에서 얼마나 양적으로 차이 나는지 보여주거나, heatmap을 사용하는 것이 대다수였다. 혹은 분석 결과를 LDA score 등 를 통해서 보여주기도 한다.
전체 샘플에서 kegg pathway에 해당하는 분포를 보고 싶을 때에는 heatmap 이 가장 적합하다. heatmap은 다른 시각화 방법보다 raw data를 반영한다.
이를 고려하여, 전반적인 raw data를 보여주기 위해서 heatmap을 사용하고, 통계적으로 유의한 feature 간의 비교를 위해서는 Error bar/ LDA score 등을 통해 보여주는 것이 나은 것 같다.
차등적으로 풍부한 개체를 판별하기 위해서는 기본적인 통계분석(wilcoxon, kruskal wallis test)과 LefSe, ANCOM-BC, ALDEx2, DESeq2등의 방법 등이 사용된다.
이 글에서는 그중에서 ALDEx2의 분석을 수행하고 결과를 시각화해 보자.
| 데이터
- 예제 데이터 : qiime2 moving picture tutorial의 데이터 (phyloseq format)
- 추가 데이터 : PICRUSt2 결과물 (pred_metagenome_unstrat.tsv)
| ggpicrsut2로 aldex2 분석하고 pathway 데이터 얻기
library(phyloseq) # 마이크로바이옴 데이터 분석 및 시각화
library(tidyverse) # R의 데이터 핸들링
library(readr) # 파일 읽어오기
library(ggpicrust2) # PICRUSt2결과 처리 및 시각화
library(ALDEx2) # Aldex2 분석
library(pheatmap) # heatmap 시각화
1. 데이터 가져오기
## phyloseq 개체 읽어오기, 위치는 본인의 파일 위치에 따라 수정
ps <- readRDS("./ps.rds")
## ggpicrust2의 패키지를 이용해서 KO abundance를 kegg pathway abundance로 변환
kegg_abundance <-
ko2kegg_abundance(
"./picrust_result/KO_metagenome_out/pred_metagenome_unstrat.tsv"
)
kegg_abundance
dim(kegg_abundance)
# 243 34
2. Gut와 tongue 데이터만 추출하기
# 혀와 장 데이터만 추출
ps_sub <- subset_samples(ps, body.site %in% c("gut", "tongue"))
# phyloseq-class experiment-level object
# otu_table() OTU Table: [ 770 taxa and 17 samples ]
# sample_data() Sample Data: [ 17 samples by 8 sample variables ]
# tax_table() Taxonomy Table: [ 770 taxa by 7 taxonomic ranks ]
# phy_tree() Phylogenetic Tree: [ 770 tips and 768 internal nodes ]
# 샘플의 메타데이터 추출
metadata <- sample_data(ps_sub) %>% data.frame()
# kegg_abundance에서도 장과 혀 샘플만 추출
metadata <- metadata %>% rownames_to_column("SampleID")
metadata <- metadata[metadata$body.site %in% c("gut", "tongue"), ]
kegg_abundance2 <- kegg_abundance[, metadata$SampleID]
dim(kegg_abundance2)
# 243 17
3. Aldex2 분석
Group <- factor(metadata$body.site)
Level <- levels(Group)
## ALDEx2 분석 ##
# 1) Aldex2분석을 위해 OTU read count 데이터를 CLR로 normalization
ALDEx2_object <- ALDEx2::aldex.clr(round(kegg_abundance2), metadata$body.site)
# clr 변환은 각 샘플의 모든 taxa값이 0이며, taxa간의 상대적 비율 값이 그대로 보존된다
# 2) Aldex 통계 분석
ALDEx2_results <- ALDEx2::aldex.ttest(ALDEx2_object, paired.test = FALSE, verbose = FALSE)
# t.test결과와 wilcoxon rank sum test결과를 반환한다
# 각 taxa가 그룹내에서 유의한 수준으로 차이가 '있는지 없는지' 판별
# 3) Effect size계산
ALDEx2_effect <- ALDEx2::aldex.effect(ALDEx2_object)
# 각 taxa가 차이가 있다면(유의하다면), '얼마나' 차이가 있는지 판별
# 4) 결과물 정리
p_values_df <- data.frame(
feature = rep(rownames(ALDEx2_results), 2),
method = c(rep("ALDEx2_Welch's t test", nrow(ALDEx2_results)),
rep("ALDEx2_Wilcoxon rank test", nrow(ALDEx2_results))),
group1 = rep(Level[1], 2 * nrow(ALDEx2_results)),
group2 = rep(Level[2], 2 * nrow(ALDEx2_results)),
p_values = c(ALDEx2_results$we.ep, ALDEx2_results$wi.ep),
effect = ALDEx2_effect$effect)
# 5) P- value값 보정 (fdr 사용시)
# BH사용시, 기존 p_values_df에 wi.eBH를 사용
adjusted_p_values <- data.frame(
feature = p_values_df$feature,
p_adjust = p.adjust(p_values_df$p_values, method = "fdr"))
# 6) 결과물 최종 정리
ko_AlDEx_result <- cbind(p_values_df, p_adjust = adjusted_p_values$p_adjust)
ko_AlDEx_result
세 군 이상에서는 method가 ALDEx2_Kruskal-Wallace test 혹은 ALDEx2_glm test 이 사용됩니다.
그러나 세 군간의 비교에서 effect 값을 산정할 수는 없습니다.
4. Pathway 정보 가져오기
# "ALDEx2_Wilcoxon rank test"결과 중에서 p-value값 보정이 유의한 것만 추출
ko_AlDEx_result2 <- ko_AlDEx_result[ko_AlDEx_result$method == "ALDEx2_Wilcoxon rank test" &
ko_AlDEx_result$p_adjust < 0.05, ]
dim(ko_AlDEx_result2) # 91 7
# ggpicrust2의 annotation함수를 통해 이 kegg pathway데이터 가져오기
ko_ann <-pathway_annotation(pathway = "KO",
daa_results_df = ko_AlDEx_result2,
ko_to_kegg = TRUE)
dim(ko_ann) # 91 11pathway_annotation 함수는 보정된 p-value값이 0.05 이하인 경우에만 kegg data를 가져온다.
kegg database에 접속하여 데이터를 다운로드하기 때문에, 만약 본인이 가진 데이터의 ko가 1000개 이상일 경우 일부만 다운로드되거나, 하루 이상의 시간이 걸릴 수 있다.
+) annotation 데이터가 아예 존재하지 않는 ko도 있다. 즉, input data와 output data의 크기가 달라질 수 있음
+) MetaCyc 데이터도 가능하다. 자세한 내용은 ggpicrust2 깃허브 참고
| Heatmap 그리기
1. 기본 heatmap 그리기
- heatmap에 그려지는 값은 각 taxa의 relative abundance값이다.
- 우리는 여기에 annotation bar로 p-value, effect size값을 추가해 보자.
# Pathway 1,2,3을 하나의 이름으로 표시하기
pathway_table <- ko_ann %>%
na.omit %>% # 44
mutate(Level_123 = paste0(pathway_class," - ", pathway_map))
rownames(pathway_table) <- pathway_table$feature
# annotation이 완료된 kegg만 추출하기
kegg_abundance3 <- kegg_abundance2[pathway_table$feature, ]
# kegg abundance를 relativer abundance로 바꾸기
col_sums <- colSums(kegg_abundance3)
relative_abundance <- t(t(kegg_abundance3) / col_sums)
# heatmap을 위한 데이터 생성
table <- merge(pathway_table, relative_abundance, by = "row.names")
rownames(table) <- table$Level_123
# heatmap 그리기
library(pheatmap)
p <- pheatmap(table[, 14:30])
p
하지만 이 그림을 보면 어떤 샘플이 gut인지 tongue인지 잘 보이지 않는다.
이를 구분하기 위해 heatmap의 annotation 기능을 추가하고, 각 kegg pathway별로 abundance값을 표준화하여서 다시 그려보자.
2. 보정한 heatmap 그리기
annotation_row <- data.frame(
row.names = rownames(table),
p_adjust = table$p_adjust,
Effect = table$effect
)
metadata$body.site %>% table
annotation_row <- annotation_row[order(annotation_row$Effect, decreasing = F), ,
drop = FALSE # rownames 사라지는걸 막음
]
annotation_col <- data.frame(
row.names = metadata$SampleID,
Site = metadata$body.site
)
annotation_col$Site <- factor(annotation_col$Site , levels = c("gut", "tongue") )
annotation_col <- annotation_col[order(annotation_col$Site, decreasing = FALSE), , drop = FALSE ]
ann_colors = list(
p_adjust = colorRampPalette(c("white", "green"))(100),
Effect = colorRampPalette(c("red", "white", "blue"))(100),
Site = c("gut" = "orange","tongue"="purple")
)
df <- table[, 14:30]
colnames(df)
df <- df[rownames(annotation_row), rownames(annotation_col)]
plot <- pheatmap::pheatmap(mat = as.matrix(df),
color = colorRampPalette(c("blue", "white", "red"))(100),
annotation_col = annotation_col,
annotation_row = annotation_row,
annotation_colors = ann_colors,
scale = "row", # 행별로 정규화
cluster_rows = F,
cluster_cols = F,
gaps_col = 8, # 8번째 열에서 heatmap 분리
gaps_row = 21,# 21번째 행에서 heatmap 분리
legend = T,
border_color=NA)
각 그룹에서 유의하게 차등적인 발현을 보이는 kegg pathway를 관찰해 보았다.
위와 같이 각 부위에서 확연하게 차이나는 PICRUSt2그룹을 관찰할 수 있다.
| Effect Size 그리기
p-value값은 임계치를 넘는지 넘지 않는지를 판별한다. 그러나 과학적으로 중요한 것은 이 값이 '얼마나' 중요한지 양적인 판단을 내리는 것이다. 이 때문에 ALDEx2의 저자는 p-value값 보다 effect size에 집중하라고 한다.
heatmap으로 aldex2 결과와 relative abundance의 차이가 동일한지 확인 한 수, effect size의 절댓값이 3 이상인 값을 확인해 보자.
# 순서 정렬
lev = pathway_table[order(pathway_table$effect), "Level_123"]
pathway_table$Level_123 <- factor(pathway_table$Level_123, level = lev)
# ggplot
pathway_table %>%
mutate(Group = if_else(effect>0, "gut", "tongue")) %>%
filter(abs(effect)>3) %>%
ggplot(aes(x = effect, y = Level_123, fill = Group)) +
geom_col() + theme_classic() +
labs(y = NULL, x = 'Effect zise')+
scale_x_continuous(limits = c(-7, 7), breaks = seq(-7,7,by=2))
이미지 출처
-Kim, HJ., Kim, J.J., Myeong, N.R. et al. Segregation of age-related skin microbiome characteristics by functionality. Sci Rep 9 , 16748 (2019).
- https://doi.org/10.1038/s41598-019-53266-3/ https://picrust.github.io/picrust/tutorials /stamp_tutorial.html
- Bailén M, Tabone M, Bressa C, Lominchar MGM, Larrosa M, González-Soltero R. Unraveling Gut Microbiota Signatures Associated with PPARD and PARGC1A Genetic Polymorphisms in a Healthy Population. Genes (Basel). 2022 Feb 1;13(2):289. doi: 10.3390/genes13020289. PMID: 35205333; PMCID: PMC8871880.
