对于许多做微生物的同学来说，基于高通量测序的16S扩增子测序是对细菌群落组成及多样性进行分析的重要手段。
细菌的16S rDNA全长约1500 bp，由V1-V9九个可变区和若干保守区组成，而16S扩增子测序的原理就是通过对细菌的16S rDNA基因上的可变区域进行扩增，然后利用变异区域的差异来对不同菌属的细菌进行分类鉴定。
16S扩增子测序结束后，我们会拿到测序的原始序列，然后经过拼接、质控等一系列流程后，需要进行OTU（Operational Taxonomic Units，操作分类单元）分析以拿到OTU丰度表，基于OTU丰度表我们可以进行的分析项目主要包括：1）Alpha多样性分析；2）Beta多样性分析；3）物种分类及丰度分析；4）物种差异及关联分析。
今天要介绍的内容是基于R语言的微生物群落组成及多样性分析——Alpha多样性分析及可视化：
1、工作目录设置及相关包安装、加载
rm(list=ls())#clear Global Environment
#设置工作目录
setwd('D:\\桌面\\α-diversity')
#安装包
install.packages("ggplot2")
install.packages("ggpubr")
install.packages("ggsignif")
install.packages("vegan")
install.packages("ggprism")
install.packages("picante")
install.packages("dplyr")
install.packages("RColorBrewer")
#加载包
library(ggplot2)
library(ggpubr)
library(ggsignif)
library(ggprism)
library(vegan)
library(picante)
library(dplyr)
library(RColorBrewer)
2、Alpha多样性指数的计算，主要包括ACE、Chao、Shannon、Simpson等指数的计算
##导入数据，所需是数据行名为样本名、列名为OTUxxx的数据表
df <- read.table("otu.txt",header = T, row.names = 1, check.names = F)
#使用vegan包计算多样性指数
Shannon <- diversity(df, index = "shannon", MARGIN = 2, base = exp(1))
Simpson <- diversity(df, index = "simpson", MARGIN = 2, base = exp(1))
Richness <- specnumber(df, MARGIN = 2)#spe.rich =sobs
###将以上多样性指数统计成表格
index <- as.data.frame(cbind(Shannon, Simpson, Richness))
tdf <- t(df)#转置表格
tdf<-ceiling(as.data.frame(t(df)))
#计算obs，chao，ace指数
obs_chao_ace <- t(estimateR(tdf))
obs_chao_ace <- obs_chao_ace[rownames(index),]#统一行名
#将obs，chao，ace指数与前面指数计算结果进行合并
index$Chao <- obs_chao_ace[,2]
index$Ace <- obs_chao_ace[,4]
index$obs <- obs_chao_ace[,1]
#计算Pielou及覆盖度
index$Pielou <- Shannon / log(Richness, 2)
index$Goods_coverage <- 1 - colSums(df ==1) / colSums(df)
#导出表格
write.table(cbind(sample=c(rownames(index)),index),'diversity.index.txt', row.names = F, sep = '\t', quote = F)
3、差异性计算及绘图
读入数据文件及分组文件
index$samples <- rownames(index)#将样本名写到文件中
#读入分组文件
groups <- read.delim('group.txt',header = T, stringsAsFactors = F)
colnames(groups)[1:2] <- c('samples','group')#改列名
#合并分组信息与多样性指数
df2 <- merge(index,groups,by = 'samples')
绘图
#Shannon
p1 <- ggplot(df2,aes(x=group,y=Shannon))+#指定数据
stat_boxplot(geom = "errorbar", width=0.1,size=0.8)+#添加误差线,注意位置，放到最后则这条先不会被箱体覆盖
geom_boxplot(aes(fill=group), #绘制箱线图函数
outlier.colour="white",size=0.8)+#异常点去除
theme(panel.background =element_blank(), #背景
axis.line=element_line(),#坐标轴的线设为显示
plot.title = element_text(size=14))+#图例位置
# scale_fill_manual(values=c("#ffc000","#a68dc8","blue"))+#指定颜色
geom_jitter(width = 0.2)+#添加抖动点
geom_signif(comparisons = list(c("A","B"),
c("A","C"),
c("B","C")),# 设置需要比较的组
map_signif_level = T, #是否使用星号显示
test = t.test, ##计算方法
y_position = c(3,3.5,3.25),#图中横线位置 设置
tip_length = c(c(0,0),
c(0,0),
c(0,0)),#横线下方的竖线设置
size=0.8,color="black")+
theme_prism(palette = "candy_bright",
base_fontface = "plain", # 字体样式，可选 bold, plain, italic
base_family = "serif", # 字体格式，可选 serif, sans, mono, Arial等
base_size = 16, # 图形的字体大小
base_line_size = 0.8, # 坐标轴的粗细
axis_text_angle = 45)+ # 可选值有 0，45，90，270
scale_fill_prism(palette = "candy_bright")+
theme(legend.position = 'none')#去除图例
p1
image.png
#Simpson
p2 <- ggplot(df2,aes(x=group,y=Simpson))+#指定数据
stat_boxplot(geom = "errorbar", width=0.1,size=0.8)+#添加误差线,注意位置，放到最后则这条先不会被箱体覆盖
geom_boxplot(aes(fill=group), #绘制箱线图函数
outlier.colour="white",size=0.8)+#异常点去除
theme(panel.background =element_blank(), #背景
axis.line=element_line(),#坐标轴的线设为显示
plot.title = element_text(size=14))+#图例位置
# scale_fill_manual(values=c("#ffc000","#a68dc8","blue"))+#指定颜色
geom_jitter(width = 0.2)+#添加抖动点
geom_signif(comparisons = list(c("A","B"),
c("A","C"),
c("B","C")),# 设置需要比较的组
map_signif_level = T, #是否使用星号显示
test = t.test, ##计算方法
y_position = c(1,1.1,1.05),#图中横线位置 设置
tip_length = c(c(0,0),
c(0,0),
c(0,0)),#横线下方的竖线设置
size=0.8,color="black")+
theme_prism(palette = "candy_bright",
base_fontface = "plain", # 字体样式，可选 bold, plain, italic
base_family = "serif", # 字体格式，可选 serif, sans, mono, Arial等
base_size = 16, # 图形的字体大小
base_line_size = 0.8, # 坐标轴的粗细
axis_text_angle = 45)+ # 可选值有 0，45，90，270
scale_fill_prism(palette = "candy_bright")+
theme(legend.position = 'none')#去除图例
p2
image.png
#Ace
p3 <- ggplot(df2,aes(x=group,y=Ace))+#指定数据
stat_boxplot(geom = "errorbar", width=0.1,size=0.8)+#添加误差线,注意位置，放到最后则这条先不会被箱体覆盖
geom_boxplot(aes(fill=group), #绘制箱线图函数
outlier.colour="white",size=0.8)+#异常点去除
theme(panel.background =element_blank(), #背景
axis.line=element_line(),#坐标轴的线设为显示
plot.title = element_text(size=14))+#图例位置
# scale_fill_manual(values=c("#ffc000","#a68dc8","blue"))+#指定颜色
geom_jitter(width = 0.2)+#添加抖动点
geom_signif(comparisons = list(c("A","B"),
c("A","C"),
c("B","C")),# 设置需要比较的组
map_signif_level = T, #是否使用星号显示
test = t.test, ##计算方法
y_position = c(55,65,60),#图中横线位置 设置
tip_length = c(c(0,0),
c(0,0),
c(0,0)),#横线下方的竖线设置
size=0.8,color="black")+
theme_prism(palette = "candy_bright",
base_fontface = "plain", # 字体样式，可选 bold, plain, italic
base_family = "serif", # 字体格式，可选 serif, sans, mono, Arial等
base_size = 16, # 图形的字体大小
base_line_size = 0.8, # 坐标轴的粗细
axis_text_angle = 45)+ # 可选值有 0，45，90，270
scale_fill_prism(palette = "candy_bright")+
theme(legend.position = 'none')#去除图例
p3
image.png
#Chao
p4 <- ggplot(df2,aes(x=group,y=Chao))+#指定数据
stat_boxplot(geom = "errorbar", width=0.1,size=0.8)+#添加误差线,注意位置，放到最后则这条先不会被箱体覆盖
geom_boxplot(aes(fill=group), #绘制箱线图函数
outlier.colour="white",size=0.8)+#异常点去除
theme(panel.background =element_blank(), #背景
axis.line=element_line(),#坐标轴的线设为显示
plot.title = element_text(size=14))+#图例位置
# scale_fill_manual(values=c("#ffc000","#a68dc8","blue"))+#指定颜色
geom_jitter(width = 0.2)+#添加抖动点
geom_signif(comparisons = list(c("A","B"),
c("A","C"),
c("B","C")),# 设置需要比较的组
map_signif_level = T, #是否使用星号显示
test = t.test, ##计算方法
y_position = c(55,65,60),#图中横线位置 设置
tip_length = c(c(0,0),
c(0,0),
c(0,0)),#横线下方的竖线设置
size=0.8,color="black")+
theme_prism(palette = "candy_bright",
base_fontface = "plain", # 字体样式，可选 bold, plain, italic
base_family = "serif", # 字体格式，可选 serif, sans, mono, Arial等
base_size = 16, # 图形的字体大小
base_line_size = 0.8, # 坐标轴的粗细
axis_text_angle = 45)+ # 可选值有 0，45，90，270
scale_fill_prism(palette = "candy_bright")+
theme(legend.position = 'none')#去除图例
p4
image.png
拼图
library("gridExtra")
library("cowplot")
plot_grid(p1,p2,p3,p4, labels=c('A','B','C','D'), ncol=2, nrow=2)#拼图及标注
ggsave('xxx.pdf',width=12,height = 4)
image.png
这里我们主要使用ggprism包美化图形，我们可以通过names(ggprism_data$themes)输出该包中的所有主题后，然后在theme_prism()及scale_fill_prism()中通过更改主题名称调整作图风格
names(ggprism_data$themes)
image.png
