R优雅的进行代谢组KEGG富集分析

R优雅的进行代谢组KEGG富集分析 by R语言数据分析指南

欢迎关注R语言数据分析指南

❝
本节来介绍如何对代谢组数据进行KEGG富集分析;数据及代码已上传到「VIP群」已经加群的观众老爷可以直接获取，小编的VIP群目前已经上传「公众号文档数据+代码约160余篇」，有需要加群的欢迎「先点击发消息」之后公众号右下角添加微信「付费99元」后邀请进群
❞

往期内容

安装R包

library(tidyverse)
library(magrittr)
library(clusterProfiler)

导入KEGG数据库注释文件

keggannotation <- read_tsv("pathway",col_names = F) %>% 
  left_join(.,read_tsv('map.txt',col_names = F),by="X1") %>% 
  select(-1) %>% set_colnames(c("pathway","ID")) %>% 
  mutate(across("ID",str_replace,"cpd:","")) %>% select(2,1) %>% 
  arrange(ID)

全部代谢物编号

allkeggid <- read_tsv("meta_intensity_neg.anno.xls") %>% select(KEGG) %>% 
  bind_rows(.,read_tsv("meta_intensity_pos.anno.xls") %>% select(KEGG)) %>% 
  arrange() %>% filter(KEGG !="_") %>% set_colnames(c("ID"))

差异代谢物

diffkeggID <- read_tsv("diff.xls") %>% select(KEGG) %>% 
  arrange() %>% filter(KEGG !="_") %>% set_colnames(c("ID"))

数据整合

total <- right_join(keggannotation,allkeggid,by="ID") %>% select(2,1)

富集分析

x <- clusterProfiler::enricher(gene = diffkeggID$ID,TERM2GENE = total,minGSSize = 1,pvalueCutoff = 1,qvalueCutoff = 1)

结果导出

write.csv(as.data.frame(x@result) %>% select(-1,-2),file="KEGG_enrichment_result.csv")

数据清洗

df <- read_csv("KEGG_enrichment_result.csv") %>%
  dplyr::rename("Description"="...1") %>% 
  arrange(desc(Count)) %>%
  select(1,2,3,4,8) %>% 
  separate(`GeneRatio`,into=c("A","B"),sep="/") %>%
  mutate(A=as.numeric(A),B=as.numeric(B)) %>%
  mutate(count=A/B) %>% head(30) %>% arrange(Count)

定义因子

df$Description <- factor(df$Description,levels = c(df$Description %>% as.data.frame() %>% pull()))

数据可视化

df %>% ggplot(aes(count,Description))+
  geom_point(aes(size=Count,color=pvalue,fill=pvalue),pch=21)+
  scale_color_gradientn(colours = (rev(RColorBrewer::brewer.pal(11,"RdBu"))))+
  scale_fill_gradientn(colours =(rev(RColorBrewer::brewer.pal(11,"RdBu"))))+
  guides(size=guide_legend(title="Count"))+
  labs(x=NULL,y=NULL)+
  theme(axis.title = element_blank(),
        axis.text.x=element_text(color="black",angle =0,hjust=0.5,vjust=0.5, margin = margin(b =5)),
        axis.text.y=element_text(color="black",angle =0,hjust=1,vjust=0.5),
        panel.background = element_rect(fill = NA,color = NA),
        panel.grid.minor= element_line(size=0.2,color="#e5e5e5"),
        panel.grid.major = element_line(size=0.2,color="#e5e5e5"),
        panel.border = element_rect(fill=NA,color="black",size=1,linetype="solid"),
        legend.key=element_blank(),
        legend.title = element_text(color="black",size=9),
        legend.text = element_text(color="black",size=8),
        legend.spacing.x=unit(0.1,'cm'),
        legend.key.width=unit(0.5,'cm'),
        legend.key.height=unit(0.5,'cm'),
        legend.background=element_blank(),
        legend.box="horizontal",
        legend.box.background = element_rect(color="black"),
        legend.position = c(1,0),legend.justification = c(1,0))+
  scale_y_discrete(labels = function(y) str_wrap(y,width=30))

数据获取

❝
本节介绍到此结束，需要获取之前绘图文档数据及代码的朋友欢迎加入小编的2022年VIP交流群，付费99元，群内会同步上传公众号文档代码；目前已上传2021-2022约160篇文档代码；添加小编微信时，请备注单位-方向-姓名以及来意以便高效处理, 本篇文档数据可从小编的「vip交流群」中获取，为了满足各位观众老爷的需求小编开通了同名淘宝店，提供数据分析服务欢迎各位观众老爷光顾
❞