14.GO富集分析R语言代码及5种图的绘制
一、举例回顾
本节使用GSE1009数据集，已经用limma包对数据集中的样本进行差异分析，现对差异基因(DEGs)做GO富集分析。
GSE1009数据集介绍：
样本量：共6个样本，其中后3个为糖尿病肾病(DN)肾小球样本，前3个为正常肾小球样本。
使用芯片：Affymetrix Human Genome U95 Version 2 Array。
平台：GPL8300。
DEGs:共有66个DEGs（diffsig）,22个上调(diffup),44个下调（diffDown）(详见上两章).
二、需要准备的文件：
包含差异基因名字+logFC值的文本文件，命名为symbol(下面有介绍详细做法。)
三、具体做法：
1. 整理symbol文件
在上一节中diffSig变量，用write.table函数保存为包含差异基因名字+logFC值等的文件，命名为diff。
将diffSig变量保存哦
保存的方法如上，如果上次已经保存了的，直接调用
保存好diff文件，是下面这样的：
保存好的diff文件
手动给基因名加上列名，为gene,如下：
加上第一列名gene
重点：新建一个文本文件，取名为symbol,将第一列gene和第二列logFC复制symbol.
symbol
2.ID转换(将基因名转换为entrezID
setwd("D:\\Rfile")
rm(list = ls())
options(stringsAsFactors=F)
#老规矩，先设置工作目录。
library("clusterProfiler")
library("org.Hs.eg.db")
library("enrichplot")
library("ggplot2")
#加载这些包，加载之前记得先安装，已经安装过的复制代码直接调用。
rt=read.table("symbol.txt",sep="\t",check.names=F,header=T)
#读取symbol文件,并赋值给rt
genes=as.vector(rt[,1])
#取rt的第一列，即基因名字，将其转换为向量，并赋值给genes变量
entrezIDs <- mget(genes, org.Hs.egSYMBOL2EG, ifnotfound=NA)
#找出基因对应的id，未找到的赋值为NA
entrezIDs <- as.character(entrezIDs)
out=cbind(rt,entrezID=entrezIDs)
#将基因ID转换为entrezIDs
write.table(out,file="id.txt",sep="\t",quote=F,row.names=F) #输出结果，结果为id文本文档
3.GO分析及绘图
##读取ID转换后文件
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变量
##GO分析
#GO富集分析
kk <- enrichGO(gene = gene,
OrgDb = org.Hs.eg.db,
pvalueCutoff =0.05,
qvalueCutoff = 0.05,
ont="all",
readable =T)
#p值和q值可以自己设置，我们这里都设置成0.05。这一步需要一点时间，耐心等待。
write.table(kk,file="GO.txt",sep="\t",quote=F,row.names = F)
#保存GO富集结果
#1.GO柱状图
pdf(file="GO柱状图.pdf",width = 10,height = 8)
barplot(kk, drop = TRUE, showCategory =10,split="ONTOLOGY") + facet_grid(ONTOLOGY~., scale='free')
dev.off()
#pdf格式
png(file="GO柱状图.png",width = 800,height = 1000)
barplot(kk, drop = TRUE, showCategory =10,split="ONTOLOGY") + facet_grid(ONTOLOGY~., scale='free')
dev.off()
#保存为png格式，图片保存格式只举例这次，后面一样的道理，也可以保存为tiff格式等
#2.GO点图
pdf(file="GO点图.pdf",width = 10,height = 8)
dotplot(kk,showCategory = 10,split="ONTOLOGY") + facet_grid(ONTOLOGY~., scale='free')
dev.off()
##3.GO浮动气泡图
library(GOplot)
ego=read.table("GO.txt", header = T,sep="\t",check.names=F)
#读取GO富集结果文件
go=data.frame(Category = ego$ONTOLOGY,ID = ego$ID,Term = ego$Description, Genes = gsub("/", ", ", ego$geneID), adj_pval = ego$p.adjust)
#读取基因的logFC文件
id.fc <- read.table("id.txt", header = T,sep="\t",check.names=F)
genelist <- data.frame(ID = id.fc$gene, logFC = id.fc$logFC)
row.names(genelist)=genelist[,1]
circ <- circle_dat(go, genelist)
#绘制GO浮动气泡图
pdf(file="GO气泡图.pdf",width = 10,height = 8)
GOBubble(circ, labels = 3,table.legend =F)
dev.off()
#4.绘制GO圆圈图
pdf(file="GO圈图.pdf",width = 14,height = 6)
GOCircle(circ,rad1=2.5,rad2=3.5,label.size=4,nsub=10)
#rad1外圈的注释；nsub=10中10代表显示GO的数据，可修改
dev.off()
#5.绘制GO热图
termNum = 20 #限定term数目
geneNum = nrow(genelist) #限定基因数目
chord <- chord_dat(circ, genelist[1:geneNum,], go$Term[1:termNum])
pdf(file="GO热图.pdf",width = 11,height = 5)
GOHeat(chord, nlfc =1, fill.col = c('red', 'white', 'blue'))
dev.off()
GO分析和GO的几种图就完了，大家可根据自身情况选择一种图，下一章是KEGG分析和绘图。
