一、舉例回顧
本節(jié)所使用GSE1009數(shù)據(jù)集，已經(jīng)用limma包進(jìn)行差異分析，現(xiàn)對DEGs做GO富集分析。
GSE1009數(shù)據(jù)集介紹：??
樣本量：共6個樣本，其中后3個為糖尿病腎病(DN)腎小球樣本，前3個為正常腎小球樣本。
使用芯片：Affymetrix Human Genome U95 Version 2 Array。
平臺：GPL8300。
DEGs:共有66個DEGs（diffsig）,22個上調(diào)(diffup),44個上調(diào)（diffDown）(詳見上兩章).
二、需要準(zhǔn)備的文件：
包含差異基因名字+logFC值的文本文件，命名為symbol(上一章有介紹詳細(xì)做法。)
[if !supportLists]三、[endif]具體步驟：
[if !supportLists]1.?[endif]ID轉(zhuǎn)換（同上一章）
setwd("D:\\Rfile")
rm(list = ls())
options(stringsAsFactors=F)
#老規(guī)矩，先設(shè)置工作目錄。
library("clusterProfiler")
library("org.Hs.eg.db")
library("enrichplot")
library("ggplot2")
#加載這些包，加載之前記得先安裝，已經(jīng)安裝過的復(fù)制代碼直接調(diào)用。
rt=read.table("symbol.txt",sep="\t",check.names=F,header=T) ???
#讀取symbol文件,并賦值給rt
genes=as.vector(rt[,1])
#取rt的第一列，即基因名字，將其轉(zhuǎn)換為向量，并賦值給genes變量
entrezIDs <- mget(genes, org.Hs.egSYMBOL2EG, ifnotfound=NA) ???
#找出基因?qū)?yīng)的id，未找到的賦值為NA
entrezIDs <- as.character(entrezIDs)
out=cbind(rt,entrezID=entrezIDs)
#將基因ID轉(zhuǎn)換為entrezIDs
write.table(out,file="id.txt",sep="\t",quote=F,row.names=F) ???#輸出結(jié)果，結(jié)果為id文本文檔
##讀取ID轉(zhuǎn)換后文件
rt=read.table("id.txt",sep="\t",header=T,check.names=F) ??????????#讀取id.txt文件
rt=rt[is.na(rt[,"entrezID"])==F,] ??????????????????????????????#去除基因id為NA的基因
gene=rt$entrezID#取entrezID賦值給gene變量
2.KEGG
kk2 <- enrichKEGG(gene = gene, organism = "hsa", pvalueCutoff =0.05, qvalueCutoff =0.05) ?
#富集分析
write.table(kk2,file="KEGGId.txt",sep="\t",quote=F,row.names = F) ?????????????????????????#保存富集結(jié)果
#KEGG柱狀圖
pdf(file="KEGG柱狀圖.pdf",width = 10,height = 7)
barplot(kk2, drop = TRUE, showCategory = 30)
dev.off()
#點圖
pdf(file="KEGG點圖.pdf",width = 10,height = 7)
dotplot(kk2, showCategory = 30)
dev.off()
#（因為我這個數(shù)據(jù)集做出來的結(jié)果不好，只有兩條通路，就不繼續(xù)做其他圖了，需要做其他圖的，代碼如下）
##KEGG氣泡圖
library(GOplot)
ego2=read.table("KEGGId.txt", header = T,sep="\t",check.names=F) ??????????#讀取kegg富集結(jié)果文件
go2=data.frame(Category = "ALL",ID = ego2$ID,Term = ego2$Description, Genes = gsub("/", ", ", ego2$geneID), adj_pval = ego2$p.adjust)
#讀取基因的logFC文件
id.fc2 <- read.table("id.txt", header = T,sep="\t",check.names=F)
genelist2 <- data.frame(ID = id.fc2$gene, logFC = id.fc2$logFC)
row.names(genelist2)=genelist2[,1]
circ2 <- circle_dat(go2, genelist2)
#繪制KEGG氣泡圖
pdf(file="KEGG氣泡圖.pdf",width = 10,height = 8)
GOBubble(circ2, labels = 3,table.legend =F)
dev.off()
#繪制KEGG圈圖
pdf(file="KEGG圈圖.pdf",width = 15,height = 6)
GOCircle(circ2,rad1=2.5,rad2=3.5,label.size=4,nsub=10) ???????????#nsub=10中10代表顯示KEGG的數(shù)據(jù)，可修改
dev.off()
#繪制KEGG熱圖
termNum =20 ????????????????????????????????????#限定term數(shù)目
geneNum = nrow(genelist2) ????????????????????????#限定基因數(shù)目
chord2 <- chord_dat(circ2, genelist2[1:geneNum,], go2$Term[1:termNum])
pdf(file="KEGGHeat.pdf",width = 9,height = 5)
GOHeat(chord2, nlfc =1, fill.col = c('red', 'white', 'blue'))
dev.off()
KEGG通路富集分析就完了，下一章是一般醫(yī)學(xué)專業(yè)才需要的DO分析。
