差异表达基因富集结果可视化

1. 导入差异表达基因

### 1. 导入差异表达基因

# if (!require("BiocManager", quietly = TRUE))
#   install.packages("BiocManager")
# 
# BiocManager::install("ReactomePA")

library(enrichplot)
library(clusterProfiler)
library(org.Hs.eg.db)
library(DOSE)
library(ReactomePA)
data(geneList)
gene <- names(geneList)[abs(geneList) > 2]
# gene:gene ids 向量，geneList：foldChanges 向量，names为基因ids

2. 差异表达基因富集分析

### 2. 差异表达基因富集分析
##ORA
res1 <- enrichPathway(gene)
# enrichPathway {ReactomePA}:Pathway Enrichment Analysis of a gene set.
res1 <- setReadable(res1, OrgDb = org.Hs.eg.db, keyType="ENTREZID")
# setReadable {DOSE}:mapping geneID to gene Symbol
res2 <- enrichGO(gene, OrgDb = org.Hs.eg.db, ont = "ALL")
res2 <- setReadable(res2, OrgDb = org.Hs.eg.db, keyType="ENTREZID")
##GSEA
res3 <- gseKEGG(geneList = geneList,organism= 'hsa')
res3 <- setReadable(res3, OrgDb = org.Hs.eg.db, keyType="ENTREZID")

3. 富集通路/基因集作图展示

### 3. 富集通路/基因集作图展示

########## ORA富集结果可视化#############
## 条形图展示，前10个
barplot(res1, showCategory=10)

##泡泡图
dotplot(res1,showCategory=20,title="enrichPathway Top20") 

##dotplot的棒棒糖图
library(ggplot2)
library(forcats)
ggplot(res1, showCategory = 20, 
       aes(GeneRatio, fct_reorder(Description, GeneRatio))) + 
  geom_segment(aes(xend=0, yend = Description)) +
  geom_point(aes(color=p.adjust, size = Count)) +
  scale_color_viridis_c(guide=guide_colorbar(reverse=TRUE)) +
  scale_size_continuous(range=c(2, 10)) +
  theme_minimal() + 
  ylab(NULL) 

## 展示感兴趣的通路
df = res1@result # 数据框，含所有富集的通路，基因等
choice_set = c("R-HSA-68886","R-HSA-69618", "R-HSA-69620")
# 找到Description对应的ID
index= which(rownames(res1@result) %in% choice_set)
Description=res1@result$Description[index]
barplot(res1, showCategory=Description) 
dotplot(res1,showCategory=Description) 

## 根据GO层次("BP","CC","MF"),分别展示
table(res2@result$ONTOLOGY)

barplot(res2, split = "ONTOLOGY") +
  facet_grid(ONTOLOGY~., scale = "free")

dotplot(res2, split = "ONTOLOGY") +
  facet_grid(ONTOLOGY~., scale = "free")

########## GSEA富集结果可视化#############
## 显示富集通路的Enrichment score及q值
library(ggplot2)
library(forcats) #  因子(factor)处理包
library(ggstance) #ggstance:实现常见图形的横向版本
ggplot(res3, aes(NES, fct_reorder(Description, NES), fill=qvalues), showCategory=20) + 
  geom_bar(stat='identity') + 
  scale_fill_continuous(low='red', high='blue', guide=guide_colorbar(reverse=TRUE)) + 
  theme_minimal() + ylab(NULL)

## 选择感兴趣的上下调通路，如下选择的是上下调的Top10（ES）
df = res3@result
# 可以跟据p.adj值过滤，去掉不显著的通路
df <- df[which(df$p.adjust>0.001),]

df = df[order(df$enrichmentScore),]
#设定分组：1--ES>0;   -1--ES<0
up = head(subset(df, enrichmentScore>0),10);up$group=1
down = tail(subset(df, enrichmentScore<0),10);down$group=-1
dat=rbind(up,down)
#dat$group = factor(dat$group);str(dat)
dat$pvalue = -log10(dat$pvalue)
dat$pvalue=dat$pvalue*dat$group 
dat=dat[order(dat$pvalue,decreasing = F),]
ggplot(dat, aes(x=reorder(Description,order(pvalue, decreasing = F)), y=pvalue, fill=group)) + 
  geom_bar(stat="identity", aes(fill=factor(group, levels = c(1,-1),  
                                            labels = c("ES>0","ES<0")))) + 
  xlab("Pathway names") +
  ylab("-log10(P.adj)") +
  coord_flip() + 
  scale_y_continuous(breaks=c(-4, -2, 0, 2, 4),
                     labels=c("4", "2", "0","2","4")) +
  scale_fill_manual(values = c("#ff6633","#34bfb5")) + 
  theme(plot.title = element_text(size = 15,hjust = 0.5),  
        axis.text = element_text(size = 12,face = 'bold'),
        panel.grid.major.x = element_blank(),
        panel.grid.minor.x = element_blank(),
        legend.title = element_blank()) +
  ggtitle("Pathway Enrichment") 

4. 富集通路与差异基因的网络关系图

### 4. 富集通路与差异基因的网络关系图
## 网络图
library(patchwork) #基于ggplot2的拼图包
# install.packages("ggnewscale")
library(ggnewscale) # ggnewscale实现在ggplot2中使用多种颜色和填充比例。
p1 = cnetplot(res1, foldChange=geneList,showCategory = 3)
p2 = cnetplot(res1, foldChange=geneList,showCategory = 3, circular = T)
p3 = cnetplot(res1, foldChange=geneList,showCategory = 3, colorEdge = T, node_label="gene")
p4 = cnetplot(res1, foldChange=geneList,showCategory = 3, layout = "gem")
(p1 + p2) / (p3 + p4)
# 选择你喜欢的layout
cnetplot(res3, foldChange=geneList,showCategory = 3)

## 热图的形式展示富集通路与差异基因
heatplot(res1, foldChange=geneList)
heatplot(res3, foldChange=geneList)

5. 富集通路之间的网络关系图

### 5. 富集通路之间的网络关系图

# emapplot {enrichplot}:Enrichment Map for enrichment result 
# of over-representation test or gene set enrichment analysis.

res1 <- enrichplot::pairwise_termsim(res1)
# pairwise_termsim {enrichplot}:Get the similarity matrix
emapplot(res1, showCategory = 15, layout="kk", cex_category=1.5,min_edge = 0.8) 

res3 <- enrichplot::pairwise_termsim(res3)
emapplot(res3, showCategory = 15, layout="kk", min_edge = 0.1, color = "NES") 

6. GSEA running score作图

### 6. GSEA running score作图
gseaplot(res3, geneSetID = 1, title = res3$Description[1])
gseaplot2(res3, geneSetID = 1:3, pvalue_table = T)
gseaplot2(res3,1:5,subplots = 1:2, base_size = 20,pvalue_table = T)

参考：

enrichplot||基因富集结果可视化解决方案 - 简书

RNAseq数据，下载GEO中的FPKM文件后该怎么下游分析