R包vegan的Mantel tests探索群落物種組成是否與環境相關

Mantel tests是確定兩組距離測度矩陣（而非兩組變量矩陣）之間相關性的相關性測試方法，用於判斷一個矩陣中的樣本距離與另一矩陣中的樣本距離是否相關。Mantel tests零假設為響應變量矩陣中對象之間的距離與解釋變量矩陣不存在相關，如果結果中p值顯著，則拒絕零假設，即存在相關性，隨著一個矩陣中樣本之間距離的增加（或減少），另一矩陣中對應樣本之間的距離也增加（或減少）。此外，Mantel方法還可用於檢驗假設或模型。在這種模型測試方法中，一個矩陣包含響應數據，另一個矩陣代表了要測試的先驗模型（檢驗的備擇假設），如果找到了重要的Mantel統計信息，它們將為模型提供一些支持。本篇重點介紹Mantel tests確定相關性的方法。例如在群落分析中，為了探索群落物種組成是否與環境相關，經常使用到Mantel tests。群落分析中通常存在兩組變量矩陣，樣方-物種多度矩陣和樣方-環境測量矩陣。首先根據兩組變量矩陣計算樣方間的相異（距離）矩陣，即分別獲得通過物種豐度計算的樣方距離（通常為Bray-Curtis距離），以及通過某幾種環境參數計算的樣方距離（通常為歐幾裡得距離）。有時也會使用樣方間真實的地理距離直接作為某種距離測度。

之後使用兩組距離測度矩陣執行Mantel tests，例如確定樣方之間的群落組成差異是否與樣方之間的溫度差異或樣方之間的物理距離相關，或者說「共變」。這些測試可用於解決環境是針對微生物群落的「選擇」，還是存在強烈的距離衰減模式，表明存在擴散限制。這些通常是生物地理學研究中的重要問題。

本篇同樣以群落分析為例，簡介R包vegan的Mantel tests。

假設存在如下數據集。第1列是樣方名稱，第2-5列為各樣方中的環境參數（即鹽度、溫度等），第6-7列為各樣方的緯度和經度，第8列及之後為各樣方中的物種及其豐度。我們期望通過Mantel tests，查看對於該數據集，作用於物種變化的最主要因素，是由環境引起的「選擇」，還是由地理因素的擴散限制所致。

 

加載R包，如上所述，首先計算兩組樣方距離測度，然後執行Mantel tests。

library(vegan)
 
#讀取上述數據集
df <- read.csv('Your_OTU_table.csv', header= TRUE)
 
##計算距離
#根據物種豐度數據，計算樣方間的 Bray-curtis 距離
abund <- df[ ,8:ncol(df)]
dist.abund <- vegdist(abund, method = 'bray')
 
#根據環境測量指標，計算樣方間的歐幾裡得距離
#這裡只選擇了其中的溫度指標，期望關注物種變化與溫度的相關性
temp <- df$Temperature
dist.temp <- dist(temp, method = 'euclidean')
 
#如果期望關注多種環境的協同作用，就選擇一個環境子集，計算樣方間的歐幾裡得距離
#例如使用 4 種環境數據，但此時需要執行數據標準化，以消除量綱差異
env <- df[ ,2:5]
scale.env <- scale(env, center = TRUE, scale = TRUE)
dist.env <- dist(scale.env, method = 'euclidean')
 
#根據經緯度，計算樣方間實際的地理距離
geo <- data.frame(df$Longitude, df$Latitude)
d.geo <- distm(geo, fun = distHaversine)       #library(geosphere)
dist.geo <- as.dist(d.geo)
 
##執行 Mantel tests，詳情 ?mantel，以下為 3 個示例
#物種豐度和溫度的相關性，以 spearman 相關係數為例，9999 次置換檢驗顯著性（Mantel tests 基於隨機置換的方法獲取 p 值）
abund_temp <- mantel(dist.abund, dist.temp, method = 'spearman', permutations = 9999, na.rm = TRUE)
abund_temp
 
#物種豐度和地理距離的相關性，以 spearman 相關係數為例，9999 次置換檢驗顯著性
abund_geo <- mantel(dist.abund, dist.geo, method = 'spearman', permutations = 9999, na.rm = TRUE)
abund_geo
 
#物種豐度和 4 種環境組合的相關性，以 spearman 相關係數為例，9999 次置換檢驗顯著性
abund_env <- mantel(dist.abund, dist.env, method = 'spearman', permutations = 9999, na.rm = TRUE)
abund_env

