跟小新老师学转录组的第五天

rm(list = ls())
options(stringsAsFactors = F)
library(clusterProfiler)
library(org.Hs.eg.db)
library(GSEABase)
library(ggplot2)
library(tidyverse)
packageVersion("clusterProfiler")

## Error in download.KEGG.Path(species)
# https://github.com/YuLab-SMU/clusterProfiler/pull/471
getOption("clusterProfiler.download.method")
#R.utils::setOption("clusterProfiler.download.method",'auto')
options(clusterProfiler.download.method = "wininet")
#options(clusterProfiler.download.method = "wget")
getOption("clusterProfiler.download.method")



# 读取差异分析结果
load(file = "data/Step03-edgeR_nrDEG.Rdata")
ls()

# 提取所有差异表达的基因名
DEG <- DEG_edgeR_symbol[DEG_edgeR_symbol$regulated!="normal", 2]
DEG <- as.character(na.omit(DEG))
head(DEG)

## ===GO数据库, 输出所有结果，后续可根据pvalue挑选结果
ego_CC <- enrichGO(gene=DEG, OrgDb= 'org.Hs.eg.db', keyType='SYMBOL', ont="CC", pvalueCutoff= 1,qvalueCutoff= 1)
ego_MF <- enrichGO(gene=DEG, OrgDb= 'org.Hs.eg.db', keyType='SYMBOL', ont="MF", pvalueCutoff= 1,qvalueCutoff= 1)
ego_BP <- enrichGO(gene=DEG, OrgDb= 'org.Hs.eg.db', keyType='SYMBOL', ont="BP", pvalueCutoff= 1,qvalueCutoff= 1)

p_BP <- barplot(ego_BP,showCategory = 10, label_format = 100) + ggtitle("Biological process")
p_CC <- barplot(ego_CC,showCategory = 10, label_format = 100) + ggtitle("Cellular component")
p_MF <- barplot(ego_MF,showCategory = 10, label_format = 100) + ggtitle("Molecular function")
plotc <- p_BP/p_CC/p_MF
plotc
ggsave('result/6.enrichGO.png', plotc, width = 10,height = 16)


## 使用ggplot2绘制
# 对富集结果使用qvalue从小到大排列，取top10
data <- ego_CC@result %>% top_n(n = 10, wt = -(qvalue))
colnames(data)
p <- ggplot(data = data, aes(x=-log10(qvalue), y=reorder(Description,-log10(qvalue)) )) + 
  geom_bar(stat = 'identity', fill="OliveDrab4", position="dodge", width =0.5) + 
  scale_x_continuous(expand = c(0,0)) + # 调整柱子底部与y轴紧贴
  xlab(label = "") + ylab(label = "") + # 设置x，y坐标轴标题
  theme_classic() +  # 设置主题，只保留x，y轴
  ggtitle(label = "-LOG(q-value)") + # 设置标题
  theme(plot.title = element_text(size = 12,color="black",hjust = 0.5), # 标题居中
        axis.title = element_text(size = 15,color ="black"), # 标题字体大小
        axis.text = element_text(size= 12,color = "black"))  # 坐标轴字体大小
  

p


ego_BP <- data.frame(ego_BP)
ego_CC <- data.frame(ego_CC)
ego_MF <- data.frame(ego_MF)
write.csv(ego_BP,'result/6.enrichGO_BP.csv')
write.csv(ego_CC,'result/6.enrichGO_CC.csv')
write.csv(ego_MF,'result/6.enrichGO_MF.csv')



## === KEGG
genelist <- bitr(gene=DEG, fromType="SYMBOL", toType="ENTREZID", OrgDb='org.Hs.eg.db')
genelist <- pull(genelist,ENTREZID)               
ekegg <- enrichKEGG(gene = genelist, organism = 'hsa', pvalueCutoff = 1, qvalueCutoff = 1)
p1 <- barplot(ekegg, showCategory=10, label_format = 100)
p2 <- dotplot(ekegg, showCategory=10, label_format = 100)
plotc = p1/p2
plotc
ggsave('result/6.enrichKEGG.png', plot = plotc, width = 8, height = 10)

ekegg <- data.frame(ekegg)
write.csv(ekegg,'result/6.enrichKEGG.csv')



## === 其他数据库通路
geneset <- read.gmt("data/MsigDB/v7.4/h.all.v7.4.symbols.gmt")
as.data.frame(table(geneset$term))
geneset$term <- gsub(pattern = "HALLMARK_","", geneset$term)
geneset$term <- str_to_title(geneset$term)

my_path <- enricher(gene=DEG, pvalueCutoff = 1, qvalueCutoff = 1, TERM2GENE=geneset)
p1 <- barplot(my_path, showCategory=15,color = "pvalue", label_format = 100)
p1
ggsave("result/6.enrich_HALLMARK.png", plot = p1, width = 8, height = 7)
  
my_path <- data.frame(my_path)
write.csv(my_path,"result/6.enrich_HALLMARK.csv")
  

# 清空当前环境变量
rm(list = ls())
options(stringsAsFactors = F)

# 加载包
library(GSEABase)
library(clusterProfiler)
library(enrichplot)
library(ggplot2)
library(stats)

# 加载数据
load("data/Step03-edgeR_nrDEG.Rdata")
DEG <- DEG_edgeR_symbol

## 构造GSEA分析数据
# 去掉没有配对上symbol的行
DEG <- DEG[!is.na(DEG$SYMBOL),]

# 去掉重复行
DEG <- DEG[!duplicated(DEG$SYMBOL),]

geneList <- DEG$logFC
names(geneList) <- DEG$SYMBOL
head(geneList)
geneList <- sort(geneList,decreasing = T)
head(geneList)
tail(geneList)


# 选择gmt文件
geneset <- read.gmt("data/MsigDB/h.all.v7.5.1.symbols.gmt")

# 运行,输出全部结果
egmt <- GSEA(geneList, TERM2GENE=geneset, pvalueCutoff = 1, verbose = T)

#出点图
dotplot(egmt,label_format = 100)

#按p值出点图
dotplot(egmt,color="pvalue",label_format = 100)   

# 单个通路图
# 按照通路名
gseaplot2(egmt, "HALLMARK_ADIPOGENESIS",  
          title = "HALLMARK_ADIPOGENESIS")
# 按照行数
gseaplot2(egmt, 5, color="red", pvalue_table = T)

#按第1到第6个出图，不显示p值
gseaplot2(egmt, 1:6, color="red") 


# 保存结果
go_gsea <- as.data.frame(egmt@result)
write.csv(go_gsea,"result/6.gsea_go_fc.csv",row.names = F)