【原】跟著iMeta學做圖｜NMDS分析展示群落beta多樣性

宏基因組 2022-11-21 發(fā)布于北京

展開全文

原始教程鏈接：https://github.com/iMetaScience/iMetaPlot/tree/main/221108NMDS

寫在前面

非度量多維尺度分析(Non-metric multidimensional scaling, NMDS)，是基于相異矩陣或距離矩陣進行排序分析的間接梯度分析方法，在微生物組研究中可以用來展示群落beta多樣性。本期我們挑選2022年2月24日刊登在iMeta上的Linking soil fungi to bacterial community assembly in arid ecosystems - iMeta：西農(nóng)韋革宏團隊焦碩等-土壤真菌驅(qū)動細菌群落的構(gòu)建，選擇文章的Figure 6C進行復現(xiàn)，基于vegan包，講解和探討和NMDS分析和可視化的方法，先上原圖：

接下來，我們將通過詳盡的代碼逐步拆解原圖，最終實現(xiàn)對原圖的復現(xiàn)。代碼編寫及注釋：農(nóng)心生信工作室。

R包檢測和安裝

安裝核心R包vegan以及ggplot2，并載入所有R包。

if (!require("vegan"))  install.packages('vegan')if (!require("ggplot2"))  install.packages('ggplot2') # 加載包library(vegan)library(ggplot2)

生成測試數(shù)據(jù)

由于缺少原始數(shù)據(jù)，因此本例使用vegan包自帶的dune數(shù)據(jù)集進行測試。dune包含了20個樣品，每個樣品有30個物種豐度，每一行是一個樣品，每一列是一個物種。

# 載入dune數(shù)據(jù)集data(dune)#載入dune包含分組信息等的元數(shù)據(jù)（即metadata），分組信息為Management列data(dune.env)

NMDS分析

獲取數(shù)據(jù)后，即可利用vegan包進行NMDS分析。

#計算bray_curtis距離distance <- vegdist(dune, method = 'bray')#NMDS排序分析，k = 2預設兩個排序軸nmds <- metaMDS(distance, k = 2)#> Run 0 stress 0.1192678 #> Run 1 stress 0.1192678 #> ... Procrustes: rmse 4.495733e-05  max resid 0.0001375161 #> ... Similar to previous best#> Run 2 stress 0.1183186 #> ... New best solution#> ... Procrustes: rmse 0.02026799  max resid 0.06495211 #> Run 3 stress 0.1183186 #> ... New best solution#> ... Procrustes: rmse 1.832694e-05  max resid 5.57604e-05 #> ... Similar to previous best#> Run 4 stress 0.1809577 #> Run 5 stress 0.1192678 #> Run 6 stress 0.1183186 #> ... New best solution#> ... Procrustes: rmse 5.582524e-06  max resid 1.803473e-05 #> ... Similar to previous best#> Run 7 stress 0.1192678 #> Run 8 stress 0.1192678 #> Run 9 stress 0.1192678 #> Run 10 stress 0.1192678 #> Run 11 stress 0.1192679 #> Run 12 stress 0.1808911 #> Run 13 stress 0.1192678 #> Run 14 stress 0.1183186 #> ... Procrustes: rmse 5.943311e-06  max resid 1.823899e-05 #> ... Similar to previous best#> Run 15 stress 0.1886532 #> Run 16 stress 0.1192678 #> Run 17 stress 0.1183186 #> ... Procrustes: rmse 3.001088e-06  max resid 9.607646e-06 #> ... Similar to previous best#> Run 18 stress 0.1192679 #> Run 19 stress 0.1808911 #> Run 20 stress 0.1183186 #> ... Procrustes: rmse 2.027412e-05  max resid 6.520856e-05 #> ... Similar to previous best#> *** Best solution repeated 4 times#查看結(jié)果#summary(nmds)

獲取可視化所需數(shù)據(jù)。

#獲得應力值（stress）stress <- nmds$stress#將繪圖數(shù)據(jù)轉(zhuǎn)化為數(shù)據(jù)框df <- as.data.frame(nmds$points)#與分組數(shù)據(jù)合并df <- cbind(df, dune.env)

