寫在前面

之前寫過一篇關于統計學軟件STAMP的教程，使用STAMP對微生物群落數據進行統計學分析還是挺方便的，尤其是對R語言不是很熟悉的朋友來說，圖形化的界面和相對簡單的操作還是挺友好的。

• STAMP——微生物組間差異統計分析 簡明教程 中文幫助文檔

我們通常使用的STAMP的結果主要就是兩組數據之間差異性檢驗的被稱作Extended error bar(擴展柱狀圖)的圖像。

由于STAMP的結果圖相對固定，可修改的圖像參數有限，經常會遇到一些問題，比如靶標物種或功能基因名字過程就會導致顯示不全，在與其它圖像拼接成一副圖的時候也會出現字號太小導致看不清楚的問題。

正好前幾天在群里有人詢問了這個圖有沒有其它的繪制辦法，今天就給大家帶來一個使用ggplot2繪制Extended error bar的方法。

數據準備

這里我將使用一套同一環境位點水體和沉積物16S擴增子測序的PICRUSt功能預測結果作為示例。

選擇的是KEGG L2水平的功能預測的相對豐度數據。

繪圖的數據文件有兩個，一個是豐度數據，另一個是樣本分組數據。

后臺回復“STAMP”獲取示例文件和完整代碼。

首先將數據導入R環境，我是首先過濾掉了平均豐度低于1%的功能分類。

data <- read.table("KEGG_L2.txt",header = TRUE,row.names = 1,sep = "t")
group <- read.table("group.txt",header = FALSE,sep = "t")
library(tidyverse)
data <- data*100
data <- data %>% filter(Apply(data,1,mean) > 1)

??如果使用不同分類學水平的微生物數據或著更深層次的功能注釋數據，由于物種或功能基因種類較多，可能會導致結果中具有差異的數目特別多，比如幾十個差異物種放在一副圖里基本上是不可能看清的，這種時候就要對數據進行過濾，去除低豐度的物種或基因，具體的過濾標準請根據自身數據情況自行確定。

統計學檢驗

在繪制Extended error bar之前首先要對數據進行差異顯著性檢驗，以選出豐度在不同組間具有顯著差異的物種或功能基因，這里以兩組數據為例，使用的檢驗方式通常為t-test和Wilcox秩和檢驗。

當分析數據符合正態分布時，使用t-test，如不符合正態分布，則使用Wilcox秩和檢驗。

t-test

首先對數據進行調整，構建用于t-test的數據框。

data <- t(data)
data1 <- data.frame(data,group$V2)
colnames(data1) <- c(colnames(data),"Group")
data1$Group <- as.factor(data1$Group)

由于R語言的t-test一次只能分析一列數據，在網上搜到了一個批量進行t-test的方法，感覺是最簡便的了。

首先使用select_if選擇格式為數字列，然后使用map_df分別對每一個列進行t-test，最后使用broom：tidy將結果整合在tidy的數據框中。

diff <- data1 %>%
select_if(is.numeric) %>%
map_df(~ broom::tidy(t.test(. ~ Group,data = data1)), .id = 'var')

最后對t-test的p值進行校正，保留校正后p值小于0.05的數據。

diff$p.value <- p.adjust(diff$p.value,"bonferroni")
diff <- diff %>% filter(p.value < 0.05)

秩和檢驗

秩和檢驗和上面的t-test一樣，只需要把代碼中的t.test換成wilcox.test就可以了。

diff1 <- data1 %>%
select_if(is.numeric) %>%
map_df(~ broom::tidy(wilcox.test(. ~ Group,data = data1)), .id = 'var')

diff1$p.value <- p.adjust(diff1$p.value,"bonferroni")
diff1 <- diff %>% filter(p.value < 0.05)

數據可視化

這個Extended error bar的結果圖整體分為兩個部分，左側為組建物種或基因豐度平均值的比較條形圖，右側為組間平均豐度及其95%置信區間的散點圖。

畫圖的思路是首先分別繪制左右兩幅圖，之后再拼接在一起，因此需要首先構建繪制兩幅圖所需的繪圖文件。

繪圖數據獲取

對于左側的組間豐度均值比較條形圖，我們首先根據差異性檢驗的結果從原始的豐度數據文件中提取具有顯著差異的列，之后按照分組計算其組內平均豐度，再轉換成ggplot繪圖所用的長格式數據框。

abun.bar <- data1[,c(diff$var,"Group")] %>%
gather(variable,value,-Group) %>%
group_by(variable,Group) %>%
summarise(Mean = mean(value))

對于右側散點圖所需要的數據，在上一步差異學檢驗的結果中已經給出了，我們只需要提取對應的數據列，之后根據差異的正負給其賦予對應的分組信息，最后按照差異豐度的大小對數據進行重新的排序

diff.mean <- diff[,c("var","estimate","conf.low","conf.high")]
diff.mean$Group <- c(ifelse(diff.mean$estimate >0,levels(data1$Group)[1],
levels(data1$Group)[2]))
diff.mean <- diff.mean[order(diff.mean$estimate,decreasing = TRUE),]

左側條形圖

在繪制條形圖之前，要把物種或基因名按照其組間豐度差異從大到小進行一個排序，以便圖像的美觀。

在畫圖之前首先要根據物種或基因豐度的差異大小，對變量的名稱設置一個因子并排序，以保證兩幅圖變量排序的一致性。

library(ggplot2)
cbbPalette <- c("#E69F00", "#56B4E9")
abun.bar$variable <- factor(abun.bar$variable,levels = rev(diff.mean$var))

