작성: 2023-08-25
Metacoder란?
- 공식 튜토리얼 : https://grunwaldlab.github.io/metacoder_documentation/workshop--05--plotting.html
- 논문 : Foster, Z. S., Sharpton, T. J., & Grünwald, N. J. (2017). Metacoder: An R package for visualization and manipulation of community taxonomic diversity data. PLoS computational biology, 13(2), e1005404. https://doi.org/10.1371/journal.pcbi.1005404
- 인용수: 498(2023.08.25 기준)
- R에서 taxonomy의 다양성 데이터를 시각화하고 관리하는 패키지를 말한다. 대표적으로 아래와 같은 이미지를 만들어낸다.
- metacoder에서 시각적으로 표현하고자 하는 바는 아래와 같다
1. Node(원)의 크기는 OTU count(혹은 abundance)
2. 색(e.g. 민트색 - 남자, 갈색 - 여자)은 각 그룹에서 차등적으로 풍부한 taxa를 나타낸다. 이는 Wilcoxon Rank sum test를 통해 정해진다.
3. 각 그룹에서 풍부한 taxa가 어떤 계통수에 위치하는지 나타낸다.
- 이와 마찬가지로 계층적 구조를 가진 데이터(KEGG pathway 등)를 시각화하기에도 적합하다.
- 어떻게 보면 Krona plot 혹은 lefse의 output으로 활용되는 Cladogram과 비슷한 용도라고 말할 수 있다.
예제 데이터
- qiime2 moving-picrutres data
- 사람의 신체 부위 4곳을 샘플링한 데이터
- phyloseq 형태로 만들어 놓음
1. import data
library(phyloseq)
library(metacoder)
library(ggplot2)
library(dplyr)
ps <- readRDS("~/ps.rds") %>%
transform_sample_counts(function(x) x/sum(x) ) %>%
subset_samples(body.site %in% c("gut", "tongue"))
ps # 770 taxa and 17 samples
# read가 0인 샘플, taxa제거
ps.f = prune_samples(sample_sums(ps)>0, ps)
ps.f = prune_taxa(rowSums(otu_table(t(ps.f))) > 0, ps.f) # 위 데이터에선 에러
Heat tree 그리는 방법
1. Down load pacakge
- CRAN에서 다운로드
install.packages("metacoder")- Github에서 개발자 버전 다운로드 (가장 최근 버전/ 버그가 있을 수도 있고, 나타난 버그가 고쳐진 것도 있음 )
install.packages("devtools")
devtools::install_github("grunwaldlab/metacoder")
2. Phyloseq to metacoder object
obj <- parse_phyloseq(ps.f)
heat_tree(obj, output_file = "~/basic1.png")
기본으로 ASV레벨의 계통수를 출력할 수 있다. 여기에 각 taxa에 따른 abundance를 node의 크기에 적용하고, 각 taxa의 이름을 띄어주자.
heat_tree(obj,
node_label = taxon_names,
node_size = n_obs,
output_file = "~/basic2.png")
이.. 이름이 너무 작으니 일단 Order level에서 살펴보자.
Order level
obj %>%
filter_taxa(taxon_ranks == "Order", supertaxa = TRUE) %>%
heat_tree(node_label = taxon_names,
node_size = n_obs,
node_color = n_obs,
node_color_range = diverging_palette(),
output_file = "~/basic3.png")- 이왕 살펴보는 김에 node color도 정해보자
One group
- 그런데 Gut에서만 풍부하게 나타나는 taxa를 강조하고 싶을 때는?
(이왕 그려주는 김에 layout을 달리해서 그려보자)
obj$data$tax_data <- calc_obs_props(obj, "tax_data") # 각 taxa 계산
obj$data$tax_abund <- calc_taxon_abund(obj, "otu_table") # 샘플의 abundance 계산
obj$data$tax_occ <- calc_n_samples(obj, "tax_abund",
groups = obj$data$sample_data$body.site,
cols = obj$data$sample_data$sample_id)
set.seed(42)
obj %>%
filter_taxa(taxon_ranks == "Order", supertaxa = TRUE) %>%
heat_tree(node_label = taxon_names,
node_size = n_obs,
node_color = gut, # tax_occ의 계산에 있는 column을 뜻함
node_size_axis_label = "OTU count",
node_color_axis_label = "Samples with reads",
layout = "davidson-harel",
initial_layout = "reingold-tilford",
output_file = "~/ex1.png")
- reingold-tilford옵션을 제외하고 그려보자
obj %>%
filter_taxa(taxon_ranks == "Order", supertaxa = TRUE) %>%
heat_tree(node_label = taxon_names,
node_size = n_obs,
node_color = tongue,
node_size_axis_label = "OTU count",
node_color_axis_label = "Samples with reads",
layout = "da", # = davidson-harel
output_file = "~/ex2.png")
Two groups
- 두 그룹 간의 abudance가 다른 taxa를 강조하고 싶을 땐? compare_groups로 abundance를 계산 후 시각화 한다.
(이왕 그릴 거 Genus level까지 나타내보자)
- 이때 wilcoxon rank sum test가 사용된다. 다른 test방법(lefse, deseq2 등)을 사용하려면, 추가적인 계산이 필요하다.
obj$data$diff_table <- compare_groups(obj, data = "tax_abund",
cols = obj$data$sample_data$sample_id,
groups = obj$data$sample_data$body.site)
print(obj$data$diff_table)
obj %>%
filter_taxa(taxon_ranks == "Genus", supertaxa = TRUE) %>%
heat_tree(node_label = taxon_names,
node_size_axis_label = "Number of OTUs",
node_size = n_obs,
node_color_axis_label = "Mean difference",
node_color = mean_diff,
node_color_range= c("#E31A1C","grey90", "#1F78B4"),
layout = "davidson-harel",
initial_layout = "reingold-tilford",
output_file = "~/ex3.png")
붉은색이 혀, 푸른색이 장에서 더 많은 taxa를 나타낸다.
DESeq2
- deseq2에서는 read count를 사용함으로, relative abundance로 바꾸는 단계를 제외한다.
ps <- readRDS("C:/Users/82105/Downloads/ps.rds") %>%
subset_samples(body.site %in% c("gut", "tongue"))
ps.f = prune_samples(sample_sums(ps)>0, ps)
obj2 = parse_phyloseq(ps.f)
obj2$data$tax_data <- calc_obs_props(obj2, "tax_data") # 각 taxa 계산
obj2$data$tax_abund <- calc_taxon_abund(obj2, "otu_table") # 샘플의 abundance 계산
obj2$data$tax_occ <- calc_n_samples(obj2, "tax_abund",
groups = obj2$data$sample_data$body.site,
cols = obj$data$sample_data$sample_id)
obj2$data$tax_data
obj2$data$tax_abund
obj2$data$tax_occ
obj2$data$diff_table <- calc_diff_abund_deseq2(obj2, data = "tax_abund",
cols = obj$data$sample_data$sample_id,
groups = obj$data$sample_data$body.site)
print(obj2$data$diff_table)
obj2 %>%
filter_taxa(taxon_ranks == "Genus", supertaxa = TRUE) %>%
heat_tree(node_label = taxon_names,
node_size = n_obs,
node_color_axis_label = "Log2 fold change",
node_size_axis_label = "Number of OTUs",
node_color = ifelse(is.na(padj) | padj > 0.05, 0, log2FoldChange),
node_color_range= c("#E31A1C","grey90", "#1F78B4"),
layout = "davidson-harel",
initial_layout = "reingold-tilford",
output_file = "~/ex5.png")
참고
- https://blog.csdn.net/woodcorpse/article/details/121867820
(중국에서 따거들이 가끔 마이크로바이옴 관련 튜토리얼을 상세하게 올려준다.. 한국보다 훨씬 많이 연구되어 있는 것 같다)
- metacoder group : https://groups.google.com/g/metacoder-discussions (혹시 버그가 생긴다면 여기로 제보/ 혹은 깃헙)
엄청 간단한 프로그램인 줄 알았는데, 생각보다 옵션이 많아서 당황했다.
혹시 나처럼 당황했던 사람들을 위해 작성해 보았다.
교수님께서 발표용으로 아주 적합한 figure라고 말씀해 주셨다. 물론 마이크로바이옴을 아는 사람들에게는 일반적인 그림이 더 잘 이해되겠지만, 대다수의 청중을 겨냥한 발표용으로는 설명하기 매우 쉬운 그림이라고 말씀하셨다.
인용수도 꽤 높은 편이라 논문에 사용하기에도 적합해 보인다. 그러나 main figure로 사용하기에는 가독성이 떨어져서 supplementary로 사용하는 것이 좋을 것 같다.
