작성: 2024-10-21
수정: 2024-01-08
side panel 은 beta diversity를 나타내는 PCA, PCoA 그림의 추가 정보를 나타냅니다. PCoA의 side panel 은 각 샘플의 dissimilarity 간의 분포를 축소해서 나타내는 PCoA 의 축 1과 2를 boxplot이나 density plot을 사용하여 나타냅니다. 5가지 방법(patchwork, MicrobiomeStat, ggpurb, cowplot, ggExtra) ggExtra)을 통해 아래와 같은 그림을 그려봅시다.
예제 데이터
언제나 qiime2의 moving picture 데이터를 사용합니다.
아래 예제 코드에서는 신체의 각 부위에 따른 마이크로바이옴의 분포를 비교합니다.
library(phyloseq)
library(tidyverse)
library(devtools)
library(vegan)
library(glue)
library(ggtext)
library(ggpubr)
library(patchwork)
1. ggplot + patchwork 패키지(수동)
- R을 사용하는 사람이라면 ggplot은 다 사용하실 겁니다. 가장 기본적인 패키지로 그림을 만들어 봅시다!
- 단점이라면, 각 PCoA의 x축과 y축을 고정해서 사용해야 한다는 것입니다. 저는 lims(x = c(-0.8, 0.8), y = c(-0.8, 0.8)) 로 고정하였습니다.
#### 1. 데이터 변수 입력 ####
phyloseq = readRDS("./ps.rds") # phyloseq
index = "bray" # distance
seed = 42 # seed값 고정
type = "body.site" # 비교 group
plot = "PCoA" # 차원축소
SampleID = "SampleID"
shap = NULL # point 모양 group
type_col = c("#E31A1C", "#1F78B4", "#4D4D4D", "green") # 색 지정
#### 2. PCoA를 위한 데이터 변환 ####
set.seed(seed)
x.dist <- phyloseq::distance(phyloseq, method = index) # distance계산
ord <- ordinate(phyloseq, plot, index)
# 3. stat - PERMANOVA
myfunc <- function(v1) {
deparse(substitute(v1))
}
dist <- myfunc(x.dist)
a <- adonis2(as.formula (glue("{dist} ~ {type}")), data=data.frame(sample_data(phyloseq)),
permutations=9999, method=index)
Perm.p <- a$`Pr(>F)`[1]
#### 4. data 변환 for PCoA ####
mat <- ord$vectors[, 1:2] %>% as.data.frame() %>% # PCoA 1, 2 추출
mutate(SampleID = rownames(.)) %>%
arrange(SampleID)
meta <- sample_data(phyloseq) %>% data.frame() %>%
arrange(SampleID)
pcoa_df <- inner_join(meta, mat, by = "SampleID") # meta data + PCoA 1,2
## pcoa_df는 아래에서 계속 사용할 데이터입니다 ##
PC1 <- round(ord$values["Relative_eig"][1,]*100, 1)
PC2 <- round(ord$values["Relative_eig"][2,]*100, 1)
#### 5. PCoA ####
## main plot
main.plot <- pcoa_df %>%
ggplot2::ggplot(aes(x = Axis.1, y=Axis.2)) +
ggplot2::geom_vline(xintercept = 0, colour = "grey80") +
ggplot2::geom_hline(yintercept = 0, colour = "grey80") +
ggplot2::geom_point(aes_string(shape = shap, color=type),alpha = 0.7, size=2.5) +
lims(x = c(-0.8, 0.8), y = c(-0.8, 0.8)) +
ggplot2::stat_ellipse(aes_string(color= type) ) +
ggplot2::theme_test() +
ggplot2::scale_color_manual(values = type_col) +
ggplot2::labs(
y =paste0("PCoA2 (", PC2, "%)"),
x =paste0("PCoA1 (", PC1, "%)")) +
ggplot2::theme(plot.caption = element_text(hjust = 0)) +
ggplot2::theme(plot.caption = element_markdown(),
# aspect.ratio=1,
legend.position = "bottom",
legend.title = element_blank())+
theme(legend.position="none",plot.margin=unit(c(0,0,0,0),"points"))+
ggtext::geom_richtext(label.color = NA, size = 4, fill = NA,
hjust = 0, vjust =0, x = -Inf, y = -Inf,
label= paste0( "**PERMANOVA** *p*-value=", Perm.p ) )
main.plot
main plot을 기준으로 side panel 을 boxplot과 density 두 종류로 만들어서 붙여 봅시다
- box plot
# theme 설정
theme0 <- function(...) theme( legend.position = "none",
panel.background = element_blank(),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
panel.margin = unit(0,"null"),
axis.ticks = element_blank(),
axis.text.x = element_blank(),
axis.text.y = element_blank(),
axis.title.x = element_blank(),
axis.title.y = element_blank(),
axis.ticks.length = unit(0,"null"),
axis.ticks.margin = unit(0,"null"),
panel.border=element_rect(color=NA),...)
# box plot
pcoa1.boxplot <-
ggplot2::ggplot(pcoa_df) +
ggplot2::geom_boxplot(ggplot2::aes(x = !!sym(type),y = Axis.1, fill = !!sym(type)), color = "black", alpha = 0.5, show.legend = F) +
ggplot2::scale_fill_manual(values = type_col) +coord_flip() +
lims(y = c(-0.8, 0.8)) +
theme_void()
pcoa2.boxplot <-
ggplot2::ggplot(pcoa_df) +
ggplot2::geom_boxplot(ggplot2::aes(x = !!sym(type),y = Axis.2, fill = !!sym(type)), color = "black", alpha = 0.5, show.legend = F) +
ggplot2::scale_fill_manual(values = type_col) +
lims(y = c(-0.8, 0.8)) +
theme_void()
pcoa2.boxplot
pcoa1.boxplot + plot_spacer() + main.plot + pcoa2.boxplot +
plot_layout(ncol = 2, nrow = 2, widths = c(4, 0.5), heights = c(0.5, 4))
- 각 PCoA에 따라 분포된 boxplot은 각 부위에 대한 dissimilarity의 축소된 데이터를 나타냅니다.
- 빨간색이 장 샘플, 초록색이 혀 샘플, 하늘색과 회색이 손바닥 샘플입니다.
- 요약하자면, PCoA1에서 장샘플의 dissimilarity은 낮은 편이며, 혀의 dissimilarity는 높은 편입니다. 이는 장샘플보다 혀의 샘플이 샘플마다 이질적이고 다양성이 높다는 것을 뜻합니다.
- density plot
# density plot
pcoa1.density <-
ggplot(pcoa_df,aes(x=Axis.1 ,fill=!!sym(type))) +
geom_density(alpha=0.5) +
scale_x_continuous(breaks=NULL,expand=c(0.02,0)) +
scale_y_continuous(breaks=NULL,expand=c(0.02,0)) +
theme_bw() +
lims(x = c(-0.8, 0.8)) +
theme0(plot.margin = unit(c(1,0,-0.48,2.2),"lines"))
pcoa1.density
pcoa2.density <-
ggplot(pcoa_df,aes(x=Axis.2 ,fill=!!sym(type))) +
geom_density(alpha=0.5) +
scale_x_continuous(breaks=NULL,expand=c(0.02,0)) +
scale_y_continuous(breaks=NULL,expand=c(0.02,0)) +
theme_bw() +
lims(x = c(-0.8, 0.8)) +
theme0(plot.margin = unit(c(0,1,1.2,-0.48),"lines")) +
coord_flip()
pcoa2.density
pcoa1.density + plot_spacer() + main.plot + pcoa2.density +
plot_layout(ncol = 2, nrow = 2, widths = c(4, 0.5), heights = c(0.5, 4))
- density는 샘플의 어디에 집중되어 있는지를 나타냅니다.
- 가장 좋은 결과는 각 부위마다 density의 정점이 나누어지는 것입니다.
- 위 데이터를 보면 PCoA 1과 2에서 장과 혀는 완전히 다른 density를 보이지만, 손바닥은 매우 겹쳐서 나타납니다.
2. MicrobiomeStat 패키지
- 통계분석부터 시각화까지 올인원을 원하신다면 MicrobiomeStat도 좋은 패키지입니다.
- #2. stat에서 pairwise 분석 통계값을 얻을 수 있습니다.
library(microbiome)
devtools::install_github("cafferychen777/MicrobiomeStat")
library(MicrobiomeStat)
library(aplot)
#### 1. import data ####
ps <- readRDS("./ps.rds")
data.obj <- mStat_convert_phyloseq_to_data_obj(ps)
#### 2. stat ####
beta.df <- generate_beta_test_single(
data.obj = data.obj,
dist.obj = NULL,
t.level = "2",
group.var = "body.site",
dist.name = c('BC', 'Jaccard')
)
beta.df$p.tab
beta.df$aov.tab
#### 3. PCoA ####
data.obj$meta.dat$SubjectID <- rownames(data.obj$meta.dat)
p <- generate_beta_ordination_single(
data.obj = data.obj,
dist.obj = NULL,
pc.obj = NULL,
subject.var = "SubjectID",
group.var = "body.site",
strata.var = "subject",
dist.name = c("BC"),
base.size = 10,
theme.choice = "bw",
custom.theme = NULL,
palette = NULL,
pdf = TRUE,
file.ann = NULL,
pdf.wid = 11,
pdf.hei = 8.5,
)
p
3. ggpurb::ggscatterhist
- ggpurb 패키지답게 깔끔한 그림을 얻을 수 있습니다.
# histogram
ggscatterhist(
pcoa_df, x = "Axis.1", y = "Axis.2",
color = "body.site", # comment out this and last line to remove the split by species
margin.plot = "histogram", # I'd suggest removing this line to get density plots
margin.params = list(fill = "body.site", color = "black", size = 0.2)
)
# density
ggscatterhist(
pcoa_df, x = "Axis.1", y = "Axis.2",
color = "body.site", # comment out this and last line to remove the split by species
margin.plot = "density", # I'd suggest removing this line to get density plots
margin.params = list(fill = "body.site", color = "black", size = 0.2)
)
4. cowplot::axis_canvas
library("cowplot")
# density
xdens <-
axis_canvas(main.plot, axis = "x") +
geom_density(data = pcoa_df, aes(x = Axis.1, fill = body.site, colour = body.site), alpha = 0.3)
ydens <-
axis_canvas(main.plot, axis = "y", coord_flip = TRUE) +
geom_density(data = pcoa_df, aes(x = Axis.1, fill = body.site, colour = body.site), alpha = 0.3) +
coord_flip()
main.plot %>%
insert_xaxis_grob(xdens, grid::unit(1, "in"), position = "top") %>%
insert_yaxis_grob(ydens, grid::unit(1, "in"), position = "right") %>%
ggdraw()
# box plot
xbox <-
axis_canvas(main.plot, axis = "x") +
geom_boxplot(data = pcoa_df, aes(x = Axis.1, fill = body.site, colour = body.site), alpha = 0.3)
ybox <-
axis_canvas(main.plot, axis = "y", coord_flip = TRUE) +
geom_boxplot(data = pcoa_df, aes(x = Axis.1, fill = body.site, colour = body.site), alpha = 0.3) +
coord_flip()
main.plot %>%
insert_xaxis_grob(xbox, grid::unit(1, "in"), position = "top") %>%
insert_yaxis_grob(ybox, grid::unit(1, "in"), position = "right") %>%
ggdraw()
5. ggExtra::ggMarginal
가장자리 plot의 가장 대표적인 패키지이죠. ggExtra입니다. 깔아 두면 언젠가는 한 번쯤 쓰는 날이 올 것입니다.
main plot의 데이터를 사용하여 그려주기 때문에, 따로 색이나 x축, y축의 값을 자동으로 지정해 줍니다.
# boxplot
ggExtra::ggMarginal(
main.plot,
type =c('boxplot'),
margins = 'both',
size = 3.5,
groupColour = F,
groupFill = T
)
# density
ggExtra::ggMarginal(
main.plot,
type =c('density'),
margins = 'both',
size = 3.5,
groupColour = F,
groupFill = T
)
상세 옵션이 적어서 처음 그리는 분도 부담이 없고, output도 가장 깔끔합니다.
여기에 dissimilarity 간의 통계 분석도 추가하고 싶다면?
ggplot으로 가장자리 plot을 먼저 그리고 ggpurb으로 통계값을 추가합니다.
library(pracma )
comp_ <- pcoa_df$body.site %>% unique() %>% as.vector()
comp_pair <- combn(comp_, 2)
combined_list <- lapply(1:ncol(comp_pair), function(i) comp_pair[, i])
ybox <-
ggplot(pcoa_df, aes(x = body.site, y = Axis.2)) +
geom_boxplot(aes(fill = body.site, colour = body.site), alpha = 0.3)+
scale_fill_manual(values = type_col) +
scale_color_manual(values = type_col) +
lims(y =c(-0.9, 0.9)) +
stat_compare_means(method = "wilcox.test", tip.length=0.02,
label = "p.signif",
comparisons = combined_list) +
theme_void() +
theme(legend.title = element_blank(),
legend.position = 'none',
strip.background = element_blank())
xbox <-
ggplot(pcoa_df, aes(x = body.site, y = Axis.1)) +
geom_boxplot(aes(fill = body.site, colour = body.site), alpha = 0.3)+
scale_fill_manual(values = type_col) +
scale_color_manual(values = type_col) +
lims(y =c(-0.9, 0.9)) +
stat_compare_means(method = "wilcox.test", tip.length=0.02,
label = "p.signif",
comparisons = combined_list) +
theme_void() +
theme(legend.title = element_blank(),
legend.position = 'none',
strip.background = element_blank())+
coord_flip()
마지막으로 cowplot으로 그림을 합쳐주면, patchwork나 ggarrange보다 더 깔끔하게 그림을 합칠 수 있습니다.
pcoa <- main.plot %>%
cowplot::insert_yaxis_grob(ybox, grid::unit(1, "in"), position = "right") %>%
cowplot::insert_xaxis_grob(xbox, grid::unit(1, "in"), position = "top") %>%
ggdraw()
pcoa
- Wilcoxon rank sum test의 p-value를 추가해 보았습니다.
- 양 손바닥을 제외하고 모두 유의한 차이를 보입니다.
교수님한테 한소리 들을 것 같은 그림이 되었네요.. 역시 figure은 깔끔한 게 답인 것 같습니다.
참고
반응형
