群落堆叠柱状图+冲击图绘制

群落堆叠柱状图示例文件
1.某一分类级别丰度文件：

图片.png

2.分组文件：

图片.png

OK，接下来开始绘制：

1.准备一下颜色信息 因为属或种水平物种比较多，可以多准备一些

#准备一下颜色信息 因为属或种水平物种比较多，可以多准备一些
library(RColorBrewer)
rm(list = ls())
m1 = brewer.pal(9,"Set1")
m2 = brewer.pal(12,"Set3")
Palette1 <- c("#B2182B","#E69F00","#56B4E9","#009E73","#F0E442","#0072B2","#D55E00","#CC79A7","#CC6666","#9999CC","#66CC99","#999999","#ADD1E5")
Palette2 <- c('blue', 'orange', 'green', 'yellow', 'red', 'hotpink', 'cyan','purple', 'burlywood1','skyblue','grey','#8b8378','#458b74','#f0ffff','#eeb422','#ee6aa7','#8b3a62','#cd5c5c','#ee6363','#f0e68c','#e6e6fa','#add8e6','#bfefff','#f08080','#d1eeee','#7a8b8b','#8b814c','#8b5f65','gray')
mix <- c(m1,Palette2,Palette1,m2)

2. 准备作图文件

  phylum <- read.delim('test.txt',row.names = 1,sep = '\t',stringsAsFactors = FALSE)
 #过滤平均丰度大于0.001的物种
  phylum_filter <- phylum[rowMeans(phylum)>0.001,] 
  phylum_filter$sum <- rowSums(phylum_filter)
  #求各类群的丰度总和，并排序
  phylum_filter <-  phylum_filter[order(phylum_filter$sum,decreasing = TRUE),]
  #删除最后一列和,并取Top10物种作图
  phylum_top10 <- phylum_filter[c(1:10),-ncol(phylum_filter)]
 # 剩余的物种合并为Others
  phylum_top10['Others',] <- 1-colSums(phylum_top10)
 
 #最后一个设置为灰色
  colour <- mix[1:nrow(phylum_top10)]
  #个人喜欢将最后一个设置为灰色，也就是Others
  colour[length(colour)] <- 'gray' 

 #设置因子，改为长格式
  phylum_top10$Taxonomy <- factor(rownames(phylum_top10),levels = rev(rownames(phylum_top10)))
  phylum_top10 <- melt(phylum_top10,id = 'Taxonomy')

  #添加分组 合并信息
  group <- read.delim('sample.txt', sep = '\t', stringsAsFactors = FALSE)
  names(group)[1] <- 'variable'
  phylum_top10 <- merge(phylum_top10, group, by = 'variable')

3. 普通柱状图绘制

#普通柱状图
  p1 <- ggplot(phylum_top10, aes(variable, 100 * value, fill = Taxon)) +
    #facet_wrap(~group,scales = 'free_x',ncol = 2)+
    geom_col(position = 'stack', width = 0.6)  +
    scale_fill_manual(values =  rev(c(colour))) +
    theme_classic()+scale_y_continuous(expand = c(0,0))+
    labs(x = '', y = 'Relative Abundance(%)')+
    theme(panel.grid = element_blank(), panel.background = element_rect(color = 'black', fill = 'transparent'), strip.text = element_text(size = 12)) +
    theme(axis.text = element_text(size = 12), axis.title = element_text(size = 13), legend.text = element_text(size = 11),legend.background = element_blank())
#小技巧若之前没有设置value*100,纵轴变为百分比形式可用scale_y_continuous(labels = scales::percent),同时添加上%
#更多theme细节可自行调节
  p1
  ggsave('Genu.pdf',p1,height = 7,width = 8)

