- 수정 : 2023-04-11
- 원문: https://stackoverflow.com/questions/14711470/plotting-envfit-vectors-vegan-package-in-ggplot2-
- 원문을 번역하여 풀어쓴 글임
오늘 분석에 사용할 데이터는 "dune"으로 네덜란드 모래언덕 초지의 식물 데이터이다. dune은 20 개의 지역에서 관찰한 30종의 식물의 개수가 적혀있는 표형태의 데이터이다. 이 데이터를 이용해 vegan 패키지의 envif를 적용하여 어떤 지역에 어떤 종이 가장 많이 자라는지를 확인해 볼 것이다.
# 데이터 불러오기
library(vegan)
data(dune)
head(dune)
# Achimill Agrostol Airaprae Alopgeni Anthodor Bellpere Bromhord Chenalbu Cirsarve Comapalu Eleopalu Elymrepe Empenigr Hyporadi Juncarti Juncbufo
# 1 1 0 0 0 0 0 0 0 0 0 0 4 0 0 0 0
# 2 3 0 0 2 0 3 4 0 0 0 0 4 0 0 0 0
# 3 0 4 0 7 0 2 0 0 0 0 0 4 0 0 0 0
# 4 0 8 0 2 0 2 3 0 2 0 0 4 0 0 0 0
# 5 2 0 0 0 4 2 2 0 0 0 0 4 0 0 0 0이 예제의 결과 그림은 아래와 같다
legend를 보면 각 샘플의 채취된 그룹을 뜻하며, 채취 장소는 총 20곳이다.
그림의 각 점은 sample을 채취한 site 20개를 뜻하고, 화살표에 적힌 이름은 초지에 자생하는 풀의 종명이다.
dune데이터를 NMDS로 분석하여, 두 개의 축(mds1, mds2)으로 나누어 데이터(각 샘플링한 지역) 간의 거리를 그림으로 나타낸 것이다. 또한 envfit(그림 속 화살표)을 이용해 site와 각 생물 종의 상관관계 또한 분석하였다.
같은 데이터를 아래와 같이 그릴 수도 있다.
두 그림을 종합해서 보면, 왼쪽 하단에 위치한 Poatri와 Elyrep는 hi그룹에 속하는 지역에 서식하는 것을 알 수 있으며.
오른쪽 상단의 Junart, Potpal 같은 식물은 md그룹에 속하는 지역에 서식하는 것을 알 수 있다.
| NMDS계산
이제 첫 번째 그림을 그려보자.
library(ggplot2)
library(grid)
# calculate distance for NMDS
NMDS.log<-log(dune+1)
sol <- metaMDS(NMDS.log)log단위로 바꿀 때 log(1)은 0 값을 가짐으로 1을 더한 뒤 log를 취해준다. 그러다 현재 분석하는 데이터인 dune 은 +1의 유무로 결과가 크게 달라지지 않는다.
## log 취하지 않고 NMDS계산 ##
sol <- metaMDS(dune)
head(sol$points)
# MDS1 MDS2
# 1 -0.84056878 -0.71580311
# 2 -0.50485241 -0.40893588
# 3 -0.08267594 -0.43666992
# 4 -0.11563020 -0.52223904
# 5 -0.62654209 -0.08668995
## log처리 한 후 NMDS계산 ##
NMDS.log<-log(dune+1)
sol <- metaMDS(NMDS.log)
head(sol$points)
# MDS1 MDS2
# 1 -0.88388170 -0.7042364543
# 2 -0.49706451 -0.4173348657
# 3 -0.06826958 -0.3994499673
# 4 -0.12820511 -0.4699872619
# 5 -0.58377203 -0.0494215159보다시피 큰 차이는 없다. 마지막에 그려지는 plot에서도 큰 차이는 없었다.
그러나 환경 데이터가 아닌 미생물 데이터 분석 시 count-table이 아니라 OTU table을 사용한다. 이때 표 안에 카운팅 되는 read의 수는 0에서 10,000의 값을 가질 수 있다. 이렇게 수의 차이가 클 때는 +1을 해준 뒤에 log상태로 계산해 준다.
계산된 NMDS의 평면 값 중 두 평면만 가져온다.
NMDS = data.frame(NMDS1 = sol$points[,1], NMDS2 = sol$points[,2])
head(NMDS)
# NMDS1 NMDS2
# 1 -0.88388118 -0.7042371356
# 2 -0.49706606 -0.4173326048
# 3 -0.06827147 -0.3994498770
# 4 -0.12820729 -0.4699866952
# 5 -0.58377296 -0.0494198342
각 값은 종의 풍부도를 기준으로 20여 개의 지역이 멀고 가까움을 거리로 표시한 것이다.
| envfit() 함수 사용
envfit이란 sol$point 가 각 지역의 멀고 가까운 관계를 나타냈다고 했을 때, 이들과 각 종 간의 관계를 또 한 번 계산하는 것이다.
vec.sp<-envfit(sol$points, dune, perm=1000)
vec.sp.df<-as.data.frame(vec.sp$vectors$arrows*sqrt(vec.sp$vectors$r))
vec.sp.df$species<-rownames(vec.sp.df)
head(vec.sp.df)
# MDS1 MDS2 species
# Achimill -0.64489339 0.008495629 Achimill
# Agrostol 0.71095363 -0.453260640 Agrostol
# Airaprae -0.16361717 0.690373490 Airaprae
# Alopgeni 0.25893419 -0.580802109 Alopgeni
# Anthodor -0.47660281 0.561990600 Anthodor
- envif함수 옵션 중 perm은 999나 1000을 쓰는데 많으면 소프트웨어에 부담 가지만 오류를 줄일 수 있다.
- 화살표의 방향은 list vector와 matrix arrow로 저장된다. 좌표를 얻기 위해선 vectors$r 에 있는 좌표를 제곱한 후에 루트씨워서 음수로 바꾸어 주자 → vec.sp.df<-as.data.frame(vec.sp$vectors$arrows*sqrt(vec.sp$vectors$r))
- 아니면 vec.sp.df <- as.data.frame(scores(vec.sp, display = "vectors")) 로 적어도 무방하다
# Create meta data for grouping
MyMeta = cbind(NMDS, Group = c("hi","hi","hi","md","lo","hi","hi","lo","md","md",
"lo","lo","hi","lo","hi","md","md","lo","hi","lo"))
head(MyMeta)
# NMDS1 NMDS2 Group
# 1 -0.84053172 -0.71584750 hi
# 2 -0.50485056 -0.40894127 hi
# 3 -0.08266363 -0.43667949 hi
# 4 -0.11561734 -0.52224112 md
# 5 -0.62654283 -0.08670242 lo
20개의 지역을 구분하기 위해 group데이터를 추가한다.
| 시각화 (ggplot)
ggplot(MyMeta) +
geom_point(mapping = aes(x = NMDS1, y = NMDS2, colour = Group))+ # NMDS for group
coord_fixed()
각 지역 간의 관계를 나타내는 nmds plot이 그려졌으며, 색을 통해 각 지역이 속하는 그룹을 구분하였다.
arrow를 그리기 위해선 일단 geom_segment () 함수를 이용한다.
ggplot(MyMeta) +
geom_point(mapping = aes(x = NMDS1, y = NMDS2, colour = Group))+ # NMDS for group
coord_fixed()+
geom_segment(data = vec.sp.df, # only line
aes(x = 0, xend = MDS1, y = 0, yend = MDS2) , #node=0 / length = NMDS /
arrow = arrow(length = unit(0.25, "cm")), # 화살표
colour = "grey")+
geom_text(data = vec.sp.df, aes(x = MDS1, y = MDS2, label = species),size = 3)이때, 0점에서 mdns값까지 뻗어나가는 형태의 화살표를 표시하기 위해 옵션은 aes(x = 0, xend = MDS1, y = 0, yend = MDS2)로 적어준다
처음 예제 그림과 같은 그림이 완성되었다.
| 깔끔하게 그리기
# 배경과 겹친 이름 삭제 하기
ggplot(MyMeta) +
geom_point(mapping = aes(x = NMDS1, y = NMDS2, colour = Group))+ # NMDS for group
coord_fixed()+
#geom_point(data = vec.sp.df, mapping = aes(x = MDS1, y = MDS2)) +
geom_segment(data = vec.sp.df, # only line
aes(x = 0, xend = MDS1, y = 0, yend = MDS2) , #node=0 / length = NMDS /
arrow = arrow(length = unit(0.25, "cm")), # 화살표
colour = "grey")+
geom_text(data = vec.sp.df, aes(x = MDS1, y = MDS2, label = species),size = 3, check_overlap=TRUE)+
theme_bw()+
theme(legend.position="none")
- ggrepel package를 참고하여 그린 그림
(ggrepel : https://mran.microsoft.com/snapshot/2017-08-20/web/packages/ggrepel/vignettes/ggrepel.html)
library('ggrepel')
ggplot(MyMeta) +
geom_point(mapping = aes(x = NMDS1, y = NMDS2, colour = Group))+ # NMDS for group
coord_fixed()+
#geom_point(data = vec.sp.df, mapping = aes(x = MDS1, y = MDS2)) +
geom_segment(data = vec.sp.df, # only line
aes(x = 0, xend = MDS1, y = 0, yend = MDS2) , #node=0 / length = NMDS /
arrow = arrow(length = unit(0.25, "cm")), # 화살표
colour = "grey")+
#geom_text_repel(data = vec.sp.df, aes(x = MDS1, y = MDS2, label = species),size = 3)+
theme_bw()+
theme(legend.position="none")+
geom_label_repel(data = vec.sp.df, aes(x = MDS1, y = MDS2, label = species, fill = factor(species)),size = 3,
fontface = 'bold', color = 'black',
box.padding = unit(0.35, "lines"),
point.padding = unit(0.5, "lines"))
geom_label_repel()을 쓰지 않고 geom_text_repel(data = vec.sp.df, aes(x = MDS1, y = MDS2, label = species), size = 3)를 사용하면 그림이 조금(?) 이상하다
화살표 초점이 안 맞는 느낌.. 이 난다
| 더보기
[Vegan] phyloseq 개체를 이용해서 biplot그리기 (2023.09.05)
| 출처
- Dune 데이터 : Jongman, R.H.G, ter Braak, C.J.F & van Tongeren, O.F.R. (1987). Data Analysis in Community and Landscape Ecology. Pudoc, Wageningen.
- 시각화 : https://stackoverflow.com/questions/14711470/plotting-envfit-vectors-vegan-package-in-ggplot2
- 같은 데이터를 사용한 비슷한 예제 : https://stackoverflow.com/questions/13794419/plotting-ordiellipse-function-from-vegan-package-onto-nmds-plot-created-in-ggplo
