群落堆叠柱状图+冲击图绘制
群落堆叠柱状图示例文件
1. 某一分类级别丰度文件：
2.分组文件：
接下来开始绘制：
准备一下颜色信息 因为属或种水平物种比较多，可以多准备一些
#准备一下颜色信息 因为属或种水平物种比较多，可以多准备一些 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)
准备作图文件
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')
普通柱状图绘制
#普通柱状图 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)
冲击图绘制
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)
合并起来作图
有时候我们不想看各个样品的组成情况，可以通过求和或者求均值的情况，将组合并只看俩组均值丰度的比较情况
#分组画图 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
除此之外还有一种方式是批量的循环作图
假设文件夹下有各个层级的丰度表如下：
写个循环一次性展示各层级平均物种丰度大于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) }
好啦 出图成功，Nice！大家也试试吧！