图片.png
冲击图绘制

  p2 <- ggplot(data = phylum_top10,aes(x = variable, y = value*100, alluvium = Taxonomy, stratum = Taxonomy))+
    geom_alluvium(aes(fill = Taxonomy),alpha = .5,width = 0.5)+ 
    geom_stratum(aes(fill = Taxonomy),width = 0.5)+scale_fill_manual(values = rev(c(colour)))+
    theme_classic()+scale_y_continuous(expand = c(0,0))+
    labs(x = '', y = 'Relative Abundance(%)')+
    geom_flow(alpha = 0.5)+ #设置透明度
    theme(panel.grid = element_blank(), panel.background = element_rect(color = 'black', fill = 'transparent'), strip.text = element_text(size = 12)) +
    theme(axis.text = element_text(size = 12), axis.title = element_text(size = 13), legend.text = element_text(size = 11),legend.background = element_blank())
  p2
ggsave('Genu.pdf',p2,height = 7,width = 12)

图片.png

合并起来作图
有时候我们不想看各个样品的组成情况，可以通过求和或者求均值的情况，将组合并只看俩组均值丰度的比较情况

  #分组画图
  phylum_top10$D <- rowMeans(phylum_top10[,1:6])
  phylum_top10$H <- rowMeans(phylum_top10[,7:12])
  phylum_top10 <- select(phylum_top10,D,H)
  colour <- mix[1:nrow(phylum_top10)]#设置对应属的个数
  colour[length(colour)] <- 'gray' #最后一个设置为灰色
  phylum_top10$Taxonomy <- factor(rownames(phylum_top10),levels = rev(rownames(phylum_top10)))
  phylum_top10 <- melt(phylum_top10,id = 'Taxonomy')

#作图 普通柱状图
  p1 <- ggplot(phylum_top10, aes(variable, 100 * value, fill = Taxonomy)) +
    #facet_wrap(~group,scales = 'free_x',ncol = 2)+
    geom_col(position = 'stack', width = 0.6)  +
    scale_fill_manual(values =  rev(c(colour))) +
    theme_classic()+scale_y_continuous(expand = c(0,0))+
    labs(x = '', y = 'Relative Abundance(%)')+
    theme(panel.grid = element_blank(), panel.background = element_rect(color = 'black', fill = 'transparent'), strip.text = element_text(size = 12)) +
    theme(axis.text = element_text(size = 12), axis.title = element_text(size = 13), legend.text = element_text(size = 11))
  p1

#冲击图
 p2 <- ggplot(data = phylum_top10,aes(x = variable, y = value*100, alluvium = Taxonomy, stratum = Taxonomy))+
    geom_alluvium(aes(fill = Taxonomy),alpha = .5,width = 0.5)+ 
    geom_stratum(aes(fill = Taxonomy),width = 0.5)+scale_fill_manual(values = rev(c(colour)))+
    theme_classic()+scale_y_continuous(expand = c(0,0))+
    labs(x = '', y = 'Relative Abundance(%)')+
    #geom_flow(alpha = 0.5)+ #设置透明度
    theme(panel.grid = element_blank(), panel.background = element_rect(color = 'black', fill = 'transparent'), strip.text = element_text(size = 12)) +
    theme(axis.text = element_text(size = 12), axis.title = element_text(size = 13), legend.text = element_text(size = 11),legend.background = element_blank())
  p2

图片.png

~~当然有时候我们不会单画某一个层级的图，很多层级都需要做的时候，一行行的回车简直累死，此时我们就可以将上边代码写在函数里运行会方便许多的
如下示例，可能还有更方便的请自行摸索：
写个函数作图：

#导入文件
input_file <- read.delim('Single_Phylum.txt',row.names = 1,sep = '\t',stringsAsFactors = FALSE)
Bar_plot <- function(input_file,top){
  #中间代码和上面一样 ‘从加载包到保存的代码’
}
Bar_plot(phylum = input_file,top = 20) #导出top前20的物种出图

除此之外还有一种方式是批量的循环作图
假设文件夹下有各个层级的丰度表如下：

图片.png
写个循环一次性展示各层级平均物种丰度大于1%的物种：

