本文使用某宏基因组测序数据中的种水平类群丰度表，并直接根据所有种相对丰度总和，提取出top10丰度的种类群，作图展示这些种的丰度组成，top10外的类群合并为“Others”。
-
作图数据说明
-
文件“species_RPM_table-NC.txt”为某宏基因组测序所得的种水平类群丰度(RPM)表格（即界门纲目科属种的“门”这一水平），第一列（Taxonomy）为每类物种的名称，第一行为样本ID（名称及顺序要与以下两个文件一致）；每一行为一类种，交叉区域为每类种在各样本中的丰度(RPM)，内容如下：
species_RPM_table-NC.txt -
文件“sample_frequency.txt”为各样本的总序列数（total mapped reads），第1列为样本ID，第2列为物种数量，第3列为总序列数，内容如下：sample_frequency.txt
-
文件“sample_group.txt”为各样本的分组信息，第一列为各样本名称；第二列（group）为各样本的分组信息，其内容如下：sample_group.txt
-
绘图脚本说明
操作环境: Rstudio
脚本保存位置：D:\ ...\new_database\S_table\RPM\cutoff-0.5\draw_taxa_barplot.r
####读取并挑选 top10 丰度种(RPM)，并绘制物种丰度堆叠柱状图
#读取数据：特征表
species <- read.delim('species_RPM_table-NC.txt', row.names = 1, sep = '\t', stringsAsFactors = FALSE, check.names = FALSE)
#求各类群的丰度总和，并排序
species$sum <- rowSums(species)
species <- species[order(species$sum, decreasing = TRUE), ]
#挑选 top10 种类群，并将 top10 外的类群合并为“Others”
species_top10 <- species[1:10, -ncol(species)]
## 读取数据：样本的总特征序列数
sample_frequency <- read.delim('sample_frequency.txt', row.names = 1, sep = '\t', stringsAsFactors = FALSE, check.names = FALSE)
## others：用每个样本对应的总特征序列数减去top10
species_top10['Others', ] <- sample_frequency$feature_count - colSums(species_top10)
## others：当使用相对丰度表时，用1减去top10物种的相对丰度
# species_top10['Others', ] <- 1 - colSums(species_top10)
#可选输出（如 csv 格式）
write.csv(species_top10, 'species_top10.csv', quote = FALSE)
####ggplot2 作图
library(reshape2) #用于排列数据
library(ggplot2) #ggplot2 作图
#调整数据布局
species_top10$Taxonomy <- factor(rownames(species_top10), levels = rev(rownames(species_top10)))
species_top10 <- melt(species_top10, id = 'Taxonomy')
#添加分组：variable列为各样本名称，Taxonomy列为各类群名称，value列为各样本中各类群的相对丰度，group列为各样本的分组信息
group <- read.delim('sample_group.txt', sep = '\t', stringsAsFactors = FALSE)
names(group)[1] <- 'variable'
species_top10 <- merge(species_top10, group, by = 'variable')
## 指定横坐标（样本ID）顺序
species_top10$variable=factor(species_top10$variable,levels = c("A0065","A4313","C0065","C4313","C4445","V0065","V4313","V4445"))
#绘制带分面的柱状图
# p <- ggplot(species_top10, aes(variable, 100 * value, fill = Taxonomy)) + ##相对丰度乘以100，求百分比
p <- ggplot(species_top10, aes(variable, value, fill = Taxonomy)) +
geom_col(position = 'stack', width = 0.6) +
facet_wrap(~group, scales = 'free_x', ncol = 3) + ### ncol 分组数量
scale_fill_manual(values = rev(c('blue', 'orange', 'green', 'yellow', 'red', 'hotpink', 'cyan','purple', 'burlywood1', 'skyblue', 'gray'))) +
# labs(x = '', y = 'Relative Abundance(%)') +
labs(x = '', y = 'RPM') +
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 = 10), axis.title = element_text(size = 13), legend.title = element_blank(), legend.text = element_text(size = 11))
## 使用geom_col()绘制柱状图，position = 'stack'，堆叠柱状图，width = 0.6，柱子宽度；labs()设置坐标轴标签；theme()中，axis.text和axis.title分别指定坐标轴刻度字体大小及标签字体大小，legend.text指定图例字体大小。
## 保存图片：结果如下图
ggsave('species_taxa_barplot.png', p, width = 8, height = 6)
#### 其他方法：barplot 作图（一个简单示例）
#定义作图整体分布
png('barplot_plot.png', width = 1000, height = 700)
par(xpd = TRUE, mar = par()$mar + c(1, 3, 1, 16))
#barplot作图，包含了颜色、字体大小、图例位置等信息
barplot(as.matrix(100 * phylum_top10), col = c('blue', 'orange', 'green', 'yellow', 'red', 'hotpink', 'cyan','purple', 'burlywood1', 'skyblue', 'gray'),
legend = rownames(species_top10),
cex.axis = 2, cex.names = 2, ylim = c(0, 100), las = 1, width = 0.5, space = 0.5, beside = FALSE,
args.legend = list(x = 'right', bty = 'n', inset = -0.18, cex = 2, y.intersp = 1.2, x.intersp = 0.7, text.width = 1))
mtext('Relative Abundance(%)', cex = 2, side = 2, line = 4)
dev.off()
species_taxa_barplot.png
