수정: 24.10.15
- 예제데이터와 분석 코드 포함
1. Beta diversity란?
마이크로바이옴 연구에서 beta diversity란 미생물 군집의 다양성을 분석하기 위한 단계이며, 이는 각 샘플 혹은 샘플이 포함된 집단(질환군 vs. 정상대조군) 사이의 차이를 측정하여 얻어진다. 가장 많이 사용되는 방법은 각 비교 집단 간 유사성 혹은 비유사성을 사용하며, 미생물의 풍부도, 미생물의 계통 간 거리가 사용되기도 한다.
2. Beta diversity 분석 방법
가장 대표적인 접근 방법은 총 4가지로 나눌 수 있다.
- Bray-Curtis dissmilarity: 샘플 간 미생물 분포의 차이에 따라 측정
- Jaccard index: 집단 간 미생물의 존재/부존재에 따라 측정
- weighted UniFac: 미생물의 계통 거리에 영향을 받으면서, 동시에 분포에 따라 측정
- unweighted UniFrac: 미생물의 계통 거리에 영향을 받으면서, 집단 간 미생물의 존재/부존재에 따라 측정
3. Beta diversity 분석 도구
beta diversity측정 도구를 사용해서 이렇게 각 샘플 간의 차이를 확인한 후, 통계 분석을 통해 이들이 유의한 차이가 나는지 알아내야 한다. 이때 우리는 PERMANOVA, ANOVA, ANOSIM 등의 도구를 사용한다.
또한 각 샘플 간의 거리를 시각화하기 위해 PCoA, PCA, NMDS 같은 차원축도 도구를 사용하여 2차원 평면에 투사하기도 한다.
4. PERMANOVA란?
PERMANOVA란 비모수적 다변수 순열검정(non-parametric multivariate statistical permutation test)을 말한다. R의 환경데이터 분석 프로그램인 "vegan" 패키지에서 `vegan::adonis2()`로 제공하고 있다.
비모수 검정임으로 정규성, 등분산성의 만족이 필요하지 않고, 순열을 사용하여 계산된다. 그러므로 샘플 간의 거리(예: Bray-Curtis 거리)를 기반으로 측정된 거리행렬에서 바로 적용이 가능하다.
1) PERMANOVA의 검정
귀무가설 (H0): 그룹 내 중심값(centroid)과 산포도(dispersion)는 모든 그룹에서 동일하다.
The centroids of the groups, as defined in the space of the chosen resemblance measure, are equivalent for all groups.
- 만약 귀무가설이 사실이라면, 데이터의 중심값의 관찰된 차이는 샘플을 각 그룹에서 무작위로 얻어낸 값과 유사하다. 즉, 그룹 간 차이는 우연으로 가정한다.
Thus, if H0 were true, any observed differences among the centroids in a given set of data will be similar in size to what would be obtained under random allocation of individual sample units to the groups (i.e., under permutation).
- 귀무가설이 기각되면, 개체의 중심값 혹은 산포도가 그룹 간에 다름을 뜻한다.
A rejection of the null hypothesis means that either the centroid and/or the spread of the objects is different between the groups.
계산은 중심값과 각 개체 간의 거리 제곱합을 개체의 수로 나눈 것과 같다.
2) PERMANOVA의 검정 수행하기
오늘도 어김없이 QIIME2의 moving-picture예제 데이터를 사용하겠습니다. 이는 사람의 혀, 양 손바닥, 장의 마이크로바이옴 분포를 담고 있다.
library(phyloseq)
library(vegan)
library(tidyverse)
ps <- readRDS("./ps.rds")
먼저, bray-curtis 비유사성을 기반으로 거리행렬을 만든 다음, PERMANOVA를 수행해 보겠습니다.
set.seed(42)
braycurtis <- phyloseq::distance(ps, method = "bray")
set.seed(42)
a <- adonis2(braycurtis~body.site,
data=data.frame(sample_data(ps)),
permutations=9999
)
a
# Permutation test for adonis under reduced model
# Terms added sequentially (first to last)
# Permutation: free
# Number of permutations: 9999
#
# adonis2(formula = distance ~ body.site, data = data.frame(sample_data(ps)), permutations = 9999)
# Df SumOfSqs R2 F Pr(>F)
# body.site 3 5.1913 0.39988 6.6633 1e-04 ***
# Residual 30 7.7909 0.60012
# Total 33 12.9823 1.00000
R2값은 각 변수가 전체 변동에 몇 퍼센트를 차지하고 있는지 나타낸다. 이는 0에서 1 사이 값을 가집니다.
위 결과를 보면 body.site가 마이크로바이옴 전체 변동에 약 40% 를 설명하고 있습니다.
이때 R 값이 0.01 이상(1% 이상)이면 효과를 미치는 변수라고 가정합니다.
또한 p-value값이 0.001 이하이므로 신체에 따라 마이크로바이옴의 변동은 유의한 차이를 나타내고 있습니다.
이때 adonis 함수에 `method = "bray"`라는 옵션이 존재합니다. 이게 기본 값이라서 input 데이터로 거리행렬을 넣으면 중복 계산되는 것이 아니냐 헷갈릴 수 있다.
실제로 관련 글이 커뮤니티에 질문으로 올라와 있습니다.
1. https://stats.stackexchange.com/questions/431433/in-vegan-adonis-what-does-method-bray-do-when-i-pass-in-a-dista nce-matrix-i
2. https://github.com/vegandevs/vegan/issues/330
결과적으로만 보면, 거리행렬을 input으로 가져가면 위 변수는 무효화됩니다.
이를 직접 비교해 봅시다.
set.seed(42)
b <- adonis2(data.frame(t(otu_table(ps)))~body.site,
data=data.frame(sample_data(ps)),
permutations=9999,method = "bray"
)
b
# Permutation test for adonis under reduced model
# Terms added sequentially (first to last)
# Permutation: free
# Number of permutations: 9999
#
# adonis2(formula = data.frame(t(otu_table(ps))) ~ body.site, data = data.frame(sample_data(ps)), permutations = 9999, method = "bray")
# Df SumOfSqs R2 F Pr(>F)
# body.site 3 5.1913 0.39988 6.6633 1e-04 ***
# Residual 30 7.7909 0.60012
# Total 33 12.9823 1.00000
OTU table을 input으로 넣고, vegan::adonis2 함수 내의 변수를 bray로 설정하면, 초기 방식과 동일한 결과를 생산합니다.
set.seed(42)
c <- adonis2(distance~body.site,
data=data.frame(sample_data(ps)),
permutations=9999,
method = "jaccard"
)
c
# Permutation test for adonis under reduced model
# Terms added sequentially (first to last)
# Permutation: free
# Number of permutations: 9999
#
# adonis2(formula = distance ~ body.site, data = data.frame(sample_data(ps)), permutations = 9999, method = "bray")
# Df SumOfSqs R2 F Pr(>F)
# body.site 3 5.1913 0.39988 6.6633 1e-04 ***
# Residual 30 7.7909 0.60012
# Total 33 12.9823 1.00000
이렇게 input을 거리행렬로 주면, adonis2 함수 내에서 method를 바꾸어도 아무런 영향이 없는 것을 볼 수 있습니다.
3) PERMANOVA 주의점
adonis의 사용법을 설명하는 영상(https://www.youtube.com/watch?v=oLf0EpMJ4yA&t=1367s&ab_channel=RiffomonasProject ) 나와있는 것처럼, permanova는 일부 단점을 가지고 있습니다.
- 산포도가 불균일할 때 영향을 받는다 :
- PERMANOVA는 그룹 간 산포도가 다를 경우, 실제로는 그룹 간 중심값(centroid)이 같아도 차이가 있다고 잘못 판단할 수 있다. 이는 산포도 차이가 PERMANOVA 결과에 영향을 미쳐 잘못된 p-value 출력하는 것이다.
- 거리 척도기반의 한계:
- PERMANOVA는 거리 행렬을 기반으로 하기 때문에, 선택한 거리 척도(예: Bray-Curtis, Jaccard)에 따라 결과가 달라진다.
- 상호작용 효과의 해석 어려움:
- 여러 요인이 동시에 영향을 미칠 때, PERMANOVA는 이들 요인 간의 상호작용(term interaction)을 포함한 모델을 처리할 수 있지만, 이러한 상호작용 효과가 실험 결과에 어떤 의미를 가지는지 해석하기 어렵다.
- 순열 검정의 민감성:
- PERMANOVA의 p-value는 순열(permutation)의 수에 민감하다. 순열 수가 충분하지 않을 경우, p-value의 정확도가 떨어질 수 있으며, 이는 분석 결과의 신뢰성에 영향을 미칠 수 있다.
- 대안 부족:
- PERMANOVA는 여러 비모수적 다변량 분석 방법 중에서 산포도 문제에 덜 민감하지만, 여전히 산포도 차이에 완전히 자유로운 방법은 아니다. 따라서 추가적인 검정 도구인 betadisper를 사용하여 산포도의 동질성을 확인하는 것이 중요하다.
- 상관성을 입증하기 어렵다.
4) ANOSIM 과의 비교
Analysis of similarities (ANOSIM) 도 유사한 비모수 순열 검정법이다. 위 두 방법들은 같이 자주 사용되지만, 목적이 다르다.
- PERMANOVA: 거리행렬을 기반으로 한 분산 분석
- ANOSIM: 개체 간 분산이 아닌 유사성에 중점을 두어 분석
# ANOSIM
set.seed(42)
Ans <- anosim(braycurtis, data.frame(sample_data(ps))$body.site, permutations = 9999)
Ans
# Call:
# anosim(x = braycurtis, grouping = data.frame(sample_data(ps))$body.site, permutations = 9999)
# Dissimilarity: bray
#
# ANOSIM statistic R: 0.6856
# Significance: 1e-04
#
# Permutation: free
# Number of permutations: 9999
공통점
- 그룹 간 효과가 강하면 두 방법 모두 유사한 결과를 도출할 수 있다.
- 만약 PERMANOVA의 F-값과 ANOSIM의 R-값이 모두 낮으면, 이는 그룹 간 차이가 약함을 나타낸다.
- you usually obtain similar results if the effect of your factor is strong
- if both values(F-statistic of PERMANOVA and R value of ANOSIM) are low, it will indicate a weak effect of your factor
- Remember that PERMANOVA and ANOSIM test if your distribution is unique or not, not directly if you factor is significant or not
- The P value for PERMANOVA and ANOSIM is sensitive to the number of permutations you do (not the case for F and R values). So I would also advise you to try higher number of permutations (and the same between R and Primer) or all the permutations possible, and see how it affect your results.
+) Edinburgh Napier University의 Aimeric Blaud 曰 :
2 factor이상이 영향을 칠 때 PERMANOVA가 더 잘 처리한다
5. Betadisper?
베타 다양성 분석에서 흔히 발생하는 문제는 그룹 간 산포도 차이가 PERMANOVA 결과에 영향을 미칠 수 있다는 점이다. 따라서 betadisper는 그룹 간 산포도의 동질성을 검정하는 데 사용된다.
braycurtis <- phyloseq::distance(ps, method = "bray")
# 산포도 동질성 검정을 위한 betadisper 계산
betadisper_res <- betadisper(braycurtis, sample_data(ps)$body.site)
# Homogeneity of multivariate dispersions
#
# Call: betadisper(d = braycurtis, group = sample_data(ps)$body.site)
#
# No. of Positive Eigenvalues: 33
# No. of Negative Eigenvalues: 0
#
# Average distance to median:
# gut left palm right palm tongue
# 0.3963 0.5087 0.5918 0.3565
#
# Eigenvalues for PCoA axes:
# (Showing 8 of 33 eigenvalues)
# PCoA1 PCoA2 PCoA3 PCoA4 PCoA5 PCoA6 PCoA7 PCoA8
# 3.4040 2.3294 0.9766 0.9574 0.7273 0.6228 0.5628 0.4793
# betadisper 결과에 대해 ANOVA 수행
anova(betadisper_res)
# Analysis of Variance Table
#
# Response: Distances
# Df Sum Sq Mean Sq F value Pr(>F)
# Groups 3 0.30376 0.101254 13.076 1.231e-05 ***
# Residuals 30 0.23230 0.007743
# ---
# Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
betadisper의 분석은 각 그룹에서 중심값과 각 그룹에 해당하는 샘플 간 거리의 중앙값을 뜻한다. 여기서 여기서 right palm 그룹이 가장 큰 변동성을 보이는 반면, tongue 그룹이 가장 작은 변동성을 보인다.
PCoA 축은 각 데이터의 변동성을 설명하는 차원이다. 이때 PCoA1과 PCoA2가 가장 큰 변동성을 나타냄으로 시각화할 때, 두 축을 사용하면 충분하다는 것을 나타낸다.
ANOVA분석 결과 각 신체 부위별 산포도 값이 유의하게 다르다는 것을 의미한다. 이는 다시 말해서 PERMANOVA의 결과를 해석할 때 주의해야 한다는 것을 말한다.
즉, 각 그룹 내 차이(e.g. 혀 샘플 내 차이)가 그룹 간 차이(e.g. 혀와 손바닥 간의 차이) 보다 클 수 있다는 것을 의미한다.
그럼에도 불구하고 PERMANOVA의 F값과 ANOSIM의 R값이 모두 높다는 것은, bray-curtis 거리행렬에서 적어도 body site에 따른 영향이 크다는 것을 나타낸다.
6. 결론
마이크로바이옴 분석에서 PERMANOVA는 그룹 간 미생물 군집의 분산을 평가하는 다변량 분석 도구이다. 그러나 그룹에 따른 산포도에 따라 영향을 받을 수 있어, betadisper를 통해 확인하는 것이 중요하다. 또한 많이 사용되는 beta diversity 검정 도구인 ANOSIM은 순위 기반 유사성을 평가하는 반면, PERMANOVA는 분산을 비교한다.
분석 결과는 신중하게 해석해야 하며, p-value값이 높다고 해고 R 값과 F 값이 낮으면 전제 변동에 큰 영향을 미치지 않는 변수일 수 있다!
7. Reference
- 포럼 내 질문 들
https://www.researchgate.net/post/Why_do_I_obtain_different_results_using_PERMANOVA_or_ANOSIM
- Paliy, O., & Shankar, V. (2016). Application of multivariate statistical techniques in microbial ecology. Molecular Ecology, 25(5), 1032-1057. https://doi.org/10.1111/mec.13536
- Anderson, M. J., & I. Walsh, D. C. (2013). PERMANOVA, ANOSIM, and the Mantel test in the face of heterogeneous dispersions: What null hypothesis are you testing? Ecological Monographs, 83(4), 557-574. https://doi.org/10.1890/12-2010.1
- https://en.wikipedia.org/wiki/Permutational_analysis_of_variance
- 실행 코드
https://uw.pressbooks.pub/appliedmultivariatestatistics/chapter/permanova/