基於物種豐度的距離矩陣與基於溫度指標的距離矩陣之間有很強的相關性（Mantel statistic R: 0.667，p value = 1e-04）。換句話說，隨著樣方在溫度方面的差異逐漸增大，它們在物種組成方面的差異也越來越大。

#物種豐度和溫度的相關性
> abund_temp

Mantel statistic based on Spearman's rank correlation rho 

Call:
mantel(xdis = dist.abund, ydis = dist.temp, method = "spearman", permutations = 9999, na.rm = TRUE) 

Mantel statistic r: 0.677 
      Significance: 1e-04 

Upper quantiles of permutations (null model):
  90%   95% 97.5%   99% 
0.148 0.198 0.246 0.290 
Permutation: free
Number of permutations: 9999

基於物種豐度的距離矩陣與樣方間的地理距離沒有顯著關係（Mantel statistic R: 0.138，p value = 0.052）。因此可知，對於該測試數據集，不存在物種豐度的距離衰減效應。

#物種豐度和地理距離的相關性
> abund_geo

Mantel statistic based on Spearman's rank correlation rho 

Call:
mantel(xdis = dist.abund, ydis = dist.geo, method = "spearman", permutations = 9999, na.rm = TRUE) 

Mantel statistic r: 0.1379 
      Significance: 0.0525 

Upper quantiles of permutations (null model):
  90%   95% 97.5%   99% 
0.107 0.140 0.170 0.204 
Permutation: free
Number of permutations: 9999

同時對於4種環境變量組合，累積的環境因素與群落物種組成高度相關（Mantel statistic r: 0.686, p value = 1e-04）。

#物種豐度和 4 種環境組合的相關性
> abund_env

Call:
mantel(xdis = dist.abund, ydis = dist.env, method = "spearman",      permutations = 9999, na.rm = TRUE) 

Mantel statistic r: 0.6858 
      Significance: 1e-04 

Upper quantiles of permutations (null model):
  90%   95% 97.5%   99% 
0.151 0.201 0.244 0.292 
Permutation: free
Number of permutations: 9999

綜上結論，對於該數據集，與地理距離相比，群落物種組成與環境參數的相關性更強。因此在該系統中，主要發生環境對群落作出的「選擇」，地理因素的擴散限制相對微弱。

最後不妨作圖觀測變量間的關係，加深對這種相關性的理解。

library(ggplot2)
 
#某物種與溫度的相關性，橫軸溫度，縱軸物種豐度，顏色表示樣方的緯度
xx = ggplot(df, aes(x = Temperature, y = Pelagibacteraceae.OTU_307744)) + 
    geom_smooth(method = 'lm', alpha = 0.2, colour = 'black') + 
    geom_point(aes(colour = Latitude), size = 4) +
    labs(y = 'Pelagibacteraceae (OTU 307744) (%)', x = 'Temperature (C)') + 
    theme( axis.text.x = element_text(face = 'bold',colour = 'black', size = 12), 
        axis.text.y = element_text(face = 'bold', size = 11, colour = 'black'), 
        axis.title= element_text(face = 'bold', size = 14, colour = 'black'), 
        panel.background = element_blank(), 
        panel.border = element_rect(fill = NA, colour = 'black'), 
        legend.title = element_text(size =12, face = 'bold', colour = 'black'),
        legend.text = element_text(size = 10, face = 'bold', colour = 'black')) +
    scale_colour_continuous(high = 'navy', low = 'salmon')
 
xx

#對於圖中的線性回歸
fit <- lm(df$Temperature~df$Pelagibacteraceae.OTU_307744)
summary(fit)

Call:
lm(formula = df$Temperature ~ df$Pelagibacteraceae.OTU_307744)

Residuals:
    Min      1Q  Median      3Q     Max 
-2.2053 -0.9336 -0.5215  0.5028  3.8232 

Coefficients:
                                Estimate Std. Error t value Pr(>|t|)    
(Intercept)                       0.4082     0.4476   0.912    0.372    
df$Pelagibacteraceae.OTU_307744   1.3008     0.1280  10.165 1.45e-09 ***
---
Signif. codes:  0 『***』 0.001 『**』 0.01 『*』 0.05 『.』 0.1 『 』 1

Residual standard error: 1.634 on 21 degrees of freedom
Multiple R-squared:  0.8311, Adjusted R-squared:  0.823 
F-statistic: 103.3 on 1 and 21 DF,  p-value: 1.454e-09

#分面圖展示多組變量的相關性，橫軸溫度，縱軸為多個物種的豐度，顏色表示樣方的緯度
library(reshape2)
 
