Editor's Note

生統和生信分析一直是咱們投必得關注的重要內容，雖然我們現在還沒有固定的編輯在寫，但是也會經常轉載我自己和其它幾個朋友在其它公眾號首發的一些相關干貨，今天這篇來源于“微生物生信”，帶大家看看堆疊柱狀圖?和沖擊圖是如何實現的，并且附代碼奧~

00. 堆疊柱狀圖和沖擊圖介紹

堆疊柱狀圖用于表征高維數據，例如微生物群落，代謝物組成是一種常見的圖表類型，所以使用的次數也非常多，理所應當被封裝成函數，一鍵繪制。

沖擊圖其實添加了堆疊柱狀圖之間的連線，這對于柱子之間的橫向比較會更加容易一些，同堆疊柱狀圖互為替代。

沖擊圖其實添加了堆疊柱狀圖之間的連線，這對于柱子之間的橫向比較會更加容易一些，同堆疊柱狀圖互為替代。

01. 一條函數和參數解釋

result = barplot(otu = "./otutab.txt",tax = "./taxonomy.txt",map = "./metadata.tsv",j = "Phylum",k= 0.01,rep = 6,axis_ord = "KO-OE-WT",label = FALSE,sd = TRUE)

otu = "./otutab.txt" ：OTU表格

tax = "./taxonomy.txt" ：物種注釋文件

map = "./metadata.tsv"：樣本分組文件

j = "Phylum" ：設置按照門為分類水平合并OTU，當然可以選擇不同的分類水平；j = "Class"；j = "Order"；j = "Family"；j = "Genus"；甚至是OTU，如果你不嫌圖例多的話。

k= 0.01：按照豐度過濾掉第豐度的物種，減少圖例信息。

rep = 6 ：測定每個處理重復數量。

axis_ord = "KO-OE-WT" ：柱子在x軸額排布順序。

label = FALSE：是否在柱子上添加標簽。

sd = TRUE :是否對堆疊柱狀圖添加標準

02. 處理本地輸出圖片外我們還可以通過result提取圖片

我設置函數返回圖形結果的目的是方便大家根據自己需求修改顏色等圖形屬性。

#提取堆疊柱狀圖
p = result[[1]]
p
#提取沖擊圖
p = result[[2]]
p

3. 附：主函數

注意一下幾個包，大家都看看裝上沒有，務必都安裝后在運行，方可成功。

library(phyloseq)
library(tidyverse)
library(vegan)
library(reshape2)
library("plyr")
library(ggalluvial)
library(ggplot2)

