GOKEGG自定义通路的富集方法

目录

前言

一、数据背景

二、使用步骤

1.加载所需的R包

2.读入数据与处理表格

3.enrichGO函数进行GO/KEGG/自定义通路的富集

4.简单的可视化

三、结论

---

前言

clusterProfiler 是业界大神Y叔写的一个R包，可以用来做各种富集分析，如GO、KEGG、以及GSEA富集分析等，并且对富集分析结果进行可视化。

这里将使用clusterProfiler包对一些数据进行GO、KEGG等富集分析，并对富集结果可视化。

一、数据背景

本次使用到两组预后不同的临床样本的基因，希望通过富集通路来比较两者的差异通路。

二、使用步骤


1.加载所需的R包

library(clusterProfiler) #clusterProfiler的主体
library(org.Hs.eg.db) #人类的参考基因组数据包
library(stringr) #处理表格数据的包
library(msigdbr) #对GSEA官网的通路数据

2.读入数据与处理表格

mutation <- read.csv('Mutation_list_Final.anno.hg19_multianno.final.xls',
                     header = T,sep = '\t') #读取数据
gene <- mutation[,c('patient','Gene_refGene')] #提取病人样本编号和基因列表
> head(gene)
  patient Gene_refGene
1  S01001       ARID1A
2  S01001         FLT3
3  S01001          RB1
4  S01001          RB1
5  S01001         ARAF
6  S01001        ERBB2

gene_id <- bitr(gene$Gene_refGene,fromType = 'SYMBOL', #数据源，fromtype是基因的id类型
           toType = c('ENTREZID'), #totype指需要转换成的ID类型
           OrgDb='org.Hs.eg.db') #基因ID转换，bitr函数可以将基因名转换为基因ID，这里的ENTREZID是clusterProfiler常用的ID。
> head(gene_id)
  SYMBOL ENTREZID
1 ARID1A     8289
2   FLT3     2322
3    RB1     5925
4   ARAF      369
5  ERBB2     2064
6 PDGFRA     5156

effect <- read.csv('effect.csv') #读取病人分组数据。
> head(effect)
     ID siteBestResponse SD.24wk C1有无血
1 01001               PD      NA       NA
2 01002               PR      NA        0
3 01003               SD       1       NA
4 01004               PR      NA        0
5 01006               PR      NA        0
6 01007               PR      NA       NA

gene$patient <- str_replace(gene$patient, "S","") #这里可以看到gene和effect的patientID是不一样的，这里将S去除。
effect <- effect[,-3:-4] #去掉后两列无用数据
PR <- effect[c(which(effect$siteBestResponse == 'PR')),] #提取出所有疗效为PR的分组
SD_PD <- rbind(effect[c(which(effect$siteBestResponse == 'SD')),],
               effect[c(which(effect$siteBestResponse == 'PD')),]) #提取出所有疗效为PD和SD的分组

colnames(gene)[2] <- 'SYMBOL' #修改gene第二列的列名为SYMBOL
gene <-  merge(gene,gene_id, by = 'SYMBOL',all.x = T) #合并gene和gene_id表格
  SYMBOL patient ENTREZID
1   AKT1   01010      207
2   AKT1   01027      207
3   AKT2   01022      208
4   AKT2   01023      208
5   AKT2   01022      208
6   AKT2   01027      208

colnames(SD_PD)[1] <- 'patient' #修改列名
SD_PD <-  merge(gene,SD_PD, by = 'patient',all.x = T) #通过patient列合并gene组和SD_PD组来标记所有gene中的SD_PD组的数据
colnames(PR)[1] <- 'patient'  #同上
PR <-  merge(gene,PR, by = 'patient',all.x = T) #同上
PR <- PR[which(PR$siteBestResponse == 'PR'),] #提取PR组的数据
SD_PD <- SD_PD[c(which(SD_PD$siteBestResponse == 'PD'),
                 which(SD_PD$siteBestResponse == 'SD')),] #提取SD_PD组的数据

3.enrichGO函数进行GO/KEGG/自定义通路的富集

#进行GO富集
PR_GO <- enrichGO(PR$ENTREZID, #数据源
                  pvalueCutoff = 1, #P值阈值
                  qvalueCutoff = 1, #qvalue是P值的校正值，P值会过滤掉很多，可以全部输出
                  OrgDb = org.Hs.eg.db, #人类参考基因组
                  ont = "ALL", #主要的分为三种，三个层面来阐述基因功能，生物学过程（BP），细胞组分（CC），分子功能（MF）
                  readable = TRUE) #是否将基因ID转换为基因名
write.csv(as.data.frame(PR_GO),"PR_GO.csv",row.names =FALSE) #保存结果
SD_PD_GO <- enrichGO(SD_PD$ENTREZID,pvalueCutoff = 1,qvalueCutoff = 1,
                     OrgDb = org.Hs.eg.db,ont = "ALL",readable = TRUE)
write.csv(as.data.frame(SD_PD_GO),"SD_PD_GO.csv",row.names =FALSE)

#进行KEGG富集
PR_KEGG <- enrichKEGG(PR$ENTREZID, #数据源
                      organism = 'hsa', #物种
                      keyType = 'kegg', #"kegg"/'ncbi-geneid'/'ncib-proteinid'/'uniprot'之一，KEGG就写kegg
                      pvalueCutoff = 1,
                      pAdjustMethod = 'BH',#P值校正方法
                      qvalueCutoff = 1)
write.csv(as.data.frame(PR_KEGG),"PR_KEGG.csv",row.names =FALSE)
SD_PD_KEGG <- enrichKEGG(SD_PD$ENTREZID, organism = 'hsa', keyType = 'kegg', pvalueCutoff = 1,
                       pAdjustMethod = 'BH',qvalueCutoff = 1)
write.csv(as.data.frame(SD_PD_KEGG),"SD_PD_KEGG.csv",row.names =FALSE)

library(msigdbr) #加载GSEA官网的通路数据包
DmGO <- msigdbr(species="Homo sapiens", #物种名
                category="C2") #选择目录，可以通过官网查询自己想富集的通路的目录
PID_pathway <- DmGO[c(which(DmGO$gs_subcat == 'CP:PID'), 
                      which(DmGO$gs_subcat == 'CP')),] #通过自己需要富集的ID号提取通路
PID <- enricher(gene$SYMBOL,TERM2GENE=PID_pathway[,c(3,7)],pvalueCutoff = 1) #富集指定通路集的通路,自定义通路富集需要使用基因名，即SYMBOL
head(PID)
write.csv(as.data.frame(PID),"PID.csv",row.names =FALSE)

4.简单的可视化

#绘制条形图
barplot(ego_ALL, #数据源
        showCategory=20,#展示前20条通路
        drop=T)
#绘制点状图
dotplot(go,showCategory=50) 

三、结论
Y数的包果然还是非常强的，除了富集功能，包中还提供了可视化功能，非常易于使用。同时还提供了自定义通路集的富集功能，适用性很广。