otus <- df[ ,1:11]
otus_melt <- melt(otus, id = c('Station', 'Salinity', 'Temperature', 'Oxygen', 'Nitrate', 'Latitude', 'Longitude'))
 
xx <- ggplot(otus_melt, aes(x = Temperature, y = value)) + 
    facet_wrap(.~variable, scales = 'free_y') +
    geom_smooth(method = 'lm', alpha = 0.2, colour = 'black') + 
    geom_point(aes(colour = Latitude), size = 4) +
    labs(y = 'Relative Abundance (%)', x = 'Temperature (C)') + 
    theme( axis.text.x = element_text(face = 'bold',colour = 'black', size = 12), 
        axis.text.y = element_text(face = 'bold', size = 10, colour = 'black'), 
        axis.title= element_text(face = 'bold', size = 14, colour = 'black'), 
        panel.background = element_blank(), 
        panel.border = element_rect(fill = NA, colour = 'black'), 
        legend.title = element_text(size =12, face = 'bold', colour = 'black'),
        legend.text = element_text(size = 10, face = 'bold', colour = 'black'), 
        legend.position = 'top', strip.background = element_rect(fill = 'grey90', colour = 'black'),
        strip.text = element_text(size = 9, face = 'bold')) +
    scale_colour_continuous(high = 'navy', low = 'salmon')
 
xx

 

上述主要展示的變量間相關性的散點圖。

接下來是對於距離測度間的相關性。

#將上文獲得的距離測度，轉化為數據框，一一對應起來
aa <- as.vector(dist.abund)
tt <- as.vector(dist.temp)
gg <- as.vector(dist.geo)
mat <- data.frame(aa, tt, gg)
 
#基於物種豐度的距離與基於溫度指標的距離之間的相關性散點圖，上文已知二者顯著相關；同時顏色表示樣方間地理距離
mm <- ggplot(mat, aes(y = aa, x = tt)) + 
    geom_point(size = 4, alpha = 0.75, colour = "black",shape = 21, aes(fill = gg/1000)) + 
    geom_smooth(method = "lm", colour = "black", alpha = 0.2) + 
    labs(x = "Difference in Temperature (C)", y = "Bray-Curtis Dissimilarity", fill = "Physical Separation (km)") + 
    theme( axis.text.x = element_text(face = "bold",colour = "black", size = 12), 
           axis.text.y = element_text(face = "bold", size = 11, colour = "black"), 
           axis.title= element_text(face = "bold", size = 14, colour = "black"), 
           panel.background = element_blank(), 
           panel.border = element_rect(fill = NA, colour = "black"),
           legend.position = "top",
           legend.text = element_text(size = 10, face = "bold"),
           legend.title = element_text(size = 11, face = "bold")) +
    scale_fill_continuous(high = "navy", low = "skyblue")
 
mm
 
#基於物種豐度的距離與樣方間地理距離之間的相關性散點圖，上文已知二者無相關性
mm <- ggplot(mat, aes(y = aa, x = gg/1000)) + 
    geom_point(size = 3, alpha = 0.5) + 
    labs(x = "Physical separation (km)", y = "Bray-Curtis Dissimilarity") + 
    theme( axis.text.x = element_text(face = "bold",colour = "black", size = 12), 
        axis.text.y = element_text(face = "bold", size = 11, colour = "black"), 
        axis.title= element_text(face = "bold", size = 14, colour = "black"), 
        panel.background = element_blank(), 
        panel.border = element_rect(fill = NA, colour = "black"))
 
mm

 

Mantel Test in R：https://jkzorz.github.io/2019/07/08/mantel-test.htmlAppendix 4: Graphical Description of MantelsTest and ANOSIM：https://www.mfe.govt.nz/publications/environmental-reporting/new-zealand-marine-environment-classification-overview/append-2The Mantel test：https://mb3is.megx.net/gustame/hypothesis-tests/the-mantel-test

  

為鼓勵讀者交流、快速解決科研困難，我們建立了「宏基因組」專業討論群，目前己有國內外5000+ 一線科研人員加入。參與討論，獲得專業解答，歡迎分享此文至朋友圈，並掃碼加主編好友帶你入群，務必備註「姓名-單位-研究方向-職稱/年級」。PI請明示身份，另有海內外微生物相關PI群供大佬合作交流。技術問題尋求幫助，首先閱讀《如何優雅的提問》學習解決問題思路，仍未解決群內討論，問題不私聊，幫助同行。