NMDS可視化

根據(jù)分組繪制一個最基礎的散點圖。

p <- ggplot(df, aes(MDS1, MDS2))+  geom_point(aes(color = Management), size = 5)

我們注意到，原圖中，每個分組被連接成不規(guī)則的多邊形并用不同顏色表示，我們可以通過ggplot2中geom_polygon()來繪制。geom_polygon()會按照數(shù)據(jù)中出現(xiàn)的順序連接觀測值，內(nèi)部可填充顏色。

p <- ggplot(df, aes(MDS1, MDS2))+  geom_point(aes(color = Management), size = 5)+  geom_polygon(aes(x = MDS1, y = MDS2, fill = Management, group = Management, color = Management),               alpha = 0.3, linetype = "longdash", linewidth = 1.5) #通過按順序連接觀測值繪制多邊形

由于geom_polygon()會按照數(shù)據(jù)中出現(xiàn)的順序連接觀測值，因此如果我們按照df自身順序來繪制多邊形，多邊形會非常奇怪，沒法代表不同分組。因此，我們需要預先處理df的順序，按合理的順序連接觀測值。

df <- df[order(df$Management), ]#先按分組排序df$Order <- c(2, 1, 3, 1, 2, 3, 4, 5, 3, 5, 1, 6, 2, 4, 1, 2, 6, 3, 5, 4)#添加一列Order，給每個分組內(nèi)觀測點的手動排序df <- df[order(df$Management, df$Order), ]#按分組和Order排序p <- ggplot(df, aes(MDS1, MDS2))+  geom_point(aes(color = Management), size = 5)+  geom_polygon(aes(x = MDS1, y = MDS2, fill = Management, group = Management, color = Management),               alpha = 0.3, linetype = "longdash", linewidth = 1.5)

分別進行Anosim分析（Analysis of similarities）和PERMANOVA（即adonis）檢驗分析。

#設置隨機種子set.seed(123)#基于bray-curtis距離進行PERMANOVA分析adonis <-  adonis2(dune ~ Management, data = dune.env, permutations = 999, method = "bray")#基于bray-curtis距離進行anosim分析anosim = anosim(dune, dune.env$Management, permutations = 999, distance = "bray")

美化圖片，并用AI微調(diào)。

# 應力值stress，Adonis R2與顯著性，Anosim R與顯著性stress_text <- paste("Stress  =", round(stress, 4))adonis_text <- paste(paste("Adonis  =", round(adonis$R2, 2)), "**")[1]anosim_text <- paste(paste("Anosim  =", round(anosim$statistic, 2)), "**")
p <- ggplot(df, aes(MDS1, MDS2))+  geom_point(aes(color = Management), size = 5)+  geom_polygon(aes(x = MDS1, y = MDS2, fill = Management, group = Management, color = Management), alpha = 0.3, linetype = "longdash", linewidth = 1.5)+  theme(plot.margin = unit(rep(1, 4), 'lines'),         panel.border = element_rect(fill = NA, color = "black", size = 0.5, linetype = "solid"),         panel.grid = element_blank(),         panel.background = element_rect(fill = 'white'))+  guides(color = "none", fill = "none")+  ggtitle(paste(paste(stress_text, adonis_text), anosim_text))

完整代碼

