基于R语言的微生物群落组成分析—差异OTU筛选及展示

安装、加载R包

rm(list=ls())#clear Global Environment
setwd('D:\\桌面\\差异OTU筛选')#设置工作路径

#安装所需R包
# install.packages("BiocManager")
# library(BiocManager)
# BiocManager::install("DESeq2")
# install.packages("ggplot2")
# install.packages('ggrepel')
#加载包
library(DESeq2)
library(ggplot2)
library(ggrepel)

加载数据

#OTU数据
otu <- read.table(file="otu.txt",sep="\t",header=T,
                  check.names=FALSE ,row.names=1)
#分组数据
group <- read.table(file="group.txt",sep="\t",
                    header=T,check.names=FALSE,row.names=1 )
group$group<-factor(group$group,levels = c('con','tes'))

image.png
image.png

OTU差异分析

1、构建DESeqDataSet对象

df_dds <- DESeqDataSetFromMatrix(countData = otu, colData = group, design = ~group)

2、差异分析

df_dds_res <- DESeq(df_dds)
suppressMessages(df_dds_res)

image.png

3、提取分析结果

#提取分析结果
df_res <- results(df_dds_res)
# 根据p-value进行重新排序
df_res = df_res[order(df_res$pvalue),]
df_res #查看结果

image.png

4、统计结果

summary(df_res)

image.png

5、合并数据

df <-  merge(as.data.frame(df_res),
                  as.data.frame(counts(df_dds_res,normalize=TRUE)),
                  by="row.names",sort=FALSE)

image.png

设置阈值并对数据进行分类

#去除空值
df1<-df[!is.na(df$padj),]
#数据分类——根据其中log2FoldChange、padj指标对OTU进行分类
df1$SG<-as.factor(ifelse(df1$padj<0.05&abs(df1$log2FoldChange)>=2,"Y","N"))
#需要注释的差异OTU标签——差异倍数大于6的进行注释，大家可自行根据需要设置
df1$label<-ifelse(df1$padj<0.05&abs(df1$log2FoldChange)>=6,"Y","N")
df1$label<-ifelse(df1$label == 'Y', as.character(df1$Row.names), '')
#图例标签
df1$SG <- factor(df1$SG, levels = c('Y', 'N'), labels = c('differences','No difference'))

image.png

使用火山图进行展示

p <- ggplot(df1, aes(log2FoldChange, -log10(padj))) +
  geom_point(aes(color = SG),alpha=0.6, size=2)+
  theme_bw()+
  theme(legend.title = element_blank(),
        panel.grid = element_blank()) +
  scale_x_continuous(breaks=seq(-10,10, 2))+
  geom_vline(xintercept = c(-2, 2), lty=3,color = 'black', lwd=0.5) + 
  geom_hline(yintercept = -log10(0.05), lty=3,color = 'black', lwd=0.5) +
  scale_color_manual(values = c( 'red','grey'))+
  labs(title="volcanoplot",
       x = 'log2 fold change',
       y = '-log10 pvalue')+
  geom_text_repel(aes(x = log2FoldChange,
                      y = -log10(padj),          
                      label=label),                       
                  max.overlaps = 1000,
                  size=3,
                  box.padding=unit(0.8,'lines'),
                  point.padding=unit(0.8, 'lines'),
                  segment.color='black',
                  show.legend=FALSE)
p

image.png

参考：https://blog.csdn.net/qq_42458954/article/details/104078845