?? 利用R語言繪制物種組成圖。本文以堆疊面積圖的方式與大家分享。

面積圖簡介

面積圖又叫區(qū)域圖。它是在折線圖的基礎之上形成的, 它將折線圖中折線與自變量坐標軸之間的區(qū)域使用顏色或者紋理填充，這樣一個填充區(qū)域我們叫面積。顏色的填充可以更好地突出趨勢信息（比如時間上的差異，分類上的差異），需要注意的是顏色要帶有一定的透明度，透明度可以很好地幫助使用者觀察不同序列之間的重疊關系，沒有透明度的面積會導致不同序列之間相互遮蓋減少可以被觀察到的信息。

與折線圖相似，面積圖可用于強調(diào)數(shù)量隨時間或分類而變化的程度，也可用于引起人們對總值趨勢的注意。他們最常用于表現(xiàn)趨勢和關系，而不是傳達特定的值。

• 所有的數(shù)據(jù)都從相同的零軸開始。
• 標準面積圖適用于展示或者比較隨著時間連續(xù)變化的定量。
• 在需要繪制大量數(shù)據(jù)系列的情況下，折線圖通常是更清晰的可視化表達方式。

準備數(shù)據(jù)集

加載數(shù)據(jù)集

otudf <- read.csv("otudf.csv", header = T, row.names = 1)
group_infor <- read.csv("group_infor.csv", header = T, row.names = 1)head(otudf[1:6, 1:6])
# CK_L_1 CK_L_2 CK_L_3 A_L_1 A_L_2 A_L_3
# otu_1    414   1371    113  1190   427  1327
# otu_2    462   1028    937   394  1266   547
# otu_3    178    193    722  1421   705   579
# otu_4    525     56    822   885   690   944
# otu_5    194   1353    365   680  1218   677
# otu_6    937   1140   1499  1098  1460   494
head(group_infor)
# Repeat Treat Depth Samples
# 1      1    CK     L  CK_L_1
# 2      2    CK     L  CK_L_2
# 3      3    CK     L  CK_L_3
# 4      1     A     L   A_L_1
# 5      2     A     L   A_L_2
# 6      3     A     L   A_L_3

數(shù)據(jù)處理

