DESeq2筛选差异OTU及绘制火山图

1.DESeq包安装

install.packages("BiocManager")
library(BiocManager)
BiocManager::install("DESeq2")
library(DESeq2)

2.数据准备

共需要两份数据文件

a.OTU丰度表格,otutab.txt;(每一列为一个样本,每一行为一种OTU,交叉区域为每种OTU在各样本中的丰度,DESeq2计算时只能识别整数,不识别小数,所以请不要使用相对丰度表格)

b.样本分组信息表格,即metadata.txt。(行为样本名称,第二列为各样本对应的分组信息)

3.常规流程

a. 读入文件

otu_file <- read.delim('otutab.txt', row.names = 1, sep = '\t', stringsAsFactors = FALSE, check.names = FALSE)
group_file <- read.delim('metadata.txt', row.names = 1, sep = '\t')

b. 使用?? DESeq2参阅官方说明文档中的常规流程

Step 1,使用DESeqDataSet(),构建DESeq2存储read数的数据集;

Step 2,使用DESeq(),基于负二项式(Negative Binomial或称 Gamma-Poisson)分布,对数据(本示例中为OTU丰度数据)进行DESeq差异分析;

Step 3,使用results(),用于在DESeq分析中提取结果。

4.png

library(DESeq2)
dds <- DESeqDataSetFromMatrix(countData = otu_file, colData = group_file, design = ~group)#构建 DESeqDataSet 对象  
dds <- DESeq(dds) #差异分析
suppressMessages(dds)

c.结果查看

接下来,我们要查看总体结果,并根据p-value进行重新排序。利用summary命令统计显示一共多少个OTU上调和下调

res <- results(dds, contrast=c('group', 'treat', 'control'))#提取分析结果
res = res[order(res$pvalue),]
res #查看结果
summary(res)  #简要统计结果
table(res$padj<0.05) #查看fdr校正后的P<0.05的个数

d.提取差异OTU并进行注释

获取padj(p值经过多重校验校正后的值)小于0.05,表达倍数取以2为对数后大于1或者小于-1的差异OTU

diff_OTU_deseq2 <-subset(res, padj < 0.05 & abs(log2FoldChange) > 1)
#或
# diff_OTU_deseq2 <-subset(res,padj < 0.05 & (log2FoldChange > 1 | log2FoldChange < -1))
dim(diff_OTU_deseq2)
head(diff_OTU_deseq2)
write.csv(diff_OTU_deseq2,file= "DEG_treat_vs_control.csv")

e.合并结果

resdata <-  merge(as.data.frame(res),as.data.frame(counts(dds,normalize=TRUE)),by="row.names",sort=FALSE)
write.table(resdata,file= "DESeq2_Group_genus.txt",sep="\t",quote=F,row.names=F)

4.火山图绘制

首先在上文得到的转化为数据框类型的DESeq2差异分析结果“resdata”的基础上,根据其中log2FoldChange、padj这两个重要的指标,对OTU进行分类.

for (i in 1:nrow(resdata)) {
    if (abs(resdata[i,'log2FoldChange']) >= 0.5) resdata[i,'select_change'] <- 'y' else resdata[i,'select_change'] <- 'n'
    if (resdata[i,'padj'] %in% NA | abs(resdata[i,'padj']) >= 0.05) resdata[i,'select_pvalue'] <- 'n' else resdata[i,'select_pvalue'] <- 'y'
    resdata[i,'select'] <- paste(resdata[i,'select_change'], resdata[i,'select_pvalue'], sep = '')
}

对于每个OTU,若log2FoldChange < 0.5,即该OTU在两组间的丰度差异倍数低于2(默认情况下,将差异倍数2判定为分界点),则在“resdata”的“select_change”列标记为“n”,反之为“y”;若padj >= 0.05,即校正后的显著性p值未低于0.05的显著性水平,则在“resdata”的“select_pvalue”列标记为“n”,反之为“y”。

同时合并“select_change”和“select_pvalue”的结果,可得到“nn”(p >= 0.05,FC < 2)、“ny”(p < 0.05,FC < 2)、“yn”(p >= 0.05,FC >= 2)、“yy”(p < 0.05,FC >= 2)四种组合。

library(ggplot2)
 
##ggplot2 差异火山图
resdata$select <- factor(resdata$select, levels = c('nn', 'ny', 'yn', 'yy'), labels = c('p >= 0.05, FC < 2', 'p < 0.05, FC < 2', 'p >= 0.05, FC >= 2', 'p < 0.05, FC >= 2'))
 