接下來繪制圖像。

p1 <- ggplot(abun.bar,aes(variable,Mean,fill = Group)) +
scale_x_discrete(limits = levels(diff.mean$var)) +
coord_flip +
xlab("") +
ylab("Mean proportion (%)") +
theme(panel.background = element_rect(fill = 'transparent'),
panel.grid = element_blank,
axis.ticks.length = unit(0.4,"lines"),
axis.ticks = element_line(color='black'),
axis.line = element_line(colour = "black"),
axis.title.x=element_text(colour='black', size=12,face = "bold"),
axis.text=element_text(colour='black',size=10,face = "bold"),
legend.title=element_blank,
legend.text=element_text(size=12,face = "bold",colour = "black",
margin = margin(r = 20)),
legend.position = c(-1,-0.1),
legend.direction = "horizontal",
legend.key.width = unit(0.8,"cm"),
legend.key.height = unit(0.5,"cm"))for (i in 1:(nrow(diff.mean) - 1))
p1 <- p1 + annotate('rect', xmin = i+0.5, xmax = i+1.5, ymin = -Inf, ymax = Inf,
fill = ifelse(i %% 2 == 0, 'white', 'gray95'))

p1 <- p1 +
geom_bar(stat = "identity",position = "dodge",width = 0.7,colour = "black") +
scale_fill_manual(values=cbbPalette)

本來是不想用循環添加矩形的方式來實現灰白較差的底紋的，但是研究了半天也沒找到更方便的辦法，各位要是有更簡便的辦法歡迎私聊～～

這里要特別注意圖層的順序，一定要在添加舉行之后在繪制條形圖，不然會覆蓋掉。

右側散點圖

同樣要先對變量名進行一個因子排序，保證兩個圖像的變量順序的一致。

此外散點圖要在右側添加校正后的p-value，因此需要先對p-value取有效數字，之后轉換成文本。

diff.mean$var <- factor(diff.mean$var,levels = levels(abun.bar$variable))
diff.mean$p.value <- signif(diff.mean$p.value,3)
diff.mean$p.value <- as.character(diff.mean$p.value)

本來是想畫散點圖的同時直接把p-value的文本加上，但是由于這幅圖p-value的文本是添加在繪圖區以外的，研究了半天沒弄明白，后來給分成了兩個圖，一副圖只有散點圖，另一副圖只有右側的p-value文本。

散點圖的繪制。

p2 <- ggplot(diff.mean,aes(var,estimate,fill = Group)) +
theme(panel.background = element_rect(fill = 'transparent'),
panel.grid = element_blank,
axis.ticks.length = unit(0.4,"lines"),
axis.ticks = element_line(color='black'),
axis.line = element_line(colour = "black"),
axis.title.x=element_text(colour='black', size=12,face = "bold"),
axis.text=element_text(colour='black',size=10,face = "bold"),
axis.text.y = element_blank,
legend.position = "none",
axis.line.y = element_blank,
axis.ticks.y = element_blank,
plot.title = element_text(size = 15,face = "bold",colour = "black",hjust = 0.5)) +
scale_x_discrete(limits = levels(diff.mean$var)) +
coord_flip +
xlab("") +
ylab("Difference in mean proportions (%)") +
labs(title="95% confidence intervals")

for (i in 1:(nrow(diff.mean) - 1))
p2 <- p2 + annotate('rect', xmin = i+0.5, xmax = i+1.5, ymin = -Inf, ymax = Inf,
fill = ifelse(i %% 2 == 0, 'white', 'gray95'))

p2 <- p2 +
geom_errorbar(aes(ymin = conf.low, ymax = conf.high),
position = position_dodge(0.8), width = 0.5, size = 0.5) +
geom_point(shape = 21,size = 3) +
scale_fill_manual(values=cbbPalette) +
geom_hline(aes(yintercept = 0), linetype = 'dashed', color = 'black')

p-value文本的繪制。

p3 <- ggplot(diff.mean,aes(var,estimate,fill = Group)) +
geom_text(aes(y = 0,x = var),label = diff.mean$p.value,
hjust = 0,fontface = "bold",inherit.aes = FALSE,size = 3) +
geom_text(aes(x = nrow(diff.mean)/2 +0.5,y = 0.85),label = "P-value (corrected)",
srt = 90,fontface = "bold",size = 5) +
coord_flip +
ylim(c(0,1)) +
theme(panel.background = element_blank,
panel.grid = element_blank,
axis.line = element_blank,
axis.ticks = element_blank,
axis.text = element_blank,
axis.title = element_blank)

這樣所有的繪圖單元就齊了。

圖像拼接

使用patchwork對3個繪圖部分進行拼接。

library(patchwork)
p <- p1 + p2 + p3 + plot_layout(widths = c(4,6,2))

最后保存圖片就齊活了。

ggsave("stamp.pdf",p,width = 10,height = 4)

后臺回復“STAMP”獲取示例文件和完整代碼。

本文中差異性檢驗的代碼參考了：https://www.it1352.com/1581321.html

另外感謝“生信小白魚”提供循環添加矩形的代碼！！

10000+：菌群分析 寶寶與貓狗 梅毒狂想曲 提DNA發Nature Cell專刊 腸道指揮大腦

系列教程：微生物組入門 BIOStar 微生物組 宏基因組

專業技能：學術圖表 高分文章 生信寶典 不可或缺的人

一文讀懂：宏基因組 寄生蟲益處 進化樹

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