## 提取物種豐度和物種注釋信息
species <- otudf[c(1:36)]
taxa <- otudf[-c(1:36)]## 確定數(shù)據(jù)正確與否
head(species)
# CK_L_1 CK_L_2 CK_L_3 A_L_1 A_L_2 A_L_3 B_L_1 B_L_2 B_L_3 C_L_1 C_L_2 C_L_3 CK_M_1 CK_M_2 CK_M_3 A_M_1 A_M_2 A_M_3
# otu_1 414 1371 113 1190 427 1327 631 207 1473 1322 1370 1476 28 1271 125 978 945 0
# otu_2 462 1028 937 394 1266 547 1144 851 365 1418 651 1329 1061 1124 332 919 1483 600
# otu_3 178 193 722 1421 705 579 1103 1398 241 1328 1286 452 263 1446 893 1098 868 661
# otu_4 525 56 822 885 690 944 86 720 1135 37 1323 792 362 696 806 1398 11 37
# otu_5 194 1353 365 680 1218 677 1129 547 829 1038 103 276 1166 630 539 1298 277 663
# otu_6 937 1140 1499 1098 1460 494 988 709 1137 251 356 130 344 727 1332 516 1154 454
# B_M_1 B_M_2 B_M_3 C_M_1 C_M_2 C_M_3 CK_H_1 CK_H_2 CK_H_3 A_H_1 A_H_2 A_H_3 B_H_1 B_H_2 B_H_3 C_H_1 C_H_2 C_H_3
# otu_1 1394 213 1449 1431 786 697 627 463 935 58 118 967 864 226 1053 834 551 961
# otu_2 985 247 1135 943 584 264 1286 838 37 1057 320 982 1351 621 1058 667 696 780
# otu_3 1117 566 508 142 1016 924 244 1426 1093 994 1166 962 632 628 163 1455 756 878
# otu_4 244 787 1261 533 1025 583 201 533 1051 651 1294 20 206 1151 369 1499 358 1083
# otu_5 716 736 1219 1411 1295 3 434 1023 406 1215 1074 994 1220 1220 528 1349 1138 906
# otu_6 1053 163 295 412 473 317 1065 981 1036 1113 435 19 1396 1151 1048 134 1199 383
head(taxa)
# phylum class order family genus species
# otu_1 Spirochaetes Gemmatimonadetes BRC1 Firmicutes Planctomycetes Nitrospirae
# otu_2 WPS_1 Euryarchaeota Cyanobacteria Armatimonadetes Spirochaetes Chloroflexi
# otu_3 Gemmatimonadetes Gemmatimonadetes Armatimonadetes Cyanobacteria WPS_2 WPS_1
# otu_4 WPS_1 Chloroflexi Acidobacteria Firmicutes Acidobacteria WPS_2
# otu_5 Acidobacteria Cyanobacteria BRC1 Chloroflexi Bacteroidetes Spirochaetes
# otu_6 Planctomycetes Gemmatimonadetes WPS_1 Chloroflexi Chloroflexi Armatimonadetes# 計算不同處理下的門水平均值
library(tidyverse)
Phylum <- species %>%group_by(taxa$phylum) %>% # 使用taxa中的門水平進行分類summarise_all(sum) %>% # 計算總的otu數(shù)量rename(Phylum = `taxa$phylum`) %>% # 修改名稱gather(key = "Samples", value = "Abundance", -Phylum) %>% # 數(shù)據(jù)形式轉(zhuǎn)換 寬-長left_join(group_infor, by = c("Samples" = "Samples")) %>%select(Phylum, Treat, Depth, Abundance) %>%group_by(Phylum, Treat, Depth) %>% # 求均值summarise_all(mean) %>%arrange(Phylum, Treat, Depth, desc(Abundance))
Phylum
# # A tibble: 204 × 4
# # Groups: Phylum, Treat [68]
# Phylum Treat Depth Abundance
# <chr> <chr> <chr> <dbl>
# 1 Acidobacteria A H 462665.
# 2 Acidobacteria A L 459384.
# 3 Acidobacteria A M 454049.
# 4 Acidobacteria B H 454310
# 5 Acidobacteria B L 455244.
# 6 Acidobacteria B M 448743.
# 7 Acidobacteria C H 449425
# 8 Acidobacteria C L 455059.
# 9 Acidobacteria C M 454707.
# 10 Acidobacteria CK H 456931
# # … with 194 more rows
# # ? Use `print(n = ...)` to see more rows## 選擇top9 及 合并剩余物種作為Other
phy_select <- unique(Phylum$Phylum)[1:9]top9 <- Phylum[Phylum$Phylum %in% phy_select, ]
other <- Phylum[!Phylum$Phylum %in% phy_select, ]
other <- other %>% group_by(Treat, Depth) %>% summarise(Abundance = mean(Abundance)) %>% cbind(Phylum = "Other") %>%select(Phylum, Treat, Depth, Abundance, everything())top10 <- rbind(top9, other)
head(top10)
# # A tibble: 6 × 4
# # Groups: Phylum, Treat [2]
# Phylum Treat Depth Abundance
# <chr> <chr> <chr> <dbl>
# 1 Acidobacteria A H 462665.
# 2 Acidobacteria A L 459384.
# 3 Acidobacteria A M 454049.
# 4 Acidobacteria B H 454310
# 5 Acidobacteria B L 455244.
# 6 Acidobacteria B M 448743.

數(shù)據(jù)可視化

# 物種組成堆疊面積圖
library(ggplot2)
library(ggalluvial)
ggplot(data = top10,aes(x = Depth, y = Abundance, fill = reorder(Phylum, -Abundance),colour = reorder(Phylum, -Abundance),stratum = reorder(Phylum, -Abundance) ,alluvium = reorder(Phylum, -Abundance))) +geom_alluvium(aes(fill = reorder(Phylum, -Abundance)), alpha = 0.7, decreasing = FALSE) +geom_stratum(aes(fill = reorder(Phylum, Abundance)), width = 0.3, size = 0.1, color = "black") +scale_y_continuous(expand = c(0, 0)) +theme_bw() +facet_grid(. ~ Treat, scales = "fixed") +scale_fill_manual(values = c("#EB7369", "#CF8B0B", "#9D9F20", "#2BB077", "#2BB077","#1BB3B7", "#29A4DE", "#8989C1", "#B174AD","#DE66A1"), name =  "Phylum") +scale_color_manual(values = c("#EB7369", "#CF8B0B", "#9D9F20", "#2BB077", "#2BB077","#1BB3B7", "#29A4DE", "#8989C1", "#B174AD","#DE66A1")) +guides(color = "none")+theme(panel.grid=element_blank(),panel.spacing.x = unit(0, units = "cm"),strip.background = element_rect(color = "white", fill = "white", linetype = "solid", size = 1),strip.placement = "outside",axis.line.y.left = element_line(color = "black", size = 0.7),axis.line.x.bottom = element_line(color = "black", size = 0.7),strip.text.x = element_text(size = 14, face = "bold"),axis.text = element_text(face = "bold", size = 12, color = "black"),axis.title = element_text(face = "bold", size = 14, colour = "black"),legend.title = element_text(face = "bold", size = 12, color = "black"),legend.text = element_text(face = "bold", size = 12, color = "black"),axis.ticks.x = element_line(size = 1),axis.ticks.y = element_line(size = 1),)+labs(x = "Depth",y= "Relative Abundance of Phylum (%)")