temp=list.files(path="./",pattern="Single.*.txt") #匹配文件
dir.create('plot',recursive = TRUE) # 创建出图的目录
#循环遍历
for (i in 1:length(temp)){
  phylum <- read.delim(temp[i],row.names = 1,sep = '\t',stringsAsFactors = FALSE)
  phylum_filter <- phylum[rowMeans(phylum)>0.01,]#筛选
  
  phylum_filter$sum <- rowSums(phylum_filter)#求各类群的丰度总和，并排序
  phylum_filter <- phylum_filter[order(phylum_filter$sum,decreasing = TRUE),]
  #phylum_filter <- phylum_filter[!grepl('unclassi',rownames(phylum_filter)),]#可以去除unclassified的物种
  phylum_top10 <- phylum_filter[1:nrow(phylum_filter),-ncol(phylum_filter)]
  phylum_top10['Others',] <- 1-colSums(phylum_top10)
  
  colour <- mix[1:nrow(phylum_top10)]#设置对应属的个数
  colour[length(colour)] <- 'gray' #最后一个设置为灰色
  
  phylum_top10$Taxonomy <- factor(rownames(phylum_top10),levels = rev(rownames(phylum_top10)))
  phylum_top10 <- melt(phylum_top10,id = 'Taxonomy')
  
  #添加分组
  group <- read.delim('sample.txt', sep = '\t', stringsAsFactors = FALSE)
  names(group)[1] <- 'variable'
  phylum_top10 <- merge(phylum_top10, group, by = 'variable')
  
  id <- strsplit(temp[i],split = '.t')[[1]][1]  #提取id 方便后面保存
  # 或者id <- unlist(strsplit(temp[i],split = '.t'))[1]
  
  #普通柱状图
  p1 <- ggplot(phylum_top10, aes(variable, 100 * value, fill = Taxonomy)) +
    #facet_wrap(~group,scales = 'free_x',ncol = 2)+
    geom_col(position = 'stack', width = 0.6)  +
    scale_fill_manual(values =  rev(c(colour))) +
    theme(panel.grid = element_blank(), panel.background = element_rect(color = 'black', fill = 'transparent'), strip.text = element_text(size = 12)) +
    theme(axis.text = element_text(size = 12), axis.title = element_text(size = 13), legend.title = element_blank(), legend.text = element_text(size = 11))+
    theme_classic()+scale_y_continuous(expand = c(0,0))+
    labs(x = '', y = 'Relative Abundance(%)')
  #facet_wrap(~Group, scales = 'free_x', ncol = 2) 
  #theme(axis.text.x=element_text(colour="black",size=12,face = "bold",angle = -90,vjust = 0.5,hjust = 0)) 
  # theme(text = element_text(family = "Times"))
  p1
  ggsave(paste('plot1/Taxa_bar1_',id,'.tiff',sep = ''),p1,height = 7,width = 8)
  #冲击图
  p2 <- ggplot(data = phylum_top10,aes(x = variable, y = value*100, alluvium = Taxonomy, stratum = Taxonomy))+
    geom_alluvium(aes(fill = Taxonomy),alpha = .5,width = 0.6)+ 
    geom_stratum(aes(fill = Taxonomy),width = 0.6)+scale_fill_manual(values = rev(c(colour)))+
    theme_classic()+scale_y_continuous(expand = c(0,0))+
    labs(x = '', y = 'Relative Abundance(%)')+
    theme(panel.grid = element_blank(), panel.background = element_rect(color = 'black', fill = 'transparent'), strip.text = element_text(size = 12)) +
    theme(axis.text = element_text(size = 12), axis.title = element_text(size = 13), legend.text = element_text(size = 11))
  p2
  ggsave(paste('plot1/Taxa_bar2_',id,'.tiff',sep = ''),p2,height = 7,width = 10)
}

图片.png

~好啦 出图成功，Nice！大家也试试吧！

下面是刚创建个人公众号，会定时更新R、linux、python，组学方面的学习内容，请多多支持呦~~

image