基于R语言的微生物群落组成多样性分析——α多样性及其可视化

点击上方蓝字关注我哦

引言

      对于许多做微生物的同学来说，基于高通量测序的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、差异性计算及绘图

1）读入数据及分组文件

#读入文件
index <- read.delim('diversity.index.txt', header = T, row.names = 1)
##figure:take shannon for example
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')

2）绘图

--Shannon

#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

--Simpson

#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

--Ace

#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

--Chao

#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

4、拼图

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)

---

END

以上就是今天所有的内容，希望对您有所帮助！！！

大家如果有什么问题请在后台留言，小编看到后会在第一时间回复

欢迎大家点赞、转发并点亮在看，让更多朋友看到！！！

也欢迎大家关注我的公众号