【富集分析】GO & KEGG

RNASeq

第二章 富集分析之GO和KEGG

---

---

前言

通过差异分析后可以对差异基因进行富集分析，常用的方式就是：

• GO（Gene Ontology）
  一个国际标准化的基因功能分类体系，提供了一套动态更新的标准词汇表 (controlled vocabulary)来全面描述生物体中基因和基因产物的属性。分别为生物学过程（biological process BP）、 细胞组分（cellular component CC）和分子功能（molecular function MF)。

• KEGG （Kyoto Encyclopedia of Genes and Genomes）
  系统分析基因产物在细胞中代谢途径的数据库，是一种最常用的代谢通路分析。

---

一、环境配置

代码如下（示例）：

Sys.setenv(language = "en") # 英文环境
options(stringsAsFactors = F) # 默认不转换为因子类型
rm(list = ls()) # 清空变量

二、加载包

代码如下（示例）：

library(DESeq2) 
library(dplyr)
library(clusterProfiler)
library(tibble)
library(AnnotationDbi)
library(org.Hs.eg.db)
library(msigdbr)
library(enrichplot)

三、差异分析

在进行富集分析前，需要先获得差异基因。RNA-seq的count数据有三个包可供使用来进行差异分析，上一篇有介绍，下面使用DESeq2包进行分析。

代码如下（示例）：

1.准备数据

# 读取数据
cds <- read.csv(file = "data/airway_scaledcounts.csv",row.names = 1) # 表达矩阵
colData <- read.csv(file = "data/airway_metadata.csv",row.names = 1) # 分组信息
annoData <- read.csv("data/annotables_grch38.csv") %>% 
	select(ensgene,entrez,symbol) # 注释信息

# 表达矩阵和分组信息匹配
if(!all(rownames(colData) == colnames(cds))){
  cds <- cds[colnames(cds),]
}

# 分组
colData$dex <- factor(colData$dex,levels = c("control","treated"))

2.分析

# 创建DESeqDataSet
dds <- DESeqDataSetFromMatrix(countData = cds,
                              colData = colData,
                              design = ~ dex)
                              # tidy = T)  # 在没有行名的情况下可以使用

# 过滤
# counts(dds) 提取原始表达矩阵
keep <- rowSums(counts(dds)) > 1  # 行和大于1的Ensembl Gene ID
dds <- dds[keep,]

# 分析
dds <- DESeq(dds)
resultsNames(dds) # 查看dds的结果名，便于指定想要提取的列
degDeSeq2 <- results(dds,name = "dex_treated_vs_control") %>% 
  as.data.frame() %>% 
  na.omit()

3. ID转换

本来是将打算将行名转换为entrez id，但是结果报错了
Error in .rowNamesDF<-(x, value = value) : missing values in 'row.names' are not allowed
后来换成symbol就没问题了，行名应该是不可以为数字。

# 通过读取的数据将ensgene转换成symbol
degDeSeq2$ensgene <- rownames(degDeSeq2)
degSymbol <- inner_join(degDeSeq2,annoData,by = "ensgene") %>%  # 通过ensgene合并
  dplyr::select(- c(ensgene,entrez)) %>% 
  distinct(symbol,.keep_all = T) %>% 
  column_to_rownames(var = "symbol")

4.差异基因筛选

# 设置差异基因筛选条件
tPValue <- 0.05
tLogFc <- 2

# 上调
pUp <- degSymbol$padj < tPValue 
upLogFC <- degSymbol$log2FoldChange > tLogFc
up <- pUp & upLogFC

# 下调
downLogFC <- degSymbol$log2FoldChange < -tLogFc
down <- pUp & downLogFC

# 合并上调和下调
degUp <- rownames(degSymbol)[up]
degDown <- rownames(degSymbol)[down]
deg <- c(degUp,degDown)

5.保存对象

# 保存对象 
save(degSymbol,degUp,degDown,deg,file = "outData/deg.Rdata")

三、GO分析

经过差异分析并且整理后的数据框行为基因名，列为差异倍数、p值等信息。通过对enrichgo()的查看可知参数gene的格式要求为entrez gene id的向量形式，所以：

1. 准备数据
   ！symbol转为为entrez gene id
2. go分析
3. 可视化

degSymbol数据框内容如下：

代码如下（示例）：

1.准备数据

load(file ="outData/deg.Rdata") # 提取deseq2分析结果
# 将symbol转entrez gene id 三种方式
## 方式一 AnnotationDb包 mapIds()
keytypes(org.Hs.eg.db)
degEntrezId <- mapIds(org.Hs.eg.db,
                      keys = deg,  # 需要转化的基因
                      column = "ENTREZID", # 转后的类型
                      keytype = "SYMBOL")  # 转换前类型
degEntrezId
class(degEntrezId)

## 方式二 AnnotationDb包 select() 返回数据框格式，column可以是多种类型
# degEntrezId2 <- select(org.Hs.eg.db,
#                       keys = deg,  # 需要转化的基因
#                       column = "ENTREZID", # 转后的类型
#                       keytype = "SYMBOL")  # 转换前类型

## 方式三 clusterprofiler包 bitr() 返回类型也是数据框
# degEntrezId3 <- bitr(deg,
#                      fromType = "SYMBOL",
#                      toType = "ENTREZID",
#                      OrgDb = org.Hs.eg.db)
# View(degEntrezId3)

2.go分析

goMF <- enrichGO(degEntrezId,
                 keyType = "ENTREZID",    # degEntrezIdl类型
                 OrgDb = "org.Hs.eg.db",  # 人的数据库
                 ont = "MF")    # 分子功能
goBP <- enrichGO(degEntrezId,
                 keyType = "ENTREZID",
                 OrgDb = "org.Hs.eg.db",
                 ont = "BP")  # 生物学过程
goCC <- enrichGO(degEntrezId,
                 keyType = "ENTREZID",
                 OrgDb = "org.Hs.eg.db",
                 ont = "CC") # 细胞组分
goAll <- enrichGO(degEntrezId,
                 keyType = "ENTREZID",
                 OrgDb = "org.Hs.eg.db",
                 ont = "ALL")  # 以上三种              

3.可视化

# 气泡图
dotplot(goBP)  # enrichplot包

# 条形图
barplot(goBP) # enrichplot包

# 有向无环图
goplot(goBP)  # enrichplot包

# 树形图
plotGOgraph(goBP) # clusterprofiler包

# 原型布局
cnetplot(goBP,colorEdge = T)

四.KEGG

kegg和go步骤一致

分析和可视化

degKEGG <- enrichKEGG(degEntrezId,
                      organism = "hsa")
dotplot(degKEGG)

---

总结

可视化结果要能看懂，看横纵坐标参数，当然最重要的还是要补充生物学知识。简单来说所有的分析就是三个步骤：

• 整理数据
  通过读取数据，然后进行清洗，至于数据格式最终整理成什么样子还是要看分析用的函数。
  通过？函数可以查看需要的数据类型从而进行转换、提取等操作。
• 分析
  每一个分析，不管是差异分析还是富集分析都可以查看R包说明文档，然后一步一步操作。
• 可视化
  比较好用的包是ggplot2，不需要参数都记得，使用的时候查看函数就行，记得汇总以备下次使用。
  ！常看R包的说明文档和?函数
！常看R包的说明文档和?函数
！常看R包的说明文档和?函数