绘制堆叠柱形图（ggplot 2）使用群落微生物物种丰度表数据

2021-05-25 本文已影响0人 fafu生信小蘑菇

群落微生物物种丰度表-绘制堆叠柱形图（ggplot 2）

今天我们先来讲一下如何利用细菌群落的高通量测序数据来绘制堆叠柱形图来展示不同样品中优势菌的相对多度。

话不多说，直接上数据。

1 加载数据及预处理

1.1查看数据

现在先解释一下，这组数据一个6个样本（sample_1 到sample_6），每个样本有3个重复，一共有18个样本进行高通量测序，高通量测序样本名称为（A1-A18）。主要有两个表，一个名称为genus的属水平的丰度表；一个名称为group的样本分组信息表。

作为例子数据是属水平的物种丰度表，如图1-1所示：

图1.1-1 属水平的物种丰度表

group的样本分组信息表，如图·1-2所示：

这里注意一下，由于前期准备原因把样本名称输成了sample-1，导致后面计算出错，所以这里修改样本名称为(sample_1到sample_6),特此纠正。

图1.1-2 group分组信息表

图1.1-3 测序样本名与样本名称的对应关系

1.2加载数据及数据预处理

#设置工作路径
setwd("C:/Users/shanpengloveforever/Desktop/图/微信") 
#加载genus 物种丰度表
data<-read.table("genus.txt",header=T,sep="\t",row.names=1)

data$sum <- rowSums(data) #求每一行的和
# 按每行的和降序排列
data1 <- data[order(data$sum, decreasing=TRUE), ]
#data2 <- data1[order(data1$sum, decreasing=FALSE), ] 按每行的和升序排列
data1 <- data1[,-19] #删除sum列，为了计算后面分组的平均值

图1.2-1

图1.2-2

#按行求指定列平均值，并且把算好的平均值添加data1数据框
data1$sample_1 <- apply(data1[,1:3], 1, mean) 
data1$sample_2 <- apply(data1[,4:6], 1, mean)
data1$sample_3 <- apply(data1[,7:9], 1, mean)
data1$sample_4 <- apply(data1[,10:12], 1, mean)
data1$sample_5<- apply(data1[,13:15], 1, mean)
data1$sample_6 <- apply(data1[,16:18], 1, mean)

#提取出已经算好的平均值到data2数据集
data2 <- data1[,19:24] 

#取出丰富度排名前10的物种，并且计算相对丰度
#由于之间已经按照每行的和进行过升序排列，所以可以直接去前10行
data3<- data2[1:10,]/apply(data2,2,sum)

data4 <- 1-apply(data3, 2, sum) #计算剩下物种的总丰度
#合并数据
data3 <- rbind(data3,data4)

图1.2-3

图1.2-4

图1.2-5

图1.2-6

图1.2-7

使用R语言将data3 数据集导出

write.table (data3, file ="data3.csv",sep =",", quote =FALSE) #将数据导出

图1.2-8

在Excel中修改data3 数据集，并且另存为genus1文本文件

图1.2-9

加载新数据集genus1

#导入修改好的数据
data3 <- read.table("genus1.txt",header=T,sep="\t",row.names=1)
#查看数据
row.names(data3)
colnames(data3)
apply(data3, 2, sum)

图1.2-10

2. 使用ggplot2 进行绘制堆叠柱形图

2.1 加载group分组信息及数据集的组合

#加载包
library(reshape2)
library(ggplot2)

#把data3 数据整理成 ggplot2 作图格式
#将菌名添加到data3里面，为了后面的数据转化
data3$Taxonomy <- factor(rownames(data3), levels = rev(rownames(data3)))
#宽数据转化为长数据
data4 <- melt(data3, id = 'Taxonomy')

图2.1-1

图2.1-2

#加载group分组信息表
group <- read.delim('group.txt', sep = '\t', stringsAsFactors = FALSE)

names(data4)[2] <- 'sample'  #修改列名
data5 <- merge(data4, group, by = 'sample')

图2.1-3

图2.1-4

2.2 使用ggplot2绘图

p<- ggplot(data5, aes(x=sample, y=100 * value, fill = Taxonomy)) +
  #数据输入：样本、物种、丰度
  geom_col(position = 'stack', width = 0.6) + # stack：堆叠图
  scale_y_continuous(expand=c(0, 0))+# 调整y轴属性，使柱子与X轴坐标接触
  scale_fill_manual(values =  rev(c('#FF0000', 
                                    '#FF88C2', '#FF00FF', '#9999FF', '#33FFFF',
                                    '#33FF33', '#D1BBFF', '#770077', '#EE7700', 
                                    '#CCEEFF', '#0000AA'))) + #手动修改颜色
  labs(x = 'Samples', y = '相对分度\n Relative Abundance(%)') + #设置X轴和Y轴的信息
  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))

p

图2.2-1

将绘制好的堆叠柱状图，保存为pdf和png格式

ggsave(filename = "genus.pdf",
       p,
       width=10,
       heigh=8)
ggsave('genus.png', p, width = 10, height = 8)


#补充
  theme(axis.text.x=element_text(angle=45, hjust=1))
# angle：调整横轴标签倾斜角度
# hjust：上下移动横轴标签

今天的内容就是这些，主要是数据处理和ggplot2 绘制堆叠柱状图，有什么不懂的可以私聊我。