DAVID富集分析
DAVID富集支持多種基因名類型破托，本次以genesymbol為例
DAVID：https://david.ncifcrf.gov/
進入之后跟著以下步驟操作即可
選擇相應(yīng)物種雄驹，以人和小鼠為例田巴，載點擊Submit
一般是進行GO（BP\CC\MF）、KEGG分析龙屉，注意選擇
右擊Download File -> 鏈接另存為嗡善，即可將富集分析結(jié)果保存至本地文檔
```
R語言可視化
setwd("G:/富集分析")? #設(shè)置工作路徑
#此處以本次分析的上調(diào)基因為例，可以分別可視化GO和KEGG的結(jié)果褐荷，本次以GO_BP和KEGG為例
#上調(diào)富集分析繪圖
DA_UP <- read.table("DAVID富集結(jié)果.txt",sep = "\t",header = T)
#select Pvalue<0.05? 將P<0.05認為有統(tǒng)計學(xué)意義
DA_UP <- subset(DA_UP,DA_UP$PValue<0.05)
DA_UP <- subset(DA_UP,DA_UP$Count > 3)
DA_UP$PValue <- -(log10(DA_UP$PValue))
UP_BP <- subset(DA_UP,DA_UP$Category == "GOTERM_BP_DIRECT")
UP_KEGG <- subset(DA_UP,DA_UP$Category == "KEGG_PATHWAY")
UP_BP$Category <- sub("GOTERM_BP_DIRECT","GOBP_PATHWAY",UP_BP$Category)
#保留繪圖所需的信息
library(tidyverse)
UP_KEGG <- UP_KEGG[,c("Category","Term","Count","PValue")]
#修改通路類別名稱格式，使結(jié)果更好看
UP_KEGG <- separate(data = UP_KEGG, col = Term, into = c("ID", "Term"), sep = ":")
UP_KEGG <- UP_KEGG[,-2]
UP_BP <- UP_BP[,c("Category","Term","Count","PValue")]
UP_BP <- separate(data = UP_BP, col = Term, into = c("ID", "Term"), sep = "~")
UP_BP <- UP_BP[,-2]
UP <- rbind.data.frame(UP_KEGG,UP_BP)
DOWN <- rbind.data.frame(DOWN_KEGG,DOWN_BP)
write.table(UP,file = "UP_enrichment.txt",sep = "\t",quote = F,row.names = F)
write.table(DOWN,file = "DOWN_enrichment.txt",sep = "\t",quote = F,row.names = F)
##繪制柱狀圖
#使用reorder對柱狀圖排序
#想讓y軸從大到小變化嘹悼，只需在Change前面加個負號
##########################################################
##上下調(diào)差異基因富集圖
########################################################
setwd("G:/科研/畢設(shè)/富集分析")
UP <- read.table("UP_enrichment.txt",header = T,sep = "\t")
DOWN <- read.table("DOWN_enrichment.txt",header = T,sep = "\t")
#因為富集結(jié)果過多叛甫，因此選擇與研究內(nèi)容相關(guān)需要的通路
UP_END <- UP[c(1,2,4,5,6,8,9,13,17,19,20,21,24,28,32,46,49,51,53,62),]
DOWN_END <- DOWN[c(1,4,6,7,10,13,14,17,20,21,31,33,34,36,40,44,45,53,59,71),]
library(ggplot2)
a <- UP_END
b <- "UP_END"
pngname <- paste(b,".png",sep = "")
pdfname <- paste(b,".pdf",sep = "")
p1 <- ggplot(data = a,
? ? ? ? ? ? aes(x = reorder(Term, PValue),?
? ? ? ? ? ? ? ? y = PValue,
? ? ? ? ? ? ? ? fill = Category)) +?
? geom_bar(stat = "identity",
? ? ? ? ? width = 0.8,
? ? ? ? ? position = position_dodge(width = 0.9))? +
? geom_text(aes(label = Count),
? ? ? ? ? ? vjust = 0.4, hjust = -0.2) +
? facet_grid(Category~., scales = 'free_y') +
? theme(panel.grid.major.x = element_line(colour = "white"),
? ? ? ? panel.background = element_blank(),
? ? ? ? axis.line.y = element_blank(),
? ? ? ? axis.line.x = element_line(color = "black"),
? ? ? ? axis.title.y = element_blank(),
? ? ? ? strip.text.y = element_blank()) +
? #scale_x_discrete(position = "top") +
? ggtitle("UP_GENE")+ #添加標(biāo)題
? ylab("-log10PValue")+
? #geom_text(aes(label = y))+
? coord_flip()
p1
ggsave("UP_gene.png",p1, units="in", dpi=300, width=6, height=6, device="png")
```
