gokegg_差异基因的GO与KEGG注释

写在前面

这个其实很简单啦&＃xff01;三个R包可以搞定的事情。

三个包分是&＃xff1a;clusterProfiler&＃xff0c;pathview&＃xff0c;org.Hs.eg.db。

clusterProfiler&＃xff0c;pathview两个包用于绘图&＃xff0c;org.Hs.eg.db是用于clusterProfiler注释的人类的数据库。

先说安装&＃xff0c;我这里是R-3.4.4的版本&＃xff0c;用如下的安装方式&＃xff0c;3.5&＃43;版本的请见官网哦&＃xff1a;

source("https://bioconductor.org/biocLite.R")

biocLite("pathview")

library(pathview)

这里拿一个包举例子&＃xff0c;其它请大家自己举一反三。要碎碎念一句&＃xff0c;biocLite安装的时候包名是必须要有双引号的哦&＃xff0c;执行library的时候就很随意啦&＃xff0c;单引号&＃xff0c;双引号&＃xff0c;不加都OK。

开始分析

做差异基因的注释&＃xff0c;只需要差异基因的entrez_ID。这一步用到的包是clusterProfiler和org.Hs.eg.db。

我的输入表格格式如图所示(数据都是随便编的&＃xff0c;毕竟不能泄露真实数据)&＃xff1a;

>head(data)

Gene ABCDEF

ENSG00000067082.14100115108300303290

ENSG00000134852.1499101100700710690

ENSG00000125538.11300032003100100012001100

ENSG00000006831.920151510010595

ENSG00000285441.199101100700710690

ENSG00000141905.18300032003100100012001100

分析的命令如下&＃xff1a;

library(clusterProfiler)

library(org.Hs.eg.db)

#Ensembl ID的后缀是版本号&＃xff0c;分析的时候不需要&＃xff0c;所以这里把后缀去掉&＃xff0c;用包内的bitr函数将Ensembl ID转换为ENTREZIDdata[,1]

eg

genelist

#GO annotationgo

pvalueCutoff &＃61; 0.05, qvalueCutoff &＃61; 0.2,keyType &＃61; &＃39;ENTREZID&＃39;)

write.csv(go,"GO-All-enrich.csv",row.names &＃61;FALSE)

#KEGG annotationkegg

pAdjustMethod &＃61; &＃39;BH&＃39;, minGSSize &＃61; 10,maxGSSize &＃61; 500,

qvalueCutoff &＃61; 0.2,use_internal_data &＃61; FALSE)

write.csv(kegg,"KEGG-enrich.csv",row.names &＃61;FALSE)

脚本中enrichGO和enrichKEGG两个函数来完成GO和KEGG的富集分析。

enrichGO中的关注一下参数"ont"&＃xff0c;这个可以设置为"ALL", "CC", "BP", "MF"四者中的一种。因为GO的注释包含后三个种类。其它参数更改的话&＃xff0c;请大家自己看包内的参数。这个主语句很简单啦请大家不要偷懒啦。

图形化展示

富集结果的图形化展示命令如下&＃xff1a;

barplot(go,showCategory&＃61;20,drop&＃61;T,title &＃61; "EnrichmentGO",font.size &＃61; 10 )

dotplot(go,showCategory&＃61;20,title&＃61;"EnrichmentGO_dot")

dotplot(kegg, showCategory&＃61;20,title&＃61;"EnrichmentKEGG")

前两行是展示GO注释的结果&＃xff0c;第三行展示KEGG的注释结果&＃xff1a;

如果你们看到了黄色的区域&＃xff0c;不要惊慌&＃xff0c;不是系统的bug&＃xff0c;是我出于自身考虑遮掉了哦。

题图展示的图片&＃xff0c;是我自己根据结果加工过的(虽然也没有好看到哪里去)。如果只用于自己了解注释信息的话&＃xff0c;直接用上面简单粗暴的语句生成图就好啦&＃xff0c;我觉得也蛮好看的呢。柱状图展示GO富集结果点图展示KEGG富集结果

KEGG会富集到很多通路上&＃xff0c;我们可以根据KEGG的注释结果导出每一个通路的图哦&＃xff0c;图片大概长这样&＃xff1a;

这是官网的图啦&＃xff0c;我自己的只捕捉到了一两个基因&＃xff0c;完全没有这么好看&＃xff0c;还是给大家放标准结果。

这个时候需要导入包pathview啦&＃xff0c;命令在这里&＃xff1a;

上文中KEGG-enrich.csv的格式如下图&＃xff0c;我们需要用到ID和geneID两列信息。

library(pathview)

pathview(gene.data &＃61; as.numeric(strsplit(as.character(test$geneID[1]),"/")[[1]]),

pathway.id &＃61; "04657", species &＃61; "hsa", kegg.native &＃61; T)

gene.data就是富集到这个通路中的geneID&＃xff0c;因为我的表格里是斜杠分割的&＃xff0c;所以我做了一下数据的提取&＃xff0c;大家如果只有一列的&＃xff0c;直接用就好了。

所有的命令行都是固定的&＃xff0c;主要是把自己的数据转换为函数所需的格式。说到这里&＃xff0c;就要结束啦。

参考资料