Bioinformatics/└ 기타

PERMANOVA (Permutational analysis of variance) 란?

김해김씨99대손 2022. 5. 23. 10:06

수정: 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의 검정 수행하기

 

ps.rds
0.08MB

 

 

오늘도 어김없이 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는 일부 단점을 가지고 있습니다. 

 

  1. 산포도가 불균일할 때 영향을 받는다 :
    • PERMANOVA는 그룹 간 산포도가 다를 경우, 실제로는 그룹 간 중심값(centroid)이 같아도 차이가 있다고 잘못 판단할 수 있다. 이는 산포도 차이가 PERMANOVA 결과에 영향을 미쳐 잘못된 p-value 출력하는 것이다.
  2. 거리 척도기반의 한계:
    • PERMANOVA는 거리 행렬을 기반으로 하기 때문에, 선택한 거리 척도(예: Bray-Curtis, Jaccard)에 따라 결과가 달라진다. 
  3. 상호작용 효과의 해석 어려움:
    • 여러 요인이 동시에 영향을 미칠 때, PERMANOVA는 이들 요인 간의 상호작용(term interaction)을 포함한 모델을 처리할 수 있지만, 이러한 상호작용 효과가 실험 결과에 어떤 의미를 가지는지 해석하기 어렵다.
  4. 순열 검정의 민감성:
    • PERMANOVA의 p-value는 순열(permutation)의 수에 민감하다. 순열 수가 충분하지 않을 경우, p-value의 정확도가 떨어질 수 있으며, 이는 분석 결과의 신뢰성에 영향을 미칠 수 있다.
  5. 대안 부족:
    • PERMANOVA는 여러 비모수적 다변량 분석 방법 중에서 산포도 문제에 덜 민감하지만, 여전히 산포도 차이에 완전히 자유로운 방법은 아니다. 따라서 추가적인 검정 도구인 betadisper를 사용하여 산포도의 동질성을 확인하는 것이 중요하다.
  6. 상관성을 입증하기 어렵다. 

 

 

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

 

 

공통점

  1. 그룹 간 효과가 강하면 두 방법 모두 유사한 결과를 도출할 수 있다.
  2. 만약 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가 더 잘 처리한다

One last thing that might explain the difference: if the "models" you did for ANOSIM and PERMANOVA are not exactly the same. If you have more than 2 factors, the interactions between factors will influence the results depending how you make your two-way ANOSIM. PERMANOVA on the other hand will deal "better" with this type of issue. 

 

 

 

 

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

https://www.researchgate.net/post/Can_anyone_help_me_in_understanding_and_clearly_interpreting_ANOSIM_Analysis_of_Similarityand_SIMPER_Similarity_percentage_analysisresults

- Paliy, O., & Shankar, V. (2016). Application of multivariate statistical techniques in microbial ecology. Molecular Ecology25(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 Monographs83(4)557-574. https://doi.org/10.1890/12-2010.1

- https://archetypalecology.wordpress.com/2018/02/21/permutational-multivariate-analysis-of-variance-permanova-in-r-preliminary/

- https://en.wikipedia.org/wiki/Permutational_analysis_of_variance

- 실행 코드

  https://uw.pressbooks.pub/appliedmultivariatestatistics/chapter/permanova/

 

 

반응형