一、举例回顾

本节使用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分析和绘图。