#纵轴为显著性 p 值
volcano_plot_pvalue <- ggplot(resdata, aes(log2FoldChange, -log(padj, 10))) +
geom_point(aes(color = select), alpha = 0.6) +
scale_color_manual(values = c('gray30', 'green4', 'red2', 'blue2')) +
theme(panel.grid = element_blank(), panel.background = element_rect(color = 'black', fill = 'transparent')) +
theme(legend.title = element_blank(), legend.key = element_rect(fill = 'transparent'), legend.background = element_rect(fill = 'transparent')) +
geom_vline(xintercept = c(-0.5, 0.5), color = 'gray', size = 0.5) + 
geom_hline(yintercept = -log(0.05, 10), color = 'gray', size = 0.5) +
labs(x = 'log2 Fold Change', y = '-log10 p-value')
 
ggsave('volcano_plot_pvalue.png', volcano_plot_pvalue, width = 6, height = 5)

10.png

第二种绘图方式

 可在火山图上添加标签,以我的数据为例。ggplot中,可使用ggreprel包中的geom_text_repel()函数进行标签的添加。

with(resdata, plot(log2FoldChange, -log10(pvalue), pch=20, main="Volcano plot", xlim=c(-2.5,4)))

# Add colored points: red if padj<0.05, orange of log2FC>1, green if both)
with(subset(resdata, padj<.05 ), points(log2FoldChange, -log10(pvalue), pch=20, col="red"))
with(subset(resdata, abs(log2FoldChange)>1), points(log2FoldChange, -log10(pvalue), pch=20, col="orange"))
with(subset(resdata, padj<.05 & abs(log2FoldChange)>1), points(log2FoldChange, -log10(pvalue), pch=20, col="green"))

# Label points with the textxy function from the calibrate plot
library(calibrate)
with(subset(resdata, padj<.05 & abs(log2FoldChange)>1.5), textxy(log2FoldChange, -log10(pvalue), labs=ID, cex=.8))

绘图方式三

# rm(list = ls()) 此步骤为清空R环境,非必要不使用
library(ggplot2)
resdata$color <- ifelse(resdata$padj<0.05 & abs(resdata$log2FoldChange)>= 1,ifelse(resdata$log2FoldChange > 1,'red','blue'),'gray')
color <- c(red = "red",gray = "gray",blue = "blue")

p <- ggplot(resdata, aes(log2FoldChange, -log10(padj), col = color)) +
  geom_point() +
  theme_bw() +
  scale_color_manual(values = color) +
  labs(x="log2 (fold change)",y="-log10 (q-value)") +
  geom_hline(yintercept = -log10(0.05), lty=4,col="grey",lwd=0.6) +
  geom_vline(xintercept = c(-1, 1), lty=4,col="grey",lwd=0.6) +
  theme(legend.position = "none",
        panel.grid=element_blank(),
        axis.title = element_text(size = 16),
        axis.text = element_text(size = 14))
p

参考来源:科学网—lyao222lll的个人资料

RNA-seq(7): DEseq2筛选差异表达基因并注释(bioMart) - 简书

关于EdgeR分析差异OTU详见博文https://www.jianshu.com/p/e6bc48d8573e

如需科研绘图,可访问我的淘宝店铺:店铺名:R语言绘图,可提供专业绘图帮助,直接SCI出版。

  • 5
    点赞
  • 63
    收藏
    觉得还不错? 一键收藏
  • 12
    评论
要用R绘制单个样本OTU注释圈,需要先安装和加载相关的R包,如ggplot2、reshape2和RColorBrewer。然后,可以按照以下步骤和代码进行绘制: 1. 准备数据 首先,需要准备好OTU注释数据,其中包括OTU ID和注释信息。可以将这些数据保存在一个.csv文件中,并使用read.csv()函数将其读入到R中。 2. 筛选数据 根据需要,可以对注释信息进行筛选,比如只保留某些层级的注释信息。可以使用dplyr包中的filter()函数实现该功能。 3. 数据整理 将筛选后的注释信息与OTU ID进行合并,并使用melt()函数将数据整理成适合绘的格式。 4. 绘 使用ggplot2包中的geom_point()函数绘制圆形,根据注释信息将圆形填充颜色,并使用coord_polar()函数将圆形排列成圆环状。 以下是示例代码: ```R #加载必要的包 library(ggplot2) library(reshape2) library(RColorBrewer) #读入OTU注释数据 otu_anno <- read.csv("otu_anno.csv", header = TRUE) #筛选注释信息 otu_anno <- filter(otu_anno, level == "phylum" | level == "class") #整理数据 otu_anno_melt <- melt(otu_anno, id.vars = "OTU_ID") #绘 ggplot(data = otu_anno_melt, aes(x = variable, y = 1, fill = value)) + geom_point(size = 10, shape = 21) + scale_fill_brewer(palette = "Set1") + theme_void() + coord_polar(theta = "x") ``` 以上代码会绘制出一个圆环状的形,圆环上的每个圆形代表一个OTU,圆形的颜色代表该OTU的注释信息。

“相关推荐”对你有帮助么?

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论 12
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值