if (!require("vegan"))  install.packages('vegan')if (!require("ggplot2"))  install.packages('ggplot2') # 加載包library(vegan)library(ggplot2) # 載入dune數(shù)據(jù)集data(dune)#載入dune包含分組信息等的元數(shù)據(jù)（即metadata），分組信息為Management列data(dune.env)#計算bray_curtis距離distance <- vegdist(dune, method = 'bray')#NMDS排序分析，k = 2預設兩個排序軸nmds <- metaMDS(distance, k = 2)#> Run 0 stress 0.1192678 #> Run 1 stress 0.1192678 #> ... Procrustes: rmse 1.505128e-05  max resid 4.673581e-05 #> ... Similar to previous best#> Run 2 stress 0.1192678 #> ... Procrustes: rmse 3.715749e-06  max resid 1.009651e-05 #> ... Similar to previous best#> Run 3 stress 0.1889642 #> Run 4 stress 0.1192679 #> ... Procrustes: rmse 0.0001542849  max resid 0.0004702712 #> ... Similar to previous best#> Run 5 stress 0.1886532 #> Run 6 stress 0.2341212 #> Run 7 stress 0.1192678 #> ... Procrustes: rmse 1.328909e-05  max resid 4.273575e-05 #> ... Similar to previous best#> Run 8 stress 0.1886532 #> Run 9 stress 0.1192678 #> ... Procrustes: rmse 1.903819e-05  max resid 5.828243e-05 #> ... Similar to previous best#> Run 10 stress 0.1192678 #> ... Procrustes: rmse 6.358457e-06  max resid 1.687026e-05 #> ... Similar to previous best#> Run 11 stress 0.119268 #> ... Procrustes: rmse 5.501506e-05  max resid 0.0001605112 #> ... Similar to previous best#> Run 12 stress 0.1192678 #> ... New best solution#> ... Procrustes: rmse 5.074111e-06  max resid 1.393603e-05 #> ... Similar to previous best#> Run 13 stress 0.1192678 #> ... Procrustes: rmse 3.160318e-05  max resid 9.85043e-05 #> ... Similar to previous best#> Run 14 stress 0.1886532 #> Run 15 stress 0.2003486 #> Run 16 stress 0.2035424 #> Run 17 stress 0.1192678 #> ... Procrustes: rmse 2.440829e-05  max resid 7.079487e-05 #> ... Similar to previous best#> Run 18 stress 0.1183186 #> ... New best solution#> ... Procrustes: rmse 0.02027171  max resid 0.06497302 #> Run 19 stress 0.1183186 #> ... New best solution#> ... Procrustes: rmse 3.78469e-06  max resid 9.699447e-06 #> ... Similar to previous best#> Run 20 stress 0.1192678 #> *** Best solution repeated 1 times#查看結(jié)果#summary(nmds)#獲得應力值（stress）stress <- nmds$stress#將繪圖數(shù)據(jù)轉(zhuǎn)化為數(shù)據(jù)框df <- as.data.frame(nmds$points)#與分組數(shù)據(jù)合并df <- cbind(df, dune.env)df <- df[order(df$Management), ]#先按分組排序df$Order <- c(2, 1, 3, 1, 2, 3, 4, 5, 3, 5, 1, 6, 2, 4, 1, 2, 6, 3, 5, 4)#添加一列Order，給每個分組內(nèi)觀測點的手動排序df <- df[order(df$Management, df$Order), ]#按分組和Order排序 #設置隨機種子set.seed(123)#基于bray-curtis距離進行PERMANOVA分析adonis <- adonis2(dune ~ Management, data = dune.env, permutations = 999, method = "bray")#基于bray-curtis距離進行anosim分析anosim = anosim(dune, dune.env$Management, permutations = 999, distance = "bray") # 應力值stress，Adonis R2與顯著性，Anosim R與顯著性stress_text <- paste("Stress  =", round(stress, 4))adonis_text <- paste(paste("Adonis  =", round(adonis$R2, 2)), "**")[1]anosim_text <- paste(paste("Anosim  =", round(anosim$statistic, 2)), "**") p <- ggplot(df, aes(MDS1, MDS2))+  geom_point(aes(color = Management), size = 5)+  geom_polygon(aes(x = MDS1, y = MDS2, fill = Management, group = Management, color = Management),               alpha = 0.3, linetype = "longdash", linewidth = 1.5)+  theme(plot.margin = unit(rep(1, 4), 'lines'),         panel.border = element_rect(fill = NA, color = "black", size = 0.5, linetype = "solid"),         panel.grid = element_blank(),         panel.background = element_rect(fill = 'white'))+  guides(color = "none", fill = "none")+  ggtitle(paste(paste(stress_text, adonis_text), anosim_text))
ggsave("Figure6C.pdf", p, height = 5.69, width = 7.42)