barplot = function(otu = "./otutab.txt",tax = "./taxonomy.txt",map = "./metadata.tsv",j = "Phylum",k= 0.01,rep = 6,axis_ord = "KO-OE-WT",label = TRUE){
library(phyloseq)
library(tidyverse)
library(vegan)
library(reshape2)
library("plyr")
library(ggalluvial)
library(ggplot2)
axis_order = strsplit(basename(axis_order), "-")[[1]]#導入otu表格
otu = read.delim(otu,row.names = 1)
head(otu)
otu = as.matrix(otu)
str(otu)#導入注釋文件
tax = read.delim(tax,row.names = 1)
head(tax)
tax = as.matrix(tax)# taxa_names(tax)#導入分組文件
map = read.delim(map,row.names = 1)
head(map)# #導入進化樹# tree = read.tree("./otus.tree")# tree
ps <- phyloseq(otu_table(otu, taxa_are_rows=TRUE),
sample_data(map) ,
tax_table(tax)# phy_tree(tree)
)
ps
ps1 = ps
i = ps1
colnames(tax_table(ps1))##這里我們過濾一定閾值的otu，會出現最后堆疊柱狀圖總體豐度高于100%的情況，這是合理的###########繪制不同分類等級的柱狀圖
Taxonomies <- i %>%
tax_glom(taxrank = j) %>% # agglomerate at Class level Class
transform_sample_counts(function(x) {x/sum(x)} )%>%# Transform to rel. abundance
psmelt() %>% # Melt to long format
filter(Abundance >= k) %>% # Filter out low abundance taxa
arrange(Phylum)# head(Taxonomies)# dim(Taxonomies)
colbar <- dim(unique(select(Taxonomies, one_of(j))))[1]
colors = colorRampPalette(c( "#CBD588", "#599861", "orange","#DA5724", "#508578", "#CD9BCD",
"#AD6F3B", "#673770","#D14285", "#652926", "#C84248",
"#8569D5", "#5E738F","#D1A33D", "#8A7C64","black"))(colbar)# mi = colorRampPalette(c( "#CBD588", "#599861", "orange","#DA5724", "#508578", "#CD9BCD",# "#AD6F3B", "#673770","#D14285", "#652926", "#C84248",# "#8569D5", "#5E738F","#D1A33D", "#8A7C64","black"))(colbar)##
Taxonomies$Abundance = Taxonomies$Abundance * 100
Taxonomies$Abundance = Taxonomies$Abundance/rep# head(Taxonomies)#按照分組求均值
colnames(Taxonomies) <- gsub(j,"aa",colnames(Taxonomies))
by_cyl <- group_by(Taxonomies, Group,aa)
zhnagxu2 = dplyr :: summarise(by_cyl, sum(Abundance))#colnames(zhnagxu2) = c("group", j,"Abundance")# head(zhnagxu2)##確定因子，這里通過求和按照從小到大的順序得到因子##長變寬# head(Taxonomies)
Taxonomies2 = dcast(Taxonomies,aa ~ Sample,value.var = "Abundance")
head(Taxonomies2)
Taxonomies2[is.na(Taxonomies2)] <- 0
aa = Taxonomies2# head(aa)
n = ncol(aa)#增加一行，為整列的均值，計算每一列的均值，2就是表示列
aa[n+1]=apply(aa[,c(2:ncol(aa))],1,sum)
colnames(aa)[n+1] <- c("allsum")# str(aa)
bb<- arrange(aa, allsum)# head(bb)
bb = bb[c(1,ncol(bb))]
cc<- arrange(bb, desc(allsum))# head(cc)##使用這個屬的因子對下面數據進行排序
head(zhnagxu2)
colnames(zhnagxu2) <- c("group","aa","Abundance")
zhnagxu2$aa = factor(zhnagxu2$aa,order = T,levels = cc$aa)
zhnagxu3 = plyr::arrange(zhnagxu2,desc(aa))# head(zhnagxu3)##制作標簽坐標，標簽位于頂端# Taxonomies_x = ddply(zhnagxu3,"group", transform, label_y = cumsum(Abundance))# head(Taxonomies_x )#標簽位于中部
Taxonomies_x = ddply(zhnagxu3,"group", transform, label_y = cumsum(Abundance) - 0.5*Abundance)
head(Taxonomies_x,20 )
Taxonomies_x$label = Taxonomies_x$aa#使用循環將堆疊柱狀圖柱子比較窄的別寫標簽，僅僅寬柱子寫上標簽
for(i in 1:nrow(Taxonomies_x)){if(Taxonomies_x[i,3] > 3){
Taxonomies_x[i,5] = Taxonomies_x[i,5]
}else{
Taxonomies_x[i,5] = NA
}
}##普通柱狀圖
p4 <- ggplot(Taxonomies_x , aes(x = group, y = Abundance, fill = aa, order = aa)) +
geom_bar(stat = "identity",width = 0.5,color = "black") +
scale_fill_manual(values = colors) +
theme(axis.title.x = element_blank()) +
theme(legend.text=element_text(size=6)) +
scale_y_continuous(name = "Abundance (%)")+
scale_x_discrete(limits = axis_order)if (label == TRUE) {
p4 = p4 +
geom_text(aes(y = label_y, label = label ),size = 4,fontface = "bold.italic")
}# print(p4)# install.packages("ggalluvial")
p4 = p4+theme_bw()+
scale_y_continuous(expand = c(0,0))+
theme(
panel.grid.major=element_blank(),
panel.grid.minor=element_blank(),
text=element_text(face = "bold"),
plot.title = element_text(vjust = -8.5,hjust = 0.1),
axis.title.y =element_text(size = 20,face = "bold",colour = "black"),
axis.title.x =element_text(size = 24,face = "bold",colour = "black"),
axis.text = element_text(size = 20,face = "bold"),
axis.text.x = element_text(colour = "black",size = 14,),
axis.text.y = element_text(colour = "black",size = 14),
legend.text = element_text(size = 15)#legend.position = "none"#是否刪除圖例
)
p4
FileName1 <- paste("./a2_",j,k,"_bar",".pdf", sep = "")
ggsave(FileName1, p4, width = 12, height =8 )##柱狀圖沖擊圖#stratum定義堆疊柱狀圖柱子內容，以weight定義柱子長度，alluvium定義連線
head(Taxonomies_x )
cs = Taxonomies_x $aa# head(cs)# as.factor(Taxonomies_x $Genus)# cs = as.character(Taxonomies_x $Genus)# cs1 = as.factor(cs)
cs1 = cs#提取真正的因子的數量
lengthfactor = length(levels(cs1))#提取每個因子對應的數量
cs3 = summary (as.factor(cs1))
cs4 = as.data.frame(cs3)
cs4$id = row.names(cs4)#對因子進行排序
df_arrange<- arrange(cs4, id)#對Taxonomies_x 對應的列進行排序
Taxonomies_x1<- arrange(Taxonomies_x , aa)
head(Taxonomies_x1)#構建flow的映射列Taxonomies_x
Taxonomies_x1$ID = factor(rep(c(1:lengthfactor), cs4$cs3))#colour = "black",size = 2,,aes(color = "black",size = 0.8)
p3 = ggplot(Taxonomies_x1,
aes(x = group, stratum = aa, alluvium = ID,
weight = Abundance,
fill = aa, label = aa)) +
geom_flow(stat = "alluvium", lode.guidance = "rightleft",
color = "black",size = 0.2,width = 0.3,alpha = .2) +
geom_bar(width = 0.45)+
geom_stratum(width = 0.45,size = 0.2) +#geom_text(stat = "stratum", size = 3,family="Times New Roman",fontface = "bold.italic") +#theme(legend.position = "none") +
scale_fill_manual(values = colors)+#ggtitle("fow_plot")+
scale_x_discrete(limits = axis_order)+
geom_text(aes(y = label_y, label = label ),size = 4,fontface = "bold.italic")+
labs(x="group",
y="Relative abundancce (%)",
title="")
# p3
if (label == TRUE) {
p3 = p3 +
geom_text(aes(y = label_y, label = label ),size = 4,fontface = "bold.italic")
}
p3 =p3+theme_bw()+
scale_y_continuous(expand = c(0,0))+
#geom_hline(aes(yintercept=0), colour="black", linetype=2) +
#geom_vline(aes(xintercept=0), colour="black", linetype="dashed") +
#scale_fill_manual(values = mi, guide = guide_legend(title = NULL))+
theme(
panel.grid.major=element_blank(),
panel.grid.minor=element_blank(),
text=element_text(face = "bold"),
plot.title = element_text(vjust = -8.5,hjust = 0.1),
axis.title.y =element_text(size = 20,face = "bold",colour = "black"),
axis.title.x =element_text(size = 24,face = "bold",colour = "black"),
axis.text = element_text(size = 20,face = "bold"),
axis.text.x = element_text(colour = "black",size = 14),
axis.text.y = element_text(colour = "black",size = 14),
legend.text = element_text(size = 15,face = "bold.italic")
#legend.position = "none"#是否刪除圖例
)
p3
FileName2 <- paste("./a2_",j,k,"_bar_flow",".pdf", sep = "")
ggsave(FileName2, p3, width = 12, height =8)
return(list(p4,p3))
}

好了本期就到這里了，歡迎大家評論區討論~

每周一篇研究生科研經驗分享，每周五篇文章教你寫SCI論文，歡迎關注專欄：

投必得科研軟件安裝使用手冊； 投必得：SCI期刊介紹與選擇； 投必得，教你寫論文；投必得統計分析大講堂；投必得科研生活解憂雜貨店

這里是論文編輯潤色專家，輸出科研干貨的投必得，我們下個回答再見ヾ(￣▽￣)Bye~Bye~

總結

以上是生活随笔為你收集整理的双柱状图柱子数量比较多_一条代码完成堆叠柱状图-冲击图的操作-终结版的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。